Peak Calling with MACS3: How ChIP-seq and ATAC-seq Peaks Are Found

A ChIP-seq or ATAC-seq experiment leaves you with the same unglamorous artifact as almost every other assay: a sorted, indexed BAM full of aligned reads. Nowhere does it say "a transcription factor bound here" or "this promoter was open." Peak calling is the step that reads the coverage landscape and decides which pileups are real enrichment and which are lumpy background from open chromatin, copy-number variation, and sequencing bias. MACS3 is the tool most people reach for, and understanding the statistics it runs — fragment size estimation, a local background model, and a Poisson enrichment test — is the difference between trusting your peak list and being fooled by it.

What a Peak Actually Is

A peak is a genomic interval where more reads pile up than you would expect by chance. That is the whole idea, but "by chance" is doing enormous work in that sentence. If you scattered the same number of reads uniformly across the genome, every position would have some expected coverage. Real signal shows up as regions where the observed coverage climbs far above that expectation.

The problem is that the genome is not a uniform place. Open chromatin, repetitive regions, mappability differences, and local amplification all create coverage bumps that have nothing to do with the protein or accessibility you care about. So a naive "count reads in a window and flag the tall ones" approach drowns in false positives. MACS3's core contribution is a local background model that adapts the null expectation to each region, so a genuine binding site stands out even when it sits on top of an already-busy stretch of chromatin.

Fragment Size Estimation and Read Shifting

Here is a subtlety that trips up beginners. In a single-end ChIP-seq library, you sequence one end of each DNA fragment. Reads from the plus strand pile up on the left side of a binding site and reads from the minus strand pile up on the right, because sequencing starts from opposite ends of fragments that straddle the protein. The result is two coverage humps flanking the true site, not one peak on top of it.

MACS3 measures the distance between those plus- and minus-strand humps to estimate the fragment size d, then shifts each read d/2 toward its 3' end so the two humps collapse into one sharp peak centered on the real site. You will see this reported in the log as the estimated fragment length. If MACS3 cannot find a clean pattern it falls back to a default, and you can override the whole process with --nomodel and --extsize, which is standard practice for ATAC-seq.

The Local Lambda Background Model

Once reads are extended to full fragments, MACS3 builds its null model. It assumes that, in the absence of real enrichment, the number of reads falling in a window follows a Poisson distribution with rate parameter lambda. The trick is that lambda is not a single genome-wide constant. MACS3 computes a local lambda for every candidate region by taking the maximum expected read count across several nested windows — typically the peak region itself, a 1 kb neighborhood, a 10 kb neighborhood, and the genome-wide average.

Taking the maximum is deliberately conservative: if a region sits on a local amplification or mappability artifact, the wider windows capture that elevated baseline and the bar a real peak must clear rises with it. This is why a control matters so much — with an input or IgG control, MACS3 uses control coverage to inform the local lambda directly, catching biases shared by both channels. Without one, it leans entirely on the treatment sample's own local windows, which works but is weaker.

The final local lambda is the maximum of whichever of these apply, so no single narrow artifact can sneak a peak through.

The Poisson Enrichment Test

With a local lambda in hand, the test itself is straightforward. For a candidate region with observed read count k, MACS3 asks: under a Poisson distribution with rate lambda, how surprising is it to see k or more reads? That upper-tail probability is the p-value. Regions where observing this much coverage would be wildly improbable under the background model become peaks.

Because MACS3 evaluates a huge number of positions across the genome, raw p-values are meaningless on their own — test enough windows and some will look enriched by luck alone. MACS3 applies a Benjamini-Hochberg correction to convert p-values into q-values (false discovery rate), and the q-value is what you actually threshold on. The default is -q 0.05, meaning you accept a peak list where roughly 5% of calls are expected to be false. It is the same multiple-testing logic that shows up anywhere you score thousands of candidates at once and need to control the flood of false positives.

Narrow vs Broad: Choosing the Right Mode

Not all enrichment looks the same, and MACS3 has two shapes in mind. Transcription factors and most ATAC-seq accessibility produce compact, punctate peaks — a few hundred base pairs of sharp signal. Some histone marks, like H3K27me3 or H3K36me3, spread across tens of kilobases as diffuse domains that a narrow peak caller chops into a confetti of tiny fragments.

For punctate signal, the default macs3 callpeak in narrow mode is correct. For diffuse marks, add --broad, which tells MACS3 to first find enriched regions at the stringent q-value, then stitch nearby weakly-enriched segments together using a second, looser cutoff set by --broad-cutoff. The output changes too: narrow mode writes a .narrowPeak file with a summit column, while broad mode writes a .broadPeak (or .gappedPeak) without a single summit, because a spread-out domain does not have one.

Picking the wrong mode is one of the most common mistakes: run a broad histone mark in narrow mode and your peak count will look impressive while being biologically wrong.

A Realistic Command Walkthrough

Before MACS3 ever runs, you prepare inputs with the standard toolkit. Align with BWA or Bowtie2 (see how read aligners work), then filter and sort with samtools so only confident, primary alignments remain:

# Keep primary, properly aligned reads with decent mapping quality
samtools view -b -q 30 -F 1804 aln.bam > filtered.bam
samtools sort -o filtered.sorted.bam filtered.bam
samtools index filtered.sorted.bam

The -q 30 keeps reads with MAPQ >= 30, and -F 1804 drops unmapped reads, secondary and supplementary alignments, duplicates, and reads failing QC. Clean input matters: duplicates and multi-mappers inflate false peaks. Then call peaks:

# Narrow peaks for a transcription factor, with an input control
macs3 callpeak \
  -t chip.sorted.bam \
  -c input.sorted.bam \
  -f BAM -g hs \
  -n mytf --outdir peaks/ \
  -q 0.05

Here -t is treatment, -c is control, -g hs sets the effective human genome size, and -n names the output. For a broad histone mark you would add --broad --broad-cutoff 0.1. MACS3 writes a .narrowPeak BED-like file, a summits file, and a bedGraph if you pass --bdg.

Reading the Output and Sanity-Checking in IGV

The .narrowPeak file is a BED6+4 format: chromosome, start, end, name, score, strand, then signal value, -log10(p), -log10(q), and the summit offset. The last three columns are what you sort and filter on. A quick pass with bedtools tells you how many peaks you called, how they distribute across chromosomes, and whether they overlap features like promoters.

But numbers lie, and the fastest reality check is visual. Load the treatment and control bedGraph or BAM tracks into IGV alongside your .narrowPeak file and look at a few called peaks. A real transcription factor peak shows a clean, focused pileup in the treatment track that is absent or flat in the control. If your "peaks" look like gentle genome-wide swells present in both channels, your q-value is too loose or your background model is being fooled.

Practical Takeaways

• MACS3 estimates fragment size d from strand-shifted read humps, then extends and shifts reads so the two humps collapse into one centered peak — override with --nomodel --extsize for ATAC-seq.
• The local lambda model computes a per-region Poisson background from nested windows (1 kb, 10 kb, genome-wide), taking the maximum so local artifacts cannot slip peaks through.
• Peaks come from an upper-tail Poisson p-value, then a Benjamini-Hochberg q-value; threshold on q with -q 0.05, not on raw p-values.
• Use narrow mode for transcription factors and sharp marks, --broad with --broad-cutoff for diffuse histone domains like H3K27me3.
• Filter inputs first with samtools (MAPQ >= 30, drop duplicates and multi-mappers), and always subtract the ENCODE blacklist with bedtools before downstream analysis.
• Sanity-check a handful of peaks in IGV: real signal is focused in treatment and flat in control.

Peak calling is where read alignment turns into biology, so it rewards care at both ends — clean, deduplicated input and a critical eye on the output. If you have not set up the surrounding toolkit yet, start with the samtools and bcftools install guide and the bedtools install guide; with those in place, MACS3 slots neatly into a reproducible epigenomics pipeline.