bio-genome-intervals-coverage-analysis

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, mosdepth 0.3+, samtools 1.19+, pybedtools 0.10+, numpy 1.26+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

`samtools depth` behaviour changed across versions: pre-1.13 capped depth at 8000 and truncated silently (`-d 0` = unlimited); 1.13+ rewrote the subcommand with NO cap and `-d/-m` deprecated/ignored. Always check `samtools --version` before trusting a max-depth number. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Coverage Analysis

**"Is my sequencing deep enough to answer the question?"** -> Measure depth as a distribution over positions, then report median, breadth at a depth threshold, and an evenness number -- never the mean alone.
- CLI: `mosdepth --by 500 prefix in.bam` (windowed depth + cumulative dist), `samtools coverage in.bam` (per-contig depth+breadth), `bedtools genomecov -ibam in.bam -bga` (bedGraph track)
- Python: `pybedtools.BedTool('in.bam').genome_coverage(bga=True)` (pybedtools); parse `prefix.mosdepth.global.dist.txt` for the breadth curve

## The Single Most Important Modern Insight -- Mean Depth Is a Budget, Not a Result; Report Breadth Off a Cumulative Curve

"30x WGS" describes what was paid for, not what was achieved. Coverage is a **distribution over positions**, and the mean is its worst summary: it is dragged **up** by a fat right tail (repeats, rDNA, mitochondria, PCR pileups, segmental dups) while staying **blind** to a hard left wall of zeros and near-zeros (GC-extreme exons, poorly-mappable regions, capture dropout). Two libraries with identical mean 30x can differ completely -- one even and callable everywhere, one spiky with 20% of the target uncallable. The mean hides both failures. Four load-bearing moves:

1. **Report MEDIAN, not mean.** The median is robust to the right tail. When mean/median exceeds ~1.1-1.2 the distribution is skewed and the mean is overstating typical depth -- that gap is a free evenness diagnostic.
2. **Report a BREADTH / cumulative-coverage curve.** "% of target >= 1x, >= 10x, >= 20x, >= 30x" is the honest summary, because adequacy is a breadth statement: a base that was not covered deeply enough is uncallable no matter how deep the rest of the genome is. mosdepth's `*.mosdepth.global.dist.txt` IS this curve. The killer question for any "mean = 30x" claim is "breadth at 20x?".
3. **Quantify EVENNESS** (CV, Fano factor, Picard fold-80, or Gini) -- an even 30x and a spiky 30x are different experiments, and a spiky library cannot be rescued by sequencing deeper (extra reads follow the same biased distribution; the holes stay holes). Fix the library (PCR-free, better capture, UMIs), not the lane count.
4. **Say WHAT WAS COUNTED.** A depth number is meaningless until the recipe is stated: duplicates dropped (only if MARKED first)? secondary/supplementary included? MAPQ filter? read span or fragment? mate-overlap corrected? per-base or per-region? The tools disagree on every one of these by default.

## Tool Taxonomy

| Tool | Counts what (defaults) | Per-base or region | When |
|------|------------------------|--------------------|------|
| mosdepth | corrects mate-overlap by default (off under `--fast-mode`/`-x`); `-Q` MAPQ filter; emits cumulative dist + summary | windowed (`--by`), per-region, or callable bins (`--quantize`) | the modern fast default for WGS/WES/targeted; gives the breadth curve directly |
| samtools coverage | per-reference summary (added 1.10); `coverage` column = breadth %, `meandepth` = depth | per-contig | quick "is this contig actually covered?" -- spots high-mean/low-breadth pileups |
| samtools depth | drops UNMAP/SECONDARY/QCFAIL/DUP by default; `-Q`/`-q` filters; `-s` de-double-counts overlap; CRAM needs `--reference` | per-base | exact per-base depth over small regions; watch the 8000-cap version trap |
| bedtools genomecov | counts READ coverage by default (double-counts mate overlap); `-pc` = fragment; `-split` for spliced | per-base / bedGraph / histogram | bedGraph tracks, genome-wide depth histogram |
| bedtools coverage | per-A-interval stats from B reads; `-a`/`-b` flipped at v2.24.0 | per-region (or `-d` per-base) | per-target counts/breadth/mean over a BED |
| Picard CollectHsMetrics | capture-kit QC | per-target panel | exome/panel uniformity: on-target %, fold-80, PCT_TARGET_BASES_20X |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| WGS / WES breadth + adequacy | `mosdepth --by` then parse `*.global.dist.txt` | emits the cumulative curve + median directly; fast |
| Quick per-contig depth & breadth glance | `samtools coverage` | one line/contig; `coverage` col = breadth, `meandepth` = depth |
| Exact per-base depth, small region | `samtools depth -a -r chr:from-to` | per-base; add `-s` for short-insert; check version for 8000 cap |
| bedGraph coverage TRACK for a browser | `bedtools genomecov -ibam -bga` (or `-bg`) | `-bga` marks zero-coverage gaps; convert to bigWig -> bigwig-tracks |
| Per-target counts/breadth/mean over a BED | `bedtools coverage -a targets.bed -b in.bam` | A = targets, B = reads (post-v2.24.0); `-mean` for mean depth |
| Callable-region BED (NO/LOW/CALLABLE/HIGH) | `mosdepth --quantize 0:1:4:150:` | lightweight CallableLoci replacement at scale |
| Target-capture uniformity QC | -> Picard CollectHsMetrics (fold-80, on-target %) | the capture-QC standard; off-target loss + bait unevenness |
| Spliced/RNA-seq depth | add `-split` (genomecov/coverage) | without it an intron (N CIGAR) is counted as covered |
| Short-insert VAF (amplicon/cfDNA) | correct mate-overlap: `samtools depth -s` / `genomecov -pc` / mosdepth default | naive per-base double-counts the overlap, corrupting VAFs |
| Normalized cross-sample track | -> chip-seq/chipseq-visualization (deepTools bamCoverage) | library-size correction (RPGC/CPM/BPM) for comparison |
| Pileup/variant evidence from BAM | -> alignment-files/pileup-generation | depth is upstream of per-call DP/AD |

## mosdepth -- The Modern Default

**Goal:** Get the median depth and the full breadth curve for a BAM in one fast pass.

**Approach:** Run mosdepth windowed (or whole-genome), then read the cumulative distribution file -- it already holds breadth at every depth threshold; no histogram integration needed.

```bash
mosdepth --by 500 -Q 20 sample in.bam     # --by 500 = 500 bp windows; -Q 20 = drop MAPQ<20 (repeat coverage collapses, intentionally)
# Outputs: sample.mosdepth.summary.txt (mean/min/max per chrom + total)
#          sample.mosdepth.global.dist.txt (cumulative: chrom, depth, proportion >= depth)
#          sample.regions.bed.gz (per-window mean depth)
```

The `*.global.dist.txt` rows are `chrom  depth  proportion_of_bases_at_least_this_depth` -- the breadth curve directly. Read median as the depth where proportion crosses 0.5. `--fast-mode`/`-x` is ~2x faster but SILENTLY disables mate-overlap correction -- fine for a rough WGS glance, wrong for VAF-sensitive short-insert data.

**Goal:** Emit a callable-region BED (NO_COVERAGE / LOW / CALLABLE / HIGH) without GATK3.

**Approach:** Use `--quantize` to bin depth and merge adjacent equal-bin runs into a compact BED.

```bash
mosdepth --quantize 0:1:4:150: callable in.bam   # bins: [0,1)=NO_COVERAGE, [1,4)=LOW, [4,150)=CALLABLE, [150,inf)=HIGH
# 4 = min callable depth (tune to caller); 150 = excessive-depth ceiling (flags rDNA/artifact pileups)
zcat callable.quantized.bed.gz | head
```

## bedtools genomecov -- Tracks and the Histogram Default

```bash
bedtools genomecov -ibam in.bam -bga > cov.bedGraph   # -bga = bedGraph INCLUDING zero-coverage runs; -bg omits zeros
bedtools genomecov -ibam in.bam -pc -bg > frag.bedGraph # -pc = FRAGMENT coverage (mate overlap counted once); default c