我们在进行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) |
运行你第二个脚本时好像不行。出现
错误于.checkKeys(value, Rkeys(x), x@ifnotfound):
value for “XXX” not found
请问是怎么回事啊？
@xqyang,
试试
e<-mget(cell.death.GO.terms,org.Hs.egGO2ALLEGS,ifnotfound=NA)
@admin,
非常感谢！
你的一个帖子《bioconductor全自动处理基因芯片CEL文件R程序》，方便给我看一下吗？好像要密码的。因为我也做芯片分析的，很想学习一下。可以的话发我Email好吗（sirxqyang@gmail.com）？谢谢！
博主您好，我在尝试第一个脚本的genelist.entrez.id提示找不到该对象，是要换成具体的id么，还是要建一个id的列表？
另外第二个脚本的e<mget(XXX,org.Hs.egGO2ALLEGS,ifnotfound=NA)时，好像用具体的GO编号可以得到结果，但是用具体的GO Term不行，不知道是什么原因，您能解答一下吗？
这是因为在org.Hs.egGO2ALLEGS中左键值是GO.ID, 右键值是其别名。在使用mget时，是以左键值为索引的。这里我文中的GO Terms可能会让人产生一些误会。
明白了，谢谢
另外，请问那个genelist.entrez.id的问题是怎么回事？
是不是应该在之前加上genelist.entrez.id <- keys(org.Hs.eg.db, keytype="ENTREZID") ?
genelist.entrez.id是指任何一串human的entrez id.
知道了，谢谢！