orthofinder

批量GO-KEGG注释何必自己写脚本

2018-12-24  本文已影响53人  因地制宜的生信达人

批量GO-KEGG注释何必自己写脚本

之前我有多个基因集的时候,比如下面这个:

> head(moduleGenes[,2:3])
     module             chr
1     brown ENSG00000175899
2      blue ENSG00000166535
3 turquoise ENSG00000256069
4      grey ENSG00000128274
5 turquoise ENSG00000205002
6      blue ENSG00000283199

简单统计得到:

  black      blue     brown     green      grey       red 
       56      1233       656       337       320       102 
turquoise    yellow 
     1783       524 

其实就是WGCNA得到的不同module的基因集。

对于这样的结果,我都是自己写脚本一个个使用 enrichKEGG 函数,然后把返回值合并。

library(clusterProfiler)
library(org.Hs.eg.db)
tmp=split(moduleGenes,moduleGenes$module)
enrich_kegg_results <- lapply(tmp,function(x){
  g_diff=unique(x[,3])
  gene.df <- bitr(g_diff, fromType = "ENSEMBL",
                  toType = c("SYMBOL", "ENTREZID"),
                  OrgDb = org.Hs.eg.db) 
  kk <- enrichKEGG(gene.df[,3],'hsa', pvalueCutoff = 0.99, qvalueCutoff=0.99 )  ## g_diff should be ENTREZID
  dat=kk@result[,c(1,2,4:6)] 
  return(dat)
  
}) 

re = do.call("rbind",enrich_kegg_results) 
re$sample=rep(names(enrich_kegg_results),sapply(enrich_kegg_results, nrow))
head(re) 
enrich_mat <- tidyr::spread(re[,c(2,5,6)],'sample','p.adjust')
rownames(enrich_mat)=enrich_mat[,1]
head(enrich_mat)
enrich_mat=enrich_mat[,-1]

但是我看clusterProfiler包自带的说明书,发现了 一个好用的函数,compareCluster ,就可以把上面的代码改写。

首先把基因ID简单转换:

> head(mydf)
     module ENTREZID
1     brown     2268
2     green     3075
3 turquoise    22875
4 turquoise     1080
5 turquoise       26
6    yellow    23072

然后两行代码即可:

formula_res <- compareCluster(ENTREZID~module, 
                              data=mydf, fun="enrichKEGG")

head(as.data.frame(formula_res))
 dotplot(formula_res)

还能顺便出图:https://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

上一篇下一篇

猜你喜欢

热点阅读