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.

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

1. 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

1. 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

1. Abdennur & Mirny 2020 - "Cooler: scalable storage for Hi-C data and

other genomically labeled arrays" (Bioinformatics) DOI: 10.1093/bioinformatics/btz540

1. 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

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

Infrastructure parameters (nextflow.config)

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.

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

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):

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

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"