/macs3-peak-calling

name: "macs3-peak-calling"
description: "Poisson-model peak caller for ChIP-seq/ATAC-seq BAMs. MACS3 callpeak finds enriched regions (TF sites or histone marks) vs input/IgG; outputs BED narrowPeak/broadPeak for motif analysis, annotation, and differential binding. Use narrow peaks for TF ChIP-seq and ATAC-seq; broad for H3K27me3, H3K9me3, and other broad marks."
license: "BSD-3-Clause"

MACS3 — ChIP-seq and ATAC-seq Peak Caller

Overview

MACS3 (Model-based Analysis of ChIP-seq) identifies regions of significant read enrichment (peaks) from ChIP-seq, ATAC-seq, CUT&RUN, and CUT&TAG experiments. It models the fragment length distribution from paired-end data or estimates it from mono-nucleosomal read shifting in single-end data, then applies a Poisson model to identify fold-enrichment over an input/IgG control. MACS3 produces BED-format narrowPeak (for transcription factors) or broadPeak (for histone marks) files with signal and q-value tracks for visualization in IGV or UCSC Genome Browser.

When to Use

• Calling transcription factor binding peaks from ChIP-seq experiments (use `--nomodel --extsize 200` or let MACS3 estimate fragment length)
• Identifying open chromatin regions from ATAC-seq experiments (use `--nomodel --shift -100 --extsize 200 -f BAMPE`)
• Calling broad histone modification peaks (H3K27me3, H3K9me3, H3K36me3) with `--broad`
• Generating peak signal tracks (bedGraph/bigWig) for genome browser visualization with `-B --SPMR`
• Performing differential binding analysis: MACS3 peaks as input to DiffBind or DESeq2
• Use **HMMRATAC** (part of MACS3) for nucleosome-resolution ATAC-seq peak calling
• Use **SPP** or **HOMER** as alternatives; MACS3 is the ENCODE-recommended standard

Prerequisites

• **Python packages**: `macs3` (Python ≥ 3.8)
• **Input**: Sorted BAM files (with index) from ChIP-seq or ATAC-seq alignment (e.g., using STAR or Bowtie2)
• **Optional**: Input/IgG control BAM for background normalization

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v macs3` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run macs3` rather than bare `macs3`.

# Install with pip or conda
pip install macs3
# or
conda install -c bioconda macs3

# Verify
macs3 --version
# macs3 3.0.2

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: assayType
    kind: required
    source: user
    ask: "What produced these reads - a transcription-factor ChIP, a histone mark, or open-chromatin (ATAC)?"
    default: null

  - id: D2
    param: controlSample
    kind: required
    source: data
    ask: "Is there a matched input or IgG control for this sample?"
    default: "none - background is estimated locally instead"

  - id: D3
    param: effectiveGenomeSize
    kind: derived
    source: upstream
    ask: "Which genome were the reads aligned to?"
    default: "carried from the alignment reference"

  - id: D4
    param: peakShape
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Are the enriched regions sharp summits or broad domains?"
    default: "narrow for TF and ATAC; broad for spreading histone marks"

  - id: D5
    param: significanceThreshold
    kind: required
    source: user
    ask: "How strong must enrichment be before a region is called a peak?"
    default: "q-value 0.05"

  - id: D6
    param: fragmentModel
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Should fragment length be modelled from the data, or fixed because the assay's signal is not paired summits?"
    default: "model for ChIP; fixed shift/extension for ATAC"

  - id: D7
    param: duplicateHandling
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Are identical read positions PCR duplicates, or real independent insertions?"
    default: "one per position for ChIP; keep all for ATAC"

  - id: D8
    param: signalTracks
    kind: optional
    source: user
    ask: "Emit normalized coverage tracks alongside the peak calls?"
    default: "peaks only"

D1 is asked first because four later settings follow from it. ATAC in particular inverts three ChIP defaults at once - no fragment model, a negative shift, and duplicates kept - so running ATAC through ChIP defaults produces peaks that look plausible and are wrong.

Quick Start

# Call peaks for TF ChIP-seq (narrow peaks, with input control)
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample_tf \
    --outdir peaks/ \
    -q 0.05

# Output: peaks/sample_tf_peaks.narrowPeak
wc -l peaks/sample_tf_peaks.narrowPeak

Workflow

Step 1: Prepare Input BAM Files

MACS3 requires sorted, indexed BAM files from genome alignment.

# Sort and index ChIP and control BAMs (if not already done)
samtools sort -@ 8 chip_raw.bam -o chip.bam
samtools sort -@ 8 input_raw.bam -o input.bam
samtools index chip.bam
samtools index input.bam

# Check read counts
echo "ChIP reads: $(samtools view -c -F 4 chip.bam)"
echo "Input reads: $(samtools view -c -F 4 input.bam)"

Step 2: Call Narrow Peaks (TF ChIP-seq)

Use the default mode for transcription factor binding site identification.

# TF ChIP-seq with input control
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n tf_chip \
    --outdir peaks/ \
    -q 0.05 \
    --keep-dup auto

echo "Peaks called: $(wc -l < peaks/tf_chip_peaks.narrowPeak)"
echo "Summit file: peaks/tf_chip_summits.bed"

# Without input control (less recommended)
macs3 callpeak \
    -t chip.bam \
    -f BAM \
    -g hs \
    -n tf_noinput \
    --outdir peaks/ \
    --nolambda

Step 3: Call Broad Peaks (Histone Marks)

Use `--broad` for spread histone modifications like H3K27me3 or H3K36me3.

# H3K27me3 broad