ENCODE Hi-C Pipeline: FASTQ to Contact Matrices and Loops

When to Use

• User wants to run a Hi-C processing pipeline from FASTQ to contact matrices and loop calls
• User asks about "Hi-C pipeline", "contact matrix", "loop calling", "Juicer", "HiCCUPS", or "TAD detection"
• User needs to process Hi-C data for 3D genome structure analysis
• Example queries: "process my Hi-C FASTQs", "generate contact matrices from Hi-C", "call chromatin loops with HiCCUPS"

Execute the ENCODE Hi-C pipeline for chromatin conformation capture data, producing multi-resolution contact matrices and loop calls.

Pipeline Overview

FASTQ -> FastQC (raw reads)
      -> bwa mem -SP5M (both mates in one call) -> {sample}.paired.bam
         -> pairtools parse -> sort -> dedup -> select UU
              |
              +-> Juicer pre -> .hic -> HiCCUPS -> loops (BEDPE)
              |
              +-> cooler cload + zoomify -> .mcool

Not run by this workflow

Adapter trimming, TAD calling, A/B compartment calling, and every cooltools analysis are outside this workflow. They are documented as manual, optional steps only: trimming in references/01-qc-trimming.md, distance decay and compartments in references/04-matrix-generation.md. cooltools and bedtools are not installed in the container image, so those manual commands need the conda environment (environments/hic-env.yml in the bioinformatics-installer skill) or a separate install.

ENCODE Repository

• GitHub: ENCODE-DCC/hic-pipeline
• Container: built from scripts/Dockerfile in this skill (docker build -t encode-toolkit/pipeline-hic:1.0.0 scripts/); override with --container
• WDL: Available for Cromwell execution
• This skill: Nextflow DSL2 reimplementation for portability

Core Tools and Versions

Versions are those installed by scripts/Dockerfile, which is what the workflow runs.

Tool	Version	Purpose	Citation
BWA-MEM	0.7.18	Alignment (both mates, -SP5M)	Li & Durbin 2009
pairtools	1.1.2	Pair classification, dedup	Open2C
Juicer tools	2.20.00	.hic generation, HiCCUPS	Durand et al. 2016
cooler	0.9.3	.cool/.mcool generation	Abdennur & Mirny 2020
samtools	1.19	BAM operations	Li et al. 2009
FastQC	0.12.1	Read quality	Andrews (Babraham)
MultiQC	1.21	Aggregated QC	Ewels et al. 2016

The conda alternative (hic-env.yml) ships no juicer_tools jar, only a JRE, so .hic generation and HiCCUPS are unavailable on that route; conversely it provides cooltools and bedtools, which the container image does not.

Key Literature

1. Rao et al. 2014 - "A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping" (Cell, ~5,000 citations) DOI: 10.1016/j.cell.2014.11.021

2. Lieberman-Aiden et al. 2009 - "Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome" (Science, ~6,000 citations) DOI: 10.1126/science.1181369

3. Durand et al. 2016 - "Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments" (Cell Systems, ~2,000 citations) DOI: 10.1016/j.cels.2016.07.002

4. Abdennur & Mirny 2020 - "Cooler: scalable storage for Hi-C data and other genomically labeled arrays" (Bioinformatics) DOI: 10.1093/bioinformatics/btz540

5. Amemiya et al. 2019 - "The ENCODE Blacklist" (Scientific Reports, ~1,372 citations) DOI: 10.1038/s41598-019-45839-z

Execution

Quick Start (Local)

nextflow run scripts/main.nf \
    -profile local \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index '/ref/bwa_index/GRCh38.fa' \
    --chrom_sizes '/ref/hg38.chrom.sizes' \
    --outdir results/ \
    -resume

SLURM HPC

nextflow run scripts/main.nf \
    -profile slurm \
    --container /path/to/pipeline-hic.sif \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index '/ref/bwa_index/GRCh38.fa' \
    --chrom_sizes '/ref/hg38.chrom.sizes' \
    --outdir results/ \
    -resume

Cloud (GCP / AWS)

# Google Cloud Batch
nextflow run scripts/main.nf -profile gcp \
    --container us-docker.pkg.dev/<project>/<repo>/pipeline-hic:1.0.0 \
    --gcp_project <project> \
    --gcp_workdir gs://<bucket>/work \
    --reads 'gs://<bucket>/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index gs://<bucket>/ref/GRCh38.fa \
    --chrom_sizes gs://<bucket>/ref/hg38.chrom.sizes \
    --outdir gs://<bucket>/results

# AWS Batch
nextflow run scripts/main.nf -profile aws \
    --container <account>.dkr.ecr.<region>.amazonaws.com/pipeline-hic:1.0.0 \
    --aws_queue <job-queue> \
    --aws_workdir s3://<bucket>/work \
    --reads 's3://<bucket>/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index s3://<bucket>/ref/GRCh38.fa \
    --chrom_sizes s3://<bucket>/ref/hg38.chrom.sizes \
    --outdir s3://<bucket>/results

--outdir only sets where results are published; Google Batch and AWS Batch stage every task through the work directory, and the workflow stops with an error if it or the project/queue is missing.

Resource Requirements

Step	CPUs	RAM	Time (2B contacts)
BWA alignment	8	16 GB	4-6 hours
pairtools parse + sort	4	16 GB	2-3 hours
pairtools dedup	4	16 GB	1-2 hours
Juicer pre + hic	4	64 GB	2-4 hours
HiCCUPS	4	16 GB	1-2 hours
Total	8	64 GB	8-16 hours

The RAM column is each step's first-attempt request. Every process asks for that much memory per attempt, so a task killed for exceeding it is retried with more (at most two retries, capped by --max_memory). Failures with any other exit status stop the run.

Pipeline Parameters

Parameter	Default	Description
--reads	required	Glob pattern to paired FASTQ files
--bwa_index	required	BWA index prefix: the genome FASTA path whose .amb .ann .bwt .pac .sa files sit beside it (every file starting with this prefix is staged)
--chrom_sizes	required	Chromosome sizes file
--outdir	./results	Output directory
--resolutions	1000,5000,10000,25000,50000,100000,250000,500000,1000000	Matrix resolutions for juicer_tools pre and cooler zoomify: a comma-separated list of positive integers. The smallest value is the cooler base bin, and every other value must be a multiple of it; the workflow stops with an error otherwise
--hiccups_resolutions	5000,10000,25000	Resolutions HiCCUPS calls loops at. Only 5000, 10000 and 25000 are accepted, and each must also be listed in --resolutions; the workflow stops with an error otherwise. Peak width (-p), window width (-i), merge radius (-d) and FDR (-f) follow Juicer's published per-resolution defaults, one value per resolution
--min_mapq	30	Minimum MAPQ passed to pairtools parse
--hiccups_gpu	false	Run HiCCUPS on an NVIDIA GPU. By default the CPU mode is used, which only searches within 8 Mb of the diagonal
--assembly	hg38	Assembly name recorded in the .mcool metadata (cooler cload --assembly). The .hic file is built from the chrom.sizes file only

Infrastructure parameters (nextflow.config)

Parameter	Default	Description
--container	encode-toolkit/pipeline-hic:1.0.0	Image built from scripts/Dockerfile. Pass a registry image for gcp/aws, or a .sif file for slurm
--max_cpus, --max_memory, --max_time	16, 128.GB, 48.h	Upper bounds applied to every process
--slurm_queue, --slurm_account	normal, none	SLURM partition and account
--gcp_project, --gcp_workdir	none (both required for -profile gcp)	Google Cloud project and gs:// work directory
--gcp_location, --gcp_disk	us-central1, 200.GB	Google Batch region and per-task disk
--aws_queue, --aws_workdir	none (both required for -profile aws)	AWS Batch job queue and s3:// work directory
--aws_region, --aws_cli_path	us-east-1, /home/ec2-user/miniconda/bin/aws	AWS region, and the AWS CLI path inside the Batch AMI

Output Files

results/
  fastqc/
    *_fastqc.html                 # Raw read quality
    *_fastqc.zip
  alignment/
    {sample}.paired.bam           # Both mates from one bwa mem -SP5M call
  pairs/
    {sample}.parse_stats.txt      # pairtools parse stats (pair-type breakdown)
    {sample}.dedup.pairs.gz       # Classified, sorted, deduplicated pairs
    {sample}.dedup_stats.txt      # pairtools dedup stats (duplication, complexity)
  matrices/
    {sample}.hic                  # Juicer .hic file (primary output)
    {sample}.mcool                # Cooler multi-resolution matrix
  loops/
    {sample}.hiccups_loops.bedpe  # HiCCUPS merged_loops.bedpe, renamed
  qc/
    {sample}.contact_stats.txt    # pairtools stats on the selected UU pairs
  multiqc/
    multiqc_report.html
  pipeline_info/
    timeline.html
    report.html
    trace.txt

The UU-selected pairs file and the pairtools stats run on it are intermediate: only qc/{sample}.contact_stats.txt is published, not the selected pairs.

.hic File Format

The .hic format (Juicer) stores multi-resolution contact matrices with normalization vectors. Can be visualized in Juicebox and loaded by hic-straw in Python/R.

.mcool File Format

The .mcool format (cooler) is an HDF5-based multi-resolution contact matrix. Widely supported by cooler, cooltools, HiGlass, and FAN-C.

QC Thresholds (ENCODE Standards)

This is the only QC threshold table for this skill; the reference files point back to it.

Metric	Pass	Warning	Fail	Computed from
Valid (UU) pair fraction	>40%	25-40%	<25%	pairs/{sample}.parse_stats.txt
Cis contacts (>20kb)	>40%	25-40%	<25%	qc/{sample}.contact_stats.txt
Cis/trans ratio	>1.5	1.0-1.5	<1.0	qc/{sample}.contact_stats.txt
Library complexity (unique/total)	>0.7	0.5-0.7	<0.5	pairs/{sample}.dedup_stats.txt

qc/{sample}.contact_stats.txt is computed after UU selection, so its pair-type breakdown is 100% UU by construction. Read pair types from pairs/{sample}.parse_stats.txt instead.

Resolution vs Depth Requirements

Resolution	Minimum Contacts Needed	Typical Depth
1 kb	>2 billion	Very deep
5 kb	>500 million	Deep
10 kb	>200 million	Standard
25 kb	>50 million	Moderate
100 kb	>10 million	Low

Pair Classification

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

Category	Description	Use
UU	Both sides uniquely mapped	Valid contact -- the only type this workflow keeps
UR / RU	One unique, one rescued	Valid but not selected here
NU	One unique, one unmapped	Not used
NM	One unmapped, one multi-mapped	Not used
MM	Both multi-mapped	Not used
WW	Complex walk (multiple ligation events), masked by --walks-policy mask	Not used
DD	Duplicate	Removed by pairtools dedup
XX	Corrupt record	Not used

This workflow selects UU only (pairtools select '(pair_type == "UU")') before matrix generation.

Critical Pitfalls

Restriction Enzyme Choice

The restriction enzyme determines fragment size and resolution:

• MboI/DpnII (GATC): 4-cutter, ~256 bp average fragment -- higher resolution
• HindIII (AAGCTT): 6-cutter, ~4 kb average fragment -- lower resolution
• Arima (proprietary): Two enzymes, ~160 bp average -- highest resolution
• Always verify which enzyme was used before interpreting resolution
• The workflow itself is enzyme-agnostic: pairtools works at read-pair level and the .hic file is built without a restriction-site file, so there is no enzyme parameter to set

Normalization Method

Different normalization methods yield different results:

• KR (Knight-Ruiz): built by juicer_tools pre -k KR,VC,VC_SQRT and used by HiCCUPS (-k KR)
• ICE (Imakaev et al.): applied to the .mcool by cooler zoomify --balance
• VC (Vanilla Coverage): simple coverage normalization, also built into the .hic
• Always document which normalization a downstream analysis read.

Resolution Depends on Depth

Do not call features at resolutions unsupported by sequencing depth:

• Calling 1 kb loops from 100M contacts will produce noise
• Check the Juicer resolution QC to determine achievable resolution
• HiCCUPS runs at the resolutions in --hiccups_resolutions (default 5 kb, 10 kb and 25 kb); each of them must also be listed in --resolutions, because HiCCUPS reads them out of the .hic file

Ligation Artifacts

Monitor the pair-type breakdown in pairs/{sample}.parse_stats.txt:

• WW pairs are complex walks (more than one ligation in a read); --walks-policy mask masks them so they never reach the contact matrix
• A large unmapped/multi-mapped fraction points at poor library or the wrong genome
• The workflow produces no re-ligation distance plot; derive one manually from the published pairs file if needed

Provenance Integration

After pipeline completion, log all outputs:

encode_log_derived_file(
    file_path="/results/matrices/sample1.hic",
    source_accessions=["ENCSR...", "ENCFF..."],
    description="Hi-C contact matrix from ENCODE Hi-C pipeline",
    file_type="hic",
    tool_used="BWA 0.7.18 + pairtools 1.1.2 + Juicer 2.20.00",
    parameters="--min_mapq 30, UU pairs only, KR/VC/VC_SQRT normalization, resolutions 1kb-1Mb"
)

Reference Files

Detailed step-by-step documentation is provided in the references/ directory:

1. 01-qc-trimming.md -- Read QC (trimming is a manual option, not run here)
2. 02-alignment.md -- BWA -SP5M alignment of both mates in one call
3. 03-pair-processing.md -- pairtools parse, sort, dedup, and select
4. 04-matrix-generation.md -- Juicer .hic and cooler .mcool generation; manual cooltools analyses
5. 05-loop-calling.md -- HiCCUPS loop detection and QC

Walkthrough: Processing ENCODE Hi-C from FASTQ to Contact Maps and Loops

Goal: Process raw Hi-C FASTQ files through the ENCODE pipeline to generate contact matrices and chromatin loop calls. Context: Hi-C captures 3D chromatin organization. The pipeline uses BWA for chimeric read alignment, pairtools for pair processing, and Juicer/HiCCUPS for loop calling.

Step 1: Find Hi-C experiment

encode_get_experiment(accession="ENCSR000AKA")

Expected output:

{
  "accession": "ENCSR000AKA",
  "assay_title": "Hi-C",
  "biosample_summary": "GM12878",
  "bio_replicate_count": 2,
  "status": "released"
}

Step 2: List FASTQ files

encode_list_files(experiment_accession="ENCSR000AKA", file_format="fastq")

Expected output (a JSON array of file records; fields abridged):

[
  {"accession": "ENCFF500HI1", "file_format": "fastq", "output_type": "reads", "biological_replicates": [1], "file_size_human": "34.2 GB", "status": "released"},
  {"accession": "ENCFF501HI2", "file_format": "fastq", "output_type": "reads", "biological_replicates": [1], "file_size_human": "35.2 GB", "status": "released"}
]

Interpretation: Hi-C paired-end reads represent chimeric ligation junctions. Each read pair captures a 3D contact.

Step 3: Name the files so a read-pair glob can find them

ENCODE FASTQs are named by accession, so the two mates of a pair share no prefix, and the workflow matches file pairs with a {1,2} glob. Which mate a file is comes from its page on encodeproject.org (paired_end 1 or 2, and paired_with naming the other accession), not from any tool here. Link the files into the shape the glob expects:

mkdir -p fastq
ln -s "$PWD/ENCFF500HI1.fastq.gz" fastq/ENCSR000AKA_R1.fastq.gz
ln -s "$PWD/ENCFF501HI2.fastq.gz" fastq/ENCSR000AKA_R2.fastq.gz

Step 4: Run the Hi-C pipeline

nextflow run scripts/main.nf \
  -profile local \
  --reads 'fastq/ENCSR000AKA_R{1,2}.fastq.gz' \
  --bwa_index '/ref/bwa_index/GRCh38.fa' \
  --chrom_sizes '/ref/hg38.chrom.sizes' \
  --outdir results/ \
  -resume

Key pipeline steps:

1. FastQC on the raw reads
2. BWA-MEM -SP5M alignment of both mates in one call (chimeric read handling)
3. pairtools parse + sort (classify pairs, MAPQ 30, mask walks)
4. pairtools dedup (remove PCR duplicates), then select UU pairs
5. Contact matrix generation (.hic via Juicer, .mcool via cooler)
6. Loop calling (HiCCUPS at the --hiccups_resolutions, by default 5 kb, 10 kb and 25 kb, merged into one BEDPE)

Step 5: Validate output quality

Use the QC threshold table above with pairs/{sample}.parse_stats.txt, pairs/{sample}.dedup_stats.txt and qc/{sample}.contact_stats.txt.

Step 6: Identify significant loops

Download loop calls for downstream analysis:

encode_list_files(experiment_accession="ENCSR000AKA", file_format="bedpe", assembly="GRCh38")

Integration with downstream skills

• Loop calls (BEDPE) feed into -> hic-aggregation for cross-tissue loop catalog
• Loop anchors feed into -> peak-annotation for enhancer-promoter assignment
• Contact data integrates with -> visualization-workflow for 3D genome display
• Pipeline provenance logged by -> data-provenance

Code Examples

1. Find Hi-C data for 3D genome analysis

encode_search_experiments(
  assay_title="Hi-C",
  organ="heart"
)

Expected output:

{
  "results": [
    {
      "accession": "ENCSR654HRT",
      "assay_title": "Hi-C",
      "biosample_summary": "heart left ventricle tissue male adult (51 years)",
      "status": "released"
    }
  ],
  "total": 4,
  "limit": 25,
  "offset": 0,
  "has_more": false,
  "next_offset": null
}

2. Get experiment details for pipeline configuration

encode_get_experiment(accession="ENCSR654HRT")

Expected output:

{
  "accession": "ENCSR654HRT",
  "assay_title": "Hi-C",
  "bio_replicate_count": 2,
  "biosample_summary": "heart left ventricle tissue male adult (51 years)",
  "assembly": ["GRCh38"],
  "audit_error_count": 0,
  "audit_warning_count": 1
}

Integration

This skill produces...	Feed into...	Purpose
Chromatin loops (BEDPE)	hic-aggregation	Cross-tissue loop catalog
Loop anchors (BED)	peak-annotation	Assign genes to loop-connected enhancers
Contact matrices (.hic / .mcool)	visualization-workflow	3D genome visualization
Loop-disrupting coordinates	variant-annotation	Identify variants breaking chromatin contacts
QC metrics	quality-assessment	Validate Hi-C library quality
Pipeline parameters	data-provenance	Record BWA/pairtools/Juicer versions
Loop anchor regions	motif-analysis	Discover CTCF motifs at loop anchors

Related Skills

• pipeline-guide -- Parent skill with compute resource assessment and cloud setup
• hic-aggregation -- Aggregate Hi-C loops across samples/tissues
• quality-assessment -- Evaluate pipeline output quality metrics
• data-provenance -- Track all pipeline inputs, outputs, and parameters
• download-encode -- Download ENCODE Hi-C FASTQ files for pipeline input
• publication-trust -- Verify literature claims backing analytical decisions

Presenting Results

When reporting Hi-C pipeline results:

• Valid pair count: Report the UU pair count and its fraction of all parsed pairs from pairs/{sample}.parse_stats.txt. UU is the only pair type this workflow carries forward
• Cis/trans ratio: Report the cis/trans contact ratio (>1.5 pass) and long-range cis fraction (>20kb, >40% expected) from qc/{sample}.contact_stats.txt. These are the primary Hi-C quality indicators
• Contact matrix resolution: Report the achievable resolution based on sequencing depth (e.g., "500M valid pairs supports 5kb resolution") and list the --resolutions actually generated
• Loop counts: Report the number of loops in loops/{sample}.hiccups_loops.bedpe (the file starts with a #chr1 ... header line, so exclude it from the count). HiCCUPS searches the resolutions in --hiccups_resolutions (default 5 kb, 10 kb and 25 kb) and merges them into that single file; per-resolution files stay in the Nextflow work directory
• Matrix paths: Provide paths to the .hic file (Juicebox-compatible) and .mcool file (cooler/HiGlass-compatible)
• Key QC metrics: Present library complexity (unique/total >0.7, from pairs/{sample}.dedup_stats.txt) and the pair-type breakdown from pairs/{sample}.parse_stats.txt in a summary table
• Normalization: Note that .hic carries KR, VC and VC_SQRT vectors (HiCCUPS uses KR) and the .mcool is ICE-balanced by cooler zoomify --balance
• Not produced here: TAD calls, A/B compartments and cooltools outputs are not generated by this workflow; say so rather than implying they are missing
• Next steps: Suggest hic-aggregation for cross-sample loop catalogs, or visualization-workflow for Juicebox/HiGlass session setup

For the request: "$ARGUMENTS"