前言
思考
该问题中卡方检验的理解
总TSS上下游1kb范围长度之和
总长度计算
卡方检验
转换为reads
注意
今天进行一个思维的简单训练。
现在有一个问题，我们来思考下：如果我们做组蛋白的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.gtf116561151总长度计算
下载染色体长度：
$ 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"|bc2971708681卡方检验
接下来我们去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=1e9m=1for (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
发现只要当比例最低达到0.03777549时就会出现显著变化了。
转换为reads既然知道了比例，那么接下来就是根据这个比例去计算了。
首先，reads必须要占满所有的TSS，即所有基因的TSS上都要有reads。
我们假设reads长度为x，那么此时需要的reads数最低需要：116561151/x
然后我们根据上述计算比例，就可以计算出此时的总reads数应该有：
116561151/x/0.03777549=4403155/x
也就是说，如果测序长度平均为30bp，那么理论上我们需要的reads数目为：4403155/30=146772条即可。
注意因为这里只是单纯从理论上去计算，read的估算其实是存在很大误差的，例如：
不可能我们在测序时仅仅要求所有TSS的覆盖度只有1x
所有TSS上下游区域的划分并不是连续的，计算reads时我们当作了连续的进行处理
其他所以为了保险，我们可以将最终结果乘以一个系数，但是这个系数具体是多少我这里也不知道，但是我也许就估算为10（即TSS区域为10x）。
那么在这个结果中，我们需要的reads则为1.5e6条。