Bigwig Tracks

---
name: bio-genome-intervals-bigwig-tracks
description: Reads, queries, and writes bigWig indexed binary signal tracks (coverage, fold-change, conservation, methylation-rate) with pyBigWig (Python) and the UCSC Kent tools (bedGraphToBigWig, bigWigToBedGraph, bigWigInfo, bigWigSummary, bigWigAverageOverBed) and deepTools (multiBigwigSummary, computeMatrix, bigwigCompare). Covers the central trap that a wide query returns a precomputed zoom-level summary (by default the mean, which annihilates narrow peaks) not per-base data, when exact=True/values() is mandatory, the NaN-not-zero gap-handling fork, choosing mean vs max vs sum vs coverage by biological question, and the sorted-bedGraph plus chrom.sizes build requirement. Use when extracting signal at regions, computing mean signal per gene/peak, building a browser track from bedGraph, comparing tracks, or building TSS/gene-body metaprofiles.
tool_type: mixed
primary_tool: pyBigWig
---

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

| Tool | Role | Mechanism | When |
|------|------|-----------|------|
| pyBigWig | Python read/write | C-extension over libBigWig; `stats`/`values`/`intervals`/`addEntries` | inside a Python pipeline; custom per-region extraction; writing a bigWig |
| bigWigAverageOverBed | mean signal per BED region | one row per feature; `name,size,covered,sum,mean0,mean` | the right tool for "average signal per gene/peak"; gives both mean0 and mean |
| bigWigSummary | region -> N equal bins | reads zoom levels (like `exact=False`); `-type=mean/min/max/std/coverage` | quick binned profile at the command line |
| bigWigInfo | header/stats sanity check | version, zoom-level count, basesCovered, min/max/mean without parsing data | first thing to run on an unfamiliar file |
| bedGraphToBigWig / wigToBigWig | build bigWig | needs sorted input + chrom.sizes | converting a coverage bedGraph/WIG to a track |
| bigWigToBedGraph / bigWigToWig | bigWig -> text | `-chrom/-start/-end` for a sub-region | exact arithmetic; inspecting values as text |
| multiBigwigSummary | score matrix across many bigWigs | mean per bin over genome `bins` or a `BED-file` | track correlation/PCA (-> `plotCorrelation`/`plotPCA`) |
| computeMatrix | signal across many regions | `reference-point` (TSS/peak center) or `scale-regions` (gene body) | metaprofiles/heatmaps (-> `plotHeatmap`/`plotProfile`) |
| bigwigCompare | combine two bigWigs bin-by-bin | `--operation log2/ratio/subtract/...` `--pseudocount` | a log2(IP/input) or (treat-control) track |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Mean signal per gene/peak (one number per BED row) | `bigWigAverageOverBed` | purpose-built; pick `mean` (covered-only) vs `mean0` (gaps as zero) deliberately |
| Peak height / "is there a binding event here?" | `stats(type='max')` or `bigWigSummary -type=max` | `mean` dilutes a narrow peak to background |
| Total signal over an exon/gene (an amount) | `stats(type='sum')` or the `sum` column | extensive quantity; do not use `mean` for a total |
| A number going into a table/threshold | `stats(..., exact=True)` or `values()` | the default `exact=False` reads zoom levels, not base data |
| Per-base values for plotting/analysis | `values(numpy=True)` | one number per base; `nan` for gaps -> `np.nanmean`, never `np.mean` |
| "Is this region even assayed/mappable?" | `stats(type='coverage')` | fraction with data; a different axis from magnitude |
| Compare many tracks (correlation/PCA) | `multiBigwigSummary` -> `plotCorrelation`/`plotPCA` | mean per bin; bin size matters |
| Metaprofile/heatmap over TSS or gene bodies | `computeMatrix reference-point`/`scale-regions` -> `plotHeatmap`/`plotProfile` | match anchored-point vs whole-body mode |
| A ratio/difference track | `bigwigCompare --operation log2 --pseudocount` | pseudocount only meaningful for log2/ratio |
| Build a normalized coverage track from a BAM | -> chip-seq/chipseq-visualization or atac-seq/footprinting (deepTools `bamCoverage`) | generation is upstream; library-size normalization lives there |
| Render the track in a browser figure | -> data-visualization/genome-tracks | pyGenomeTracks/IGV; zoom-out IS the summary trap made visual |
| Discrete features (peaks/genes), not signal | bigBed, not bigWig | continuous-vs-interval; one interval per base defeats the format |

## Inspect a File Before Trusting It

```bash
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.

```bash
# 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:

```python
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)

```python
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)

```python
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).

```bash
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:

```python
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)

```bash
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; `nan_to_num`/`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

| Threshold | Source | Rationale |
|-----------|--------|-----------|
| `exact=True` when a number enters a result | pyBigWig design | default `exact=False` reads zoom levels (up to ~16x coarser than data); fine for exploration only |
| Zoom ladder: smallest bin ~16x mean interval size, each level 4x the previous | Kent 2010 convention | the resolution at which a wide query is answered; `bigWigInfo -zooms` shows the actual levels |
| Index < ~1% of data; ~10x smaller than bedGraph | Kent 2010 | order-of-magnitude; exact ratio is data-dependent (sparse vs dense) |
| bedGraph must be sorted `-k1,1 -k2,2n`, non-overlapping | bedGraphToBigWig requirement | signal is a function (one value per base); unsorted/overlapping input errors out |
| computeMatrix flank `-b/-a` 2000-3000 bp at TSS | metaprofile convention | captures promoter-proximal signal; widen for distal features; state the value used |
| bin size (e.g. 10-50 bp) | resolution vs file size | finer bins preserve narrow features but enlarge the file; state the bin when reading values back |

## Common Errors

| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| Region reads flat but browser shows a peak | wide `mean` query diluted the peak | use `type='max'`, more bins, or zoom to feature resolution |
| Per-region numbers a reviewer cannot reproduce | `exact=False` zoom approximation | re-extract with `exact=True` / `bigWigAverageOverBed` |
| `np.mean` returns NaN | gaps in the track (NaN, not 0) | `np.nanmean`, or `np.nan_to_num` if gaps are biologically zero |
| `mean` and `mean0` differ a lot in `bigWigAverageOverBed` | track is sparsely covered | pick the column by biology (depth -> mean0; rate -> mean) |
| `bedGraphToBigWig` errors / drops intervals | unsorted input or chrom-name/length mismatch | `sort -k1,1 -k2,2n`; match chrom.sizes to the reference |
| IGV will not render zoom-out | bigWig built with `maxZooms=0` | rebuild with zoom levels (default 10) |
| `bigwigCompare` rejects `--ratio` | flag renamed | use `--operation log2` |

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