Bio Atac Seq Atac Peak Calling

Version Compatibility

Reference examples tested with: MACS3 3.0.2+, MACS2 2.2.9+, Genrich 0.6.1+, HMMRATAC 1.2+ (now bundled in MACS3 as macs3 hmmratac), samtools 1.19+, bedtools 2.31+, IDR 2.0.4+.

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

• CLI: <tool> --version then <tool> --help to confirm flags

If code throws unexpected errors, introspect the installed binary (<tool> -h) and adapt the example to match the actual CLI rather than retrying.

ATAC-seq Peak Calling

"Call accessible regions from my ATAC-seq BAM" -> Identify Tn5-hypersensitive open chromatin, treating fragments as point insertion events (not protein-bound regions as in ChIP-seq) and accounting for the lack of input control.

• CLI (canonical, ENCODE 4): macs2 callpeak -t atac.bam -f BAMPE -g hs -n sample --nomodel --shift -75 --extsize 150 --keep-dup all -B --SPMR -p 0.01
• CLI (HMM-based, single sample): macs3 hmmratac -i atac.bam -n sample --outdir hmm_out
• CLI (joint replicates): Genrich -j -t rep1.bam,rep2.bam -o peaks.narrowPeak -e chrM -E blacklist.bed

The -p 0.01 (loose) plus IDR is the ENCODE pattern: low stringency increases peak overlap between replicates, and IDR rescues the reproducible set. Single-sample workflows usually swap to -q 0.05 instead.

Algorithmic Taxonomy

Methodology evolves; verify the current ENCODE ATAC-seq Standards (encodeproject.org pipelines/atac-seq) before locking parameters. ENCODE 4 still defaults to MACS2 (not MACS3) at time of writing; macs3 callpeak is API-compatible for ATAC parameters but not yet the official ENCODE binary.

Shift-Extend vs BAMPE: The Critical Choice

Two valid ways to feed paired-end ATAC into MACS:

Pattern A (ENCODE / "single-end-ified"): -f BAMPE actually IGNORES --shift/--extsize. To activate them, use -f BAM and treat each end independently. ENCODE's pipeline uses -f BAM --shift -75 --extsize 150 to model each Tn5 cut as a 150 bp window centered on the insertion site, ignoring fragment lengths.

Pattern B (paired-fragment): -f BAMPE uses the full paired-end fragment span as the signal interval. Best when fragment lengths are biologically meaningful (e.g., NFR-only peak calling at 38/75 bp). In BAMPE mode, do NOT set --shift/--extsize (silently ignored, but confusing).

For most bulk ATAC, Pattern A matches ENCODE convention and is reproducible against published peak sets. Pattern B can be more sensitive at narrow regulatory elements but does not match ENCODE outputs.

Effective Genome Size

-g hs and -g mm are MACS shorthands for old defaults. Modern values:

Wrong size shifts every q-value but rarely changes peak ranks. Use unique-kmers.py (khmer) or the deepTools tabulated values for exact sizes; the shorthand is a decade-old approximation.

Effective Genome Size: When It Matters

Trigger: Comparing peaks across genome builds or species; reproducing published q-value cutoffs; hi-resolution lambda estimation.

Mechanism: MACS estimates genome-wide lambda as total_reads / effective_size. Wrong size -> wrong null -> shifted q-values, especially at the marginal cutoff.

Symptom: Peak counts diverge ~10-20% from published numbers when re-running an old dataset.

Fix: Pull the read-length-matched value from deepTools effectiveGenomeSize table. For pipelines, parameterize this; never inline the shorthand for cross-study comparisons.

Per-Tool Failure Modes

MACS2/MACS3 -- Confounded NFR + broad accessibility

Trigger: Cell type with extended open domains (e.g., active super-enhancers, MYOD1 regulons, locus-control regions).

Mechanism: Default narrow-peak mode segments wide accessible domains into multiple smaller peaks at local lambda spikes; --broad --broad-cutoff 0.1 merges them but inflates total length and breaks IDR comparability.

Symptom: Peak count >> 200k for human bulk ATAC at ENCODE depth; mean peak width < 200 bp; visual inspection in IGV shows 3-5 calls under one continuous accessibility block.

Fix: Run both narrow and broad; use narrow for differential analysis, broad for domain-level enrichment (e.g., super-enhancer overlap). Do NOT use --call-summits for broad mode.

Genrich -- Replicate weighting and chrM exclusion

Trigger: Replicates with very different library sizes; high-mitochondrial samples not pre-filtered.

Mechanism: Genrich's joint mode pools reads via Fisher's method. Library-size imbalance dominates the joint p-value; chrM reads inflate background unless -e chrM is set.

Symptom: Most-significant peaks cluster on chrM or on the largest-library replicate's high-coverage regions.

Fix: Always pass -e chrM (Genrich 0.6+) and -E blacklist.bed. Down-sample BAMs to common depth (samtools view -s) before joint calling if libraries differ >2x. Add -r to remove PCR duplicates inside Genrich, OR pre-deduplicate (do not do both).

MACS3 hmmratac (HMMRATAC) -- Depth and fragment-size dependence

Trigger: Library < 25M nuclear reads, or libraries with degraded chromatin and flat fragment-size distribution.

Mechanism: The 3-state HMM is trained from fragment-size classes (NFR ~50 bp, mono ~200 bp, di ~400 bp peaks). Without periodicity the emission distributions collapse and the HMM cannot separate states.

Symptom: Output BED is empty, or all peaks are tiny (~150 bp) with no nucleosome flanks called; runtime explodes (>24h) on shallow data.

Fix: Verify fragment-size periodicity in QC first (atac-qc skill). If flat, fall back to MACS3 callpeak. HMMRATAC needs >= 30M deduplicated nuclear reads per ENCODE recommendation.

HOMER findPeaks -- Window-size sensitivity

Trigger: Default -style dnase uses 75 bp peaks; ATAC peaks are 250-500 bp typically.

Mechanism: HOMER's window-based caller does not auto-fit width to ATAC.

Fix: Use -style factor -size 150 for narrow ATAC peaks, or skip HOMER for peak calling and use it only for downstream motif analysis on MACS peaks.

Aligner choice -- chromap vs bwa-mem2 vs bowtie2 affects peak shape

Trigger: Switching aligners between datasets and expecting reproducible peaks.

Mechanism: chromap (Zhang 2021) applies its own ATAC-specific 4 bp / -5 bp Tn5 shift before fragment output; bwa-mem2 and bowtie2 do not. Downstream --shift -75 --extsize 150 parameters are calibrated for unshifted bwa/bowtie BAMs; applying them to chromap output double-shifts the signal.

Symptom: Peaks called from chromap output are shifted by ~5-10 bp relative to bwa output at the same locus.

Fix: When using chromap, drop --shift and --extsize (chromap's pre-shift is sufficient) OR use chromap's --no-correction flag to disable Tn5 shift and proceed with standard MACS parameters. Document the aligner version and any shift choices in methods. Within a project, pin the aligner.

Single-sample (no replicate) -- Rotation / circular-shift permutation

Trigger: Single biological sample without any replicate for IDR.

Mechanism: IDR requires two replicates by construction. For n=1, statistical confidence per peak comes from local background (Poisson p-value) but reproducibility cannot be assessed.

Fix: Apply a stricter -q 0.01 (vs ENCODE -p 0.01 + IDR pattern) and additionally apply rotation/circular-shift permutation: shift the BAM cuts by a random distance modulo each chromosome and re-call peaks; the per-peak persistence rate across rotations is a non-parametric reproducibility proxy. Document this is a single-sample heuristic, not ENCODE-compliant.

ENCODE 3 vs ENCODE 4 Differences

To reproduce a published ENCODE 3 dataset, pin the original pipeline and threshold exactly. ENCODE 4 results are not directly numerically comparable to ENCODE 3 even on the same input BAM.

Super-Enhancer Detection

For active super-enhancer (SE) annotation alongside narrow-peak workflow, ROSE (Whyte 2013) and LILY (Boeva 2017) stitch ATAC or H3K27ac peaks separated by < 12.5 kb and rank by signal:

# ROSE expects H3K27ac BAM but works on ATAC narrowPeak with care
ROSE_main.py -g hg38 -i atac_peaks.gff -r atac.bam -o rose_out/ -t 2500

ROSE-style stitching is complementary to MACS3 narrow peaks: narrow peaks for differential analysis; SE annotation for biology interpretation. SE calls require H3K27ac input for definitive annotation; ATAC alone produces "stretch enhancers" that overlap but are not identical to H3K27ac SE.

ENCODE 4 ATAC-seq Pipeline (Reference Implementation)

The exact ENCODE pattern produces the most-comparable peak sets:

# Per-replicate peak calling (loose threshold)
macs2 callpeak \
    -t rep1.filt.dedup.bam \
    -f BAM -g hs \
    -n rep1 --outdir peaks/rep1/ \
    --nomodel --shift -75 --extsize 150 \
    --keep-dup all \
    -B --SPMR \
    -p 0.01

# Pooled (all replicates)
macs2 callpeak \
    -t rep1.filt.dedup.bam rep2.filt.dedup.bam \
    -f BAM -g hs -n pooled --outdir peaks/pooled/ \
    --nomodel --shift -75 --extsize 150 --keep-dup all -B --SPMR -p 0.01

# Pseudoreplicates (split each rep BAM in half)
samtools view -b -h -s 1.5 rep1.filt.dedup.bam > rep1.psr1.bam   # seed.fraction
samtools view -b -h -s 2.5 rep1.filt.dedup.bam > rep1.psr2.bam   # different seed
# (call peaks on each pseudoreplicate the same way)

--SPMR writes signal as Signal Per Million Reads (normalized bedGraph). -p 0.01 is intentionally loose; IDR will tighten to a reproducible set.

IDR for Reproducible Peaks

Goal: Find peaks reproducible across biological replicates at controlled IDR.

Approach: Score paired peak lists by signalValue, fit IDR's two-component mixture (reproducible + noise), threshold at IDR <= 0.05 (true reps) or 0.10 (pseudoreplicates).

# Sort peaks by p-value (column 8) so IDR scores by significance
sort -k8,8nr rep1_peaks.narrowPeak > rep1.sorted.narrowPeak
sort -k8,8nr rep2_peaks.narrowPeak > rep2.sorted.narrowPeak

# True replicates -- threshold IDR <= 0.05
idr --samples rep1.sorted.narrowPeak rep2.sorted.narrowPeak \
    --input-file-type narrowPeak --rank p.value \
    --output-file true_reps.idr \
    --idr-threshold 0.05 --plot --log-output-file idr.log

# Pseudoreplicates -- threshold IDR <= 0.10 (looser, ENCODE Nself <= 2 rule)
idr --samples psr1_peaks.narrowPeak psr2_peaks.narrowPeak \
    --input-file-type narrowPeak --rank p.value \
    --output-file psr.idr --idr-threshold 0.10 --plot

ENCODE consistency rules: Nt = peaks passing IDR on true reps; Nself = peaks passing IDR on pseudoreps. Library passes if max(Nt, Nself) / min(Nt, Nself) <= 2. If both ratios > 2, the library is rejected.

IDR fails when: Ranking column choice matters. --rank p.value (column 8) is robust; --rank signal.value (column 7) breaks if MACS pile-up scaling differs between replicates.

Decision Tree by Experimental Scenario

Reconciliation: When Callers Disagree

Operational rule for high-confidence reporting: Require a peak to pass IDR <= 0.05 on true replicates AND survive blacklist/greylist filtering AND have mean signalValue >= 5 across reps. Two callers from different families (MACS + Genrich) agreeing within 250 bp is acceptable evidence when IDR is unavailable.

Blacklist and Greylist

# ENCODE blacklist (Amemiya 2019) -- always remove
wget https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz
gunzip hg38-blacklist.v2.bed.gz
bedtools intersect -v -a peaks.narrowPeak -b hg38-blacklist.v2.bed > peaks.no_blacklist.narrowPeak

# Sample-specific greylist (input-derived high-signal regions; rarely available for ATAC)
# For ATAC, ENCODE recommends pooling all samples' top-percentile signal and removing
# regions exceeding 100x median coverage as a "soft greylist"

Blacklist is mandatory; greylist is optional and most useful when the same library prep produces consistent artifact regions across samples.

NFR-Only Peak Calling

Goal: Call peaks using only sub-nucleosomal fragments (<100 bp) for sharper TF-binding-relevant accessibility.

Approach: Pre-filter BAM to short fragments, then call peaks with parameters scaled to the smaller fragment length.

samtools view -h sample.dedup.bam | \
    awk 'substr($0,1,1)=="@" || ($9 > 0 && $9 < 100) || ($9 < 0 && $9 > -100)' | \
    samtools view -b > nfr.bam
samtools index nfr.bam

macs2 callpeak -t nfr.bam -f BAM -g hs -n sample_nfr \
    --nomodel --shift -37 --extsize 75 \
    --keep-dup all -p 0.01

--shift -37 --extsize 75 halves both parameters to match shorter fragments; this is what TOBIAS recommends for footprinting input.

Output Files (narrowPeak)

Convert to bigWig for browsers: sort -k1,1 -k2,2n sample_treat_pileup.bdg | bedGraphToBigWig - chrom.sizes sample.bw.

Common Errors

References

• Buenrostro JD et al 2013 Nat Methods 10:1213 (ATAC-seq protocol)
• Corces MR et al 2017 Nat Methods 14:959 (Omni-ATAC, fixed-width peaks)
• Tarbell ED & Liu T 2019 Nucleic Acids Res 47:e91 (HMMRATAC)
• Gaspar JM 2018 bioRxiv 459545 (Genrich)
• Li Q et al 2011 Ann Appl Stat 5:1752 (IDR framework)
• Landt SG et al 2012 Genome Res 22:1813 (ENCODE/modENCODE peak calling guidelines, IDR Nself rule)
• Amemiya HM et al 2019 Sci Rep 9:9354 (ENCODE blacklist v2)
• ENCODE ATAC-seq Standards (encodeproject.org/atac-seq) -- canonical pipeline parameters

Related Skills

• atac-seq/atac-qc - Verify TSS enrichment, FRiP, and fragment periodicity before calling
• atac-seq/consensus-peakset - Combine per-sample peaks into a fixed-width differential-ready set
• atac-seq/single-cell-atac - Pseudobulk peak calling per cluster
• atac-seq/differential-accessibility - Downstream DiffBind/csaw/DESeq2 testing
• atac-seq/deep-learning-atac - chromBPNet bias-corrected per-base profiles as alternative input
• read-alignment/bowtie2-alignment - Upstream ATAC alignment
• alignment-files/duplicate-handling - Pre-call dedup with Picard MarkDuplicates
• chip-seq/peak-calling - ChIP-seq comparison (uses input control)
• chip-seq/super-enhancers - ROSE / LILY for super-enhancer annotation
• genome-intervals/bed-file-basics - Peak file manipulation