Bio Genome Intervals Bigwig Tracks

Version Compatibility

Reference examples tested with: pyBigWig 0.3.22+, numpy 1.26+, ucsc-bedgraphtobigwig/ucsc-tools 469+, deeptools 3.5+.

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version (or bigWigInfo with no args for usage) then <tool> --help to confirm flags
• Python: pip show pyBigWig then help(pyBigWig.bigWigFile.stats) to check signatures

Building any bigWig needs a chrom.sizes file (name<TAB>length) and a coordinate-sorted bedGraph; pyBigWig's numpy return path requires numpy present at compile time. If code throws an error, introspect the installed tool and adapt rather than retrying.

BigWig Tracks

"Get the signal from my bigWig over these regions / build a browser track." -> Query an indexed binary signal track, choosing the summary statistic and exactness that match the biological question, or build one from a sorted bedGraph + chrom.sizes.

• CLI: bigWigAverageOverBed in.bw regions.bed out.tab, bigWigSummary in.bw chr s e N -type=max, bedGraphToBigWig in.sorted.bedGraph chrom.sizes out.bw, bigWigInfo in.bw
• Python: bw=pyBigWig.open('x.bw'), bw.stats(chr,s,e,type='max',exact=True), bw.values(chr,s,e,numpy=True), bw.intervals(chr,s,e) (pyBigWig)

The Single Most Important Modern Insight -- A Wide Query Returns a Zoom-Level Summary, Not the Underlying Data

bigWig is fast (Kent 2010) because it stores, alongside base-resolution values, a ladder of precomputed zoom levels holding per-bin sum/sumSquared/min/max/nBasesCovered. A B+ tree resolves the chromosome, an R-tree (cirTree) finds the data blocks in O(log n), per-block zlib keeps it ~10x smaller than bedGraph, and the zoom ladder answers a wide region in near-constant time by reading a precomputed summary instead of the base data. That speed is bought with two stacked approximations, both ON by default, and a third trap at the moment the signal is reduced to a single number:

1. WHICH statistic. Over one wide bin mean (the default) dilutes a narrow tall feature toward background: a 200 bp ChIP summit of 500 in a 1 Mb sea of 1 averages to ~1.1 -- indistinguishable from background, while type='max' returns 500. Same file, same coordinates, opposite conclusions, decided by the named statistic. mean is faithful for broad features (domains, gene-body coverage) and a lie for narrow ones. max=peak height, sum=total amount (scales with width), coverage=fraction of bases with any data (ignores magnitude), std=variability.
2. WHERE the number comes from. exact=False (the pyBigWig default, and what bigWigSummary and every zoomed-out browser do) computes from the nearest zoom level, not base data. Fine for exploration and broad features; exact=True (or values()) is mandatory whenever a number enters a result -- a per-region average in a table, a threshold call, anything a reviewer recomputes. Plausible-but-zoom-approximated is the worst failure: it does not error, it rounds the biology.
3. NaN is NOT zero. Uncovered positions are no-data, surfaced as NaN in values() and as gaps between intervals() runs -- never 0. On a region 30% covered at signal 10: np.mean -> NaN (poisons), np.nanmean -> 10 (covered-only, = bigWigAverageOverBed mean column), gaps-as-zero (np.nan_to_num().mean(), deepTools --missingDataAsZero, the mean0 column) -> 3.0. A >3x swing in the headline number, and which is correct is biological: coverage/read-depth tracks -> gaps are zero (mean0); rate/ratio tracks (methylation %, log2FC, conservation) -> gaps are undefined (mean/nanmean).

Name the biological question first; the statistic, the exact flag, and the gap-handling then follow deterministically. Left on default, all three conspire to hand back a fast, confident, wrong answer.

Tool Taxonomy

Decision Tree by Scenario

Inspect a File Before Trusting It

bigWigInfo coverage.bw                  # version, # zoom levels, basesCovered, min/max/mean/std
bigWigInfo -chroms coverage.bw          # chrom names + lengths (the file carries its own chrom list)

A zero or low zoom-level count means a zoomed-out browser will read base data slowly (or, with maxZooms=0, IGV breaks). basesCovered far below the genome size means most positions are no-data (NaN), which makes the mean-vs-mean0 choice below load-bearing.

Mean Signal per Region (the most common task)

Goal: Compute one signal number per gene/peak, choosing covered-only vs gaps-as-zero by the track's biology.

Approach: Use the purpose-built bigWigAverageOverBed (BED needs a unique name column) and read the right output column -- mean (covered bases only) for rate/ratio tracks, mean0 (uncovered counted as zero) for coverage/depth tracks.

# BED4+ with a UNIQUE name in column 4; output columns: name size covered sum mean0 mean
bigWigAverageOverBed coverage.bw genes.bed signal_per_gene.tab
# -> read $6 (mean, covered-only) for methylation/log2FC; $5 (mean0, gaps=0) for read depth

The pyBigWig equivalent, when the extraction is inside a Python pipeline -- note exact=True because these numbers enter a result, and an explicit gap decision:

import pyBigWig
import numpy as np

bw = pyBigWig.open('coverage.bw')
GAPS_ARE_ZERO = False   # True for read-depth/coverage tracks; False for rate/ratio (methylation, log2FC)

def region_signal(chrom, start, end):
    v = bw.values(chrom, start, end, numpy=True)                         # per-base, nan for gaps
    if GAPS_ARE_ZERO:
        return float(np.nan_to_num(v).mean())                           # mean0: gaps counted as 0 (= bigWigAverageOverBed mean0)
    return np.nanmean(v) if not np.all(np.isnan(v)) else float('nan')    # covered-only (= bigWigAverageOverBed mean); stats(type='mean') is also covered-only, NOT mean0

Peak Height vs Total vs Coverage (statistic = question)

import pyBigWig
bw = pyBigWig.open('chip.bw')
region = ('chr1', 1_000_000, 2_000_000)

peak  = bw.stats(*region, type='max', exact=True)[0]        # binding-event height; mean would dilute it
total = bw.stats(*region, type='sum', exact=True)[0]        # total signal (amount; scales with width)
assayed = bw.stats(*region, type='coverage', exact=True)[0] # FRACTION of bases with any data (0..1), ignores magnitude
profile = bw.stats(*region, type='max', nBins=1000)         # 1000-bin max profile; nBins keeps narrow features visible

stats() returns a list of length nBins (default 1). type is one of mean(default)/max/min/coverage/std/sum. Use max with nBins>1 to see narrow features across a wide window; a single-bin mean over a megabase buries every peak.

Per-Base Values (NaN is not zero)

import pyBigWig
import numpy as np
bw = pyBigWig.open('coverage.bw')

v = bw.values('chr1', 1_000_000, 1_001_000, numpy=True)   # list by default; numpy=True -> ndarray, nan for gaps
covered_mean = np.nanmean(v)                               # ignores gaps (= bigWigAverageOverBed mean)
depth_mean = np.nan_to_num(v).mean()                       # gaps counted as zero (= mean0); only for depth tracks
raw = bw.intervals('chr1', 1_000_000, 1_001_000)          # [(start,end,value),...] the unresampled stored runs
bw.close()

Build a Valid bigWig

Goal: Turn a coverage bedGraph into an indexed, browser-ready bigWig.

Approach: Coordinate-sort the bedGraph, supply a chrom.sizes whose names and lengths match the bedGraph exactly, and run bedGraphToBigWig (which builds the index + zoom levels).

sort -k1,1 -k2,2n coverage.bedGraph > coverage.sorted.bedGraph   # bedGraphToBigWig REQUIRES sorted, non-overlapping input
cut -f1,2 reference.fa.fai > chrom.sizes                          # or fetchChromSizes hg38 > chrom.sizes
bedGraphToBigWig coverage.sorted.bedGraph chrom.sizes coverage.bw

Writing directly with pyBigWig -- addHeader (ordered chrom list) MUST precede addEntries, and entries must be added in sorted (chrom, start) order matching the header:

import pyBigWig
bw = pyBigWig.open('out.bw', 'w')
bw.addHeader([('chr1', 248956422)])   # ordered (name,length); maxZooms default 10; maxZooms=0 disables zoom and breaks IGV
bw.addEntries(['chr1'], [0], ends=[100], values=[1.5])           # mode (a) variable intervals
# mode (b) variableStep: bw.addEntries('chr1', [0,100], values=[1.5,2.3], span=20)
# mode (c) fixedStep:    bw.addEntries('chr1', 0, values=[1.5,2.3], span=20, step=30)
bw.close()                            # close() builds the R-tree index + zoom ladder