前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >R包是否应该每次使用都联网?以及生信必备统计学实例推荐

R包是否应该每次使用都联网?以及生信必备统计学实例推荐

作者头像
生信技能树
发布2018-11-08 15:45:18
1.3K0
发布2018-11-08 15:45:18
举报
文章被收录于专栏:生信技能树

在墙内搞生物信息学偶尔会遇到一些莫名其妙的问题,比如突然间有天就没办法使用下面这个函数:

代码语言:javascript
复制
  ###   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函数了:

代码语言:javascript
复制
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),]

写完我就又思考了,这个统计学应该是生信工程师的必备技能,那么除了我演示的超几何分布检验,还有哪些统计学实例是一定要掌握的呢?

大家可以留言推荐一下!

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2018-10-19,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档