/bio-hi-c-analysis-hic-differential

name: bio-hi-c-analysis-hic-differential
description: Compare Hi-C contact matrices between conditions to identify differential chromatin interactions. Compute log2 fold changes, statistical significance, and visualize differential contact maps. Use when comparing Hi-C contacts between conditions.
tool_type: python
primary_tool: cooltools

Version Compatibility

Reference examples tested with: cooler 0.9+, cooltools 0.6+, matplotlib 3.8+, numpy 1.26+, pandas 2.2+, scipy 1.12+, statsmodels 0.14+

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

• Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Hi-C Differential Analysis

**"Compare Hi-C contacts between my conditions"** → Compute log2 fold-change contact maps, identify statistically significant differential interactions, and visualize changes in 3D genome organization.

• Python: `cooltools` for expected values, custom differential analysis with `scipy.stats`

Compare Hi-C contact matrices between conditions.

Required Imports

import cooler
import cooltools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from scipy import stats
import bioframe

Load Two Conditions

# Load balanced cooler files at same resolution
clr1 = cooler.Cooler('condition1.mcool::resolutions/10000')
clr2 = cooler.Cooler('condition2.mcool::resolutions/10000')

print(f'Condition 1: {clr1.info["sum"]:,} contacts')
print(f'Condition 2: {clr2.info["sum"]:,} contacts')

Compute Log2 Fold Change

def log2_fold_change(clr1, clr2, region, pseudocount=1):
    '''Compute log2(condition2/condition1) for a region'''
    mat1 = clr1.matrix(balance=True).fetch(region)
    mat2 = clr2.matrix(balance=True).fetch(region)

    # Add pseudocount and compute log2 ratio
    log2fc = np.log2((mat2 + pseudocount) / (mat1 + pseudocount))
    log2fc[np.isinf(log2fc)] = np.nan

    return log2fc

region = 'chr1:50000000-60000000'
log2fc = log2_fold_change(clr1, clr2, region)
print(f'Log2FC range: {np.nanmin(log2fc):.2f} to {np.nanmax(log2fc):.2f}')

Plot Differential Contact Map

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Condition 1
mat1 = clr1.matrix(balance=True).fetch(region)
im1 = axes[0].imshow(np.log2(mat1 + 1), cmap='Reds', vmin=-10, vmax=-3)
axes[0].set_title('Condition 1')
plt.colorbar(im1, ax=axes[0])

# Condition 2
mat2 = clr2.matrix(balance=True).fetch(region)
im2 = axes[1].imshow(np.log2(mat2 + 1), cmap='Reds', vmin=-10, vmax=-3)
axes[1].set_title('Condition 2')
plt.colorbar(im2, ax=axes[1])

# Log2 fold change (diverging colormap)
norm = TwoSlopeNorm(vmin=-2, vcenter=0, vmax=2)
im3 = axes[2].imshow(log2fc, cmap='coolwarm', norm=norm)
axes[2].set_title('Log2(Cond2/Cond1)')
plt.colorbar(im3, ax=axes[2])

plt.tight_layout()
plt.savefig('differential_hic.png', dpi=150)

Split View Comparison

def plot_split_view(mat1, mat2, title=''):
    '''Upper triangle: condition1, Lower triangle: condition2'''
    combined = np.triu(mat1) + np.tril(mat2, k=-1)

    fig, ax = plt.subplots(figsize=(8, 8))
    im = ax.imshow(np.log2(combined + 1), cmap='Reds', vmin=-10, vmax=-3)
    ax.axline((0, 0), slope=1, color='black', linewidth=0.5)
    ax.set_title(f'{title}\nUpper: Cond1, Lower: Cond2')
    plt.colorbar(im, ax=ax)
    return fig

mat1 = clr1.matrix(balance=True).fetch(region)
mat2 = clr2.matrix(balance=True).fetch(region)
fig = plot_split_view(mat1, mat2)
plt.savefig('split_view.png', dpi=150)

Depth Normalization

def depth_normalize(clr, target_depth=None):
    '''Normalize matrix to target sequencing depth'''
    total = clr.info['sum']
    if target_depth is None:
        return 1.0
    return target_depth / total

# Normalize both samples to same depth
target = min(clr1.info['sum'], clr2.info['sum'])
scale1 = depth_normalize(clr1, target)
scale2 = depth_normalize(clr2, target)

mat1_norm = clr1.matrix(balance=True).fetch(region) * scale1
mat2_norm = clr2.matrix(balance=True).fetch(region) * scale2

Statistical Testing (Per-Pixel)

**Goal:** Identify individual contact pixels that are statistically significantly different between two conditions using biological replicates.

**Approach:** For each pixel position, collect values across replicates in both conditions, apply a per-pixel t-test or Mann-Whitney U test, then correct for multiple testing with FDR.

def differential_test(matrices1, matrices2, method='ttest'):
    '''
    Test for differential contacts between replicates.
    matrices1/2: lists of numpy arrays (replicates)
    '''
    n1, n2 = len(matrices1), len(matrices2)
    shape = matrices1[0].shape

    pvalues = np.ones(shape)
    log2fc = np.zeros(shape)

    for i in range(shape[0]):
        for j in range(shape[1]):
            vals1 = [m[i, j] for m in matrices1 if not np.isnan(m[i, j])]
            vals2 = [m[i, j] for m in matrices2 if not np.isnan(m[i, j])]

            if len(vals1) >= 2 and len(vals2) >= 2:
                if method == 'ttest':
                    _, p = stats.ttest_ind(vals1, vals2)
                elif method == 'mannwhitneyu':
                    _, p = stats.mannwhitneyu(vals1, vals2, alternative='two-sided')
                pvalues[i, j] = p
                log2fc[i, j] = np.log2((np.mean(vals2) + 1) / (np.mean(vals1) + 1))

    return log2fc, pvalues

# Example with replicates
rep1_cond1 = [clr.matrix(balance=True).fetch(region) for clr in condition1_reps]
rep1_cond2 = [clr.matrix(balance=True).fetch(region) for clr in condition2_reps]

log2fc, pvalues = differential_test(rep1_cond1, rep1_cond2)

FDR Correction

from statsmodels.stats.multitest import multipletests

# Flatten p-values, apply FDR
pval_flat =