ChIP

吴小邪 2023-05-23 发布于湖南

展开全文

前言

今天进行一个思维的简单训练。

现在有一个问题，我们来思考下：如果我们做组蛋白的ChIP-seq分析，那么对于我们组蛋白数据，最少应该有多少的测序深度，才能够有效的进行数据分析呢？

这里涉及到一系列的统计学知识。

我这里的思考不一定完全正确，仅仅写出来让大家一起来思考。

思考

什么才是一个合格的组蛋白ChIP-seq数据呢？那么肯定是在有限的reads中，这些reads大都可以富集集中到一些基因组上的区域中。

但是结合不同组蛋白自身特点，所以这里我们将数据分为2种：

sharp marker

broad marker

之所以做这种分类，主要考虑到，基因组上broad marker多呈现为大片段连续分布，而sharp的marker则多呈现为在TSS区域的富集。

为了简化我们的思考，我们这次主要考虑sharp的marker。也就是考虑需要多少reads时，数据可以表现出在TSS有富集的信号。

为了对TSS范围做出定义，我以TSS上下游1kb标准，需要考虑2kb的范围内有富集。

描述到这里，其实如果有一定统计基础的同学，差不多已经知道用什么统计方法来进行计算了——超几何检验。

但是这里不能用超几何检验，因为在超几何检验中，球的总数需要是固定的，并且白球和黑球总数也是固定的，超几何检验计算的是每次抽取的球中包含白球总数的概率。

而在我们的实例中，落入到TSS上下游1kb范围的reads以及落入到其他区域的reads总数是不固定的。

所以我们不能用超几何检验。那我们应该用什么统计方式呢？

我们可以用卡方检验。

该问题中卡方检验的理解

在这个问题中，因为每个read都可以随机分布在整个基因组上，而我们关心的只是，这些reads能不能主要落在TSS上下游1kb的范围，所以我们需要计算出：

整个基因组的长度

所有TSS上下游1kb的长度之和

基因组上除去TSS上下游1kb范围之外的长度

根据卡方检验的原理，我们可以把思路转变一下：去比较落在TSS上下游1kb范围内的reads数和落在其他区域内的reads数什么时候可以有显著差异。只要当两者比例存在显著差异时，我们就可以得出结论了。

既然比较差异，那么我们可以先假设没有富集现象的存在，这是reads是均匀分布的，所以基于该假设的结果即为：所有TSS上下游1kb的长度之和比上基因组上除去TSS上下游1kb范围之外的长度。

如果reads的分布存在富集，那么在该假设下，我们应该存在另外一个比例值。

而我们的标准是在P<0.05时，计算出这个比例值。

这样我们的思路就确定了，我们需要计算的结果有：

总TSS上下游1kb范围长度之和

基因组上除去TSS上下游1kb范围之外的长度

而我们的结果则是计算出在卡方检验P<0.05时，在reads分布存在富集的假设下，最小的比例值。
考虑到不同基因如果距离较近，那么他们之间会存在有overlap，所以这里我们先提取出所有基因TSS上下游1kb范围的坐标，然后对其取并集：

$ grep -v ^# gencode.v40.annotation.gtf |awk -F'\t' -v OFS='\t' '{if($3~/gene/) {if($4-1000<0) print $1,0,$4 1000;else print $1,$4-1000,$4 1000}}'|bedtools merge -i - >gene.gtf

然后我们计算所有并集的区域之和：

$ awk -F '\t' 'BEGIN{sum=0} {b=$3-$2;sum=sum b} END{print sum}' gene.gtf
116561151

总长度计算

下载染色体长度：

$ wget https://hgdownload-test.gi.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes

计算常见的染色体长度之和：

$ egrep 'chr[0-9XY]{1,2}.[0-9]{1,}' hg38.chrom.sizes | awk -F'\t' 'BEGIN{sum=0} {sum=sum $2} END{print sum}'
3088269832

这样我们就知道了基因组上除去TSS上下游1kb范围之外的长度：

$ echo '3088269832-116561151'|bc
2971708681

卡方检验

接下来我们去R中进行计算。

首先计算出均匀分布时，总TSS上下游1kb范围长度之和占总长度的比例

116561151/3088269832
# 0.03774319

然后写循环计算：

rm(list = ls())
options(stringsAsFactors = F)
library(stringr)
library(ggplot2)
dat <- data.frame(proportion=NA,
                  p_value=NA)
base=1e9
m=1
for (i in seq(0.03774319,0.03774319 0.01,length.out=1e5)) {
  mytable <- data.frame(men = c(0.03774319*base,(1-0.03774319)*base),
                        women = c(i*base,(1-i)*base))
  p <- chisq.test(mytable)$p.value
  dat[m,'proportion']=i
  dat[m,'p_value']=p
  m <- m 1
}
dat$fdr <- p.adjust(dat$p_value)
head(dat[dat$fdr<0.05,])
#     proportion      p_value        fdr
# 324 0.03777549 0.0001512562 0.04900699
# 325 0.03777559 0.0001442724 0.04688854
# 326 0.03777569 0.0001375932 0.04485539
# 327 0.03777579 0.0001312061 0.04290439
# 328 0.03777589 0.0001250991 0.04103250
# 329 0.03777599 0.0001192608 0.03923679