在墙内搞生物信息学偶尔会遇到一些莫名其妙的问题，比如突然间有天就没办法使用下面这个函数：
### over-representation test
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
问题是我一直在大力宣传这个clusterProfiler包的这种enrich方式，如果我一直推荐的是一个错误的代码，那该多尴尬呀！
怪不得总是有些人问到使用它的各种失败，各种报错，因为我大部分时间都是在墙外所以根本就没办法重复出求助者的错误。
简单看了加源代码，发现是download_KEGG('hsa') 这个函数的问题，这个函数每次使用都会去下载这个数据，如果某些地方某些电脑无法访问这个KEGG的官方网站，那么这个包就没办法使用了。
所以我的第一个问题来了？
一个主打统计学功能函数的R包需要每次都联网吗？
毕竟很多工作场景是不允许联网的，先不说墙内墙外的问题。
虽然无法使用这个包来做注释，但是其统计学原理我是懂的，只好自己写enrich函数了：
library(ggplot2)
library(KEGG.db)
library(org.Hs.eg.db)
load('for_annotation.Rdata')
## KEGG pathway analysis
library(KEGG.db)
tmp=toTable(org.Hs.egPATH)
minS=10;
maxS=500;
GeneID2Path<<- tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
# 超几何分布检验是统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的个数（不归还）
gened=intersect(names(GeneID2Path),gene_up) ## all genes been choosed
geneb=intersect(names(GeneID2Path),gene_all) ## all genes as backgroud
N=length(geneb);N ## number of all the genes
M=length(gened);M ## number of all the choosed genes
kegg_r <- do.call(rbind,lapply(names(Path2GeneID),function(thisP){
# thisP=names(Path2GeneID)[1];thisP
n=length(intersect(geneb,Path2GeneID[[thisP]]));n
if(n>maxS) return(NULL);
if(n<minS) return(NULL);
( exp_count=n*M/N)
k=length(intersect(gened,Path2GeneID[[thisP]]));k
OddsRatio=k/exp_count
p=phyper(k-1,M, N-M, n, lower.tail=F)
#print(paste(i,p,OddsRatio,exp_count,k,M,sep=" "))
return(c(thisP,p,OddsRatio,exp_count,k,n,M,N))
}))
colnames(kegg_r)=c('path_id','p','OddsRatio','exp_count','k','n','M','N')
kegg_r=merge(kegg_r,toTable(KEGGPATHID2NAME),by='path_id')
kegg_r=kegg_r[order(kegg_r$p),]
写完我就又思考了，这个统计学应该是生信工程师的必备技能，那么除了我演示过的超几何分布检验，还有哪些统计学实例是一定要掌握的呢？
大家可以留言推荐一下！
网友评论