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.

--broad

-f BAMPE

-style super

-C

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:

--broad

--broad

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

-style histone

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

--nomodel

--extsize 150 --shift -75

-f BAMPE

When using

--nomodel

--extsize

--mfold

--mfold 3 50

--nomodel

3. Fragment Size and Extension

When using

--nomodel

--extsize

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.

--broad

The ENCODE pipeline always uses

--nomodel --shift 0 --extsize {fraglen}

4. Input Format

TagAlign is a BED6 format widely used by ENCODE. MACS3 reads it with

-f BED

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

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:

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

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

--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

[5, 50]

[3, 50]

--nomodel

--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

-F 4

-L 4

-C 2

Histone Marks (H3K4me3, H3K27ac, H3K27me3, H3K36me3)

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

-style histone

-L 0

-L 0

Converting HOMER Output to BED

pos2bed.pl peaks.txt > peaks.bed

HOMER peak files use 1-indexed inclusive coordinates;

pos2bed.pl

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

HOMER Key Parameters

auto

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:

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

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

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

--keep-dup all

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