Bio Hi C Analysis Compartment Analysis

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

cooltools had a major API shift around 0.5 -> 0.7+ (functions standardized on view_df/viewframe arguments; eigs_cis, expected_cis, saddle signatures changed). The cooler MUST be balanced before any compartment analysis: clr.matrix(balance=True) requires a stored weight column. A .mcool is multi-resolution -- pass a single-resolution URI (file.mcool::/resolutions/100000), not the bare .mcool. The phasing_track MUST share the cooler's exact binning or phasing silently no-ops. If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

A/B Compartment Analysis

"Which regions of my genome are in the active vs inactive compartment?" -> Distance-normalize the cis matrix, take an eigenvector of its Pearson correlation matrix, then orient it by GC/gene density so positive = A (active) -- but verify the kept eigenvector is the compartment one, not an arm gradient.

• Python: cooltools.eigs_cis(clr, gc_track, view_df=arms, n_eigs=3, sort_metric='pearsonr')

The Single Most Important Modern Insight -- E1 Is a Choice, Not an Output, and Its Sign Is Arbitrary Until Phased

The two most damaging beginner assumptions are "E1 = compartments" and "positive E1 = active." Both are false out of the box, and both fail silently -- the pipeline runs, returns a track, and is wrong.

1. E1 is not guaranteed to be the compartment track. On a whole-chromosome O/E correlation matrix the largest eigenvalue very often belongs to a smooth p-arm-vs-q-arm or centromere-to-telomere GRADIENT, not the plaid A/B checkerboard -- the real compartment signal then lands in E2 or E3. cooltools' own docs concede the first eigenvector "occasionally describes chromosomal arms or translocation blowouts." The compartment eigenvector is the one with the largest |correlation| to an activity track (GC, gene density, H3K27ac), not the one with the largest eigenvalue. The structural fix removes the gradient at the source: run eigs_cis per chromosome ARM (a view_df split at centromeres, from bioframe.make_chromarms), so the arm gradient is never in the within-arm matrix. Set sort_metric='pearsonr' so the returned eigenvectors are ordered by GC correlation, not eigenvalue -- otherwise the arm gradient is reported as "E1." A monotonic "compartment track" with no sign flips across a chromosome is the failure signature of a captured arm gradient.

2. The sign is arbitrary until phased. Eigenvectors are defined up to sign; the positive lobe is meaningless and can differ per chromosome AND per sample. The eigenvector MUST be oriented with an external active-chromatin track via the phasing_track argument so A = positive. GC content is the field default (it needs no extra assay and tracks compartment A; Lieberman-Aiden 2009 Science 326:289) -- compute it with bioframe.frac_gc at the compartment resolution, exactly matching the cooler's binning. Wrong/weak phasing flips A<->B silently, and every downstream saddle, switch call, and differential result inverts with no error. This is a classic source of irreproducible compartment papers.

3. Compartments are an equilibrium phenomenon decoupled from TADs/loops. Compartments = microphase separation of A/B chromatin states (cohesin-independent; survive cohesin loss, Schwarzer 2017 Nature 551:51, and CTCF loss, Nora 2017 Cell 169:930). TADs/loops = ATP-driven loop extrusion stalled at CTCF. Removing cohesin reinforces compartments while erasing TADs (Schwarzer 2017 reports reinforced compartmentalization on Nipbl loss; Rao 2017 Cell 171:305 eliminates all loop domains with compartments retained) -- loop extrusion actively mixes chromatin across compartment boundaries, so removing the extruder lets microphase separation run to completion (Nuebler 2018 PNAS 115:E6697). A preserved-or-stronger saddle after a cohesin/Nipbl/RAD21 perturbation is the EXPECTED result, not a bug; compartment-strength and TAD-strength are antagonistic. If a CTCF/cohesin perturbation makes compartments vanish, suspect a phasing artifact, not biology.

Method / Output Taxonomy

Decision Tree by Scenario

Per-Arm Eigenvector with GC Phasing

Goal: Assign each genomic bin to the active (A) or inactive (B) compartment with a non-arbitrary sign, avoiding the centromere arm-gradient artifact.

Approach: Build a per-arm view_df (split at centromeres) so the arm gradient never enters the matrix; compute a GC-content phasing track at the cooler's exact binning; run eigs_cis with the GC track and sort_metric='pearsonr' so eigenvectors are ordered by GC correlation; then take the GC-correlated eigenvector as the compartment track.

import cooler
import cooltools
import bioframe

clr = cooler.Cooler('matrix.mcool::/resolutions/100000')   # 100kb: compartments are coarse-scale
chromsizes = clr.chromsizes
cens = bioframe.fetch_centromeres('hg38')
arms = bioframe.make_chromarms(chromsizes, cens)            # per-arm view removes the centromere gradient
arms = arms[arms.chrom.isin(clr.chromnames)].reset_index(drop=True)

genome = bioframe.load_fasta('hg38.fa')                     # FASTA index (.fai) must exist
bins = clr.bins()[:][['chrom', 'start', 'end']]
gc = bioframe.frac_gc(bins, genome)                         # phasing track at the cooler's exact binning

eigvals, eigvecs = cooltools.eigs_cis(clr, gc, view_df=arms, n_eigs=3, sort_metric='pearsonr')
eigvecs['compartment'] = ['A' if e > 0 else 'B' for e in eigvecs['E1']]   # GC-phased: positive E1 = A

After phasing, sanity-check that E1 correlates with gc['GC'] (sign and magnitude). If the strongest correlation is in E2/E3, that component -- not E1 -- is the compartment track; re-derive the call from it.

Compartment Strength via Saddle Plot

Goal: Quantify how strongly the genome demixes into A and B with a single comparable number across conditions.

Approach: Compute the distance-decay expected; pass the cooler, the expected, and the phased E1 eigenvector to saddle, which digitizes E1 into quantile groups internally (via qrange) and aggregates O/E into a 2D table of same-vs-cross-compartment interactions; then read saddle_strength (the (AA+BB)/(AB+BA) corner ratio) at one fixed extent, applied identically to every sample compared.

N_GROUPS = 38            # quantile groups for digitizing E1; ~30-50 is conventional (cooltools tutorial)
Q_LO, Q_HI = 0.025, 0.975   # trim the extreme 2.5% tails before digitizing to resist outlier bins

expected = cooltools.expected_cis(clr, view_df=arms)        # has the 'balanced.avg' column saddle needs
track = eigvecs[['chrom', 'start', 'end', 'E1']]           # the SAME phased E1 used for the A/B call
interaction_sum, interaction_count = cooltools.saddle(
    clr, expected, track, 'cis', n_bins=N_GROUPS, qrange=(Q_LO, Q_HI), view_df=arms
)
strength = cooltools.api.saddle.saddle_strength(interaction_sum, interaction_count)   # 1D array; lives in cooltools.api.saddle, not top level
EXTENT = N_GROUPS // 5   # read strength at the top/bottom ~20% of bins; pick one extent, use it everywhere
score = strength[EXTENT]

saddle_strength returns an ARRAY (cumulative corner ratio over increasing extent), not a scalar -- there is no canonical single number, so choose an extent and apply it identically across compared samples. Mismatched n_bins, resolution, qrange, or extent make strengths incomparable. Remember: a preserved-or-higher strength after cohesin/Nipbl/RAD21 loss is the expected result.

Compare Compartments Across Conditions

Goal: Find bins that switch A<->B between two conditions without being fooled by per-sample sign flips.

Approach: Do NOT independently phase two eigenvectors and diff them bin-by-bin -- a weak per-chromosome GC correlation can flip the sign in one sample only, manufacturing fake "switches." Use dcHiC, which computes eigenvectors on quantile-normalized scores in a shared framework (sign-coherent across the cohort) and reports a multivariate significance per bin. Route this to hic-differential.

# Hand-diffing is only safe when you have CONFIRMED both eigenvectors are sign-coherent (same arms phased
# to the same GC track with strong correlation). Otherwise use dcHiC -- see hic-differential.
import pandas as pd
merged = eig1.merge(eig2, on=['chrom', 'start', 'end'], suffixes=('_1', '_2'))
merged['switch'] = (merged['E1_1'] > 0) != (merged['E1_2'] > 0)   # only meaningful if both are phased coherently

Per-Method Failure Modes

Eigenvector captured the arm gradient

Trigger: eigs_cis run per whole chromosome (no per-arm view_df) and/or sort_metric=None. Mechanism: the largest eigenvalue belongs to the smooth p-vs-q arm / centromere gradient, not the A/B checkerboard. Symptom: a monotonic "compartment track" across a chromosome with no sign flips; weak correlation of E1 with GC. Fix: run per chromosome arm (bioframe.make_chromarms); set sort_metric='pearsonr'; pick the eigenvector with the largest |corr| to GC.