bio-chipseq-peak-calling

name	bio-chipseq-peak-calling
description	ChIP-seq peak calling using MACS3 and HOMER findPeaks. Call narrow peaks for transcription factors or broad peaks for histone modifications. Supports single-caller and multi-caller consensus approaches, input control, fragment size modeling, and various output formats. Use when calling peaks from ChIP-seq alignments.
tool_type	cli
primary_tool	macs3

Version Compatibility

Reference examples tested with: MACS2 2.2+, MACS3 3.0+, HOMER 4.11+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Peak Calling

"Call peaks from my ChIP-seq data" → Identify significantly enriched regions (narrow peaks for TFs, broad peaks for histone marks) by comparing IP signal to input control.

• CLI (MACS3): macs3 callpeak -t chip.bam -c input.bam -f BAM -g hs -n sample
• CLI (HOMER): makeTagDirectory tags/ chip.bam then findPeaks tags/ -style factor -i input_tags/ -o peaks.txt

Choosing a Peak Caller

MACS3 and HOMER use fundamentally different statistical approaches. MACS3 builds a dynamic local Poisson model (taking the maximum of genome-wide, 1kb, 5kb, and 10kb background estimates). HOMER applies three independent sequential filters: control enrichment, local enrichment, and clonal signal complexity. Neither is universally superior — the choice depends on the analysis context.

Context	Recommended	Why
ENCODE compliance / IDR workflows	MACS3	ENCODE standard; IDR designed around MACS2/SPP output
Broad marks (H3K27me3, H3K36me3)	MACS3 --broad	Two-tier Poisson produces gapped peaks with internal structure
Paired-end data	MACS3 with -f BAMPE	Native fragment pileup without read extension
Integrated motif + annotation workflow	HOMER	Peak calling, motif discovery, and annotation share one tag directory
Super-enhancer identification	HOMER -style super	Native ROSE-like algorithm from H3K27ac data
High-confidence peak sets	Both + consensus	Multi-caller intersection removes tool-specific false positives
Repetitive genomes or PCR artifacts	HOMER	Clonal filter (-C) catches false positives over repeats

When both tools are available and runtime permits, running both and taking the intersection produces a higher-confidence peak set than either caller alone. See Multi-Caller Consensus below.

MACS3 is the actively developed successor to MACS2. Commands are identical except the binary name. MACS2 is in maintenance mode.

Decision Framework

Before running any peak calling command, assess the data and choose parameters accordingly. The three key decisions are peak mode, model building strategy, and genome size.

1. Peak Mode: Narrow vs Broad

The choice depends on the biology of the target, not a preference:

Target type	Examples	Mode	Why
Transcription factors	CTCF, p53, GATA1	Narrow (default)	TFs bind discrete motif sites, producing sharp peaks
Active promoter marks	H3K4me3, H3K27ac	Narrow (default)	These marks localize to narrow regulatory elements
Elongation/body marks	H3K36me3, H3K79me2	--broad	Deposited across gene bodies, forming wide domains
Repressive marks	H3K27me3, H3K9me3	--broad	Spread across large chromatin domains

If uncertain about the target's enrichment pattern, check published ENCODE data for that mark or search current literature before proceeding.

HOMER note: The MACS3 modes above do not map directly to HOMER styles. For HOMER, use -style factor only for transcription factors. All histone marks — including narrow marks like H3K4me3 and H3K27ac — perform better with -style histone (benchmarked in Omnipeak, NAR 2025). See the HOMER section below.

2. Model Building vs --nomodel

MACS3 estimates fragment size by finding paired plus/minus strand peaks. This requires at least 100 such pairs within the --mfold enrichment range. Use --nomodel when model building is expected to fail:

Condition	Use --nomodel?	Why
Whole-genome, >1M treatment reads	No — let MACS3 model	Enough data for reliable fragment size estimation
Single chromosome or small region	Yes	Too few enriched regions for 100 paired peaks
Low read count (<500k treatment)	Yes	Sparse signal makes modeling unreliable
ATAC-seq or DNase-seq	Yes, with --extsize 150 --shift -75	Open chromatin has no directional shift to model
Paired-end data (-f BAMPE)	N/A	Fragment size comes from mate pairs, no modeling needed

When using --nomodel, set --extsize based on the mark type (see table below). When model building fails unexpectedly, try widening --mfold (e.g., --mfold 3 50) before falling back to --nomodel.

3. Fragment Size and Extension

When using --nomodel, the --extsize value should come from data when possible, not a generic default:

1. Best: Estimate from strand cross-correlation using phantompeakqualtools (SPP). This is the ENCODE standard — the same analysis produces QC metrics (NSC/RSC). See chipseq-qc for cross-correlation details.
2. Good: Estimate using macs3 predictd -i chip.bam -g hs --outdir . and read the fragment size from stderr. Requires sufficient enriched regions (same constraint as model building).
3. Fallback: Use mark-type defaults from the table below.

Fragment sizes vary across experiments — ENCODE H3K4me3 datasets show cross-correlation estimates from 95-245bp within a single experiment. Using an estimated value is always preferable to a default.

Mark type	Fallback --extsize	Rationale
TFs (CTCF, p53)	150-200	Typical ChIP fragment size range
H3K4me3, H3K27ac	147	Nucleosome core particle is 147bp; these marks are nucleosome-proximal
H3K4me1	200	Slightly broader enhancer marks
H3K36me3, H3K79me2	200-300	Gene body marks, broader signal
H3K27me3, H3K9me3	200	Extension less critical since --broad links subpeaks

The ENCODE pipeline always uses --nomodel --shift 0 --extsize {fraglen} where fraglen comes from cross-correlation, even when model building would succeed. This ensures consistency between peak calling and signal track generation.

4. Input Format

File type	Extension	-f flag	Notes
Aligned BAM (SE)	.bam	BAM	Most common
Aligned BAM (PE)	.bam	BAMPE	Uses actual fragment sizes from mate pairs
tagAlign / BED6	.tagAlign.gz, .bed.gz	BED	ENCODE standard format; stores chr, start, end, name, score, strand
SAM	.sam	SAM	Rarely used directly

TagAlign is a BED6 format widely used by ENCODE. MACS3 reads it with -f BED (not a special tagAlign flag).

Basic Peak Calling

macs3 callpeak -t chip.bam -c input.bam -f BAM -g hs -n sample --outdir peaks/

Without Input Control

macs3 callpeak -t chip.bam -f BAM -g hs -n sample --outdir peaks/

Calling without input control uses genomic background only. This is less accurate — always prefer matched input/IgG when available.

Narrow Peaks (TF, H3K4me3, H3K27ac)

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample_narrow \
    --outdir peaks/ \
    -q 0.05

Broad Peaks (H3K36me3, H3K27me3, H3K9me3)

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample_broad \
    --outdir peaks/ \
    --broad \
    --broad-cutoff 0.1

Paired-End Data

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAMPE \
    -g hs \
    -n sample_pe \
    --outdir peaks/

BAMPE uses actual fragment sizes from mate pairs — no model building or --extsize needed.

tagAlign / BED Input

macs3 callpeak \
    -t treatment.tagAlign.gz \
    -c control.tagAlign.gz \
    -f BED \
    -g hs \
    -n sample \
    --outdir peaks/

MACS3 accepts gzipped tagAlign directly. No decompression needed.

Multiple Replicates

macs3 callpeak \
    -t rep1.bam rep2.bam rep3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n pooled \
    --outdir peaks/

Genome Size

Use built-in shortcuts for model organisms, or provide a numeric value for non-model organisms, single-chromosome data, or targeted captures:

Genome	Flag	Effective Size
Human (whole)	hs	2.7e9
Mouse (whole)	mm	1.87e9
C. elegans	ce	9e7
D. melanogaster	dm	1.2e8
Human chr21 only	46700000	~46.7M mappable bases
Custom/targeted	numeric	Sum of mappable bases in target regions

For single-chromosome or targeted data, always provide the numeric effective genome size rather than the whole-genome shortcut. MACS computes genome-wide background as lambda_BG = (control_reads x fragment_size) / genome_size. Using hs (2.9e9) on chr21-only data (~46.7M) deflates lambda_BG by ~62x, making the background floor artificially low and allowing false positives in regions where local lambda is also low.

Fixed Fragment Size (--nomodel)

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --nomodel \
    --extsize 150 \
    -n sample \
    --outdir peaks/

See the decision framework above for when to use --nomodel and how to choose --extsize.

Rescuing Failed Model Building

If model building fails with "needs at least 100 paired peaks":

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --mfold 3 50 \
    -n sample \
    --outdir peaks/

Widening --mfold from the default [5, 50] to [3, 50] lowers the enrichment threshold for selecting model-building regions, which can recover enough paired peaks. If this still fails, use --nomodel with an appropriate --extsize.

Generate Signal Tracks

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    --outdir peaks/ \
    -B \
    --SPMR

sort -k1,1 -k2,2n peaks/sample_treat_pileup.bdg > peaks/sample.sorted.bdg
bedGraphToBigWig peaks/sample.sorted.bdg chrom.sizes peaks/sample.bw

Cutoff Analysis

macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --cutoff-analysis \
    -n sample \
    --outdir peaks/

Generates a table of peak counts at various q-value cutoffs. Useful for choosing a threshold when expected peak counts are unclear.

HOMER Peak Calling

HOMER peak calling requires two steps: creating a tag directory (which estimates fragment size via strand autocorrelation), then calling peaks against it.

Tag Directories

makeTagDirectory chip_tags/ chip.bam
makeTagDirectory input_tags/ input.bam

For tagAlign/BED input, specify format explicitly:

makeTagDirectory chip_tags/ treatment.tagAlign.gz -format bed
makeTagDirectory input_tags/ control.tagAlign.gz -format bed

Tag directories store pre-processed alignment data and compute QC metrics (fragment length, strand enrichment) automatically.

Transcription Factors

findPeaks chip_tags/ -style factor -i input_tags/ -gsize 2.7e9 -o peaks.txt

-style factor uses fixed-width peaks with all three filters active (control enrichment -F 4, local enrichment -L 4, clonal filtering -C 2). Peak width is auto-estimated from tag autocorrelation. Use for point-source TF binding (CTCF, p53, GATA1).

Histone Marks (H3K4me3, H3K27ac, H3K27me3, H3K36me3)

findPeaks chip_tags/ -style histone -i input_tags/ -gsize 2.7e9 -o regions.txt

-style histone uses variable-width region stitching (500bp building blocks stitched within 1000bp) and disables local enrichment filtering (-L 0). Use for all histone marks, including narrow marks like H3K4me3 and H3K27ac. Benchmarking (Omnipeak, NAR 2025) shows histone mode outperforms factor mode for H3K4me3 because it captures the variable-width enrichment around modified nucleosomes rather than clipping to a fixed-width center. For broad marks (H3K27me3, H3K36me3), -L 0 is essential since local enrichment filtering would eliminate the signal being sought.

Converting HOMER Output to BED

pos2bed.pl peaks.txt > peaks.bed

HOMER peak files use 1-indexed inclusive coordinates; pos2bed.pl converts to 0-indexed BED format. For manual conversion:

awk 'BEGIN{OFS="\t"} !/^#/ && NF>=5 {print $2, $3, $4, $1, $8}' peaks.txt > peaks.bed

HOMER Key Parameters

Parameter	Default	Description
-style	factor	Peak mode: factor, histone, super, groseq, tss, dnase
-i	none	Control tag directory
-gsize	2e9	Effective genome size
-size	auto	Peak size (auto for factor; 500 for histone)
-F	4.0	Fold enrichment over control
-L	4.0	Fold enrichment over local region (0 in histone mode)
-C	2.0	Clonal signal filter (0 to disable)
-fdr	0.001	FDR threshold
-o	stdout	Output file (auto writes to tag directory)

Multi-Caller Consensus

When both MACS3 and HOMER are available, running both and intersecting results is the recommended approach for final peak sets. MACS3's local Poisson model and HOMER's three-filter approach have different error profiles, so their intersection removes tool-specific false positives. For quick exploration with a single tool, MACS3 alone is sufficient.

bedtools intersect -a macs3_peaks.narrowPeak -b homer_peaks.bed -wa -u > consensus_peaks.bed

For lenient matching (peaks within 500bp), use bedtools window:

bedtools window -a macs3_peaks.narrowPeak -b homer_peaks.bed -w 500 | cut -f1-10 | sort -k1,1 -k2,2n | uniq > consensus_peaks.bed

Single-caller vs consensus: Single-caller (MACS3 alone) is appropriate for exploratory analysis, quick surveys, or when only one tool is installed. Consensus is preferred for final peak sets used in downstream motif analysis, functional enrichment, cross-condition comparison, or publication.

Consensus vs IDR vs naive overlap: These address different questions:

• IDR measures replicate reproducibility (same caller, multiple replicates). ENCODE standard for TF ChIP-seq.
• Naive overlap identifies peaks found across replicates. ENCODE standard for histone marks — IDR is too conservative for the variable dynamic range of histone signal.
• Multi-caller consensus measures algorithmic agreement (multiple callers, same data). Widely used for higher-confidence peak sets.

For publication-quality analysis, consider both replicate thresholding (IDR or naive overlap) and multi-caller consensus.

Sanity-Checking Results

After peak calling, verify that results are biologically plausible before proceeding to downstream analysis:

Target	Typical peak count (whole genome)	Typical peak width
TFs (CTCF, p53)	10,000-80,000	200-500 bp
H3K4me3	20,000-50,000	500-2,000 bp
H3K27ac	30,000-80,000	500-3,000 bp
H3K27me3	5,000-30,000 broad domains	10-100+ kb
H3K36me3	10,000-30,000 broad domains	5-50 kb

For subset data (e.g., single chromosome), scale expectations proportionally — chr21 is ~1.5% of the human genome, so expect ~1.5% of whole-genome peak counts.

Red flags that indicate parameter problems:

• Zero peaks: Check genome size flag (using hs on subset data?), input format (-f), and whether input control is properly matched
• Orders of magnitude too many peaks: Control sample may be swapped with treatment, or q-value cutoff is too permissive
• Orders of magnitude too few peaks: Genome size may be too large for the data, or model building failed silently (check stderr)

When results are unexpected, re-read the MACS3 stderr output — it reports tag counts, redundancy rates, fragment size estimates, and model building status. Search current MACS3 documentation or literature if the issue is not covered here.

Output Files

File	Description
*_peaks.narrowPeak	Peak coordinates (BED6+4)
*_peaks.broadPeak	Broad peak coordinates
*_summits.bed	Peak summit positions
*_model.r	R script for model visualization
*_treat_pileup.bdg	Treatment signal (with -B)
*_control_lambda.bdg	Control signal (with -B)

narrowPeak Format

chr1  100  200  peak_1  100  .  5.2  10.5  8.3  50

Columns: chr, start, end, name, score, strand, signalValue, pValue, qValue, peak

Convert narrowPeak to BED

Many downstream tools and pipelines expect simple BED format (chr, start, end, name, score). Extract the first 5 columns from narrowPeak:

cut -f1-5 peaks.narrowPeak > peaks.bed

Filter Peaks

awk '$9 > 2' peaks.narrowPeak > peaks.filtered.narrowPeak  # -log10(q) > 2 means q < 0.01
sort -k7,7nr peaks.narrowPeak > peaks.sorted.narrowPeak

Key Parameters

Parameter	Default	Description
-t	required	Treatment file(s)
-c	none	Control file(s)
-f	AUTO	Format: BAM, BAMPE, BED, SAM
-g	hs	Effective genome size (shortcut or numeric)
-n	NA	Output prefix
-q	0.05	Q-value cutoff
-p	none	P-value cutoff (overrides -q)
--broad	false	Broad peak calling
--broad-cutoff	0.1	Q-value cutoff for broad regions
--nomodel	false	Skip fragment size modeling
--extsize	200	Read extension size (with --nomodel)
--mfold	5 50	Enrichment range for model building
-B	false	Generate bedGraph
--SPMR	false	Signal per million reads
--nolambda	false	Use local background only (skip genome-wide lambda)

ENCODE Pipeline Reference

The ENCODE ChIP-seq pipeline uses these MACS2 parameters as its standard:

macs2 callpeak -t chip.tagAlign.gz -c input.tagAlign.gz \
    -f BED -g {gsize} -p 1e-2 \
    --nomodel --shift 0 --extsize {fraglen} \
    --keep-dup all -B --SPMR

Key differences from default MACS3 settings: uses -p 0.01 (p-value, not q-value) because IDR or naive overlap thresholding is applied downstream to filter peaks across replicates; uses --keep-dup all because duplicates are handled upstream during BAM processing; fragment length comes from phantompeakqualtools cross-correlation, not MACS model building. For histone marks, ENCODE applies naive overlap across replicates rather than IDR.

Related Skills

• peak-annotation - Annotate peaks to genes (ChIPseeker or HOMER annotatePeaks.pl)
• motif-analysis - De novo and known motif enrichment (HOMER findMotifsGenome.pl, MEME)
• differential-binding - Compare peaks between conditions
• chipseq-qc - QC metrics including FRiP, NSC/RSC, IDR
• super-enhancers - Super-enhancer identification with ROSE or HOMER -style super
• alignment-files/sam-bam-basics - BAM file preparation
• chipseq-visualization - Visualize peaks and signal tracks
• genome-intervals/interval-arithmetic - BED intersection, merging, and window operations