/bio-hi-c-analysis-loop-calling

name: bio-hi-c-analysis-loop-calling
description: Detect chromatin loops and point interactions from Hi-C data using cooltools, chromosight, and HiCCUPS-like methods. Identify CTCF-mediated loops and enhancer-promoter contacts. Use when detecting chromatin loops from Hi-C data.
tool_type: mixed
primary_tool: cooltools

Version Compatibility

Reference examples tested with: bedtools 2.31+, cooler 0.9+, cooltools 0.6+, matplotlib 3.8+, numpy 1.26+, pandas 2.2+, pybedtools 0.9+

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.

Chromatin Loop Calling

**"Call chromatin loops from my Hi-C data"** → Detect point enrichments in contact matrices representing CTCF-mediated loops and enhancer-promoter interactions.

• Python: `cooltools.dots()` or `chromosight detect --pattern=loops`

Detect chromatin loops and point interactions from Hi-C data.

Required Imports

import cooler
import cooltools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import bioframe

Call Loops with cooltools (Dot Calling)

clr = cooler.Cooler('matrix.mcool::resolutions/10000')
view_df = bioframe.make_viewframe(clr.chromsizes)

# Compute expected values
expected = cooltools.expected_cis(clr, view_df=view_df, ignore_diags=2)

# Call dots (loops)
dots = cooltools.dots(
    clr,
    expected=expected,
    view_df=view_df,
    max_loci_separation=2000000,  # Max loop size (2Mb)
    max_nans_tolerated=0.5,
)

print(f'Found {len(dots)} loops')
print(dots.head())

Using chromosight (CLI)

# Call loops with chromosight
chromosight detect \
    --pattern loops \
    --min-dist 20000 \
    --max-dist 2000000 \
    matrix.cool \
    loops_output

# Output: loops_output.tsv with loop coordinates and scores

Parse chromosight Output

loops = pd.read_csv('loops_output.tsv', sep='\t')
print(f'Found {len(loops)} loops')
print(loops.head())

# Columns: chrom1, start1, end1, chrom2, start2, end2, score, etc.

Using HiCExplorer hicDetectLoops

# Call loops with HiCExplorer
hicDetectLoops \
    -m matrix.cool \
    -o loops.bedgraph \
    --maxLoopDistance 2000000 \
    --windowSize 10 \
    --peakWidth 6 \
    --pValuePreselection 0.05 \
    --pValue 0.05

Loop Statistics

# Calculate loop sizes
loops['size'] = abs(loops['end2'] - loops['start1'])

print('Loop size statistics:')
print(f'  Mean: {loops["size"].mean() / 1000:.0f} kb')
print(f'  Median: {loops["size"].median() / 1000:.0f} kb')
print(f'  Min: {loops["size"].min() / 1000:.0f} kb')
print(f'  Max: {loops["size"].max() / 1000:.0f} kb')

# Size distribution
plt.hist(loops['size'] / 1000, bins=50)
plt.xlabel('Loop size (kb)')
plt.ylabel('Count')
plt.savefig('loop_sizes.png', dpi=150)

Filter Loops by Score

# Keep high-confidence loops
score_threshold = loops['score'].quantile(0.75)
high_conf_loops = loops[loops['score'] >= score_threshold]
print(f'High confidence loops: {len(high_conf_loops)}')

Annotate Loops with Features

import pybedtools

# Convert loop anchors to BED
anchor1 = loops[['chrom1', 'start1', 'end1']].copy()
anchor1.columns = ['chrom', 'start', 'end']
anchor2 = loops[['chrom2', 'start2', 'end2']].copy()
anchor2.columns = ['chrom', 'start', 'end']

# Load CTCF peaks
ctcf_peaks = pybedtools.BedTool('ctcf_peaks.bed')

# Intersect anchors with CTCF
anchor1_bed = pybedtools.BedTool.from_dataframe(anchor1)
anchor1_ctcf = anchor1_bed.intersect(ctcf_peaks, wa=True, u=True)

print(f'Anchors with CTCF: {len(anchor1_ctcf)} / {len(anchor1)}')

Compare Loops Between Conditions

# Load loops from two conditions
loops1 = pd.read_csv('condition1_loops.bedpe', sep='\t')
loops2 = pd.read_csv('condition2_loops.bedpe', sep='\t')

# Find overlapping loops
tolerance = 20000  # 20kb

def loops_overlap(l1, l2, tol):
    return (l1['chrom1'] == l2['chrom1'] and
            l1['chrom2'] == l2['chrom2'] and
            abs(l1['start1'] - l2['start1']) <= tol and
            abs(l1['start2'] - l2['start2']) <= tol)

shared = []
for _, loop1 in loops1.iterrows():
    for _, loop2 in loops2.iterrows():
        if loops_overlap(loop1, loop2, tolerance):
            shared.append(loop1)
            break

print(f'Shared loops: {len(shared)}')
print(f'Condition 1 specific: {len(loops1) - len(shared)}')
print(f'Condition 2 specific: {len(loops2) - len(set(range(len(loops2))) - set([]))}')

Aggregate Peak Analysis (APA)

**Goal:** Assess the overall strength and validity of called loops by stacking contact sub-matrices centered on loop anchors and averaging the signal.

**Approach:** For each loop, extract a fixed-size snippet from the contact matrix centered on the loop anchor pair, then compute the element-wise mean across all snippets to produce an aggregate enrichment map.

# Stack loops and compute average signal
from cooltools.lib import snip

def compute_apa(clr, loops, window=100000, resolution=10000):
    '''Compute average peak analysis'''
    flank = window // resolution

    stacks = []
    for _, loop in loops.iterrows():
        try:
            # Get region around loop
            snippet = clr.matrix(balance=True).fetch(
                f"{loop['chrom1']}:{loop['start1']-window}-{loop['end1']+window}",
                f"{loop['chrom2']}:{loop['start2']-window}-{loop['end2']+window}"
            )
            if snippet.shape[0] == snippet.shape[1]:
                stacks.append(snippet)
        except:
            continue

    if len(stacks) > 0:
        apa = np.nanmean(stacks, axis=0)
        return apa
    return None

apa_matrix = compute_apa(clr, loops.hea