bio-hi-c-analysis-tad-detection

What this skill does

## Version Compatibility

Reference examples tested with: cooler 0.10+, cooltools 0.7+, bioframe 0.7+, HiCExplorer 3.7+

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.

cooltools changed its API around 0.5 -> 0.7+ (functions standardized on `view_df`/viewframe arguments; `insulation` returns the `boundary_strength_{W}`/`is_boundary_{W}` columns). A `.cool` MUST be balanced (a stored `weight` column) before insulation; `clr.matrix(balance=True)` on an unbalanced cooler returns all-NaN. A `.mcool` is multi-resolution: pass a single-resolution URI (`file.mcool::/resolutions/10000`), never the bare `.mcool`.

# TAD Detection

**"Where are the reproducible domain boundaries in my Hi-C matrix, and how strong?"** -> Compute the diamond-window insulation score on the balanced matrix, take valley minima as boundaries and their prominence as strength, and report across a LIST of window sizes rather than a single magic scale.
- Python: `cooltools.insulation(clr, [3*res, 5*res, 10*res, 25*res])` then rank by `boundary_strength_{W}`
- CLI: `hicFindTADs -m corrected.cool --outPrefix tads --correctForMultipleTesting fdr` (sweep `--minDepth/--maxDepth/--step`)

## The Single Most Important Modern Insight -- The Boundary Is Real; the Domain Is Mostly an Averaging Artifact

A population Hi-C "TAD" is the ensemble average over a heterogeneous mixture of cell-specific, stochastic domains. Single-cell imaging (Bintu 2018 *Science* 362:eaau1783) shows individual cells DO have sharp domains, but the boundary POSITION varies cell to cell - the population boundary is a *preferred* position, not a wall. Cohesin depletion abolishes population TADs while leaving single-cell domains intact, removing only the preferred-position bias. Three consequences govern every decision in this skill:

1. **"How many TADs are there" is the wrong question.** TAD number and size vary 2-5x across caller, resolution, and normalization, with NO ground truth (Forcato 2017 *Nat Methods* 14:679; Zufferey 2018 *Genome Biol* 19:217). Insulation valleys and directionality-index sign-changes give DIFFERENT boundary sets on the same matrix. The quoted "average TAD size ~880kb" is an artifact of one caller at one resolution - Zufferey explicitly states there is no average TAD size.
2. **The boundary is the reproducible, mechanistic unit.** Strong, CTCF-backed boundaries survive caller and resolution swaps; weak boundaries and exact domain extents do not. Prefer the continuous insulation/boundary-strength track over a hard domain partition for any downstream claim.
3. **The diamond window IS the analysis.** Small window (~3x bin) -> sub-TAD/fine boundaries; large window (~25x bin) -> compartment-scale domains - same matrix, totally different "TADs." Always run a list of windows and report multi-scale; never report a single-window partition as ground truth.

## TAD-Caller Taxonomy

| Method | Role | Mechanism | When |
|--------|------|-----------|------|
| cooltools `insulation` | boundary score + strength | diamond-window valleys; prominence = strength; Li/Otsu threshold | cooler-native pipeline, multi-scale, the modern default |
| HiCExplorer `hicFindTADs` | domains + boundaries + FDR | multi-window TAD-separation score with per-bin multiple-testing | CLI workflow, hierarchical sweep, FDR-controlled boundaries |
| directionality index (DI) | boundary direction | HMM on up/downstream interaction bias (Dixon 2012) | classic comparison, legacy reproducibility, gives a partition |
| Arrowhead (Juicer) | corner-score domains | arrowhead transform on the `.hic`; loop-anchored contact domains | Juicer/.hic ecosystems; calls fewer, sharper domains (Rao 2014, median ~185kb) |
| OnTAD / TADtree / rGMAP | nested/hierarchical | explicitly models meta-TADs > TADs > sub-TADs | when the question is about hierarchy or sub-TAD insulation (An 2019) |
| Stripenn / JOnTADS | stripe-aware | calls asymmetric stripes as first-class objects | when the map shows flames/stripes a flat caller mis-segments |

Insulation/DI/hicFindTADs are blind to stripes (asymmetric one-sided extrusion). Arrowhead "contact domains" are corner-anchored and categorically different from track-based boundaries - do NOT cross-compare their counts naively.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Matrix not yet balanced | `cooler balance` / `cooler.balance_cooler` first | unbalanced insulation = coverage-driven garbage valleys |
| Reproducible boundaries, one sample | `insulation` multi-window, rank by `boundary_strength_{W}` | strength is continuous and comparable; partition is brittle |
| "What scale of domain?" | run windows `[3,5,10,25]x` bin; ~10x is the mammalian sweet spot | the window sets sub-TAD vs TAD vs compartment-domain |
| Need FDR-controlled domains | `hicFindTADs` with `--correctForMultipleTesting fdr`, sweep depths | per-bin multiple-testing on the TAD-separation score |
| Hierarchical/nested structure | OnTAD/TADtree (or compare windows) | a flat caller picks ONE level set by its window |
| Asymmetric stripes/flames present | Stripenn/JOnTADS | insulation-only pipelines are blind to stripes |
| Two conditions, boundary change | differential SCORE at matched bins, NOT intersected domain BEDs | partitions are unstable; set-differencing manufactures spurious gain/loss |
| Annotate boundaries with CTCF | -> chip-seq/peak-annotation, genome-intervals/overlap-significance | ~76-85% of boundaries are convergent CTCF + cohesin |
| Overlap boundaries with features | -> genome-intervals/interval-arithmetic | boundary BED set operations live there |
| Render domains on the matrix | -> hic-visualization | the TAD square is a colormap/resolution choice as much as a measurement |

## Insulation Score and Boundaries (multi-scale)

**Goal:** Produce a continuous boundary-strength track and threshold-flagged boundaries at several scales, so the analysis reports where insulation reproducibly dips rather than a single brittle partition.

**Approach:** Run `cooltools.insulation` on the balanced cooler with a LIST of window sizes (3-25x the bin). Each window appends its own `log2_insulation_score_{W}` (valleys = boundaries), `boundary_strength_{W}` (valley prominence - the quantitative, comparable strength), and `is_boundary_{W}` (the prominence passed through a Li histogram threshold). Rank and compare on `boundary_strength`, not on the boolean flag.

```python
import cooler
import cooltools

clr = cooler.Cooler('matrix.mcool::/resolutions/10000')   # single-resolution URI, must be balanced
res = clr.binsize
windows = [3 * res, 5 * res, 10 * res, 25 * res]   # 30k,50k,100k,250k: sub-TAD -> compartment-domain
ins = cooltools.insulation(clr, windows, verbose=True)   # clr_weight_name='weight' default -> needs ICE balancing

strong = ins[ins[f'is_boundary_{10 * res}']]   # 100kb window: ~10x bin, mammalian interphase sweet spot
ranked = ins.dropna(subset=[f'boundary_strength_{10 * res}']).sort_values(f'boundary_strength_{10 * res}', ascending=False)
```

`boundary_strength_{W}` is the scipy-style PROMINENCE of the insulation valley - continuous, quantitative, and comparable across samples; use it for ranking and cross-condition deltas. `is_boundary_{W}` is just that prominence passed through `threshold='Li'` (skimage threshold_li, an Otsu-like histogram split that is MORE PERMISSIVE than Otsu). Because the Li cutoff is fit per dataset, `is_boundary` is dataset-dependent and NOT directly comparable across samples - compare `boundary_strength`, then threshold consistently. `min_frac_valid_pixels` (default 0.66) and `min_dist_bad_bin` gate which bins get a score; sparse/