bio-genome-intervals-bedgraph-handling

What this skill does

## Version Compatibility

Reference examples tested with: deeptools 3.5+, bedtools 2.31+, ucsc-bedgraphtobigwig 445+, pyBigWig 0.3.22+.

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

bedGraphToBigWig has a hard, under-advertised input contract: the bedGraph must be `LC_COLLATE=C`-sorted by chrom then start, contain non-overlapping intervals, and ship with a chrom.sizes derived from the exact assembly the reads were aligned to. deepTools effective-genome-size tables are occasionally updated between releases - re-check the installed version's table. If code throws an error, introspect the installed tool and adapt the example to match the actual API rather than retrying.

# bedGraph Handling

**"Make me a coverage/signal track I can compare across samples and load in a browser"** -> Generate a per-bin signal track, normalize it onto a common scale (or decide a spike-in is required), then convert the text bedGraph to an indexed bigWig under the strict sort/overlap/chrom.sizes contract.
- CLI: `bamCoverage -b s.bam -o s.bw --normalizeUsing RPGC --effectiveGenomeSize <N>`; `bedtools genomecov -ibam s.bam -bga`; `LC_COLLATE=C sort -k1,1 -k2,2n in.bdg | bedGraphToBigWig /dev/stdin chrom.sizes out.bw`
- Python: `pyBigWig.open('s.bw')` to read/extract; `bw.intervals(chrom, start, end)` returns the bedGraph rows

## The Single Most Important Modern Insight -- A Raw Coverage bedGraph Is a Library-Size Artifact, and the Wrong Normalization Is Worse Than None

Column 4 of a raw coverage bedGraph is not biology - it is sequencing depth. Two libraries of *identical* biology sequenced to different depths produce different heights, so any cross-sample statement ("more signal at this promoter in treatment") on un-normalized tracks is a category error. The modern path skips the text intermediate entirely: **deepTools bamCoverage takes BAM -> normalized bigWig in one step**, because bigWig is indexed, binary, random-access and bedGraph is flat text. Three load-bearing moves:

1. **Every library-size normalization (CPM/RPKM/BPM/RPGC=1x) assumes total signal is conserved across samples.** They all just rescale each library to a common total (per-million reads, or to 1x genome coverage). That model is correct when signal only *redistributes* locally - the usual case - and **actively wrong** when the perturbation changes *global* levels (histone-mark KD, BET-bromodomain inhibitor, global pol-II collapse). A genuine 3-fold global increase becomes, after CPM/RPGC, *no change* - the extra signal is spread thin and rescaled away. The model is unfalsifiable from the normalized data: forcing both libraries to the same total defines away any global difference. **Library-size normalization assumes the very thing under measurement does not happen.**
2. **There is no computational rescue for a global change after the fact.** The only fix is an external ruler decided AT THE BENCH - a spike-in of fixed foreign chromatin per cell (ChIP-Rx, Orlando 2014; defined reference epigenome, Bonhoure 2014) - scaled by the spike-in reads, not the sample reads. The wet-lab decision had to be made before sequencing; with no spike-in, the global scale is unrecoverable. The mechanics live in chip-seq/spike-in-normalization; the *decision* (could this perturbation change global levels?) belongs here, up front.
3. **bedGraph is scratch; bigWig is the artifact.** The text bedGraph is the last human-readable checkpoint - `awk '$4 > 1000'` to find blacklist pileups, confirm the sort/overlap invariants - before opaque binary. Inspect it, then ship bigWig. Never distribute a bedGraph as a final product: it is unindexed, so a browser reads the whole file to render any region.

## Normalization Taxonomy

| Method | What it assumes | When to use | When WRONG |
|--------|-----------------|-------------|------------|
| None | nothing (raw counts) | single-sample inspection only | any cross-sample comparison - depth confounds it |
| CPM | total mapped reads is the right denominator; total signal conserved | depth-only normalization; quick cross-sample on a common assay | a few high-coverage bins dominate (composition skew); global change |
| RPKM | as CPM plus bin length matters; total signal conserved | legacy default; depth + bin-length normalized | composition skew; global change; superseded by BPM for tracks |
| BPM (TPM-analog) | sum over all bins fixed at 1e6; total signal conserved | composition-aware cross-sample default; robust to a few dominant bins | global change (still a conserved-total rescale) |
| RPGC (1x) | mean genome-wide coverage = 1x; correct effective-genome-size; total signal conserved | field-standard ChIP/ATAC browser viewing; most interpretable height | wrong effective-genome-size (linear scaling error); global change |
| spike-in (external) | spike-in amount is constant per cell (a ruler that does not move) | global-level change plausible or under test | nothing computational - requires a bench step before sequencing |

All five library-size methods share one axiom: **total signal is conserved**. The decision is not which library-size method, it is whether library-size normalization is legitimate at all (see Decision Tree).

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| One BAM -> browser track, local redistribution | `bamCoverage --normalizeUsing RPGC --effectiveGenomeSize <N>` | one-step BAM->normalized bigWig; RPGC is the interpretable ChIP/ATAC standard |
| Cross-sample, composition skew likely | `bamCoverage --normalizeUsing BPM` | bins-per-million fixes the per-bin sum; robust to dominant bins |
| Global-level change plausible (KD/KO of a chromatin modifier, BET inhibitor) | spike-in -> chip-seq/spike-in-normalization | library-size normalization erases the global change by construction |
| RNA-seq coverage track | `bamCoverage --filterRNAstrand` or `genomecov -bga -split` | `-split`/strand handling so spliced reads do not paint introns |
| ChIP/ATAC track | add `--extendReads` (and `--centerReads` for footprints) | a read is a fragment END; raw read-end coverage is double-humped and wrong |
| Treatment vs input from raw BAMs | `bamCompare -b1 chip.bam -b2 input.bam --operation log2` | normalizes depth THEN does the arithmetic |
| Two already-normalized bigWigs | `bigwigCompare --operation log2` | arithmetic only - feeding un-normalized tracks manufactures a fake change |
| Stack N samples into a value matrix | `bedtools unionbedg -header -names ...` | union interval partition; feed the matrix to R/Python for testing |
| Sample-relatedness QC | `multiBigwigSummary bins` -> `plotCorrelation`/`plotPCA` | genome-wide value matrix for correlation/PCA |
| Need exact per-base arithmetic (not a browser) | keep bedGraph (`genomecov -bga`) | bedGraph is exact text; bigWig is binned/lossy |
| Convert finished bedGraph -> bigWig | `LC_COLLATE=C sort` then `bedGraphToBigWig` + matched chrom.sizes | the strict contract; inspect the text first |

## Generate a Normalized Track with bamCoverage (the modern default)

**Goal:** Turn one BAM into a normalized, browser-ready bigWig in a single command.

**Approach:** Let bamCoverage bin, normalize, and write bigWig directly; pick the normalization from the taxonomy, supply the effective-genome-size for RPGC, extend reads for ChIP/ATAC, and exclude chrX/chrM (and any spike-in contigs) from the scale-factor calculation.

```bash
BIN_SIZE=25                  # bp; smaller = finer + noisier + bigger. Match to feature width (sharp TF/ATAC 10-25; broad marks 50-200)
EFFGENOME=2913022398         # GRCh38 non-N length (faCount); use ONLY if multimappers were kept (see Effective Genome Size)

bamCoverage -b sample.bam -o sample.bw \
  --binSize $BIN_SIZE --normalizeUsing RPGC --effectiveGenomeSize $EFFGENOME \
  --extendReads --i