bio-hi-c-analysis-loop-calling

What this skill does

## 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

| Tool | Philosophy | Mechanism | When |
|------|-----------|-----------|------|
| cooltools `dots` | local enrichment (CPU HiCCUPS) | pixel must beat 4 local backgrounds (donut/horizontal/vertical/lower-left); Poisson p; lambda-binned BH-FDR | cooler/.mcool pipelines, the modern default; pure-CPU |
| Juicer HiCCUPS | local enrichment (GPU original) | same 4-kernel model on `.hic`; CUDA-bound | `.hic`/Juicer ecosystem with a GPU available |
| chromosight `detect` | template correlation | Pearson correlation of a loop/border/stripe kernel vs each window | want loops AND borders AND stripes from one engine; Micro-C-friendly |
| Mustache | scale-space blobs | Difference-of-Gaussians across scales; multi-scale catches loops of different sizes | mixed loop sizes, kb-resolution Micro-C, recovers more E-P/ChIA-PET loops |
| SIP | image processing | Gaussian blur + regional-max + watershed | `.hic` image-based alternative |
| cooltools `pileup` (APA) | CONFIRMATION, not discovery | aggregate snippets centered on an anchor set; measure center vs corner | validate/quantify a loop SET; works on shallow maps |

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

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Shallow map (tens of M pairs) | APA/pileup on KNOWN anchors (CTCF/cohesin ChIP or reference loops) -- STOP de-novo | callers return noise below ~5-10kb resolution |
| Deep cooler/.mcool, CPU only | `cooltools.dots` (default) | pure-CPU HiCCUPS reimplementation on balanced cooler |
| Deep `.hic` with a GPU | Juicer HiCCUPS | CUDA original built for billion-contact `.hic` scans |
| Mixed loop sizes / kb Micro-C | Mustache | scale-space natively spans loop sizes; sub-5kb-friendly |
| Want stripes/borders too | chromosight (swap `--pattern`) | same template engine; stripes are a SEPARATE class, not loops |
| Validate a call set | `cooltools.pileup` -> APA score vs corner control | aggregate enrichment + visual QC of the dot |
| Confirm anchors are real loops | convergent-CTCF check -> chip-seq/peak-annotation, atac-seq/footprinting | extrusion loops carry convergent CTCF motifs |
| Annotate loop anchors | -> chip-seq/peak-annotation, atac-seq/enhancer-gene-linking | E-P / TF context lives there |
| Anchor-overlap enrichment p-value | -> genome-intervals/overlap-significance | turn an anchor-overlap count into a permutation test |
| Two conditions, loop strength shift | union anchors -> chromosight `quantify` per condition -> test delta (or `diff_mustache`) | NO bin-level DESeq for loops; quantify a fixed coordinate set |
| HiChIP / PLAC-seq / PCHi-C | FitHiChIP / MAPS / HiC-DC+ -> chip-seq/peak-calling | protein-anchored, coverage-biased; HiCCUPS null is wrong |

## 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`.

```python
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