R获取KEGG/GO中某一通路所有基因

1、通过KEGGREST包来探索KEGG数据库的12大代谢通路

#BiocManager安装"KEGGREST"，
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("KEGGREST")

#加载"KEGGREST"
library(KEGGREST) #用于提取通路及基因信息

#获取KEGG数据库收录的所有物种的清单
org <- keggList('organism') 
##小鼠为mmu
#获取小鼠的KEGG数据库的全部通路及基因集
mmu_path <- keggLink("pathway","mmu") 
#得到字符型向量。元素名为基因id，元素为通路名
length(unique(names(mmu_path)))
#[1] 9394 即涉及到的基因数量是 9394个
length(unique(mmu_path))
#[1] 352  即KEGG通路有 352 条

通过观察KEGG官网 ：https://www.genome.jp/kegg/pathway.html 的12大代谢通路，可以看到通路的id都是00开头（大多数，还是有些以01开头）

meta= unique(mmu_path)[grepl('mmu00',unique(mmu_path))]
mmu_info <- lapply(meta, keggGet)  

nm=unlist(lapply( mmu_info , function(x) x[[1]]$NAME))
library(stringr)
genes = unlist(lapply(mmu_info , function(x) {
  g = x[[1]]$GENE
  paste(str_split(g[seq(2,length(g),by=2)],';',simplify = T)[,1],collapse =';')
}))
df =data.frame(
  mmu= meta,
  nm=nm,
  genes =genes
)

1701916503072.png

以上方法，仅能够提取00开头的代谢通路基因，01还需要其他方法

2、找到KEGG特定通路的所有基因

例如，提取牛的TCA，如图，查到KEGG基因号为bta00020

1701916742833.png
1701916806581.png

使用包装函数

getKEGG <- function(ID){
  
  library("KEGGREST")
  
  gsList = list()
  for(xID in ID){
    
    gsInfo = keggGet(xID)[[1]]
    if(!is.null(gsInfo$GENE)){
      geneSetRaw = sapply(strsplit(gsInfo$GENE, ";"), function(x) x[1])   
      xgeneSet = list(geneSetRaw[seq(2, length(geneSetRaw), 2)])          
      NAME = sapply(strsplit(gsInfo$NAME, " - "), function(x) x[1])
      names(xgeneSet) = NAME
      gsList[NAME] = xgeneSet 
    } else{
      cat(" ", xID, "No corresponding gene set in specific database.\n")
    }
  }
  return(gsList)
genelist <- as.data.frame(getKEGG('bta00020'))
write.csv(genelist,"bta00020.csv")

3、R获取GO中某一通路所有基因

#根据物种预处理GO数据库
#需要下载getGoTerm.R脚本https://github.com/JiahaoWongg/Files/blob/main/getGoTerm.R
source("getGoTerm.R")
GO_DATA <- get_GO_data("org.Mm.eg.db", "ALL", "SYMBOL")
save(GO_DATA, file = "GO_DATA_mmu.RData")
#上述处理时间较长，以鼠为例，保存为Rdata后续直接加载

##下面是两个包装函数，建议写为R直接调用
findGO <- function(pattern, method = "key"){
  
  if(!exists("GO_DATA"))
    load("GO_DATA.RData")
  if(method == "key"){
    pathways = cbind(GO_DATA$PATHID2NAME[grep(pattern, GO_DATA$PATHID2NAME)])
  } else if(method == "gene"){
    pathways = cbind(GO_DATA$PATHID2NAME[GO_DATA$EXTID2PATHID[[pattern]]])
  }
  
  colnames(pathways) = "pathway"
  
  if(length(pathways) == 0){
    cat("No results!\n")
  } else{
    return(pathways)
  }
}

getGO <- function(ID){
  
  if(!exists("GO_DATA"))
    load("GO_DATA.RData")
  allNAME = names(GO_DATA$PATHID2EXTID)
  if(ID %in% allNAME){
    geneSet = GO_DATA$PATHID2EXTID[ID]
    names(geneSet) = GO_DATA$PATHID2NAME[ID]
    return(geneSet)     
  } else{
    cat("No results!\n")
  }
}

load("GO_DATA.RData") # 载入数据 GO_DATA
findGO("insulin") # 寻找含有指定关键字的 pathway name 的 pathway
findGO("INS", method = "gene") # 寻找含有指定基因名的 pathway
getGO("GO:1901142") # 获取指定 GO ID 的 gene set