bio-hi-c-analysis-matrix-operations

What this skill does

## Version Compatibility

Reference examples tested with: cooler 0.10+, cooltools 0.7+, bioframe 0.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 standardized its API around 0.7 (functions take a `view_df` viewframe; `expected_cis` defaults to `smooth=True, aggregate_smoothed=True`). `.mcool` is multi-resolution: analysis functions take a single-resolution URI (`file.mcool::/resolutions/10000`), never the bare `.mcool`. A matrix must be balanced (a stored `weight` column) before O/E, compartments, insulation, or dots; `clr.matrix(balance=True)` on an unbalanced cooler returns all-NaN.

# Hi-C Matrix Operations

**"Make the pixels of my Hi-C matrix comparable to each other."** -> Balance (remove per-bin coverage bias under equal-visibility), then divide by distance-matched expected (remove the polymer P(s) background) to get O/E.
- Python: `cooler.balance_cooler(clr, cis_only=True, store=True)`, then `cooltools.expected_cis(clr)` and divide observed by per-diagonal expected.

## The Single Most Important Modern Insight -- Balancing Makes ONE Matrix Self-Consistent; It Does NOT Make Two Matrices Comparable

Balancing (ICE/KR) is a *within-matrix* operation: it solves for per-bin bias weights so every bin has equal genome-wide visibility, making a single map internally consistent. It does nothing to relate map A to map B. Two balanced matrices at different sequencing depth still differ in absolute magnitude, dynamic range, and noise floor -- and `rescale_marginals` makes the absolute balanced values arbitrary-scaled anyway. "I balanced both, now I'll subtract/log2-ratio them" is the single most common error in the field: the depth difference is read as biology. Cross-sample comparison requires downsampling to equal valid-pair count, distance-matched O/E, and a replicate-aware differential tool (multiHiCcompare, HiCcompare loess-over-distance, dcHiC) -- route to hic-differential.

Two corollaries that flow from the same equal-visibility model:

1. **CNV silently breaks balancing.** The premise is that every bin *should* make the same number of contacts; any deviation is technical bias. That is true for a diploid uniform-copy genome and FALSE for tumors/aneuploids -- a 3-copy region genuinely contacts ~3x more. ICE forces equal marginals and ERASES that real copy-number, then redistributes it perversely (post-ICE high-copy regions go cis-depleted, trans-enriched; Servant 2018). Use **raw counts** for CNV/SV calling (coverage IS the signal -> copy-number); plain ICE/KR on an aneuploid is a category error (use CNV-aware LOIC/CAIC for 3D structure).
2. **A balanced cis map is still dominated by distance-decay.** The A/B plaid and focal loops are a faint modulation under the P(s) background. Dividing by distance-matched expected (O/E) before eigendecomposition is mandatory, or the top eigenvector is just the decay curve, not compartments.

## Normalization-Method Taxonomy

| Method | What it does | Mechanism | When |
|--------|-------------|-----------|------|
| ICE (cooler native) | true matrix balancing | iterative proportional fitting (Sinkhorn); equalizes all marginals | default; robust, converges on sparse/low-depth where KR fails |
| KR (juicer) | true matrix balancing | Knight-Ruiz Newton solver; SAME fixed point as ICE | fast (few iterations); fails to converge on sparse/high-res maps |
| SCALE (juicer) | true matrix balancing | modern KR-family solver, more robust | juicer's default; converges where KR diverges on sparse maps |
| VC / vanilla coverage | NOT true balancing | single pass: divide by row-coverage * col-coverage (one ICE iteration) | fast robust fallback; leaves residual bias |
| VC_SQRT | NOT true balancing | divide by sqrt of coverage product (gentler than VC) | very sparse data where full balancing overfits |
| LOIC / CAIC (Servant 2018) | CNV-aware balancing | condition on copy-number; LOIC keeps the CN effect, CAIC removes it | aneuploid/tumor genomes (plain ICE is wrong here) |

KR and ICE reach the same balanced map -- choose by **convergence, not quality**: KR is faster but blows up on sparse/low-depth/high-resolution matrices; ICE is the robust default; SCALE is the juicer-side answer when KR fails.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Need to balance a diploid map | `cooler.balance_cooler(cis_only=True, store=True)` (ICE) | robust default; cis-only is the analysis convention |
| KR failed to converge (sparse/high-res) | fall back to ICE, or SCALE on the juicer side | same fixed point; ICE/SCALE are the robust solvers |
| Tumor / aneuploid genome, 3D structure | CNV-aware LOIC/CAIC (Servant 2018) | plain ICE erases real copy-number |
| CNV / SV calling from Hi-C | RAW counts (no balancing) | coverage is the signal -> copy-number |
| Compartments at 100kb-1Mb | balance -> `expected_cis` -> O/E -> Pearson -> eigenvector | O/E removes P(s) so the plaid is visible -> compartment-analysis |
| Focal loops / TADs | balance -> `expected_cis` -> O/E | local enrichment needs the distance background removed -> loop-calling, tad-detection |
| P(s) / polymer-state diagnostic | `expected_cis(smooth=True)` -> log-derivative | the derivative reads out loop-extrusion machinery |
| Imported a juicer KR/VC weight column | check `divisive_weights` before applying | cooler weights are multiplicative, juicer's are divisive |
| Compare two conditions | downsample to equal depth, then -> hic-differential | balancing is within-matrix, not a cross-sample normalizer |
| Bare `.mcool` passed and KeyError | use `file.mcool::/resolutions/<bp>` URI | `.mcool` is a container of resolutions |

## Balance a Matrix (ICE)

**Goal:** Remove one-dimensional per-bin coverage bias so every bin has equal genome-wide visibility within this single map.

**Approach:** Mask low-coverage and blacklisted bins FIRST (mad_max on log-marginals + explicit blacklist of centromere/rDNA/unmappable), drop the first two diagonals (ligation chemistry, not 3D contact), then run cis-only ICE; the multiplicative weight vector is stored in the `weight` column.

```python
import cooler

clr = cooler.Cooler('matrix.mcool::/resolutions/10000')
bias, stats = cooler.balance_cooler(clr, cis_only=True, mad_max=5, ignore_diags=2, blacklist=None, store=True)
print('converged:', stats['converged'], 'scale:', stats['scale'])   # stats also reports var, divisive_weights

clr = cooler.Cooler('matrix.mcool::/resolutions/10000')              # re-open to see the stored weights
balanced = clr.matrix(balance=True).fetch('chr1')                    # raw[i,j] * w[i] * w[j]; masked bins -> NaN
```

`cis_only=True` is the convention for compartment/TAD/loop work -- trans signal is weak, noisy ambient ligation that pulls the bias estimates toward trans noise. `ignore_diags=2` drops the main diagonal (self-ligation/dangling ends) and first off-diagonal (undigested/religated fragments): huge untrustworthy counts that would dominate the marginals. `mad_max=5` filters bins whose **log**-marginal is >5 MAD below the median; without it a near-empty unmappable/centromeric bin gets a gigantic weight and ICE diverges. Masking (mad_max + `blacklist`) MUST precede balancing -- balancing cannot rescue a no-signal bin, it amplifies it.

CLI equivalent:

```bash
cooler balance --cis-only --mad-max 5 --ignore-diags 2 matrix.mcool::/resolutions/10000
```

## Expected: cis P(s) Curve and trans Scalar

**Goal:** Build the distance-decay background (the denominator for O/E) and the P(s) curve for diagnostics.

**Approach:** cis expected is a per-diagonal curve (one value per separation s -- this IS P(s)); trans expected