Peaks 筛选过滤策略

生信菜鸟团

发布于 2025-04-09 10:26:08

2580

ChIP-Seq 实验中，不同的 callpeak 策略和过滤选择能有多少影响呢？我们这里探索一下。

这里选择的数据来源是这篇文章：

《Genome-Wide Mapping of Binding Sites Reveals Multiple Biological Functions of the Transcription Factor Cst6p in Saccharomyces cerevisia》

GSM1975223 Cst6p_glucose_rep1
GSM1975224 Cst6p_glucose_rep2
GSM1975225 Cst6p_ethanol_rep1
GSM1975226 Cst6p_ethanol_rep2

探究的是 Cst6p 转录因子在 glucose 和 ethanol 两个条件下的结合。为了比较两种条件下转录因子结合的不同，glucose 作为对照组，ethanol 作为实验组。

数据下载，质控，比对，这些通用步骤在很多教程都有讲，不再多少。现在我们得到了 bam 文件。怕不好找，这里贴一个我的吧！

ChIP-Seq 分析流程-上游

callpeak 策略

合并还是单独？

我们现在每一组都有两个重复，我们是直接合并同组 bam, 然后进行callpeak，还是单样本 callpeak，然后取交集。

合并样本 callpeak

先合并文件后 IGV 看一下

samtools merge -r bam/glucose.bam bam/SRR3033154.bam bam/SRR3033155.bam
samtools merge -r bam/ethanol.bam bam/SRR3033156.bam bam/SRR3033157.bam

samtools index bam/glucose.bam
samtools index bam/ethanol.bam

我们可以清楚地看到，在我们可视化的区域中，glucose 样本的覆盖范围更均匀，而 ethanol 影响的样本则显示出不同的结合事件。

callpeak

Macs2 支持输入多个实验样本和对照样本，大多时候对照样本我们使用的是input 或者 igg做空白对照。macs2 先将同组多个文件进行 merge 合并成一个文件，和我们上面的步骤是一样的（我验证了！结果是一样的）。

GLU1=bam/SRR3033154.bam
GLU2=bam/SRR3033155.bam

ETH1=bam/SRR3033156.bam
ETH2=bam/SRR3033157.bam

macs2 callpeak -t $ETH $ETH2 -c $GLU $GLU --gsize 1E7  --name ETH1 --outdir ETH/

$ wc -l *bed
139 ETH_summits.bed

可以看出我们得到了139个peaks。

看一下我们发现的第一个peak。

还不错，看起来像是一个真实的结果。narrowPeak 每一列意义如下：

单独样本callpeak

我们一对一进行 callpeak

GLU1=bam/SRR3033154.bam
GLU2=bam/SRR3033155.bam

ETH1=bam/SRR3033156.bam
ETH2=bam/SRR3033157.bam

macs2 callpeak -t $ETH1 -c $GLU1 --gsize 1E7  --name ETH1 --outdir ETH1/ 
macs2 callpeak -t $ETH2 -c $GLU2 --gsize 1E7  --name ETH2 --outdir ETH2/

$ wc -l ETH1_summits.bed
90 ETH1_summits.bed

$ wc -l ETH2_summits.bed 
87 ETH2_summits.bed

继续看第一个peak

可以看到单独样本 peak 的峰值位点都不太一样，也可以看出合并样本的包容性似乎是最强的。

评估重复性

单样本分开分析就可以评估样本的重复性。推荐两个工具吧。

bedtools jaccard

$ bedtools jaccard -a ETH1/*narrow* -b ETH2/*narrow*
intersection    union   jaccard n_intersections
10090   17339   0.581925        58

指标	含义
intersection	两重复峰区域重叠的碱基数（bp）
union	两重复峰区域合并后的总碱基数（去重）
jaccard	重复性指标 = intersection / union，范围[0,1]，值越高重复性越好
n_intersections	实际重叠的峰数量（非碱基数）

IDR

Peaks 排序

sort -k8,8nr ETH1/ETH1_peaks.narrowPeak > IDR/ETH1_sorted.narrowPeak
sort -k8,8nr ETH2/ETH2_peaks.narrowPeak > IDR/ETH2_sorted.narrowPeak

IDR 分析

cd IDR

idr --samples ETH1_sorted.narrowPeak ETH2_sorted.narrowPeak \
--input-file-type narrowPeak \
--rank p.value \
--output-file ETH \
--plot \
--log-output-file ETH.idr.log

详细解读看推文：评估 ChIP-Seq 重复性

我们有两个重复样本，我们先看 IDR 考虑了多少重叠的 peaks。

$ wc -l ETH
58 ETH

这两个重复样本，重叠区域有多少区域 IDR < 0.05 呢。

$ awk '{if($5 >= 540) print $0}' ETH | wc -l
36

那我们也可以选择这些 peaks 为我们关注的 peaks.

合并多个 narrow peak 文件

我们也可以将多个 narrow peak 文件，进行合并，先单文件分析，再合并整合也是没有问题的。

合并两个 narrowPeak 文件中的峰值

cat ETH1/ETH1_peaks.narrowPeak ETH2/ETH2_peaks.narrowPeak | sort -k1,1 -k2,2n > merge/merged_sorted.narrowPeak

将重叠或相邻的区域合并为单个连续的峰值区域

bedtools merge -i merged_sorted.narrowPeak > final_union_peaks.bed

查看结果

$ wc -l merged_sorted.narrowPeak 
177 merged_sorted.narrowPeak

可以看到走合并样本 callpeak 依旧是最宽容的。

过滤掉基因组中重复或低复杂度的区域

由于测序方案的偏倚性，ChIP-Seq 数据通常会受到大量假阳性的影响。过滤掉基因组中重复或低复杂度的区域，来自这些区域的数据几乎总是显示错误信号。这些区域的 peaks 应该被过滤掉。

peaks 注释，GO和KEGG富集分析，我们一般会使用峰顶文件，即仅包含峰顶坐标，我们可以筛选剔除重复或低复杂度的区域的峰顶文件。

拓展：可以先根据峰值先过滤一下

cat ETH/ETH_summits.bed | awk ' $5> 10 { print $0 }' > highscore.bed

sdust https://github.com/lh3/minimap 程序可用于创建低复杂度区域的区间文件

sdust refs/saccer3.fa > lowcomplexity.bed

只保留不属于低复杂度区域的高分值

bedtools intersect -v -a highscore.bed -b lowcomplexity.bed > highcomplexity.bed

查看结果

$ wc -l *bed
    28 highcomplexity.bed
    54 highscore.bed
 15603 lowcomplexity.bed
 15685 total

只保留基因起始上游的峰值（根据需求）

基因起始位点上游（如-2kb至+0.5kb）是核心启动子区域，富集转录因子结合位点（TFBS）、RNA聚合酶结合位点和表观调控元件（如CpG岛）。
调控机制：保留这些区域可筛选出可能直接参与基因转录起始调控的峰（如转录因子结合或组蛋白修饰H3K4me3）。

提取基因上游区域 (1kb）

bedtools flank -g refs/sc.fa.fai -l 1000 -r 0 -s -i refs/genes.gff > flanked.gff

取交集

bedtools intersect -u -a ETH/ETH_summits.bed -b flanked.gff > upstream-peaks.bed

观点

说了这么多探索性分析，那我们到底应该怎么分析呢？（个人观点）

方法	优势	局限性
合并BAM后Callpeak	提升低丰度信号检测能力	可能引入系统性误差
单独Callpeak+IDR	严格评估可重复性	仅支持两两比较
单独Callpeak+bedtools	支持多重复直接交集	忽略信号强度差异