1）对每个基因进行p-value的计算
假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布。已知样本一中唯一比对到基因组的总reads数为N1，样本二中唯一比对到基因组的总reads数为N2，样本一中唯一比对到基因A的总reads数为x，样本二中唯一比对到基因A的总reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算：
图片1.png
2）用FDR错误控制法对p-value作多重假设检验校正
FDR错误控制法是Benjamini于1995年提出的一种方法，通过控制FDR（False Discovery Rate）来决定P值的域值。假设你挑选了R个差异表达的基因，其中有S个是真正有差异表达的，另外有V个其实是没有差异表达的，是假阳性的。实践中希望错误比例Q＝V/R平均而言不能超过某个预先设定的值（比如0.05），在统计学上，这也就等价于控制FDR不能超过5％。
对所有候选基因的p值进行从小到大排序，则若想控制FDR不能超过q，则只需找到最大的正整数i，使得 p(i)≤ (i*q)/m.然后，挑选对应p(1)，p(2)，...，p(i)的基因做为差异表达基因，这样就能从统计学上保证FDR不超过q。 因此，FDR的计算公式如下：
q-value(i)=p(i)*length(p)/rank(p)
首先将n个基因的原始P值由大到小排序，将最大P值赋值为n，最小P值赋值为1。校正P值=原始P值（n/i）。其中n表示所有的基因个数，i表示从小到大第i个P值。例如下表中，共6个基因，因此n=6，从大到小排序后，第1大的原始P值0.0687，其校正P值为0.06873（6/6）=0.0687；第2大的原始P值0.0235，其校正P值为0.0235（6/5）=0.0282；……；第6大的原始P值0.0003，其校正P值为0.0003（6/1）=0.002。
但是要注意第3大的原始P值0.0192，如果按公式计算其校正P值为0.0195*（6/4）=0.0288，但是FDR在计算校正P值时，需要将当前计算值与上一个计算值比较，取二者中最小值。比如，第3大原始P值所计算的校正P值为0.0288，与第2大校正P值0.0282相比，0.0282更小，因此这里的值不是0.0288，而是0.0282。
参考文献：
1.Audic, S. and J. M. Claverie (1997). The significance of digital gene expression profiles. Genome Res 7(10): 986-95.
2.Benjamini, Y. and D. Yekutieli (2001). The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics. 29: 1165-1188.