批量GO-KEGG富集分析注释何必自己写脚本,看Y叔的神器
之前我有多个基因集的时候,比如下面这个:
> 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)
还能顺便出图,这里不方便用我自己的图,就直接给Y叔的样图给大家看:
来自于:
https://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html