GPTomics/bio-workflows-atacseq-pipeline

Version Compatibility

Reference examples tested with: Bowtie2 2.5.3+, MACS3 3.0+, Genrich 0.6+, bedtools 2.31+, deepTools 3.5+ (alignmentSieve), fastp 0.23+, samtools 1.19+, DiffBind 3.12+, TOBIAS 0.16+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: macs3 callpeak -f BAMPE uses real fragment lengths and IGNORES --shift/--extsize/--nomodel; the cut-site style needs -f BAM/-f BED on Tn5-shifted reads. alignmentSieve --ATACshift applies the +4/-5 shift once. ENCODE ATAC-seq v3 and v4 QC thresholds are not interchangeable. Confirm in-tool before quoting.

ATAC-seq Pipeline

"Run ATAC-seq from FASTQ to differential accessibility and footprints" -> Chain QC/trim, alignment, chrM removal, dedup, a single Tn5 shift, peak calling, fixed-width consensus, differential accessibility, and footprinting.

• CLI + R: fastp -> bowtie2 -> drop chrM -> markdup -> Tn5 shift (once) -> macs3 -> Corces consensus -> DiffBind/csaw -> TOBIAS

This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.

The governing principle

ATAC-seq differs from ChIP-seq at four seams, and each is where the analysis goes wrong.

1. There is NO input control -- the shift-extend cut-site model IS the background. ATAC has no matched IP/input; peak significance comes from local lambda over the Tn5 insertion signal. Do not invent a "control"; commit instead to the build + ENCODE blacklist (removed before calling) as the coordinate frame.
2. The Tn5 +4/-5 shift is applied EXACTLY ONCE, after dedup and chrM removal. alignmentSieve --ATACshift (or one bedtools awk) applies it. Applying it twice, or combining -f BAMPE with --shift/--extsize (silently ignored), misplaces every cut site. Pick ONE calling mode: cut-site (-f BAM/-f BED + --nomodel --shift -75 --extsize 150) OR fragment (-f BAMPE on shifted reads, NO --shift).
3. chrM is removed BEFORE peak calling. Mitochondrial reads dominate ATAC libraries (often 20-50%, less with Omni-ATAC); leaving them in inflates depth and distorts FRiP and normalization.
4. Differential accessibility requires a FIXED-WIDTH consensus peakset. Variable-width MACS peaks make per-sample counts non-comparable. Build the Corces 501 bp iterative-overlap consensus (Corces 2018) so every region is the same width before counting; DiffBind/csaw then count into uniform intervals.

Reporting corollary: ENCODE ATAC v3 and v4 define TSS-enrichment/FRiP thresholds differently -- pick one standard and state which; do not mix rows across versions.

Pipeline map

FASTQ (paired, Nextera)
  | [1] QC & trim -----------------> fastp (Nextera adapters)   (read-qc/fastp-workflow)
  v
  | [2] Align ---------------------> bowtie2 --very-sensitive -X 2000   (read-alignment/bowtie2-alignment)
  v     ^-- commitment: build + ENCODE blacklist (NO input control)
  | [3] Drop chrM (BEFORE dedup/peaks) -> mito can be 20-50% of reads
  v
  | [4] Dedup --------------------> markdup -r   (alignment-files/duplicate-handling)
  v
  | [5] Tn5 shift ONCE ------------> alignmentSieve --ATACshift (+4/-5)
  v     ^-- pick ONE calling mode; never BAMPE + --shift
  | [6] Peak calling -------------> macs3 (cut-site -f BAM --shift/--extsize | -f BAMPE)  (atac-seq/atac-peak-calling)
  v
  | [7] Fixed-width consensus -----> Corces 501 bp iterative overlap   (atac-seq/consensus-peakset)
  v
  | [8] QC + differential + footprints -> TSS/FRiP/fragment; DiffBind/csaw; TOBIAS  (atac-seq/atac-qc, differential-accessibility, footprinting)
  v
Accessibility peaks + differential regions + TF activity

Made-once commitments

| Commitment | Choice | Consequence inherited downstream | |------------|--------|----------------------------------| | Build + blacklist | One build; ENCODE blacklist (removed before calling) | ATAC has no input, so the blacklist + shift-extend model ARE the background control | | Tn5 shift | Applied ONCE (--ATACshift), then ONE calling mode | Double-shift or BAMPE+--shift misplaces cut sites | | chrM handling | Removed before dedup/peaks | Mito reads (20-50%) inflate depth, FRiP, normalization | | Differential interval | Fixed-width Corces 501 bp consensus | Variable-width peaks make per-sample counts non-comparable |

The canonical order and why

1. QC/trim with Nextera adapters (CTGTCTCTTATACACATCT).
2. Align (bowtie2 --very-sensitive -X 2000) so the full nucleosome-spanning fragment distribution is captured.
3. Remove chrM, then compute NRF/PBC, then dedup -- order-trap on both ends: markdup -r physically removes duplicates, so NRF/PBC1 computed afterwards are identically 1.0; and mito reads are over-amplified, so computing them before chrM removal measures chrM chemistry, not nuclear-library complexity. The binding constraint is PRE-DEDUP. Mito must go before peak calling regardless.
4. Dedup (collate -> fixmate -m -> sort -> markdup -r).
5. Tn5 shift ONCE (alignmentSieve --ATACshift).
6. Call peaks in ONE mode -- order-trap: -f BAMPE + --shift/--extsize silently drops the flags.
7. Build the fixed-width consensus (Corces 501 bp) -- order-trap: differential on variable-width peaks is not comparable.
8. QC, differential (DiffBind/csaw on the consensus), footprinting (TOBIAS).

Choosing the caller and calling mode

Pipeline-level selection only; mechanism lives in the component skills.

| Fork | Lean toward | Hand off to | |------|-------------|-------------| | Caller | MACS3 (standard); Genrich (-j ATAC mode: handles replicates + chrM + blacklist in one pass); HMMRATAC (nucleosome-aware HMM) | atac-seq/atac-peak-calling | | Calling mode | Cut-site -f BAM/-f BED + --nomodel --shift -75 --extsize 150 (ENCODE smoothing window on shifted reads) vs fragment -f BAMPE on shifted reads (no --shift) | atac-seq/atac-peak-calling | | Consensus | Corces 2018 iterative-overlap fixed-width 501 bp | atac-seq/consensus-peakset | | Differential | DiffBind / csaw / DESeq2 on the fixed-width count matrix; spike-in for global shifts | atac-seq/differential-accessibility |

Primary path: Bowtie2 + Tn5 shift + MACS3

Goal: turn Nextera FASTQ into shifted, chrM-free peaks ready for a fixed-width consensus.

Approach: align with a wide insert window, drop chrM, dedup, Tn5-shift once, then call in ONE mode. Full runnable script: examples/atacseq_workflow.sh; differential: examples/differential_atac.R.

bowtie2 -p 8 -x bt2_index/genome -1 trimmed/${s}_R1.fq.gz -2 trimmed/${s}_R2.fq.gz \
    --very-sensitive --no-mixed --no-discordant -X 2000 2> aligned/${s}.log \
  | samtools view -@4 -bS -q 30 -f 2 - | samtools sort -@4 -o aligned/${s}.sorted.bam
samtools index aligned/${s}.sorted.bam

# Drop chrM BEFORE dedup/peaks (mito dominates ATAC), then dedup
samtools idxstats aligned/${s}.sorted.bam | cut -f1 | grep -v -e '^chrM$' -e '^MT$' \
  | xargs samtools view -b aligned/${s}.sorted.bam > aligned/${s}.noMT.bam
samtools collate -@8 -O -u aligned/${s}.noMT.bam | samtools fixmate -m -u - - \
  | samtools sort -@8 -u - | samtools markdup -r -@8 - aligned/${s}.dedup.bam
samtools index aligned/${s}.dedup.bam            # alignmentSieve needs an indexed input BAM

# Tn5 +4/-5 shift ONCE
alignmentSieve -b aligned/${s}.dedup.bam -o aligned/${s}.shifted.bam --ATACshift -p 8
samtools index aligned/${s}.shifted.bam

# Remove ENCODE blacklist regions BEFORE calling (the made-once commitment above; see the example script)
# Everything downstream (peaks, counts, footprints) consumes ${s}.filt.bam, never ${s}.shifted.bam.
# NOTE: examples/atacseq_workflow.sh names its blacklist-FILTERED output `.shifted.bam`; same reads,
# different name. Match on the step, not the suffix.
bedtools intersect -v -a aligned/${s}.shifted.bam -b "$BLACKLIST" > aligned/${s}.filt.bam
samtools index aligned/${s}.filt.bam

# Cut-site calling on the shifted, blacklist-filtered reads (ONE mode; do NOT also use -f BAMPE with these flags)
macs3 callpeak -t aligned/${s}.filt.bam -f BAM -g hs -n ${s} --outdir peaks \
    --nomodel --shift -75 --extsize 150 --keep-dup all -q 0.01

For the ENCODE 4 IDR + pseudoreplicate pipeline and the Corces 501 bp iterative-overlap consensus, see atac-seq/atac-peak-calling and atac-seq/consensus-peakset.

Differential accessibility and footprinting

Goal: compare accessibility across conditions on comparable intervals, then read TF activity.

Approach: count into the fixed-width consensus with DiffBind (or csaw), then run the TOBIAS three-step (ATACorrect -> ScoreBigwig -> BINDetect) for footprints.

library(DiffBind)                                  # counts into the fixed-width consensus
dba <- dba(sampleSheet = samples)                  # bamReads = shifted BAMs, Peaks = per-sample narrowPeak
dba <- dba.count(dba)                              # use summits/consensus for uniform width
dba <- dba.normalize(dba); dba <- dba.contrast(dba, categories = DBA_CONDITION)
dba <- dba.analyze(dba); report <- dba.report(dba)

# peaks/consensus.bed is the Corces 501 bp FIXED-WIDTH consensus from atac-seq/consensus-peakset (step 7).
# It is NOT peaks/consensus_peaks.narrowPeak, which is the variable-width pooled MACS3 call; build the
# fixed-width set first or these three commands have no input.
# TOBIAS three-step: bias-correct -> score -> detect bound motifs (differential across two conditions).
# Footprint on the BLACKLIST-FILTERED reads (${s}.filt.bam), the same reads MACS3 called peaks from --
# blacklist regions are artifact pileups, and bias-correcting over them corrupts the footprint scores.
TOBIAS ATACorrect -b aligned/${s}.filt.bam -g genome.fa -p peaks/consensus.bed --outdir foot --cores 8
TOBIAS ScoreBigwig --signal foot/${s}_corrected.bw --regions peaks/consensus.bed --output foot/${s}.bw --cores 8
TOBIAS BINDetect --motifs motifs.jaspar --signals foot/ctrl.bw foot/treat.bw --genome genome.fa \
    --peaks peaks/consensus.bed --outdir foot/bindetect --cores 8

QC checkpoints between steps

| After | Gate | Interpretation | |-------|------|----------------| | Alignment | Mapping >80%, mito <20% (Omni-ATAC lower) | High mito = suboptimal lysis; drop before calling | | PRE-dedup | NRF >0.8, PBC1 >0.8 | Low complexity = over-amplification/low input; compute before dedup | | Peaks | FRiP >0.2, TSS enrichment >5 (v3) | Low TSS/FRiP = over/under-digestion or degraded chromatin (atac-seq/atac-qc) | | Fragment size | NFR <100 bp, mono ~200 bp, di ~400 bp periodicity | Loss of nucleosome periodicity = over-digestion (Tn5:DNA too high) | | Consensus | Fixed-width (501 bp) built before counting | Variable-width peaks make counts non-comparable |

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | Depth/FRiP dominated by one contig; few real peaks | chrM not removed before calling | Drop chrM/MT before dedup and peak calling | | Cut sites offset / footprints smeared | Tn5 shift applied twice, or -f BAMPE used with --shift/--extsize | Shift ONCE; pick ONE calling mode (cut-site -f BAM OR fragment -f BAMPE) | | Differential counts not comparable across samples | Counted into variable-width MACS peaks | Build the Corces 501 bp fixed-width consensus first | | Looked for an input/IgG track and found none | ATAC has no input control | Use the blacklist + shift-extend model as background; do not fabricate a control | | QC numbers disagree with a reference | Mixed ENCODE v3 and v4 thresholds | Pick one ENCODE version and report which |

Pipeline map (hand-offs)

• read-qc/fastp-workflow - Nextera adapter trimming
• read-alignment/bowtie2-alignment - the aligner, wide insert window
• alignment-files/duplicate-handling - collate/fixmate/sort/markdup order
• atac-seq/atac-peak-calling - MACS3/Genrich/HMMRATAC, ENCODE 4 IDR, calling modes
• atac-seq/atac-qc - TSS enrichment, FRiP, NRF/PBC, fragment periodicity
• atac-seq/consensus-peakset - Corces 2018 iterative-overlap fixed-width consensus
• atac-seq/differential-accessibility - DiffBind/csaw/DESeq2 on the consensus
• atac-seq/footprinting - TOBIAS three-step and per-TF failure modes
• atac-seq/nucleosome-positioning - V-plot, NucleoATAC, +1 nucleosome

The complete runnable scripts are in this skill's examples/ (atacseq_workflow.sh, differential_atac.R).

Related Skills

• database-access/sra-data - Pull ATAC-seq FASTQ from SRA / ENA
• database-access/geo-data - Resolve GEO accessions for ATAC datasets
• read-qc/fastp-workflow - Nextera adapter trimming and quality filtering
• read-alignment/bowtie2-alignment - Standard ATAC-seq aligner
• alignment-files/duplicate-handling - MarkDuplicates pre-peak-calling
• atac-seq/atac-peak-calling - MACS3 / Genrich / HMMRATAC details, ENCODE 4 IDR
• atac-seq/atac-qc - TSS enrichment, FRiP, NRF/PBC1/PBC2 details
• atac-seq/consensus-peakset - Corces 2018 iterative-overlap fixed-width consensus
• atac-seq/differential-accessibility - DiffBind / csaw / DESeq2; spike-in normalization
• atac-seq/footprinting - TOBIAS three-step; per-TF failure modes
• atac-seq/motif-deviation - chromVAR for motif accessibility variability
• atac-seq/nucleosome-positioning - V-plot, NucleoATAC, +1 nucleosome
• atac-seq/single-cell-atac - For scATAC instead of bulk
• atac-seq/co-accessibility - Cicero cis-regulatory inference
• atac-seq/enhancer-gene-linking - ABC, ENCODE-rE2G enhancer-gene mapping
• atac-seq/deep-learning-atac - chromBPNet variant-effect prediction
• atac-seq/allele-specific-accessibility - WASP + caQTL mapping
• chip-seq/peak-annotation - Annotate ATAC peaks to genes

References

• Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ (2013) Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nature Methods 10:1213-1218. DOI 10.1038/nmeth.2688. (original ATAC-seq.)
• Corces MR, Trevino AE, Hamilton EG, et al (2017) An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nature Methods 14:959-962. DOI 10.1038/nmeth.4396. (Omni-ATAC.)
• Corces MR, Granja JM, Shams S, et al (2018) The chromatin accessibility landscape of primary human cancers. Science 362:eaav1898. DOI 10.1126/science.aav1898. (fixed-width iterative-overlap consensus peakset.)
• Bentsen M, Goymann P, Schultheis H, et al (2020) ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. Nature Communications 11:4267. DOI 10.1038/s41467-020-18035-1. (TOBIAS.)