Calls ChIP-seq peaks with MACS3, MACS2, HOMER, or SPP across narrow (TF) and broad (histone) modes. Handles input control matching, fragment-size modeling vs --nomodel, effective genome size, ENCODE-style IDR vs naive overlap, hyper-ChIPable artifacts, and aligner-specific shifts. Use when calling peaks from ChIP-seq alignments, choosing between narrow vs broad mode for a histone mark, deciding model vs nomodel for low-depth data, applying ENCODE pseudoreplicate IDR, or reconciling MACS vs HOMER vs SPP results.
Instalar com Codex ou Claude Copie este prompt, cole no Codex, Claude ou outro assistente e deixe que ele revise a página da skill e instale para você.
Um comando direto ignora o prompt de revisão. Verifique a origem antes de executá-lo.
Instruções da origem · Visualização somente leitura
name
bio-chipseq-peak-calling
description
Calls ChIP-seq peaks with MACS3, MACS2, HOMER, or SPP across narrow (TF) and broad (histone) modes. Handles input control matching, fragment-size modeling vs --nomodel, effective genome size, ENCODE-style IDR vs naive overlap, hyper-ChIPable artifacts, and aligner-specific shifts. Use when calling peaks from ChIP-seq alignments, choosing between narrow vs broad mode for a histone mark, deciding model vs nomodel for low-depth data, applying ENCODE pseudoreplicate IDR, or reconciling MACS vs HOMER vs SPP results.
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 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:
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.
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.
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
Tool
Model
Treats fragments as
Strength
Fails when
MACS3/MACS2 callpeak
Dynamic local Poisson (max of genome-wide, 1kb, 5kb, 10kb lambda) + BH-FDR
Single-end shifts; PE fragments via BAMPE
Mature, fast, ENCODE-default, narrow + broad modes, integrated signal tracks
Confounds NFR with broad accessible domains; default narrow mode segments broad enrichment; assumes most genome NOT enriched (breaks for genome-wide marks)
ENCODE TF caller; integrated NSC/RSC QC; robust for sharp TF peaks
Underperforms for broad marks; older R codebase; phantompeakqualtools wrapper has R-version compatibility issues
HOMER -style factor
Fixed-width peaks + three sequential filters (control / local / clonal)
Tag positions; auto-estimated width
Fast on tag directories; clonal filter -C removes PCR-artifact peaks
Less calibrated p-values; fixed width clips variable-width factor binding
HOMER -style histone
Variable-width region stitching (500 bp blocks, 1000 bp gap merging); L=0 (no local enrichment)
Tag positions
Captures variable-width histone enrichment; Omnipeak benchmark (Shpynov & Artyomov 2026): outperforms -style factor for histone marks including H3K4me3
Less sensitive than MACS for very sharp TF binding
Genrich -y (ChIP mode)
q-value on log-transformed p-value, joint replicate model
Whole fragments (PE intervals)
Joint replicate analysis; chrM exclusion via -e chrM; auto blacklist via -E
Less peer-reviewed than MACS/SPP; thin literature; control handling less mature
MACS3 hmmratac
3-state HMM on fragment-size signal
Fragment-size classes
Best for ATAC, not ChIP
Wrong tool for ChIP; ChIP fragment-size distribution doesn't drive useful HMM states
SEACR (Meers 2019)
Empirical threshold on signal block totals
Bedgraph signal blocks
Designed for sparse CUT&RUN/CUT&Tag data; "stringent" mode with IgG strongly preferred
Not for traditional ChIP-seq (assumes near-zero background); see cut-and-run-tag
LanceOtron (Hentges 2022)
CNN trained on ENCODE peaks
bigWig signal
Competitive for both narrow and broad without parameter tuning
Newer; less validated; web-only or pip install
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.
Target
Mode
Why
Transcription factors (CTCF, p53, GATA1, FOXA1)
Narrow (default)
Discrete motif binding produces sharp peaks
H3K4me3, H3K27ac at promoters/enhancers
Narrow
Localized at regulatory elements
H3K4me1 at enhancers
Narrow or broad-cutoff 0.1
Variable; check published data for the cell type
H3K36me3, H3K79me2 (elongation)
Broad
Deposited across active gene bodies (5-50 kb domains)
H3K27me3, H3K9me3 (repressive)
Broad
Spread across 10-100+ kb domains
H4K20me3 (constitutive het)
Broad
Heterochromatin domains
Pol II (RNAPII)
Narrow at promoter + broad option for elongation profile
Two separate analyses if doing elongation biology
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.
Condition
Model?
Fallback
Whole-genome, ≥1M treatment reads, narrow TF
Yes
--mfold 3 50 if fails
Paired-end with -f BAMPE
N/A
Fragment size from mate pairs
Single chromosome or targeted capture
No
--nomodel --extsize <data-derived or mark default>
Low read count (<500k)
No
Same
Broad histone mark
Either
Mark-type default if no estimate available
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):
Genome
Read length
Effective size
hg38
50 bp
2.701e9
hg38
75 bp
2.748e9
hg38
100 bp
2.806e9
hg38
150 bp
2.862e9
mm10
50 bp
2.308e9
mm10
100 bp
2.467e9
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)
# 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).
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
Feature
ENCODE 3
ENCODE 4
TF peak ranker
SPP
SPP (unchanged)
Histone caller
MACS2
MACS2 (MACS3 not yet adopted)
Aligner
bwa-mem
bwa-mem (chromap evaluated; not yet swapped)
Blacklist
v1 (ENCODE DAC, unpublished resource)
v2 (Amemiya 2019)
TF significance
-p 1e-2 + IDR @ 0.05
Same
Histone significance
-p 1e-2 + naive overlap
Same
Effective genome size
hs/mm shorthand
deepTools read-length-tabulated
Pseudoreplicate IDR threshold
0.10 self-consistency
0.10 self-consistency
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).
MACS2/3 -- Confounded narrow vs broad on intermediate marks
Trigger: Marks of intermediate breadth (H3K4me1, H3K9ac) called with default narrow mode.
Symptom: Peak count 3-5× higher than published for same cell type; mean peak width < 200 bp at known enhancer regions.
Fix: For H3K4me1, try --broad --broad-cutoff 0.1 and compare; for H3K9ac, narrow mode typically OK. Always cross-reference published peak counts for the cell type and antibody lot.
MACS2/3 -- --call-summits double-counts
Trigger: Narrow mode + --call-summits flag.
Mechanism: MACS adds sub-peak summits at multi-mode pile-ups; broad-shouldered peaks get split into 2-3 entries.
Symptom: Peak count inflated; same genomic region appears as 2-3 adjacent peaks in narrowPeak output.
Fix: Drop --call-summits unless deliberately analyzing multi-mode binding (rare); merge bedtools merge -d 200 if needed post-hoc.
HOMER -- Wrong style for histones
Trigger:-style factor used for histone marks.
Mechanism: Factor mode uses fixed-width peaks with local enrichment filter -L 4 that eliminates broad signal.
Symptom: Far fewer peaks than expected for H3K4me3/H3K27ac/H3K27me3; missed enrichment at known regions.
Fix: Use -style histone for ALL histone marks (Omnipeak benchmark, Shpynov & Artyomov 2026); reserve -style factor for TFs only.
SPP / phantompeakqualtools -- R version incompatibility
Trigger: Running phantompeakqualtools wrapper script with R ≥ 4.0.
Mechanism: spp R package has unmaintained dependencies; some functions silently fail or return NaN for NSC/RSC.
Fix: Use conda env pinned to R 3.6 + spp 1.16; or use kundajelab/phantompeakqualtools fork (current); or substitute deepTools plotFingerprint for QC and MACS-derived fragment length.
Trigger: Using chromap (fast aligner) output as MACS input with --shift -75 --extsize 150.
Mechanism: chromap pre-applies a Tn5/cut-site shift before fragment output (designed for ATAC); ChIP cut-site reasoning doesn't apply but the shift still happens silently.
Symptom: Peaks shifted ~5-10 bp from bwa-mem output at the same locus.
Fix: When using chromap, drop downstream shift OR use chromap's --no-correction. For ChIP, bwa-mem or bowtie2 are safer defaults until ENCODE switches.
Reconciliation: When Callers Disagree
Pattern
Likely cause
Action
MACS finds peak; HOMER misses
HOMER local-enrichment filter (-L 4) removed it at low-signal regions; or -style factor clipped a histone peak
Re-run HOMER with -style histone -L 0 for histones; if persists, trust MACS