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.oup.com/bioinformatics/article/36/7/2251/5631907
DOI:
https://doi.org/10.1093/bioinformatics/btz859
第一作者 :
Takuya Aramaki
通讯作者:
Hiroyuki Ogata
其他作者:
Romain Blanc-Mathieu, Hisashi Endo, Koichi Ohkubo, Minoru Kanehisa, Susumu Goto
单位:
1 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.genome.jp/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.genome.jp/pub/db/kofam/ko_list.gz # download the ko list
wget ftp://ftp.genome.jp/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