Chatomics Field GuideWhat They Don't Teach You

Glossary · Genomics and Variants

Peak calling

The MACS2 output you trust without checking known biology is the fastest way to publish noise as a binding site.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 3 min read

Also: MACS2, MACS3, narrowPeak

Definition

Peak calling is the step that converts aligned reads from a ChIP-seq, ATAC-seq, CUT&RUN, or CUT&TAG experiment into a discrete list of genomic intervals where the signal is significantly enriched above local background. MACS2 (and its newer rewrite, MACS3) is the field's default tool: it models the expected bimodal read pileup around a true binding site, estimates background with a dynamic Poisson distribution that adjusts to local genome complexity, and reports each enriched region with a p-value and a Benjamini-Hochberg-corrected q-value. The standard output format is narrowPeak, a 10-column BED6+4 file with one row per peak: coordinates, signal strength, -log10(p-value), -log10(q-value), and an optional summit offset.

You hit peak calling right after alignment and duplicate marking, before you do anything that treats a genomic interval as biologically meaningful: motif enrichment, annotation to the nearest gene, differential binding between conditions, or overlap with a public ChIP-seq track. Whatever MACS2 emits becomes your working definition of "binding site" or "accessible region" for every step downstream, so a wrong parameter here doesn't stay local, it propagates through the whole analysis.

The decision that actually matters isn't which tool to run (most labs default to MACS2 without debate); it's whether the parameters match the assay. ChIP-seq, ATAC-seq, and CUT&RUN produce fragments with different size distributions and different strand-shift signatures, and MACS2's peak model assumes you told it which one you're looking at.

Why it matters

Get the parameters wrong and you don't get an error message, you get a peak file that looks plausible and is wrong. The clearest check is biological: if you run ER (estrogen receptor) ChIP-seq, you should see peaks at TFF1 and GREB1. If you don't, the problem is the ChIP itself (failed pulldown), the antibody (nonspecific), or the sample (no real ER binding), and no amount of tuning MACS2's q-value cutoff fixes any of those three. Peak calling can't tell you which failure mode you're in; matching your peak list against loci you already know the biology of can.

The same logic applies to parameter choice. Feed ATAC-seq reads to MACS2 with ChIP-seq defaults, letting it build an automatic shift model instead of setting --nomodel, and you corrupt the fragment-length estimate, because Tn5 insertion produces a different fragment-size distribution than sonicated ChIP-seq chromatin. The ENCODE3 pipeline fixes this with --nomodel --shift -37 --extsize 73, tuned to the nucleosome-free footprint Tn5 cuts around. Skip that and your peaks shift off the true cut sites, merging into broad blobs or splitting a single accessible region into several fragments.

Where people get it wrong

The mistake is treating peak calling as assay-agnostic: running the same MACS2 command on ChIP-seq and ATAC-seq data because "it's the same tool." MACS2's automatic shift-model building exists to detect the bimodal strand-shift pattern of ChIP-seq fragments (forward-strand reads pile up upstream of the binding site, reverse-strand reads pile up downstream, and MACS uses the offset between them to estimate fragment length). ATAC-seq fragments don't follow that pattern; Tn5 transposase inserts directly into accessible chromatin, so the automatic shift model either fails to converge or converges on a meaningless number, and every downstream peak inherits that error. The fix isn't a different tool, it's --nomodel plus assay-appropriate --shift/--extsize values, set explicitly instead of left to MACS2's defaults.

A concrete example

Calling peaks from an ER ChIP-seq BAM against an input control with pre-estimated fragment length, then a separate ATAC-seq call using ENCODE3's nucleosome-adjusted parameters. The last step is the one people skip: load the narrowPeak file in IGV and confirm a peak actually sits at TFF1 or GREB1. If neither shows up, the q-value cutoff isn't the problem, the antibody or the ChIP itself is.

bash
# ChIP-seq: sharp, narrow peaks (transcription factor), fragment length pre-estimated
macs2 callpeak -t IP.bam -c Input.bam -n ER_rep1 \
  -p 0.01 --nomodel --extsize 146 --keep-dup all -g hs

# ATAC-seq: nucleosome-adjusted parameters (ENCODE3), never the ChIP-seq shift model
macs2 callpeak -t ATAC.bam -c Input.bam -n atac_rep1 \
  --nomodel --shift -37 --extsize 73 -q 0.05 -g hs

# sanity check: does the narrowPeak file hit the loci you expect?
grep -E "TFF1|GREB1" ER_rep1_peaks.narrowPeak

Related terms

Questions people ask

What is peak calling in ChIP-seq and ATAC-seq analysis?

It's the step that identifies genomic regions with significantly enriched read coverage after alignment, representing protein-DNA binding sites in ChIP-seq or accessible chromatin in ATAC-seq. MACS2 is the standard tool; it outputs a narrowPeak file listing each region's coordinates, signal strength, and statistical significance.

What's the difference between MACS2 narrow peaks and broad peaks?

Narrow peaks use a Poisson p-value threshold (commonly 0.01) and suit sharp, punctate signal like transcription-factor ChIP-seq. Broad peaks (--broad --broad-cutoff 0.1) merge nearby enriched regions into domains and suit diffuse marks like H3K27me3 or H3K36me3; MACS2 can also report gapped peaks, broad domains that contain at least one narrow sub-peak.

What q-value cutoff should I use in MACS2?

MACS2's default is 0.05. For narrow peaks (transcription-factor ChIP-seq), 0.01 is the tighter, more commonly recommended threshold; broad marks typically stay at 0.05. For single-replicate ATAC-seq without IDR filtering, some labs relax to 0.1 to avoid discarding real signal, but that trade-off increases false positives, so pair it with a biological sanity check.

What columns are in a narrowPeak file?

narrowPeak is a 10-column BED6+4 format: chrom, chromStart, chromEnd, name, score, strand, signalValue, pValue (-log10), qValue (-log10, FDR-corrected), and peak (summit offset from chromStart, or -1 if not called). The pValue and qValue columns are what you filter or rank peaks on, not the raw score column.

MACS2 or MACS3, which should I use?

MACS2 (currently v2.2.7.1) is still the version most published pipelines and tutorials build around, and it's the safer default if you need to match an existing protocol exactly. MACS3 is the actively maintained successor; check its release notes against your assay before swapping it into an established pipeline, since exact parameter parity with MACS2 isn't guaranteed.

Related pages

Related reading on the blog

Sources

  1. Identifying peaks in *-seq data using shape information — Peak calling methodology, Poisson modeling, and sources of false peaks
  2. MACS2 peak calling details and parameters — MACS2 parameter explanations, p-value/q-value thresholds, narrow vs broad peak calling, --keep-dup
  3. MACS2 parallel peak calling tutorial — Practical MACS2 commands and fragment-size parameter usage referenced in the example
  4. Investigating ChIP-seq: Peak calling with MACS — Shift model, FDR correction, narrow/broad/gapped peaks, ENCODE3 ATAC-seq parameters
  5. UCSC Genome Browser narrowPeak Format — narrowPeak BED6+4 column specification