ENCODE ATAC-seq Pipeline

When to Use

• User wants to run an ATAC-seq processing pipeline from FASTQ to peaks and signal tracks
• User asks about "ATAC-seq pipeline", "Tn5 shift", "chromatin accessibility pipeline", or "Bowtie2 for ATAC"
• User needs to process ATAC-seq data with proper Tn5 insertion site correction
• Example queries: "process my ATAC-seq FASTQs", "run ENCODE ATAC-seq pipeline", "call accessibility peaks from ATAC-seq"

Execute the ENCODE ATAC-seq processing pipeline from raw FASTQ files through Tn5 offset correction, peak calling, IDR analysis, and signal track generation. This skill provides a Nextflow DSL2 implementation following ENCODE uniform analysis standards.

TSS enrichment scoring is a manual post-processing step; the workflow does not compute it (see "Manual QC steps" below and references/05-qc-metrics.md).

Overview

ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) uses the Tn5 transposase to probe open chromatin regions. This pipeline processes ATAC-seq data through quality control, alignment with Bowtie2, mitochondrial read removal, duplicate removal, Tn5 insertion site correction (+4/-5 bp offset), blacklist filtering, nucleosome-free fragment selection, MACS2 peak calling, FRiP calculation, and an IDR comparison for every pair of replicates.

Key differences from ChIP-seq: Bowtie2 aligner (optimized for short fragments), Tn5 transposase shift correction, mitochondrial read filtering (chrM can be 30-80% of reads), and no input control.

The workflow is paired-end only. Passing --single_end stops the run with an error, because Tn5 shifting, nucleosome-free selection and BAMPE peak calling all depend on fragment length.

Key Literature

Pipeline Stages

FASTQ ──> FastQC / Trim Galore ──> Bowtie2 ──> Mito Removal ──> Picard MarkDuplicates
  │                                            (chrM dropped)   (duplicates REMOVED)
  │                                                                       │
  │           ┌───────────────────────────────────────────────────────────┘
  │           v
  │     Tn5 Shift (alignmentSieve --ATACshift) ──> Blacklist Filter ──> Size Selection
  │                                                       │                    │
  │                                                       │        ┌───────────┴────────┐
  │                                                       v        v                    v
  │                                               Signal Track   NFR (<150 bp)   Mono-nucleosome
  │                                                (all frags)      │              (150-300 bp)
  │                                                                 v
  │                                                   MACS2 Peak Calling ──> IDR (every pair)
  │                                                                 │
  │                                                                 v
  │                                                   FRiP (NFR peaks vs final BAM)
  v
 QC reports ────────────────────────────────────────────────────────────────> MultiQC

The Tn5 shift runs after duplicate removal, and peaks are called on the nucleosome-free BAM only. The signal track is built from all fragments in the blacklist-filtered BAM, not from the NFR BAM.

Stage Summary

Input Requirements

Required

• ATAC-seq FASTQ (--reads): paired-end reads, gzipped. A Nextflow file-pair glob,

e.g. 'fastq/*_R{1,2}.fq.gz'.

• Bowtie2 index directory (--bowtie2_index): Bowtie2 is invoked as

bowtie2 ... -x <dir>/<genome>, so the directory must hold index files named after the 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 with bowtie2-build GRCh38.fa GRCh38_bowtie2_index/GRCh38. If the flag is omitted, the workflow looks for ./<genome>_bowtie2_index in the launch directory. The workflow does not build or download the index.

Optional

• Blacklist (--blacklist): defaults to the ENCODE Blacklist v2 URL for --genome.

There is no sample sheet and no input control. Inputs are globs, and every sample in a run shares one --genome. Unlike ChIP-seq, ATAC-seq does not need a separate input or IgG control; MACS2 calls peaks against a local background model.

Tn5 Transposase Offset Correction

The Tn5 transposase inserts sequencing adapters with a 9-bp duplication. To center reads on the actual cut site:

• Forward strand (+): shift +4 bp
• Reverse strand (-): shift -5 bp

The workflow applies this with alignmentSieve --ATACshift (deeptools) after duplicate removal and before blacklist filtering. The correction is essential for footprinting and motif analysis.

Fragment Size Distribution

ATAC-seq produces a characteristic nucleosomal ladder pattern:

The workflow calls peaks on the nucleosome-free BAM. The NFR/mono-nucleosome boundary is --nfr_max (default 150); the mono-nucleosome selection is --nfr_max to 300 bp. The workflow does not plot the fragment size distribution.

Parameters

Pipeline parameters (main.nf)

Infrastructure parameters (nextflow.config)

Profiles are local, slurm, gcp and aws.

QC Thresholds

The workflow computes only the metrics marked "workflow" below. TSS enrichment, NRF/PBC, fragment-size plots and ataqv are manual post-processing steps documented in references/05-qc-metrics.md.

qc/<sample>.idxstats.txt is samtools idxstats of the BAM before mitochondrial reads are removed (contig, length, mapped, unmapped): 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 that fraction.

TSS Enrichment Score (manual)

The TSS enrichment score measures the fold enrichment of ATAC-seq signal at transcription start sites compared to flanking regions. It is the single most informative QC metric for ATAC-seq, but this workflow does not compute it: there is no TSS BED input and no computeMatrix/plotProfile step. Run it manually against signal/<sample>.signal.bw with a TSS BED for your assembly; the commands are in references/05-qc-metrics.md.

Execution

Quick Start (Local Docker)

nextflow run scripts/main.nf \
  -profile local \
  --reads 'fastq/*_R{1,2}.fq.gz' \
  --genome GRCh38 \
  --bowtie2_index GRCh38_bowtie2_index \
  --blacklist hg38-blacklist.v2.bed.gz \
  --outdir results/

--blacklist is optional; without it the workflow downloads the ENCODE Blacklist v2 for --genome. Give the glob at least two replicates if you want the IDR step to run.

SLURM HPC

The slurm profile runs through Singularity, so pass a local image file rather than the default Docker image name:

singularity build pipeline-atacseq.sif docker-daemon://encode-toolkit/pipeline-atacseq:1.0.0

nextflow run scripts/main.nf \
  -profile slurm \
  --container /path/to/pipeline-atacseq.sif \
  --slurm_queue normal \
  --reads 'fastq/*_R{1,2}.fq.gz' \
  --genome GRCh38 \
  --bowtie2_index GRCh38_bowtie2_index \
  --outdir results/

Cloud

# Google Cloud Batch
nextflow run scripts/main.nf -profile gcp \
    --container us-docker.pkg.dev/<project>/<repo>/pipeline-atacseq:1.0.0 \
    --gcp_project <project> \
    --gcp_workdir gs://<bucket>/work \
    --reads 'gs://<bucket>/fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --bowtie2_index gs://<bucket>/reference/GRCh38_bowtie2_index \
    --outdir gs://<bucket>/results

# AWS Batch
nextflow run scripts/main.nf -profile aws \
    --container <account>.dkr.ecr.<region>.amazonaws.com/pipeline-atacseq:1.0.0 \
    --aws_queue <job-queue> \
    --aws_workdir s3://<bucket>/work \
    --reads 's3://<bucket>/fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --bowtie2_index s3://<bucket>/reference/GRCh38_bowtie2_index \
    --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.

Cloud Cost Estimates

Output Directory Structure

results/
  fastqc/                   # FastQC reports for raw and trimmed reads (.html, .zip)
  trimmed/                  # Trimmed FASTQ (*_val_1.fq.gz / *_val_2.fq.gz) + trimming reports
  aligned/                  # <sample>.bam, .bam.bai, <sample>.flagstat.txt,
                            #   <sample>.bowtie2.log
  filtered/                 # <sample>.dup_metrics.txt, <sample>.final.bam(.bai),
                            #   <sample>.final.flagstat.txt
    shifted/                # <sample>.shifted.bam(.bai) -- Tn5-corrected, pre-blacklist
    nfr/                    # <sample>.nfr.bam(.bai) and <sample>.mononuc.bam
  peaks/
    narrow/                 # <sample>_peaks.narrowPeak, _summits.bed, _peaks.xls,
                            #   _treat_pileup.bdg, _control_lambda.bdg
    idr/                    # <sampleA>_vs_<sampleB>.idr_peaks.txt (+ .png), one file per
                            #   replicate pair
  signal/                   # <sample>.signal.bw (all fragments, RPKM)
  qc/
    <sample>.idxstats.txt   # samtools idxstats before chrM removal:
                            #   contig / length / mapped / unmapped
    <sample>.frip_mqc.tsv   # Peak set / FRiP / reads_in_peaks / total_reads
    multiqc/                # multiqc_report.html, multiqc_data/
  pipeline_info/            # timeline.html, report.html, trace.txt

The nucleosome-free and mono-nucleosome BAMs are both written to filtered/nfr/; there is no filtered/mononuc/ directory, and mononuc.bam has no index.

Common Pitfalls

1. High Mitochondrial Read Fraction

Mitochondrial DNA lacks chromatin and is highly accessible, often capturing 30-80% of reads. This is the most common ATAC-seq quality issue. The workflow removes --mito_name reads and publishes the per-contig counts it used, qc/<sample>.idxstats.txt: divide the mapped count on the --mito_name row by the sum of the mapped column to get the fraction, or read it from the samtools section of the MultiQC report. If >50% mito, consider optimizing the cell lysis step.

2. Wrong mitochondrial contig name

--mito_name defaults to chrM. Assemblies that call the contig MT need --mito_name MT. Confirm the name with samtools idxstats on a BAM from your index before running.

3. Using BWA Instead of Bowtie2

Bowtie2 handles the short fragments from ATAC-seq (especially NFR <150bp) better than BWA-MEM. The workflow uses Bowtie2 with --very-sensitive.

4. Only one replicate in the --reads glob

IDR needs a pair of samples. With one sample there is no pair, so IDR is skipped silently and peaks/idr/ is never created. Every sample matched by --reads is treated as a replicate of the same experiment, so unrelated samples in one glob produce meaningless pairwise comparisons.

5. TSS enrichment is not in the output

TSS enrichment is the most informative single metric for ATAC-seq, but the workflow does not compute it. Run the manual computeMatrix/plotProfile step in references/05-qc-metrics.md before judging a library.

Pipeline Scripts

The image is pinned to linux/amd64; on an arm64 host it runs under emulation.

Tool versions in the image: Bowtie2 2.5.4, samtools 1.19, bedtools 2.31.0, Picard 3.1.1 (Java 17), Trim Galore 0.6.10 with cutadapt 4.6, FastQC 0.12.1, MACS2 2.2.9.1, IDR 2.0.4.2, deepTools 3.5.5, MultiQC 1.21. The conda environment bioinformatics-installer/environments/atacseq-env.yml pins the same versions of the tools it lists.

ENCODE Data Integration

After running on your own data, compare with ENCODE reference:

# Find matching ENCODE ATAC-seq experiments
encode_search_experiments(
    assay_title="ATAC-seq",
    organ="pancreas",
    biosample_type="tissue"
)

# Download ENCODE peaks for comparison
encode_batch_download(
    download_dir="/data/encode_reference/",
    output_type="IDR thresholded peaks",
    assay_title="ATAC-seq",
    organ="pancreas",
    assembly="GRCh38"
)

Pitfalls & Edge Cases

• Tn5 shift is critical: ATAC-seq reads must be shifted +4/-5 bp to center on the Tn5

insertion site. The workflow does this with alignmentSieve --ATACshift after duplicate removal. Without the correction, footprinting is offset by ~5 bp.

• Mitochondrial reads dominate: Expect 30-80% mitochondrial reads. The workflow drops

them right after alignment, before duplicate removal and peak calling. >80% chrM indicates dead/dying cells or poor nuclei isolation.

• Fragment size distribution is diagnostic: a nucleosomal ladder (sub-nucleosomal

<150bp, mono-nucleosomal ~200bp, di-nucleosomal ~400bp) confirms successful transposition. The workflow does not plot it; use bamPEFragmentSize manually.

• TSS enrichment threshold: ENCODE requires TSS enrichment >=5 (GRCh38), >=6 (hg19), or

>=10 (mm10) for ATAC-seq (ENCODE data standards). Values below 4 indicate poor signal-to-noise. Computed manually, not by this workflow.

• MACS2 `--shift`/`--extsize` do not apply here: the workflow calls peaks in -f BAMPE

mode, where MACS2 takes fragment coordinates from read pairs, forces --nomodel and neutralises --shift internally. Shift/extension values only matter when calling peaks on BED or single-end input. See references/04-peak-calling.md.

• Paired-end only: --single_end is rejected with an error. Single-end ATAC-seq cannot

distinguish nucleosome-free from nucleosomal fragments.

Walkthrough: Processing ENCODE ATAC-seq from FASTQ to Accessible Chromatin Peaks

Goal: Process raw ATAC-seq FASTQ files through this pipeline to generate nucleosome-free region peaks and a signal track. Context: Bowtie2 alignment, chrM removal, duplicate removal, Tn5 shift (+4/-5), blacklist filtering, NFR selection and MACS2 peak calling.

Step 1: Find ATAC-seq experiment

encode_get_experiment(accession="ENCSR637ENO")

Expected output (fields abridged):

{
  "accession": "ENCSR637ENO",
  "assay_title": "ATAC-seq",
  "biosample_summary": "GM12878",
  "assembly": ["GRCh38"],
  "bio_replicate_count": 2,
  "tech_replicate_count": 2,
  "status": "released"
}

Step 2: List FASTQ files

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

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

[
  {"accession": "ENCFF100ATQ", "file_format": "fastq", "output_type": "reads", "biological_replicates": [1], "file_size_human": "1.8 GB"},
  {"accession": "ENCFF101ATQ", "file_format": "fastq", "output_type": "reads", "biological_replicates": [1], "file_size_human": "1.9 GB"},
  {"accession": "ENCFF102ATQ", "file_format": "fastq", "output_type": "reads", "biological_replicates": [2], "file_size_human": "1.7 GB"},
  {"accession": "ENCFF103ATQ", "file_format": "fastq", "output_type": "reads", "biological_replicates": [2], "file_size_human": "1.8 GB"}
]

The listing does not say which file of a pair is read 1 and which is read 2 -- no encode_* tool reports that. Open each file's page on encodeproject.org, where paired_end is 1 or 2 and paired_with names the other accession. Both replicates are needed for the IDR step.

Step 3: Download and name the FASTQs so a read-pair glob can find them

encode_download_files(file_accessions=["ENCFF100ATQ", "ENCFF101ATQ", "ENCFF102ATQ", "ENCFF103ATQ"], download_dir="/data/atacseq/fastq")

ENCODE FASTQs are named by accession (ENCFF123ABC.fastq.gz) with no _R1/_R2 in the name, so the --reads glob (*_R{1,2}.fq.gz) cannot pair them. Take the mate assignment from each file's page on encodeproject.org (paired_end is 1 or 2, paired_with names the other accession), then link them into the shape the glob expects:

cd /data/atacseq/fastq
ln -s ENCFF100ATQ.fastq.gz gm12878_rep1_R1.fq.gz
ln -s ENCFF101ATQ.fastq.gz gm12878_rep1_R2.fq.gz
ln -s ENCFF102ATQ.fastq.gz gm12878_rep2_R1.fq.gz
ln -s ENCFF103ATQ.fastq.gz gm12878_rep2_R2.fq.gz

Step 4: Run the ATAC-seq pipeline

nextflow run scripts/main.nf \
  -profile local \
  --reads '/data/atacseq/fastq/gm12878_*_R{1,2}.fq.gz' \
  --genome GRCh38 \
  --bowtie2_index /data/reference/GRCh38_bowtie2_index \
  --blacklist /data/reference/hg38-blacklist.v2.bed.gz \
  --mito_name chrM \
  --outdir /data/atacseq/results

Pipeline steps, in the order the workflow runs them:

1. FastQC on raw reads
2. Adapter trimming with Trim Galore (--nextera), plus FastQC on the trimmed reads
3. Alignment (Bowtie2 --very-sensitive, MAPQ 30, properly paired only)
4. Mitochondrial read removal, with samtools idxstats of the pre-removal BAM published as qc/<sample>.idxstats.txt
5. Duplicate removal (Picard REMOVE_DUPLICATES=true)
6. Tn5 shift correction (+4/-5, alignmentSieve --ATACshift)
7. Blacklist filtering of the BAM
8. Nucleosome-free (<150 bp) and mono-nucleosome (150-300 bp) selection
9. Peak calling on the NFR BAM (MACS2 -f BAMPE --nomodel --keep-dup all --call-summits --qvalue 0.05 -B)
10. IDR on every pair of replicates, signal track from all fragments, FRiP, MultiQC

Step 5: Validate output quality

From the workflow:

With the two replicates above there is one IDR file; a third replicate would add gm12878_rep1_vs_gm12878_rep3 and gm12878_rep2_vs_gm12878_rep3.

Manual follow-ups (not run by this workflow): TSS enrichment, fragment-size distribution plots, NRF/PBC and ataqv. Commands are in references/05-qc-metrics.md.

Step 6: Track and log provenance

encode_track_experiment(accession="ENCSR637ENO", notes="GM12878 ATAC-seq processed through the pipeline-atacseq skill")

Integration with downstream skills

• Accessible chromatin peaks feed into -> accessibility-aggregation for cross-experiment union merge
• Peak regions feed into -> motif-analysis for TF motif enrichment
• Signal tracks feed into -> visualization-workflow for browser display
• Peaks feed into -> regulatory-elements for cCRE classification
• QC metrics validated by -> quality-assessment

Code Examples

1. Find ATAC-seq data for processing

encode_search_experiments(
  assay_title="ATAC-seq",
  organ="pancreas"
)

Expected output:

{
  "results": [
    {
      "accession": "ENCSR789PAN",
      "assay_title": "ATAC-seq",
      "biosample_summary": "pancreas tissue male adult (44 years)",
      "status": "released"
    }
  ],
  "total": 8,
  "limit": 25,
  "offset": 0,
  "has_more": false,
  "next_offset": null
}

2. Check file details before download

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

Expected output:

[
  {
    "accession": "ENCFF100ATQ",
    "file_format": "fastq",
    "output_type": "reads",
    "biological_replicates": [1],
    "file_size_human": "3.1 GB",
    "status": "released"
  }
]

Integration

Related Skills

• pipeline-guide (parent): General pipeline selection and resource assessment
• accessibility-aggregation: Merge ATAC-seq peaks across samples
• quality-assessment: Deep-dive QC analysis beyond basic metrics
• regulatory-elements: Annotate peaks with regulatory element classifications
• compare-biosamples: Compare accessibility profiles across cell types
• pipeline-chipseq: Sibling pipeline for ChIP-seq data
• publication-trust: Verify literature claims backing analytical decisions

Presenting Results

When reporting ATAC-seq pipeline results:

• Mitochondrial fraction: Report it from qc/<sample>.idxstats.txt -- mapped reads on

the --mito_name row over the sum of the mapped column -- or from the samtools section of the MultiQC report (ideal <5%, acceptable <20%)

• Key QC metrics from the run: mapping rate (bowtie2 log, samtools flagstat),

duplication rate (Picard dup_metrics.txt), and read counts, all aggregated in qc/multiqc/multiqc_report.html

• Peak counts: Report the per-replicate MACS2 peak count and, for every replicate pair,

the IDR peak count at the 0.05 threshold. IDR is skipped when only one sample was processed

• TSS enrichment: State plainly that the workflow does not compute it. Report it only

if the user ran the manual step, with the quality tier (Excellent >=7, Good 5-7, Marginal 3-5, Poor <3)

• FRiP: Report the value from qc/<sample>.frip_mqc.tsv (>=0.3 for ATAC-seq). It is

computed from the blacklist-filtered BAM (all fragments) against the peaks called on the nucleosome-free fragments

• Fragment size distribution / NFR fraction / ataqv: also manual; do not report values

the run did not produce

• Output paths: List key outputs (peaks/narrow/, peaks/idr/, signal/,

filtered/nfr/, qc/, pipeline_info/)

• Next steps: Suggest motif-analysis for TF footprinting and de novo motif discovery,

or visualization-workflow for genome browser session generation

For the request: "$ARGUMENTS"