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 BAM -g hs -n sample --nomodel --shift -75 --extsize 150 --keep-dup all -B --SPMR -p 0.01 (use -f BAM, not -f BAMPE -- BAMPE reads true fragment ends and ignores --shift/--extsize)
• 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 combines the per-replicate p-values at each position 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 deep coverage (a practical minimum around 30M deduplicated nuclear reads).

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 omit chromap's --Tn5-shift (the shift is opt-in, applied only when that flag or an ATAC preset is set) so it is not double-applied, then 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.