8 比较分析
- ChIP-seq 数据的比较分析对于在重复实验之间或在不同的生理条件下比较感兴趣的蛋白质的结合是必不可少的。此章节将介绍通过
Peaks之间取交集来快速比较和如何使用基于Peakreads 密度的定量的方法来研究差异结合。
8.1 Peak 区域的交集
-
比较样本间
Peak的一种简单的方法是简单的取交集。这导致了在两个样品中是否发现峰的binary view，提供了对不同样品之间的Peak区域的相似性的粗略估计。然而，这种方法本质上低估了相似性，因为Peak calling依赖于应用于区域排序列表上的置信阈值（ 见章节 6）。因此，一个样品中超过阈值的Peak在第二个样品中可能刚好低于阈值，即使它显示出相当的富集（ref1）。另一种选择是将来自第一个样本的高置信度Peaks与第二个样本中用较低置信阈值识别的所有Peak取交集（ref1/ ref2）。虽然在较小的程度上，Peak的简单取交集也低估了样品之间的差异，因为即使两个样品中的一个Peak高于阈值，它仍然可以显示出非常不同的富集。 -
下面的代码显示如何计算所有成对的交集。
-
shell计算交集百分比
for sample1 in NRF1_CHIP_WT_1 NRF1_CHIP_WT_2 NRF1_CHIP_TKO_1 NRF1_CHIP_TKO_2 do L=$sample1"\t"$(cat peaks/${sample1}_peaks_peakzilla.bed | wc -l) for sample2 in NRF1_CHIP_WT_1 NRF1_CHIP_WT_2 NRF1_CHIP_TKO_1 NRF1_CHIP_TKO_2 do L=$L"\t"$(intersectBed -u -a peaks/${sample1}_peaks_peakzilla.bed -b peaks/${sample2}_peaks_peakzilla.bed | wc -l) done echo -e $L done > change/peak_overlap_table.txt # Table with: sample / total_number_of_peaks / overlap_for_each_sample cat change/peak_overlap_table.txt | column -t # Transformation into percentages cat changes/peak_overlap_table.txt | awk '{L=$1"\t"$2;for(i=3;i<=NF;i++){L=L"\t"$i*100/$2};print L}' > changes/peak_overlap_percent_table.txt cat changes/peak_overlap_percent_table.txt | column -t-
R中对样品聚类以及通过热图形式展示。# Clustering of samples and representation in a heatmap # Open R and load library R library(NMF) # load table x = read.table("changes/peak_overlap_percent_table.txt") data = x[, 3:6] rownames(data) = colnames(data) = x[, 1] # Generate a heatmap that clusters the samples pdf("changes/peak_overlap_percent_heatmap.pdf") aheatmap(data, breaks = 50, annRow = data.frame(condition = c("WT", "WT", "TKO", "TKO"), replicate = as.factor(c(1, 2, 1, 2))), annColors = list('condition' = c("blue", "red"), 'replicate' = c("balck", "grey")) ) dev.off() # Cloes R q()
-
-
重复样本数不同，使用默认
Peak阈值标识的Peak数量不同（ 7167 Vs 10232 ）。因此，分析是不对称的：WT_1中 98% 的Peak与 WT_2 重合，但 WT_2 中只有 67% 的Peak与 WT_1 重合。然而，如果考虑到附加的富集区，WT_2 中的附加Peak可能已经很好地存在于 WT_1 中。当样品未饱和时，更深的测序也会增加重复交集数目。Peak区域长度的差异会进一步扭曲结果，特别是如果接受任何重叠，比如 1bp 。注：由于这些偏差，通常在 Venn 图中显示的峰重叠可能具有误导性，因为不重叠的峰不应被解释为特定于样品的峰。
-
尽管存在这些限制，我们可以在本示例中得出结论，WT 和 TKO 细胞中 NRF1 的两个重复实验似乎共享它们的大部分
Peak区域。正如预期的那样，我们观察到 WT 和 TKO 之间重叠的峰值较少。然而，仍然有很大的重叠（ WT 76% 或者 TKO 45% )，这表明许多峰值是在不同条件之间共享的。 -
所有成对比较的结果热图允许我们检查哪些样本比其他样本更相似（ 图 8.1A ）。正如预期的那样，样本按条件进行聚类。此外，WT
Peak似乎更经常与 TKOPeak共享，而反之亦然，尽管这可能是由于 TKO 样本中的较高峰数而造成的伪像。
8.2 IRREPRODUCIBLE DISCOVERY RATE （ IDR ）
-
如上所述，ChIP-seq 分析中的每一个比较都强烈依赖于在
Peak calling步骤中选择的阈值。因此，ENCODE开发了不可复制的发现率（IDR），作为一种基于它们在重复之间的可重复性来识别真正Peak的标准（ref3）。其基本思想是，使用线性阈值生成的Peak list将包含真正的Peak和noise。当列表中的Peak被排序时，例如在它们的富集倍数或显著性上，对于真正的Peak，这些排序将很好地相关，而noise将显示不相关。IDR 使用统计方法来找到曲线中的点，在该点上，重复rank之间的关联的heterogeneity急剧增加。注意：此方法也可用于来自不同
Peak callers生成的相同样本的Peak list。在没有重复的情况下，这种方法可以帮助在单个样品中定义可靠的峰。注意：IDR 不适用于组蛋白修饰数据，因为宽峰的交集不明确。
8.2.1 Peak calling for IDR
-
只要评分不会产生太多的关联，从而导致排名不明确。IDR 可以与任何对
Peak进行排名的Peak caller的输出一起使用。在这里，我们展示了如何将 IDR 与peakzilla峰值一起使用。如前所述，为了使 IDR 工作，峰值列表必须同时包含真正的峰值和噪声，因此必须放宽峰值调用参数。 -
当比较重复时，需要一个统一的
Peak集合作为参考集合。这既可以由 IDR 软件从单独识别的Peak区域创建，也可以由用户提供。对于后者，在Peak calling之前合并重复应该有助于识别所有可能的Peak区域，来自合并样本的Peak list通常用于 IDR 分析。# Call peaks on each replicate with a relaxed threshold (-s 0.1 instead of default 1) peakzilla.py -s 0.1 <(bamToBed -i reads/NRF1_CHIP_WT_1.bam) <(bamToBed -i reads/NRF1_INPUT_WT.bam) -l peaks/NRF1_CHIP_WT_1_all_peakzilla_report.txt |\ awk -v OFS="\t" '(NR>1&&$1!="chrM"){print $1,$2-1,$3,$5,$6,"+"}' |\ intersectBed -v -a stdin -n mm10_black.bed.gz >peaks/NRF1_CHIP_WT_1_all_peaks_peakzilla.bed peakzilla.py -s 0.1 <(bamToBed -i reads/NRF1_CHIP_WT_2.bam) <(<(bamToBed -i reads/NRF1_INPUT_WT.bam) -l peaks/NRF1_CHIP_WT_2_all_peakzilla_report.txt |\ awk -v OFS="\t" '(NR>1&&$1!="chrM"){print $1,$2-1,$3,$5,$6,"+"}' |\ intersectBed -v -a stdin -b mm10_blacklist.bed.gz >peaks/NRF1_CHIP_WT_2_all_peaks_peakzilla.bed # Pool BAM files and call peaks on merged file with a relaxed threshold. samtools merge reads/NRF1_CHIP_WT.bam reads/NRF1_CHIP_WT_1.bam reads/NRF1_CHIP_WT_2.bam peakzilla.py -s 0.1 <(bamToBed -i reads/NRF1_CHIP_WT.bam) <(bamToBed -i reads/NRF1_INPUT_WT.bam) -l peaks/NRF1_CHIP_WT_1_all_peakzilla_report.txt |\ awk -v OFS="\t" '(NR>1&&$1!="chrM"){print $1,$2-1,$3,$5,$6,"+"}' |\ intersectBed -v -a stdin -b mm10_blacklist.bed.gz >peaks/NRF1_CHIP_WT_all_peaks_peakzilla.bed
8.2.2 计算 IDR
-
下一步是运行 IDR 软件以鉴定可重现的
Peak。对于生物学重复，通常的阈值是 0.05，即最后list中高达 5% 不能被重现出来。技术重复应使用较低的阈值，因为它们的总体可变性较低。idr --samples peaks/NRF1_CHIP_WT_1_all_peaks_peakzilla.bed peaks/NRF1_CHIP_WT_2_all_peaks_peakzilla.bed --peak-list peaks/NRF1_CHIP_WT_all_peaks_peakzilla.bed --input-file-type bed --rank score --idr-threshold 0.05 --output-file peaks/NRF1_CHIP_WT_idr --plot -
IDR 的概念在很大程度上依赖于有两个好的重复。如果其中一个重复显示质量较差，例如如果
IP不是有效的，IDR 将只记录非常少的可重现峰。对于这种情况，ENCODE 开发了一种拯救策略，方法是将两个重复汇集在一起，然后随机地将 reads 分成两个伪重复。这些并不代表真正的生物或实验变异，但用于对来自DNA片段群体的reads 采样中的随机噪声进行建模。对于伪复制的IDR比较，为了降低了噪声，建议使用较低的阈值 （ 0.0025 ）。
8.2.3 reads 密度的比较
-
在 reads 密度的水平上，评估样本之间的总体相似性的最简单的方法是全局地关联它们的 reads 密度。这可以通过计算整个基因组的标准化 reads 密度上的
Pearson相关系数（PCC）来实现，无论是对于每个单独的碱基对（ref3），还是在沿着基因组的滑动窗口中。然而，由于Peak区域仅代表基因组的一小部分， reads 密度的高相关性将主要反映一致的背景信号，如 ChIP 和 Input 样本之间的高相关性所显示的（见表1）。表1 reads 密度的
Pearson相关系数（ PCC ）Sample1 Sample2 PCC NRF1_CHIP_WT_1 NRF1_CHIP_WT_2 0.98 NRF1_CHIP_WT_1 NRF1_CHIP_TKO_1 0.97 NRF1_CHIP_WT-1 NRF1_INPUT_WT 0.96 沿着基因组的每个碱基对的 reads 密度的PCC（ 不包括两个样本中具有零 reads 的位置 ）。
-
为了具体比较
Peak区域中的 reads 密度，可以仅在至少一个样本中包含Peak的区域内计算PCC值。还可以使用散点图在视觉上比较每个区域的标准化平均 reads 密度。这代表了Peak区域中信号的更定量比较，而不是 8.1节 中解释的重叠Peak区域的二元方法。
8.3.1 合并 Peak 区间
- 为了聚焦于包含高信号并且可能在不同样本之间存在差异的基因组窗口，我们在
Peak区域内执行 reads 密度的比较。为此，我们将来自所有实验的Peak区域合并（ 可以针对所有可能的成对比较单独执行 ）。尽管这也依赖于Peak caller阈值，但它允许对所有Peak进行定量比较，包括仅存在于一个样本中的那些。
# 连接所有样本的 Peak，然后通过染色体和起始坐标进行排序
cat peaks/NRF1_CHIP_WT_1_peaks_peakzilla.bed peaks/NRF1_CHIP_WT_2_peaks_peakzilla.bed peaks/NRF1_CHIP_TKO_1_peaks_peakzilla.bed peaks/NRF1_CHIP_TKO_2_peaks_peakzilla.bed |\
sort -k1,1 -k2,2n > changes/NRF_all_regions_tmp.txt
# 合并所有样本的 Peak 区间
mergeBed -i change/NRF1_all_regions_tmp.txt > changes/NRF1_all_regions.txt
# 删除中间文件
rm changes/NRF1_all_regions_tmp.txt
# 一行命令解决
cat peaks/NRF1_CHIP_WT_1_peaks_peakzilla.bed peaks/NRF1_CHIP_WT_2_peaks_peakzilla.bed peaks/NRF1_CHIP_TKO_1_peaks_peakzilla.bed peaks/NRF1_CHIP_TKO_2_peaks_peakzilla.bed |\
sort -k1,1 -k2,2n |\
mergeBed -i stdin > changes/NRF1_all_regions.txt
8.3.2 计算每个样本的 reads 数
-
我们将合并区域与原始 reads 重叠，以计算每个样本中落入其中的reads 数。这些原始 reads 数将直接用作下一节中识别差异
Peak区域的输入。# 对于一个样品 intersectBed -c -sorted -a changes/NRF1_all_regions.txt -b reads/NRF1_CHIP_WT_1.bam > changes/NRF1_all_regions_count1.txt # 对于所有样本 intersectBed -c -sorted -a changes/NRF1_all_regions.txt -b reads/NRF1_CHIP_WT_1.bam |\ intersectBed -c -sorted -a stdin -b reads/NRF1_CHIP_WT_2.bam |\ intersectBed -c -sorted -a stdin -b reads/NRF1_CHIP_TKO_1.bam |\ intersectBed -c -sorted -a stdin -b reads/NRF1_CHIP_TKO_2.bam > changes/NRF1_all_regions_count.txt cat changes/NRF1_all_regions_count.txt |\ head | column -t
8.3.3 标准化 reads 数
-
由于 Peak 区域具有不同的大小，并且样本包含不同数量的比对 reads ，因此将每个 reads 计数标准化为该区域的大小和样本中比对 reads 的总数，以获得
RPKM值（ reads per kilobase per million）。# 将 reads 数标准化 （ RPKM ） # RPKM = reads * （1000000 / 总的 reads 数）* （ 1000 / Peak 区域大小） TOT1=$(samtools idxstats reads/NRF1_CHIP_WT_1.bam | awk '($1!˜"*"){t=t+$3}END{print t}') TOT2=$(samtools idxstats reads/NRF1_CHIP_WT_2.bam | awk '($1!˜"*"){t=t+$3}END{print t}') TOT3=$(samtools idxstats reads/NRF1_CHIP_TKO_1.bam | awk '($1!˜"*"){t=t+$3}END{print t}') TOT4=$(samtools idxstats reads/NRF1_CHIP_TKO_2.bam | awk '($1!˜"*"){t=t+$3}END{print t}') cat changes/NRF1_all_regions_count.txt | awk -v OFS="\t" -v TOT1=$TOT1 -v TOT2=$TOT2 -v TOT3=$TOT3 -v TOT4=$TOT4 '{size=$3-$2; print $1,$2,$3, $4*(1000000/TOT1)*(1000/size), $5*(1000000/TOT2)*(1000/size), $6*(1000000/TOT3)*(1000/size), $7*(1000000/TOT4)*(1000/size)}’ > changes/NRF1_all_regions_rpkm.txt cat changes/NRF1_all_regions_rpkm.txt |\ head | column -t# 使用循环计算 reads 数和产出 RPKM 值 cat changes/NRF1_all_regions.txt > changes/NRF1_all_regions_count.txt cat changes/NRF1_all_regions.txt > changes/NRF1_all_regions_rpkm.txt for sample in NRF1_CHIP_WT_1 NRF1_CHIP_WT_2 NRF1_CHIP_TKO_1 NRF1_CHIP_TKO_2 do TOT=$(samtools idxstats reads/${sample}.bam |\ awk '($1!˜"*"){t=t+$3}END{print t}') intersectBed -c -sorted \ -a changes/NRF1_all_regions_count.txt\ -b reads/${sample}.bam > tmp mv tmp changes/NRF1_all_regions_count.txt paste changes/NRF1_all_regions_rpkm.txt\ <(cat changes/NRF1_all_regions_count.txt |\ awk -v OFS="\t" -v TOT=$TOT ’{size=$3-$2; print $NF*(1000000/TOT)*(1000/size)}’) > tmp mv tmp changes/NRF1_all_regions_rpkm.txt done cat changes/NRF1_all_regions_count.txt |\ head | column -t cat changes/NRF1_all_regions_rpkm.txt |\ head | column -t
8.3.4 reads 数的比较
-
我们现在可以可视化散点图中跨样本的
Peak区域中的标准化 reads 数，并计算相关的 PCC 值。散点图为探索数据和得出结论提供了一种不带偏见的方式。为了更好地可视化数据的分布，RPKM 值以 log2 标准化。为了确认 WT 和 TKO 细胞之间 NRF1 结合的变化是可重复的，我们比较了重复之间的 TKO 和 WT 变化倍数。为此，我们向所有数据点添加一个伪计数（ 这里是0.1 ），以避免被0除，并以 log2 标准化，以获得以 0 为中心的正态分布。# 打开 R R # 加载 RPKM 矩阵 data = read.table("changes/NRF1_all_regions_rpkm.txt") colnames(data)=c("chr", "start", "end", "NRF1_CHIP_WT_1", "NRF1_CHIP_WT_2", "NRF1_CHIP_TKO_1", "NRF1_CHIP_TKO_2") # 散点图两两比较可视化 pdf("changes/NRF1_all_regions_rpkm_cor_scatter.pdf") par(bg="white",mfrow=c(2,2)) smoothScatter(log2(data$NRF1_CHIP_WT_1), log2(data$NRF1_CHIP_WT_2), xlim=c(-4,10), ylim=c(-4,10), xlab="NRF1_CHIP_WT_1(log2)", ylab="NRF1_CHIP_WT_2(log2)", main=paste("PCC=", signif(cor(data$NRF1_CHIP_WT_1,data$NRF1_CHIP_WT_2),2))) smoothScatter(log2(data$NRF1_CHIP_TKO_1), log2(data$NRF1_CHIP_TKO_2), xlim=c(-4,10), ylim=c(-4,10), xlab="NRF1_CHIP_TKO_1(log2)", ylab="NRF1_CHIP_TKO_2(log2)", main=paste("PCC=", signif(cor(data$NRF1_CHIP_TKO_1,data$NRF1_CHIP_TKO_2),2))) smoothScatter(log2(data$NRF1_CHIP_WT_1), log2(data$NRF1_CHIP_TKO_1), xlim=c(-4,10), ylim=c(-4,10), xlab="NRF1_CHIP_WT_1(log2)", ylab="NRF1_CHIP_TKO_1(log2)", main=paste("PCC=", signif(cor(data$NRF1_CHIP_WT_1,data$NRF1_CHIP_TKO_1),2))) smoothScatter(log2(data$NRF1_CHIP_WT_2), log2(data$NRF1_CHIP_TKO_2), xlim=c(-4,10), ylim=c(-4,10), xlab="NRF1_CHIP_WT_2(log2)", ylab="NRF1_CHIP_TKO_2(log2)", main=paste("PCC=", signif(cor(data$NRF1_CHIP_WT_2,data$NRF1_CHIP_TKO_2),2))) dev.off()# Compare the RPKM changes between WT and TKO samples across replicated samples pdf("changes/NRF1_all_regions_rpkm_cor_delta_scatter.pdf") par(bg="white") smoothScatter( log2((data$NRF1_CHIP_TKO_1+0.1)/(data$NRF1_CHIP_WT_1+0.1)), log2((data$NRF1_CHIP_TKO_2+0.1)/(data$NRF1_CHIP_WT_2+0.1)), xlim=c(-6,10), ylim=c(-6,10), xlab="NRF1_CHIP_TKO_1 / NRF1_CHIP_WT_1(log2)", ylab="NRF1_CHIP_TKO_2 / NRF1_CHIP_WT_2(log2)", main=paste("PCC=", signif( cor( (data$NRF1_CHIP_TKO_1+0.1)/(data$NRF1_CHIP_WT_1+0.1),(data$NRF1_CHIP_TKO_2+0.1)/(data$NRF1_CHIP_WT_2+0.1)), 2 ) ) ) dev.off() # 关闭 R q() -
我们观察到，WT 和 TKO 的重复样本的 reads 密度沿对角线排列，并且相关性很好（ PCC > 0.9 ）（ 图8.1 C ）。当比较 WT 和 TKO 样本时，reads 密度仍然相关，但小于重复之间的相关性（ PCC = 0.74 或 0.81 ）（图8.1D）。
-
此外，我们发现在 WT 样品中识别的所有峰在两种条件下都显示出相似的 reads 密度，因为 WT 样品中具有高 reads 密度的所有区域在 TKO 样品中也具有高 reads 密度并沿对角线排列（ 图8.1 D ）。相反，在 TKO 样本中识别的许多峰具有低 reads 数或不存在于 WT 样本中，通过图左上部分的数据点群体可视化。
注意：在曲线图的左下部分有很少或没有具有低 reads 数密度的点，这是由于我们只选择了在至少一个样本中被称为
Peak的区域，因此显示了 reads 计数的最小富集。在图 8.1 C 中，对于 TKO 重复样品，左下角的点表示在 WT 样品中识别但在 TKO 样品中具有背景 reads 密度的峰。如果从成对比较中合并Peak区域，则它们不会出现。 -
由于两个重复都显示了在 TKO 样本中获得的
Peak，因此检查重复 1和重复 2中的这些Peak是否相同是很有趣的。这是通过比较delta-delta曲线图中的变化倍数来确认的，这表明在 WT 和 TKO 之间观察到的 reads 密度变化在两个重复之间高度一致（ PCC = 0.67 ）（图 8.1 E）。
图 8.1 Peak 区间比较
（ A ）WT 和 TKO 两样本的相关性热图
（ B ）
Peakrank heterogeneity （ -log10IDR ）表示为peak秩的函数。Boxplots 显示每个包含 1000 个Peak的 bin 的总结（ C ）比较 TKO 重复样品之间
Peakreads 密度的散点图（ D ）比较 TKO 样品与 WT 样本之间
Peakreads 密度的散点图（ E ）两个重复的 TKO 和 WT 之间
Peak区域的 reads 密度的散点图比较比率image
8.4 差异结合分析
- 一旦散点图确认样品间存在差异
Peak，统计方法可以定义共享Peak或差异Peak的组，用于进一步分析。差异结合分析主要有两种类型的工具（ref4）。第一种类型采用基于 reads 计数数据的定量方法来比较一种条件下的结合强度与另一种条件下的结合强度。第二种类型使用隐马尔可夫模型将基因组分割成丢失、不变或获得的区域。然而，这些工具不允许在这三种截然不同的状态之外进行定量描述。在这里，我们介绍了使用DESeq2和DiffBind进行定量分析的典型分析流程，DiffBind为 ChIP-seq 分析提供了专门围绕DESeq2的封装。 - 在以下部分中，如果显示条件之间差异结合的
Peak区域分别在 WT或 TKO 细胞中显示更多的 NRF1 结合，则它们被称为 “WT-specific ” 或 “ TKO-specific ”。相反，显示条件之间的结合（ 在任一方向上 ）变化小于2倍的Peak被称为 “shared Peaks” 。
8.4.1 DESeq2
-
具有差异富集的
Peak区域的鉴定在概念上类似于差异表达基因的鉴定，因为两者都依赖于 reads 数的比较。这使我们能够采用最初为RNA-seq数据分析而设计的成熟的统计方法，例如R/Bioconductor软件包DESeq2（ref5）和edgeR（ref6）。DESeq2 使用基于负二项的广义线性模型来检验零假设，即两个条件之间 reads 数的log2FC等于零。它可以分解为四个主要步骤（ 包含在 DESeq() 函数中）：-
Size factors的估计，即将原始 reads 数标准化为Peak区域中的 reads 总数（代表测序深度 ）。Size factors被计算为每个Peak区域的 reads 数的几何平均值的中位数。 - 离散度的估计，即重复之间的变异性。由于少量的重复（ 通常是 2-3 个）会导致较大的方差，DESeq2 使用具有相似 reads 数的区域来更好地估计离散度（ 使用经验贝叶斯收缩方法 ）。
- 变化倍数的估计。由于小数字的比率本质上是
noisy，DESeq2 降低了低 reads 数的变化倍数估计（ 使用经验贝叶斯收缩方法 ）。因此，这些所谓的适当的变化倍数不同于根据标准化 reads 数数据直接计算比率。 - 在两个条件之间 reads 数的 log2 变化倍数为零的零假设下的差异富集检验。P 值通过
Wald 检验计算，并通过Benjamini-Hochberg进行校正。附加过滤使用所有区域的标准化 reads 数的平均值，以去除给定 FDR 阈值 以下的校正 p 值（ 默认 0.1；引入 NA ）。
library(DESeq2) # 加载章节 8.3.2 中的 reads 数矩阵 x = read.table("changes/NRF1_all_regions_count.txt") colnames(x) = c("chr","start","end", "NRF1_CHIP_WT_1", "NRF1_CHIP_WT_2", "NRF1_CHIP_TKO_1", "NRF1_CHIP_TKO_2") d = [, 4:7] # conditions colData = data.frame(condition = c("WT","WT","TKO","TKO")) rownames(colData) = colnames(d) # DESeq2 matrix dds = DESeqDataSetFromMatrix(countData = d, colData=colData, design=˜condition) # 过滤没有 reads 的区域（ 或者 reads 较低的 ） dds = dds[rowSums(counts(dds))>1,] # 转换数据进行方差分析 # #当样本数 < 30 时使用 regularised-logarithm transformation (RLD) rld=rlog(dds,blind=F) # 否则使用 variance-stabilising transformation (VST) vsd=vst(dds,blind=F) # 主成分分析（ PCA ） pdf("changes/NRF1_all_regions_count_pca.pdf") plotPCA(rld,intgroup=c("condition")) dev.off() # 差异分析（ 没有转换的数据 ） dds=DESeq(dds) # TKO versus WT res = results(dds, contrast = c("condition","TKO","WT")) # 火山图绘制 （ log2FC vs adjusted p-value ） pdf("changes/NRF1_all_regions_count_TKO_vs_WT_volcano.pdf") par(bg="white") plot(res$log2FoldChange, -log10(res$padj), xlim=c(-10,10), ylim=c(0,83), xlab="TKO / WT (log2)", ylab="P-value(-log10)", pch=20) lines(c(-10,10),c(10,10)) dev.off() # Table chr / start / end / NRF1_CHIP_WT_1 / NRF1_CHIP_WT_2 / NRF1_CHIP_TKO_1 / NRF1_CHIP_TKO_2 / log2FoldChange / padj y=cbind(x[,1:3],counts(dds,normalized=T), as.matrix(res)[,c(2,6)]) write.table(y, "changes/NRF1_all_regions_count_TKO_vs_WT_table.txt", quote=F, sep="\t",row.names=F,col.names=T) q()# 选择差异基因 # 由于校正后的 p 值与 FC 值成比例，因此我们只需要一个 p 值阈值 # 比如 10^-10 cat changes/NRF1_all_regions_count_TKO_vs_WT_table.txt | head | column -t # TKO-specific peaks: 2212 # No header, p-value (adjusted) <= 1e-10, log2 fold change >0 ( TKO > WT ) cat changes/NRF1_all_regions_count_TKO_vs_WT_table.txt | awk '(NR>1&&$9<=1e-10&&$8>0)' > changes/NRF1_all_regions_count_TKO_spec_table.txt # WT-specific peaks: 13 # No header, adjusted p-value <= 1e-10, log2 fold change < 0 ( WT > TKO ) cat changes/NRF1_all_regions_count_TKO_vs_WT_table.txt | awk ’(NR>1 && $9<=1e-10 && $8<0)’ > changes/NRF1_all_regions_count_WT_spec_table.txt # Shared peaks: 7068 # No header, absolute log2 fold change < 1 (less than 2 fold change in either direction) cat changes/NRF1_all_regions_count_TKO_vs_WT_table.txt | awk ’(NR>1 && $8<1 && $8>-1)’ > changes/NRF1_all_regions_count_shared_table.txt -
-
DESeq2 中的一个常见的
diagnostic plot是主成分分析（ PCA ），它允许我们检查多个样本之间的总体相似性（ 默认情况下，使用样本间方差最高的500个Peak区域 ）。在 NRF1 ChIP-seq 数据的情况下，这确认了重复聚集在一起，并且 99% 的方差是由 WT 和 TKO 条件之间的变化来解释的（图 8.2A）。火山图清楚地表明，在 TKO 条件下，大多数差异Peak都增加了 NRF1 结合（图 8.2B）。对于第 9 章中描述的下游分析，我们输出TKO 特异的和共享Peak list（ 由于 WT 特定peak的数目较低，因此忽略了它们）。为此，根据严格的显著性阈值（校正后的 p 值 ≤ 10^-10 ）选择显著的Peak，而要求共享峰显示小于两倍的变化。这些选择标准在这些阈值之间留下了一个灰色的峰值区域，允许两个组之间的清晰分离。
8.4.2 Diffbind
-
DiffBind 是一个封装工具，它将
R/Bioconductor软件包DESeq、DESeq2或edgeR应用于 ChIP-seq 数据（ 默认：DESeq2 ）。它提供了一个简单的流程，并在几个步骤中进行数据可视化，这允许检测重复一致性和条件之间的总体差异。它需要一个类似于ChIPQC的样本表 （ 参见附件中的NRF1_Sample_Sheet.csv），该样本表以data.frame或CSV格式总结有关样本所需的信息。library(DiffBind) # 加载数据 NRF1 = dba(sampleSheet="NRF1_sample_sheet.csv") -
导入数据后，通过生成显示成对欧几里德距离的热图以及由此产生的样品的层次聚类，可以基于
Peak位置检查样品的总体相似性（ 类似于8.1 节中生成的热图 ）。 -
DiffBind 中的下一步定义将用于比较的一致
Peak值集。函数dba.count()中的minOverlapp参数设置了Peak必须出现在给定数量的样本中才能包含到一致集合中的要求（ 默认值：2 ）。使用一致Peak集合处的富集（ 输入标准化 reads数 ）值重复成对距离的热图可视化通常将改善按样本类型的聚类。minOverlay可以降低到 1，以包括所有Peak，这可能会增加噪音，但会减少假阴性。函数dba.contrast()的作用是：定义使用样本表中的哪一列进行比较。通过在block参数中指定混杂参数的列，可以考虑实验设置中的混杂参数。minMembers参数设置每个比较组中所需的最小唯一样本数（ 默认：2 ）。运行差异结合分析的最后一个函数是dba.Analyze()。这里需要考虑的一个重要参数是bFullLibrarySize，它决定是否将整个文库的大小用于标准化（ 默认值：TRUE ）。这对于预期全局变化的比较是可取的。相反，仅考虑峰值区域内的 reads 数（bFullLibrarySize = false）适用于预期不到一半的峰值将发生变化的比较。注意：全局变化可能需要使用
spike-in进行标准化（ref7 / ref8）。
# 差异结合分析
plot(NRF1)
NRF1 = dba.count(NRF1, minOverlap=2)
NRF1 = dba.contrast(NRF1, categories=DBA_CONDITION, minMembers=2)
NRF1 = dba.analyze(NRF1, method=c(DBA_DESEQ2))
-
DiffBind 可以使用函数
dba.report()导出差异Peak并可视化结果。NRF1.DB=dba.report(NRF1) dba.plotPCA(NRF1, contrast=1, label=DBA_CONDITION) dba.plotMA(NRF1) dba.plotBox(NRF1) dba.plotHeatmap(NRF1, contrast=1, correlations=FALSE) -
从DiffBind 获得的差异
Peak上的 PCA 图 再次表明 WT 和 TKO 样品在重复之间紧密聚类，并且在条件之间很好地分离。MA图显示了分析中每个Peak的 log2 转换的富集倍数变化与 log2 转换的平均富集（图 8.2 C）。在 FDR < 5% 的默认阈值下，DiffBind 鉴定到了 6946 个差异结合Peak；其中大多数在 TKO 细胞中表现出更强的结合。请注意，这个阈值比我们在第8.4.1节中的 DESeq2 分析中的阈值更宽松，反映在鉴定到更多数量的差异Peak。箱式图显示了与那些在 WT 细胞中显示明显的更多 （ + ）或更少（ - ）结合的Peak相比， log2 转化的富集在所有Peak中的分布（图 8.2 D）。在 NRF1 数据中，我们观察到所有Peak的 reads 密度都有很大的变化，表明全库大小标准化更适合于此数据集。最后，可以使用热图对每个差异Peak的每个重复的标准化富集进行可视化和聚类（图 8.2E）。这再次证实，大多数Peak在 TKO 中显示出更高的富集，并且重复之间具有相似的水平，因此聚在一起。
图 8.2 差异结合分析
（ A ）WT 和 TKO ChIP-seq 样本的主成分分析
（ B ）合并后的 Peaks 中倍数与显著性的火山图
（ C ）MA-plot
（ D ）箱式图
（ E ）在 WT 和 TKO 细胞中差异结合的 Peak 的聚类热图
image
进一步阅读
- Bardet, AF, He, Q, Zeitlinger, J, Stark, A (2011). A computational
pipeline for comparative ChIP-seq analyses. Nat Protoc, 7(1):
45–61. - Steinhauser S. et al. (2016). A comprehensive comparison of tools for differential ChIP-seq analysis. Brief Bioinform, 17(6): 953–966.
ref1：A computational pipeline for comparative ChIP-seq analyses.
ref3：Measuring reproducibility of high-throughput experiments
ref4：An introduction to computational tools for differential binding analysis with chip-seq data.
ref5：Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
ref6：edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
网友评论