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

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.

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.