一、MetaboAnalysis
网站功能是基于R语言实现，也有同名的R包可供使用。
网址：https://www.metaboanalyst.ca/
这种方式比较简单粗暴，支持的ID类型有化合物名称, HMDB ID, KEGG ID, PubChem CID等等，但其弊端是富集分析基于的通路为人类的代谢通路，因此最后可能会富集到一些奇怪的结果。
二、massdatabase包
massDatabase一共支持12个在线数据库, 支持下载物种的化合物/通路数据库并将它们转换为特定的数据结构，然后我们结合clusterprofiler包从而非常方便的进行代谢物富集分析。这里以研究物种—弓形虫(Toxoplasma gondii，缩写为tgo)的代谢通路和化合物数据库为例。
KEGG Organisms缩写查询：https://www.genome.jp/kegg/catalog/org_list.html
下载tgo的kegg通路数据库
主要利用3个函数功能 ：1. download_kegg_pathway 2. read_kegg_pathway 3. convert_kegg2metpath ，分别进行下载，读取和转换作用。
#安装R包 remotes::install_github("tidymass/massdatabase", dependencies = TRUE) ## 下载tgo的通路代谢物关系 download_kegg_pathway(path = "kegg_tgo_pathway", sleep = 1, organism = "tgo") ## 然后读取所下载的数据 tgo_data <- read_kegg_pathway(path = "kegg_tgo_pathway") class(tgo_data) ## 转变数据库形式 kegg_pathway_database <- convert_kegg2metpath(data = tgo_data, path = ".") ##查看有多少通路 length(unlist(kegg_pathway_database@pathway_name)) [1] 104
观察一下kegg_pathway_database的结果，共有104个通路，我们只需要关注这5个内容即可，一个最简单的思路就是通过for循环遍历将各个通路的信息提取出来，当然也可以用lapply等函数解析列表。
但是要注意一个细节，就是有的通路是没有gene或者compound对应, 也没有前两层级的KEGG信息，因此这部分需要我们for循环中多添加一个条件，如果遍历到这几行，我们赋值为None。
提取通路与gene对应信息
path_num <- c(68:76) ## 排除没有对应gene的几个通路 ## 提取1到67 77 到104的通路 result1 <- NULL for (i in 1:104) { a <- tgo_data[[i]]$pathway_id # 提取id b <- tgo_data[[i]]$pathway_name # 提取名称 e <- tgo_data[[i]]$gene_list # 提取gene对应关系 c <- tgo_data[[i]]$describtion if ( i %in% path_num){ d <- 'None;None' }else { d <- tgo_data[[i]]$pathway_class # 提取前两层级注释 Pathway_re <- e%>% mutate(pathway_id = a, pathway_name = b, pathway_class= d) %>% separate(pathway_class,into =c("Pathway1","Pathway2"),sep = ';',remove = FALSE) ##前两层分开 result1 <- rbind(result1,Pathway_re)} } write.table(result1, 'tgo_KEGG_path.tsv', sep = '\t',row.names = FALSE,quote = FALSE)
gene、注释及前几层通路信息如下：
提取通路与化合物对应信息
result2 <- NULL com_num <- c(38,68:76,78:94,96,98,103) ## 排除没有对应化合物的几个通路 for (i in 1:104) { a <- tgo_data[[i]]$pathway_id # 提取id b <- tgo_data[[i]]$pathway_name # 提取名称 f <- tgo_data[[i]]$compound_list # 提取化合物对应关系 # c <- tgo_data[[i]]$describtion if ( i %in% com_num){ d <- 'None;None' }else { d <- tgo_data[[i]]$pathway_class Compound_re <- f %>% mutate(pathway_id = a, pathway_name = b, pathway_class= d)%>% separate(pathway_class,into =c("Pathway1","Pathway2"),sep = ';',remove = FALSE) result2 <- rbind(result2,Compound_re) } } write.table(result2, 'tgo_KEGG_com.tsv', sep = '\t',row.names = FALSE,quote = FALSE)
化合物、名称及通路对应信息如下：
代谢物富集分析
有了上面的compound对应信息，我们只需提取其中几列作为背景文件轻松用clusterprofiler包进行富集。
library(clusterProfiler) ##制作背景文件 com_ano <- result2 %>% select(KEGG.ID,pathway_id,pathway_name) ## 富集ID gene_select <- sample(com_ano$KEGG.ID,100) #KEGG 富集分析 #默认以所有注释到 KEGG 的基因为背景集，也可通过 universe 参数指定其中的一个子集作为背景集 #默认以 p<0.05 为标准，Benjamini 方法校正 p 值，q 值阈值 0.2 #默认输出 top500 富集结果 kegg_rich <- enricher(gene = gene_select, TERM2GENE = com_ano[c('pathway_id', 'KEGG.ID')], TERM2NAME = com_ano[c('pathway_id', 'pathway_name')], pvalueCutoff = 0.05, pAdjustMethod = 'BH', qvalueCutoff = 1, maxGSSize = 500) dotplot(kegg_rich) #富集气泡图 cnetplot(kegg_rich) #网络图展示富集功能和基因的包含关系 emapplot(kegg_rich) #网络图展示各富集功能之间共有基因关系 heatplot(kegg_rich)
