Bio Hi C Analysis Loop Calling

Version Compatibility

Reference examples tested with: cooltools 0.7+, cooler 0.10+, bioframe 0.7+, chromosight 1.6+, mustache 1.3+

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.

The .cool must be BALANCED before calling loops -- dots/pileup read the weight column and raw counts are unsupported. An .mcool is multi-resolution; pass a single-resolution URI (file.mcool::/resolutions/10000), not the bare .mcool. cooltools changed signatures around 0.5 -> 0.7 (view_df/expected_value_col conventions); verify help(cooltools.dots) for the installed version. The view_df passed to expected_cis MUST be the same one passed to dots/pileup.

Chromatin Loop Calling

"Where are the focal loops (CTCF/cohesin corner-dots, E-P contacts) in my Hi-C map?" -> Test each off-diagonal pixel for focal enrichment against its local background (on a balanced, expected-normalized matrix), control FDR, then validate the set by aggregation and orthogonal support.

• Python: cooltools.dots(clr, expected=cooltools.expected_cis(clr, view_df=arms), view_df=arms)
• CLI: chromosight detect --pattern loops --min-dist 20000 --max-dist 2000000 sample.cool::/resolutions/5000 out

The Single Most Important Modern Insight -- Loop Calling Is Depth-Limited, Not Algorithm-Limited; the First Question Is "How Deep Is the Map?"

The caller choice is second-order. The dominant variable in whether loops are found at all is sequencing depth / map resolution. Rao 2014 needed ~4.9 BILLION contacts in GM12878 to reach 1kb bins and call ~10,000 loops; robust de-novo calling realistically wants 5-10kb resolution, which is hundreds of millions to billions of valid cis pairs. Below that, every caller returns near-nothing or noise, and tuning the FDR will not rescue it. So the workflow forks on depth before any tool is chosen:

1. Deep map (>=~500M-1B valid pairs, 5-10kb resolution): de-novo calling is licensed. Run cooltools dots (or chromosight / Mustache), then validate (see below).
2. Shallow map: do NOT de-novo call. Run APA / pileup on a KNOWN anchor set -- loops imported from a deep reference map, or anchor pairs built from CTCF/cohesin ChIP-seq peaks. This is the single most important practical reframe in the skill: shallow data can still confirm and quantify a hypothesized loop set even when it cannot discover one.

Two corollaries that follow directly:

De-novo calling DISCOVERS; APA CONFIRMS -- never conflate them. APA aggregates many putative loops to surface mean signal no individual loop could pass FDR for. An enriched APA center pixel proves "this SET of pairs is enriched on average"; it does NOT prove any single pair is a loop and it cannot discover new loops. Presenting an APA pileup as evidence that "these loops exist" is the classic abuse. And the APA score is meaningless without a corner control -- center pixel divided by an off-diagonal corner block of the flank is the on-vs-off measurement; the bare center value alone says nothing.

Loops form between convergent CTCF motifs -- biology AND a validation filter. Loops preferentially link two CTCF motifs in CONVERGENT orientation (Rao 2014 observation; de Wit 2015, Sanborn 2015 extrusion mechanism; proven by CTCF-site inversion experiments that kill or reroute the loop). A called corner-dot whose two anchors carry convergent CTCF motifs is high-confidence; one with no CTCF/anchor support on a shallow map is likely a false positive. Not all loops are CTCF loops (E-P and polycomb loops exist), so convergent-CTCF is a strong positive filter, not a universal requirement.

Loop-Caller Taxonomy

Forcato 2017 (Nat Methods 14:679) is the canonical finding that loop callers show LOW pairwise overlap and poor replicate reproducibility -- far worse than TAD callers. Practical consequence: a loop called by only one tool is suspect. Trust comes from consensus across >=2 callers plus orthogonal support (convergent CTCF, ChIA-PET/HiChIP), not from any single tool's list length.

Decision Tree by Scenario

De-Novo Loop Calling with cooltools dots

Goal: Discover focal loops genome-wide on a deep, balanced map with honest FDR control.

Approach: Build chromosome-arm regions, compute the distance-decay expected on those arms, then run dots -- which convolves four local-background kernels and runs Benjamini-Hochberg FDR independently within geometrically-spaced lambda-bins of locally-adjusted expected. The arms view_df must be identical for expected_cis and dots.

import cooler, cooltools, bioframe

clr = cooler.Cooler('matrix.mcool::/resolutions/10000')   # 10kb: a realistic de-novo floor; finer needs more depth
arms = bioframe.make_viewframe(clr.chromsizes)   # or cooltools.lib.read_viewframe_from_file('hg38_arms.bed', clr) for per-arm
expected = cooltools.expected_cis(clr, view_df=arms, nproc=4)   # distance-matched background; same view as dots
loops = cooltools.dots(
    clr, expected=expected, view_df=arms,
    max_loci_separation=10_000_000,   # ignore pixels farther than 10Mb from the diagonal
    n_lambda_bins=40, lambda_bin_fdr=0.1,   # FDR run independently per geometric lambda-bin (HiCCUPS default)
    clustering_radius=20_000,   # merge called pixels within 20kb into one loop
    nproc=4,
)

The four backgrounds, and why lower-left is the clever one. A pixel must beat ALL four local-background kernels, not one. Donut = is it a focal enrichment at all. Horizontal and vertical = is it actually a STRIPE pixel masquerading as a dot (these kernels exist to NOT call architectural stripes as loops). Lower-left = is it just a TAD/contact-domain CORNER -- a domain corner is enriched vs the donut but NOT vs its lower-left neighborhood, so requiring the pixel to also beat lower-left separates a genuine point loop from a generic domain corner. Skipping lower-left inflates calls with domain corners.

Lambda-chunking is why HiCCUPS FDR is honest. Contact counts span orders of magnitude with genomic distance, so a single genome-wide BH-FDR would be dominated by the high-count near-diagonal regime and over-call. dots bins pixels by their locally-adjusted expected into geometrically-spaced lambda-bins (n_lambda_bins=40) and runs BH-FDR independently within each (lambda_bin_fdr=0.1), so low-count and high-count regimes are each thresholded correctly.

Template-Matching with chromosight

# detect: <contact_map> <prefix> are positional and come LAST
chromosight detect --pattern loops --threads 8 \
  --min-dist 20000 --max-dist 2000000 --pearson 0.4 \
  sample.cool::/resolutions/5000 sample_loops
# output: sample_loops.tsv -> chrom1,start1,end1,chrom2,start2,end2,bin1,bin2,score,pvalue,qvalue

The score is a Pearson correlation (-1..1) between a loop kernel and each windowed submatrix. The same engine finds borders and stripes by swapping --pattern (loops, loops_small, borders, hairpins, centromeres, stripes_left, stripes_right) -- but stripes are a separate feature class, NOT loops. --pearson is the correlation cutoff; raise it for fewer, higher-confidence calls.

Scale-Space with Mustache

mustache -f sample.mcool -r 5000 -o loops.tsv -pt 0.1 -st 0.88 -norm weight -p 8
# output: BIN1_CHR BIN1_START BIN1_END BIN2_CHR BIN2_START BIN2_END FDR DETECTION_SCALE

-pt is the FDR/p-value threshold (default 0.1), -st the sparsity filter (default 0.88), -norm weight for a balanced .cool (KR for .hic). Mustache spans loop sizes natively via Difference-of-Gaussians across scales, which is why it adapts to kb-resolution Micro-C better than fixed-kernel HiCCUPS.

Aggregate Peak Analysis (APA) -- Confirm, Don't Discover

Goal: Quantify whether a loop SET is enriched on average and visually QC the call set (a clean aggregate dot = mostly real; a smeared/absent center = contaminated).

Approach: Compute expected, pile up observed/expected snippets centered on each anchor pair, average across the stack, then report the APA score = center pixel divided by an off-diagonal corner-control block. Pass expected_df so snippets are O/E and comparable across genomic separations.

import numpy as np
import cooltools

expected = cooltools.expected_cis(clr, view_df=arms, nproc=4)
stack = cooltools.pileup(clr, loops, view_df=arms, expected_df=expected, flank=100_000, nproc=4)   # bedpe two-anchor features
apa = np.nanmean(stack, axis=0)   # pileup returns (n_snippets, D, D); average over axis 0 -> 2D aggregate

center = apa.shape[0] // 2
corner = 3   # 3x3 corner-control block (Rao 2014 lower-left convention)
apa_score = apa[center, center] / np.nanmean(apa[-corner:, :corner])   # center vs lower-left corner; >1 = enriched

Differential Loops -- Union Anchors, Not a Bin-Level Tool

Goal: Find loops whose strength changes between conditions.

Approach: There is NO DESeq-for-loops. Build a UNION anchor set across conditions, then score each loop's strength per condition at a FIXED coordinate set (chromosight quantify, which is purpose-built for this, or APA per condition), then test the strength delta.

# quantify scores a FIXED coordinate set; arg order: <bed2d> <contact_map> <prefix>
chromosight quantify --pattern loops union_anchors.bed2d condA.cool condA_q
chromosight quantify --pattern loops union_anchors.bed2d condB.cool condB_q
# compare the per-loop score columns; or diff_mustache.py -f1 A -f2 B -pt 0.05 -pt2 0.1 -r 5000 -o diff

diffHic, multiHiCcompare, and dcHiC operate on BINS or COMPARTMENTS, not focal loops -- do NOT use them as a loop-differential tool. Cross-reference hic-differential for the bin/compartment regime.

Per-Method Failure Modes

De-novo calling on a shallow map

Trigger: running dots/chromosight/Mustache on tens of millions of pairs or >=25kb bins. Mechanism: focal signal is below the noise floor without depth. Symptom: zero or a handful of scattered, irreproducible calls. Fix: STOP de-novo; run APA on a known anchor set (CTCF/cohesin ChIP or reference loops).

APA reported without a corner control

Trigger: quoting the aggregate center-pixel value as the loop "strength." Mechanism: without an off-diagonal corner the number has no on-vs-off baseline. Symptom: a "high" APA that reflects distance-decay, not looping. Fix: APA score = center / corner-control block (Rao 2014 lower-left convention).