Bio applied coverage tracks

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

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.

# 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

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.

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

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