R语言获取人类KEGG通路的所有基因或代谢物

宏发饰业

自行写代码进行KEGG富集分析时需要制作富集背景文件，将所有代谢通路中包含的所有基因（蛋白质）或代谢物提取出来，生成背景文件。
KEGGREST：A package that provides a client interface to the Kyoto Encyclopedia of Genes and Genomes (KEGG) REST server.
该软件包为使用者提供了KEGG数据库的API接口，首先简单看一下它的基本使用：
# 查看KEGG包含的子数据库
> listDatabases()
[1] "pathway"  "brite"    "module"   "ko"       "genome"   "vg"       "ag"       "compound" "glycan"   "reaction"
[11] "rclass"   "enzyme"   "disease"  "drug"     "dgroup"   "environ"  "genes"    "ligand"   "kegg"

# 获取KEGG数据库中某个物种的所有通路（如人类）
> keggList("pathway","hsa")

# 获取某一条KEGG通路的全部信息，部分结果如图1所示。
> keggGet("hsa00020")
# KEGG通路的基因、代谢物等信息就包含在keggGet()函数获得的结果中


图1 通路hsa00020的部分信息

一、提取某条KEGG通路中的所有基因
# 获取某一条KEGG通路的全部信息
path <- keggGet("hsa00010")
head(path)
[[1]]
[[1]]$ENTRY
   Pathway
"hsa00010"
……

# 获取通路中的基因信息
gene.info <- path[[1]]$GENE
gene.info
  [1] "3101"                                                                                         
  [2] "HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]"                                                   
  [3] "3098"                                                                                         
  [4] "HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]"                                                   
……

# 提取gene symbol和Entrez ID
genes <- unlist(lapply(gene.info,function(x) strsplit(x,";")))
gene.symbol <- genes[1:length(genes)%%3 == 2]
gene.id <- genes[1:length(genes)%%3 == 1]
gene.symbol
[1] "HK3"     "HK1"     "HK2"     "HKDC1"   "GCK"     "GPI"     "PFKM"    "PFKP"    "PFKL"    "FBP1"    "FBP2"    "ALDOC"  ……
gene.id
[1] "3101"   "3098"   "3099"   "80201"  "2645"   "2821"   "5213"   "5214"   "5211"   "2203"   "8789"   "230"    "226"    "229"   ……

# 生成gene symbol和Entrez ID匹配的数据框
gene.df <- data.frame(gene.symbol = gene.symbol,gene.id = gene.id)
head(gene.df)
  gene.symbol gene.id
1         HK3    3101
2         HK1    3098
3         HK2    3099

# 提取通路ID和通路名称（有需要的自己添加到数据框gene.df中即可）
> path[[1]]$ENTRY
   Pathway
"hsa00010"
> path[[1]]$NAME
[1] "Glycolysis / Gluconeogenesis - Homo sapiens (human)"提取某条KEGG通路的基因（名称和ID）时有一个需要注意的点，就是有的基因可能有Entrez ID，但没有gene symbol，此时提取会出错。
二、提取某条KEGG通路中的所有代谢物
# 获取某一条KEGG通路的全部信息
path <- keggGet("hsa00010")

# 获取通路中的代谢物信息
compound <- path[[1]]$COMPOUND
compound
          C00022
       "Pyruvate"
          C00024
     "Acetyl-CoA"

# 获取代谢物名称和CPD ID
cpd.name <- as.character(compound)
cpd.id <- names(compound)

# 生成代谢物名称和CPD ID匹配的数据框
metabolism.df <- data.frame(cpd.name = cpd.name,cpd.id = cpd.id)
head(metabolism.df)
             cpd.name cpd.id
1            Pyruvate C00022
2          Acetyl-CoA C00024
3           D-Glucose C00031
# 如果需要通路的ID和名称，操作同上三、提取所有人类KEGG通路中的所有代谢物
# BiocManager::install("KEGGREST",force = TRUE)
# devtools::install_github("https://github.com/cran/RbioRXN.git")
# BiocManager::install("fmcsR")
library(KEGGREST)
library(fmcsR)
library(plyr)
library(devtools)
library(RbioRXN)
library(stringr)
library(dplyr)

hsa_pathway <- keggList("pathway","hsa") # 获取KEGG数据库中所有人类通路
hsa_path <- data.frame(hsa_pathway) # 转成数据框,方便后续分析
hsa_path$pathID <- substr(rownames(hsa_path),6,nchar(rownames(hsa_path)[1])) # 提取pathway ID

# 提取通路中的所有代谢物ID及名称
match.df <- vector()
for (i in 1:nrow(hsa_path)) {
  hsa_info <- keggGet(hsa_path[i,"pathID"])
  hsa_compound <- hsa_info[[1]]$COMPOUND
  path_name <- hsa_info[[1]]$NAME
  if(length(hsa_compound)>0)
  {
    cpd <- names(hsa_compound[1])
    cpd_name <- as.character(hsa_compound[1])
    for (j in 1:(length(hsa_compound)-1)) {
      cpd <- paste(cpd,names(hsa_compound)[j+1],sep = ";")
      cpd_name <- paste(cpd_name,as.character(hsa_compound)[j+1],sep = ";")
    }
    match.ver <- c(hsa_path[i,"pathID"],path_name,cpd,cpd_name)
    match.df <- rbind(match.df,match.ver)
  }
  if(length(hsa_compound)==0){
    match.df <- rbind(match.df,c(hsa_path[i,"pathID"],path_name,"",""))
  }
}

rownames(match.df) <- match.df[,1]
colnames(match.df) <- c("Pathway_ID","Pathway_Name","Compound_ID","Metabolism_Name")


图2 match.df数据框部分截图



图3 代谢物背景文件

自写代码进行富集分析时，需要制作富集背景文件。“提取所有人类KEGG通路中的所有代谢物”部分的R代码提供了制作代谢物背景文件的一种方法，生成的数据框中，每一行表示一条KEGG通路，获得的代谢物名称和ID分别存储在一列中。如果是需要如图3所示的格式（每一行为一个代谢物），对循环稍作修改即可。
制作基因富集背景文件同理，写循环完成即可。
2022.04.21 Thursday 更新
最近KEGG富集分析的背景基因更新了，我在重新制作背景文件时，发现有更简单的获取KEGG通路包含的所有基因相关信息。KEGGREST包提供了KEGG数据库的API接口，但实际上KEGG数据库网页上也是直接有API的，之前我倒是知道，但使用不熟练，也用的比较少，没想到今天真香了。直接从网站下载数据比写循环来得更简单。
KEGG API：https://www.kegg.jp/kegg/rest/keggapi.html
关于KEGG API的其他使用方法，网站是给出了示例的，感兴趣的可自行研究。



图4 KEGG API提供的人类通路和基因的信息

/link/pathway/hsa和/link/hsa/pathway中提供的信息是相同的，即所有参与人类KEGG通路的基因，选择其中一个下载即可，这两个文件只提供KEGG通路ID和基因ID，如果想获得通路名称，就需要另一个文件。如果不想获取通路名称，那这个文件已经足够用来进行富集分析了。



图5 KEGG API提供的人类通路ID和名称的对应关系

/list/pathway/hsa文件提供了人类KEGG通路ID和名称的对应关系。
有了这两个文件，再使用R语言合并数据不能更简单了。比起使用R包，效率会更高。
2022.12.08 Thursday
回看自己半年前写的代码，这什么玩意儿？！太辣鸡了……
对原本的代码进行了优化，就不再贴出来了。

郭亚军

请问图3格式那种循环是怎么个修改思路呀，改了几次都运行不成功[捂脸]

广西百色李祥贵

我是将代谢物名称、ID、通路名称和ID各设为四个vector，最外层的for循环是相同的，里面的if语句也是这两种情况。当hsa_compound的长度大于0时，使用for循环分别提取代谢物名称和ID，hsa_compound的长度是多少，通路名和ID就重复几次（写在提取代谢物名称和ID的for循环外）。当hsa_compound的长度等于0时，通路名和ID只重复一次，代谢物名称和ID用NA填充。最后在最外层for循环外面将四个character类型的变量cbind成一个数据框。

追随昨天

非常感谢！我马上去试试！！[爱]

越天之阶

请问下这个接下来要做代谢物的富集分析，该如何做呢[赞]

一土翁

可以参考ALIEz：R语言进行富集分析及画图，原理一样

后天也不差

为啥我的没有出现数据框啊，无事发生

欧古橱卫全屋定制

emmm，这我不知道如何回答

吉生

为什么我的获取KEGG数据库中所有人类通路这一步 网页连接不上 不能接收数据；如果直接从网页下载那个地方也没有下载地方 是直接全部复制粘贴吗[流泪]

阿里先人

Microsoft Edge可以直接右键另存为