# QC and Trimming for Hi-C Data

Hi-C reads contain chimeric sequences from ligation junctions. Ligation
junctions are handled downstream by pairtools.

**This workflow does not trim.** `main.nf` runs FastQC on the raw reads and
passes the same raw reads to `bwa mem -SP5M`, which soft-clips adapter
sequence. Trim Galore is present in the container image but no process calls
it. The trimming commands below are a manual, optional pre-processing step.

## Pre-Alignment QC with FastQC

```bash
fastqc --threads 4 --outdir fastqc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz
```

Key checks:
- Per-base quality (expect Phred >28 across most positions)
- Adapter content (Illumina adapters)
- Sequence length distribution (Hi-C reads are typically 50-150 bp)
- Duplication level (Hi-C libraries often show high duplication at low depths)

**Note**: Hi-C data will show unusual insert size distributions because reads
originate from ligation junctions, not contiguous fragments. This is expected.

## Optional: Adapter Trimming with Trim Galore (not run by this workflow)

If you choose to trim before running the workflow, keep it light: do NOT
aggressively trim quality, since chimeric reads may have lower quality at the
junction. Feed the resulting `*_val_{1,2}.fq.gz` files to `--reads`.

```bash
trim_galore \
    --paired \
    --quality 15 \
    --phred33 \
    --length 30 \
    --cores 4 \
    --fastqc \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

### Parameter Rationale

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--quality 15` | Phred 15 | Lenient -- chimeric reads have junction artifacts |
| `--length 30` | 30 bp | Short reads still carry valid contact information |
| `--cores 4` | 4 | Trim Galore uses ~3x threads internally |

### Default here: No Trimming

Like Juicer, this workflow skips trimming entirely and relies on the aligner
to handle adapter contamination via soft-clipping, which BWA-MEM does:

```bash
# No trim_galore step; raw FASTQs go straight to bwa mem -SP5M,
# which soft-clips adapter sequence.
```

## Restriction Enzyme Verification

Before processing, verify which restriction enzyme was used:

| Enzyme | Recognition Site | Ligation Junction | Average Fragment |
|--------|-----------------|-------------------|------------------|
| MboI/DpnII | GATC | GATCGATC | ~256 bp |
| HindIII | AAGCTT | AAGCTAGCTT | ~4 kb |
| Arima | Two sites | Multiple | ~160 bp |
| NcoI | CCATGG | CCATGCATGG | ~2 kb |

```bash
# Check for ligation junction sequence in reads (MboI example)
zcat sample_R1.fastq.gz | head -10000 | grep -c 'GATCGATC'
```

If the junction sequence appears frequently (>1% of reads), the enzyme
assignment is confirmed.

## Post-Trimming Summary (only if you trimmed manually)

After trimming, verify:
- >95% reads pass filters
- Median read length >40 bp
- Adapter detection rate is reasonable (typically 5-20% for Hi-C)

# BWA `-SP5M` Alignment for Hi-C

Hi-C reads are chimeric: each mate can originate from a different genomic
location due to proximity ligation. Standard paired-end alignment fails
because mates are not from a contiguous fragment. The solution is `bwa mem
-SP5M`, which skips pairing and mate rescue while still writing both mates to
one BAM.

## Genome Index Preparation

```bash
# Build BWA index (one-time step, ~1 hour for human genome)
bwa index -a bwtsw genome.fa
```

Requires ~8 GB disk and ~8 GB RAM for human genome.

## Alignment Strategy Used by This Workflow

Align both mates in a single `bwa mem -SP5M` call and write one BAM:

```bash
bwa mem -t 8 -SP5M genome.fa \
    sample_R1.fastq.gz sample_R2.fastq.gz \
    | samtools view -@ 4 -bhS - \
    > sample.paired.bam
```

`main.nf` publishes this file as `alignment/{sample}.paired.bam`. It keeps the
pair information in one file, which pairtools parses directly, and it is the
ENCODE-recommended approach.

### BWA-MEM Flags for Hi-C

| Flag | Meaning | Reason |
|------|---------|--------|
| `-S` | Skip mate rescue | Mates are on different chromosomes |
| `-P` | Skip pairing | Do not try to pair mates |
| `-5` | Split alignment: primary = 5' end | Ensures primary alignment is the 5' portion |
| `-M` | Mark shorter split as secondary | Compatible with downstream tools |

These flags are CRITICAL for Hi-C. Without `-SP`, BWA will try to pair
mates expecting a standard insert size, which fails for Hi-C contacts.

## Per-Mate Alignment (Alternative, not used here)

Some pipelines align R1 and R2 independently as single-end reads and merge the
two BAMs before parsing:

```bash
bwa mem -t 8 -SP5M genome.fa sample_R1.fastq.gz \
    | samtools view -@ 4 -bS - > sample_R1.bam

bwa mem -t 8 -SP5M genome.fa sample_R2.fastq.gz \
    | samtools view -@ 4 -bS - > sample_R2.bam
```

This workflow does not do this and never writes per-mate BAMs.

## Alignment QC (manual)

`main.nf` runs no flagstat on the Hi-C BAM; run it yourself on the published
file if you need it:

```bash
samtools flagstat alignment/sample.paired.bam
```

Expected metrics:
- **Mapped**: >80% (both mates combined)
- **Properly paired**: Low percentage is NORMAL for Hi-C (mates are on different chromosomes)
- **Supplementary alignments**: 5-20% (chimeric reads split across junction)

### Chimeric Read Handling

BWA-MEM produces supplementary alignments for reads spanning ligation junctions.
The `-5` flag ensures the 5' portion is reported as primary:

```
Read spans junction:
  [ChromA:100-200]---GATCGATC---[ChromB:500-600]

BWA reports:
  Primary:      ChromA:100-200  (5' end)
  Supplementary: ChromB:500-600  (3' end)
```

pairtools uses the primary alignment for contact calling.

## MAPQ Filtering

Filter low-quality alignments after pairing (in pairtools step), not here.
Keeping all alignments allows pairtools to classify pair types correctly.

Typical MAPQ threshold: 30 (applied during pairtools parse).

## Memory and Time Estimates

| Genome | Threads | RAM | Time (500M reads) |
|--------|---------|-----|---------------------|
| Human (hg38) | 8 | 12 GB | 3-5 hours |
| Mouse (mm10) | 8 | 10 GB | 2-4 hours |

BWA-MEM memory usage scales with index size, not read count. 8 GB is
sufficient for human genome alignment.

# Pair Processing with pairtools

pairtools is the ENCODE-standard tool for Hi-C pair classification,
filtering, and deduplication. It processes aligned BAMs into sorted,
deduplicated .pairs files ready for matrix generation.

## Parse: Classify Read Pairs

Convert aligned BAM to .pairs format with pair type classification:

```bash
# pairtools sort hands --tmpdir to GNU sort, which does not create it:
# the directory must already exist.
mkdir -p tmp

pairtools parse \
    --chroms-path chrom.sizes \
    --min-mapq 30 \
    --walks-policy mask \
    --max-inter-align-gap 30 \
    --nproc-in 4 \
    --nproc-out 4 \
    --output-stats sample.parse_stats.txt \
    sample.paired.bam \
    | pairtools sort \
        --nproc 4 \
        --tmpdir $PWD/tmp \
        -o sample_parsed_sorted.pairs.gz
```

`main.nf` publishes the statistics file as `pairs/{sample}.parse_stats.txt`.
It holds the pair-type breakdown for the whole library, which is the only
place those counts are available (see below).

### Parse Parameters

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--min-mapq 30` | MAPQ 30 | Filter low-confidence alignments |
| `--walks-policy mask` | Mask walks | Handle complex ligation events |
| `--max-inter-align-gap 30` | 30 bp | Maximum gap between split alignments |

### Pair Types Output

pairtools assigns each pair a two-letter code, one letter per side: U unique,
R rescued, M multi, N null (unmapped), W walk, D duplicate, X corrupt.

| Code | Meaning | Use |
|------|---------|-----|
| UU | Both uniquely mapped | Primary contacts |
| UR/RU | One unique, one rescued | Valid with caution |
| MU | One multi-mapped, one unique | Ambiguous, excluded |
| MM | Both multi-mapped | Excluded |
| NU | One unique, one unmapped | Excluded |
| NM | One unmapped, one multi-mapped | Excluded |
| NN | Both unmapped | Excluded |
| WW | Complex walk (multiple ligations), masked by `--walks-policy mask` | Excluded |
| DD | Duplicate | Removed in dedup step |
| XX | Corrupt record | Excluded |

## Sort Pairs

Pairs must be sorted by genomic position for deduplication:

```bash
pairtools sort \
    --nproc 4 \
    --tmpdir /tmp/ \
    sample_parsed.pairs.gz \
    -o sample_sorted.pairs.gz
```

## Deduplicate

Remove PCR/optical duplicates based on alignment positions:

```bash
pairtools dedup \
    --nproc-in 4 \
    --nproc-out 4 \
    --mark-dups \
    --output-stats sample.dedup_stats.txt \
    -o sample.dedup.pairs.gz \
    sample_sorted.pairs.gz
```

Both files are published: `pairs/{sample}.dedup.pairs.gz` and
`pairs/{sample}.dedup_stats.txt`.

### Dedup Statistics

The stats file reports:
- Total pairs processed
- Unique pairs retained
- PCR duplicate pairs removed
- Optical duplicate pairs removed
- Complexity estimate

Expected duplication rate: 10-40% depending on library complexity and depth.

## Filter for Valid Contacts

Select only UU pairs for contact matrix generation:

```bash
pairtools select \
    '(pair_type == "UU")' \
    sample.dedup.pairs.gz \
    -o sample_valid.pairs.gz
```

This is what `main.nf` does; the selected pairs file itself is an intermediate
and is not published.

For higher sensitivity (at cost of some noise), you can include rescued pairs
manually. The workflow does not offer this as an option:

```bash
pairtools select \
    '(pair_type == "UU") or (pair_type == "UR") or (pair_type == "RU")' \
    sample.dedup.pairs.gz \
    -o sample_valid_rescued.pairs.gz
```

## Pair Statistics

Generate detailed contact statistics:

```bash
pairtools stats \
    sample_valid.pairs.gz \
    -o sample_contact_stats.txt
```

`main.nf` runs exactly this on the UU-selected pairs and publishes it as
`qc/{sample}.contact_stats.txt`.

Key metrics from the stats output:
- **cis contacts**: Same chromosome
- **trans contacts**: Different chromosomes
- **cis >20kb**: Long-range cis contacts (biologically meaningful)
- **cis <20kb**: Short-range, often ligation artifacts

Because the input is UU-only, the pair-type distribution in this file is 100%
UU by construction. Read pair types from `pairs/{sample}.parse_stats.txt`
instead.

## Cis/Trans Ratio

The cis/trans ratio is a key QC metric. Match whole keys: `pairtools stats`
also emits `cis_1kb+` ... `cis_40kb+` rows, and a regex like `/cis/` would
pick up the last of those instead of the `cis` total.

```bash
awk '$1=="cis" {c=$2} $1=="trans" {t=$2} END {
    print "Cis:", c;
    print "Trans:", t;
    print "Cis/Trans ratio:", c/t
}' sample_contact_stats.txt
```

For thresholds, use the single QC table in `SKILL.md` (cis/trans >1.5 pass,
1.0-1.5 warning, <1.0 fail).

# Contact Matrix Generation with Juicer and Cooler

Generate multi-resolution contact matrices in .hic (Juicer) and .mcool
(cooler) formats from deduplicated .pairs files.

The "Matrix Quality Assessment" commands that use cooltools (distance decay,
A/B compartments) are manual steps: the workflow does not run them, and
cooltools is not installed in the container image. Use the conda environment
(`hic-env.yml`) or install cooltools separately for those.

## Juicer Tools: Generate .hic File

The .hic format is the ENCODE standard for Hi-C contact matrices.

### Convert pairs to Juicer medium format

```bash
# pairtools output -> Juicer medium format
# Format: strand1 chr1 pos1 frag1 strand2 chr2 pos2 frag2
pairtools select \
    '(pair_type == "UU")' \
    sample_dedup.pairs.gz \
    | awk 'BEGIN{OFS="\t"} !/^#/ {
        s1 = ($6 == "+") ? 0 : 16;
        s2 = ($7 == "+") ? 0 : 16;
        print s1, $2, $3, 0, s2, $4, $5, 1
    }' > sample_juicer_medium.txt
```

### Build .hic file

```bash
java -Xmx54g -jar juicer_tools.jar pre \
    --threads 4 \
    -r 1000,5000,10000,25000,50000,100000,250000,500000,1000000 \
    -k KR,VC,VC_SQRT \
    sample_juicer_medium.txt \
    sample.hic \
    hg38.chrom.sizes
```

### Key Parameters

| Parameter | Value | Reason |
|-----------|-------|--------|
| `-Xmx54g` | JVM heap | Large matrices need significant memory. `main.nf` gives the heap 85% of the task's memory allocation (54 GB on the first attempt of the 64 GB request), leaving the rest for the JVM itself |
| `--threads` | 4 | Threads for the normalization step only. Without `--mndindex`, `pre` builds the matrix single-threaded and says so on stderr ("Using single threaded preprocessor") |
| `-r` | Multiple resolutions | Enables multi-scale analysis (`--resolutions`) |
| `-k` | KR,VC,VC_SQRT | Generate multiple normalization vectors |

### Normalization Methods

| Method | Description | Use Case |
|--------|-------------|----------|
| KR | Knight-Ruiz balancing | Default, most common |
| VC | Vanilla Coverage | Simple, transparent |
| VC_SQRT | Square root VC | Reduces extreme values |

ENCODE standard: KR normalization for primary analysis.

## Cooler: Generate .mcool File

The .mcool format is preferred by Python/R analysis tools.

### Load pairs into cooler

cooler bins once at a base resolution and coarsens that into the others, so
`main.nf` uses the smallest value in `--resolutions` as the base bin (1000 in
the default list). Every other resolution has to be a multiple of it, which the
workflow checks before it starts.

```bash
# Create single-resolution .cool files, then zoom
cooler cload pairs \
    --chrom1 2 --pos1 3 --chrom2 4 --pos2 5 \
    --assembly hg38 \
    hg38.chrom.sizes:1000 \
    sample_valid.pairs.gz \
    sample_1kb.cool

# Generate multi-resolution mcool
cooler zoomify \
    --balance \
    --resolutions 1000,5000,10000,25000,50000,100000,250000,500000,1000000 \
    --nproc 4 \
    sample_1kb.cool \
    -o sample.mcool
```

### ICE Balancing

cooler uses ICE (Iterative Correction and Eigenvector decomposition) by
default via the `--balance` flag. This is equivalent to the `cooler balance`
command:

```bash
cooler balance --nproc 4 sample_1kb.cool
```

## Matrix Quality Assessment

### Contact Count per Resolution

```bash
# Check contacts at each resolution
for res in 1000 5000 10000 25000; do
    cooler info sample.mcool::resolutions/${res} | \
        python3 -c "import json,sys; d=json.load(sys.stdin); print(f'${res}bp: {d[\"sum\"]:,} contacts')"
done
```

### Distance Decay Curve (manual, not run by this workflow)

The contact frequency vs distance curve is a fundamental QC metric:

```bash
# Using cooltools (not in the container image)
cooltools expected-cis \
    --nproc 4 \
    sample.mcool::resolutions/10000 \
    -o sample_expected_cis.tsv
```

The slope should be approximately -1 on a log-log plot (power law decay).
Deviations indicate:
- Steeper slope: Over-digestion or poor ligation
- Shallower slope: Incomplete digestion

### Compartment Calling (A/B) -- manual, not run by this workflow

The workflow has no compartment step. To call A/B compartments yourself, run
eigenvector decomposition on the published `.mcool` at 100 kb resolution:

```bash
# cooltools is not in the container image
cooltools eigs-cis \
    --n-eigs 3 \
    --phasing-track /ref/gene_density.bedGraph \
    sample.mcool::resolutions/100000 \
    -o sample_compartments
```

The first eigenvector (E1) separates active (A) from inactive (B) compartments.
Use gene density as phasing track to orient E1 correctly (positive = A compartment).

## Visualization

### Juicebox (for .hic files)
```bash
# Launch Juicebox desktop application
java -jar Juicebox.jar sample.hic
```

### HiGlass (for .mcool files)
```bash
# Ingest into HiGlass server
higlass-manage ingest sample.mcool --assembly hg38
```

# Loop Calling with HiCCUPS

HiCCUPS (Hi-C Computational Unbiased Peak Search) identifies chromatin loops
from contact matrices. It is the ENCODE-standard loop caller, part of
Juicer tools.

## HiCCUPS Loop Calling

```bash
java -Xmx13g -jar juicer_tools.jar hiccups \
    --cpu \
    --threads 4 \
    -k KR \
    -r 5000,10000,25000 \
    -f 0.1,0.1,0.1 \
    -p 4,2,1 \
    -i 7,5,3 \
    -d 20000,20000,50000 \
    sample.hic \
    loops_output/
```

This is the command `main.nf` runs with the default `--hiccups_resolutions`
(`5000,10000,25000`). The heap is 85% of the task's memory allocation, so it is
13 GB on the first attempt of the 16 GB request in `nextflow.config` and grows
with each retry, rather than being a fixed number. `--cpu` is
required with the container image, which has no CUDA runtime; drop it only with
`--hiccups_gpu` on a host with an NVIDIA GPU. CPU mode restricts the search to a
band along the diagonal (8 Mb by default), so very long-range loops are not
reported.

### Key Parameters

`-r`, `-f`, `-p`, `-i` and `-d` take one value per resolution, in the order given to
`--hiccups_resolutions`; the peak widths, window widths and merge radii are Juicer's
published defaults for 5 kb, 10 kb and 25 kb. With `--hiccups_resolutions 10000` the
workflow runs `-r 10000 -f 0.1 -p 2 -i 5 -d 20000`.

| Parameter | Value | Meaning |
|-----------|-------|---------|
| `--cpu` | flag | CPU mode; needed without a CUDA runtime, searches 8 Mb from the diagonal |
| `-k` | KR | Normalization vector to read from the .hic file |
| `-r` | 5000,10000,25000 | Resolutions to search for loops (`--hiccups_resolutions`) |
| `-f` | 0.1,0.1,0.1 | FDR threshold per resolution |
| `-p` | 4,2,1 | Peak width (pixels) per resolution |
| `-i` | 7,5,3 | Window width (pixels) of the local background region per resolution |
| `-d` | 20000,20000,50000 | Merge radius (bp) around a loop centroid per resolution: 20 kb at 5 kb and 10 kb, 50 kb at 25 kb. juicer_tools 2.20.00 reads one value per `-r` resolution (`HiCCUPSConfiguration.extractIntegerValues(..., resolutions.length)`); its usage text still says "three values", but a list of any other length (except a single value, which is applied to every resolution) stops HiCCUPS with "Must pass N parameters" and exit code 30 |

`juicer_tools pre` must have written the `-k` vector into the .hic file:
`main.nf` builds KR, VC and VC_SQRT.

### Resolution Selection

| Resolution | Detects | Minimum Contacts |
|------------|---------|------------------|
| 5 kb | Fine-scale loops | >1 billion |
| 10 kb | Standard loops | >500 million |
| 25 kb | Large-scale loops | >100 million |

## HiCCUPS Output Format

HiCCUPS produces a BEDPE-like file with loop anchors. juicer_tools 2.20.00
writes the BEDPE core columns first, then the loop attributes in alphabetical
order, on a header line that starts with `#`:

```
#chr1  x1  x2  chr2  y1  y2  name  score  strand1  strand2  color  centroid1  centroid2  expectedBL  expectedDonut  expectedH  expectedV  fdrBL  fdrDonut  fdrH  fdrV  numCollapsed  observed  radius
```

Key columns (1-based, in that order):
- `chr1 x1 x2` (1-3) -- Upstream anchor
- `chr2 y1 y2` (4-6) -- Downstream anchor
- `centroid1`, `centroid2` (12-13) -- Centroid of the merged pixel cluster on each side
- `expectedDonut` (15) -- Expected count from the donut background model
- `fdrDonut` (19) -- FDR from the donut model (primary significance)
- `numCollapsed` (22) -- Enriched pixels merged into this call
- `observed` (23) -- Observed contact count
- `radius` (24) -- Radius of the merged cluster

The header line is the only `#` line; skip it before any arithmetic on the
file.

## Merge Loops Across Resolutions

HiCCUPS calls loops at each resolution independently and merges them itself:
it writes `merged_loops.bedpe` into the output directory alongside the
per-resolution files. There is no separate merge command (`juicer_tools` has
no `hiccups_merge` tool).

`main.nf` copies that merged file to `loops/{sample}.hiccups_loops.bedpe`,
which is the only loop file published; the per-resolution files stay in the
Nextflow work directory.

## Alternative: Mustache Loop Caller (manual, not in the image)

Mustache (Roayaei Ardakany 2020, ~165 citations) uses a scale-space
representation for loop detection. It is not installed in the container
image and the workflow does not run it:

```bash
mustache \
    -f sample.mcool \
    -r 10000 \
    -ch hg38.chrom.sizes \
    -o mustache_loops.bedpe \
    -pt 0.05 \
    -st 0.8
```

### HiCCUPS vs Mustache

| Feature | HiCCUPS | Mustache |
|---------|---------|----------|
| Background model | Donut + 3 others | Scale-space Gaussian |
| GPU support | Yes (CUDA; the image ships CPU mode only) | No |
| Speed | Faster with GPU | Moderate |
| Sensitivity | Standard | Higher (more loops) |
| ENCODE standard | Yes | Alternative |
| Concordance | ~50% overlap between callers (Wolff 2022) |

## Loop QC Metrics

### Loop Count by Resolution

HiCCUPS writes two files per resolution: `enriched_pixels_<res>.bedpe`
(pre-filter candidates) and `postprocessed_pixels_<res>.bedpe` (the final
calls that go into `merged_loops.bedpe`). Count the post-filter file:

```bash
# One file per resolution in --hiccups_resolutions (default: 5000 10000 25000)
for res in 5000 10000 25000; do
    count=$(grep -vc '^#' loops_output/postprocessed_pixels_${res}.bedpe)
    echo "Resolution ${res}: ${count} loops"
done
```

These per-resolution files are not published by the workflow; they remain in
the task's work directory. From a finished run, count the merged file instead,
excluding its header line:

```bash
grep -vc '^#' results/loops/sample.hiccups_loops.bedpe
```

Expected loop counts (human cell line, >1B contacts):
- 5 kb: 5,000-15,000 loops
- 10 kb: 3,000-10,000 loops
- 25 kb: 1,000-5,000 loops

### Loop Size Distribution

```bash
# Loop size = y1 - x1, i.e. column 5 minus column 2; skip the header line
awk '!/^#/ {print $5 - $2}' results/loops/sample.hiccups_loops.bedpe | \
    sort -n | \
    awk '{a[NR]=$1} END {
        print "Median loop size:", a[int(NR/2)];
        print "Min:", a[1];
        print "Max:", a[NR]
    }'
```

Typical loop sizes:
- Median: 200-400 kb
- Range: 50 kb to 5 Mb
- Loops <50 kb may be artifacts at lower resolutions

### CTCF Enrichment at Anchors (manual, not run by this workflow)

True loops are enriched for CTCF binding at anchors. bedtools is not in the
container image, so run this outside the workflow:

```bash
bedtools intersect \
    -a loop_anchors.bed \
    -b CTCF_peaks.bed \
    -u | wc -l
```

Expect >60% of loop anchors to overlap CTCF peaks for convergent CTCF loops.

# Hi-C Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-hic skill — papers defining ENCODE Hi-C processing standards, chromatin conformation capture methodology, contact matrix generation, loop calling, and 3D genome organization tools.

---

## Hi-C Method Development

---

### Lieberman-Aiden et al. 2009 — Comprehensive mapping of long-range interactions

- **Citation:** Lieberman-Aiden E, van Berkum NL, Williams L, Imakaev M, Raber T, Lajoie BR, Dekker J, Bhatt DM, Nussbaum MC, Anton IM, Ahmed K, Gnirke A, Lander ES. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science, 326(5950):289-293, 2009.
- **DOI:** [10.1126/science.1181369](https://doi.org/10.1126/science.1181369)
- **PMID:** 19815776 | **PMC:** PMC2858594
- **Citations:** ~6,000
- **Key findings:** Introduced the Hi-C method, the first technique to map genome-wide chromatin interactions in an unbiased manner. Discovered that the human genome is organized into two spatial compartments: A compartments (open, gene-rich, active chromatin) and B compartments (closed, gene-poor, inactive chromatin), visible as a plaid/checkerboard pattern in contact matrices at megabase resolution. Demonstrated that chromatin conformation is consistent with a "fractal globule" model where the genome folds without knots, enabling maximum accessibility. This foundational paper established Hi-C as the standard method for studying 3D genome organization.

---

### Rao et al. 2014 — A 3D map of the human genome at kilobase resolution

- **Citation:** Rao SSP, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, Sanborn AL, Machol I, Omer AD, Lander ES, Aiden EL. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell, 159(7):1665-1680, 2014.
- **DOI:** [10.1016/j.cell.2014.11.021](https://doi.org/10.1016/j.cell.2014.11.021)
- **PMID:** 25497547 | **PMC:** PMC5635824
- **Citations:** ~5,000
- **Key findings:** Generated the highest-resolution Hi-C maps to date using the in situ Hi-C protocol, which performs ligation inside intact nuclei rather than in dilute solution. Identified ~10,000 chromatin loops in human cells, the majority of which are anchored by convergent CTCF/cohesin binding sites (the "convergent rule" of loop formation). Discovered six subcompartments (A1, A2, B1, B2, B3, B4) refining the original A/B compartment classification. Introduced HiCCUPS (Hi-C Computational Unbiased Peak Search) for systematic loop detection and Arrowhead for TAD/domain annotation. This paper established the in situ Hi-C protocol and loop calling methodology used in the ENCODE pipeline.

---

## ENCODE Pipeline Standards

---

### ENCODE Project Consortium 2020 — Expanded encyclopaedias of DNA elements

- **Citation:** ENCODE Project Consortium et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature, 583(7818):699-710, 2020.
- **DOI:** [10.1038/s41586-020-2493-4](https://doi.org/10.1038/s41586-020-2493-4)
- **PMID:** 32728249 | **PMC:** PMC7410828
- **Citations:** ~2,500
- **Key findings:** ENCODE Phase 3 paper establishing Hi-C as a core 3D genome assay. Hi-C contact matrices are used to identify topologically associating domains (TADs), chromatin loops, and A/B compartments that organize gene regulation. Established quality standards: cis/trans ratio >1.5, long-range cis contacts (>20kb) >40%, and sufficient depth for the target resolution. Assembly standard: GRCh38/mm10.

---

### Hitz et al. 2023 — The ENCODE Uniform Analysis Pipelines

- **Citation:** Hitz BC, Lee JW, Jolanki O, et al. The ENCODE Uniform Analysis Pipelines. Research Square (preprint), 2023.
- **DOI:** [10.21203/rs.3.rs-311127/v1](https://doi.org/10.21203/rs.3.rs-311127/v1)
- **Citations:** ~84
- **Key findings:** Documents the official ENCODE Hi-C pipeline specification: BWA-MEM with `-SP5M` for chimeric Hi-C alignment, pairtools for pair classification and deduplication, Juicer tools for .hic file generation with KR normalization, cooler for .mcool generation, and HiCCUPS for loop calling. Specifies MAPQ >= 30 filtering, restriction site-aware pair classification, and multi-resolution matrix generation. Pipeline infrastructure uses Docker containers and WDL for reproducible execution.

---

### Yardimci et al. 2019 — Measuring the reproducibility and quality of Hi-C data

- **Citation:** Yardimci GG, Ozadam H, Sauria MEG, Ursu O, Yan KK, Yang T, Chakraborty A, Kaul A, Lajoie BR, Song F, Zhan Y, Ay F, Gerstein M, Kundaje A, Li Q, Taylor J, Yue F, Dekker J, Noble WS. Measuring the reproducibility and quality of Hi-C data. Genome Biology, 20:57, 2019.
- **DOI:** [10.1186/s13059-019-1658-7](https://doi.org/10.1186/s13059-019-1658-7)
- **PMID:** 30890173 | **PMC:** PMC6425651
- **Citations:** ~250
- **Key findings:** Systematic evaluation of Hi-C data quality metrics and reproducibility. Established that cis/trans ratio, long-range cis fraction, and stratum-adjusted correlation coefficient (SCC) are the most informative quality metrics. Demonstrated that loop calling requires substantially deeper sequencing than compartment or TAD analysis. Recommended minimum depths for different analyses: 10M contacts for compartments, 50M for TADs, 500M+ for loops at 5-10kb resolution. These thresholds inform the QC standards in this pipeline.

---

## Core Pipeline Tools

---

### Durand et al. 2016 — Juicer: one-click Hi-C analysis

- **Citation:** Durand NC, Shamim MS, Machol I, Rao SSP, Huntley MH, Lander ES, Aiden EL. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Systems, 3(1):95-98, 2016.
- **DOI:** [10.1016/j.cels.2016.07.002](https://doi.org/10.1016/j.cels.2016.07.002)
- **PMID:** 27467249 | **PMC:** PMC5846465
- **Citations:** ~2,000
- **Key findings:** Introduced the Juicer pipeline and juicer_tools for Hi-C data processing. The ENCODE pipeline uses juicer_tools pre for generating .hic files from sorted pairs, with Knight-Ruiz (KR) balanced normalization as the default. Juicer_tools also provides HiCCUPS for loop calling and Arrowhead for domain annotation. The .hic format stores multi-resolution contact matrices with embedded normalization vectors, enabling interactive exploration in Juicebox.

---

### Abdennur & Mirny 2020 — Cooler: scalable storage for Hi-C data

- **Citation:** Abdennur N, Mirny LA. Cooler: scalable storage for Hi-C data and other genomically labeled arrays. Bioinformatics, 36(1):311-316, 2020.
- **DOI:** [10.1093/bioinformatics/btz540](https://doi.org/10.1093/bioinformatics/btz540)
- **PMID:** 31290943
- **Citations:** ~300
- **Key findings:** Introduced the cooler format (.cool/.mcool), an HDF5-based storage format for Hi-C contact matrices that supports arbitrary genomic resolution and out-of-core operations on matrices too large for memory. The .mcool format stores multiple resolutions in a single file. Cooler integrates with cooltools for compartment analysis, insulation score calculation, and other downstream analyses. Used in this pipeline alongside Juicer to generate .mcool files as an alternative to .hic, with broader support in the Python scientific computing ecosystem.

---

### Open2C — pairtools: pair-level Hi-C data processing

- **Citation:** Open2C, Abdennur N, Abraham S, Fudenberg G, Flyamer IM, Galitsyna AA, Goloborodko A, Imakaev M, Oksuz BA, Venev SV. Pairtools: from sequencing data to chromosome-level 3D maps. PLoS Computational Biology, 20(5):e1012164, 2024.
- **DOI:** [10.1371/journal.pcbi.1012164](https://doi.org/10.1371/journal.pcbi.1012164)
- **PMID:** 38753885
- **Citations:** ~50
- **Key findings:** Pairtools provides a suite of command-line tools for processing Hi-C read pairs: parse (classify pairs from SAM), sort, dedup (remove PCR duplicates from pairs), select (filter by pair type), and stats (generate QC statistics). The ENCODE pipeline uses pairtools parse to classify read pairs into UU (both uniquely mapped), UR (one rescued), WW (walk/same-strand artifact), and other categories, followed by pairtools dedup to remove PCR duplicates. The .pairs format output is the standard intermediate format for Hi-C data before matrix generation.

---

### Imakaev et al. 2012 — Iterative correction of Hi-C data

- **Citation:** Imakaev M, Fudenberg G, McCord RP, Naumova N, Goloborodko A, Lajoie BR, Dekker J, Mirny LA. Iterative correction of Hi-C data reveals hallmarks of chromosome organization. Nature Methods, 9(10):999-1003, 2012.
- **DOI:** [10.1038/nmeth.2148](https://doi.org/10.1038/nmeth.2148)
- **PMID:** 22941365 | **PMC:** PMC3816492
- **Citations:** ~1,500
- **Key findings:** Introduced ICE (Iterative Correction and Eigenvector decomposition), a matrix balancing algorithm for removing systematic biases in Hi-C data including GC content, mappability, and restriction fragment length. ICE is the default normalization in the cooler/cooltools ecosystem (alternative to KR normalization in Juicer). Also introduced eigenvector decomposition of the normalized contact matrix for A/B compartment identification. The pipeline supports both KR (Juicer default) and ICE (cooler default) normalization methods.

---

## Supplementary Tools (Non-Hi-C-Specific)

See pipeline-chipseq/references/literature.md for detailed descriptions of shared tools.

---

### Li & Durbin 2009 — BWA: Burrows-Wheeler Aligner

- **DOI:** [10.1093/bioinformatics/btp324](https://doi.org/10.1093/bioinformatics/btp324) | **PMID:** 19451168 | **Citations:** ~25,000
- **Hi-C role:** Per-mate alignment of Hi-C reads. Unlike standard paired-end alignment, Hi-C reads are aligned independently because mates can map to distant genomic locations (chimeric ligation products). BWA-MEM is used with -SP flags to prevent rescue of discordant pairs. MAPQ >= 30 filtering removes ambiguously mapped reads.

---

### Li et al. 2009 — SAMtools

- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352) | **PMID:** 19505943 | **Citations:** ~53,700
- **Hi-C role:** BAM operations: converting the `bwa mem -SP5M` stream to the BAM that pairtools parse reads for pair classification.

---

### Ewels et al. 2016 — MultiQC

- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354) | **PMID:** 27312411 | **Citations:** ~6,800
- **Hi-C role:** Aggregates QC metrics from FastQC, BWA alignment, pairtools stats, and contact statistics into a unified HTML report. Provides batch-level assessment of pair type distributions and library quality.

---

### Amemiya et al. 2019 — ENCODE Blacklist

- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z) | **PMID:** 31249361 | **Citations:** ~1,372
- **Hi-C role:** Blacklist v2 filtering removes artifact-prone genomic regions from Hi-C contact matrices. Blacklisted regions can create spurious high-contact areas in matrices due to multi-mapping artifacts, affecting loop calling and compartment analysis.

---

### Andrews 2010 — FastQC

- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **Hi-C role:** Raw read quality assessment. Hi-C libraries typically show broader insert size distributions than standard WGS libraries and may have elevated adapter contamination rates from short ligation products.