原来本菜鸟做玩差异分析后通路注释都是喜欢跑到DAVID里面的网页搞定,感觉还要整理表格之类,贼烦。最近翻了翻群主大大的富集分析帖子,发现了该bioconductor神包clusterprifiler,于是参照群主自己探索了一下用法,确实太好用了,代码非常简单和简洁。据说写该包的作业也是NB人物。 接下来切入正题:先要转化成entrezID,这个参考博客中的http://www./thread-404-1-1.html,大神写很清楚了 1.做GO和KEGG: 接下来就是clusterprofiler的函数: library(clusterProfiler) ego=enrichGO(OrgDb="org.Hs.eg.db", gene = gene1,pvalueCutoff = 0.01,readable=TRUE)#一行代码完成GO富集,gene1就是上图的entrezID的list Pvalue可自己调节 接下来输出结果: write.csv(summary(ego),"G-enrich.csv",row.names =F) 同样的KEGG: ekk <- enrichKEGG(gene = gene1, organism = 'hsa', pvalueCutoff = 0.05) write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F); 那么这样我们就用4行代码完成了GO和kegg的富集分析。 后面我在bioconductor的Rscript里还发现这包还能玩GSEA富集.....(TM好强啊我去) 代码也是极其简单的: gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler") c5 <- read.gmt(gmtfile)#前面俩行大概就是来选择GSEA的数据库用的,c5.cc.v5.0.entrez.gmt这里他事例用的是GO的cc数据库,这个可以调整的,参考GSEA的网站来选择 egmt <- enricher(gene1, TERM2GENE=c5) write.table()write.csv(summary(egmt),"Gsea-enrich.csv",row.names =F)#后面俩行就调整你的GENElist,我这里是上述的gene1。 到此,GO,Kegg,GSEA完成。这包真的是强。。。。 |
|