Bio core hic analysis

Hi-C Analysis: 3D Genome Organization

When to Use

• Loading, inspecting, or balancing a Hi-C contact matrix stored as .cool/.mcool
• Visualizing raw or log-transformed contact matrices to spot TAD-like blocks
• Computing the P(s) contact-decay curve for Hi-C library QC
• Calling A/B compartments (active/inactive chromatin) via eigenvector decomposition of the O/E matrix
• Detecting TAD boundaries with the insulation score, or building aggregate pileup/APA plots over loop anchors

Version Compatibility

cooler ≥0.9, cooltools ≥0.6, numpy ≥1.24, pandas ≥2.0, matplotlib ≥3.7, Python ≥3.10.

Prerequisites

pip install cooler cooltools numpy pandas matplotlib

Assumes familiarity with genomic coordinates/bins and basic pandas. For raw read pair → .cool matrix generation, see bio-hi-c-analysis-hic-data-io; for matrix balancing details, see bio-hi-c-analysis-matrix-operations.

Background

DNA folds into a hierarchy of 3D structures: compartments (~1–10 Mb, A=active/B=inactive, found by eigenvector decomposition), TADs (~100 kb–3 Mb, found by insulation score), and loops (~10–300 kb, CTCF-anchored, found by pileup/APA). The Hi-C protocol cross-links cells, digests DNA with a restriction enzyme, proximity-ligates, and sequences — read pairs mapping far apart in the linear genome but close in 3D space appear as off-diagonal contacts. Contacts are binned into a matrix where entry (i, j) counts ligations between bins i and j; the cooler format stores this sparse matrix (plus optional multi-resolution zoom levels in .mcool) in HDF5.

Goal: Load a Hi-C contact matrix and inspect its resolution, chromosomes, and shape. Approach: Use cooler.Cooler() for a single-resolution .cool, or "file.mcool::resolutions/25000" for one resolution inside a multi-resolution .mcool. If no real file is available, generate a synthetic distance-decay matrix for a runnable demo.

import numpy as np
import pandas as pd
import cooler

def make_demo_cooler(path, n_bins=200, binsize=25_000, seed=42):
    """Create a minimal demo .cool file with synthetic distance-decay contacts."""
    chroms = pd.DataFrame({"name": ["chr1"], "length": [n_bins * binsize]})
    bins = pd.DataFrame({
        "chrom": ["chr1"] * n_bins,
        "start": np.arange(n_bins) * binsize,
        "end":   np.arange(1, n_bins + 1) * binsize,
    })
    rng = np.random.default_rng(seed)
    rows, cols, vals = [], [], []
    for i in range(n_bins):
        for j in range(i, min(i + 80, n_bins)):
            w = np.exp(-0.1 * (j - i)) * rng.poisson(20)
            if w > 0:
                rows.append(i); cols.append(j); vals.append(int(w))
    pixels = pd.DataFrame({"bin1_id": rows, "bin2_id": cols, "count": vals})
    cooler.create_cooler(path, bins=bins, pixels=pixels, dtypes={"count": np.int32})
    return path

COOL_FILE = "demo_hic_25kb.cool"
make_demo_cooler(COOL_FILE)

clr = cooler.Cooler(COOL_FILE)                 # for real data: cooler.Cooler("file.mcool::resolutions/25000")
print(f"Resolution: {clr.binsize:,} bp | Chromosomes: {clr.chromnames} | Shape: {clr.shape}")
mat = clr.matrix(balance=False).fetch("chr1:0-5000000").astype(float)  # balance=True once ICE weights exist

Goal: Assess Hi-C data quality with the contact-decay curve P(s) and detect A/B compartments. Approach: cooltools.expected_cis averages contacts by genomic distance to give the P(s) curve; cooltools.eigs_cis runs PCA on the observed/expected matrix, whose first eigenvector (E1) separates A (active, gene-dense) from B (inactive) chromatin. E1's sign is arbitrary — flip using GC content or gene density so positive = A.

import cooltools

def hic_view(clr):
    """Build the whole-genome view DataFrame cooltools needs for expected/eigs/insulation."""
    return pd.DataFrame({
        "chrom": clr.chromnames,
        "start": [0] * len(clr.chromnames),
        "end":   list(clr.chromsizes.values),
        "name":  clr.chromnames,
    })

view_df = hic_view(clr)

expected = cooltools.expected_cis(clr, view_df=view_df, ignore_diags=2)  # P(s): exclude near-diagonal artifacts
dist_bp = expected["dist"] * clr.binsize
count_col = expected.filter(like="avg").columns[0]

eigvals, eigvecs = cooltools.eigs_cis(clr, view_df=view_df, n_eigs=3, ignore_diags=2)  # A/B via O/E eigendecomposition
ev = eigvecs[eigvecs["chrom"] == "chr1"]
compartment = np.where(ev["E1"] > 0, "A", "B")     # sign is arbitrary — verify against GC/gene density
print(eigvecs[["chrom", "start", "end", "E1"]].head())

Goal: Call TAD boundaries from the insulation score. Approach: The insulation score at bin i averages contacts within a sliding square window centered on the diagonal at i; local minima are boundaries where cross-boundary contacts are depleted. Run cooltools.insulation with one or more window sizes (typically 100–500 kb) to probe TAD hierarchy at different scales.

def call_tad_boundaries(clr, view_df, window_bp):
    """Compute insulation score and boundary calls for one window size (bp)."""
    ins = cooltools.insulation(clr, window_bp=[window_bp], view_df=view_df, ignore_diags=2)
    score_col = f"log2_insulation_score_{window_bp}"
    boundary_col = f"is_boundary_{window_bp}"
    return ins, score_col, boundary_col

window = 10 * clr.binsize  # 10 bins ~ 250 kb at 25 kb resolution
insulation, score_col, boundary_col = call_tad_boundaries(clr, view_df, window)
boundaries = insulation[insulation.get(boundary_col, False)]
print(f"Boundaries found: {len(boundaries)} at window {window // 1000} kb")

Pitfalls

• ICE balancing: raw counts are biased by GC content, mappability, and fragment density. Use balance=True for compartment/insulation analysis, balance=False only for raw visualization — balanced weights must already exist in the cooler (from cooler balance) or these calls fail.
• Diagonal artifacts: the first few diagonals reflect unligated/self-ligated fragments, not real contacts — always pass ignore_diags=2 (or more) to expected_cis, eigs_cis, and insulation.
• E1 sign is arbitrary: don't assume positive = A compartment; verify against a GC-content or gene-density track and flip if needed.
• Resolution-dependent TAD calls: TADs appear at 25–40 kb; 5–10 kb shows sub-TADs, 100 kb shows compartment-scale domains. Always report the resolution and insulation window size used.
• .mcool vs .cool: .mcool holds multiple resolutions — you must select one with "file.mcool::resolutions/25000"; opening it bare raises an error.

See Also

• bio-hi-c-analysis-hic-data-io
• bio-hi-c-analysis-compartment-analysis
• bio-hi-c-analysis-tad-detection
• bio-hi-c-analysis-loop-calling