8. Peaks定量与差异分析

东船西舫悄无言
271 阅读4分钟

image.png

  • 单样本：比如有一个脱落酸处理的样本，想研究ZAT6这个转录因子在脱落酸处理的样本中是一个什么情况，它有哪些结合位点，这些结合位点位于哪些基因的附近。
  • 多样本：有时候样本可能是分组的，有一组样本是用脱落酸处理的，另外一组样本是没有用脱落酸处理的。这种情况下就需要进行差异分析。两组样本之间进行比较。

我们现在想知道的是：哪些peak在样本A和样本B中是有明显的差异的。

什么叫做差异？

  • 就是说有一个peak在样本A中有这个peak，而在样本B中没有这个peak；【有和无的差异】
  • 还有一种情况，在样本A中这个peak的丰度，就是那个峰特别特别的高，而在样本B的这个peak中，那个峰特别特别的低。也就是有一个量度上的差异。【峰的高和低的差异】

类似于转录组学里的差异表达分析。

image.png

  1. 在做基因的转录差异分析的时候： 比如在测序数据量都一样的情况下，就不涉及到测序深度导致的影响。

前提：geneA在样本A和样本B中是统一的，样本AB中的geneA的位置都是从同一个地方到同一个地方，中间的外显子的位置都是确定的。

即：都有这个基因，且基因长得一样。

image.png

  1. peak差异分析

image.png

  • 合并生成所有样本的共识性的peak。

image.png

image.png

操作：生成共识peaks

image.png

  • merge的算法

image.png

  • 然后就可以计算落在这些个位置的reads的数目

image.png

定量

image.png

  • 参考课程

image.png

image.png

image.png

2. 转化为 saf 格式，用于 featureCounts 定量

# 方法1
awk 'OFS="\t" {print $1":"$2+1"-"$3, $1, $2+1, $3, "."}' \
  clean_peaks/merge/merged_peaks.bed \
  > clean_peaks/merge/merged_peaks.saf

# 方法2
awk '{print $1":"$2"_"$3"\t"$1"\t"$2"\t"$3"\t."}' clean_peaks/merge/merged_peaks.bed > clean_peaks/merge/merged_peaks.saf

Q：在生物信息学中有2种格式：0-based, 1-based。

image.png

因此，该代码应该改为：

awk '{print $1":"$2+1"_"$3"\t"$1"\t"$2+1"\t"$3"\t."}' clean_peaks/merge/merged_peaks.bed > clean_peaks/merge/merged_peaks.saf

代码如下：

Rscript software/RunFeatureCounts/run-featurecounts.R \
  -b clean_bams/ZAT6_ABA_rep1_final.bam \
  --saf clean_peaks/merge/merged_peaks.saf \
  --isPairedEnd FALSE \
  -o counts/ZAT6_ABA_rep1 \
  1>logs/ZAT6_ABA_rep1_run_feature_counts.log 2>&1

image.png

结果：

image.png

合并丰度矩阵，并进行 TMM：

ls counts/H3K4me3_WT_rep1.count \
  counts/H3K4me3_WT_rep2.count \
  counts/ZAT6_ABA_rep1.count \
  counts/ZAT6_ABA_rep2.count >counts/count.fofn
  
perl software/RunFeatureCounts/abundance_estimates_to_matrix.pl \
  --est_method featureCounts \
  --quant_files counts/count.fofn \
  --out_prefix counts/merged_peaks \
  1>logs/create_count_matrix.log 2>&1

image.png

结果：这步操作就是合并了不同的样本到一个表格。

image.png

后续可以做PCA分析

要用TMM校正后的矩阵。

image.png

需要准备的文件：

1. 样本信息表

image.png

2. merged_peaks.TMM.TPM.matrix矩阵

做PCA和样本相关性分析。

定完量之后就可以进行差异peak的分析。【和转录组里完全一样】

image.png

image.png

结果： image.png

  • 文章
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02686-y#Abs1

image.png

软件测评：

image.png

总结起来就4个软件： Diffbind, MACS2, edgeR, DESeq2.

其实Diffbind就是整合了edgeR+DESeq2。

MACS2就是call peak那一步用的，自带一个差异分析的程序；但这个程序不支持生物学重复。

所以我们就用：edgeR和DESeq2就好了。

  • 存在一个问题：

    • 跑完之后，会给你算一个Fold Change出来，转录组中习惯叫Fold Change，ChIP-seq里面叫Fold Enrichment。
    • 还会给你算一个pvalue, qvalue出来。

  • Fold Change是可靠的 √

  • pvalue, qvalue没那么可靠，因为edgeR, DESeq2等给转录组开发的软件，计算p值，q值的时候，依据的数据的模型是转录组的模型。转录组的模型应用到chip-seq中来还是有差异的。

  • 所以算出来的pvalue, qvalue仅有参考价值，没有那么重要。

    • 比如：peakA的p值是0.001，peakB的p值是0.01，peakA就是比peakB差异更明显，这个没什么问题；
    • 但当我们筛选有意义的peak的时候，peakA是0.001，peakB是0.06，按照统计学，必须卡0.05的peak，那么peakB就认为不是了，这也没有必要。假如peakB调控的区域正好位于基因上游的promoter区域，而且这个基因又特别特别重要，你想在文章里说一说这个基因，但发现pvalue是0.06，不具有统计学的显著性，我觉得大可不必，该说就说，没什么问题。
avatar
文章
阅读
粉丝