Bio Genome Intervals Bedgraph Handling

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

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

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.

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 --ignoreForNormalization chrX chrM -p 8

Defaults to verify: --binSize 50, --normalizeUsing None, --scaleFactor 1.0, --extendReads off, --centerReads off. For single-end ChIP supply the fragment length (--extendReads 200); paired-end infers it. --scaleFactor with --scaleFactorsMethod None is the hook for a bench-derived spike-in factor. --outFileFormat bedgraph writes the text form when the raw numbers are needed.

Generate with bedtools genomecov (text, flexible, no normalization)

-bg collapses equal-coverage runs but omits zero-coverage regions; -bga additionally tiles zeros (use when downstream tools need explicit 0s). -split is mandatory for RNA-seq so spliced reads do not paint introns. -scale 1000000/<nreads> is a crude manual RPM; deepTools is preferred for real normalization.

bedtools genomecov -ibam sample.bam -bga -split > sample.bedgraph

Convert bedGraph -> bigWig (the silent-corruption trap)

Goal: Produce a valid bigWig from a finished bedGraph without shipping a file that loads but lies.

Approach: C-locale-sort, guarantee non-overlapping intervals, derive chrom.sizes from the exact aligned-to FASTA, inspect the text, then convert.

samtools faidx ref.fa && cut -f1,2 ref.fa.fai > chrom.sizes   # chrom.sizes from the SAME FASTA the reads aligned to
LC_COLLATE=C sort -k1,1 -k2,2n sample.bedgraph > sample.sorted.bedgraph   # C locale: locale-aware sort triggers "is not case-sensitive sorted"
bedGraphToBigWig sample.sorted.bedgraph chrom.sizes sample.bw

If concatenation/merging introduced overlaps, collapse with an explicit aggregation BEFORE converting - and note max vs mean vs sum are different signals, there is no safe default:

bedtools merge -i sample.sorted.bedgraph -d 0 -c 4 -o max > sample.nonoverlap.bedgraph

bigWig round-trips losslessly: bigWigToBedGraph sample.bw out.bedgraph (optionally -chrom=chr1 -start=1000 -end=2000).

Multi-Sample Arithmetic

Goal: Compare two tracks (treatment/input, two conditions) without letting a depth difference masquerade as biology.

Approach: From raw BAMs use bamCompare, which normalizes depth THEN applies the operation; only use bigwigCompare on bigWigs that are already on a common scale.

bamCompare -b1 chip.bam -b2 input.bam -o log2ratio.bw \
  --operation log2 --pseudocount 1 --binSize 25 --scaleFactorsMethod readCount

--operation (NOT --ratio) chooses log2/ratio/subtract/add/mean/reciprocal_ratio/first/second; default log2. --scaleFactorsMethod readCount (the default) scales by library size; --scaleFactorsMethod SES (signal-extraction scaling, Diaz 2012) instead estimates the factor from the shared background bins and is more robust than readCount for SHARP/punctate marks and TF ChIP where enrichment is a small genomic fraction; it DEGRADES for broad marks (H3K27me3/H3K9me3) where the diffuse enrichment cannot be cleanly separated from background, so use readCount (or spike-in) there. --pseudocount (default 1) prevents divide-by-zero in log2/ratio but pulls low-coverage bins toward 0 - a log2 track's apparent dynamic range is partly a pseudocount+bin-size artifact, do not read fold-changes off a browser track as measured. bigwigCompare --skipZeroOverZero drops bins that are 0 in both rather than flooding the output with log2(1)=0. Stack many samples and QC relatedness:

bedtools unionbedg -i s1.bdg s2.bdg s3.bdg -header -names s1 s2 s3 > matrix.txt   # inputs must be coordinate-sorted
multiBigwigSummary bins -b s1.bw s2.bw s3.bw -o scores.npz && plotCorrelation -in scores.npz --corMethod spearman --whatToPlot heatmap -o corr.png

Read/Extract Signal with pyBigWig

import pyBigWig

bw = pyBigWig.open('sample.bw')
mean_over_region = bw.stats('chr1', 1_000_000, 1_010_000, type='mean')[0]   # binned summary, not per-base
rows = bw.intervals('chr1', 1_000_000, 1_010_000)   # the underlying bedGraph rows: (start, end, value)
bw.close()

bw.stats()/bw.values() return what the bin resolution preserved, not a faithful per-base record - coarse bins silently change the values read back.