p.adjust()
library("fdrtool")
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
http://www.360doc.com/content/17/1228/11/50153987_717073620.shtml
http://www.360doc.com/content/18/0914/23/52043738_786760347.shtml
在生信分析中,隔三差五地就需要和p值探讨是否显著差异,还要搬出FDR对p值进行校正。让每个基因根据p值大小从小到大排个队,拿个号牌,然后把自己的p值乘上总基因数,再除以自己号牌上的数就是FDR校正后的p值啦。这个过程用数学语言表示为:
其中,q-valuei是校正后的p值,length(p)是总基因数,rank(p)是每个基因排队的号牌。
那么,为什么要用这个公式完成校正呢?且听小锐从原理说起。
寻找差异基因的漫漫长路,从假设开始……
我们在处理宏基因组差异基因的选择时,需要对两个样本的每个基因进行一次假设检验。如果我们有m个基因,那么我们就要做m次假设检验。每一次的假设检验的零假设H0为:两个样本的这个基因没有显著性差异。其中有m0个零假设是正确的,即这个基因在两个样本中确实没有显著性差异;但有m1=m-m0个零假设是错误的,即两个样本的这个基因是有显著性差异。m次检验之后,被拒绝的零假设的个数记为R。为了方便记忆,可用一张表格来表示假设检验的结果:
m次假设检验的结果
此时,我们需要多考虑一点:找到的“差异基因”里面,是否有找错了的呢?
为了解决这一问题,统计学家提出了控制错误率的标准。其中,Benjamini和Hochberg(1995)提出的FDR准则最为流行,它的标准相对来说更宽泛,在实际应用中能够获得更大的功效。也就是说,通过这种方式,我们可以在错误率很低的情况下找到真正有差异的基因。
根据“m次假设验证的结果”(记性不好的童鞋可以翻看前面的表格),我们可以得到FDR准则,即要求控制错误拒绝率。令Q=V / (V+S),它表示被错误拒绝的零假设数目占所有被拒绝的零假设数目的比例。Q也是一个不可观测的随机变量。定义FDR如下:
当R=0时,定义Q=0
Benjamini和Hochberg提出了以FDR作为多重检验的准则,但是其检验的方法采用的是Simes(1986)提出的算法。
怎么修正?
重回开始的例子:在做差异基因的筛选时,有两个样本,m个基因。根据上述,我们需要做m次假设检验,会得到m个p值,分别为p1, p2,…, pm。将这m个p值由小到大排序,得到另外一组p(1)≤p(2)≤…≤p(m),其中下标表示p值的排序位置。做这m次假设检验的控制错误率为
上式表示期望,也就是均值,即平均控制错误率。
假设我们按照每个P值小于0.05这个标准,我们一共可以挑出V+S=R个差异基因,这R个差异基因中有V个其实是没有差异表达的,S个是真正有差异表达的基因。如果有10,000个基因接受显著性检验,其中有90%,也就是9000个(即m0)满足原假设,那我们依然会从这些满足原假设的基因中找到450(即V)个有差异的,即使剩下1,000个(S)全部都被认为是差异基因,如果不利用FDR值进行修正,我们找到的差异基因中450/(450+1000)=31%(即Q=V/(V+S))的差异基因是错误的。
为了使得这个错误率很低,我们要求平均错误控制率要小于一个预先设定的q值(比如0.05)即E(Q)≤q*
在统计学上,这也就等价于控制FDR不能超过5%。将E(Q)≤q*进行拆分:
因为V代表实际的差异基因,是固定不变的;要想使得控制错误率达到最小,就必须要求分母(V+S)最大化。
我们做了次假设检验,得到了V+S个差异基因(V+S个基因的P值≤0.05),V+S也表示排序后的p值小于0.05中的最大的那个p值对应的排序号(秩),即最大的那个p值为上述的pi,排序号为上述i。而此时V等于pi*m,所以E(Q)≤q*可以写成
修正的q*值取最小值,即为我们平时工作中用到的修正的FDR值q-valuei。这里的公式即为Bonferroni型多重检验过程中的公式。也是开始FDR的计算公式:
最后是FDR校正后的p值计算的一个小例子。
大家可以移步该网页查看
http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html