macs3-peak-calling

MACS3 — ChIP-seq 和 ATAC-seq 峰识别工具

概述

MACS3（基于模型的 ChIP-seq 分析）可从 ChIP-seq、ATAC-seq、CUT&RUN 和 CUT&TAG 实验中识别显著富集的读段区域（峰）。它会根据双端数据对片段长度分布建模，或根据单端数据中单核小体读段的偏移估算片段长度，然后应用泊松模型，识别相对于 input/IgG 对照的富集倍数。MACS3 会生成 BED 格式的 narrowPeak（适用于转录因子）或 broadPeak（适用于组蛋白标记）文件，其中包含信号和 q 值轨道，可用于在 IGV 或 UCSC Genome Browser 中可视化。

适用场景

• 从 ChIP-seq 实验中识别转录因子结合峰（使用 --nomodel --extsize 200，或让 MACS3 估算片段长度）
• 从 ATAC-seq 实验中识别开放染色质区域（使用 --nomodel --shift -100 --extsize 200 -f BAMPE）
• 使用 --broad 识别宽广的组蛋白修饰峰（H3K27me3、H3K9me3、H3K36me3）
• 使用 -B --SPMR 生成峰信号轨道（bedGraph/bigWig），用于基因组浏览器可视化
• 进行差异结合分析：将 MACS3 峰作为 DiffBind 或 DESeq2 的输入
• 使用 HMMRATAC（MACS3 的一部分）进行核小体分辨率的 ATAC-seq 峰识别
• 也可使用 SPP 或 HOMER；MACS3 是 ENCODE 推荐的标准工具

前置条件

• Python 包：macs3（Python ≥ 3.8）
• 输入：来自 ChIP-seq 或 ATAC-seq 比对的已排序 BAM 文件（带索引）（例如使用 STAR 或 Bowtie2）
• 可选：用于背景归一化的 input/IgG 对照 BAM 文件

安装前请检查：当前环境中可能已经有此工具（例如位于 pixi / conda 环境中）。先运行 command -v macs3；如果返回路径，则跳过以下安装命令。在 pixi 项目中运行时，请使用 pixi run macs3 调用工具，而不要直接运行 macs3。

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

# Verify
macs3 --version
# macs3 3.0.2

预先确认

编写任何分析代码前，请先与用户确认以下事项。

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，因为后续四项设置都取决于它。ATAC 尤其需要同时反转三项 ChIP 默认设置——不使用片段模型、采用负向偏移、保留重复读段——因此用 ChIP 默认设置运行 ATAC 会产生看似合理、实际却错误的峰。

快速开始

# 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

工作流

步骤 1：准备输入 BAM 文件

MACS3 需要来自基因组比对的已排序、带索引 BAM 文件。

# 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)"

步骤 2：识别窄峰（转录因子 ChIP-seq）

使用默认模式识别转录因子结合位点。

# 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

步骤 3：识别宽峰（组蛋白标记）

对于 H3K27me3 或 H3K36me3 等分布较广的组蛋白修饰，使用 --broad。

# H3K27me3 broad histone mark
macs3 callpeak \
    -t h3k27me3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n h3k27me3 \
    --outdir peaks/ \
    --broad \
    --broad-cutoff 0.1 \
    -q 0.05

echo "Broad peaks: $(wc -l < peaks/h3k27me3_peaks.broadPeak)"

# H3K4me3 (sharp mark — use narrow peaks)
macs3 callpeak \
    -t h3k4me3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n h3k4me3 \
    --outdir peaks/ \
    -q 0.05

步骤 4：识别 ATAC-seq 峰

ATAC-seq 需要针对 Tn5 插入位点进行特殊处理。

# ATAC-seq with paired-end BAM (recommended)
macs3 callpeak \
    -t atac.bam \
    -f BAMPE \
    -g hs \
    -n atac_sample \
    --outdir peaks/ \
    --nomodel \
    --nolambda \
    -q 0.05 \
    --keep-dup all

echo "ATAC peaks: $(wc -l < peaks/atac_sample_peaks.narrowPeak)"

# Single-end ATAC-seq: shift reads to center on Tn5 cut site
macs3 callpeak \
    -t atac_se.bam \
    -f BAM \
    -g hs \
    -n atac_se \
    --outdir peaks/ \
    --nomodel \
    --shift -100 \
    --extsize 200 \
    --keep-dup all

步骤 5：生成用于可视化的信号轨道

生成 bedGraph 和 bigWig 文件，用于基因组浏览器可视化。

# Generate bedGraph normalized to million reads (SPMR)
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n chip_track \
    --outdir tracks/ \
    -B \
    --SPMR \
    --keep-dup auto

# Convert bedGraph to bigWig for IGV/UCSC
# Requires bedGraphToBigWig and chrom.sizes
sort -k1,1 -k2,2n tracks/chip_track_treat_pileup.bdg > tracks/chip_sorted.bdg
bedGraphToBigWig tracks/chip_sorted.bdg genome/hg38.chrom.sizes tracks/chip.bw

echo "BigWig track: tracks/chip.bw"

步骤 6：注释并分析峰

解析 narrowPeak 输出，并将峰注释到基因组特征。

import pandas as pd

# Load narrowPeak file
# Columns: chrom, start, end, name, score, strand, signalValue, pValue, qValue, peak
cols = ["chrom", "start", "end", "name", "score", "strand",
        "signalValue", "pValue", "qValue", "peak"]
peaks = pd.read_csv("peaks/tf_chip_peaks.narrowPeak", sep="\t",
                    header=None, names=cols)

print(f"Total peaks: {len(peaks)}")
print(f"Peaks on chr1: {(peaks['chrom'] == 'chr1').sum()}")
print(f"Median peak width: {(peaks['end'] - peaks['start']).median():.0f} bp")
print(f"Peaks with q-value < 0.01: {(peaks['qValue'] > 2).sum()}")  # -log10(q) > 2

# Filter high-confidence peaks
high_conf = peaks[peaks["qValue"] > 2].copy()  # q < 0.01
high_conf["width"] = high_conf["end"] - high_conf["start"]
print(f"\nHigh-confidence peaks: {len(high_conf)}")
high_conf.to_csv("high_confidence_peaks.bed", sep="\t", index=False, header=False,
                 columns=["chrom", "start", "end", "name", "score", "strand"])

关键参数

参数	默认值	范围/选项	作用
-t / --treatment	必需	BAM/BED/SAM	ChIP 或 ATAC 处理文件
-c / --control	—	BAM/BED/SAM	Input/IgG 对照；若无对照，请省略 --nolambda
-g / --gsize	必需	hs、mm、ce、dm 或整数	有效基因组大小；hs=2.7e9（人类），mm=1.87e9（小鼠）
-q / --qvalue	0.05	0–1	峰调用的 FDR 阈值
-p / --pvalue	—	0–1	P 值截止值（需要严格控制时使用此选项，而非 q 值）
--broad	关闭	标志	为弥散的组蛋白标记调用宽峰
--broad-cutoff	0.1	0–1	宽区域合并的 Q 值截止值
--nomodel	关闭	标志	跳过片段长度建模；ATAC-seq 必需
--extsize	200	50–1000	设置 --nomodel 时的片段延伸长度
--shift	0	-500–500	读取位置偏移，单位为 bp；ATAC-seq 使用 --extsize 200 时设为 -100
--keep-dup	1	auto、all、整数	重复项处理方式；auto 使用泊松模型，all 保留全部（ATAC-seq）
-B / --bdg	关闭	标志	写入 bedGraph 信号轨道
--SPMR	关闭	标志	将 bedGraph 归一化为每百万次读取的信号

常用方案

方案 1：对多个样本批量调用峰

#!/bin/bash
# Call peaks for multiple ChIP-seq samples with the same input
INPUT="input.bam"
GENOME="hs"
OUTDIR="peaks"
mkdir -p "$OUTDIR"

SAMPLES=(H3K4me3 H3K27ac H3K27me3 CTCF)
MODES=(narrow narrow broad narrow)

for i in "${!SAMPLES[@]}"; do
    sample="${SAMPLES[$i]}"
    mode="${MODES[$i]}"
    echo "Calling peaks: $sample ($mode)"
    
    if [ "$mode" == "broad" ]; then
        BROAD_FLAG="--broad --broad-cutoff 0.1"
    else
        BROAD_FLAG=""
    fi
    
    macs3 callpeak \
        -t "${sample}.bam" \
        -c "$INPUT" \
        -f BAM \
        -g "$GENOME" \
        -n "$sample" \
        --outdir "$OUTDIR" \
        $BROAD_FLAG \
        -q 0.05 \
        --keep-dup auto
    
    echo "$sample: $(wc -l < $OUTDIR/${sample}_peaks.*Peak) peaks"
done

方案 2：使用 IDR（不可重复发现率）获取可重复的峰

# Call peaks on individual replicates (lenient thresholds for IDR)
for rep in rep1 rep2; do
    macs3 callpeak \
        -t "chip_${rep}.bam" \
        -c input.bam \
        -f BAM \
        -g hs \
        -n "tf_${rep}" \
        --outdir peaks/ \
        -p 0.1 \
        --keep-dup auto
done

# Run IDR to find reproducible peaks
# pip install idr
idr --samples peaks/tf_rep1_peaks.narrowPeak peaks/tf_rep2_peaks.narrowPeak \
    --input-file-type narrowPeak \
    --output-file peaks/tf_idr_peaks.txt \
    --idr-threshold 0.05 \
    --plot

echo "IDR peaks: $(wc -l < peaks/tf_idr_peaks.txt)"

预期输出

输出	格式	描述
*_peaks.narrowPeak	BED6+4	窄峰，包含信号、p 值、q 值和峰顶偏移量
*_peaks.broadPeak	BED6+3	宽峰（使用 --broad 时）：chrom、start、end、signal、p-val、q-val
*_summits.bed	BED3+2	峰顶位置（1 bp）及得分；用于基序分析
*_treat_pileup.bdg	bedGraph	处理信号轨道（使用 -B 时）
*_control_lambda.bdg	bedGraph	对照/局部 lambda 轨道（使用 -B 时）
*_model.r	R 脚本	片段长度模型；运行 Rscript *_model.r 生成图表

故障排除

问题	原因	解决方案
调用出的峰很少	q 值过严格或读取深度较低	放宽至 -p 1e-3；检查测序深度（建议至少有 ≥10M 条比对后的读取）
峰过多（>100k）	阈值过宽松或没有输入对照	添加 --control input.bam；使用 -q 0.01；按 signalValue 过滤
峰调用失败并提示“no reads”	BAM 文件未排序或未建立索引	在运行 MACS3 前执行 samtools sort 和 samtools index
ATAC-seq 峰出现在细胞线粒体中	mtDNA 含量较高	过滤：`samtools view -h chip.bam
片段模型失败	读取数量过少或读取长度异常	添加 --nomodel --extsize 200 以跳过建模
bedGraph 输出过大	高覆盖度数据未归一化	添加 --SPMR，将信号归一化为每百万次读取
--broad 漏掉窄峰	信号实际呈尖锐峰形	检查 ChIP 靶标：TF 和 H3K4me3 需要使用窄峰模式
gsize 不匹配	组装版本对应的基因组大小错误	hg19/hg38 使用 hs，mm9/mm10 使用 mm；或提供准确的整数值

参考资料

• MACS3 GitHub：macs3-project/MACS — 源代码、文档和变更日志
• Zhang Y 等（2008）“基于模型的 ChIP-Seq 分析（MACS）” — Genome Biology 9:R137。DOI:10.1186/gb-2008-9-9-r137
• ENCODE ATAC-seq 流水线 — 使用 MACS3 的 ENCODE 标准化 ATAC-seq 工作流
• IDR 框架 — 用于获得可重复峰调用的不可重复发现率