ATAC-seq Peak Calling
Basic MACS3 for ATAC-seq
# Standard ATAC-seq peak calling
macs3 callpeak \
-t sample.bam \
-f BAMPE \
-g hs \
-n sample \
--outdir peaks/ \
-q 0.05 \
--nomodel \
--shift -75 \
--extsize 150 \
--keep-dup all \
-B
Key ATAC-seq Parameters
# Explained parameters
macs3 callpeak \
-t sample.bam \ # Treatment BAM
-f BAMPE \ # Paired-end BAM (uses fragment size)
-g hs \ # Genome size: hs (human), mm (mouse)
-n sample \ # Output name prefix
--nomodel \ # Don't build shifting model
--shift -75 \ # Shift reads to center on Tn5 cut site
--extsize 150 \ # Extend reads to this size
--keep-dup all \ # Keep duplicates (ATAC has natural duplicates)
-B \ # Generate bedGraph for visualization
--call-summits # Call peak summits
Why These Parameters?
| Parameter |
Reason |
| --nomodel |
ATAC doesn't have control, can't build model |
| --shift -75 |
Centers on Tn5 insertion site |
| --extsize 150 |
Smooths signal around cut sites |
| --keep-dup all |
Tn5 creates duplicate cuts at accessible sites |
| -f BAMPE |
Uses actual fragment size from paired-end |
Paired-End vs Single-End
# Paired-end (recommended for ATAC)
macs3 callpeak -f BAMPE -t sample.bam ...
# Single-end (less common)
macs3 callpeak -f BAM -t sample.bam \
--nomodel --shift -75 --extsize 150 ...
Call Peaks on NFR Only
# First, filter to nucleosome-free reads (<100bp fragments)
samtools view -h sample.bam | \
awk 'substr($0,1,1)=="@" || ($9>0 && $9<100) || ($9<0 && $9>-100)' | \
samtools view -b > nfr.bam
# Call peaks on NFR
macs3 callpeak \
-t nfr.bam \
-f BAMPE \
-g hs \
-n sample_nfr \
--nomodel \
--shift -37 \
--extsize 75 \
--keep-dup all \
-q 0.01
Broad Peaks (Optional)
# For broader accessible regions
macs3 callpeak \
-t sample.bam \
-f BAMPE \
-g hs \
-n sample_broad \
--nomodel \
--shift -75 \
--extsize 150 \
--broad \
--broad-cutoff 0.1
Batch Processing
#!/bin/bash
GENOME=hs # hs for human, mm for mouse
OUTDIR=peaks
mkdir -p $OUTDIR
for bam in *.bam; do
sample=$(basename $bam .bam)
echo "Processing $sample..."
macs3 callpeak \
-t $bam \
-f BAMPE \
-g $GENOME \
-n $sample \
--outdir $OUTDIR \
--nomodel \
--shift -75 \
--extsize 150 \
--keep-dup all \
-q 0.05 \
-B \
--call-summits
done
Output Files
| File |
Description |
| _peaks.narrowPeak |
Peak locations (BED-like) |
| _summits.bed |
Peak summit positions |
| _peaks.xls |
Peak statistics (Excel format) |
| _treat_pileup.bdg |
Signal track (bedGraph) |
| _control_lambda.bdg |
Background (if control provided) |
narrowPeak Format
chr1 100 500 peak1 500 . 10.5 50.2 45.1 200
Columns: chrom, start, end, name, score, strand, signalValue, pValue, qValue, summit_offset
Convert to BigWig
# Sort bedGraph
sort -k1,1 -k2,2n sample_treat_pileup.bdg > sample.sorted.bdg
# Convert to BigWig
bedGraphToBigWig sample.sorted.bdg chrom.sizes sample.bw
Merge Replicates
# Pool BAMs before peak calling (recommended for final peaks)
samtools merge -@ 8 merged.bam rep1.bam rep2.bam rep3.bam
# Call peaks on merged
macs3 callpeak -t merged.bam -f BAMPE -g hs -n merged ...
IDR for Replicate Consistency
# Call peaks on each replicate
macs3 callpeak -t rep1.bam -f BAMPE -g hs -n rep1 ...
macs3 callpeak -t rep2.bam -f BAMPE -g hs -n rep2 ...
# Run IDR
idr --samples rep1_peaks.narrowPeak rep2_peaks.narrowPeak \
--input-file-type narrowPeak \
--output-file idr_peaks.txt \
--plot
# Filter by IDR threshold
awk '$5 >= 540' idr_peaks.txt > reproducible_peaks.bed
Related Skills
- read-alignment/bowtie2-alignment - Align ATAC-seq reads
- atac-seq/atac-qc - Quality control
- chip-seq/peak-calling - ChIP-seq comparison
- genome-intervals/bed-file-basics - Work with peak files
diegosouzapw