利用KEGGREST包提取所有代谢相关的基因

写笔记主要是为了自己进步，也是方便自己查看。R语言新手，如有不对请指正。

直接上代码吧

rm(list = ls())
library(KEGGREST)
library(stringr)

##提取KEGG所有的通路
hsa_path <- keggLink("pathway","hsa")

#KEGG共收集8443个相关通路的基因
length(unique(names(hsa_path)))

#涉及通路353个
length(unique(hsa_path))

#代谢相关的通路在KEGG中以“00……”开头，提取”00……“开头通路共84个
pathway= unique(hsa_path)[grepl('hsa00',unique(hsa_path))];length(pathway)


#提取这84个通路的基因
hsa_info <- lapply(pathway, keggGet)

分别把基因名字和基因的ID提取出来，合成一个数据框
gene_symbol = unlist(lapply( hsa_info , function(x) {
  g = x[[1]]$GENE
  str_split(g[seq(2,length(g),by=2)],';',simplify = T)[,1]
}))

gene_ID = unlist(lapply( hsa_info , function(x) {
  g = x[[1]]$GENE
  paste0("hsa:",g[seq(1,length(g),by=2)])
}))

genelist <- data.frame(gene_ID,gene_symbol)

#去重复,提取84个代谢通路中1728个基因
genelist <- genelist[!duplicated(genelist$gene_symbol),]
length(genelist$gene_symbol)

head(genelist)

#保存为Rdata，随时可以加载使用
save(genelist，file = "metabolism.Rdata")

最终结果长这样子：
##   gene_ID gene_symbol
1  hsa:3101         HK3
2  hsa:3098         HK1
3  hsa:3099         HK2
4 hsa:80201       HKDC1
5  hsa:2645         GCK
6  hsa:2821         GPI

KEGG上的基因数量和通路数量都是不断更新的，所以KEGGREST包能提取的基因数量也在不断变化。我去年运行的时候总基因8192个，现在已经有8443个了，以后也会不断增加。

KEGG官网上也有少量的代谢通路不以00开头，以01开头。数量太少，且不是重要的通路，所以我省略了。想要提取更全的基因基的朋友可以自己加上。

最后，本文参考自生信技能树曾老师的代码，感谢曾老师的分享。