Bio Hi C Analysis Hic Differential

Version Compatibility

Reference examples tested with: cooler 0.10+, cooltools 0.7+, bioframe 0.7+, multiHiCcompare 1.20+, diffHic 1.34+, dcHiC (2022 release), hicrep (Bioconductor) 1.26+, edgeR 4.0+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• 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.

Note: cooltools provides the Python feature-extraction parts (expected, eigenvectors, insulation, pileups) but NO turnkey two-condition test -- the differential statistics live in R/Bioconductor (multiHiCcompare, diffHic, dcHiC, hicrep). The get.scc signature changed between the TaoYang-dev and current Bioconductor/qunhualilab releases; verify with ?get.scc before calling. A .mcool is multi-resolution: pass a single-resolution URI (file.mcool::/resolutions/10000), not the bare file.

Hi-C Differential Analysis

"What changed in 3D genome organization between my conditions?" -> Equalize the two maps with a distance-stratified between-sample normalization, then test at the SCALE of the question (compartment, TAD boundary, loop, or bin-pair) with a replicate-aware method, not pixel-wise log2 subtraction.

• Python (features): cooltools.expected_cis(clr), cooltools.eigs_cis(clr), cooltools.insulation(clr)
• R (bin-pair test): make_hicexp(...) |> cyclic_loess() |> hic_exactTest() |> results() (multiHiCcompare)
• R (compartments): Rscript dchicf.r --pcatype cis|select|analyze (dcHiC)

The Single Most Important Modern Insight -- balancing makes a matrix self-consistent, NOT cross-comparable

ICE/KR/SCALE balancing equalizes the marginals WITHIN one map (it removes per-bin visibility bias). It says nothing about whether map A and map B are on the same footing. Two balanced maps still differ in (a) total sequencing depth and (b) cis/trans ratio, and a naive log2(A/B) is dominated by those two nuisances plus the shared distance-decay P(s) -- with biology buried underneath. "I balanced both, so I can subtract them" is the single most common error in differential Hi-C. The fix is a BETWEEN-sample, DISTANCE-STRATIFIED normalization (multiHiCcompare's cyclic loess on the M-D plot, or diffHic's trended loess offsets) before any difference is interpretable.

The M-D plot is the RNA-seq MA-plot's distance-aware cousin: M = log2(IF1/IF2), but plotted against genomic DISTANCE D instead of mean abundance. The loess fit is done PER distance stratum because both bias and variance depend on distance -- a sparser library loses long-range pairs faster, so the depth bias is itself distance-dependent and a single global size-factor cannot fix it. After normalization M should center on 0 at every D; a residual M-trend at large D means normalization failed at long range. Inspect MD_hicexp() -- do not trust the call set blind.

Differential-Method Taxonomy (scale-matched -- one tool cannot do all four)

Using a bin-pair tool (multiHiCcompare/diffHic) to "find compartment changes" is a category error: compartments are an eigenvector property of the whole chromosome, not a sum of independent bin-pairs, so the result is a noisy bin-pair list, not coherent compartment calls. Match the tool to the scale of the question.

Decision Tree by Scenario

Replicate QC Gate with HiCRep SCC (do this FIRST)

Goal: Decide whether replicates are reproducible enough that a differential call is meaningful at all.

Approach: Plain Pearson on Hi-C always looks reproducible -- the shared P(s) decay alone drives r > 0.9 even between unrelated maps. HiCRep's SCC stratifies by distance (removing the decay) and smooths for sparsity, then variance-weights the strata. Compute SCC for within-condition replicate pairs and between-condition pairs; within must clearly exceed between or there is nothing to call.

library(hicrep)

# Current Bioconductor/qunhualilab interface: dat is a 4-column table
# (mid1, mid2, IF_A, IF_B); resol = bin size; max = max distance considered.
# Verify with ?get.scc -- the older TaoYang-dev interface is get.scc(mat1, mat2, resol, h, lbr, ubr).
scc_out <- get.scc(dat_repA_vs_repB, resol = 50000, max = 5000000)
scc_out$scc   # stratum-adjusted correlation coefficient in [-1, 1]

A differential claim is only meaningful when within-condition SCC clearly exceeds between-condition SCC. If they overlap, the "differential" signal is replicate noise.

Differential Bin-Pair Contacts with multiHiCcompare

Goal: Find individual bin-pairs whose contact frequency changes between conditions, with a calibrated FDR from replicate variance.

Approach: Build a Hi-C experiment from per-replicate sparse upper-triangular tables, normalize jointly across all samples with cyclic loess on the M-D plot, then run edgeR's exact test (2 groups) or GLM (covariates). Filtering on mean abundance happens at make_hicexp time -- it is independent of the contrast, so it shrinks the multiple-testing burden without inflating FDR.

library(multiHiCcompare)

# Each replicate is a 4-column sparse table: chr, region1(bp), region2(bp), IF
# chr coded 1-22, 23=X, 24=Y.
hicexp <- make_hicexp(c1_r1, c1_r2, c2_r1, c2_r2,
                      groups = c(0, 0, 1, 1),
                      zero.p = 0.8,                 # drop bin-pairs >80% zero across samples
                      A.min  = 5,                   # drop bin-pairs with low mean IF (independent filter)
                      filter = TRUE,
                      remove.regions = hg19_cyto)   # blacklist centromeres/telomeres
hicexp <- cyclic_loess(hicexp, span = NA)           # span=NA -> GCV chooses the loess span
hicexp <- hic_exactTest(hicexp)                     # 2-group; use hic_glm(hicexp, design) for covariates
res <- results(hicexp)                              # chr, region1, region2, D, logFC, logCPM, p.value, p.adj
sig <- topDirs(hicexp, logfc_cutoff = 1, logcpm_cutoff = 1, p.adj_cutoff = 0.1, return_df = 'pairedbed')
MD_hicexp(hicexp)                                   # diagnostic: M should center on 0 at every D

diffHic is the alternative when full edgeR control is wanted: squareCounts -> filterDirect -> normOffsets(type='loess') -> asDGEList -> estimateDisp -> glmQLFit -> glmQLFTest. Same NB-GLM engine as edgeR/csaw; requires biological replicates for dispersion.

Differential A/B Compartments with dcHiC

Goal: Detect compartment changes between conditions, including graded shifts that never cross the A/B boundary, with cross-sample-comparable eigenvectors.

Approach: Everyone computes PC1, but PC1's sign is arbitrary per chromosome per sample and sometimes the compartment signal is in PC2 -- naive multi-sample comparison silently compares flipped or mismatched axes. dcHiC anchors the sign (GC/gene-density) and selects the correct PC, quantile-normalizes the scores, then uses a multivariate Mahalanobis distance per bin to flag outliers across all samples. Run as a staged CLI; the input file lists <matrix> <bed> <replicate_prefix> <experiment_prefix> per replicate (no dashes/dots in prefixes).

# Staged CLI (dchicf.r). cis = per-sample compartments; select = pick PC per chr;
# analyze = differential PCA (Mahalanobis); viz = IGV-style browser.
Rscript dchicf.r --file input.txt --pcatype cis    --dirovwt T --cthread 2 --pthread 4
Rscript dchicf.r --file input.txt --pcatype select --dirovwt T --genome hg38
Rscript dchicf.r --file input.txt --pcatype analyze --dirovwt T --diffdir cond1_vs_cond2
Rscript dchicf.r --file input.txt --pcatype viz    --diffdir cond1_vs_cond2 --genome hg38
# Differential calls land in DifferentialResult/<diffdir>/fdr_result/ ; optional: --pcatype subcomp (HMM) and dloop.

For a quick Python eyeball of compartment switching (NOT a replicate-aware test -- use dcHiC for that), difference the phased E1 from cooltools per condition and flag sign flips. State it as exploratory.

Differential TAD Boundaries via Delta Insulation

Goal: Find boundaries that strengthen, weaken, or appear/disappear between conditions.

Approach: Boundary changes are about local insulation DEPTH, not bin-pair counts, so compute the Crane-style insulation score per condition with cooltools at a fixed window, then difference the scores; attach significance with replicate insulation profiles. The window must match the bin size (typically 5-25x the bin).

import cooltools

ins1 = cooltools.insulation(clr1, window_bp=[200000], ignore_diags=2)   # window ~5-25x bin size
ins2 = cooltools.insulation(clr2, window_bp=[200000], ignore_diags=2)
merged = ins1.merge(ins2, on=['chrom', 'start', 'end'], suffixes=('_1', '_2'))
merged['delta_insulation'] = merged['log2_insulation_score_200000_2'] - merged['log2_insulation_score_200000_1']

GENOVA and FAN-C also compute the insulation score for differencing; significance is added separately (paired test across replicate insulation tracks).

CNV Correction Before Cancer Comparisons

Goal: Avoid calling copy-number differences as 3D-structure differences in tumor-vs-normal data.

Approach: Contact count scales with copy number, and balancing assumes equal visibility -- so ICE on aneuploid data SHIFTS contacts between amplified and deleted regions instead of removing the artifact, lighting up enormous spurious "differential interaction" blocks. Either regress out the marginal (1D) coverage log-ratio as a covariate (diffHic marginCounts + normalizeCNV, or OneD's 1D GAM), OR call CNV first (HiNT/HiCnv) and interpret rearrangement blocks separately. Never run vanilla balanced-matrix differential on cancer data.

library(diffHic)

margins <- marginCounts(data)              # 1D marginal coverage per bin (a RangedSummarizedExperiment)
nb.off  <- normalizeCNV(data, margins)     # 2D loess on (abundance, marginal log-ratio) -> GLM offsets; matches margins internally
y <- asDGEList(data)
y$offset <- nb.off                         # attach the CNV offsets; asDGEList does not carry them automatically
# then estimateDisp -> glmQLFit -> glmQLFTest as usual

Per-Method Failure Modes

Naive log2 of two balanced maps

Trigger: log2((mat2+1)/(mat1+1)) on two balanced coolers of different depth. Mechanism: balancing is within-map only; depth + cis/trans + P(s) differences dominate. Symptom: a smooth distance-dependent gradient in the "differential" map, strongest off-diagonal. Fix: distance-stratified between-sample normalization (multiHiCcompare cyclic_loess / diffHic loess offsets) before differencing.