--- name: bio-hi-c-analysis-hic-visualization description: Renders Hi-C contact matrices honestly and reproducibly with matplotlib, cooltools, HiCExplorer, pyGenomeTracks, FAN-C, CoolBox, and plotgardener. Covers the raw/ICE-balanced/observed-over-expected transform choice, LogNorm vs symmetric-diverging colormaps with vmax/percentile clipping, resolution-to-feature matching (compartments 100-500kb, TADs 10-40kb, loops 5-10kb), square vs rotated-triangle track-stacking, NaN/white-stripe handling, virtual 4C, APA/saddle/on-diagonal pileups, two-condition side-by-side and log2-ratio maps, and interactive (HiGlass) vs scripted-static publication figures. Use when plotting a contact matrix, choosing a normalization or color scale, building a multi-track Hi-C figure, making a virtual 4C profile, piling up loops/boundaries, or comparing two conditions. tool_type: python primary_tool: matplotlib --- ## 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 ` then `help(module.function)` to check signatures - CLI: ` --version` then ` --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 Track-Stacking **Goal:** Hang the contact map above aligned genome-browser tracks (genes, ChIP, insulation) on a shared x-axis. **Approach:** Prefer a library that owns the 45deg shear and the `depth` crop -- pyGenomeTracks/HiCExplorer (`file_type = hic_matrix`), FAN-C, plotgardener, or CoolBox -- because hand-rolling the `Affine2D` shear is where `extent`/`aspect` alignment bugs live. The matplotlib reference below is for a single panel; for a real stack, route to data-visualization/genome-tracks. ```bash # HiCExplorer / pyGenomeTracks: config-driven, reproducible. depth = how far up the diagonal. # A 2 Mb TAD needs depth >= ~2_000_000 or it is silently truncated. hicPlotTADs --tracks tracks.ini --region chr1:50000000-60000000 -o stack.png ``` ```python from matplotlib.transforms import Affine2D # matplotlib single-panel reference: shear the square map onto the diagonal. t = Affine2D().rotate_deg(45) + ax.transData im = ax.pcolormesh(np.log2(oe_mtx), cmap='coolwarm', vmin=-v, vmax=v) im.set_transform(t) ax.set_ylim(0, m.shape[0]) # crop the depth; the y-axis is genomic SEPARATION, not a 2nd coordinate ``` ## Virtual 4C from a Viewpoint **Goal:** Extract a 1D contact profile from one viewpoint bin to compare against a real 4C experiment. **Approach:** Take the viewpoint row from the balanced chromosome matrix; the near-cis distance-decay spike swamps distal signal on a linear axis, so plot on log-y (or mask the +/- few bins around the viewpoint) and say which. Cross-condition profiles must be balanced, depth-matched, and on identical y-axes. ```python res = clr.binsize vp_bin = (55_000_000 // res) - (50_000_000 // res) # viewpoint index within the region profile = clr.matrix(balance=True).fetch(region)[vp_bin, :] fig, ax = plt.subplots(figsize=(11, 2.5)) ax.semilogy(np.arange(len(profile)) * res / 1e6 + 50, profile) # log-y: the near-cis spike lies on linear ax.axvline(55, color='red', ls='--') ``` ## APA Pileup over a Loop Set **Goal:** Validate a loop call set genome-wide by averaging the contact signal centered on every loop's anchor pair. **Approach:** Pass BEDPE features and the cis expected to `cooltools.pileup` for an observed/expected stack, average over the feature axis (`axis=0`), and plot `log2` with a symmetric cmap; the center pixel is the loop, the APA score is center / a corner-background patch (Rao 2014 lower-left 3x3 convention). ```python import pandas as pd loops = pd.read_csv('loops.bedpe', sep='\t') # chrom1,start1,end1,chrom2,start2,end2 stack = cooltools.pileup(clr, loops, view_df=view_df, expected_df=expected, flank=100_000) apa = np.nanmean(stack, axis=0) # stack is (n_features, D, D) -> average over features c = apa.shape[0] // 2 apa_score = apa[c, c] / np.nanmean(apa[-3:, :3]) # center / lower-left 3x3 background fig, ax = plt.subplots(figsize=(5, 5)) im = ax.matshow(np.log2(apa), cmap='coolwarm', vmin=-1, vmax=1) ax.set_title(f'APA score {apa_score:.2f}') ``` For on-diagonal pileups over CTCF sites, strand-orient before averaging (`stack[mask] = stack[mask][:, ::-1, ::-1]` for `-` strand) or convergent/divergent signals cancel. For HiChIP/PLAC-seq anchored loops, route to loop-calling and the peak-anchored pileup conventions there. ## Two-Condition Comparison **Goal:** Show a contact change between two conditions without it being a depth/coverage artifact. **Approach:** Both matrices must be ICE-balanced AND depth-matched (downsample the deeper library to equal valid pairs) BEFORE ratioing. Then either side-by-side panels on an IDENTICAL cmap/norm/vmin/vmax/resolution, or a single `log2(cond1/cond2)` divergent map with symmetric limits and white = no change; grey out very-distal noise-amplified bins. ```python m1 = clr1.matrix(balance=True).fetch(region) m2 = clr2.matrix(balance=True).fetch(region) # clr1, clr2 already depth-equalized upstream ratio = np.log2((m1 + 1e-5) / (m2 + 1e-5)) # pseudocount tames divide-by-small off-diagonal fig, ax = plt.subplots(figsize=(7, 7)) im = ax.matshow(ratio, cmap='RdBu_r', vmin=-2, vmax=2) # symmetric, white=no change fig.colorbar(im, ax=ax, fraction=0.046, label='log2(cond1/cond2)') ``` Quantitative replicate-aware differential testing (not just a figure) lives in hic-differential. ## Per-Method Failure Modes ### Negative claim on a balanced map **Trigger:** "no compartments/loops" read off a balanced (not O/E) map. **Mechanism:** the polymer distance-decay gradient dominates balanced data and hides checkerboard/dots. **Symptom:** features absent that O/E would reveal. **Fix:** replot as `log2(O/E)` with a symmetric cmap before making any negative claim. ### Asymmetric divergent limits **Trigger:** `vmin != -vmax` on an O/E or log2-ratio map. **Mechanism:** zero (no change/enrichment) is no longer the white midpoint. **Symptom:** neutral regions read as enriched or depleted; reviewers flag it. **Fix:** `vmin=-v, vmax=v` (or `TwoSlopeNorm(vcenter=0)`). ### vmax = data max **Trigger:** letting the heavy-tailed diagonal set vmax. **Mechanism:** a few super-bins are orders of magnitude above the bulk. **Symptom:** the whole map looks empty/dark; over-clipping the other way fabricates structure. **Fix:** vmax = a stated high off-diagonal percentile (95th-99.5th). ### Resolution mismatched to feature **Trigger:** plotting at whatever the `.mcool` defaults to. **Mechanism:** loops at 100kb are averaged away (oversmoothing); compartments at 5-10kb mix in TAD/loop noise. **Symptom:** vanished dots or a noisy checkerboard. **Fix:** compartments 100-500kb, TADs 10-40kb, loops 5-10kb (Micro-C 1-2kb). ### Interpolated matrix quantified **Trigger:** `interp_nan`/`adaptive_coarsegrain` (or scHi-C smoothing) then measuring on the result. **Mechanism:** interpolation/imputation fabricates contacts for DISPLAY. **Symptom:** a filled centromere looks like contiguous chromatin; "structure" that is the smoother's prior. **Fix:** fill for display only; quantify on the raw/balanced matrix and disclose the smoother. ### Triangle depth too small **Trigger:** `depth` < the largest feature. **Mechanism:** the triangle crop truncates separations above `depth`. **Symptom:** a 2 Mb TAD silently cut off. **Fix:** set `depth >= ~feature size`; remember the y-axis is genomic separation. ### Unequal-depth ratio **Trigger:** `log2(cond1/cond2)` on un-depth-matched or unbalanced maps. **Mechanism:** a global depth difference is a uniform multiplicative offset. **Symptom:** a whole-map color shift read as biology. **Fix:** ICE-balance and downsample to equal valid pairs first. ## Quantitative Thresholds | Threshold | Source | Rationale | |-----------|--------|-----------| | Compartment resolution 100-500kb | compartment scale (Lieberman-Aiden 2009) | checkerboard is Mb-scale; finer bins add noise, not detail | | TAD resolution 10-40kb | domain scale (Dixon 2012) | insulation/boundary structure lives at sub-Mb | | Loop resolution 5-10kb (Micro-C 1-2kb) | focal-contact scale (Rao 2014) | a loop is a ~10kb focal pixel; coarse bins blur it, too-fine buries it in Poisson noise | | vmax = 95th-99.5th off-diagonal percentile | heavy-tailed counts | data-max vmax leaves the map dark; report the percentile | | Divergent limits symmetric `vmin=-vmax` | zero must be the midpoint | asymmetric limits misplace the white neutral point | | APA flank +/- 100kb | corner-background convention | too small contaminates the corner; too large averages in neighbors | | APA score = center / lower-left 3x3 | Rao 2014 | standard center-to-background loop enrichment ratio | | HiGlass zoom levels < ~5x apart | Kerpedjiev 2018 | adjacent resolutions must be close for smooth multires rendering | ## Common Errors | Error / symptom | Cause | Solution | |-----------------|-------|----------| | `matrix(balance=True)` all NaN | cooler not balanced | run `cooler balance` / `cooler.balance_cooler` first | | Empty / wrong-resolution result | bare `.mcool` passed | use `file.mcool::/resolutions/` | | White stripes mistaken for "no contact" | NaN bins left at default | `cmap.set_bad('lightgray')` to render masked bins | | `cmap='fall'` KeyError | colormap not registered | `import cooltools.lib.plotting` first | | `ImportError` on `import cooltools.lib.plotting` | matplotlib >= 3.9 dropped `register_cmap` (cooltools 0.7.x) | pin matplotlib < 3.9, or use a stock cmap (`'afmhot_r'`) | | Pileup looks averaged-out / scrambled | wrong nanmean axis | aggregate over `axis=0` (stack is `(n_features, D, D)`) | | Empty region / no overlap | chrom naming (`chr1` vs `1`) | harmonize names across cooler, BED/BEDPE, fasta | | `AttributeError` on cooltools fn | pre-0.7 vs 0.7+ API | `help(cooltools.)`; update to the `view_df`/`expected_df` signature | ## References - cooler: Abdennur N, Mirny LA. Cooler: scalable storage for Hi-C data and other genomically labeled arrays. *Bioinformatics* 2020;36(1):311-316. - cooltools: Open2C, Abdennur N, Abraham S, Fudenberg G, et al. Cooltools: enabling high-resolution Hi-C analysis in Python. *PLoS Comput Biol* 2024;20(5):e1012067. - HiCExplorer: Ramirez F, Bhardwaj V, Arrigoni L, et al. High-resolution TADs reveal DNA sequences underlying genome organization in flies. *Nat Commun* 2018;9(1):189. - pyGenomeTracks: Lopez-Delisle L, Rabbani L, Wolff J, et al. pyGenomeTracks: reproducible plots for multivariate genomic datasets. *Bioinformatics* 2021;37(3):422-423. - FAN-C: Kruse K, Hug CB, Vaquerizas JM. FAN-C: a feature-rich framework for the analysis and visualisation of chromosome conformation capture data. *Genome Biol* 2020;21(1):303. - CoolBox: Xu W, Zhong Q, Lin D, et al. CoolBox: a flexible toolkit for visual analysis of genomics data. *BMC Bioinformatics* 2021;22(1):489. - plotgardener: Kramer NE, Davis ES, Wenger CD, et al. Plotgardener: cultivating precise multi-panel figures in R. *Bioinformatics* 2022;38(7):2042-2045. - HiGlass: Kerpedjiev P, Abdennur N, Lekschas F, et al. HiGlass: web-based visual exploration and analysis of genome interaction maps. *Genome Biol* 2018;19(1):125. - coolpup.py: Flyamer IM, Illingworth RS, Bickmore WA. Coolpup.py: versatile pile-up analysis of Hi-C data. *Bioinformatics* 2020;36(10):2980-2985. - A/B compartments: Lieberman-Aiden E, van Berkum NL, Williams L, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. *Science* 2009;326(5950):289-293. - TAD directionality index: Dixon JR, Selvaraj S, Yue F, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. *Nature* 2012;485(7398):376-380. - Loops/HiCCUPS/APA: Rao SSP, Huntley MH, Durand NC, et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. *Cell* 2014;159(7):1665-1680. ## Related Skills - hic-data-io - Load the cooler/.mcool files and zoomify for multires HiGlass tilesets - matrix-operations - Balancing and O/E that the divergent map depends on - compartment-analysis - Eigenvector phasing behind the saddle plot - tad-detection - Insulation/boundary tracks stacked under the triangle - loop-calling - Loop calls and peak-anchored pileup conventions visualized here - hic-differential - Replicate-aware testing behind the two-condition comparison - data-visualization/genome-tracks - Config-driven multi-track stacks (pyGenomeTracks hic_matrix) - genome-intervals/bigwig-tracks - Export eigenvector/insulation as bigWig for the track stack - data-visualization/heatmaps-clustering - General heatmap color/norm conventions