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.

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

QC/trim with Nextera adapters (CTGTCTCTTATACACATCT).
Align (bowtie2 --very-sensitive -X 2000) so the full nucleosome-spanning fragment distribution is captured.
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.
Dedup (collate -> fixmate -m -> sort -> markdup -r).
Tn5 shift ONCE (alignmentSieve --ATACshift).
Call peaks in ONE mode -- order-trap: -f BAMPE + --shift/--extsize silently drops the flags.
Build the fixed-width consensus (Corces 501 bp) -- order-trap: differential on variable-width peaks is not comparable.
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.)