Bio Applied Coverage Tracks

---
name: bio-applied-coverage-tracks
description: "Generate normalized bigWig coverage tracks from BAM with deepTools bamCoverage/bamCompare (RPKM/CPM/RPGC), summarize with multiBamSummary, and plot TSS/region signal with computeMatrix + plotHeatmap/plotProfile; pyBigWig for programmatic access. Use when normalizing BAM to bigWig, computing ChIP/input log2 ratio tracks, making TSS metagene heatmaps, or querying bigWig values in Python."
tool_type: bash
primary_tool: deepTools
---

# Coverage Tracks with deepTools

## When to Use

- Converting a sorted, indexed BAM into a normalized signal track (bigWig) for IGV/UCSC browser display
- Normalizing for sequencing depth across samples with RPKM, CPM, or RPGC (1x genome coverage) before comparing tracks
- Computing a ChIP-vs-input or treatment-vs-control log2 ratio track for enrichment visualization
- Building TSS-centered or peak-centered heatmaps/profile plots to compare signal shape across samples
- Summarizing multi-sample bigWig/BAM coverage over bins or regions for correlation/PCA QC, or pulling per-base values programmatically with pyBigWig

## Version Compatibility

deepTools ≥3.5.6, pyBigWig ≥0.3.25, samtools ≥1.19 (BAM must be coordinate-sorted and `.bai`-indexed). deepTools requires Python ≥3.9 and numpy/scipy/matplotlib on the same environment.

## Prerequisites

```bash
pip install deeptools pyBigWig
# or: conda install -c bioconda deeptools pybigwig
samtools sort -o sample.sorted.bam sample.bam && samtools index sample.sorted.bam
```

Prior concepts: coordinate-sorted/indexed BAM (`bio-alignment-files-alignment-sorting`, `bio-alignment-files-alignment-indexing`), BED/bigWig interval basics (`bio-genome-intervals-bed-file-basics`, `bio-genome-intervals-bigwig-tracks`).

## Step 1: BAM to Normalized bigWig

**Goal:** turn a BAM into a depth-normalized bigWig for browser display or downstream matrix building.
**Approach:** `bamCoverage` with `--normalizeUsing` set to RPKM/CPM for general use, or RPGC (needs `--effectiveGenomeSize`) to target 1x genome coverage; `--extendReads` for ChIP/ATAC fragment-level signal.

```bash
# RPKM-normalized, 10-bp bins, extend single-end reads to fragment size
bamCoverage -b sample.sorted.bam -o sample.rpkm.bw \
    --normalizeUsing RPKM --binSize 10 --extendReads --numberOfProcessors 8

# RPGC (1x genome coverage) — requires effective genome size (hg38: 2913022398)
bamCoverage -b sample.sorted.bam -o sample.rpgc.bw \
    --normalizeUsing RPGC --effectiveGenomeSize 2913022398 \
    --binSize 10 --extendReads --ignoreDuplicates

# ChIP-vs-input log2 ratio track
bamCompare -b1 chip.sorted.bam -b2 input.sorted.bam -o chip_over_input_log2.bw \
    --operation log2 --scaleFactorsMethod readCount --binSize 25 --numberOfProcessors 8
```

```python
import subprocess
from typing import Optional


def run_bam_coverage(bam: str, out_bw: str, normalize: str = "RPKM",
                      bin_size: int = 10, effective_genome_size: Optional[int] = None) -> str:
    """Wrap deepTools bamCoverage to produce a normalized bigWig.

    normalize: one of 'RPKM', 'CPM', 'BPM', 'RPGC', 'None'.
    RPGC requires effective_genome_size (e.g. 2913022398 for hg38).
    Assumes bam is coordinate-sorted with a .bai index alongside it.
    """
    if normalize == "RPGC" and not effective_genome_size:
        raise ValueError("RPGC normalization requires effective_genome_size")
    cmd = ["bamCoverage", "-b", bam, "-o", out_bw,
           "--normalizeUsing", normalize, "--binSize", str(bin_size),
           "--numberOfProcessors", "4"]
    if effective_genome_size:
        cmd += ["--effectiveGenomeSize", str(effective_genome_size)]
    subprocess.run(cmd, check=True)
    return out_bw


if __name__ == "__main__":
    # demo: verify command construction without requiring a real BAM
    import unittest.mock as mock
    with mock.patch("subprocess.run") as m:
        run_bam_coverage("sample.bam", "sample.bw", normalize="CPM")
        assert "--normalizeUsing" in m.call_args[0][0] and "CPM" in m.call_args[0][0]
        try:
            run_bam_coverage("sample.bam", "sample.bw", normalize="RPGC")
            raise AssertionError("expected ValueError for missing genome size")
        except ValueError:
            pass
    print("bamCoverage command construction OK")
```

## Step 2: Multi-Sample Summary and TSS/Region Matrices

**Goal:** QC many samples together (correlation/PCA over genome-wide bins) and visualize average signal shape around TSS or a BED of regions.
**Approach:** `multiBamSummary bins` (or `BED-file` mode) for correlation matrices; `computeMatrix reference-point`/`scale-regions` feeding `plotHeatmap`/`plotProfile`.

```bash
# Genome-wide 10kb-bin correlation matrix across samples (QC: replicate concordance)
multiBamSummary bins --bamfiles sampleA.bam sampleB.bam input.bam \
    --binSize 10000 --labels sampleA sampleB input -o readCounts.npz \
    --numberOfProcessors 8
plotCorrelation -in readCounts.npz --corMethod spearman --whatToPlot heatmap \
    -o figures/sample_correlation.png --plotNumbers

# TSS-centered matrix (+/-3kb) over a gene BED, then heatmap + profile
computeMatrix reference-point -S sample.rpkm.bw input.rpkm.bw \
    -R genes_hg38.bed --referencePoint TSS -b 3000 -a 3000 --binSize 10 \
    --skipZeros --numberOfProcessors 8 -o tss_matrix.gz \
    --outFileSortedRegions tss_regions_sorted.bed

plotHeatmap -m tss_matrix.gz -out figures/tss_heatmap.png \
    --colorMap RdBu_r --whatToShow 'heatmap and colorbar' \
    --samplesLabel ChIP Input --zMin 0 --zMax 5
plotProfile -m tss_matrix.gz -out figures/tss_profile.png \
    --samplesLabel ChIP Input --perGroup --plotTitle 'Signal +/-3kb around TSS'
```

## Step 3: Programmatic Access with pyBigWig

**Goal:** pull per-base or per-bin values, and header stats, directly in Python without going through deepTools plotting.
**Approach:** `pyBigWig.open()` gives header/chrom info; `.values()`/`.stats()` extract signal over an interval, optionally binned with `nBins`.

```python
import numpy as np
import pyBigWig


def mean_signal_over_region(bw_path: str, chrom: str, start: int, end: int, n_bins: int = 50) -> np.ndarray:
    """Return a fixed-length binned mean-signal profile over one region from a bigWig.

    Coordinates are 0-based half-open (BED convention), matching pyBigWig's API.
    """
    bw = pyBigWig.open(bw_path)
    try:
        if chrom not in bw.chroms():
            raise KeyError(f"{chrom} not found in {bw_path}")
        vals = bw.stats(chrom, start, end, type="mean", nBins=n_bins)
        return np.array([v if v is not None else 0.0 for v in vals])
    finally:
        bw.close()


def genome_wide_mean(bw_path: str) -> float:
    """Whole-genome mean signal, useful as a quick sanity check on normalization."""
    bw = pyBigWig.open(bw_path)
    try:
        return bw.header()["sumData"] / bw.header()["nBasesCovered"]
    finally:
        bw.close()


if __name__ == "__main__":
    # demo: create a tiny bigWig and verify round-trip reads
    tmp_bw = pyBigWig.open("/tmp/_demo.bw", "w")
    tmp_bw.addHeader([("chr1", 1000)])
    tmp_bw.addEntries("chr1", [0, 100, 500], values=[1.0, 2.0, 3.0], span=100, step=100)
    tmp_bw.close()

    profile = mean_signal_over_region("/tmp/_demo.bw", "chr1", 0, 600, n_bins=6)
    assert profile.shape == (6,)
    assert genome_wide_mean("/tmp/_demo.bw") > 0
    print("pyBigWig round-trip OK:", profile)
```

## Pitfalls

- **Normalization mismatch across samples**: comparing an RPKM track to a CPM or raw-coverage track invisibly biases fold-change interpretation — pick one method per comparison and apply it to every sample
- **Missing `--effectiveGenomeSize`**: RPGC silently uses the wrong scale factor if the genome size doesn't match the actual reference build (hg38 ≠ hg19 ≠ mm10 effective sizes)
- **Fragment vs. read signal**: forgetting `--extendReads` on single-end ChIP/ATAC data plots raw read length instead of fragment-length coverage, understating peak width
- **BAM not indexed/sorted**: `bamCoverage`/`multiBamSummary` fail or silently misbehave on coordinate-unsorted or unindexed BAM — always `samtools sort` + `samtools index` first
- **`computeMatrix` region file ordering**: `--outFileSortedRegions` reorders regions by clustering/sorting choice in `plotHeatmap`, so downstream region-to-signal joins must use that sorted BED, not the original
- **bigWig coordinates are 0-based half-open**: mixing this with 1-based GTF/VCF coordinates when building region BEDs causes off-by-one shifts in TSS-centered plots

## See Also

- `bio-applied-chipseq-pipeline` — full ChIP-seq pipeline including this deepTools track-generation step
- `chipseq-epigenomics` — peak calling, differential binding, and annotation upstream of coverage tracks
- `atac-seq-analysis` — Tn5-corrected fragment QC and footprinting on ATAC BAMs before generating tracks
- `bio-alignment-files-pileup-generation` — alternative per-base coverage via samtools/pysam pileups