文章阅读目录大纲
KEGG 里面目前并没有“现成的每个 KO 一条代表性序列 FASTA”这种官方序列数据库，假若我们需要基于KEGG数据库中的KO信息的注释，那我们一般会需要自己从 KEGG GENES 里面把每个 KO 对应的基因/蛋白序列抓出来，再按 KO 编号组织成 fasta 集合构建出对应的数据库。基于所建立好的KEGG基因序列数据库，我们就可以实现下面的一些基因注释工作：
- 在全基因组规模代谢网络重建工作中，进行我们的目标基因组中的代谢网络中的酶节点的直系同源推断，从而将我们的目标基因组中的基因映射到具体的KEGG代谢网络上的节点位置，从而重建出代谢网络模型（使用带有KO编号的蛋白序列做比对注释）
- 假若我们在进行宏基因组的基因丰度的计算，则可以基于所建立的KEGG基因序列数据库作为参考库，进行宏基因组测序数据中的KO基因丰度的计算（使用带有KO编号的基因序列做比对注释）
KEGG Orthology介绍
对于KO直系同源（KEGG Orthology）而言，其是一种跨物种的基因功能ortholog 组，即意味着每一个KO编号并不是对应着一条单独的序列，而是把不同物种里功能类似的基因归到一个 K 编号下面。具体而言，KO编号就是一种在KEGG代谢网络中的节点模型概念，而具体的物种中具有对应的KO编号的序列就是该代谢网络节点模型的具体对象实例。即一个 KO（概念） 会对应多个基因（这些基因是来源于不同物种、不同菌株的，并且在同一个基因组内，可能会出现具有相同的KO编号的多个重复的基因），而一个基因序列通常只属于一个 KO。
因此，现在我们假若要构建“每个 KO 对应的基因序列集合”，本质上就是：
- 建一个映射：KO ID → {gene ID1, gene ID2, …}
- 然后对每个 gene ID 抓取它的序列（核苷酸或氨基酸），然后按 KO 归类即可。
从零开始建立KEGG Orthology数据库
KEGG数据库中并没有提供直接的下载链接来进行数据库的下载，但是KEGG数据库提供了相应的REST API编程接口来让我们进行使用。那在这里我们实际上就可以基于此来编写出一些自动化的脚本来实现自动化的建立其KO序列数据库，用于上述的两种场景下的注释计算工作。
今天在这里整理了一下我们如何手动从零开始建立起针对KEGG的KO序列的数据库构建的流程，简单的总结以下，我们大致会需要通过下面的三个步骤来完成数据库的构建：
- 用 KEGG REST API 把 KO 和所属基因 ID 的对应关系拿下来；
- 对每个基因 ID 用 REST API 拿它的核苷酸/蛋白序列；
- 按 KO 号分组，整理成每个 KO 对应的多序列 FASTA（如果嫌多，可以在每个 KO 里只保留一条代表性序列）。
使用到的KEGG REST API介绍
首先，我们需要拿到在KEGG数据库中所有的KO编号的集合，这个可以直接通过ko枚举的api（https://rest.kegg.jp/list/ko）来实现，例如我们在浏览器中访问这个api的url，可以得到：
# https://rest.kegg.jp/list/ko
K00001 E1.1.1.1, adh; alcohol dehydrogenase [EC:1.1.1.1]
K00002 AKR1A1, adh; alcohol dehydrogenase (NADP+) [EC:1.1.1.2]
K00003 hom; homoserine dehydrogenase [EC:1.1.1.3]
K00004 BDH, butB; (R,R)-butanediol dehydrogenase / meso-butanediol dehydrogenase / diacetyl reductase [EC:1.1.1.4 1.1.1.- 1.1.1.303]
K00005 gldA; glycerol dehydrogenase [EC:1.1.1.6]
K00006 GPD1; glycerol-3-phosphate dehydrogenase (NAD+) [EC:1.1.1.8]
K00007 dalD; D-arabinitol 4-dehydrogenase [EC:1.1.1.11]
K00008 SORD, gutB; L-iditol 2-dehydrogenase [EC:1.1.1.14]接着，针对上面的操作得到的每一个KO编号，就可以从kegg数据库中基于基因id的枚举api得到对应关系，例如，我们在浏览器中访问https://rest.kegg.jp/link/genes/K00001，可以得到输出结果：
# https://rest.kegg.jp/link/genes/K00001
ko:K00001 dme:Dmel_CG3481
ko:K00001 der:6540581
ko:K00001 dse:6611292
ko:K00001 dsi:Dsimw501_GD23968
ko:K00001 dya:Dyak_GE10683
ko:K00001 dya:Dyak_GE19037
ko:K00001 dan:6492907
ko:K00001 dsr:110186092通过循环，吧所有的KO编号都处理一遍，这样子我们就可以得到了一个KO编号对应基因编号的数据表。最后，我们就可以通过相应的KEGG数据库中的get工具得到具体的基因序列或者蛋白质序列，例如我们想要针对大肠杆菌的b0001基因拿到对应的序列，而可以访问对应的api：
# gene nucleotide sequence
# https://rest.kegg.jp/get/eco:b0001/ntseq
>eco:b0001 K08278 thr operon leader peptide | (RefSeq) thrL; thr operon leader peptide (N)
atgaaacgcattagcaccaccattaccaccaccatcaccattaccacaggtaacggtgcg
ggctga
# protein sequence
# https://rest.kegg.jp/get/eco:b0001/aaseq
>eco:b0001 K08278 thr operon leader peptide | (RefSeq) thrL; thr operon leader peptide (A)
MKRISTTITTTITITTGNGAG将从API处理得到的序列保存为fasta文件，即可用于下游的数据库分析中建立相应的序列数据库。
GCModeller自动化流程脚本编程示例
在这里我们以构建微生物的GEMs模型所需要的KO注释数据库为例，进行相关的脚本的编写。首先我们需要限定一下fasta序列下载的范围为原核生物内，不然将KEGG数据库中所有的的物种的序列都下载下来，时间成本太高了。我们可以通过下面的GCModeller程序包的脚本编程代码进行KEGG数据库中的原核生物的基因组获取：
require(GCModeller);
require(kegg_api);
imports "dbget" from "kegg_api";
imports "kegg_api" from "kegg_api";
let genomes = kegg_api::listing("organism");
# split text file into multiple columns
# column 1 - kegg organism species code
# column 2 - kegg organism name
# column 3 - taxonomy tree lineage
genomes <- lapply(genomes, str -> strsplit(str, "\t"));
genomes <- data.frame(
row.names = names(genomes),
kegg_code = genomes@{1},
organism = genomes@{2},
kingdom = sapply(genomes@{3},str -> strsplit(str,";")[1]),
lineage = genomes@{3}
);
let bacterial = genomes[genomes$kingdom == "Prokaryotes", ];
print(bacterial, max.print = 6);
# kegg_code organism kingdom lineage
# ----------------------------------------------------------------------------------------------------------------
# <mode> <string> <string> <string> <string>
# T00007 "eco" "Escherichia coli K-12 MG1655" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
# T00068 "ecj" "Escherichia coli K-12 W3110" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
# T00666 "ecd" "Escherichia coli K-12 DH10B" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
# T00913 "ebw" "Escherichia coli K-12 BW2952" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
# T02541 "ecok" "Escherichia coli K-12 MDS42" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
# T09611 "ecoc" "Escherichia coli K-12 C3026" "Prokaryotes" "Prokaryotes;Bacteria;Enterobacteria;Escherichia"
#
# [ reached 'max' / getOption("max.print") -- omitted 10166 rows ]
write.csv(bacterial, file = relative_work("bacterial.csv"));从KEGG数据库上所获取得到的物种列表信息来看，目前在KEGG数据库中，原核生物大概包含有10172个基因组。那后续的fasta序列下载中，基因的物种范围将会被限定在这约1万多个基因组中。从上面所描述的第一步可以知道，我们可以直接通过一个API就可以从KEGG数据库中请求得到数据库中所有记录下的KO编号列表。在下面的代码中，我们进行了相关的API请求，然后将结果解析为一个数据框，方便做数据查看：
# cache table if table file is missing
let index = REnv::getHtml("https://rest.kegg.jp/list/ko", interval = 3, filetype = "txt");
let ko_index = textlines(index)
|> strsplit("(\t+)|(;\s+)")
|> tqdm()
|> lapply(t -> {
# ko_id \t+ gene_name ;\s+ description
# 1 2 3 4 5
# K00001 E1.1.1.1, adh; alcohol dehydrogenase [EC:1.1.1.1]
# K00002 AKR1A1, adh; alcohol dehydrogenase (NADP+) [EC:1.1.1.2]
# K00003 hom; homoserine dehydrogenase [EC:1.1.1.3]
if (length(t) == 3) {
# this KO has no gene name
# due to the gene name is missing, so that
# only has 3 elements, gene name is empty
c(t[1],"",t[3]);
} else {
# pick the 1,3,5
# 1. KO id
# 3. gene name
# 5. gene function text
t[c(1,3,5)];
}
})
;
ko_index = data.frame(
row.names = ko_index@{1},
ko = ko_index@{1},
gene_names = ko_index@{2},
description = ko_index@{3}
);
print("request ko index from KEGG web:");
print(ko_index, max.print = 6);
# [1] "load ko index table for make request data:"
# ko gene_names description
# ---------------------------------------------------------------------------------------------------------------------------------------
# <mode> <string> <string> <string>
# K00001 "K00001" "E1.1.1.1, adh" "alcohol dehydrogenase [EC:1.1.1.1]"
# K00002 "K00002" "AKR1A1, adh" "alcohol dehydrogenase (NADP+) [EC:1.1.1.2]"
# K00003 "K00003" "hom" "homoserine dehydrogenase [EC:1.1.1.3]"
# K00004 "K00004" "BDH, butB" "(R,R)-butanediol dehydrogenase / meso-butanediol dehydrogenase / diacetyl reduc...|36 chars truncated
# K00005 "K00005" "gldA" "glycerol dehydrogenase [EC:1.1.1.6]"
# K00006 "K00006" "GPD1" "glycerol-3-phosphate dehydrogenase (NAD+) [EC:1.1.1.8]"
# K00007 "K00007" "dalD" "D-arabinitol 4-dehydrogenase [EC:1.1.1.11]"
# K00008 "K00008" "SORD, gutB" "L-iditol 2-dehydrogenase [EC:1.1.1.14]"
# K00009 "K00009" "mtlD" "mannitol-1-phosphate 5-dehydrogenase [EC:1.1.1.17]"
# K00010 "K00010" "iolG" "myo-inositol 2-dehydrogenase / D-chiro-inositol 1-dehydrogenase [EC:1.1.1.18 1....|9 chars truncated
# K00011 "K00011" "AKR1B" "aldehyde reductase [EC:1.1.1.21]"
# K00012 "K00012" "UGDH, ugd" "UDPglucose 6-dehydrogenase [EC:1.1.1.22]"
# K00013 "K00013" "hisD" "histidinol dehydrogenase [EC:1.1.1.23]"
#
# [ reached 'max' / getOption("max.print") -- omitted 28015 rows ]
write.csv(ko_index, file = file.allocate("/ko.csv", fs = db));可以看得到，目前在KEGG数据库之中，记录了多于两万8千个直系同源组的保守基因簇的信息。那现在有了KO的id列表后，我们所要做的就是针对每一个KO编号，从KEGG数据库中拿到对应的基因ID列表，然后根据基因ID列表下载对应的fasta序列即可。下面的代码实现了提取某一个KO编号所对应的基因ID列表：
let geneId_file = `/ko/${ko_id}.txt`;
let ko_genes = REnv::getHtml(`https://rest.kegg.jp/link/genes/${ko_id}`, interval = 3, filetype = "txt");
# split each text line with white space
# ko:K00001 dme:Dmel_CG3481
# ko:K00001 der:6540581
# ko:K00001 dse:6611292
# ko:K00001 dsi:Dsimw501_GD23968
# ko:K00001 dya:Dyak_GE10683
#
# first column is the KO id
# use the second column as gene id
ko_genes = ko_genes |> textlines() |> strsplit("\s+");
ko_genes = ko_genes@{2};
writeLines(ko_genes, con = file.allocate(geneId_file, fs = db));在之前我们所提到，我需要下载的是原核生物的基因序列用于构建GEMs，所以在这里我们需要根据前面所得到的bacterial数据框中的kegg_code列中的基因组代码，来做这里得到的ko_genes的过滤，仅保留下原核生物的基因id列表：
# split the kegg gene id by use `:` symbol as delimiter
# eco:b0001
# first column is the kegg organism species code
# second column is the gene locus tag
ko_genes = readLines(file.allocate(geneId_file, fs = db));
ko_genes = strsplit(ko_genes,"[:]");
ko_genes = data.frame(species = ko_genes@{1}, gene_id = ko_genes@{2});
print("genes that belongs to KO group:");
print(ko_genes, max.print = 6);
if (length(species) > 0) {
let i = ko_genes$species in species;
if (sum(i) > 0) {
ko_genes = ko_genes[i,];
} else {
ko_genes = NULL;
}
}最后，基于上面所得到的基因ID列表，就可以创建url进行对应的fasta序列的下载了。
let seqfile = `/fasta/${ko_id}.txt`;
let fasta = NULL;
# split the ko_genes id character vector into multiple parts
# each part has 20 elements
ko_genes = sprintf("%s:%s", ko_genes$species, ko_genes$gene_id);
ko_genes = split(ko_genes, size = 20);
print("download gene sequence with data task blocks:");
print(length(ko_genes));
str(ko_genes);
for(block in ko_genes) {
let url = `https://rest.kegg.jp/get/${paste(block, sep="+")}/${seqtype}`;
let seqs = REnv::getHtml(url, interval = 3, filetype = "txt");
fasta = c(fasta, seqs);
}
writeLines(fasta, con = file.allocate(seqfile, fs = db));在KEGG数据库中，get工具是可以一次性的批量请求多个基因ID的序列信息的。在上面的函数中，我们按照20个id为一组，使用加号连接这些gene id，生成一个url，进行fasta序列文本数据的请求，将所有的分片请求得到的fasta序列文本数据合并到一块后，保存到对应的文件中。我们就可以构建出对应的KO编号的fasta序列数据库了。在上面所建立的url中，有一个seqtype参数，这个参数可以取两个值：ntseq为基因序列，aaseq为蛋白序列。我们从KEGG数据库中获取基因序列，可以用于构建宏基因组测序进行基因注释和基因丰度计算所需要的序列数据库。而使用aaseq获取得到的蛋白序列，则可以能够让我们进行直系同源注释，进行GEM模型的建立。
- 建立KEGG的KO序列数据库 - 2026年1月4日
- 环境微生物群落GEMs建模综述 - 2025年12月29日
- 【虚拟细胞】翻译事件建模 - 2025年12月21日
No responses yet