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.nantonum().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

Compare and Profile Tracks (deepTools)

bigwigCompare -b1 treat.bw -b2 control.bw -o log2ratio.bw --operation log2 --pseudocount 1   # NOT --ratio (older flag name)
multiBigwigSummary BED-file -b a.bw b.bw -o scores.npz --BED regions.bed                      # then plotCorrelation/plotPCA
computeMatrix reference-point -S signal.bw -R tss.bed -b 2000 -a 2000 -o matrix.gz            # anchored on TSS
plotHeatmap -m matrix.gz -o heatmap.png

Both bigwigCompare and multiBigwigSummary use mean-per-bin, so the zoom-level dilution caveat above applies; computeMatrix --missingDataAsZero is the same NaN-vs-zero fork inside deepTools.

Per-Method Failure Modes

Wide mean read as the peak

Trigger: bw.stats(chrom, start, end) (default type='mean') over a region wide relative to the feature. Mechanism: the mean dilutes a narrow tall peak toward background. Symptom: "no signal here" that a browser zoom-in contradicts. Fix: use type='max' (or nBins>1, or values()) for narrow features.

exact=False leaks zoom approximations into a result

Trigger: shipping default-exact stats() numbers into a table/threshold. Mechanism: the value is computed from the nearest zoom level, not base data, at up to 16x coarser granularity. Symptom: plausible numbers a reviewer cannot reproduce. Fix: pass exact=True (or use values()/bigWigAverageOverBed) whenever a number enters a result.

Averaging NaN as zero (or poisoning to NaN)

Trigger: np.mean(bw.values(...)) over a track with gaps, or mean0 on a rate track. Mechanism: np.mean poisons to NaN; nantonum/mean0/--missingDataAsZero averages real gaps as zeros. Symptom: a >3x swing or a NaN where a number was expected. Fix: decide biologically -- coverage track -> mean0/zero; rate/ratio track -> mean/np.nanmean.

addEntries before addHeader (or out of order)

Trigger: writing entries before the header, or in non-sorted order. Mechanism: the chrom list and offsets must exist and be ordered before data is appended. Symptom: runtime error or a corrupt file. Fix: addHeader([(chrom,length),...]) first, add entries in (chrom, start) order matching the header; close() to finalize.

chrom.sizes / naming mismatch on build

Trigger: bedGraphToBigWig with a chrom.sizes from a different assembly or naming (chr1 vs 1). Mechanism: the builder validates intervals against chrom lengths. Symptom: chromosome not found, or silently dropped/truncated intervals. Fix: derive chrom.sizes from the same reference (cut -f1,2 ref.fa.fai); harmonize naming; sort first.

Forcing peaks into a bigWig (or dense signal into bigBed)

Trigger: storing called peaks as a bigWig. Mechanism: bigWig is continuous signal; discrete features with per-feature metadata belong in bigBed. Symptom: lost boundaries/names, or an enormous one-interval-per-base file. Fix: signal -> bigWig; intervals/features -> bigBed.

Quantitative Thresholds

Common Errors

References

• Kent WJ, Zweig AS, Barber G, Hinrichs AS, Karolchik D. 2010. BigWig and BigBed: enabling browsing of large distributed datasets. Bioinformatics 26:2204-2207.
• Ramirez F, Ryan DP, Gruning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dundar F, Manke T. 2016. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44:W160-W165.
• pyBigWig (Devon Ryan / deepTools project) - C-extension wrapping libBigWig; no journal paper, see https://github.com/deeptools/pyBigWig
• UCSC Kent utilities (bedGraphToBigWig, bigWigInfo, bigWigSummary, bigWigAverageOverBed) - https://github.com/ucscGenomeBrowser/kent; cite Kent 2010 for the format.

Related Skills

• bedgraph-handling - The text bedGraph this skill converts to/from, and exact-arithmetic alternative
• coverage-analysis - Generates the per-base depth/bedGraph that becomes a bigWig
• bed-file-basics - The region BED files passed to bigWigAverageOverBed/computeMatrix
• chip-seq/chipseq-visualization - Generates normalized tracks (bamCoverage) and renders computeMatrix metaprofiles
• atac-seq/footprinting - Consumes bigWig signal over motif sites
• data-visualization/genome-tracks - Renders the bigWig in a browser figure (where zoom-out is the summary trap made visual)