diff --git a/CLAUDE.md b/CLAUDE.md index dffdae54..cb84d95b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -91,7 +91,10 @@ API points worth knowing: case rather than crashing). - `GenomicArray.autosomes()` falls back to returning the whole array (with a warning) when no autosomes are recognized. Be permissive on - unfamiliar assemblies; don't silently drop data. + unfamiliar assemblies; don't silently drop data. Its mask, + `GenomicArray.autosome_filter()`, deliberately has no fallback: use it + where "no autosomes" must mean "nothing to assess" (e.g. the + low-coverage warning in `segmentation`). - PAR coordinates live in `skgenome.genomebuild`; `cnvlib.params.PSEUDO_AUTSOMAL_REGIONS` is a back-compat re-export. diff --git a/cnvlib/batch.py b/cnvlib/batch.py index c7c50d6d..85c04280 100644 --- a/cnvlib/batch.py +++ b/cnvlib/batch.py @@ -77,7 +77,7 @@ def batch_make_reference( is_haploid_x_reference : bool True if reference samples are male (haploid X chromosome). diploid_parx_genome : str, optional - Reference genome name (e.g., 'hg19', 'hg38') for pseudo-autosomal region + Reference genome name (e.g., 'grch37', 'grch38') for pseudo-autosomal region handling on X/Y chromosomes. fasta : str Path to reference genome FASTA file. Required for GC/RepeatMasker @@ -178,7 +178,7 @@ def batch_make_reference( ... target_bed='baits.bed', ... antitarget_bed=None, # Auto-generate ... is_haploid_x_reference=False, - ... diploid_parx_genome='hg38', + ... diploid_parx_genome='grch38', ... fasta='hg38.fa', ... annotate='refFlat.txt', ... short_names=True, @@ -200,7 +200,7 @@ def batch_make_reference( ... target_bed=None, # Auto-generate from FASTA ... antitarget_bed=None, ... is_haploid_x_reference=True, - ... diploid_parx_genome='hg19', + ... diploid_parx_genome='grch37', ... fasta='hg19.fa', ... annotate='refFlat.txt', ... short_names=False, diff --git a/cnvlib/call.py b/cnvlib/call.py index 50f6e481..fdd9191d 100644 --- a/cnvlib/call.py +++ b/cnvlib/call.py @@ -88,7 +88,7 @@ def do_call( Default: False diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). Treats PAR regions as diploid + (e.g., 'grch37', 'grch38'). Treats PAR regions as diploid even in male samples. Default: None filters : list of str, optional diff --git a/cnvlib/cli/genome_instability_index.py b/cnvlib/cli/genome_instability_index.py index f59d4470..ec27ecbe 100755 --- a/cnvlib/cli/genome_instability_index.py +++ b/cnvlib/cli/genome_instability_index.py @@ -32,7 +32,7 @@ def argument_parsing(): AP.add_argument( "--diploid-parx-genome", type=str, - help="Considers the given human genome's PAR of chromosome X as autosomal. Example: 'grch38'", + help="Considers the given human genome's PAR of chromosome X as autosomal. Example: 'grch38' (or the UCSC name 'hg38')", ) AP.add_argument("-o", "--output", help="Output filename (default: stdout).") return AP.parse_args() diff --git a/cnvlib/cnary.py b/cnvlib/cnary.py index 9d7d20be..b6684cea 100644 --- a/cnvlib/cnary.py +++ b/cnvlib/cnary.py @@ -341,24 +341,29 @@ def center_all( logging.info("Shifting log2 values by %f", shift) self.data["log2"] += shift - def drop_low_coverage(self, verbose: bool = False) -> CopyNumArray: - """Drop bins with extremely low log2 coverage or copy ratio values. + def low_coverage_mask(self) -> Series: + """Flag bins with extremely low log2 coverage or copy ratio values. These are generally bins that had no reads mapped due to sample-specific issues. A very small log2 ratio or coverage value may have been substituted to avoid domain or divide-by-zero errors. NaN log2/depth values (from degenerate or malformed inputs) are also - dropped here: a bare ``log2 < min_cvg`` comparison is False for NaN, so + flagged: a bare ``log2 < min_cvg`` comparison is False for NaN, so without this NaN bins would slip past into segmentation and trigger an 'invalid value encountered in greater' warning or a CBS crash (#521). """ min_cvg = params.NULL_LOG2_COVERAGE - params.MIN_REF_COVERAGE log2 = self.data["log2"] - drop_idx = (log2 < min_cvg) | log2.isna() + mask = (log2 < min_cvg) | log2.isna() if "depth" in self: depth = self.data["depth"] - drop_idx |= (depth == 0) | depth.isna() + mask |= (depth == 0) | depth.isna() + return mask + + def drop_low_coverage(self, verbose: bool = False) -> CopyNumArray: + """Drop the bins flagged by `low_coverage_mask`.""" + drop_idx = self.low_coverage_mask() if verbose and drop_idx.any(): logging.info("Dropped %d low-coverage bins", drop_idx.sum()) return self[~drop_idx] diff --git a/cnvlib/commands.py b/cnvlib/commands.py index 9a1226c6..615b1e6d 100644 --- a/cnvlib/commands.py +++ b/cnvlib/commands.py @@ -105,7 +105,7 @@ def add_diploid_parx_genome(P: argparse._ActionsContainer) -> None: P.add_argument( "--diploid-parx-genome", type=str, - help="Considers the given human genome's PAR of chromosome X as autosomal. Example: 'grch38'", + help="Considers the given human genome's PAR of chromosome X as autosomal. Example: 'grch38' (or the UCSC name 'hg38')", ) @@ -2349,7 +2349,7 @@ def do_sex( Whether the reference sample is male (haploid X chromosome). diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). + (e.g., 'grch37', 'grch38'). Returns ------- diff --git a/cnvlib/fix.py b/cnvlib/fix.py index 370624c3..9983a86f 100644 --- a/cnvlib/fix.py +++ b/cnvlib/fix.py @@ -52,7 +52,7 @@ def do_fix( Contains expected neutral coverage and bias correction features. diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). Ensures PAR regions on X/Y chromosomes + (e.g., 'grch37', 'grch38'). Ensures PAR regions on X/Y chromosomes are treated as diploid during normalization. Default: None do_gc : bool, optional @@ -140,7 +140,7 @@ def do_fix( For male sample with PAR handling: >>> fixed = do_fix( ... target_raw, antitarget_raw, reference, - ... diploid_parx_genome='hg38' + ... diploid_parx_genome='grch38' ... ) """ # Load, recenter and GC-correct target & antitarget probes separately diff --git a/cnvlib/import_rna.py b/cnvlib/import_rna.py index 08b064f8..0fba4c75 100644 --- a/cnvlib/import_rna.py +++ b/cnvlib/import_rna.py @@ -47,7 +47,7 @@ def do_import_rna( Maximum log2 ratio to cap values at. Default is 3. diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38'). + (e.g., 'grch37', 'grch38'). min_sample_fraction : float, optional Minimum fraction of samples in which a gene must be expressed (count >= 1) for it to be retained. Default 0.5 preserves the legacy diff --git a/cnvlib/params.py b/cnvlib/params.py index 0797d5d4..2e21e337 100644 --- a/cnvlib/params.py +++ b/cnvlib/params.py @@ -18,11 +18,12 @@ ANTITARGET_NAME = "Antitarget" ANTITARGET_ALIASES = (ANTITARGET_NAME, "Background") -# PAR coordinates and the list of supported genome builds are owned by -# skgenome.genomebuild; re-exported here for back-compat with code that -# imports them from cnvlib.params. Inner coordinate values are returned -# as mutable lists (matching the historical type) so any caller that -# treated them as such continues to work. +# PAR coordinates and the canonical names of supported genome builds are +# owned by skgenome.genomebuild; re-exported here for back-compat with code +# that imports them from cnvlib.params. Aliases such as "hg38" are accepted by +# skgenome.genomebuild.get_genome_build but are not keys here. Inner coordinate +# values are returned as mutable lists (matching the historical type) so any +# caller that treated them as such continues to work. from skgenome.genomebuild import REGISTERED_BUILDS as _REGISTERED_BUILDS PSEUDO_AUTSOMAL_REGIONS = { diff --git a/cnvlib/reference.py b/cnvlib/reference.py index 777a89da..843e5fb2 100644 --- a/cnvlib/reference.py +++ b/cnvlib/reference.py @@ -48,7 +48,7 @@ def do_reference_flat( Default is False. diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). + (e.g., 'grch37', 'grch38'). Returns ------- @@ -131,7 +131,7 @@ def do_reference( Default: False (assumes female reference or mixed) diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). PAR regions on X/Y are treated as + (e.g., 'grch37', 'grch38'). PAR regions on X/Y are treated as diploid even in male samples. Default: None female_samples : bool, optional @@ -211,7 +211,7 @@ def do_reference( ... target_fnames=['M1.cnn', 'M2.cnn'], ... fa_fname='hg38.fa', ... female_samples=False, - ... diploid_parx_genome='hg38' + ... diploid_parx_genome='grch38' ... ) """ if antitarget_fnames: diff --git a/cnvlib/reports.py b/cnvlib/reports.py index a25cd97f..c1ec27cb 100644 --- a/cnvlib/reports.py +++ b/cnvlib/reports.py @@ -170,7 +170,7 @@ def do_genemetrics( Whether sample is female. If None, inferred from data. diploid_parx_genome : str, optional Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). + (e.g., 'grch37', 'grch38'). location_stats : list of str, optional Location statistics to compute: 'mean', 'median', 'mode', 'p_ttest'. Default is empty tuple. diff --git a/cnvlib/segmentation/__init__.py b/cnvlib/segmentation/__init__.py index a3acad83..a6af41e6 100644 --- a/cnvlib/segmentation/__init__.py +++ b/cnvlib/segmentation/__init__.py @@ -27,6 +27,12 @@ # HMM method variants, for dispatch (replaces scattered method.startswith("hmm")) _HMM_METHODS = frozenset({"hmm", "hmm-tumor", "hmm-germline"}) +#: Warn when at least this fraction of autosomal bins have no or near-zero +#: coverage. At 1%, placeholder values injected into a clean WGS profile already +#: shift CBS segment means by about 0.2 log2 and multiply the haar segment count +#: about 80-fold, while the targeted test fixtures stay below 0.3%. +LOW_COVERAGE_WARN_FRACTION = 0.01 + def do_segmentation( cnarr: CNA, @@ -52,8 +58,8 @@ def do_segmentation( Segmentation algorithm: 'cbs', 'haar', 'hmm', 'hmm-tumor', or 'hmm-germline'. diploid_parx_genome : str, optional - Reference genome name for pseudo-autosomal region handling - (e.g., 'hg19', 'hg38', 'mm10'). + Reference genome name for pseudo-autosomal region handling, as + registered in `skgenome.genomebuild` (e.g., 'grch38'). threshold : float, optional Significance threshold (for CBS/haar) or smoothing window size (for HMM methods). If None, uses method-specific defaults. @@ -105,6 +111,8 @@ def do_segmentation( "Input has no sample_id set; the segmented output will be " "unlabeled (CLI usage derives sample_id from the input filename)." ) + if not skip_low: + _warn_low_coverage(cnarr) if variants is not None and len(variants): # Process variants in genomic order. by_ranges/baf_by_ranges slice via @@ -200,6 +208,48 @@ def _to_str(s, enc=locale.getpreferredencoding()): # noqa: B008 return s +def _warn_low_coverage(cnarr: CNA) -> None: + """Warn when zero-coverage placeholder bins will distort segmentation. + + Without ``skip_low``, bins with no reads are segmented as though their + placeholder log2 values (near ``params.NULL_LOG2_COVERAGE``) were + measurements. Only autosomal bins are counted, because a female sample's + chrY is legitimately at the floor; an input with no recognizable autosomes + is not assessed. Called once per sample, before the per-arm split, so the + warning is not repeated for each arm. + """ + log2 = cnarr["log2"].to_numpy() + # Non-finite bins are dropped before segmentation whether or not skip_low + # is set, so they do not distort it. + is_counted = cnarr.autosome_filter().to_numpy() & np.isfinite(log2) + is_low = cnarr.low_coverage_mask().to_numpy() & is_counted + n_low = int(is_low.sum()) + n_counted = int(is_counted.sum()) + # ``not n_low`` also covers n_counted == 0 + if not n_low or n_low < LOW_COVERAGE_WARN_FRACTION * n_counted: + return + msg = ( + f"WARNING: {n_low} of {n_counted} autosomal bins " + f"({n_low / n_counted:.1%}) have no or near-zero coverage, so their " + "placeholder log2 values will be segmented as measurements." + ) + if n_low < n_counted: + is_ok = is_counted & ~is_low + shift = n_low / n_counted * (log2[is_ok].mean() - log2[is_low].mean()) + if round(shift, 2) > 0: + msg += ( + f" This lowers segment means by about {shift:.2f} log2 units " + "on average, which can report copy-neutral regions as losses." + ) + msg += ( + " Recomputing coverage with larger bins (see the 'autobin' command) " + "leaves fewer bins without reads. The --drop-low-coverage option " + "excludes these bins, but also removes genuine homozygous deletions " + "in germline and high-purity tumor samples." + ) + logging.warning(msg) + + def _ds( args: tuple[ CNA, diff --git a/doc/germline.rst b/doc/germline.rst index de5778db..1dba3b98 100644 --- a/doc/germline.rst +++ b/doc/germline.rst @@ -26,7 +26,8 @@ samples: - The ``--drop-low-coverage`` option (see :doc:`tumor`) should not be used; it will typically remove germline deep deletions altogether, which is not - desirable. + desirable. If :ref:`segment` warns that many bins have no or near-zero + coverage, recompute coverage with larger bins (see :ref:`autobin`) instead. A worked example diff --git a/doc/pipeline.rst b/doc/pipeline.rst index e35d5629..d5df0383 100644 --- a/doc/pipeline.rst +++ b/doc/pipeline.rst @@ -921,7 +921,17 @@ Bins with a weight of 0 are dropped before segmentation. Additional filters: - ``--drop-low-coverage`` -- drop bins with a read depth of 0 or very close to - 0, i.e. a log2 value suggesting the same. Use with tumor samples. + 0, i.e. a log2 value suggesting the same. Typically used with tumor samples; + see :doc:`tumor` for exceptions. + Without this option, such bins are segmented as though their placeholder + log2 values (below -15) were measurements, which lowers the means of the + segments containing them and fragments haar segmentation. When at least 1% + of autosomal bins are affected, the ``segment`` command logs a warning that + also estimates the resulting downward shift in segment means where it is + appreciable. The usual cause is bins too small for the sequencing depth, as + in shallow whole-genome sequencing; recomputing coverage with larger bins + (see :ref:`autobin`) avoids it without discarding data, whereas this option + also removes genuine homozygous deletions (see :doc:`germline`). - ``--drop-outliers`` -- drop bins with log2 value too far from a rolling average, taking local variability into account. Applied by default. diff --git a/doc/sex.rst b/doc/sex.rst index 62fa260b..8ad3e9f3 100644 --- a/doc/sex.rst +++ b/doc/sex.rst @@ -225,11 +225,20 @@ seem to exhibit a copy number loss there (they should have doubled coverage but haven't). In order to avoid this, the option ``--diploid-parx-genome `` can -be used for all affected commands. The genome build should be "grch38" or -similar. This will cause cnvkit to treat the PAR1/2 on chromosome X as -effectively autosomal (i.e. the expected copy number is the same for females -and males, the sex prediction expects the same coverage on X in PAR1/2 for both -sexes, etc.) +be used for all affected commands. This will cause cnvkit to treat the PAR1/2 +on chromosome X as effectively autosomal (i.e. the expected copy number is the +same for females and males, the sex prediction expects the same coverage on X +in PAR1/2 for both sexes, etc.) + +The supported builds are ``grch37`` and ``grch38``; the UCSC names ``hg19`` and +``hg38`` are accepted as aliases, since their chromosome X and Y coordinates +are identical to GRCh37 and GRCh38, and names are not case-sensitive. The +chromosome naming style (``chrX`` or ``X``) does not matter. The build should +be the one the bin coordinates come from: CNVkit cannot detect a mismatch, and +the PAR coordinates differ between the two builds. PAR1 differs only near its +boundaries, but the two builds' PAR2 regions do not overlap at all, so with the +wrong build the real PAR2 is not treated as autosomal and a stretch of +X-specific sequence is. When ``--male-reference`` (see above) is also used, the reference copy number for PAR1/2 on X is still 2 (same as for the autosomes). @@ -307,8 +316,11 @@ CNVkit is not detecting my sample's sex correctly. What can I do? In lower-quality samples, particularly tumor samples analyzed without a robust reference (see :doc:`tumor`), there may be many bins with no coverage which bias -the segment means. Try repeating the :ref:`segment` command with the -``--drop-low-coverage`` option if you did not do so originally. +the segment means; the :ref:`segment` command warns when this is the case. For +a bulk tumor sample, try repeating the :ref:`segment` command with the +``--drop-low-coverage`` option if you did not do so originally. Otherwise, +recompute coverage with larger bins (see :ref:`autobin`), since that option +also removes genuine homozygous deletions. If a small number of samples in your reference pool look like they were inferred incorrectly (e.g. three males among forty samples are reported as diff --git a/doc/tumor.rst b/doc/tumor.rst index 402f01e9..2b00550a 100644 --- a/doc/tumor.rst +++ b/doc/tumor.rst @@ -9,13 +9,13 @@ are available to support this use case: platform, see the :ref:`reference` command for strategies. - Use ``--drop-low-coverage`` to ignore bins with log2 normalized coverage - values below -15. Virtually all tumor samples, even cancer cell lines, are - not completely homogeneous. Even in regions of homozygous deletion in the - largest tumor-cell clonal population, some sequencing reads will be obtained - from contaminating normal cells without the deletion. Therefore, extremely low - log2 copy ratio values do not indicate homozygous deletions but failed - sequencing or mapping in all cells regardless of copy number status at that - site, which are not informative for copy number. This option in the + values below -15. Most tumor samples are not completely homogeneous. Even + in regions of homozygous deletion in the largest tumor-cell clonal + population, some sequencing reads will be obtained from contaminating normal + cells without the deletion. Therefore, extremely low log2 copy ratio values + usually do not indicate homozygous deletions but failed sequencing or mapping + in all cells regardless of copy number status at that site, which are not + informative for copy number. This option in the :ref:`batch` command applies to segmentation; the option is also available in the :ref:`segment`, :ref:`metrics`, :ref:`segmetrics`, :ref:`genemetrics` and :doc:`heterogeneity` commands. @@ -28,6 +28,12 @@ are available to support this use case: probably based on dummy values corresponding to zero-coverage (perhaps unmappable) bins, and not real observations. + - The rationale above assumes contaminating normal cells. In a sample of + very high tumor purity, such as a flow-sorted population or some cell + lines, a homozygous deletion may yield no reads at all, and this option + then removes it from the segmentation; leave the option off for such + samples. + - The :ref:`batch` command emits default-threshold integer copy number calls in each sample's ``.call.cns`` file, applying no tumor purity or ploidy correction (see :ref:`batch`). For purity-aware integer copy number, re-run diff --git a/skgenome/gary.py b/skgenome/gary.py index 2f5d540b..803296c3 100644 --- a/skgenome/gary.py +++ b/skgenome/gary.py @@ -271,8 +271,8 @@ def sample_id(self) -> pd.Series: # Traversal - def autosomes(self, also: str | list[str] | pd.Series | None = None) -> Self: - """Select autosomal rows, excluding sex chromosomes and alternative contigs. + def autosome_filter(self) -> pd.Series: + """Flag autosomal rows, excluding sex chromosomes and alternative contigs. Accepts both arabic-numeral chromosome names (e.g. ``chr1``) and Roman-numeral names (e.g. ``chrI``, ``chrXVI``). Sex chromosomes are @@ -280,20 +280,28 @@ def autosomes(self, also: str | list[str] | pd.Series | None = None) -> Self: as a sex chromosome in human data and as autosome 10 in yeast data. Mitochondrial and alternative/unplaced contigs are always excluded. - - If no autosomes can be recognized (e.g. an unfamiliar custom assembly), - emits a warning and returns the whole array unchanged. + Unlike `autosomes`, an unfamiliar assembly yields an all-False mask, and + pseudo-autosomal regions are never included (``CopyNumArray.autosomes`` + adds them back when given ``diploid_parx_genome``). """ - unique = self.chromosome.unique() - sex_x, sex_y = infer_sex_chrom_labels(unique) - is_auto = self.chromosome.map( + sex_x, sex_y = infer_sex_chrom_labels(self.chromosome.unique()) + return self.chromosome.map( lambda n: ( is_autosome(n, sex_x, sex_y) and not is_mitochondrial(n) and not is_alternative_contig(n) ) ) + + def autosomes(self, also: str | list[str] | pd.Series | None = None) -> Self: + """Select the rows flagged by `autosome_filter`, plus any in `also`. + + If no autosomes can be recognized (e.g. an unfamiliar custom assembly), + emits a warning and returns the whole array unchanged. + """ + is_auto = self.autosome_filter() if not is_auto.any(): + unique = self.chromosome.unique() logging.warning( "GenomicArray.autosomes(): no autosomes recognized in %d " "chromosomes (sample: %s); returning the array unchanged.", diff --git a/skgenome/genomebuild.py b/skgenome/genomebuild.py index 23a4bb81..1568f645 100644 --- a/skgenome/genomebuild.py +++ b/skgenome/genomebuild.py @@ -7,7 +7,8 @@ For callers that already pass a genome name as a string (e.g. ``diploid_parx_genome="grch38"``), :func:`get_genome_build` performs the -lookup and produces a clear ``KeyError`` for unknown names. +lookup, by canonical name or alias, and produces a clear ``KeyError`` for +unknown names. """ from __future__ import annotations @@ -29,10 +30,17 @@ class GenomeBuild: backing mapping is wrapped in ``MappingProxyType`` on construction so that the module-level singletons (``GRCH37``, ``GRCH38``) cannot be mutated in place by callers. + + ``aliases`` lists further names :func:`get_genome_build` accepts for this + build: other names for the assembly whose chromosome X and Y coordinates + are identical, so the same PAR coordinates apply. UCSC hg19 and hg38 are + aliases of GRCh37 and GRCh38 in this sense. Aliases must be lowercase, + since lookup lowercases its argument. """ name: str par_regions: Mapping[str, tuple[int, int]] + aliases: tuple[str, ...] = () def __post_init__(self) -> None: # Defend the singleton against in-place mutation of the backing dict. @@ -49,6 +57,7 @@ def __post_init__(self) -> None: "PAR1Y": (10000, 2649520), "PAR2Y": (59034049, 59363566), }, + aliases=("hg19",), ) GRCH38 = GenomeBuild( @@ -59,11 +68,21 @@ def __post_init__(self) -> None: "PAR1Y": (10000, 2781479), "PAR2Y": (56887902, 57217415), }, + aliases=("hg38",), ) REGISTERED_BUILDS: dict[str, GenomeBuild] = {b.name: b for b in (GRCH37, GRCH38)} +# Every name lookup accepts: canonical names plus aliases. REGISTERED_BUILDS +# keeps one entry per assembly, because cnvlib.params re-exports its keys as +# the historical PSEUDO_AUTSOMAL_REGIONS table. +_BUILDS_BY_NAME: dict[str, GenomeBuild] = { + name: build + for build in REGISTERED_BUILDS.values() + for name in (build.name, *build.aliases) +} + def get_genome_build(name: str) -> GenomeBuild: """Look up a registered genome build by name (case-insensitive). @@ -72,13 +91,13 @@ def get_genome_build(name: str) -> GenomeBuild: """ key = name.lower() try: - return REGISTERED_BUILDS[key] + return _BUILDS_BY_NAME[key] except KeyError: raise KeyError( - f"Unknown genome build {name!r}; supported: {sorted(REGISTERED_BUILDS)}" + f"Unknown genome build {name!r}; supported: {sorted(_BUILDS_BY_NAME)}" ) from None def is_supported_build(name: str) -> bool: - """Return True if *name* is a registered genome build.""" - return name.lower() in REGISTERED_BUILDS + """Return True if *name* is a registered genome build or alias.""" + return name.lower() in _BUILDS_BY_NAME diff --git a/test/test_genomebuild.py b/test/test_genomebuild.py index 18985629..a763787d 100644 --- a/test/test_genomebuild.py +++ b/test/test_genomebuild.py @@ -32,17 +32,41 @@ def test_get_genome_build_case_insensitive(self): self.assertIs(get_genome_build("grch38"), GRCH38) self.assertIs(get_genome_build("GRCH38"), GRCH38) + def test_ucsc_aliases_resolve_to_canonical_build(self): + """UCSC hg19/hg38 share GRCh37/GRCh38 chrX and chrY coordinates.""" + self.assertIs(get_genome_build("hg19"), GRCH37) + self.assertIs(get_genome_build("hg38"), GRCH38) + self.assertIs(get_genome_build("HG38"), GRCH38) + self.assertEqual(get_genome_build("hg38").name, "grch38") + def test_get_genome_build_unknown_raises(self): + # hg18 is NCBI36, whose PAR coordinates differ; it must not resolve. with self.assertRaises(KeyError) as cm: get_genome_build("hg18") - # Error message lists supported names so the user can self-correct. + # Error message lists accepted names, aliases included, so the user can + # self-correct. self.assertIn("supported", str(cm.exception)) + self.assertIn("hg38", str(cm.exception)) def test_is_supported_build(self): self.assertTrue(is_supported_build("grch37")) self.assertTrue(is_supported_build("GRCh38")) + self.assertTrue(is_supported_build("HG19")) self.assertFalse(is_supported_build("hg18")) + def test_alias_table_is_unambiguous(self): + """Every name is lowercase and resolves to the build that declares it. + + A collision between any two names, canonical or alias, would silently + resolve to whichever build came last and apply the wrong PAR + coordinates; lookup lowercases its argument, so an uppercase alias + could never match. + """ + for build in REGISTERED_BUILDS.values(): + for name in (build.name, *build.aliases): + self.assertEqual(name, name.lower()) + self.assertIs(get_genome_build(name), build) + def test_frozen(self): with self.assertRaises(FrozenInstanceError): GRCH37.name = "other" # type: ignore[misc] @@ -54,8 +78,10 @@ def test_par_regions_singleton_protected_from_mutation(self): def test_back_compat_via_params(self): """cnvlib.params.PSEUDO_AUTSOMAL_REGIONS preserves the historical API.""" - # Every registered build appears as a top-level key + # Every registered build appears as a top-level key, and aliases do not: + # the historical keys were exactly the canonical names. self.assertEqual(set(params.PSEUDO_AUTSOMAL_REGIONS), set(REGISTERED_BUILDS)) + self.assertNotIn("hg38", params.PSEUDO_AUTSOMAL_REGIONS) # All four PAR keys round-trip equally for every build for build in (GRCH37, GRCH38): for key in ("PAR1X", "PAR2X", "PAR1Y", "PAR2Y"): diff --git a/test/test_haar.py b/test/test_haar.py index 5152ebb8..e04e6b81 100644 --- a/test/test_haar.py +++ b/test/test_haar.py @@ -302,9 +302,11 @@ class HaarNoiseFloorTests(unittest.TestCase): FIXTURE = "formats/regression/p2-9_2.cnr" + def _load(self): + return cnvlib.read(self.FIXTURE) + def _n_segments(self, fdr_q): - cnarr = cnvlib.read(self.FIXTURE) - return len(segment_haar(cnarr, fdr_q)) + return len(segment_haar(self._load(), fdr_q)) def test_not_grossly_over_segmented(self): """Default-q haar averages many bins per segment, not ~7. @@ -313,7 +315,7 @@ def test_not_grossly_over_segmented(self): haar produces a few hundred (>=15 bins/seg). Assert the average segment spans at least 15 bins so the ~7-bin fragmentation cannot return. """ - cnarr = cnvlib.read(self.FIXTURE) + cnarr = self._load() n_bins = len(cnarr) n_segs = len(segment_haar(cnarr, 1e-4)) self.assertLessEqual( @@ -342,6 +344,24 @@ def test_fdr_q_controls_sensitivity(self): ) +class HaarWgsNoiseFloorTests(HaarNoiseFloorTests): + """The same noise-floor assertions at WGS bin density. + + Fixture: test/formats/wgs-chr17.cnr, 165,626 bins of median width 267 bp. + Both assertions hold once the zero-coverage placeholder bins (12.9% of + this fixture) are dropped, so the over-segmentation haar shows on the raw + array comes from those bins rather than from bin density. On the raw array + the level-1 MAD is blind to the placeholder values while the coarser + windowed means are not, the FDR threshold falls to about the dispersion of + the statistics it thresholds, and both assertions fail. + """ + + FIXTURE = "formats/wgs-chr17.cnr" + + def _load(self): + return cnvlib.read(self.FIXTURE).drop_low_coverage() + + class SegmentByPeaksTests(unittest.TestCase): """Tests for weighted segment-mean computation.""" diff --git a/test/test_segmentation.py b/test/test_segmentation.py index 0313b2a3..81acb7d8 100644 --- a/test/test_segmentation.py +++ b/test/test_segmentation.py @@ -305,6 +305,111 @@ def test_haar_vcf_without_allele_info_falls_back(self): self.assertTrue((seg.start < seg.end).all()) +class FloorBinWarningTests(unittest.TestCase): + """Zero-coverage placeholder bins are segmented as measurements by default. + + Bins with no reads carry log2 near params.NULL_LOG2_COVERAGE. Unless + skip_low is set they reach the segmenter, where at WGS bin density they + pull CBS segment means down by several log2 units and fragment haar. + Dropping them by default would also remove genuine homozygous deletions + in high-purity samples, so do_segmentation warns instead. + """ + + WGS = "formats/wgs-chr17.cnr" + + @staticmethod + def _synthetic(chrom_sizes, log2, depth=None): + """A CopyNumArray of 1-kb bins with the given per-chromosome bin counts.""" + chroms = [c for c, n in chrom_sizes.items() for _ in range(n)] + starts = np.concatenate([np.arange(n) for n in chrom_sizes.values()]) * 1000 + data = pd.DataFrame( + { + "chromosome": chroms, + "start": starts, + "end": starts + 1000, + "gene": ["-"] * len(chroms), + "log2": log2, + "weight": np.ones(len(chroms)), + } + ) + if depth is not None: + data["depth"] = depth + return cnary.CopyNumArray(data, {"sample_id": "synthetic"}) + + def _warnings(self, cnarr, skip_low=False): + with self.assertLogs(level="INFO") as cm: + segmentation.do_segmentation(cnarr, "none", skip_low=skip_low) + return [line for line in cm.output if "--drop-low-coverage" in line] + + def _assert_silent(self, cnarr, skip_low=False): + self.assertEqual(self._warnings(cnarr, skip_low), []) + + def test_warns_once_on_wgs_dropout(self): + """12.9% of wgs-chr17 bins are at the floor; warn exactly once. + + chr17 has two arms, so a warning emitted per arm would appear twice. + """ + warned = self._warnings(cnvlib.read(self.WGS)) + self.assertEqual(len(warned), 1, warned) + self.assertIn("21361 of 165626", warned[0]) + self.assertIn("lowers segment means", warned[0]) + + def test_shift_stated_only_when_means_drop(self): + """The expected shift is stated only when it is a visible downward shift. + + Bins flagged by zero depth alone keep ordinary log2 values; here they + sit 0.01 below the rest, a shift of 0.002 that would print as 0.00. + When every bin is at the floor there are no measured bins to compare. + """ + depth_only = self._synthetic( + {"chr1": 100}, + np.concatenate([np.full(20, -0.01), np.zeros(80)]), + depth=[0.0] * 20 + [30.0] * 80, + ) + all_floor = self._synthetic({"chr1": 100}, np.full(100, -19.5)) + for cnarr, counts in ((depth_only, "20 of 100"), (all_floor, "100 of 100")): + with self.subTest(counts=counts): + warned = self._warnings(cnarr) + self.assertEqual(len(warned), 1, warned) + self.assertIn(counts, warned[0]) + self.assertNotIn("lowers segment means", warned[0]) + + def test_silent_when_dropping_low_coverage(self): + """With skip_low the floor bins are excluded, so there is nothing to report.""" + self._assert_silent(cnvlib.read(self.WGS), skip_low=True) + + def test_silent_on_targeted_fixtures(self): + """Ordinary panels have at most a few autosomal floor bins.""" + for fname in ( + "formats/amplicon.cnr", + "formats/p2-20_1.cnr", + "formats/regression/p2-9_2.cnr", + ): + with self.subTest(fname=fname): + self._assert_silent(cnvlib.read(fname)) + + def test_sex_chromosome_floor_is_not_counted(self): + """A female sample's chrY sits at the floor legitimately. + + Counting it would warn on every female sample. The second case has no + autosomes, where GenomicArray.autosomes() would return the whole array. + """ + rng = np.random.default_rng(0) + normal = rng.normal(0, 0.1, 80) + for sizes in ({"chr1": 60, "chrX": 20, "chrY": 20}, {"chrX": 80, "chrY": 20}): + with self.subTest(chromosomes=list(sizes)): + log2 = np.concatenate([normal, np.full(20, -19.5)]) + self._assert_silent(self._synthetic(sizes, log2)) + + def test_non_finite_bins_are_not_counted(self): + """Non-finite bins are dropped before segmentation regardless of skip_low. + + -inf is below the floor threshold, so counting it would report 10% here. + """ + log2 = np.concatenate([np.full(10, -np.inf), np.zeros(90)]) + self._assert_silent(self._synthetic({"chr1": 100}, log2)) + + class TransferFieldsTests(unittest.TestCase): """Tests for transfer_fields and do_segmentation NaN handling."""