02 Peak 注释与多样本比较
02 Peak 注释与多样本比较
拿到 peak 文件之后，第一个问题是"这些 peak 落在基因组的什么位置"。ChIPseeker 把 peak 注释到最近的基因、标注它在启动子 / 内含子 / 基因间等区域，一步出图。
本章用 ChIPseeker 自带的 AR（雄激素受体）ChIP-seq 数据演示：3 个剂量（0M / 1nM / 100nM）的 peak 文件，看剂量增加时 peak 数量、位置分布和关联基因如何变化。
核心流程
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
peak <- readPeakFile("peaks.narrowPeak")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
anno <- annotatePeak(peak, TxDb = txdb, tssRegion = c(-3000, 3000))
plotAnnoPie(anno)
plotDistToTSS(anno)
tssRegion = c(-3000, 3000) 定义"启动子"的范围：TSS 上下游 3kb 以内的 peak 算启动子区域。
真实示例
配套脚本 epi02_chipseeker_sci.R 在 ChIPseeker 内置的 AR ChIP-seq 数据上跑完整流程：
Rscript scripts/epigenomics/epi02_chipseeker_sci.R
每张图看什么
图 1：AR 100nM 的 peak 落在哪些基因组区域。启动子占比越高，说明这个 TF 越倾向于结合在基因起始位点附近。AR 是经典的启动子 + 增强子结合 TF，所以启动子和远端基因间区域都有不少。
图 2：三个剂量的 peak 到最近 TSS 的距离分布。剂量越高，靠近 TSS 的 peak 比例越大 —— 说明高剂量下 AR 更多地占据启动子。
图 3：三个剂量的基因组区域分布对比。堆叠条形图一眼看出"启动子占比随剂量增加"。
图 4：三个剂量的 peak 关联基因重叠。"三个剂量都有"的基因是 AR 的核心靶基因；"只在 100nM 出现"的是高剂量特异的。
图 5：AR 100nM peak 关联基因的 GO BP 富集。应该能看到雄激素响应、细胞增殖调控等通路。
图 6：每条染色体上的 peak 数量。分布大致和染色体大小成正比，但某些染色体可能因为基因密度高而偏多。
套到自己数据上
把 getSampleFiles() 换成自己的 .narrowPeak 文件路径即可。注意：
- 参考基因组版本要和 peak 文件一致（hg19 / hg38 / mm10）
- 对应的 TxDb 包：
TxDb.Hsapiens.UCSC.hg38.knownGene（hg38）或TxDb.Mmusculus.UCSC.mm10.knownGene（小鼠） - 如果是 broadPeak（组蛋白修饰），
annotatePeak的参数不用改，但解读时"启动子"的含义要注意
常见坑
坑 1：tssRegion 范围拍脑袋定
默认 tssRegion = c(-3000, 3000)，但组蛋白活化标记 H3K4me3 通常在 ±1kb 内最强，K27ac 可以扩到 ±5kb。范围拍错了会让"启动子" peak 占比虚高或虚低。根据自己实验类型调，TF 一般 ±2kb 够用。
坑 2：peak 注释到"最近基因"当成调控关系
peak 离哪个基因最近 ≠ peak 调控哪个基因。增强子可以跨数十 kb 调控其他基因，严格的调控关系需要 Hi-C / 4C / ABC model 等额外数据。把"最近基因"当作初筛候选可以，写在文章里要明确说"基于 proximity 的关联"。
坑 3：多样本对比时 peak 数量直接比
样本 A 在更深的测序下 call 出更多 peak，不代表 A 比 B"信号更强"。用 consensus peak set 上的 reads 数（DiffBind 思路）做基于强度的对比，光比较 peak 数量是误导。
坑 4：GO 富集背景错
ChIPseeker 的 enrichGO 默认背景是物种全部基因，但你的 ChIP/ATAC 只在某些基因上有信号，应该用"所有可被检测到 peak 的基因"做 universe。enrichGO(universe = all_peak_genes) 显式传背景。
坑 5：染色体名前缀不一致
peak 文件可能用 "chr1"，但有些 GTF 用 "1"。ChIPseeker 不会自动统一，要么用 seqlevelsStyle() 改一致，要么 manually seqlevels(peak) <- paste0("chr", seqlevels(peak))。否则 annotatePeak 会找不到任何基因。
下一步
接着深入：
- 03 DiffBind 差异结合分析 — 在 consensus peak 上做差异强度分析
- 05 Motif 富集与 HOMER — 知道 peak 在哪后下一步：哪些 TF motif 在这些 peak 富集
横向延伸：
- 04 Peak 可视化与多样本比较 — 注释结果做成可发表的图
- bulk RNA-seq 03 GO/KEGG/GSEA — peak 关联基因的富集分析思路完全一致