bio-hi-c-analysis-hic-visualization

What this skill does

## Version Compatibility

Reference examples tested with: cooler 0.10+, cooltools 0.7+, matplotlib 3.8+, bioframe 0.7+, HiCExplorer 3.7+, pyGenomeTracks 3.9+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

Notes specific to this skill: `.mcool` is multi-resolution -- pass a single-resolution URI (`file.mcool::/resolutions/10000`), never the bare path. A cooler must be balanced (`cooler balance`) before `matrix(balance=True)` returns anything but NaN. `cooltools.pileup` returns a stack of shape `(n_features, D, D)` -- aggregate over `axis=0`. cooltools standardised on `view_df`/`expected_df` arguments around 0.7+; verify with `help(cooltools.pileup)` before chaining.

# Hi-C Visualization

**"Plot my Hi-C contact matrix"** -> Choose a transform (raw / ICE-balanced / observed-over-expected), a matched colormap+norm (LogNorm for counts, symmetric-diverging for O/E and ratios), and a resolution that fits the feature; render as a square map or a rotated triangle for track-stacking, with NaN bins shown explicitly.
- Python: `clr.matrix(balance=True).fetch(region)` then `ax.matshow(m, norm=LogNorm(...), cmap='fall')`
- CLI: `hicPlotMatrix --matrix m.cool --region chr1:50-60Mb --log1p --colorMap fall -o out.png`

## The Single Most Important Modern Insight -- The Colorscale Is Where Hi-C Figures Lie

The same matrix under raw / ICE-balanced / observed-over-expected tells three *different biological stories*, and the choice is not cosmetic -- it decides which biology is legible. A balanced map is still dominated by the polymer distance-decay gradient (the bright diagonal falling off as ~P(s)); it shows TADs but washes out compartments and loops. Dividing by the distance-matched expected and taking `log2(O/E)` with a SYMMETRIC diverging cmap (`coolwarm`/`RdBu_r`, `vmin=-vmax`) removes that gradient and makes the compartment checkerboard and loop corner-dots suddenly visible -- they were always in the data. The corollary is a reviewer's reflex: a "no compartments / no loops" claim plotted on a balanced (not O/E) map is unsupportable. A reviewer-grade figure is one that can be reconstructed from the legend -- it states (1) the normalization (raw / ICE-balanced / O/E / log2-ratio), (2) the color scale (LogNorm vs symmetric-diverging) with its limits or clip percentile, and (3) the bin resolution. If those three are absent, the figure is neither interpretable nor reproducible.

## Transform Taxonomy

| Transform | Norm + cmap | What it shows | When |
|-----------|-------------|---------------|------|
| Raw counts | `LogNorm`, sequential (`fall`) | depth + per-bin coverage bias; white stripes are artifacts | QC sanity check only -- almost never the science figure |
| ICE-balanced | `LogNorm` vmin~1e-4..1e-1, `fall` | TADs + the distance-decay gradient; loops/compartments washed out | the honest "raw structure" map; track-stacking context |
| Observed/Expected | `log2`, symmetric `coolwarm`/`RdBu_r`, `vmin=-vmax` | compartment checkerboard + loop corner-dots | compartments, loops, any focal-enrichment claim |
| log2(cond1/cond2) | symmetric `RdBu_r`, `vmin=-vmax`, white=0 | gained/lost contacts | two conditions (balance + depth-match FIRST) |

## Layout Taxonomy

| Layout | Tool | Mechanism | When |
|--------|------|-----------|------|
| Square map | matplotlib `matshow`/`pcolormesh`, FAN-C `HicPlot2D` | symmetric 2D matrix | matrix itself is the result; inter-region rectangle; difference map |
| Rotated triangle | pyGenomeTracks/HiCExplorer `hic_matrix`, FAN-C `HicPlot`, plotgardener `plotHicTriangle`, CoolBox `style='triangular'` | 45deg shear, keep upper half, diagonal on top | STACKING genome-browser tracks below on a shared x-axis |
| Pileup (APA/saddle/on-diagonal) | `cooltools.pileup`/`saddle`, coolpup.py | average snippets over a feature set | the only honest genome-wide claim from sparse data |
| Virtual 4C | one matrix row, FAN-C `HicSlicePlot` | 1D profile from a viewpoint bin | compare against a real 4C anchor |
| Interactive | HiGlass | multires `.mcool` pan/zoom | exploration -- NOT a reproducible figure |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Show compartments / loops | `log2(O/E)`, symmetric `coolwarm`, `vmin=-vmax` | balanced map's gradient hides them |
| Show TADs / domains | balanced `LogNorm`, `fall`, 10-40kb | domain insulation lives at sub-Mb scale |
| Matrix + genes + ChIP + insulation stack | -> data-visualization/genome-tracks (pyGenomeTracks `hic_matrix`) | config = reproducible provenance; library does the shear |
| Quantify compartment strength | saddle plot (phase E1 first) -> compartment-analysis | corners give the single strength number |
| Validate a loop SET genome-wide | APA pileup, `log2(O/E)`, symmetric | one loop is invisible; 10k averaged is solid |
| Compare two conditions | side-by-side same-scale OR log2-ratio | balance + depth-match both FIRST |
| Export eigenvector / insulation as a track | -> genome-intervals/bigwig-tracks | bigWig feeds the track stack |
| Explore to find a region/resolution | HiGlass, then reproduce in a scripted tool | interactive != publication |

## Square Contact Map (Balanced, Log-Scaled, NaN Shown)

**Goal:** Render a balanced cis matrix for one region with honest dynamic range and visible masked bins.

**Approach:** Fetch the balanced matrix at a single-resolution URI, set `vmax` from a high off-diagonal percentile (report it), use `LogNorm`, and explicitly color NaN bins with `set_bad` so masked regions read as gray rather than as "zero contact".

```python
import cooler, numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import cooltools.lib.plotting   # registers the 'fall' cmap; needs matplotlib < 3.9 with cooltools 0.7.x (else use a stock cmap like 'afmhot_r')

clr = cooler.Cooler('matrix.mcool::/resolutions/10000')
region = ('chr1', 50_000_000, 60_000_000)
m = clr.matrix(balance=True).fetch(region)

vmax = np.nanpercentile(m[m > 0], 99.5)   # report this percentile in the legend
cmap = plt.get_cmap('fall').copy(); cmap.set_bad('lightgray')   # NaN bins shown, not white
fig, ax = plt.subplots(figsize=(7, 7))
im = ax.matshow(m, norm=LogNorm(vmin=vmax * 1e-3, vmax=vmax), cmap=cmap)
fig.colorbar(im, ax=ax, fraction=0.046, label='balanced (ICE)')
```

## Observed/Expected Divergent Map

**Goal:** Make the compartment checkerboard and loop corner-dots legible by removing the polymer distance-decay background.

**Approach:** Compute the cis expected with cooltools, fetch the matched-distance expected per pixel, divide observed by expected and take `log2`, then plot with a symmetric diverging cmap centered at 0 -- asymmetric limits move the white midpoint off zero and make neutral regions read as enriched.

```python
import cooltools, bioframe

view_df = bioframe.make_viewframe(clr.chromsizes)
expected = cooltools.expected_cis(clr, view_df=view_df, nproc=4)

chrom = region[0]
exp_by_diag = expected.query('region1 == @chrom')['balanced.avg'].to_numpy()   # expected per genomic separation
m = clr.matrix(balance=True).fetch(region)
i, j = np.indices(m.shape)
oe_mtx = m / exp_by_diag[np.abs(i - j)]   # divide each pixel by its distance-matched expected

v = 2.0   # symmetric clip; |log2(O/E)| up to ~2 is the usual readable range
fig, ax = plt.subplots(figsize=(7, 7))
im = ax.matshow(np.log2(oe_mtx), cmap='coolwarm', vmin=-v, vmax=v)   # vmin=-vmax mandatory
fig.colorbar(im, ax=ax, fraction=0.046, label='log2(obs/exp)')
```

`expected['balanced.avg']` is the per-diagonal expected; indexing it by `|i-j|` broadcasts it to a full per-pixel expected matrix.

## Rotated Triangle for