Estimate autobin depths from bedGraph input

cnvlib/autobin.py

@@ -1,4 +1,8 @@
-"""Estimate reasonable bin sizes from BAM read counts or depths."""
+"""Estimate reasonable bin sizes from read counts or depths.
+
+The sample may be a BAM or CRAM file, or a tabix-indexed bedGraph (.bed.gz) of
+per-base read depth.
+"""
 
 from __future__ import annotations
 
@@ -20,6 +24,12 @@
 from .antitarget import compare_chrom_names
 from .descriptives import weighted_median
 
+# A bedGraph is read line by line, so each region sampled from it is trimmed to a
+# centered window of at most this many bases. Accessible and antitarget regions
+# span up to megabases each, and 100 of them untrimmed would mean reading tens of
+# megabases of per-base depth; targets are almost always shorter than this.
+BEDGRAPH_SAMPLE_WIDTH = 10000
+
 
 def midsize_file(fnames: Sequence[str]) -> str:
     """Select the median-size file from several given filenames.
@@ -42,19 +52,22 @@ def do_autobin(
     antitarget_max_size: int = 1000000,
     fasta: str | None = None,
 ) -> tuple[tuple[float, int | None], tuple[float | None, int | None]]:
-    """Quickly calculate reasonable bin sizes from BAM read counts.
+    """Quickly calculate reasonable bin sizes from read counts or depths.
 
     Parameters
     ----------
     bam_fname : str
-        Path to BAM file.
+        Path to a BAM or CRAM file, or to a tabix-indexed bedGraph file
+        (.bed.gz) of per-base read depth.
     method : str
         Sequencing method: 'wgs' (whole-genome sequencing), 'amplicon' (targeted
         amplicon capture), or 'hybrid' (hybridization capture).
     targets : GenomicArray, optional
         Targeted genomic regions (required for 'hybrid' and 'amplicon').
     access : GenomicArray, optional
-        Sequencing-accessible regions of the reference genome (for 'hybrid' and 'wgs').
+        Sequencing-accessible regions of the reference genome (for 'hybrid' and
+        'wgs'). Required for 'wgs' on a bedGraph, and for a bedGraph's 'hybrid'
+        antitarget depth.
     bp_per_bin : float, optional
         Desired number of sequencing read nucleotide bases mapped to each bin.
         Default is 100000.0.
@@ -78,7 +91,21 @@ def do_autobin(
     Raises
     ------
     ValueError
-        If targets are required for the method but not provided or empty.
+        If targets are required for the method but not provided or empty; if
+        method 'wgs' is given a bedGraph without accessible regions; or if no
+        read depth is found in the target regions.
+
+    Notes
+    -----
+    A BAM or CRAM file is summarized by its index: per-chromosome read counts
+    and the estimated read length give the 'wgs' depth and the 'hybrid'
+    antitarget depth. A bedGraph records neither read counts nor contig
+    lengths, so every depth is instead measured in a seeded random sample of
+    regions -- the targets, the accessible regions ('wgs'), or the accessible
+    regions on targeted chromosomes minus the targets ('hybrid' antitargets) --
+    as the target depth is for any input. The 'hybrid' target depth of a
+    bedGraph covers all targets, whereas for a BAM it covers only those on
+    accessible chromosomes.
     """
     if method in ("amplicon", "hybrid"):
         if targets is None:
@@ -113,6 +140,43 @@ def depth2binsize(depth: float | None, min_size: int, max_size: int) -> int | No
             bin_size = max_size
         return bin_size
 
+    if access is not None and not len(access):
+        # Empty accessible regions mean none given, as on the command line
+        access = None
+    if coverage.is_bedgraph_format(bam_fname):
+        tgt_depth, anti_depth = bedgraph_depths(bam_fname, method, targets, access)
+    else:
+        tgt_depth, anti_depth = bam_depths(bam_fname, method, targets, access, fasta)
+
+    # A non-finite or non-positive target depth means no reads were found in
+    # the sampled/accessible regions -- give an actionable error rather than a
+    # cryptic "cannot convert float NaN to integer" downstream (#421).
+    if tgt_depth is None or not np.isfinite(tgt_depth) or tgt_depth <= 0:
+        raise ValueError(
+            f"Could not estimate read depth in the target regions of "
+            f"{bam_fname!r} (method {method!r}): got {tgt_depth!r}. This "
+            "usually means no reads mapped to those regions -- most often "
+            "because the chromosome names in the sample file don't match those "
+            "in the reference, --access, or --targets regions (e.g. 'chr1' "
+            "vs. '1'), or the regions are empty. Check that the sample file, "
+            "reference FASTA, and any --access/--targets files use consistent "
+            "chromosome names."
+        )
+
+    # Clip bin sizes to specified ranges
+    tgt_bin_size = depth2binsize(tgt_depth, target_min_size, target_max_size)
+    anti_bin_size = depth2binsize(anti_depth, antitarget_min_size, antitarget_max_size)
+    return ((tgt_depth, tgt_bin_size), (anti_depth, anti_bin_size))
+
+
+def bam_depths(
+    bam_fname: str,
+    method: str,
+    targets: GA | None,
+    access: GA | None,
+    fasta: str | None,
+) -> tuple[float, float | None]:
+    """Target and antitarget read depths of a BAM or CRAM file."""
     samutil.ensure_bam_index(bam_fname)
     rc_table = samutil.idxstats(bam_fname, drop_unmapped=True, fasta=fasta)
     read_len = samutil.get_read_length(bam_fname, fasta=fasta)
@@ -121,10 +185,6 @@ def depth2binsize(depth: float | None, min_size: int, max_size: int) -> int | No
     # Dispatch
     match method:
         case "amplicon":
-            # From BAM index
-            # rc_table = update_chrom_length(rc_table, targets)
-            # tgt_depth = average_depth(rc_table, read_len)
-            # By sampling
             assert targets is not None
             tgt_depth = sample_region_cov(bam_fname, targets, fasta=fasta)
             anti_depth: float | None = None
@@ -139,30 +199,55 @@ def depth2binsize(depth: float | None, min_size: int, max_size: int) -> int | No
                 fasta,
             )
         case "wgs":
-            if access is not None and len(access):
-                rc_table = update_chrom_length(rc_table, access)
+            rc_table = update_chrom_length(rc_table, access)
             tgt_depth = average_depth(rc_table, read_len)  # type: ignore[arg-type]
-            anti_depth = None  # type: ignore[assignment]
+            anti_depth = None
+        case _:
+            raise ValueError(f"Unknown sequencing method {method!r}")
+    return tgt_depth, anti_depth
 
-    # A non-finite or non-positive target depth means no reads were found in
-    # the sampled/accessible regions -- give an actionable error rather than a
-    # cryptic "cannot convert float NaN to integer" downstream (#421).
-    if tgt_depth is None or not np.isfinite(tgt_depth) or tgt_depth <= 0:
-        raise ValueError(
-            f"Could not estimate read depth in the target regions of "
-            f"{bam_fname!r} (method {method!r}): got {tgt_depth!r}. This "
-            "usually means no reads mapped to those regions -- most often "
-            "because the chromosome names in the BAM file don't match those "
-            "in the reference, --access, or --targets regions (e.g. 'chr1' "
-            "vs. '1'), or the regions are empty. Check that the BAM, reference "
-            "FASTA, and any --access/--targets files use consistent chromosome "
-            "names."
-        )
 
-    # Clip bin sizes to specified ranges
-    tgt_bin_size = depth2binsize(tgt_depth, target_min_size, target_max_size)
-    anti_bin_size = depth2binsize(anti_depth, antitarget_min_size, antitarget_max_size)
-    return ((tgt_depth, tgt_bin_size), (anti_depth, anti_bin_size))
+def bedgraph_depths(
+    bg_fname: str, method: str, targets: GA | None, access: GA | None
+) -> tuple[float, float | None]:
+    """Target and antitarget read depths of a bedGraph, by region sampling."""
+    match method:
+        case "amplicon":
+            assert targets is not None
+            return sample_region_cov(bg_fname, targets), None
+        case "hybrid":
+            assert targets is not None
+            tgt_depth = sample_region_cov(bg_fname, targets)
+            if access is None:
+                logging.warning(
+                    "Antitarget depth can't be estimated from bedGraph %s without "
+                    "sequencing-accessible regions (--access); a bedGraph does "
+                    "not record chromosome lengths",
+                    bg_fname,
+                )
+                return tgt_depth, None
+            # Off-target regions on the targeted chromosomes, as hybrid() takes
+            # for BAM. Unlike hybrid(), which subtracts target reads from each
+            # chromosome's total, the target depth needs no such restriction.
+            antitargets = shared_chroms(access.subtract(targets), targets)[0]
+            if not len(antitargets):
+                logging.warning(
+                    "No accessible regions lie off-target on the targeted "
+                    "chromosomes; check that the --access and --targets "
+                    "chromosome names match. Antitarget depth can't be estimated."
+                )
+                return tgt_depth, None
+            return tgt_depth, sample_region_cov(bg_fname, antitargets)
+        case "wgs":
+            if access is None:
+                raise ValueError(
+                    f"Method 'wgs' needs sequencing-accessible regions (--access) "
+                    f"to estimate the read depth of bedGraph {bg_fname!r}, since a "
+                    "bedGraph does not record chromosome lengths."
+                )
+            return sample_region_cov(bg_fname, access), None
+        case _:
+            raise ValueError(f"Unknown sequencing method {method!r}")
 
 
 def hybrid(
@@ -222,15 +307,30 @@ def idxstats2ga(table: pd.DataFrame, bam_fname: str) -> GA:
 def sample_region_cov(
     bam_fname: str, regions: GA, max_num: int = 100, fasta: str | None = None
 ) -> float:
-    """Calculate read depth in a randomly sampled subset of regions."""
+    """Calculate read depth in a randomly sampled subset of regions.
+
+    `bam_fname` may be a BAM or CRAM file, or a bedGraph (.bed.gz), whose
+    sampled regions are first trimmed to `BEDGRAPH_SAMPLE_WIDTH` bases.
+    """
+    if not (regions.end > regions.start).any():
+        # No targets on the accessible chromosomes, or only zero-width regions:
+        # depth is undefined, and do_autobin raises an actionable error.
+        return float("nan")
     midsize_regions = sample_midsize_regions(regions, max_num)
-    bam_chroms = samutil.get_bam_chroms(bam_fname, fasta)
+    is_bedgraph = coverage.is_bedgraph_format(bam_fname)
+    if is_bedgraph:
+        midsize_regions = center_windows(midsize_regions, BEDGRAPH_SAMPLE_WIDTH)
+    else:
+        bam_chroms = samutil.get_bam_chroms(bam_fname, fasta)
     with tempfile.NamedTemporaryFile(suffix=".bed", mode="w+t") as f:
         tabio.write(regions.as_dataframe(midsize_regions), f, "bed4")
         f.flush()
-        table = coverage.bedcov(
-            f.name, bam_fname, 0, fasta, allow_empty=True, bam_chroms=bam_chroms
-        )
+        if is_bedgraph:
+            table = coverage.bedgraph_to_basecount(f.name, bam_fname)
+        else:
+            table = coverage.bedcov(
+                f.name, bam_fname, 0, fasta, allow_empty=True, bam_chroms=bam_chroms
+            )
     if not len(table):
         # No sampled region overlaps a BAM contig -- depth is undefined; let
         # the caller (do_autobin) raise an actionable chrom-name-mismatch error.
@@ -240,15 +340,28 @@ def sample_region_cov(
 
 
 def sample_midsize_regions(regions: GA, max_num: int) -> pd.DataFrame:
-    """Randomly select a subset of up to `max_num` regions."""
+    """Randomly select a subset of up to `max_num` regions.
+
+    Regions between the 25th and 75th percentile of size are preferred. When no
+    region falls in that range, as with two regions of very different size, all
+    regions of nonzero size are eligible instead.
+    """
     sizes = regions.end - regions.start
     lo_size, hi_size = np.percentile(sizes[sizes > 0], [25, 75])
-    midsize_regions = regions.data[(sizes >= lo_size) & (sizes <= hi_size)]
+    in_range = (sizes >= lo_size) & (sizes <= hi_size)
+    midsize_regions = regions.data[in_range if in_range.any() else sizes > 0]
     if len(midsize_regions) > max_num:
         midsize_regions = midsize_regions.sample(max_num, random_state=0xA5EED)
     return midsize_regions
 
 
+def center_windows(regions: pd.DataFrame, width: int) -> pd.DataFrame:
+    """Trim each region longer than `width` to its centered `width`-base window."""
+    lengths = regions.end - regions.start
+    starts = regions.start + (lengths - width).clip(lower=0) // 2
+    return regions.assign(start=starts, end=starts + lengths.clip(upper=width))
+
+
 def shared_chroms(*tables) -> list:
     """Intersection of DataFrame .chromosome values."""
     chroms = tables[0].chromosome.drop_duplicates()

cnvlib/batch.py

@@ -258,7 +258,7 @@ def batch_make_reference(
         # Tweak default parameters
         if not target_avg_size:
             if normal_fnames:
-                # Calculate bin size from .bai & access
+                # Calculate bin size from the BAM index or bedGraph, & access
                 if fasta and not access_arr:
                     # Calculate wgs depth from all
                     # sequencing-accessible area (it doesn't take that long
@@ -293,11 +293,15 @@ def batch_make_reference(
 
     elif method == "hybrid":
         if (not target_avg_size or not antitarget_avg_size) and normal_fnames:
-            # Calculate bin sizes from BAM read counts
+            # Calculate bin sizes from read depths in the median-size normal
             assert target_bed is not None
             bait_arr = tabio.read_auto(target_bed)
             access_arr = tabio.read_auto(access_bed) if access_bed else None
             bam_fname = autobin.midsize_file(normal_fnames)
+            if not access_arr and fasta and coverage.is_bedgraph_format(bam_fname):
+                # A bedGraph has no chromosome lengths, which bound the
+                # off-target regions of a BAM; take accessible regions from FASTA
+                access_arr = access.do_access(fasta)
             (tgt_depth, tgt_bin_size), (anti_depth, anti_bin_size) = autobin.do_autobin(
                 bam_fname,
                 "hybrid",

cnvlib/commands.py

@@ -749,7 +749,7 @@ def _cmd_antitarget(args: argparse.Namespace) -> None:
 
 
 def _cmd_autobin(args: argparse.Namespace) -> None:
-    """Quickly calculate reasonable bin sizes from BAM or CRAM read counts."""
+    """Quickly calculate reasonable bin sizes from BAM, CRAM, or bedGraph depths."""
     if args.method in ("hybrid", "amplicon") and not args.targets:
         raise RuntimeError(f"Sequencing method {args.method!r} requires targets (-t)")
     if args.method == "wgs":
@@ -830,8 +830,8 @@ def read_regions(bed_fname):
 P_autobin.add_argument(
     "bams",
     nargs="+",
-    help="""Sample BAM or CRAM file(s), or bedGraph file(s) (.bed.gz), to test for
-            target coverage""",
+    help="""Sample BAM or CRAM file(s), or bedGraph file(s) (.bed.gz with tabix index
+            .tbi or .csi), to test for target coverage""",
 )
 add_fasta(P_autobin)
 P_autobin.add_argument(
@@ -848,7 +848,9 @@ def read_regions(bed_fname):
     "--access",
     metavar="FILENAME",
     help="""Sequencing-accessible genomic regions, or exons to use as possible targets
-            (e.g. output of refFlat2bed.py)""",
+            (e.g. output of refFlat2bed.py). For bedGraph input, read depth is
+            measured in these regions for 'wgs' (required) and outside the targets
+            for the 'hybrid' antitarget bin size.""",
 )
 P_autobin.add_argument(
     "-t",

cnvlib/coverage.py

@@ -1,4 +1,4 @@
-"""Supporting functions for the 'antitarget' command."""
+"""Supporting functions for the 'coverage' command."""
 
 from __future__ import annotations
 
@@ -46,19 +46,21 @@ def is_bedgraph_format(fname: str) -> bool:
     return fname.endswith(".bed.gz")
 
 
-def validate_bedgraph_index(fname: str) -> None:
-    """Ensure bedGraph file has a tabix index.
+def find_bedgraph_index(fname: str) -> str:
+    """Locate a bedGraph file's tabix index, preferring .tbi over .csi.
 
     Raises
     ------
     FileNotFoundError
         If the index file (.tbi or .csi) is not found.
     """
-    if not (os.path.exists(f"{fname}.tbi") or os.path.exists(f"{fname}.csi")):
-        raise FileNotFoundError(
-            f"bedGraph file {fname} requires a tabix index (.tbi or .csi). "
-            f"Create one with: tabix -p bed {fname}"
-        )
+    for index_fname in (f"{fname}.tbi", f"{fname}.csi"):
+        if os.path.exists(index_fname):
+            return index_fname
+    raise FileNotFoundError(
+        f"bedGraph file {fname} requires a tabix index (.tbi or .csi). "
+        f"Create one with: tabix -p bed {fname}"
+    )
 
 
 def bedgraph_to_basecount(
@@ -90,14 +92,15 @@ def bedgraph_to_basecount(
             f"pysam is required for reading bedGraph files. {PYSAM_INSTALL_MSG}"
         ) from None
 
-    validate_bedgraph_index(bedgraph_fname)
+    index_fname = find_bedgraph_index(bedgraph_fname)
 
     # Read the BED file to get regions
     bed_data = tabio.read(bed_fname, "bed")
 
     # Open tabix file for random access with guaranteed cleanup
     try:
-        tbx = pysam.TabixFile(bedgraph_fname)
+        # pysam looks only for a .tbi index unless the index is named
+        tbx = pysam.TabixFile(bedgraph_fname, index=index_fname)
     except OSError as exc:
         raise ValueError(
             f"Failed to open bedGraph file {bedgraph_fname!r}. Error: {exc}"
@@ -196,14 +199,15 @@ def do_coverage(
     fasta: str | None = None,
     no_overlap: bool = False,
 ) -> CNA:
-    """Calculate coverage in the given regions from BAM read depths.
+    """Calculate coverage in the given regions from BAM read depths or a bedGraph.
 
     Parameters
     ----------
     bed_fname : str
         Path to BED file defining regions to measure coverage.
     bam_or_bg_fname : str
-        Path to BAM file containing aligned reads or bedGraph file (.bed.gz).
+        Path to a BAM or CRAM file of aligned reads, or a tabix-indexed bedGraph
+        file (.bed.gz) of per-base read depth.
     by_count : bool, optional
         Calculate coverage by read count instead of read depth. Default is False.
         Ignored for bedGraph input.