git:20260921.ed45b8d to git:20260921.4ba10c6

1 added, 1 removed. Audit A to A.

---

name: pipeline-hic

description: "Execute ENCODE Hi-C pipeline from FASTQ to contact matrices and loop calls. Child of pipeline-guide. Provides Nextflow execution with Docker and cloud deployment. Use when processing Hi-C data, generating contact matrices, or calling loops. Trigger on: Hi-C pipeline, chromatin conformation, contact matrix, loop calling, TAD detection, Juicer, HiCCUPS, 3D genome."

---

# 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

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)

```bash

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

```bash

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)

```bash

# 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

|------|------|-----|---------------------|

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

+ | `--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.

|--------|------|---------|------|---------------|

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

```python

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:

```json

{

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

```json

[

{"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:

```bash

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

```bash

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:

```json

{

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

```json

{

"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"