如何获取KEGG数据库中与代谢相关的全部基因呢？

曾老板发来了一篇文献，里面有一个分析方法很有意思，就是对COVID-19病人的PBMC数据仅使用代谢相关的基因做降维聚类分群，居然也能很好的对细胞亚群进行分类，文章里面使用的代谢相关基因来自KEGG数据库，共有1387个基因。文献标题《Integrated analysis of plasma and single immune cells uncovers metabolic changes in individuals with COVID-19》，于2021年9月6号发表在nature biotechnology 杂志上面。

左边就是基于代谢相关基因做的UMAP，是不是很神奇：

但是呢，作者没有提供这个基因列表的表格，那就自己去KEGG搞出来！

简单搜了一下方法，发现 KEGGREST 这个R包可以。

先来看看一个通路的情况

KEGGREST 这个R包可以获取KEGG数据库上面的很多信息，如下面：

rm(list=ls())
library(clusterProfiler) #功能富集分析包，统计学原理累计超几何分布
library(org.Hs.eg.db) # 数据包，做基因ID转换，是物种特异性包Hs表示人
library(GSEABase) # 读取gmt格式数据
library(ggplot2)
library(tidyverse)
library(KEGGREST)

# 查看KEGG包含的子数据库
listDatabases()

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

# 获取某一条KEGG通路的全部信息
keggGet("hsa00020")

# 看下细节
path <- keggGet("hsa00020")
path[[1]]$ENTRY
path[[1]]$NAME
path[[1]]$DESCRIPTION
path[[1]]$CLASS
path[[1]]$PATHWAY_MAP
path[[1]]$GENE
head(path[[1]]$GENE)

这里面最重要的就是 CLASS 和 GENE：

获取所有通路的信息

有了上面的信息之后，现在就可以直接对所有通论走一个循环提取每个通路里面的信息，然后整理为表格：

# 获取KEGG数据库中某个物种的所有通路（如人类）
path_all <- keggList("pathway","hsa")
res <- list()
for(i in 1:length(path_all) ) {
# i <- 1
# 获取某一条KEGG通路的全部信息
# path <- keggGet("hsa00630")
  id <- names(path_all[i])
  path <- keggGet(id)

  ENTRY = if_else(is.na(as.character(path[[1]]$ENTRY)), "NA", as.character(path[[1]]$ENTRY))
  NAME = if_else(is.na(as.character(path[[1]]$NAME)), "NA", as.character(path[[1]]$NAME))
#DESCRIPTION = ifelse(is.null(path[[1]]$DESCRIPTION), "NA", as.character(path[[1]]$DESCRIPTION))
  CLASS = ifelse(is.null((path[[1]]$CLASS)), "NA", as.character(path[[1]]$CLASS))
  GENE = as.character(path[[1]]$GENE)
if(length(GENE)!=0){
    GENE_ID <- GENE[seq(1,length(GENE),by=2)]
    GENE_SYMBOL <- GENE[seq(2,length(GENE),by=2)]
    GENE_SYMBOL <- str_split(GENE_SYMBOL,";",n=2,simplify = T)[,1]
    GENE_SYMBOL <- str_trim(GENE_SYMBOL,side = "both")
  } else{
    GENE_ID="NA"
    GENE_SYMBOL="NA"
  }

# 合并
  res[[id]] <- data.frame(ENTRY=ENTRY, 
                          NAME=NAME,
                          #DESCRIPTION= DESCRIPTION,
                          CLASS=CLASS,
                          GENE_ID=GENE_ID,
                          GENE_SYMBOL=GENE_SYMBOL
                          )
print(id)
}

提取代谢相关通路信息

上面得到的res是一个list，转为df：

library(data.table)
df <- rbindlist(res, fill = TRUE)
save(df,file = "df.Rdata")

# 提取代谢相关通路信息
temp <- as.data.frame(table(df$CLASS))
temp[grepl("Metabolism",temp$Var1),]

df_metabolist <- df[grepl("Metabolism",df$CLASS), ]
gene_metabolist <- unique(df_metabolist$GENE_SYMBOL)
gene_metabolist

这个里面的 GENE_SYMBOL 有KEGG数据库自己的基因命名：

所以就用上面提取的ID，然后转为symbol吧：

genelist <- bitr(gene=unique(df_metabolist$GENE_ID), fromType="ENTREZID", toType="SYMBOL", OrgDb='org.Hs.eg.db')
head(genelist)
gene_metabolist <- unique(genelist$SYMBOL)
length(gene_metabolist)
head(gene_metabolist)

共有1712个，后面再去掉免疫细胞相关marker，应该可以跟文献对应上吧！

是不是很简单！

如果上面的内容对你有用，欢迎一键三连~

转发：