注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Bioinformatics home

 
 
 

日志

 
 

富集性分析  

2012-06-23 01:14:43|  分类: 生物信息编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

转自 http://ygc.name/2012/04/28/enrichment-analysis/#comment-3188

经常看到一些饼图,描述某些事物的组成,比如说有钱人的学历分布,然后我们可以看到高学历所占比例并不高,根据这个比例下结论通常是错的,这些比例说明不了问题,如果把各种学历在总体人口中的分布做为背景进行考虑的话,你就会发现学历还是有点用的。

当我们用组学测定了一大堆分子之后,我们希望站在更高的角度去看这些分子和那些生物学过程相关。那么通常各种注释,对这些基因/蛋白进行分类,那么从分类的比例上,是不能草率下结论,正如上面有钱人学历分布的例子一样。我们需要把总体的分布考虑进去。

和某个注释/分类是否有相关性,把基因分成属于这一类,和不属于这一类两种,这就好比经典统计学中的白球和黑球的抽样问题。也可以列一个2x2的表,进行独立性分析。

以文章Gene Expression in Ovarian Cancer Reflects Both Morphology and Biological Behavior, Distinguishing Clear Cell from Other Poor-Prognosis Ovarian Carcinomas所鉴定的差异基因为例。

73个差异基因的Symbol,我把它转为 entrezgene ID得到57个(漏掉的不管它,只是做为一个例子):

> eg  [1] "7980"   "3081"   "3162"   "3059"   "1545"   "1917"   "6696"   "5797"  [9] "6648"   "10397"  "6781"   "5817"   "1282"   "1284"   "6948"   "7077" [17] "5744"   "8566"   "1368"   "1474"   "11015"  "3306"   "728441" "2678" [25] "4286"   "3929"   "5095"   "2064"   "1428"   "6590"   "3569"   "2745" [33] "3912"   "978"    "5950"   "6539"   "9445"   "5004"   "9971"   "7453" [41] "2719"   "1410"   "1311"   "4653"   "4162"   "5358"   "3484"   "3486" [49] "2261"   "307"    "1672"   "4837"   "22795"  "486"    "4118"   "3915" [57] "10140" 

下面测试一下这些基因和化学刺激响应的相关性。

goid <- “GO:0042221” # response to chemical stimulus 

那么做为背景,总体基因为N,属于“化学刺激响应”这个分类的基因有M个。

allgeneInCategory <- unique(get(goid, org.Hs.egGO2ALLEGS)) M <- length(allgeneInCategory) N <- length(mappedkeys(org.Hs.egGO)) 

样本的大小是n,属于“化学刺激响应”这个分类的基因有k个。

n <- length(eg) k <- sum(eg %in% allgeneInCategory) 

白球黑球问题,最简单的莫过于用二项式分布,从总体上看,要拿到一个基因属于“化学刺激响应”这个分类的概率是M/N。那么现在抽了n个基因,里面有k个基于这个分类,p值为:

> 1-sum(sapply(0:k-1, function(i) choose(n,i) * (M/N)^i * (1-M/N)^(n-i))) [1] 8.561432e-10 

二项式分布,是有放回的抽样,你可以多次抽到同一基因,这是不符合的。所以这个计算只能说是做为近似的估计值,无放回的抽样,符合超几何分布,通过超几何分布的计算,p值为:

> phyper(k-1,M, N-M, n, lower.tail=FALSE) [1] 7.879194e-10 

如果用2x2表做独立性分析,表如下:

> d <- data.frame(gene.not.interest=c(M-k, N-M-n+k), gene.in.interest=c(k, n-k)) > row.names(d) <- c("In_category", "not_in_category") > d                 gene.not.interest gene.in.interest In_category                  2613               28 not_in_category             15310               29 

这个也有很多方法可以做检验,经典的有卡方检验和fisher's exact test。

如果用卡方检验来做。

> chisq.test(d, )   Pearson's Chi-squared test with Yates' continuity correction  data:  d X-squared = 51.3849, df = 1, p-value = 7.592e-13 

对于2x2表来说,卡方检验通常也只能做为近似估计值,特别是当sample size或expected ell count比较小的时候,计算并不准确。

fisher's exact test,名副其实,真的就比较exact,因为它使用的是超几何分布来计算p值。这也是为什么fisher's exact test和超几何模式计算的p-值是一样的,

> fisher.test(d)   Fisher's Exact Test for Count Data  data:  d p-value = 7.879e-10 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval:  0.1013210 0.3089718 sample estimates: odds ratio  0.1767937  

通常各种软件做GO富集性分析,都是使用超几何分布进行计算。

IPA软件则是使用fisher's exact test来检验基因在某个网络中是否富集。

从这个例子上看,chi-squared test的表现最差,binomial做为p值的近似估计还是不错的,因为计算较为简单。

富集性分析应用范围非常广,从Disease OntologyGene Ontology, KEGG, 到Reactome Pathway等等

  评论这张
 
阅读(1632)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017