R语言中phyper做超几何检验
文章目录
phyper
基因集分析有两种,一种是GSEA(gene set enrichment analysis),需要根据所有基因logFC排序,根据rank来算enrichment score,还有一种是ORA(Over-representation analysis),看选出的显著的基因集是否和已知的基因集显著相关。对于ORA分析,常用超几何分布来检验。在R语言中,用的函数是phyper。
|
|
$$ phyper(…, lower.tail = FALSE)计算的是Pr(X>x)在计算p值时,x-1计算的才是Pr(X≥x) $$
case study
我有一个队列的数据,每个病人的分子分型和药物反应的信息如下表格。
Subtype1 | Subtype2 | Subtype3 | Subtype4 | |
---|---|---|---|---|
not response | 5 | 0 | 9 | 3 |
response | 1 | 3 | 1 | 6 |
Q1:亚型3是否和not response相关?
|
|
写到这里,很有意思的一件事,我想到了另外检验方法,用U(Z)检验来做test。首先把表格合并,得到
non-Subtype3 | Subtype3 | |
---|---|---|
not response | 8 | 9 |
response | 10 | 1 |
零假设:not response在3型和非3型中的分布频率p1和p2是一致的
$$ p_1=8/18,p_2=9/10,加权平均为\overline{p}={(8+9)}/(8+10+9+1)=17/28 $$ $$ 频数的标准差为S_{p1,-p2}=\sqrt{p_1p_2({{1\over{n_1}} + {1\over{n_2}} })}=\sqrt{{8\over{18}}*{9\over{10}}({{1\over{18}} + {1\over{10}} })}=0.2494 $$ $$ 因为上面第一行的频数小于30,需要进行连续校正,得到u值为 $$ $$ u_c = {{{|p_1 - p_2|} - {0.5\over{n_1}} - {0.5\over{n_2}}}\over{S_{p1,-p2}}} = {|{{8\over{18}} - {9\over{10}} | - {0.5\over{18}} - {0.5\over{10}} }\over0.2494} = 1.514 $$ 经查表可知p值在0.0643和0.0655之间,我们取0.06445,这是双侧检验,如果是单侧检验的话,p值应该是0.06445/2=0.032,和我们用超几何检验的出来的类似。
因为样本数小于40,我们还可以用fisher精准检验来算P值,p值为0.0407
|
|
####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################
文章作者 zzx
上次更新 2023-04-11