ChIP-seq 峰值调用分析¶

一句话概述¶

ChIP-seq（染色质免疫沉淀测序）用于检测蛋白质（如转录因子、组蛋白修饰）在基因组上的结合位点，MACS2 是最主流的峰调用工具。

核心知识点表格¶

知识点	说明
ChIP-seq	用抗体"钓"出与特定蛋白结合的 DNA，然后测序找到结合位点
Input 对照	未经免疫沉淀的 DNA，作为背景噪声参考
Narrow Peak	尖锐峰（转录因子，如 CTCF），通常 <500bp
Broad Peak	宽峰（组蛋白修饰，如 H3K27me3），可达数十 kb
MACS2/MACS3	最常用的峰调用工具，用泊松分布建模背景
FRiP	Fraction of Reads in Peaks，>1% 合格，>5% 良好
IDR	不可重复性发现率，用于多重复合并（ENCODE 标准）
Motif 分析	在峰区域中寻找富集的 DNA 序列模式（如转录因子结合位点）
HOMER	综合性 ChIP-seq 分析工具，擅长 motif 发现
ENCODE	DNA 元件百科全书项目，制定了 ChIP-seq 分析标准

白话解释原理¶

想象一下： 你想知道一个蛋白质（比如转录因子 p53）在基因组哪些位置"坐着"。

ChIP-seq 就像一个"钓鱼"实验： 1. 用甲醛把蛋白质和 DNA "焊"在一起（交联） 2. 把 DNA 打碎成小片段 3. 用 p53 的抗体当"鱼钩"，把带着 p53 的 DNA 片段"钓"出来 4. 把钓出来的 DNA 测序 5. 看看哪些基因组位置的测序信号特别高 → 那就是 p53 的结合位点

MACS2 做的事情：在背景噪声中找出信号显著高于周围的"山峰"（peak）。

各步骤详解¶

第一步：数据预处理¶

白话解释： 从原始 FASTQ 到比对好的 BAM 文件。

# 1. 质量控制
fastqc chip_sample.fastq.gz -o qc_results/  # 检查原始数据质量
fastqc input_control.fastq.gz -o qc_results/  # 检查input对照质量

# 2. 去接头和低质量碱基
trim_galore \
    --quality 20 \               # 质量阈值Q20
    --length 36 \                # 最短读长36bp
    --fastqc \                   # 自动运行FastQC
    chip_sample.fastq.gz         # 输入文件

# 3. 比对到参考基因组（Bowtie2）
bowtie2 \
    -x /data/genome/hg38_index \    # Bowtie2索引
    -U chip_sample_trimmed.fq.gz \  # 输入trimmed reads
    --very-sensitive \               # 高灵敏度模式
    --no-mixed \                     # 不允许不配对比对
    -p 8 \                           # 8个线程
    -S chip_sample.sam               # 输出SAM文件

# 4. 转换排序和去重
samtools view -bSq 20 chip_sample.sam | \   # 过滤低质量比对(MAPQ<20)
    samtools sort -@ 8 -o chip_sample.sorted.bam  # 排序

samtools index chip_sample.sorted.bam  # 建索引

# 去除PCR重复
picard MarkDuplicates \
    I=chip_sample.sorted.bam \         # 输入BAM
    O=chip_sample.dedup.bam \          # 输出去重BAM
    M=dup_metrics.txt \                # 重复统计
    REMOVE_DUPLICATES=true             # 移除重复reads

samtools index chip_sample.dedup.bam  # 重建索引

# 5. 去除黑名单区域
bedtools intersect \
    -v \                               # 排除模式
    -abam chip_sample.dedup.bam \      # 输入BAM
    -b hg38-blacklist.v2.bed \         # ENCODE黑名单
    > chip_sample.clean.bam            # 输出清洗后BAM

第二步：MACS2 峰调用¶

白话解释： MACS2 会扫描整个基因组，找出信号显著高于背景的区域（峰）。

# 安装 MACS2
pip install macs2  # 安装（Python包）
# 或 conda install -c bioconda macs2

# === 窄峰调用（转录因子，如 CTCF、p53） ===
macs2 callpeak \
    -t chip_sample.clean.bam \     # 处理组（ChIP样本）
    -c input_control.clean.bam \   # 对照组（Input）
    -f BAM \                        # 输入格式
    -g hs \                         # 基因组大小（hs=人类，mm=小鼠）
    -n my_chip \                    # 输出文件前缀
    --outdir peaks/ \               # 输出目录
    -q 0.05 \                       # FDR阈值0.05
    --keep-dup auto \               # 自动处理重复reads
    --bdg                           # 输出bedGraph信号轨迹

# 输出文件:
# my_chip_peaks.narrowPeak    — 峰坐标和统计信息
# my_chip_summits.bed         — 峰顶点坐标
# my_chip_treat_pileup.bdg    — 处理组信号
# my_chip_control_lambda.bdg  — 背景信号

# === 宽峰调用（组蛋白修饰，如 H3K27me3、H3K36me3） ===
macs2 callpeak \
    -t chip_h3k27me3.clean.bam \   # H3K27me3 ChIP样本
    -c input_control.clean.bam \   # Input对照
    -f BAM \                        # 格式
    -g hs \                         # 基因组大小
    -n h3k27me3 \                   # 前缀
    --broad \                       # 宽峰模式
    --broad-cutoff 0.1 \            # 宽峰FDR阈值
    --outdir peaks/                 # 输出目录

# === 无对照峰调用（不推荐但有时必须） ===
macs2 callpeak \
    -t chip_sample.clean.bam \     # 只有ChIP样本
    -f BAM \                        # 格式
    -g hs \                         # 基因组大小
    -n no_input \                   # 前缀
    --nolambda \                    # 不使用局部lambda模型
    -q 0.01                         # 用更严格的阈值补偿

第三步：质量控制¶

# 1. 计算 FRiP（Fraction of Reads in Peaks）
total_reads=$(samtools view -c chip_sample.clean.bam)  # 总reads数
reads_in_peaks=$(bedtools intersect \
    -a chip_sample.clean.bam \       # BAM文件
    -b peaks/my_chip_peaks.narrowPeak \  # 峰文件
    -u -bed | wc -l)                 # 落在峰内的reads数

frip=$(echo "scale=4; $reads_in_peaks / $total_reads" | bc)  # 计算FRiP
echo "FRiP = $frip"  # 打印FRiP值
# FRiP > 0.01 (1%): 合格
# FRiP > 0.05 (5%): 良好

# 2. IDR 分析（多重复实验一致性）
# 安装 IDR
pip install idr

# 对两个生物学重复分别调峰
macs2 callpeak -t rep1.bam -c input.bam -f BAM -g hs -n rep1 --outdir peaks/
macs2 callpeak -t rep2.bam -c input.bam -f BAM -g hs -n rep2 --outdir peaks/

# 运行 IDR
idr \
    --samples peaks/rep1_peaks.narrowPeak peaks/rep2_peaks.narrowPeak \  # 两个重复的峰
    --input-file-type narrowPeak \     # 输入格式
    --output-file idr_results.txt \    # 输出文件
    --plot                              # 生成IDR图
# IDR < 0.05 的峰被认为是可重复的

# 3. 使用 deeptools 生成信号热图
# 安装: pip install deeptools

# 生成 bigWig 信号轨迹
bamCoverage \
    -b chip_sample.clean.bam \       # 输入BAM
    -o chip_signal.bw \              # 输出bigWig
    --normalizeUsing RPKM \          # RPKM标准化
    --binSize 10                      # 10bp分辨率

# 在 TSS 周围绘制信号热图
computeMatrix reference-point \
    -S chip_signal.bw \               # 输入信号
    -R genes.bed \                     # 基因注释
    --referencePoint TSS \             # 以TSS为中心
    -a 3000 -b 3000 \                  # 上下游各3kb
    -o matrix_tss.gz                   # 输出矩阵

plotHeatmap \
    -m matrix_tss.gz \                 # 输入矩阵
    -o tss_heatmap.pdf \               # 输出图片
    --colorMap Blues                    # 蓝色系配色

第四步：Motif 分析（HOMER）¶

白话解释： 在峰区域中寻找反复出现的 DNA 序列模式，推测哪些转录因子可能在这些位置结合。

# 安装 HOMER
# conda install -c bioconda homer

# 1. 准备峰文件
# 将 narrowPeak 转为 HOMER 格式
awk 'OFS="\t" {print $1":"$2"-"$3, $1, $2, $3, "+"}' \
    peaks/my_chip_peaks.narrowPeak > peaks_homer.txt  # 转换格式

# 2. 运行 motif 发现
findMotifsGenome.pl \
    peaks/my_chip_summits.bed \    # 峰顶点文件（推荐用summit）
    hg38 \                          # 基因组版本
    motif_results/ \                # 输出目录
    -size 200 \                     # 峰顶点上下各100bp
    -mask                           # 屏蔽重复序列

# 输出关键文件:
# knownResults.html  — 已知motif富集结果
# homerResults.html  — 从头发现的motif
# knownResults.txt   — 文本格式结果

# 3. 峰注释（HOMER）
annotatePeaks.pl \
    peaks/my_chip_peaks.narrowPeak \  # 峰文件
    hg38 \                             # 基因组
    -go go_results/ \                  # GO富集分析
    -genomeOntology genome_ont/ \      # 基因组特征注释
    > annotated_peaks.txt              # 输出注释结果

# 注释结果包括: 每个峰最近的基因、距TSS距离、基因组特征（启动子/增强子/内含子等）

第五步：差异结合分析¶

# 使用 DiffBind 做差异结合分析（R语言）
# BiocManager::install("DiffBind")

library(DiffBind)  # 加载差异结合包

# 1. 创建样本信息表
samples <- data.frame(
    SampleID = c("WT_1", "WT_2", "KO_1", "KO_2"),  # 样本ID
    Condition = c("WT", "WT", "KO", "KO"),            # 条件
    Replicate = c(1, 2, 1, 2),                         # 重复
    bamReads = c("wt1.bam", "wt2.bam", "ko1.bam", "ko2.bam"),  # BAM文件
    ControlID = c("Input", "Input", "Input", "Input"),  # 对照
    bamControl = c("input.bam", "input.bam", "input.bam", "input.bam"),  # 对照BAM
    Peaks = c("wt1.narrowPeak", "wt2.narrowPeak",      # 峰文件
              "ko1.narrowPeak", "ko2.narrowPeak"),
    PeakCaller = rep("narrow", 4)                       # 峰调用器类型
)

# 2. 构建 DiffBind 对象
dba_obj <- dba(sampleSheet = samples)  # 创建对象
dba_obj <- dba.count(dba_obj)          # 计数reads
dba_obj <- dba.contrast(dba_obj,        # 设置对比
                        categories = DBA_CONDITION)
dba_obj <- dba.analyze(dba_obj)        # 差异分析

# 3. 提取结果
results <- dba.report(dba_obj, th = 0.05)  # FDR < 0.05
dba.plotMA(dba_obj)                         # MA图
dba.plotVolcano(dba_obj)                    # 火山图

常见报错与解决¶

报错	原因	解决方案
`Too few peaks called (0 or <100)`	数据质量差或参数不对	检查FRiP，放宽q值阈值(如0.1)，确认input是否正确
`macs2: command not found`	安装路径未配置	确认conda环境已激活，或检查PATH
`Redundant rate too high`	PCR重复过多	重新做去重步骤，检查文库复杂度
`Tag size is 0`	无法自动检测读长	手动指定 `--tsize 150`
`Unable to open BAM file`	BAM文件损坏或无索引	重新排序+建索引: `samtools sort && samtools index`
`HOMER: genome not installed`	缺少基因组数据	运行 `perl configureHomer.pl -install hg38`

速查表¶

┌──────────────────────────────────────────────────────────┐
│                ChIP-seq 分析速查                          │
├──────────────────────────────────────────────────────────┤
│ 流程: FASTQ → Trim → Bowtie2 → 去重 → MACS2 → Motif     │
│                                                          │
│ MACS2 关键参数:                                           │
│   -t: ChIP BAM    -c: Input BAM                         │
│   -g: 基因组大小 (hs/mm/ce/dm)                           │
│   -q: FDR阈值 (默认0.05)                                 │
│   --broad: 宽峰模式 (组蛋白修饰用)                        │
│   --nomodel --extsize 200: 跳过模型建立                  │
│                                                          │
│ 峰类型选择:                                               │
│   Narrow: TF (CTCF, p53, GATA) → 默认模式               │
│   Broad: 组蛋白 (H3K27me3, H3K36me3) → --broad          │
│   Mixed: H3K4me1 → --broad 或默认都可以                  │
│                                                          │
│ 质控标准:                                                 │
│   比对率 > 70%                                           │
│   FRiP > 1% (最低), > 5% (良好)                         │
│   峰数量: TF 通常 1k-50k, 组蛋白 10k-100k+              │
│   IDR < 0.05 (重复一致性)                                │
│                                                          │
│ Motif工具: HOMER / MEME-ChIP / DREME                    │
│ 差异分析: DiffBind / MAnorm                              │
│ 可视化: IGV / deeptools / pyGenomeTracks                 │
└──────────────────────────────────────────────────────────┘

面试高频问题¶

MACS2 是怎么调用峰的？它的统计原理是什么？ 答：MACS2 首先通过链偏移模型估算 DNA 片段大小（fragment size），然后将 reads 延伸到片段长度。接着用滑动窗口扫描基因组，在每个窗口中用泊松分布比较 ChIP 信号与局部 lambda 背景（取局部 1kb、5kb、10kb 和全基因组的最大值作为背景估计），计算 p 值，再通过 BH 方法控制 FDR。最后合并相邻显著窗口得到峰。
Narrow peak 和 Broad peak 在生物学上有什么区别？分析方法有何不同？ 答：Narrow peak 对应点状结合的蛋白（如转录因子），峰宽通常 <500bp，有明确的 summit；Broad peak 对应大范围修饰的组蛋白标记（如 H3K27me3 覆盖整个被沉默的基因），可达数十 kb。分析差异：MACS2 默认 narrow 模式；宽峰需用 --broad 参数，使用更宽松的合并策略；SICER/EPIC2 也是专门为宽峰设计的工具。
为什么 ChIP-seq 需要 Input 对照？没有 Input 能做吗？ 答：Input 对照的目的是校正基因组中的技术偏差，包括：染色质开放度不均匀导致某些区域更容易被超声打断、基因组重复区域产生假信号、GC 偏差等。没有 Input 也能调峰（--nolambda），但假阳性率会大幅升高，不推荐在正式分析中使用。ENCODE 标准强制要求有 Input 或 IgG 对照。