#会了它就可以知道哪个基因属于哪个通路，以及不同通路中包含的基因
#参考https://www.jieandze1314.com/post/cnposts/185/
# if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# BiocManager::install("KEGGREST")
library(KEGGREST)
library(tidyverse)
# 1.1先看看能获取到数据库的什么内容
listDatabases()
# > listDatabases()
# [1] "pathway" "brite" "module" "ko" "genome" "vg" "ag" "compound" "glycan"
# [10] "reaction" "rclass" "enzyme" "disease" "drug" "dgroup" "environ" "genes" "ligand"
# [19] "kegg"
#获取全部物种 耗时可根据提共的链接访问
#keggList("organism")
# http://rest.kegg.jp/list/organism
#如果直接用keggList("pathway") 得到的是map的通路
head(keggList("pathway"))
# > head(keggList("pathway"))
# path:map00010 path:map00020
# "Glycolysis / Gluconeogenesis" "Citrate cycle (TCA cycle)"
# path:map00030 path:map00040
# "Pentose phosphate pathway" "Pentose and glucuronate interconversions"
# path:map00051 path:map00052
# "Fructose and mannose metabolism" "Galactose metabolism"
## 总共有多少条map通路
length(keggList("pathway"))
# [1] 551
#如果想获得某物种的全部通路（比如人）
res <- keggList("pathway", "hsa")
length(res)
# [1] 347
# 347条人的通路
# 如果想根据一些关键词挑选通路【比较重要】
# 比如我现在想看人类的通路中和细胞凋亡、分化有关的
options(stringsAsFactors = F)
pat=keggList("pathway", "hsa")
pat=as.data.frame(res)
pat$id=stringr::str_split(rownames(pat),":",simplify = T)[,2]
# 看到我们这里设置了不区分大小写（这里写apoptosis，但识别出来的是Apoptosis）；并且不用写全单词（这里写的是different，但会识别出differentiation）
p=pat[grepl('different|apoptosis',pat[,1],ignore.case=T),]
p
# > p
# res id
# path:hsa04210 Apoptosis - Homo sapiens (human) hsa04210
# path:hsa04215 Apoptosis - multiple species - Homo sapiens (human) hsa04215
# path:hsa04380 Osteoclast differentiation - Homo sapiens (human) hsa04380
# path:hsa04658 Th1 and Th2 cell differentiation - Homo sapiens (human) hsa04658
# path:hsa04659 Th17 cell differentiation - Homo sapiens (human) hsa04659
# 2.3 利用keggGet()提取指定信息
# 来查看一个人类的基因和一个大肠杆菌的基因
query <- keggGet(c("hsa:10458", "ece:Z5100"))
length(query)
names(query[[1]])
# > names(query[[1]])
# [1] "ENTRY" "SYMBOL" "NAME" "ORTHOLOGY" "ORGANISM" "PATHWAY" "NETWORK" "BRITE"
# [9] "POSITION" "MOTIF" "DBLINKS" "STRUCTURE" "AASEQ" "NTSEQ"
# 看看这个基因涉及的通路
query[[1]]$PATHWAY
# > query[[1]]$PATHWAY
# hsa04520 hsa04810
# "Adherens junction" "Regulation of actin cytoskeleton"
# hsa05130 hsa05135
# "Pathogenic Escherichia coli infection" "Yersinia infection"
#和网页一致：https://www.kegg.jp/dbget-bin/www_bget?hsa:10458
# 获得蛋白序列
keggGet(c("hsa:10458", "ece:Z5100"), "aaseq")
# > keggGet(c("hsa:10458", "ece:Z5100"), "aaseq")
# AAStringSet object of length 2:
# width seq names
# [1] 552 MSLSRSEEMHRLTENVYKTIMEQFNPSLRNFIAMGK...VQLKPTVTNDRCDLSAQGPEGREHGDGSARTLAGR hsa:10458 K05627 ...
# [2] 248 MLNGISNAASTLGRQLVGIASRVSSAGGTGFSVAPQ...PARQAPPPPTGPSGLPPLAQALKDHLAAYEQSKKG ece:Z5100 K12786 ...
# 获得DNA序列
keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")
# > keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")
# DNAStringSet object of length 2:
# width seq names
# [1] 1659 ATGTCTCTGTCTCGCTCAGAGGAGATGCACCGGCTC...GGGATGGGAGCGCCCGCACCCTGGCTGGAAGATGA hsa:10458 K05627 ...
# [2] 747 ATGCTTAATGGAATTAGTAACGCTGCTTCTACACTA...ATTTAGCTGCCTATGAGCAATCGAAGAAAGGGTAA ece:Z5100 K12786 ...
# 如果想提取序列可以直接as.character()，例如
as.character(keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")[1])
as.character(keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")[2])
# > as.character(keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")[1])
# hsa:10458 K05627 BAI1-associated protein 2 | (RefSeq) BAIAP2, BAP2, FLAF3, IRSP53, WAML; BAR/IMD domain containing adaptor protein 2 (N)
# "ATGTCTCTGTCTCGCTCAGAGGAGATGCACCGGCTCACGGAAAATGTCTATAAGACCATCATGGAGCAGTTCAACCCTAGCCTCCGGAACTTCATCGCCATGGGGAAGAATTACGAGAAGGCACTGGCAGGTGTGACGTATGCAGCCAAAGGCTACTTTGACGCCCTGGTGAAGATGGGGGAGCTGGCCAGCGAGAGCCAGGGCTCCAAAGAACTCGGAGACGTTCTCTTCCAGATGGCTGAAGTCCACAGGCAGATCCAGAATCAGCTGGAAGAAATGCTGAAGTCTTTTCACAACGAGCTGCTTACGCAGCTGGAGCAGAAGGTGGAGCTGGACTCCAGGTATCTGAGTGCTGCGCTGAAGAAATACCAGACTGAGCAAAGGAGCAAAGGCGACGCCCTGGACAAGTGTCAGGCTGAGCTGAAGAAGCTTCGGAAGAAGAGCCAGGGCAGCAAGAATCCTCAGAAGTACTCGGACAAGGAGCTGCAGTACATCGACGCCATCAGCAACAAGCAGGGCGAGCTGGAGAATTACGTGTCCGACGGCTACAAGACCGCACTGACAGAGGAGCGCAGGCGCTTCTGCTTCCTGGTGGAGAAGCAGTGCGCCGTGGCCAAGAACTCCGCGGCCTACCACTCCAAGGGCAAGGAGCTGCTGGCGCAGAAGCTGCCGCTGTGGCAACAGGCCTGTGCCGACCCCAGCAAGATCCCGGAGCGCGCGGTGCAGCTCATGCAGCAGGTGGCCAGCAACGGCGCCACCCTCCCCAGCGCCCTGTCGGCCTCCAAGTCCAACCTGGTCATTTCCGACCCCATTCCGGGGGCCAAGCCCCTGCCGGTGCCCCCCGAGCTGGCACCGTTCGTGGGGCGGATGTCTGCCCAGGAGAGCACACCCATCATGAACGGCGTCACAGGCCCGGATGGCGAGGACTACAGCCCGTGGGCTGACCGCAAGGCTGCCCAGCCCAAATCCCTGTCTCCTCCGCAGTCTCAGAGCAAGCTCAGCGACTCCTACTCCAACACACTCCCCGTGCGCAAGAGCGTGACCCCAAAAAACAGCTATGCCACCACCGAGAACAAGACTCTGCCTCGCTCGAGCTCCATGGCAGCCGGCCTGGAGCGCAATGGCCGTATGCGGGTGAAGGCCATCTTCTCCCACGCTGCTGGGGACAACAGCACCCTCCTGAGCTTCAAGGAGGGTGACCTCATTACCCTGCTGGTGCCTGAGGCCCGCGATGGCTGGCACTACGGAGAGAGTGAGAAGACCAAGATGCGGGGCTGGTTTCCCTTCTCCTACACCCGGGTCTTGGACAGCGATGGCAGTGACAGGCTGCACATGAGCCTGCAGCAAGGGAAGAGCAGCAGCACGGGCAACCTCCTGGACAAGGACGACCTGGCCATCCCACCCCCCGATTACGGCGCCGCCTCCCGGGCCTTCCCCGCCCAGACGGCCAGCGGCTTCAAGCAGAGGCCCTACAGTGTGGCCGTGCCCGCCTTCTCCCAGGGCCTGGATGACTATGGAGCGCGGTCCATGAGCAGGAATCCCTTTGCCCACGTCCAGCTGAAGCCGACAGTGACCAACGACAGGTGTGATCTGTCCGCCCAAGGGCCAGAAGGCCGGGAGCACGGGGATGGGAGCGCCCGCACCCTGGCTGGAAGATGA"
# > as.character(keggGet(c("hsa:10458", "ece:Z5100"), "ntseq")[2])
# ece:Z5100 K12786 LEE-encoded effector EspF | (GenBank) espF; espF (N)
# "ATGCTTAATGGAATTAGTAACGCTGCTTCTACACTAGGGCGGCAGCTTGTAGGTATCGCAAGTCGAGTGAGCTCTGCGGGGGGAACTGGATTTTCTGTAGCCCCTCAGGCCGTGCGTCTTACTCCGGTGAAAGTTCATTCCCCTTTTTCTCCAGGCTCGTCGAATGTTAATGCGAGAACGATTTTTAATGTGAGCAGCCAGGTGACTTCATTTACTCCCTCTCGTCCGGCACCGCCGCCACCGACAAGTGGACAGGCATCCGGGGCATCCCGACCTTTACCGCCCATTGCACAGGCATTAAAAGAGCACTTGGCTGCCTATGAAAAATCGAAAGGTCCTGAGGCTTTAGGTTTTAAGCCCGCCCGTCAGGCACCGCCGCCACCGACAAGTGGACAGGCATCCGGGGCATCCCGACCTTTACCGCCCATTGCACAGGCATTAAAAGAGCACTTGGCTGCCTATGAAAAATCGAAAGGTCCTGAGGCTTTAGGTTTTAAGCCCGCCCGTCAGGCACCGCCGCCACCGACAAGTGGACAGGCATCCGGGGCATCCCGACCTTTACCGCCCATTGCACAGGCATTAAAAGAGCACTTGGCTGCCTATGAAAAATCGAAAGGTCCTGAGGCTTTAGGTTTTAAGCCCGCCCGTCAGGCACCACCGCCACCGACAGGGCCTAGTGGACTACCGCCCCTTGCACAGGCATTAAAAGATCATTTAGCTGCCTATGAGCAATCGAAGAAAGGGTAA"
#如果想看单个通路的图片
# png <- keggGet("hsa05130", "image")
# t <- tempfile()
# # library(png)
# # writePNG(png, t)
# # if (interactive()) browseURL(t)
# 根据通路找基因【常用】
# 先看单个通路的操作
# 首先得到全部的基因结果（其中会包含Entrez ID和Symbol以及详细的基因信息）
all=keggGet('hsa04210')[[1]]$GENE
# [271] "5595"
# [272] "MAPK3; mitogen-activated protein kinase 3 [KO:K04371] [EC:2.7.11.24]"
entrez=all[seq(1,length(all),2)]
symbol=stringr::str_split(all[seq(2,length(all),2)],";",simplify = T)[,1]
id=data.frame(entrezID=entrez,symbolID=symbol)
# 根据一些关键词挑选了通路，并保持为p
p=pat[grepl('different|apoptosis',pat[,1],ignore.case=T),]
p
# > p
# res id
# path:hsa04210 Apoptosis - Homo sapiens (human) hsa04210
# path:hsa04215 Apoptosis - multiple species - Homo sapiens (human) hsa04215
# path:hsa04380 Osteoclast differentiation - Homo sapiens (human) hsa04380
# path:hsa04658 Th1 and Th2 cell differentiation - Homo sapiens (human) hsa04658
# path:hsa04659 Th17 cell differentiation - Homo sapiens (human) hsa04659
# 之前p是字符串，现在变成列表，用于存储一会得到的各个结果
pp=lapply(p$id, function(x){
# 其实还是像单个处理一样
all=keggGet(x)[[1]]$GENE
entrez=all[seq(1,length(all),2)]
symbol=stringr::str_split(all[seq(2,length(all),2)],";",simplify = T)[,1]
id=data.frame(entrezID=entrez,symbolID=symbol)
})
names(pp)=as.character(p$id)
# 搜索基因信息中的两个词
head(keggFind("genes", c("shiga", "toxin")),3)
#搜索化学公式
head(keggFind("compound", "C7H10O5", "formula"))
# 利用keggConv()转换KEGG ID
keggConv("ncbi-proteinid", c("hsa:10458", "ece:Z5100"))
head(keggLink("pathway", "hsa"))
# 或者指定一个或多个基因（基因支持跨物种）
keggLink("pathway", c("hsa:10458", "ece:Z5100"))
####################总结#########
# 找全部通路
# keggList("pathway")
#
# 根据基因找通路
# keggLink("pathway", c("hsa:10458", "ece:Z5100"))
#
# 根据关键词找通路
# options(stringsAsFactors = F)
# pat=keggList("pathway", "hsa")
# pat=as.data.frame(res)
# pat$id=stringr::str_split(rownames(pat),":",simplify = T)[,2]
# # 看到我们这里设置了不区分大小写（这里写apoptosis，但识别出来的是Apoptosis）；并且不用写全单词（这里写的是different，但会识别出differentiation）
# p=pat[grepl('different|apoptosis',pat[,1],ignore.case=T),]
# 根据通路找基因
# # 这个p可以自己定义，也可以根据上面关键词找通路
# p=pat[grepl('different|apoptosis',pat[,1],ignore.case=T),]
# # 之前p是字符串，现在变成列表，用于存储一会得到的各个结果
# p=as.list(p)
# pp=lapply(p, function(x){
# # 其实还是像单个处理一样
# all=keggGet(x)[[1]]$GENE
# entrez=all[seq(1,length(all),2)]
# symbol=stringr::str_split(all[seq(2,length(all),2)],";",simplify = T)[,1]
# id=data.frame(entrezID=entrez,symbolID=symbol)
# })
# names(pp)=as.character(p)
1 前言
在正式学习这个包之前，还是先补充一点KEGG数据库的知识吧
KEGG数据库是一个综合数据库，其中我们最常用的就是通路（PATHWAY）了，主要介绍了细胞生化过程（还有图注），如细胞代谢、细胞周期、膜转运、信号传递。但PATHWAY只能算上是一个小分支，全部的数据库可以看下面：
1.1 首先，KEGG主要有4大分类
系统信息 Systems information
基因组信息 Genomic information
化学信息 Chemical information
健康信息Health information
后来健康信息+药物标签 Drug labels = KEGG MEDICUS
然后每个分类下面又有自己的子数据库，比如常用的pathway数据库就是在系统信息之下的
1.2 每个分类都有各自的颜色
1.3 每个分支有各自的作用
来自官网的描述：https://www.kegg.jp/kegg/kegg1a.html
KEGG数据库是生物系统的数据重现，由基因、蛋白质（基因组信息）和化学物质（化学信息）的分子模块组建而成。这些分子模块之间构成了相互作用的关系网络（系统信息）。另外还包含了对生命活动有影响的疾病和药物资源（健康信息和药物信息）。
1.4 比较常用的子数据库
1.4.1 KEGG PATHWAY代谢通路数据库
来自：https://www.kegg.jp/kegg/kegg3a.html
包含了许多手工绘制或计算机模拟的代谢通路，包含了细胞过程、遗传信息、新陈代谢、疾病、药物开发等的分子互作和反应网络图。
每个pathway map都是通过2-4个字母的前缀 + 5位数字组成
其中前缀有5种类型：
map - Reference pathway：手工绘制；根据已有的知识绘制的全面的供参考的代谢图，也是为后面几种的基础。图中的小框都是白色（后期可以填充）；每个节点表示某个基因、该基因编码的酶或者酶参与的反应，可以点击其中的链接查看详细信息
org - Organism-specific pathway map：计算机模拟；物种特有的代谢通路；绿色小框表示物种特有的基因或酶，点击绿色框的基因会看到更详细的信息
ko - Reference pathway (KO)：计算机模拟，节点只代表基因；点击蓝色框会看到直系同源基因信息
ec - Reference pathway (EC)：计算机模拟，节点只代表相关的酶
rn - Reference pathway (Reaction)：计算机模拟，节点只代表参与的某个反应、反应物及反应类型
1.4.2 KEGG ORTHOLOGY 直系同源数据库
这个数据库是基于：同源基因功能相似。会把每个功能已知的基因和其他同源的基因归为一类，即一个KO号，它不分物种，比如ko00010（https://www.kegg.jp/kegg-bin/show_pathway?ko00010）就是糖酵解
图中又这么几个形状：
圆角矩形：其他的KO通路
**方框：**如果把鼠标悬停在方框上，就会弹出K+数字。比如这里的K02777、K02779等，表示所有物种的一个同源基因
至于这个同源基因在各个物种中分别是什么，可以看详情，例如K02777（h**小圆圈：**如果把鼠标悬停在方框上，会弹出化合物信息（C符号）ttps://www.kegg.jp/dbget-bin/www_bget?ko:K02777）
