Signal heatmaps

Turning bigWig coverage into a region x track x bin signal cube, averaging replicates, and drawing grouped heatmaps of AR+F vs AR−F.

Table of contents

1. Peaks vs signal
2. The signal cube
3. Grouping replicates by averaging columns
4. Drawing the heatmap
5. Reading the result

---

Peaks vs signal

Peaks (last page) are discrete calls — where a factor binds. bigWig files hold the continuous read coverage genome-wide — how strong the signal is at every base. Heatmaps and browser tracks are built from bigWigs.

The signal cube

genomeblocks.signal.signal() (also Loci.signal) extracts binned coverage for a set of loci across a list of bigWigs:

import numpy as np

regions = ARpF + ARmF                       # rows: AR+F peaks then AR-F peaks
bw_order = ["ATAC_0h_r1", "ATAC_0h_r2", "ATAC_4h_r1", "ATAC_4h_r2",
            "AR_0h", "AR_4h", "FOXA1_0h", "FOXA1_4h"]
bigwigs = [BW[k] for k in bw_order]

S = np.nan_to_num(regions.signal(bigwigs, n_bins=200, flank=2000, workers=4))

The return value S is a NumPy array of shape (n_loci, n_tracks, n_bins):

• n_loci rows — one per row of regions, in its order: every AR+F peak followed by every AR−F peak.
• n_tracks — one slab per bigWig, in the order you passed them (a {name: path} dict works too, tracks in its order).
• n_bins columns — each window is split into n_bins equal bins.

Key arguments:

arg	meaning
n_bins=200	bins per region (heatmap column resolution)
flank=2000	window is the locus centre ± flank bp (set span=True to bin the whole interval instead)
workers=4	parallel extraction across processes into a shared-memory cube; the default 1 is one sequential native pass
backend=	the bigWig reader: pybigtools (default), pybigwig or python — same numbers from each

bigWigs can have gaps (no coverage); np.nan_to_num turns those NaNs into 0 so the heatmap and any averaging behave.

Grouping replicates by averaging columns

The two ATAC replicates are two adjacent track slabs. Averaging them is just a mean over the track axis — exactly what the browser does for a list of bigWigs. We collapse 8 tracks into 6 conditions:

S6 = np.stack([
    S[:, [0, 1], :].mean(1),   # ATAC 0h  (rep1 + rep2)
    S[:, [2, 3], :].mean(1),   # ATAC 4h  (rep1 + rep2)
    S[:, 4, :], S[:, 5, :],    # AR 0h, AR 4h
    S[:, 6, :], S[:, 7, :],    # FOXA1 0h, FOXA1 4h
], axis=1)
samples = ["ATAC 0h", "ATAC 4h", "AR 0h", "AR 4h", "FOXA1 0h", "FOXA1 4h"]

S[:, [0, 1], :].mean(1) averages the two ATAC-0h slabs over axis 1, leaving (n_loci, n_bins); np.stack(..., axis=1) reassembles the six condition slabs into a (n_loci, 6, n_bins) cube.

Drawing the heatmap

plot_heatmap (in genomeblocks.signal_draw, also Loci.plot_heatmap) renders one heatmap column per track, optionally split into row groups:

from genomeblocks.signal_draw import plot_heatmap

vmax = float(np.percentile(S6, 99))         # robust colour cap
fig = plot_heatmap(regions, S6,
                   groups={"AR+F": ARpF, "AR-F": ARmF},
                   sets=["AR+F", "AR-F"],
                   samples=samples,
                   vmax=vmax, ymax=vmax,
                   height=2000)                 # the flank used above: ticks read -2kb / +2kb

How the arguments shape the figure:

• groups — a plain dict {name: rows}. Each group becomes a block of rows. A group is a Loci (rows of regions with the same chromosome, start and end), a boolean mask over regions, or row numbers — so {"AR+F": np.arange(len(ARpF)), "AR-F": np.arange(len(ARpF), len(regions))} draws the same figure.
• sets — the order (top → bottom) the groups are drawn in.
• samples — column labels, one per track in S6.
• vmax / ymax — colour-scale cap and profile-track y-max (a scalar, or one per track). Using the 99th percentile keeps a few hot bins from washing out the rest.

To make tracks comparable across the figure, TMM-normalise the cube first with from genomeblocks import tmm; S = tmm(S) before grouping. The notebook leaves this commented out so the raw coverage is shown; flip it on when comparing libraries of different depth. genomeblocks.interop.cube_to_xarray(S6, regions, samples, flank=2000) labels the cube’s axes if you would rather slice it in xarray.

Reading the result

Rows are AR peaks, split into the FOXA1-dependent (AR+F) and FOXA1-independent (AR−F) blocks; columns are the six conditions. The interesting contrasts:

• ATAC 0 h vs 4 h — is the site already open before androgen, or does it open on stimulation? FOXA1-pioneered (AR+F) sites tend to be pre-accessible.
• AR 0 h vs 4 h — confirms the DHT-induced AR gain that defines the cistrome.
• FOXA1 0 h vs 4 h — by construction, strong in AR+F and weak in AR−F.

Next: Genomic annotation →