name: bio-workflows-atacseq-pipeline description: Orchestrates the end-to-end bulk ATAC-seq pipeline from FASTQ to differential accessibility and TF footprints, chaining Nextera-aware fastp QC, Bowtie2 alignment, chrM removal, dedup, a single Tn5 +4/-5 shift, MACS3 peak calling, Corces fixed-width consensus, DiffBind/csaw differential accessibility, and TOBIAS footprinting. Use when committing the reference build + blacklist once, recognizing ATAC has NO input control (the shift-extend model IS the background), applying the Tn5 shift exactly once (never combining -f BAMPE with --shift/--extsize), removing chrM before calling, building a fixed-width consensus so per-sample counts are comparable, or choosing MACS3 vs Genrich vs HMMRATAC. Hands mechanism to the atac-seq component skills; not a re-teach of any single step. workflow: true depends_on: - read-qc/fastp-workflow - read-alignment/bowtie2-alignment - alignment-files/duplicate-handling - atac-seq/atac-peak-calling - atac-seq/atac-qc - atac-seq/consensus-peakset - atac-seq/differential-accessibility - atac-seq/footprinting - atac-seq/motif-deviation - atac-seq/nucleosome-positioning qc_checkpoints: - after_qc: "Q30 >85%, adapter content <5% (Nextera)" - after_alignment: "Mapping rate >80%, mitochondrial <20% (Omni-ATAC lower)" - before_dedup: "NRF >0.8, PBC1 >0.8 (computed PRE-dedup)" - after_peaks: "FRiP >0.2, TSS enrichment >5 (ENCODE v3; v4 thresholds differ, do not mix)" - after_consensus: "Fixed-width (Corces 501 bp) consensus built before counting for differential" origin: openai4s category: bioskills/workflows metadata: tool_type: mixed primary_tool: MACS3 third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT
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.
"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.
ATAC-seq differs from ChIP-seq at four seams, and each is where the analysis goes wrong.
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).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.
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
| 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 |
CTGTCTCTTATACACATCT).--very-sensitive -X 2000) so the full nucleosome-spanning fragment distribution is captured.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.alignmentSieve --ATACshift).-f BAMPE + --shift/--extsize silently drops the flags.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 |
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.
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
| 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 |
| 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 |
The complete runnable scripts are in this skill's examples/ (atacseq_workflow.sh, differential_atac.R).