Warn when zero-coverage bins will bias segmentation

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.
 

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,

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

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()

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]

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
     -------

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

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

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 = {

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:

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.

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,

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

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.
 

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 <genome-build>`` 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

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