【原】kofamscan：基于HMM模型注释KEGG通路 再也不需要考虑数据库过时了

写在前面

KEGG Ortholog （KO）是最早的一批构建基因家族用于自动化基因功能注释的，随着基因组序列越来越多，虽然有加速版BLAST（Diamond）存在，对氨基酸序列库进行序列相似性搜索也很吃力。KEGG一直都没有提供基于HMM的相似性搜索库，现在有了，KEGG提供了Kofam/ KofamScan 进行KEGG KO的自动化注释，工具和库免费，再也不担心下载不到最新KEGG的序列库问题了。

当前KEGG Kofam（2019-09-09版本）估计包含了 23086的pHMM谱， 原核包含了10252的pHMM谱， 真核包含13666个pHMM谱。

阅读原文：

基于HMM模型注释KEGG通路

KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold

• 原文链接：

  https://academic./bioinformatics/article/36/7/2251/5631907

• DOI：

  https:///10.1093/bioinformatics/btz859

• 第一作者 ：

  Takuya Aramaki

• 通讯作者：

  Hiroyuki Ogata

• 其他作者：

  Romain Blanc-Mathieu, Hisashi Endo, Koichi Ohkubo, Minoru Kanehisa, Susumu Goto

单位：

• ¹ Bioinformatics Center, Institute for Chemical Research, Kyoto University, Gokasho, Uji, Kyoto 611-0011, 2Hewlett-Packard Japan Ltd.,

  
  Koto-ku, Tokyo 136-8711 and 3Database Center for Life Science, Research Organization of Information and Systems, Kashiwa, Chiba

  
  277-0871, Japan

  
  1引言

自动基因功能注释是解释基因组数据的重要第一步。京都基因与基因组百科全书（KEGG）是一个广泛使用的参考知识库，它通过将基因与诸如代谢途径和分子网络等生物学知识联系起来，帮助研究基因组功能（Kanehisa和Goto，2000年）。在KEGG中，KEGG Orthology（KO）数据库-手动管理的大量蛋白质家族（即KO家族）数据库-用作通过K编号标识符将基因与其他KEGG资源（例如代谢图）链接的基线参考。当前，KO被分配给KEGG GENES数据库中的12 934 525（48％）蛋白序列（27 173 868蛋白）。现有的三个工具BlastKOALA，GhostKOALA（Kanehisa 等（2016年）和KAAS（Moriya 等人，2007年）可用于将KO分配给蛋白质序列。这些工具使用同源性搜索软件，例如BLAST（Altschul 等，1997）和GHOSTX（Suzuki 等，2014）来搜索针对GENES的氨基酸序列。为了减少多次成对序列比较所需的冗长计算时间，这些工具使用了从GENES中选择的代表性序列来构建目标数据库。在这项研究中，我们建议采用配置文件隐马尔可夫模型（pHMM）来压缩数据库并为相似性评分定义自适应阈值，这些阈值可用于可靠的KO分配。

2 
对于以给定KO注释的GENES中的每组蛋白质序列，我们生成如下的pHMM。首先，具有100％序列同一性聚类截止值的CD-HIT减少了序列集中的序列冗余。接下来，分别使用MAFFT和HMMER / hmmbuild来比对序列并生成pHMM。

如下为每个KO家庭计算特定于家庭的适应性分数阈值。对于给定的KO家族，属于该家族的非冗余序列被随机分为三组。一组中的一组用作正训练数据集，而其余两组中的序列用于生成pHMM。属于其余KO家族的序列用作所考虑的KO的负面训练数据集。使用HMMER / hmmsearch计算正/负训练数据集中的序列与pHMM之间的序列相似性得分。根据正值和负值数据集中序列的两组位分数的分布，我们确定阈值分数T，该阈值分数使F度量最大化（补充数据）。通过替换三组中的积极训练数据集，此过程重复三遍。最后，将KO家族的自适应分数阈值定义为 T的平均值。用作KO家族分配给新序列的标准。

具有自适应得分阈值的HMM数据库被称为KOfam。当使用KEGG 88.2版进行本研究时，KOFam含有20654个pHMM。我们开发了KofamScan软件，并在KofamKOALA Web服务器中使用了该软件，以使用KOfam注释基因，并通过K编号将其与其他KEGG资源链接，以进行多功能功能研究。

3评估与讨论

为了比较KofamScan与BlastKOALA，GhostKOALA和KAAS的性能，我们使用了从GENES记录的6030个基因组中随机选择的40个基因组（20个真核生物和20个原核生物；补充表S1）作为测试查询。该测试集包含对应于16166个不同KO的383202个序列（有KO的143662个序列）。从GENES中，我们删除了属于所选40个测试查询基因组属的所有基因组。然后，使用剩下的带有KO注释的GENES序列，我们为该评估生成了一个测试的KOfam数据库。至于BlastKOALA，GhostKOALA和KAAS，在从测试查询所代表的属中删除基因组后，我们使用了各自的Web服务器中使用的默认目标数据库。

为此评估创建的KOfam数据库包含20 394个pHMM，其中9414由原核生物序列表示。对于构成我们测试集的40个基因组，KofamScan（0.866），BlastKOALA（0.889）和GhostKOALA（0.862）的性能（F- measure）相当，而KAAS的F- measure （0.810）较低（表1）。为了仅使用20个原核基因组作为测试查询进行另一项测试，我们通过排除仅由真核序列组成的pHMM（对于KofamScan）或通过对原核生物使用目标数据库（对于BlastKOALA，GhostKOALA和KAAS）来简化目标数据库。同样，KofamScan的性能（˚F = 0.875）相当于BlastKOALA（0.846）和GhostKOALA（0.886），而KAAS的F值最低（0.786）。

网页工具注释

如果基因数量不多，我们可以使用网页工具注释：

https://www./tools/kofamkoala/

安装

依赖

从这些依赖看，目前只有linux版本，还是离不开linux，好好学习linux吧。

• Linux

• Ruby >= 2.4

• HMMER >= 3.1

  
  -0 GNU Parallel

数据库安装

cd ~/db
mkdir kofamscan
cd kofamscan

下载数据库，这里一共有1G多一点，还是很小的。

wget ftp://ftp./pub/db/kofam/ko_list.gz # download the ko list
wget ftp://ftp./pub/db/kofam/profiles.tar.gz # download the hmm profiles

从GitHub上：github.com/rotheconrad/KEGGDecoder-binder，下载安装包

解压数据库

gunzip ko_list.gz
tar xf profiles.tar.gz

这里构建虚拟环境安装工具

conda create -n kofamscan hmmer parallel
conda activate kofamscan
conda install -c conda-forge ruby

代码运行

这里我使用命令行运行，不用设置那个配置文件，直接指定数据库位置，方便大家运行。

输入翻译好的氨基酸序列，使用mapper模式匹配数据库；

• p KEGG数据库比对的hmm文件。

• k KO匹配文件。

~/db/kofamscan/kofamscan/exec_annotation -f mapper -o cs.txt ~/Desktop/Shared_Folder/Helper_project/cs.faa -p ~/db/kofamscan/profiles/ -k ~/db/kofamscan/ko_list