# Stage 1: QC and Trimming

## Input
- Raw paired-end FASTQ files (this workflow is paired-end only)
- Adapter sequences: Nextera transposase adapters (not TruSeq)

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

## Key Difference from ChIP-seq
ATAC-seq uses **Nextera** transposase adapters, not Illumina TruSeq. The workflow passes
`--nextera` explicitly rather than relying on auto-detection.

## Parameters

| Parameter | Default | Notes |
|-----------|---------|-------|
| Quality cutoff | 20 | Phred score minimum |
| Min length | 20 | Shorter than ChIP-seq due to NFR fragments |
| Adapter | Nextera (`--nextera`) | Tn5 transposase adapters |
| Stringency | 1 | Overlap with adapter sequence required |

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 with Nextera adapters
trim_galore --paired --nextera --quality 20 --length 20 --fastqc \
  --cores 4 -o trimmed/ sample_R1.fastq.gz sample_R2.fastq.gz
```

## Expected Output
- `trimmed/*_trimming_report.txt` -- trimming statistics
- `trimmed/*_val_1.fq.gz`, `trimmed/*_val_2.fq.gz` -- trimmed paired-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` in the file name.

All of these feed the MultiQC report.

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Per-base quality | >Q20 after trimming | Check sequencing run quality |
| Adapter content | <5% after trimming | Verify Nextera adapter detection |
| GC content | Unimodal, matching genome | Check for contamination |
| Insert size | Nucleosomal ladder visible | Inspect library prep |
| Read count | >=50M total recommended | May need more sequencing |

## Troubleshooting
- **High adapter content**: ATAC-seq libraries with short inserts have more adapter
  read-through. This is normal for NFR fragments; the fixed `--length 20` keeps them.
- **No nucleosomal pattern in insert sizes**: May indicate failed transposition or
  over-transposition. Check Tn5:cell ratio.

# Stage 2: Alignment

## Tools
- **Bowtie2 2.5.4** (image version): Primary aligner for ATAC-seq (Langmead & Salzberg 2012, ~30,000 citations)
- **Samtools 1.19** (image version): BAM conversion, sorting, indexing, and statistics

## Why Bowtie2 Instead of BWA
Bowtie2 is preferred for ATAC-seq because:
1. Better handling of short fragments (NFR <150 bp)
2. `--very-sensitive` mode provides optimal alignment for ATAC-seq read characteristics
3. Concordant paired-end alignment with fragment size constraints
4. Standard in the ENCODE ATAC-seq pipeline

## Reference Genome

The workflow does not download or build an index. It expects a directory
(`--bowtie2_index`, default `./<genome>_bowtie2_index`) whose index files are named after
the genome, because Bowtie2 is invoked as `bowtie2 ... -x <dir>/<genome>`:

```
GRCh38_bowtie2_index/
  GRCh38.1.bt2  GRCh38.2.bt2  GRCh38.3.bt2  GRCh38.4.bt2
  GRCh38.rev.1.bt2  GRCh38.rev.2.bt2
```

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_bowtie2_index
bowtie2-build --threads 8 GRCh38.fa GRCh38_bowtie2_index/GRCh38
```

Use the ENCODE "no alt" analysis set (excludes alternate haplotypes and decoys).

## Parameters

| Parameter | Value | Notes |
|-----------|-------|-------|
| Bowtie2 mode | --very-sensitive | Optimal sensitivity for ATAC-seq |
| Threads | 8 | `BOWTIE2_ALIGN` is configured with 8 CPUs in `nextflow.config` |
| Max fragment size (-X) | 2000 | Accommodate di/tri-nucleosomal fragments |
| No mixed (--no-mixed) | yes | Require both mates to align |
| No discordant (--no-discordant) | yes | Require concordant alignment |
| MAPQ filter | 30 | Remove multi-mappers |
| Proper pairs only | `-f 2` | Applied in the same samtools step |
| 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 with Bowtie2
bowtie2 --very-sensitive -X 2000 --no-mixed --no-discordant \
  --threads 8 -x GRCh38_bowtie2_index/GRCh38 \
  -1 trimmed_R1.fq.gz -2 trimmed_R2.fq.gz 2> sample.bowtie2.log | \
  samtools view -@ 8 -bS -q 30 -f 2 - | \
  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` is run on this BAM in Stage 3 and published as
`qc/<sample>.idxstats.txt`, which is where the mitochondrial fraction comes from.

## Expected Output
- `aligned/<sample>.bam` + `.bam.bai` -- coordinate-sorted, indexed BAM
- `aligned/<sample>.flagstat.txt` -- alignment summary
- `aligned/<sample>.bowtie2.log` -- Bowtie2 alignment rate summary (parsed by MultiQC)

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Overall alignment rate | >80% | Check genome build, contamination |
| Concordant pair rate | >90% | Check library prep |
| MAPQ>=30 fraction | >60% of mapped | Expected lower than ChIP due to open chromatin |
| Mitochondrial reads | <20% | Measured and removed in Stage 3 |
| Total mapped reads | >=25M after mito removal | May need deeper sequencing |

## Notes
- The `-f 2` flag retains only properly paired reads; combined with `-q 30` this is the
  only BAM flag filtering the workflow performs.
- `-X 2000` allows Bowtie2 to map large fragments from di/tri-nucleosomal DNA.
- Mitochondrial read fraction varies by cell type and preparation; it is measured and
  filtered in Stage 3.

# Stage 3: Filtering, Tn5 Shift, and Fragment Selection

## Tools
- **Samtools 1.19** (image version): Mitochondrial read removal, sorting, indexing
- **Picard MarkDuplicates 3.1.1** (image version, Java 17): PCR duplicate removal
- **deeptools `alignmentSieve` 3.5.5** (image version): Tn5 shift and fragment size selection
- **bedtools 2.31.0** (image version): Blacklist filtering

## Order of operations in the workflow

```
aligned.bam -> mito removal -> Picard MarkDuplicates (REMOVE_DUPLICATES=true)
            -> alignmentSieve --ATACshift -> bedtools blacklist filter
            -> final.bam -> alignmentSieve size selection (NFR, mono-nucleosome)
```

The Tn5 shift runs **after** duplicate removal, not before. Fragment size selection is the
last step, and peaks are called on the NFR BAM.

## Tn5 Transposase Offset Correction

The Tn5 transposase creates a 9-bp target site duplication during insertion. To
accurately represent the cut site, reads must be shifted:
- **Forward strand (+)**: shift +4 bp
- **Reverse strand (-)**: shift -5 bp

This correction is critical for motif footprinting and accurate cut-site analysis.

## Commands

The workflow runs the equivalent of:

```bash
# Step 1: Record the per-contig read counts, then remove chrM reads
samtools idxstats aligned.bam > sample.idxstats.txt
samtools view -@ 4 -b aligned.bam $(samtools idxstats aligned.bam | \
  awk '$1 != "chrM" && $1 != "*" {print $1}' | tr '\n' ' ') > no_mito.bam

# Step 2: Mark and remove PCR duplicates
samtools sort -@ 4 -o sorted.bam no_mito.bam
picard MarkDuplicates \
  INPUT=sorted.bam OUTPUT=dedup.bam \
  METRICS_FILE=dup_metrics.txt \
  REMOVE_DUPLICATES=true VALIDATION_STRINGENCY=LENIENT
samtools index dedup.bam

# Step 3: Apply Tn5 shift (+4/-5)
alignmentSieve --bam dedup.bam --outFile shifted_unsorted.bam \
  --ATACshift --numberOfProcessors 4
samtools sort -@ 4 -o shifted.bam shifted_unsorted.bam
samtools index shifted.bam

# Step 4: Remove reads overlapping blacklist regions
bedtools intersect -v -abam shifted.bam -b hg38-blacklist.v2.bed.gz > final.bam
samtools index final.bam
samtools flagstat final.bam > final.flagstat.txt

# Step 5: Separate nucleosome-free and mono-nucleosomal fragments
alignmentSieve --bam final.bam --outFile nfr_unsorted.bam \
  --maxFragmentLength 150 --numberOfProcessors 4
samtools sort -@ 4 -o nfr.bam nfr_unsorted.bam && samtools index nfr.bam

alignmentSieve --bam final.bam --outFile mononuc_unsorted.bam \
  --minFragmentLength 150 --maxFragmentLength 300 --numberOfProcessors 4
samtools sort -@ 4 -o mononuc.bam mononuc_unsorted.bam
```

The fragment length boundary (150) is `--nfr_max`, and the mitochondrial contig name is
`--mito_name`. There is no `samtools view -F 1804` step in this workflow: flag filtering
(`-q 30 -f 2`) happens once, during alignment.

## Expected Output
- `qc/<sample>.idxstats.txt` -- `samtools idxstats` of the BAM **before** mitochondrial
  reads are removed: one row per contig with name, length, mapped and unmapped counts. The
  mitochondrial fraction is the mapped count on the `--mito_name` row divided by the sum of
  the mapped column; MultiQC's samtools module reads the same file and reports it. The step
  also works when the genome has no such contig -- nothing is removed. No `bc` is involved.
- `filtered/<sample>.dup_metrics.txt` -- Picard duplication metrics
- `filtered/shifted/<sample>.shifted.bam` + `.bai` -- Tn5-corrected, pre-blacklist
- `filtered/<sample>.final.bam` + `.final.bam.bai` -- blacklist-filtered (all fragments)
- `filtered/<sample>.final.flagstat.txt` -- read count after all filtering
- `filtered/nfr/<sample>.nfr.bam` + `.bai` -- nucleosome-free fragments
- `filtered/nfr/<sample>.mononuc.bam` -- mono-nucleosome fragments (no index is written)

Both size-selected BAMs land in `filtered/nfr/`.

## Fragment Size Classes

| Class | Size (bp) | Use | Produced here? |
|-------|-----------|-----|----------------|
| NFR (nucleosome-free) | <150 (`--nfr_max`) | Peak calling, TF footprinting | yes |
| Mono-nucleosome | 150-300 | Nucleosome positioning | yes |
| Di-nucleosome | 300-500 | Chromatin architecture | no |

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Mitochondrial fraction (`qc/<sample>.idxstats.txt`) | <20% (ideal <5%) | Optimize cell lysis |
| Duplication rate | <30% (Picard `PERCENT_DUPLICATION`) | Low complexity library |
| NFR fraction | >40% of fragments <150bp (manual) | Check transposition efficiency |
| Post-filter reads | >=25M | May need deeper sequencing |

## Notes
- `alignmentSieve --ATACshift` applies the +4/-5 offset automatically. It emits an unsorted
  BAM, so every `alignmentSieve` call is followed by `samtools sort`.
- NRF and PBC are not computed by this workflow; see `references/05-qc-metrics.md`.
- Blacklist: Amemiya et al. 2019 (hg38-blacklist.v2.bed.gz, ~900 regions, ~40 Mb). It is
  applied to the BAM, so peaks are already blacklist-clean.

# Stage 4: Peak Calling and IDR

## Tools
- **MACS2 2.2.9.1** (image version): Peak caller (Zhang et al. 2008)
- **IDR 2.0.4.2** (image version): Irreproducible Discovery Rate (Li et al. 2011)
- **bedtools 2.31.0 + samtools 1.19** (image versions): FRiP counting after peak calling

## MACS2 Parameters for ATAC-seq

ATAC-seq peak calling differs from ChIP-seq in several key ways:
1. **No control/input file** -- peaks are called against the local background model
2. **Use the NFR BAM only** -- nucleosome-free fragments for accessibility peaks
3. **BAMPE mode** -- fragment coordinates come from the read pairs

| Parameter | Value | Notes |
|-----------|-------|-------|
| `--format` | BAMPE | Use actual fragment sizes |
| `--gsize` | hs (2.7e9) | `mm` when `--genome mm10` |
| `--nomodel` | yes | Do not build a shifting model |
| `--qvalue` | 0.05 | FDR threshold |
| `--keep-dup` | all | Duplicates were already removed in Stage 3 |
| `--call-summits` | yes | Identify sub-peak summits |
| `-B` | yes | Generate bedGraphs |

**`--shift` and `--extsize` are not used, and would have no effect here.** In `-f BAMPE`
mode MACS2 takes fragment coordinates from the read pairs: it forces `nomodel = True` and
sets `shift = 0` internally, without warning. Shift/extension values such as
`--shift -100 --extsize 200` or `--shift -75 --extsize 150` apply only when calling peaks
on BED or single-end input, which this workflow does not do.

## Commands

The workflow runs the equivalent of:

```bash
# Peak calling on NFR fragments
macs2 callpeak -t sample.nfr.bam \
  -f BAMPE -g hs -n sample \
  --nomodel --keep-dup all --call-summits \
  --qvalue 0.05 -B

# IDR, once per pair of replicates, the two names sorted alphabetically
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
```

Calling peaks on the blacklist-filtered all-fragment BAM (`filtered/<sample>.final.bam`)
instead of the NFR BAM is a reasonable manual variant, but the workflow always uses NFR.

## What the workflow publishes

- `peaks/narrow/<sample>_peaks.narrowPeak`
- `peaks/narrow/<sample>_peaks.xls`
- `peaks/narrow/<sample>_treat_pileup.bdg` and `<sample>_control_lambda.bdg`
- `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/`
alongside the peak files. Only narrow peaks are called; there is no broad mode.

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

- 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 and `peaks/idr/` is
  not created.
- There are **no** pooled-replicate calls, **no** pseudoreplicates, and **no** optimal /
  conservative peak sets.
- Rescue ratio and self-consistency ratio are **not computed**; they require pooled and
  pseudoreplicated peak calls that this workflow does not produce.

`--skip_idr` turns the step off entirely.

## IDR Interpretation

| Metric | Expected (Good) | Concern | Computed here? |
|--------|-----------------|---------|----------------|
| IDR peaks (0.05 threshold) | 50,000-150,000 | <30,000 suggests poor signal | yes |
| 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.3 (ATAC-seq standard) | Poor accessibility signal |
| Peak count | >50,000 (IDR filtered) | Low enrichment |
| Peak width distribution | Median 200-500 bp | Check if calling mode correct |
| Peaks at TSS (manual) | Enrichment visible | Fundamental ATAC-seq signal |

## Notes

- ATAC-seq FRiP is typically much higher than ChIP-seq (0.3-0.6 vs 0.01-0.1)
  because open chromatin is a large fraction of the genome. The workflow computes it after
  peak calling, counting the peaks called here against the blacklist-filtered all-fragment
  BAM; see `references/05-qc-metrics.md`.
- The workflow always calls peaks on NFR fragments for accessibility analysis.
- For nucleosome positioning, use `filtered/nfr/<sample>.mononuc.bam` separately.
- IDR is standard for ATAC-seq with biological replicates.

# Stage 5: Signal Tracks and QC Report

Stage 5 of the workflow does three things: it builds one bigWig signal track per sample, it
computes FRiP, 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)

The workflow runs `bamCoverage` once per sample, on the blacklist-filtered BAM
(`filtered/<sample>.final.bam`, all fragments), and publishes
`signal/<sample>.signal.bw`:

```bash
bamCoverage -b final.bam -o sample.signal.bw \
  --normalizeUsing RPKM --binSize 10 \
  --numberOfProcessors 4 --extendReads
```

There is no NFR-only signal track. If you want one, run the same command manually against
`filtered/nfr/<sample>.nfr.bam`:

```bash
bamCoverage -b results/filtered/nfr/sample.nfr.bam -o sample_nfr_signal.bw \
  --normalizeUsing RPKM --binSize 10 \
  --numberOfProcessors 8 --extendReads
```

## FRiP (workflow)

The workflow computes the fraction of reads in peaks for every sample and publishes
`qc/<sample>.frip_mqc.tsv`. The denominator is the blacklist-filtered BAM (all fragments),
and the peaks are the ones called on the nucleosome-free fragments:

```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`; no `bc` is involved. 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 ATAC-seq expectation of >=0.3 is read
directly off the `FRiP` column. To score the nucleosome-free fragments alone instead, run
the same commands against `results/filtered/nfr/sample.nfr.bam`.

## 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
- The Bowtie2 alignment log
- `qc/<sample>.idxstats.txt`, read by MultiQC's samtools module, which reports the
  mitochondrial fraction from it
- 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 | >=50M recommended | yes |
| Mapping rate | Bowtie2 log, samtools flagstat | >80% | yes |
| Mitochondrial fraction | `qc/<sample>.idxstats.txt` | <20% | yes |
| Duplication rate | Picard | <30% | yes |
| IDR peaks | IDR | >50,000 | yes (one file per replicate pair) |
| TSS enrichment | deeptools + TSS BED | >=5 | no |
| FRiP | bedtools + samtools (`qc/<sample>.frip_mqc.tsv`) | >=0.3 | yes |
| NFR fraction | fragment size distribution | >40% | no |
| NRF / PBC1 | pre-dedup BAM | >=0.8 | no |

## Manual step: TSS Enrichment Score

The most important ATAC-seq quality metric, and **not computed by this workflow**: there is
no TSS BED parameter and no `computeMatrix`/`plotProfile` step. deeptools 3.5.5 is in the
image, so this can be run against the published bigWig with a TSS BED you supply (for
example, a GENCODE TSS BED for your assembly):

```bash
computeMatrix reference-point -S results/signal/sample.signal.bw \
  -R tss.bed -a 2000 -b 2000 \
  --referencePoint TSS -o tss_matrix.gz

plotProfile -m tss_matrix.gz -o tss_enrichment.pdf \
  --perGroup --refPointLabel TSS

# TSS enrichment score = max(TSS signal) / mean(flanking signal)
```

| TSS Score | Quality | ENCODE Standard |
|-----------|---------|-----------------|
| >=7 | Excellent | Pass |
| 5-7 | Good | Pass (ENCODE minimum = 5 for GRCh38) |
| 3-5 | Marginal | Borderline |
| <3 | Poor | Fail |

## Manual step: Fragment Size Distribution

**Not computed or plotted by this workflow.** deeptools is in the image; Picard
`CollectInsertSizeMetrics` also needs R, which the image does not have, so prefer
`bamPEFragmentSize`:

```bash
bamPEFragmentSize -b results/filtered/sample.final.bam -o fragment_sizes.pdf \
  --maxFragmentLength 1000 --numberOfProcessors 4
```

Expected pattern: peaks at ~200 bp (NFR), ~400 bp (mono-nuc), ~600 bp (di-nuc).

## Manual step: ataqv (ATAC-seq QC)

Comprehensive ATAC-seq-specific QC tool from Parker Lab (Orchard et al. 2020). **ataqv is
not in the pipeline image, not in `atacseq-env.yml`, and is never invoked by the
workflow.** Install it separately, for example with `conda install -c bioconda ataqv`:

```bash
ataqv --peak-file results/peaks/narrow/sample_peaks.narrowPeak --tss-file tss.bed \
  --name sample_name --metrics-file sample.ataqv.json \
  human results/filtered/sample.final.bam

# Generate HTML report
mkarv my_ataqv_report/ sample.ataqv.json
```

ataqv reports TSS enrichment, fragment length distribution, peak metrics,
mitochondrial fraction, and duplicate rate in a single interactive report.

# ATAC-seq Pipeline — Literature References

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

---

## ATAC-seq Method Development

These papers establish the ATAC-seq assay and its key methodological improvements.

---

### Buenrostro et al. 2013 — Transposition of native chromatin for ATAC-seq

- **Citation:** Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nature Methods, 10(12):1213-1218, 2013.
- **DOI:** [10.1038/nmeth.2688](https://doi.org/10.1038/nmeth.2688)
- **PMID:** 24097267 | **PMC:** PMC3959825
- **Citations:** ~5,600
- **Key findings:** Original ATAC-seq method paper. Described the use of the hyperactive Tn5 transposase to simultaneously fragment DNA and insert sequencing adapters into accessible chromatin regions in a single enzymatic step. The method requires only 500-50,000 cells and a simple two-step protocol (transposition + PCR amplification), making it substantially faster and simpler than DNase-seq or FAIRE-seq. Demonstrated that ATAC-seq captures open chromatin sites, DNA-binding protein footprints, individual nucleosome positions, and chromatin compaction at nucleotide resolution. Identified classes of DNA-binding factors that strictly avoided, could tolerate, or tended to overlap with nucleosomes. The characteristic nucleosomal ladder pattern in fragment size distribution (sub-nucleosomal < 150 bp, mono-nucleosome 150-300 bp, di-nucleosome 300-500 bp) serves as a key QC metric in this pipeline.

---

### Corces et al. 2017 — Omni-ATAC: Improved ATAC-seq protocol

- **Citation:** Corces MR, Trevino AE, Hamilton EG, Greenside PG, Sinnott-Armstrong NA, Vesuna S, Satpathy AT, Rubin AJ, Montine KS, Wu B, Kathiria A, Cho SW, Mumbach MR, Carter AC, Kasowski M, Orloff LA, Risca VI, Kundaje A, Khavari PA, Montine TJ, Greenleaf WJ, Chang HY. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nature Methods, 14(10):959-962, 2017.
- **DOI:** [10.1038/nmeth.4396](https://doi.org/10.1038/nmeth.4396)
- **PMID:** 28846090 | **PMC:** PMC5623106
- **Citations:** ~1,925
- **Key findings:** Introduced the Omni-ATAC protocol with three key improvements over the original method: (1) use of digitonin-based lysis with NP-40 and Tween-20 detergents to reduce mitochondrial DNA contamination, (2) ability to process archival frozen tissue samples including 50-micron sections, and (3) substantial improvement in signal-to-background ratio. The reduced mitochondrial contamination directly decreases sequencing cost by redirecting reads to informative nuclear DNA. Omni-ATAC enabled the first chromatin accessibility maps from distinct human brain structures, revealing activities of disease-associated DNA elements. This protocol is now the basis for most ENCODE ATAC-seq experiments and informs the aggressive mitochondrial filtering in this pipeline.

---

### Grandi et al. 2022 — Chromatin accessibility profiling by ATAC-seq (Nature Protocols)

- **Citation:** Grandi F, Modi H, Kampman L, Corces MR. Chromatin accessibility profiling by ATAC-seq. Nature Protocols, 17(6):1518-1552, 2022.
- **DOI:** [10.1038/s41596-022-00692-9](https://doi.org/10.1038/s41596-022-00692-9)
- **PMID:** 35478247
- **Citations:** ~287
- **Key findings:** Comprehensive step-by-step protocol for the Omni-ATAC workflow covering sample preparation, transposition, library preparation, sequencing, and data analysis. Documents that ATAC-seq libraries for 12 samples can be generated in approximately 10 hours, with downstream computational analysis implementable using benchmarked pipelines. Provides detailed recommendations for troubleshooting common issues including high mitochondrial fraction, low library complexity, and poor TSS enrichment. Serves as the current definitive experimental protocol reference for ATAC-seq.

---

## ENCODE Pipeline Standards

---

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

- **Citation:** ENCODE Project Consortium et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature, 583(7818):699-710, 2020.
- **DOI:** [10.1038/s41586-020-2493-4](https://doi.org/10.1038/s41586-020-2493-4)
- **PMID:** 32728249 | **PMC:** PMC7410828
- **Citations:** ~2,500
- **Key findings:** ENCODE Phase 3 paper establishing the uniform processing framework for all ENCODE assays including ATAC-seq. Defined the Registry of candidate cis-Regulatory Elements (cCREs) using integrated ATAC-seq and histone ChIP-seq data. ATAC-seq is a primary input for identifying open chromatin elements in the cCRE classification system, alongside DNase-seq and H3K4me3/H3K27ac ChIP-seq. Established assembly standards (GRCh38/mm10), IDR-based replicate consistency assessment, and preferred_default file selection criteria used in this pipeline.

---

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

- **Citation:** Hitz BC, Lee JW, Jolanki O, et al. The ENCODE Uniform Analysis Pipelines. Research Square (preprint), 2023.
- **DOI:** [10.21203/rs.3.rs-311127/v1](https://doi.org/10.21203/rs.3.rs-311127/v1)
- **Citations:** ~84
- **Key findings:** Documents the official ENCODE ATAC-seq pipeline specification including Bowtie2 alignment (vs BWA-MEM for ChIP-seq), Tn5 offset correction (+4/-5 bp), mitochondrial read removal, nucleosome-free fragment selection (<150 bp), MACS2 peak calling without input control, and IDR-based replicate consistency. Specifies the pipeline infrastructure using Docker containers and WDL workflow language, enabling execution on local machines, HPC clusters, or cloud environments. This paper is the authoritative reference for the pipeline architecture implemented in this skill.

---

## ATAC-seq Analysis Best Practices

---

### Yan et al. 2020 — From reads to insight: a hitchhiker's guide to ATAC-seq data analysis

- **Citation:** Yan F, Powell DR, Curtis DJ, Wong NC. From reads to insight: a hitchhiker's guide to ATAC-seq data analysis. Genome Biology, 21(1):22, 2020.
- **DOI:** [10.1186/s13059-020-1929-3](https://doi.org/10.1186/s13059-020-1929-3)
- **PMID:** 32014034 | **PMC:** PMC6996192
- **Citations:** ~328
- **Key findings:** Comprehensive review of ATAC-seq computational analysis covering pre-analysis (quality check, alignment), core analysis (peak calling), and advanced analysis (differential accessibility, motif enrichment, footprinting, nucleosome positioning). Discussed the critical importance of the Tn5 offset correction for accurate TF footprinting and the distinction between nucleosome-free region (NFR) analysis and full-fragment analysis. Highlighted that MACS2 should be run with --nomodel --shift -100 --extsize 200 for ATAC-seq (or alternatively on NFR fragments only) and that peak calling should use the --keep-dup all flag to retain all unique fragments. Reviewed the reconstruction of transcriptional regulatory networks with multi-omics data integration.

---

### Orchard et al. 2020 — ataqv: quality control for ATAC-seq data

- **Citation:** Orchard P, Kyono Y, Hensley J, Kitzman JO, Parker SCJ. Quantification, Dynamic Visualization, and Validation of Bias in ATAC-Seq Data with ataqv. Cell Systems, 10(3):298-306.e4, 2020.
- **DOI:** [10.1016/j.cels.2020.02.009](https://doi.org/10.1016/j.cels.2020.02.009)
- **PMID:** 32213349 | **PMC:** PMC7138743
- **Citations:** ~62
- **Key findings:** Introduced ataqv, a computational toolkit for ATAC-seq quality control that computes, visualizes, and compares QC metrics across samples. Analyzed 2,009 public ATAC-seq datasets revealing a 10-fold range in QC metrics. Demonstrated through Tn5 dosage experiments that technical variation in the Tn5-to-nuclei ratio and sequencing flowcell density induces systematic bias by changing enrichment of reads across functional genomic annotations including promoters, enhancers, and TF-bound regions (with CTCF as a notable exception). Established TSS enrichment as a primary quality indicator and showed that the Tn5 dosage effect is the dominant source of technical bias in ATAC-seq. ataqv is used in this pipeline's QC stage alongside deepTools and MultiQC.

---

## Core Pipeline Tools

---

### Langmead & Salzberg 2012 — Bowtie 2

- **Citation:** Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods, 9(4):357-359, 2012.
- **DOI:** [10.1038/nmeth.1923](https://doi.org/10.1038/nmeth.1923)
- **PMID:** 22388286 | **PMC:** PMC3322381
- **Citations:** ~47,300
- **Key findings:** Introduced Bowtie 2, which combines a full-text minute index with hardware-accelerated dynamic programming for fast, sensitive, and accurate gapped alignment. Unlike BWA-MEM (used for ChIP-seq), Bowtie 2 is particularly well-suited for ATAC-seq because of its superior handling of short fragments and the very short reads produced by nucleosome-free regions. The ENCODE ATAC-seq pipeline uses Bowtie 2 with --very-sensitive mode and maximum fragment length set to accommodate the nucleosomal ladder (typically --maxins 2000). Bowtie 2 is used in this pipeline's alignment stage.

---

### Zhang et al. 2008 — MACS peak caller

- **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:** MACS2 is the ENCODE-standard peak caller for both ChIP-seq and ATAC-seq. For ATAC-seq, MACS2 is run differently than for ChIP-seq: no input control is used (peaks are called against the local background model), the --nomodel flag is set because the Tn5 insertion creates a distinct fragment distribution, and the --shift and --extsize parameters are adjusted for the Tn5 offset. Peak calling is performed only on nucleosome-free fragments (<150 bp) with --call-summits for precise accessible region identification. The narrow peak mode is always used for ATAC-seq (unlike ChIP-seq which uses broad mode for some histone marks).

---

### Li et al. 2011 — IDR (Irreproducible Discovery Rate)

- **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:** IDR framework for replicate consistency assessment. In ATAC-seq, IDR is applied to MACS2 peak calls from biological replicates to produce a consensus peak set at IDR threshold 0.05. The same methodology as ChIP-seq but particularly important for ATAC-seq because the open chromatin signal is typically broader and less punctate than TF ChIP-seq peaks, making replicate concordance assessment critical for distinguishing true accessible regions from noise.

---

## Supplementary Tools (Non-ATAC-Specific)

These tools are shared with the ChIP-seq pipeline and other genomics workflows. See pipeline-chipseq/references/literature.md for detailed descriptions.

---

### Li et al. 2009 — SAMtools

- **Citation:** Li H, et al. 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 | **Citations:** ~53,700
- **ATAC-seq role:** Used for coordinate sorting, MAPQ filtering, mitochondrial read removal (samtools view filtering chrM), fragment size filtering for NFR selection, and BAM indexing. The mitochondrial filtering step is especially critical for ATAC-seq where mitochondrial reads can constitute 30-80% of total reads.

---

### Amemiya et al. 2019 — 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 | **Citations:** ~1,372
- **ATAC-seq role:** Blacklist v2 filtering is applied after Tn5 offset correction and before peak calling. Essential for ATAC-seq because Tn5 can insert preferentially into certain repetitive regions, creating artifact peaks that overlap blacklisted regions.

---

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

- **Citation:** Ramírez F, et al. 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 | **Citations:** ~6,100
- **ATAC-seq role:** Generates normalized signal tracks (bigWig), fingerprint plots, and TSS enrichment heatmaps. The TSS enrichment score (computed via computeMatrix + plotProfile at RefSeq TSS annotations) is the primary ENCODE QC metric for ATAC-seq, with assembly-dependent thresholds: >= 5 GRCh38 / >= 6 hg19 / >= 10 mm10 (ENCODE data standards).

---

### 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 | **Citations:** ~6,800
- **ATAC-seq role:** Aggregates QC metrics from FastQC, Bowtie2, Picard, ataqv, and MACS2 into a unified HTML report for batch-level quality assessment.

---

### Broad Institute — Picard MarkDuplicates

- **Citation:** Broad Institute. Picard toolkit. GitHub, 2019.
- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **ATAC-seq role:** PCR duplicate marking and library complexity estimation. Particularly important for ATAC-seq where low-input protocols may have higher duplication rates. NRF, PBC1, PBC2 metrics are computed from Picard output.

---

### 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 | **Citations:** ~12,000
- **ATAC-seq role:** Used for blacklist region removal, fragment size filtering, FRiP calculation, and peak annotation with gene models.

---

### 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
- **ATAC-seq role:** Trim Galore (Cutadapt wrapper) performs adapter trimming. ATAC-seq libraries frequently contain adapter sequences because nucleosome-free fragments are short (< 150 bp) and paired-end reads can read through the insert into the adapter on the other end.

---

### 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/)
- **ATAC-seq role:** Raw read quality assessment. Note that ATAC-seq libraries often trigger FastQC warnings for per-sequence GC content (due to Tn5 insertion bias) and sequence duplication levels (expected for low-input libraries) — these warnings are typically not concerning for ATAC-seq.