# Stage 1: QC and Trimming

## Input
- Raw FASTQ files (single-end or paired-end)
- Adapter sequences (auto-detected by Trim Galore or specify: Illumina TruSeq)

## Tools
- **FastQC 0.12.1** (image version): Per-base quality, adapter content, duplication rates, GC content
- **Trim Galore 0.6.10** (image version, wraps Cutadapt 4.6): Adapter trimming + quality filtering

## Parameters

| Parameter | Default | Notes |
|-----------|---------|-------|
| Quality cutoff | 20 | Phred score minimum |
| Min length | 36 | Discard reads shorter than this after trimming |
| Adapter | auto-detect | TruSeq for most ENCODE libraries |
| Stringency | 1 | Overlap with adapter sequence required |
| Error rate | 0.1 | Maximum allowed error rate in adapter detection |

The workflow does not expose these as parameters; they are fixed in `main.nf`.

## Commands

The workflow runs the equivalent of:

```bash
# Raw QC
fastqc -t 4 -o qc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz

# Paired-end trimming
trim_galore --paired --quality 20 --length 36 --fastqc \
  --cores 4 -o trimmed/ sample_R1.fastq.gz sample_R2.fastq.gz

# Single-end trimming (--single_end)
trim_galore --quality 20 --length 36 --fastqc \
  --cores 4 -o trimmed/ sample.fastq.gz
```

## Expected Output
- `trimmed/*_trimming_report.txt` -- trimming statistics (reads processed, trimmed, removed)
- `trimmed/*_val_1.fq.gz`, `trimmed/*_val_2.fq.gz` -- trimmed paired-end reads
- `trimmed/*_trimmed.fq.gz` -- trimmed single-end reads
- `fastqc/*.html`, `fastqc/*.zip` -- FastQC reports for raw reads (FASTQC process) and for
  the trimmed reads (`trim_galore --fastqc`). Both land in the same `fastqc/` directory;
  the trimmed-read reports carry `_val_1`/`_val_2`/`_trimmed` in the file name.

All of these feed the MultiQC report.

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Per-base quality | >Q20 across all positions after trimming | Check sequencing run quality |
| Adapter content | <5% after trimming | Verify trimming parameters |
| GC content | Unimodal, matching expected genome GC | Check for contamination |
| Sequence duplication | <50% at this stage | May indicate low complexity library |
| Read count | Record for downstream normalization | No hard threshold at this stage |

## Troubleshooting

- **High adapter content after trimming**: the workflow does not expose `--stringency` or
  `--adapter`; run Trim Galore manually with those flags if auto-detection fails
- **Bimodal GC distribution**: Indicates possible contamination; run FastQ Screen to identify
  organism of origin
- **Very short reads after trimming**: Library insert size may be too short; consider
  investigating library preparation (the 36 bp minimum is fixed in `main.nf`)

# Stage 2: Alignment

## Tools
- **BWA-MEM 0.7.18** (image version): Primary aligner for ChIP-seq (Li & Durbin, 2009)
- **Samtools 1.19** (image version): BAM conversion, sorting, indexing, and statistics

## Reference Genome

The workflow does not download or build an index. It expects a directory (`--bwa_index`,
default `./<genome>_index`) that already contains the FASTA and its BWA index files,
because BWA is invoked as `bwa mem ... <dir>/<genome>.fa`:

```
GRCh38_index/
  GRCh38.fa
  GRCh38.fa.amb
  GRCh38.fa.ann
  GRCh38.fa.bwt
  GRCh38.fa.pac
  GRCh38.fa.sa
```

Build it once from the ENCODE "no alt" analysis-set FASTA:

| Organism | Assembly | FASTA source |
|----------|----------|--------------|
| Human | GRCh38 (hg38) | `https://www.encodeproject.org/files/GRCh38_no_alt_analysis_set_GCA_000001405.15/` |
| Mouse | mm10 (GRCm38) | `https://www.encodeproject.org/files/mm10_no_alt_analysis_set_ENCODE/` |

```bash
mkdir -p GRCh38_index
cp GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta GRCh38_index/GRCh38.fa
bwa index GRCh38_index/GRCh38.fa

# Chromosome sizes for --chrom_sizes (Stage 5 needs them)
samtools faidx GRCh38_index/GRCh38.fa
cut -f1,2 GRCh38_index/GRCh38.fa.fai > GRCh38.chrom.sizes
```

**Important**: Use the ENCODE "no alt" analysis set which excludes alternate haplotype
contigs and decoy sequences. This prevents ambiguous multi-mapping to alternate loci.

## Parameters

| Parameter | Value | Notes |
|-----------|-------|-------|
| BWA algorithm | mem | Recommended for reads >70bp |
| BWA threads (`-t`) | 8 | `BWA_MEM` is configured with 8 CPUs in `nextflow.config` |
| BWA `-M` flag | yes | Mark shorter split hits as secondary (Picard compatible) |
| MAPQ filter | 30 | Remove multi-mappers (MAPQ<30) |
| Sort order | coordinate | Required for Picard and peak calling |
| Sort memory | 2G per thread | `samtools sort -m 2G`; 8 threads x 2G fits the 32 GB process limit |

## Commands

The workflow runs the equivalent of:

```bash
# Paired-end alignment
bwa mem -t 8 -M GRCh38_index/GRCh38.fa trimmed_R1.fq.gz trimmed_R2.fq.gz | \
  samtools view -@ 8 -bS -q 30 - | \
  samtools sort -@ 8 -m 2G -o aligned.bam -
samtools index aligned.bam

# Single-end alignment (--single_end)
bwa mem -t 8 -M GRCh38_index/GRCh38.fa trimmed.fq.gz | \
  samtools view -@ 8 -bS -q 30 - | \
  samtools sort -@ 8 -m 2G -o aligned.bam -
samtools index aligned.bam

# Alignment statistics (the workflow runs flagstat only)
samtools flagstat aligned.bam > flagstat.txt
```

`samtools idxstats` and `samtools stats` are useful manual follow-ups but are not run by
the workflow.

## Expected Output
- `aligned/<sample>.bam` + `.bam.bai` -- coordinate-sorted, indexed BAM
- `aligned/<sample>.flagstat.txt` -- alignment summary (total, mapped, paired, properly paired)

When `--control` is given, the control libraries appear here too, prefixed `CONTROL_`.

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Mapping rate | >80% | Check genome build match, contamination |
| Mapped reads | >=20M (TF), >=45M (histone) | Sequence more or pool |
| MAPQ>=30 fraction | >70% of mapped reads | Check for repetitive regions enrichment |
| Properly paired (PE) | >90% of mapped | Check library preparation |
| Mitochondrial reads | <5% (manual, via `samtools idxstats`) | Not measured or filtered by this workflow |

## Notes

- The MAPQ 30 filter removes reads mapping to multiple locations. For repetitive element
  analysis, consider relaxing this threshold (it is fixed in `main.nf`).
- BWA-MEM is preferred over BWA-ALN for reads longer than 70bp. For older datasets with
  shorter reads, BWA-ALN may be more appropriate.
- The `-M` flag ensures compatibility with Picard MarkDuplicates in Stage 3.
- This workflow performs no mitochondrial filtering at any stage.

# Stage 3: Filtering and Deduplication

## Tools
- **Samtools 1.19** (image version): Flag and MAPQ filtering
- **Picard MarkDuplicates 3.1.1** (image version, Java 17): Remove PCR duplicates
- **bedtools 2.31.0** (image version): Blacklist region filtering

## Order of operations in the workflow

```
aligned.bam -> samtools flag/MAPQ filter -> Picard MarkDuplicates (REMOVE_DUPLICATES=true)
            -> bedtools blacklist filter -> final.bam -> MACS2
```

Duplicates are removed and the blacklist is applied **before** peak calling, not after.

## Blacklist

The default blacklist is downloaded per `--genome`; override it with `--blacklist`:
- **Human (hg38)**: `https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz`
- **Mouse (mm10)**: `https://github.com/Boyle-Lab/Blacklist/raw/master/lists/mm10-blacklist.v2.bed.gz`
- Reference: Amemiya et al. 2019 (Scientific Reports, ~1,372 citations)

The hg38 blacklist contains ~900 regions covering ~40 Mb of problematic sequence
including high-signal artifacts, satellite repeats, and assembly gaps. Gzipped BED is
accepted directly.

## Samtools Flag Filtering

For paired-end data the workflow uses `-F 1804`, which removes:
- Bit 4: read unmapped
- Bit 8: mate unmapped
- Bit 256: secondary alignment
- Bit 512: read fails quality checks
- Bit 1024: PCR duplicate (any already flagged by the aligner)

For `--single_end` it uses `-F 1028` (unmapped + duplicate). Neither flag set removes
mitochondrial reads; this workflow does not filter chrM.

## Commands

The workflow runs the equivalent of:

```bash
# Flag and MAPQ filter (paired-end)
samtools view -@ 4 -b -F 1804 -q 30 aligned.bam | \
  samtools sort -@ 4 -o filtered.bam -

# Mark and remove PCR duplicates
picard MarkDuplicates \
  INPUT=filtered.bam \
  OUTPUT=dedup.bam \
  METRICS_FILE=dup_metrics.txt \
  REMOVE_DUPLICATES=true \
  VALIDATION_STRINGENCY=LENIENT

# Remove reads overlapping blacklist regions
bedtools intersect -v -abam dedup.bam -b hg38-blacklist.v2.bed.gz > final.bam
samtools index final.bam

# Record final read count
samtools flagstat final.bam > final_flagstat.txt
```

## Expected Output
- `filtered/<sample>.dup_metrics.txt` -- Picard duplication metrics
- `filtered/<sample>.final.bam` + `.final.bam.bai` -- blacklist-filtered, ready for peak calling
- `filtered/<sample>.final.flagstat.txt` -- final read count after all filtering

The intermediate `filtered.bam` and `dedup.bam` stay in the Nextflow work directory.

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Duplication rate | <30% (Picard `PERCENT_DUPLICATION`) | Low-input library; consider re-prep |
| Post-filter read count | >=20M TF / >=45M histone | May need deeper sequencing |
| Blacklist overlap | <1% of reads | Expected; higher suggests artifacts |

## Library Complexity Metrics (manual)

NRF, PBC1 and PBC2 are **not computed by this workflow**, and Picard MarkDuplicates does
not report them: its metrics file gives duplicate counts and `PERCENT_DUPLICATION`, while
PBC needs per-position read counts. Compute them from the pre-deduplication BAM if needed:

- **NRF** = distinct read start positions / total reads (Non-Redundant Fraction)
- **PBC1** = positions with exactly 1 read / distinct positions
- **PBC2** = positions with exactly 1 read / positions with exactly 2 reads

```bash
# From the published pre-deduplication BAM (already MAPQ >= 30 from the alignment step)
samtools view -F 1804 -q 30 results/aligned/sample.bam | \
  awk 'BEGIN{OFS="\t"}{print $3, $4, ($2%32>=16 ? "-" : "+")}' | \
  sort | uniq -c | \
  awk '{total += $1; distinct += 1; if ($1 == 1) one += 1; if ($1 == 2) two += 1}
       END {printf "NRF=%.3f PBC1=%.3f PBC2=%.3f\n", distinct/total, one/distinct, one/two}'
```

This is approximate; the ENCODE reference implementation derives the 5' positions with
`bedtools bamtobed` instead of the raw POS field.

These metrics quantify PCR amplification bottleneck severity.

# Stage 4: Peak Calling and IDR

## Tools
- **MACS2 2.2.9.1** (image version): Model-based Analysis of ChIP-Seq (Zhang et al. 2008, ~7,000 citations)
- **IDR 2.0.4.2** (image version): Irreproducible Discovery Rate (Li et al. 2011, ~1,500 citations)
- **bedtools 2.31.0 + samtools 1.19** (image versions): FRiP counting after peak calling

## MACS2 Parameters

The workflow builds one MACS2 command per ChIP sample. `--peak_type` selects the narrow or
broad column for the whole run.

| Parameter | Narrow (TF) | Broad (Histone) | Notes |
|-----------|-------------|-----------------|-------|
| `--format` | BAMPE | BAMPE | `BAM` when `--single_end` is set |
| `--gsize` | hs (2.7e9) | hs (2.7e9) | `mm` when `--genome mm10` |
| `--qvalue` | 0.05 | 0.05 | FDR threshold |
| `--broad` | no | yes | Broad peak mode for repressive marks |
| `--broad-cutoff` | n/a | 0.1 | Linking threshold for broad peaks |
| `--nomodel` | yes | yes | Fragment sizes come from the BAM |
| `--keep-dup` | all | all | Duplicates were already removed in Stage 3 |
| `--call-summits` | yes | no | Subpeak summit positions |
| `-B` | yes | yes | bedGraphs for the Stage 5 signal tracks |
| `-c` | only with `--control` | only with `--control` | All control BAMs are passed together; MACS2 pools them |

## Narrow vs Broad Mark Decision

| Peak Type | Targets | Rationale |
|-----------|---------|-----------|
| **Narrow** | H3K4me3, H3K4me1, H3K27ac, H3K9ac, all TFs, CTCF | Punctate binding pattern |
| **Broad** | H3K27me3, H3K36me3, H3K9me3, H3K79me2 | Diffuse domain spreading |

## Commands

The workflow runs the equivalent of:

```bash
# Narrow peaks (TF and active histone marks)
macs2 callpeak -t treatment.bam -c control.bam \
  -f BAMPE -g hs -n sample \
  --qvalue 0.05 --nomodel --keep-dup all --call-summits -B

# Broad peaks (repressive histone marks)
macs2 callpeak -t treatment.bam -c control.bam \
  -f BAMPE -g hs -n sample \
  --qvalue 0.05 --nomodel --keep-dup all --broad --broad-cutoff 0.1 -B

# IDR, once per pair of replicates, the two names sorted alphabetically (narrow runs only)
idr --samples rep1_peaks.narrowPeak rep2_peaks.narrowPeak \
  --input-file-type narrowPeak \
  --rank p.value \
  --output-file rep1_vs_rep2.idr_peaks.txt \
  --plot \
  --idr-threshold 0.05
```

## FRiP (workflow)

After peak calling the workflow computes the fraction of reads in peaks for every
treatment sample, in both narrow and broad mode, and publishes
`qc/<sample>.frip_mqc.tsv`. The command and the output columns are in
`references/05-qc-metrics.md`.

## What the workflow publishes

- `peaks/<narrow|broad>/<sample>_peaks.narrowPeak` (or `_peaks.broadPeak` and
  `_peaks.gappedPeak` in broad mode)
- `peaks/<narrow|broad>/<sample>_peaks.xls`
- `peaks/<narrow|broad>/<sample>_treat_pileup.bdg` and `<sample>_control_lambda.bdg`
  (consumed by Stage 5)
- `peaks/idr/<sampleA>_vs_<sampleB>.idr_peaks.txt` and, when IDR emits it,
  `<sampleA>_vs_<sampleB>.idr_peaks.txt.png`, one pair of files per replicate pair
- `qc/<sample>.frip_mqc.tsv` (see `references/05-qc-metrics.md`)

`<sample>_summits.bed` is produced by `--call-summits` and published to `peaks/narrow/`;
broad runs do not produce it.

## What the IDR step does and does not do

All samples matched by `--reads` are treated as replicates of one experiment. The workflow
sorts them by sample name and runs `idr` once for every pair. Consequently:

- Only narrow runs get IDR (`--peak_type broad` produces no `peaks/idr/`).
- Two samples give one comparison, three give three, four give six. No replicate is
  dropped.
- With a single peak file there is no pair, so IDR is skipped silently.
- There are **no** pooled-replicate calls, **no** pseudoreplicates, and **no** optimal /
  conservative peak sets.
- Rescue ratio and self-consistency ratio are **not computed**. Computing them requires
  pooled and pseudoreplicated peak calls that this workflow does not produce.

## Manual follow-ups (not run by this workflow)

```bash
# Self-pseudoreplicate IDR: split one replicate's reads in half, call peaks on each half,
# then compare. Requires re-running MACS2 outside the workflow.
idr --samples pr1_peaks.narrowPeak pr2_peaks.narrowPeak \
  --input-file-type narrowPeak \
  --output-file self_idr_peaks.txt \
  --plot
```

Pooled peak calls across replicates are likewise outside this workflow.

## IDR Interpretation

| Metric | Expected (Good) | Concern | Computed here? |
|--------|-----------------|---------|----------------|
| IDR peaks (0.05 threshold) | 50,000-200,000 (TF) | <20,000 suggests poor enrichment | yes |
| IDR / individual rep ratio | 0.3-0.7 | <0.3 very stringent; >0.7 very lenient | manual (count lines) |
| Rescue ratio | <2 | >2 suggests replicate discordance | no |
| Self-consistency ratio | <2 | >2 suggests noisy data | no |

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| FRiP (`qc/<sample>.frip_mqc.tsv`) | >=0.01 (1%) | Poor enrichment; check antibody |
| IDR peaks | >20,000 (TF) | Low enrichment or poor replicates |
| Peak count ratio | IDR = 30-70% of rep peaks | Replicate consistency issue |

## Notes

- Peaks are called on individual replicates; the IDR step then compares every pair of them.
- For broad marks, IDR is not standard practice, which matches the workflow: broad runs
  skip it. Use pooled replicate peak calls instead, produced manually.
- The `--call-summits` flag identifies subpeak summits within broader peak regions,
  useful for motif analysis downstream.
- FRiP is calculated by the workflow as: alignments of `filtered/<sample>.final.bam` that
  overlap a peak / all alignments of that BAM. Control libraries are not peak-called, so
  they get no FRiP row.

# Stage 5: Signal Tracks and QC Report

Stage 5 of the workflow does three things: it converts the MACS2 bedGraphs into bigWig
signal tracks, it computes FRiP for every treatment sample, and it runs MultiQC over the
logs collected in earlier stages. The sections marked "manual step" below are **not run by
the workflow**.

## Signal Track Generation (workflow)

Signal tracks provide normalized coverage for genome browser visualization. The workflow
runs the equivalent of the following per sample, using the `_treat_pileup.bdg` and
`_control_lambda.bdg` files MACS2 wrote with `-B`. MACS2 produces `_control_lambda.bdg`
even without `-c`, so the tracks exist with or without `--control`; without a control the
"fold change" track is fold enrichment over the local lambda model.

```bash
# Fold enrichment over the MACS2 control/lambda track
macs2 bdgcmp -t sample_treat_pileup.bdg -c sample_control_lambda.bdg \
  -o fc.bdg -m FE
sort -k1,1 -k2,2n fc.bdg > fc.sorted.bdg
bedGraphToBigWig fc.sorted.bdg GRCh38.chrom.sizes sample.fc.bw

# Signal p-value track (statistical significance)
macs2 bdgcmp -t sample_treat_pileup.bdg -c sample_control_lambda.bdg \
  -o pval.bdg -m ppois
sort -k1,1 -k2,2n pval.bdg > pval.sorted.bdg
bedGraphToBigWig pval.sorted.bdg GRCh38.chrom.sizes sample.pval.bw
```

Published as `signal/<sample>.fc.bw` and `signal/<sample>.pval.bw`. `--chrom_sizes` is
required for this stage, which is why the workflow stops at startup when it is missing.

## FRiP (workflow)

The workflow computes the fraction of reads in peaks for every treatment sample, in both
narrow and broad mode, and publishes `qc/<sample>.frip_mqc.tsv`. It runs the equivalent of:

```bash
TOTAL_READS=$(samtools view -c results/filtered/sample.final.bam)
READS_IN_PEAKS=$(bedtools intersect -u -a results/filtered/sample.final.bam \
  -b results/peaks/narrow/sample_peaks.narrowPeak | samtools view -c -)
awk -v a=$READS_IN_PEAKS -v b=$TOTAL_READS 'BEGIN{printf "%.4f\n", a/b}'
```

The division uses `awk`. The published file carries a MultiQC header and one row per peak
file, so MultiQC renders it as the section "Fraction of reads in peaks":

| Column | Meaning |
|--------|---------|
| `Peak set` | Name of the peak file the row was computed from |
| `FRiP` | `reads_in_peaks / total_reads`, printed to four decimals |
| `reads_in_peaks` | Alignments of the final BAM overlapping a peak (`bedtools intersect -u`) |
| `total_reads` | All alignments of the final BAM (`samtools view -c`) |

FRiP here is a fraction, not a percentage: the ENCODE >=1% standard is >=0.01 in this
file. Control libraries are not peak-called, so no `CONTROL_<name>.frip_mqc.tsv` is
written.

## MultiQC Aggregated Report (workflow)

```bash
multiqc . -o . -f
```

Published as `qc/multiqc/multiqc_report.html` with `qc/multiqc/multiqc_data/`. The workflow
feeds it exactly these inputs:

- FastQC reports for raw reads
- FastQC reports for trimmed reads (from `trim_galore --fastqc`)
- Trim Galore trimming reports
- `samtools flagstat` from the alignment step
- Picard MarkDuplicates metrics
- `samtools flagstat` from the blacklist-filtered BAM
- `qc/<sample>.frip_mqc.tsv` from the FRiP step, shown as "Fraction of reads in peaks"

MACS2 output is not passed to MultiQC.

## Comprehensive QC Metrics Table

| Metric | Tool | Threshold | Computed by the workflow? |
|--------|------|-----------|---------------------------|
| Total reads | FastQC, samtools flagstat | Record | yes |
| Mapped reads | samtools flagstat | >=20M TF / >=45M histone | yes |
| Mapping rate | samtools flagstat | >80% | yes |
| Duplication rate | Picard | <30% | yes |
| IDR peaks | IDR | >20,000 (TF) | yes (narrow runs, one file per replicate pair) |
| NRF | manual (see 03-filtering.md) | >=0.8 | no |
| PBC1 | manual (see 03-filtering.md) | >=0.8 | no |
| PBC2 | manual (see 03-filtering.md) | >=3 | no |
| NSC | phantompeakqualtools | >1.05 | no |
| RSC | phantompeakqualtools | >0.8 | no |
| FRiP | bedtools + samtools (`qc/<sample>.frip_mqc.tsv`) | >=0.01 (1%) | yes |
| Mitochondrial fraction | `samtools idxstats` | <5% | no |

## Manual step: Strand Cross-Correlation (phantompeakqualtools)

Computes NSC and RSC, which measure ChIP enrichment quality independent of peak calling.
Reference: Kharchenko et al. 2008; Landt et al. 2012.

**Not run by this workflow, and phantompeakqualtools is not in the pipeline image.** Run it
in the conda environment `bioinformatics-installer/environments/chipseq-env.yml`, which
includes `phantompeakqualtools=1.2.2`.

```bash
conda env create -f chipseq-env.yml && conda activate encode-chipseq
run_spp.R -c=results/filtered/sample.final.bam -savp=cc_plot.pdf -out=cc_scores.txt -tmpdir=tmp/
# Output columns: filename, numReads, estFragLen, corr_estFragLen,
#                 phantomPeak, corr_phantomPeak, argmin_corr, min_corr,
#                 NSC, RSC, QualityTag
```

## Manual step: Fingerprint Plot (deeptools)

**Not run by this workflow**, though deeptools 3.5.5 is installed in the image:

```bash
plotFingerprint -b results/filtered/chip.final.bam results/filtered/CONTROL_input.final.bam \
  --labels ChIP Input \
  --plotFile fingerprint.pdf \
  --outRawCounts fingerprint_counts.txt
```

The fingerprint plot visually shows enrichment: a good ChIP sample curves away
from the diagonal (uniform coverage), while input stays close to the diagonal.

# ChIP-seq Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-chipseq skill — papers defining ENCODE ChIP-seq processing standards and the tools used at each pipeline stage.

---

## ENCODE Pipeline Standards

These papers define the official ENCODE consortium standards that this pipeline implements.

---

### Landt et al. 2012 — ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia

- **Citation:** Landt SG, Marinov GK, Kundaje A, Kheradpour P, Pauli F, Batzoglou S, Bernstein BE, Bickel P, Brown JB, Cayting P, Chen Y, DeSalvo G, Epstein C, Fisher-Aylor KI, Euskirchen G, Gerstein M, Gertz J, Hartemink AJ, Hoffman MM, Iyer VR, Jung YL, Karmakar S, Kellis M, Kharchenko PV, Li Q, Liu T, Liu XS, Ma L, Milosavljevic A, Myers RM, Park PJ, Pazin MJ, Perry MD, Raha D, Reddy TE, Rozowsky J, Shoresh N, Sidow A, Slattery M, Stamatoyannopoulos JA, Tolstorukov MY, White KP, Xi S, Farnham PJ, Lieb JD, Wold BJ, Snyder M. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Research, 22(9):1813-1831, 2012.
- **DOI:** [10.1101/gr.136184.111](https://doi.org/10.1101/gr.136184.111)
- **PMID:** 22955991 | **PMC:** PMC3431496
- **Citations:** ~4,000
- **Key findings:** Foundational ENCODE standards paper establishing quality metrics and experimental guidelines for ChIP-seq. Defined minimum sequencing depth requirements (20M uniquely mapped reads for TF, 45M for histone modifications), antibody validation criteria (immunoprecipitation followed by western blot, mass spectrometry, or immunofluorescence), and replicate concordance requirements (at least 2 biological replicates). Introduced the FRiP (Fraction of Reads in Peaks) metric with a minimum threshold of 1% and the cross-correlation analysis framework (NSC > 1.05, RSC > 0.8) for assessing enrichment quality. Established the traffic-light QC system where no single metric is sufficient — quality must be interpreted collectively. These standards remain the basis for all ENCODE ChIP-seq data production and are implemented directly in this pipeline's QC stage.

---

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

- **Citation:** ENCODE Project Consortium, Moore JE, Purcaro MJ, Pratt HE, Epstein CB, Shoresh N, Adrian J, Kawli T, Davis CA, Dobin A, Kaul R, Halow J, Van Nostrand EL, Freese P, Gorkin DU, Shen Y, He Y, Mackiewicz M, Pauli-Behn F, Williams BA, Mortazavi A, Keller CA, Zhang XO, Amin SI, Hardy M, Mber F, Sandstrom R, Bernstein BE, Wold BJ, Kundaje A, Stam M, Partridge EC, Bristow CA, Gerstein M, Gingeras TR, Stamatoyannopoulos JA, Weng Z, Snyder M, Birney E, Myers RM, Hardison RC, Ren B, Cherry JM, Bernstein BE. 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:** Phase 3 ENCODE paper describing the expanded encyclopedia of regulatory elements. Introduced the Registry of candidate cis-Regulatory Elements (cCREs), which classified 926,535 human and 339,815 mouse cCREs into promoter-like, enhancer-like, and CTCF-only categories using integrated ChIP-seq and chromatin accessibility data. Established the uniform processing pipeline framework where all ENCODE data are processed through standardized pipelines to ensure cross-experiment comparability. Defined updated quality standards including the use of IDR (Irreproducible Discovery Rate) for replicate consistency assessment, preferred_default file selection criteria, and assembly standards (GRCh38 for human, mm10 for mouse). This paper provides the scientific basis for the pipeline's parameter choices and output format standards.

---

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

- **Citation:** Hitz BC, Lee JW, Jolanki O, Kagda MS, Graham K, Sud P, Gabdank I, Strattan JS, Sloan CA, Dreszer T, Rowe LD, Podduturi NR, Malladi V, Chan ET, Davidson JM, Ho M, Miyasato S, Simison M, Tanaka F, Luo Y, Whaling I, Hong EL, Lee BT, Sandstrom R, Rynes E, Nelson J, Nishida A, Ingersoll A, Buckley M, Frerker M, Kim DS, Boley N, Trout D, Dobin A, Rahmanian S, Wyman D, Balderrama-Gutierrez G, Reese F, Durand NC, Dudchenko O, Weisz D, Rao SSP, Blackburn A, Gkountaroulis D, Sadr M, Olshansky M, Eliaz Y, Nguyen D, Bochkov ID, Shamim M, Mahajan R, Lieberman Aiden E, Gingeras TR, Heath S, Hirst M, Kent WJ, Kundaje A, Mortazavi A, Wold BJ, Cherry JM. 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:** Comprehensive description of all ENCODE uniform analysis pipelines including ChIP-seq, ATAC-seq, RNA-seq, WGBS, Hi-C, and others. Documents the infrastructure for pipeline distribution using Docker containers and Workflow Description Language (WDL), enabling execution on local machines, HPC clusters, or cloud environments via Cromwell. Details the quality control metrics and thresholds applied uniformly across all ENCODE data, including library complexity metrics (NRF, PBC1, PBC2), enrichment metrics (NSC, RSC, FRiP), and reproducibility metrics (IDR). Establishes the principle that standardized computational methodologies are prerequisites for successful integrative analyses across different ENCODE collections. This paper serves as the authoritative reference for the pipeline architecture and parameter choices implemented in this skill.

---

## Core Pipeline Tools

Papers for the primary algorithms used at each pipeline stage. These tools are selected to match ENCODE uniform pipeline specifications.

---

### Zhang et al. 2008 — Model-based Analysis of ChIP-Seq (MACS)

- **Citation:** Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. Model-based Analysis of ChIP-Seq (MACS). Genome Biology, 9(9):R137, 2008.
- **DOI:** [10.1186/gb-2008-9-9-r137](https://doi.org/10.1186/gb-2008-9-9-r137)
- **PMID:** 18798982 | **PMC:** PMC2592715
- **Citations:** ~7,000
- **Key findings:** Introduced the MACS peak-calling algorithm, which models the shift size of ChIP-seq tags empirically from the data and uses a dynamic local Poisson model to capture local biases in the genome. MACS accounts for both the DNA fragment size and the local background to improve peak detection sensitivity and specificity. The algorithm operates in two modes: narrow peak calling (default, suitable for TF and sharp histone marks like H3K4me3) with q-value cutoff of 0.05, and broad peak calling (for diffuse histone marks like H3K27me3 and H3K36me3) using --broad with a relaxed cutoff of 0.1. MACS2 remains the ENCODE-standard peak caller for ChIP-seq and is the primary peak caller in this pipeline.

---

### Li et al. 2011 — Measuring reproducibility of high-throughput experiments (IDR)

- **Citation:** Li Q, Brown JB, Huang H, Bickel PJ. Measuring reproducibility of high-throughput experiments. Annals of Applied Statistics, 5(3):1752-1779, 2011.
- **DOI:** [10.1214/11-AOAS466](https://doi.org/10.1214/11-AOAS466)
- **Citations:** ~1,500
- **Key findings:** Developed the Irreproducible Discovery Rate (IDR) framework for assessing the reproducibility of high-throughput experiments, specifically designed for peak-based genomic assays like ChIP-seq. IDR uses a copula mixture model to distinguish signal from noise by comparing ranked peak lists from biological replicates, classifying each peak as either belonging to the reproducible signal component or the irreproducible noise component. The IDR threshold of 0.05 is used by ENCODE to define the final peak set from replicated experiments. IDR is preferred over simple overlap-based approaches because it uses the quantitative ranking information (e.g., peak signal or p-value) rather than just peak presence/absence. This pipeline applies IDR analysis as the final stage before generating the consensus peak set.

---

### Li 2013 — BWA-MEM alignment algorithm

- **Citation:** Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997, 2013.
- **DOI:** [10.48550/arXiv.1303.3997](https://doi.org/10.48550/arXiv.1303.3997)
- **Citations:** ~11,600
- **Key findings:** Introduced BWA-MEM, an alignment algorithm for mapping sequence reads against large reference genomes. BWA-MEM automatically chooses between local and end-to-end alignments, supports paired-end reads, and performs chimeric alignment. The algorithm is robust to sequencing errors and applicable to a wide range of sequence lengths from 70bp to several megabases. For mapping 100bp sequences typical of ChIP-seq, BWA-MEM outperforms other contemporary aligners in both speed and accuracy. BWA-MEM is the ENCODE-standard aligner for ChIP-seq data and is used in this pipeline's alignment stage with default parameters plus MAPQ filtering at Q30.

---

### Li et al. 2009 — The Sequence Alignment/Map format and SAMtools

- **Citation:** Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16):2078-2079, 2009.
- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352)
- **PMID:** 19505943 | **PMC:** PMC2723002
- **Citations:** ~53,700
- **Key findings:** Defined the SAM/BAM format specification, now the universal standard for storing read alignments against reference sequences. SAMtools provides utilities for post-processing alignments: indexing, sorting, merging, filtering by mapping quality, and extracting statistics. In the ChIP-seq pipeline, samtools is used for coordinate sorting, MAPQ-based filtering (removing reads with MAPQ < 30), removal of unmapped reads, mitochondrial read filtering, and BAM indexing. The SAM/BAM format enables interoperability between all downstream tools in the pipeline.

---

## Quality Control & Visualization

---

### Ramírez et al. 2016 — deepTools2

- **Citation:** Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dündar F, Manke T. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Research, 44(W1):W160-W165, 2016.
- **DOI:** [10.1093/nar/gkw257](https://doi.org/10.1093/nar/gkw257)
- **PMID:** 27079975 | **PMC:** PMC4987876
- **Citations:** ~6,100
- **Key findings:** Comprehensive suite for processing and visualizing deep-sequencing data. Provides tools for generating normalized coverage tracks (bamCoverage for RPKM/CPM bigWig files), computing enrichment heatmaps (computeMatrix + plotHeatmap), fingerprint plots for assessing IP enrichment (plotFingerprint), and correlation analyses between replicates (multiBamSummary + plotCorrelation). In this pipeline, deepTools generates the signal tracks (bigWig files) using RPGC normalization and produces QC visualizations including fingerprint plots that distinguish well-enriched from failed ChIP experiments.

---

### Ewels et al. 2016 — MultiQC

- **Citation:** Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19):3047-3048, 2016.
- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354)
- **PMID:** 27312411 | **PMC:** PMC5039924
- **Citations:** ~6,800
- **Key findings:** Tool for aggregating QC results from multiple bioinformatics tools into a single interactive HTML report. Parses output from FastQC, Picard, Samtools, MACS2, deepTools, and many other tools. Enables rapid identification of batch effects, outlier samples, and systematic quality issues across an entire project. In this pipeline, MultiQC is the final QC step that aggregates metrics from all preceding stages into a comprehensive report, providing the overview needed for the traffic-light QC assessment.

---

### Ramachandran et al. 2013 — phantompeakqualtools

- **Citation:** Kharchenko PV, Tolstorukov MY, Park PJ. Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nature Biotechnology, 26(12):1351-1359, 2008. (Algorithm basis); Landt et al. 2012 (ENCODE implementation); phantompeakqualtools software by Anshul Kundaje.
- **DOI:** [10.1038/nbt.1508](https://doi.org/10.1038/nbt.1508) (Kharchenko et al.)
- **PMID:** 19029915
- **Citations:** ~2,500 (Kharchenko et al.)
- **Key findings:** phantompeakqualtools computes strand cross-correlation metrics that quantify ChIP enrichment independent of peak calling. The Normalized Strand Coefficient (NSC) measures the ratio of the cross-correlation peak to the background minimum — values > 1.05 indicate acceptable enrichment. The Relative Strand Correlation (RSC) measures the ratio of the fragment-length cross-correlation to the read-length peak — values > 0.8 indicate the ChIP signal exceeds the phantom peak artifact. These metrics were adopted by ENCODE as mandatory QC measures (Landt et al. 2012) and are computed in this pipeline's QC stage for every ChIP-seq sample.

---

## Supplementary Tools (Non-ENCODE-Specific)

These tools are widely used in genomics pipelines and are incorporated here but are not specific to the ENCODE ChIP-seq standard.

---

### Andrews 2010 — FastQC

- **Citation:** Andrews S. FastQC: A quality control tool for high throughput sequence data. Babraham Bioinformatics, 2010.
- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **Citations:** Widely cited (no formal publication; referenced via URL)
- **Key findings:** De facto standard tool for raw sequencing data quality assessment. Generates per-base and per-sequence quality scores, GC content distribution, sequence length distribution, adapter content detection, sequence duplication levels, and overrepresented sequence identification. FastQC's modular report format integrates directly with MultiQC for batch-level quality assessment. Used in this pipeline's first stage to assess raw FASTQ quality before trimming.

---

### Martin 2011 — Cutadapt (basis for Trim Galore)

- **Citation:** Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 17(1):10-12, 2011.
- **DOI:** [10.14806/ej.17.1.200](https://doi.org/10.14806/ej.17.1.200)
- **Citations:** ~13,000
- **Key findings:** Introduced Cutadapt, an adapter trimming tool that finds and removes adapter sequences, primers, poly-A tails, and other unwanted sequences from high-throughput sequencing reads. Supports color-space data, allows multiple adapter types per read, and handles partial adapter matches. Trim Galore, used in this pipeline, is a wrapper around Cutadapt that adds automatic adapter detection and integrates quality trimming with a default Phred score cutoff of 20. Trim Galore also runs FastQC on trimmed output for post-trimming quality verification.

---

### Broad Institute — Picard MarkDuplicates

- **Citation:** Broad Institute. Picard toolkit. GitHub, 2019.
- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **Citations:** Widely cited (no formal publication; referenced via URL)
- **Key findings:** Picard MarkDuplicates identifies and flags PCR and optical duplicate reads in aligned BAM files. Duplicates arise from PCR amplification during library preparation where multiple copies of the same template molecule are sequenced. In ChIP-seq, duplicate removal is essential because PCR duplicates inflate peak signal intensity and can create false positive peaks. Picard uses read alignment coordinates (and optionally molecular barcodes) to identify duplicates, keeping the read with the highest base quality sum. Library complexity metrics (NRF, PBC1, PBC2) are derived from Picard's duplication statistics. ENCODE recommends NRF >= 0.8, PBC1 >= 0.8, PBC2 >= 3 for acceptable library complexity.

---

### Quinlan & Hall 2010 — BEDTools

- **Citation:** Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 26(6):841-842, 2010.
- **DOI:** [10.1093/bioinformatics/btq033](https://doi.org/10.1093/bioinformatics/btq033)
- **PMID:** 20110278 | **PMC:** PMC2832824
- **Citations:** ~12,000
- **Key findings:** Suite of utilities for genomic interval arithmetic: intersection, union, subtraction, merging, and complement operations on BED, BAM, VCF, and GFF files. In the ChIP-seq pipeline, bedtools is used for blacklist region removal (bedtools intersect -v), peak annotation (bedtools intersect with gene models), and computing the Fraction of Reads in Peaks (FRiP) metric. The tool's efficiency with large genomic datasets and streaming I/O support make it suitable for integration into automated pipeline workflows.

---

### Amemiya et al. 2019 — The ENCODE Blacklist

- **Citation:** Amemiya HM, Kundaje A, Boyle AP. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Scientific Reports, 9:9354, 2019.
- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z)
- **PMID:** 31249361 | **PMC:** PMC6597582
- **Citations:** ~1,372
- **Key findings:** Defined comprehensive blacklists for human (hg19/GRCh38), mouse (mm9/mm10), worm, and fly genomes containing regions that produce anomalous, unstructured, or high signal in next-generation sequencing experiments regardless of cell type or experiment. These regions include collapsed repeats, satellite sequences, and other assembly artifacts. Removing blacklisted regions is described as an essential quality measure for all functional genomics analyses. The ENCODE blacklist v2 for GRCh38 contains 910 regions (totaling ~36 Mb). In this pipeline, blacklist filtering is applied immediately after duplicate marking and before peak calling, ensuring that artifact peaks are excluded from all downstream analyses.

---

## Additional Relevant Papers

---

### Bailey et al. 2013 — Practical Guidelines for Comprehensive ChIP-seq Analysis

- **Citation:** Bailey T, Krajewski P, Ladunga I, Lefebvre C, Li Q, Liu T, Madrigal P, Taslim C, Zhang J. Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data. PLoS Computational Biology, 9(11):e1003326, 2013.
- **DOI:** [10.1371/journal.pcbi.1003326](https://doi.org/10.1371/journal.pcbi.1003326)
- **PMID:** 24244136 | **PMC:** PMC3828144
- **Citations:** ~274
- **Key findings:** Step-by-step computational guidelines complementing the ENCODE standards, covering sequencing depth selection, quality checking, mapping, data normalization, reproducibility assessment, peak calling, differential binding analysis, FDR control, peak annotation, visualization, and motif analysis. Provides a decision tree for selecting appropriate peak callers based on the type of ChIP-seq experiment (TF vs. histone) and discusses the impact of control sample selection on peak calling accuracy. Particularly useful for understanding the rationale behind the parameter choices in this pipeline.

---

### Danecek et al. 2021 — Twelve years of SAMtools and BCFtools

- **Citation:** Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, Li H. Twelve years of SAMtools and BCFtools. GigaScience, 10(2):giab008, 2021.
- **DOI:** [10.1093/gigascience/giab008](https://doi.org/10.1093/gigascience/giab008)
- **PMID:** 33590861 | **PMC:** PMC7931819
- **Citations:** ~9,400
- **Key findings:** Retrospective on 12 years of SAMtools development documenting major improvements including multi-threaded I/O, CRAM format support, improved variant calling, and the modular HTSlib library. Both SAMtools and BCFtools have been installed over 1 million times via Bioconda. Documents the samtools flagstat, samtools idxstats, and samtools stats utilities used in this pipeline for alignment QC metrics including mapping rate, properly paired percentage, and mitochondrial fraction calculation.