Bio Chipseq Peak Calling

Version Compatibility

Reference examples tested with: MACS3 3.0.4+, MACS2 2.2.9+, HOMER 4.11+, SPP 1.16+, samtools 1.19+, bedtools 2.31+, IDR 2.0.4+.

Before running, verify versions: <tool> --version and <tool> --help to confirm flags. If a flag is missing, check the changelog — MACS2->MACS3 is API-compatible for callpeak but predictd, bdgpeakcall, and hmmratac differ.

ChIP-seq Peak Calling

"Identify protein-DNA binding sites from ChIP-seq alignments" -> Detect statistically enriched genomic regions by comparing IP signal to input control (or genomic background), with peak shape (narrow/broad) determined by target biology (TF vs histone mark).

• CLI (ENCODE TF default): macs2 callpeak -t chip.bam -c input.bam -f BAM -g hs -n sample --keep-dup all -p 1e-2
• CLI (ENCODE histone default): same with --broad --broad-cutoff 0.1 for H3K27me3, H3K9me3, H3K36me3
• CLI (alternative): macs3 callpeak ... (API-identical, active development), HOMER findPeaks tags/ -style histone -i input_tags/, SPP via phantompeakqualtools wrapper

ENCODE TF pipeline still uses SPP for peak ranking + IDR, with MACS2 producing the signal tracks. Histone pipeline uses MACS2 + naive overlap (IDR is too conservative for histone signal dynamic range). MACS3 is the actively maintained successor; MACS2 receives only bug fixes.

Critical Pre-Call Validation

Before any peak calling, three things must be true or the output is unreliable:

1. Antibody validated — KO/KD orthogonal control, peptide-array specificity for histone modifications, or vendor-provided CRISPR-validated lot (Epicypher, CST). "ChIP-grade" marketing is not validation. See chipseq-qc.
2. Fragment-size distribution is sane — TF ChIP should show sub-nucleosomal (~50-100 bp) enrichment; histone ChIP should show clean mono- (~150) and di-nucleosomal (~300) peaks. Flat distribution = over-sonication; rescue is impossible. Check via samtools view -f 0x2 sample.bam | awk '{print $9}' | sort | uniq -c.
3. Input control matches — Sonicated input is biased toward open chromatin; MNase input toward nucleosomes. Input from a different library prep batch or fragmentation method introduces bias that subtraction cannot fix.

Algorithmic Taxonomy

For CUT&RUN / CUT&Tag specifically, see chip-seq/cut-and-run-tag — protocol differences (lower depth, IgG-only control, E. coli spike-in carryover) drive different caller choice (MACS2 + SEACR consensus, not MACS3 alone).

Decision: Narrow vs Broad

Driven by target biology, not preference. Calling broad mode does not make a sharp signal broad; it changes how MACS stitches adjacent enrichment.

For HOMER: use -style histone for ALL histone marks (Omnipeak benchmark, Shpynov & Artyomov 2026 NAR 54:gkaf1454); -style factor ONLY for transcription factors.

Decision: Model vs --nomodel

MACS2/3 fragment-size modeling needs ≥100 paired plus/minus enrichment regions within --mfold (default [5, 50]). Silent failure produces wrong fragment size and warped peaks — always inspect _model.r output.

When --nomodel is required, choose --extsize in priority order: (1) cross-correlation estimate from phantompeakqualtools (ENCODE standard, gives NSC/RSC simultaneously); (2) macs3 predictd -i chip.bam -g hs and read stderr; (3) mark-type fallback (147 for nucleosome-proximal marks, 200 for broader marks).

Effective Genome Size — Often Wrong, Always Matters

-g hs (2.7e9) and -g mm (1.87e9) are decade-old approximations. Modern read-length-matched values (deepTools effectiveGenomeSize table):

Wrong size shifts every q-value but rarely peak ranks. For subset data (single chromosome, targeted), provide numeric -g <bp>; the shorthand inflates lambda_BG by 60× and produces false positives at low-signal regions.

Hyper-ChIPable Regions Are a Persistent Artifact

Teytelman 2013 (PNAS) and Park 2013 (PLoS One) demonstrated that highly-transcribed genes (rRNA, tRNA, histone gene cluster, snoRNA hosts, mitochondrial-encoded genes, abundant housekeeping loci) appear "bound" in ChIP-seq with untagged GFP, no antibody, or non-existent targets. ENCODE blacklist v2 catches repeat-driven artifacts but NOT these hyper-ChIPable transcribed regions.

Always interpret peaks at rRNA loci, tRNA clusters, replication-dependent histone genes (HIST1/2 clusters), mitochondrial DNA, and the top-1% input-signal regions with skepticism. For rigorous claims: (1) require motif enrichment at the peak (artifact has no motif); (2) require KO/KD signal loss; (3) build a cell-type-specific blacklist from the top 1% of input signal and intersect-out.

Pipeline Reference: ENCODE TF vs Histone

TF pipeline (uses SPP for peak ranking):

# Per-replicate (loose) — IDR tightens downstream
macs2 callpeak -t rep1.tagAlign.gz -c input.tagAlign.gz \
    -f BED -g hs -n rep1 \
    --nomodel --shift 0 --extsize {fraglen_from_xcor} \
    --keep-dup all -B --SPMR -p 1e-2

# Repeat for rep2, pooled, and pseudoreplicates (split each rep into halves)
# Score peaks by signalValue, sort, run IDR (see Replicate Handling below)

Histone pipeline (uses MACS2 broad / narrow + naive overlap):

# Broad marks: H3K27me3, H3K9me3, H3K36me3
macs2 callpeak -t rep1.tagAlign.gz -c input.tagAlign.gz \
    -f BED -g hs -n rep1 \
    --broad --broad-cutoff 0.1 \
    --nomodel --shift 0 --extsize {fraglen} \
    --keep-dup all -B --SPMR -p 1e-2

# Naive overlap: a peak passes if it appears in ≥2 of N replicates
# with ≥40% reciprocal overlap (ENCODE default, often misquoted as 50%)
bedtools intersect -a rep1.broadPeak -b rep2.broadPeak -f 0.40 -r -u > naive_overlap.bed

--keep-dup all is intentional in the ENCODE pattern: duplicates were already filtered upstream by MarkDuplicates + samtools view -F 1804 -q 30. -p 1e-2 is permissive because IDR (TF) or overlap (histone) tightens downstream.

Replicate Handling: IDR vs Naive Overlap

ENCODE rules (Landt 2012 Genome Res):

TFs use IDR. Run on signal-ranked peaks (sort by -k8,8nr p-value; -k7,7nr signal works for SPP but breaks for MACS pile-up if libraries differ).

sort -k8,8nr rep1.narrowPeak > rep1.sorted
sort -k8,8nr rep2.narrowPeak > rep2.sorted

idr --samples rep1.sorted rep2.sorted \
    --input-file-type narrowPeak --rank p.value \
    --idr-threshold 0.05 --output-file true_reps.idr --plot

ENCODE Nself/Nt consistency rule (often misremembered):

• Nt = IDR-passing peaks across true biological replicates (threshold 0.05)
• Nself (per rep) = IDR-passing peaks across pseudoreplicates of one library (threshold 0.10)
• Library passes if max(N1self, N2self) / min(N1self, N2self) ≤ 2 AND max(Nt, max(Nself)) / min(Nt, min(Nself)) ≤ 2
• Both ratios > 2: library rejected

Histones use naive overlap. IDR's high-vs-low-rank assumption breaks for histone dynamic range. Naive overlap: pool peaks, require each to appear in ≥2 replicates with ≥40% reciprocal overlap.

ENCODE 3 vs ENCODE 4 Differences

ENCODE 4 outputs are NOT numerically comparable to ENCODE 3 on the same BAM (blacklist change + genome size update shift peak counts ~3-10%).

Per-Tool Failure Modes

MACS2/3 -- Silent fragment-size model failure

Trigger: Sparse signal, low replicate depth, or saturated samples; _model.r plot never inspected.

Mechanism: Model needs ≥100 paired plus/minus enriched regions in --mfold range. Below threshold, MACS picks an arbitrary fragment size (often 50 or 1000 bp), producing miscentered or oversized peaks. Stderr shows a warning that gets ignored.

Symptom: Peak summits shifted relative to known motif positions by hundreds of bp; visual inspection in IGV shows peaks displaced from pile-up centers.

Fix: Inspect <sample>_model.r — if peaks look reasonable, accept; if degenerate, widen with --mfold 3 50 or switch to --nomodel --extsize <data-derived>. For consistency across samples in a study, always use --nomodel --extsize {fraglen} with cross-correlation-derived fraglen (ENCODE pattern).