MACS3 Peak Calling — ChIP-seq 和 ATAC-seq 峰值检测标准工具¶

一句话说明¶

MACS3 是 MACS2 的继任版本，是 ChIP-seq 和 ATAC-seq 最广泛使用的 peak calling 工具，通过模拟泊松分布检测 DNA 蛋白结合位点或开放染色质区域，支持 Python 3.9+，MACS2 已停止维护。

安装与配置¶

# 创建专用环境（MACS3 需要 Python 3.9+）
conda create -n macs3 python=3.11 -y
conda activate macs3

# 安装 MACS3（最新稳定版 3.0+）
pip install macs3

# 验证安装
macs3 --version

# 查看帮助（各子命令功能）
macs3 --help
macs3 callpeak --help

# 安装 samtools（BAM 文件操作，经常配合使用）
conda install -c bioconda samtools -y

核心用法¶

# ── ChIP-seq Peak Calling（最基础用法） ──────────────
macs3 callpeak \
  -t ChIP.bam \                        # 处理组 BAM（target）
  -c Input.bam \                       # 对照组 BAM（control/input）
  -f BAM \                             # 输入文件格式
  -g hs \                              # 基因组大小（hs=人类，mm=小鼠）
  -n sample_name \                     # 输出文件名前缀
  --outdir peaks/ \                    # 输出目录
  2>&1 | tee macs3.log                 # 保存运行日志

# ── ATAC-seq Peak Calling（推荐参数） ────────────────
macs3 callpeak \
  -t ATAC.bam \                        # ATAC-seq 比对 BAM
  -f BAM \                             # 文件格式
  -n atac_sample \                     # 输出前缀
  --outdir atac_peaks/ \               # 输出目录
  -g hs \                              # 人类基因组
  --nomodel \                          # 不建核小体移位模型（ATAC 用此参数）
  --extsize 200 \                      # 片段延伸长度（ATAC 用 200bp）
  --shift -100 \                       # 偏移量（ATAC 用 -100，代表 cut site）
  --keep-dup all \                     # 保留所有 reads（单细胞 ATAC 不去重）
  -B \                                 # 输出 bedgraph 文件（用于可视化）
  --SPMR                               # 信号标准化为每百万 reads

参数详解¶

参数	说明	ChIP-seq	ATAC-seq
`-t`	处理组/样本 BAM 文件	实验组	ATAC BAM
`-c`	对照组 BAM（input）	input/IgG	可不加
`-g`	基因组大小（`hs`/`mm`/数字）	hs	hs
`--nomodel`	不建片段大小模型	窄峰不用	必须加
`--extsize`	片段延伸长度（nomodel 时用）	200	200
`--shift`	reads 偏移量	0	-100
`--broad`	宽峰模式（H3K27me3 等组蛋白修饰）	宽峰时用	不用
`--keep-dup`	PCR 重复处理（`all`/`1`/`auto`）	1 或 auto	all
`-q`	FDR 显著性阈值	0.05	0.05
`-p`	p 值阈值（不与 -q 同用）	0.01	0.01
`-B`	输出 bedgraph 覆盖度文件	可选	推荐

实战案例¶

# ── 实战1：ChIP-seq 窄峰（转录因子结合位点） ──────────
macs3 callpeak \
  -t CTCF_ChIP.bam \                   # CTCF ChIP-seq 实验组
  -c Input.bam \                       # Input 对照
  -f BAM \
  -g hs \
  -n CTCF_peaks \
  --outdir ChIP_results/ \
  -q 0.01 \                            # FDR < 1%（比默认 5% 更严格）
  2>&1 | tee ChIP_macs3.log

# ── 实战2：ChIP-seq 宽峰（H3K27me3 等组蛋白修饰） ───
macs3 callpeak \
  -t H3K27me3.bam \
  -c Input.bam \
  -f BAM \
  -g hs \
  -n H3K27me3_peaks \
  --outdir broad_peaks/ \
  --broad \                            # 开启宽峰模式
  --broad-cutoff 0.1 \                 # 宽峰合并阈值（宽峰比窄峰宽松）
  2>&1 | tee H3K27me3_macs3.log

# ── 实战3：单细胞 ATAC-seq 伪 bulk 分析 ──────────────
# 先合并同一 cluster 所有细胞的 BAM 文件（伪 bulk）
samtools merge \
  cluster0_merged.bam \                # 输出合并 BAM
  cell1.bam cell2.bam cell3.bam \      # 同一 cluster 的细胞 BAM 文件
  --threads 8

# 排序和索引
samtools sort -@ 8 -o cluster0_sorted.bam cluster0_merged.bam
samtools index cluster0_sorted.bam

# 对合并后的伪 bulk 做 peak calling
macs3 callpeak \
  -t cluster0_sorted.bam \
  -f BAM \
  -g hs \
  -n cluster0 \
  --outdir scATAC_peaks/ \
  --nomodel \                          # ATAC 不建模
  --extsize 200 \                      # 延伸 200bp
  --shift -100 \                       # cut site 偏移
  --keep-dup all \                     # 单细胞数据不去重复
  -q 0.05 \                            # FDR < 5%
  2>&1 | tee cluster0_macs3.log

# ── Python 分析 MACS3 输出结果 ───────────────────────
import pandas as pd                    # 数据处理
import numpy as np
import matplotlib.pyplot as plt

# 读取 MACS3 输出的 narrowPeak 文件
# 列名：chr, start, end, name, score, strand, fc, -log10(pvalue),
#       -log10(qvalue), summit_offset
narrowpeak_cols = [
    'chrom', 'start', 'end', 'name', 'score',
    'strand', 'fold_change', 'neg_log10_pvalue',
    'neg_log10_qvalue', 'summit_offset'
]

peaks_df = pd.read_csv(
    'peaks/sample_name_peaks.narrowPeak',
    sep='\t',
    header=None,
    names=narrowpeak_cols
)

# 统计 peak 信息
print(f"总 peak 数量：{len(peaks_df)}")
print(f"Peak 宽度统计：\n{peaks_df.eval('end - start').describe()}")
print(f"\n各染色体 peak 数量（top 10）：")
print(peaks_df['chrom'].value_counts().head(10))

# 过滤掉线粒体和低质量染色体的 peak
peaks_filtered = peaks_df[
    ~peaks_df['chrom'].isin(['chrM', 'chrUn']) &  # 去除线粒体
    ~peaks_df['chrom'].str.contains('_')           # 去除 scaffold 染色体
].copy()
print(f"\n过滤后 peak 数量：{len(peaks_filtered)}")

# 按 FDR 过滤（-log10 qvalue > 2 对应 FDR < 1%）
peaks_highconf = peaks_filtered[peaks_filtered['neg_log10_qvalue'] > 2]
print(f"高置信度 peak（FDR<1%）：{len(peaks_highconf)}")

# 可视化 peak 信号分布
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.hist(peaks_filtered['fold_change'], bins=50, color='steelblue', edgecolor='white')
plt.xlabel('Fold Change over Input')
plt.ylabel('Count')
plt.title('Peak 信号强度分布')

plt.subplot(1, 2, 2)
plt.hist(peaks_filtered.eval('end - start'), bins=100,
         color='coral', edgecolor='white')
plt.xlabel('Peak Width (bp)')
plt.ylabel('Count')
plt.title('Peak 宽度分布')
plt.tight_layout()
plt.savefig('peak_statistics.pdf', dpi=150)

# 保存过滤后的 peak 文件
peaks_filtered.to_csv(
    'filtered_peaks.bed',
    sep='\t',
    header=False,
    index=False
)
print("结果保存完成！")

常见报错与解决¶

报错	原因	解决方法
`MACS3 command not found`	未安装或环境未激活	`conda activate macs3`
`No peaks found`	数据质量差或参数错误	降低 `-q 0.1` 或检查比对率
`BAM not sorted`	BAM 文件未排序	`samtools sort` 先排序
`Too many duplicates`	高 PCR 重复率	加 `--keep-dup 1` 去重或检查文库
宽峰结果 peak 太少	阈值太严	`--broad-cutoff 0.15` 放宽

速查表¶

# ChIP-seq 窄峰（TF）
macs3 callpeak -t chip.bam -c input.bam -g hs -n prefix -q 0.05

# ATAC-seq peak calling
macs3 callpeak -t atac.bam -g hs -n prefix --nomodel --extsize 200 --shift -100 --keep-dup all

# 宽峰（组蛋白修饰）
macs3 callpeak -t chip.bam -c input.bam -g hs -n prefix --broad

# 输出文件
# prefix_peaks.narrowPeak  → peak 坐标（窄峰）
# prefix_peaks.broadPeak   → peak 坐标（宽峰）
# prefix_summits.bed       → 每个 peak 的最高点（summit）
# prefix_treat_pileup.bdg  → 覆盖度 bedgraph（可转 bigwig 用于 IGV）