使用R进行Gene Ontology(GO)富集分析过程中如何导出某一特定GO词条下所有的基因 11

我们在进行GO分析时,时常需要得到某一特定的GO term下所有涉及到的基因,但是R当中并没有明确的工具。下面的脚本片断可以解决这个问题。

> library("goProfiles")
> GOTermsList <- as.GOTerms.list(as.character(genelist.entrez.id),"Entrez","org.Hs.eg.db",onto="BP")
> ancestorsList<-getAncestorsLst(GOTermsList,onto="BP")
> yapply <- function(X,FUN, ...) {
+   index <- seq(length.out=length(X))
+   namesX <- names(X)
+   if(is.null(namesX))
+     namesX <- rep(NA,length(X))
+ 
+   FUN <- match.fun(FUN)
+   fnames <- names(formals(FUN))
+   if( ! "INDEX" %in% fnames ){
+     formals(FUN) <- append( formals(FUN), alist(INDEX=) )
+   }
+   if( ! "NAMES" %in% fnames ){
+     formals(FUN) <- append( formals(FUN), alist(NAMES=) ) 
+   } 
+   mapply(FUN,X,INDEX=index, NAMES=namesX,MoreArgs=list(...)) 
+ }
> ancestors<-do.call(rbind,yapply(ancestorsList,function(.ele){cbind(rep(NAMES,length(.ele)),.ele)}))
> colnames(ancestors)<-c("entrez_id","go_term")
> library("ChIPpeakAnno") #使用用其中addGeneIDs可以很方便地在种类ID之间进行转换
> anno<-addGeneIDs(as.character(genelist.entrez.id),"org.Hs.eg.db",c("symbol","genename"),"entrez_id")
> go.anno<-merge(anno,ancestors,by="entrez_id")
> c.cell.death<-lapply(cell.death.GO.terms,function(.ele,x){as.character(x[x$go_term==.ele,"symbol"])},go.anno)

今天又得到高人指点,有了新的更简单的脚本

> library(org.Hs.eg.db)
> help("org.Hs.egGO2ALLEGS")
> e<-mget(cell.death.GO.terms,org.Hs.egGO2ALLEGS)
> t<-lapply(e,function(.ele,x){intersect(as.character(x),.ele)},genelist.entrez.id)

11 thoughts on “使用R进行Gene Ontology(GO)富集分析过程中如何导出某一特定GO词条下所有的基因

  1. Reply xqyang 3月 15,2012 8:19 下午

    运行你第二个脚本时好像不行。出现
    错误于.checkKeys(value, Rkeys(x), x@ifnotfound):
    value for “XXX” not found
    请问是怎么回事啊?

  2. Reply xqyang 3月 28,2012 10:30 下午

    你的一个帖子《bioconductor全自动处理基因芯片CEL文件R程序》,方便给我看一下吗?好像要密码的。因为我也做芯片分析的,很想学习一下。可以的话发我Email好吗(sirxqyang@gmail.com)?谢谢!

  3. Reply hufanglq 12月 6,2014 6:22 上午

    博主您好,我在尝试第一个脚本的genelist.entrez.id提示找不到该对象,是要换成具体的id么,还是要建一个id的列表?
    另外第二个脚本的e<mget(XXX,org.Hs.egGO2ALLEGS,ifnotfound=NA)时,好像用具体的GO编号可以得到结果,但是用具体的GO Term不行,不知道是什么原因,您能解答一下吗?

  4. Reply hufanglq 12月 9,2014 8:12 上午

    是不是应该在之前加上genelist.entrez.id <- keys(org.Hs.eg.db, keytype="ENTREZID") ?

Leave a Reply

  

  

  

%d 博主赞过: