KEGG的BHR评分注释直系同源

在基因组学研究中，将新测序的基因或者针对目标基因组进行基于KEGG代谢通路体系的虚拟细胞建模，都会需要将目标基因组与已知功能基因进行比对注释。KEGG（Kyoto Encyclopedia of Genes and Genomes）数据库通过其KEGG Orthology (KO)系统，为基因功能注释提供了一个强大的平台。KO系统将功能上保守的直系同源基因归为一类，每个KO条目（K编号）代表一个直系同源基因群，这些基因在不同物种中通常执行相似的生物学功能。因此，将新基因的序列与KEGG数据库中的已知基因进行比对，可以推断其可能的KO编号，从而将其功能映射到KEGG通路图或功能层级中。

然而，基因组中存在大量基因，如何准确地将每个基因分配到正确的KO组是注释过程的核心挑战之一。直系同源基因（orthologs）是指由共同祖先经物种分化形成的基因，它们在不同物种中往往保留相似的功能。与之相对的是旁系同源基因（paralogs），由基因复制产生，功能可能发生分化。在功能注释中，我们希望优先识别直系同源关系，因为直系同源基因更可能具有保守的功能。为了在大规模基因组注释中自动识别直系同源基因，KEGG开发了Bidirectional Hit Rate (BHR)评分这一指标，用于衡量候选基因与参考基因之间互为最佳匹配的程度。BHR评分在KEGG的注释流程（如KAAS自动注释服务器）中发挥着关键作用，是决定基因功能注释准确性的重要参数。

https://gcmodeller.org/

KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from genomic and molecular-level information. It is a computer representation of the biological system, consisting of molecular building blocks of genes and proteins (genomic information) and chemical substances (chemical information) that are integrated with the knowledge on molecular wiring diagrams of interaction, reaction and relation networks (systems information). It also contains disease and drug information (health information) as perturbations to the biological system.

https://www.kegg.jp/kegg/kegg1a.html

在之前所记录的《记一次浅尝辄止的FBA代谢网络分析》工作之中：对于FBA模型所需要的代谢网络的重建，我们会需要将目标基因组中的蛋白序列集合与从UniProt序列库中抽取出来的具有KO编号的参考序列库做比较。基于KO编号的注释结果得到对应的KEGG代谢反应结果数据的注释，从而重建出整个代谢网络。

关于进行非参考基因组的KEGG代谢途径功能注释的操作，大家也可以观看之前的课程视频《【GCModeller教程】KEGG代谢途径注释原理》

对于酶来说，40-70%的序列相似性对于功能的预测有90%的准确性（Tian,W）。直系同源基于是来自于相同的祖先的基因分化，保存在不同的物种中的功能基因。在实际操作中，他们能够通过BBH（bi-directional best hit）来推测出来。因此，对在许多物种中的直系同源基因的鉴定是进行批量大规模的细胞模型构建，或者对新测序的基因功能预测的最便捷的途径。

蛋白质同源性注释的两个等级

大家在使用KEGG数据库之中的KAAS系统进行蛋白序列的KO编号的注释分析的时候，KAAS系统提供了两种等级的注释方法供大家进行选择：

• single-directional best hit(SBH method)
• bi-directional best hit(BBH method)

简单的说起来，上面的两个注释等级第一个就是只做query vs subject单向注释，取出blastp结果中的每一个query的第一条同源序列作为结果；第二个等级BBH要严格一些：就是除了forward方向的比较，还存在有subject vs query的reverse方向的比较，然后都同时取query结果中的两个方向的第一条序列结果的交集作为注释结果。

上面提到的两个方法非常的简单，naive。但是也会存在一些问题，例如：

• SBH方法会得到很多重复的注释结果
• 而BBH方法则可能会出现注释结果偏少的问题，因为在BBH等级中，序列必须要符合不管是正向比对还是反向比对，双方都必须要处于对方的得分最高的第一条的位置上

BHR评分注释直系同源

BHR（Bidirectional Hit Rate，双向命中率）是一种基于序列比对结果的度量，BHR评分是在BBH的方法上改进而来的，这个评分目前广泛的应用于KO编号的注释之上。其在KEGG数据库中用于评估两个基因在不同物种基因组中互为最佳匹配的程度。其计算基于双向比对（即正反两个方向的BLAST搜索）的得分，具体原理如下：

# Bi-directional hit rate
BHR = Rf * Rr

1. 双向比对搜索：首先，将待注释的基因（查询序列）与KEGG参考基因组数据库中的所有基因进行BLAST比对，记录每个参考基因的比对得分（bit score）。然后，反过来，将每个参考基因作为查询序列，与包含待注释基因的基因组进行BLAST比对，记录待注释基因的比对得分。这样，对于每一对可能的基因对（查询基因A与参考基因B），我们都有两个方向的BLAST得分：A对B的得分和重复B对A的得分。
2. 计算相对得分：由于两个方向的得分可能不同（因为BLAST是启发式算法，且参考基因组大小不同会影响得分），BHR引入了相对得分的概念。具体而言，对于基因A对基因B的比对，我们计算正向相对得分 Rf = S_AB / S_A^max，其中S_AB是A对B的BLAST得分，S_A^max是A在该参考基因组中所有比对中的最高得分（即A在该基因组中的最佳匹配得分）。同理，对于反向比对，计算反向相对得分 Rr = S_BA / S_B^max，其中S_BA是B对A的得分，S_B^max是B在查询基因组中所有比对中的最高得分（即B在该基因组中的最佳匹配得分）。这两个相对得分都在0到1之间，表示各自方向上匹配的相对强度。
3. 计算BHR：BHR定义为正向相对得分与反向相对得分的乘积，即 BHR = Rf × Rr。这个公式同时考虑了两个方向的匹配情况。只有当两个方向的匹配都接近最佳时，BHR才会接近1；如果任一方向匹配不佳，BHR会显著降低。因此，BHR可以看作是对两个基因互为最佳匹配程度的综合度量。
4. 评分范围：由于Rf和Rr均在0~1之间，BHR的取值范围也在0到1之间。BHR=1表示两个基因在彼此基因组中都是唯一的最佳匹配（即互为最佳命中），而BHR接近0则表示至少在一个方向上匹配很差。通常，BHR越高，表示两个基因越可能是直系同源基因；BHR越低，则更可能是旁系同源或无功能关联的基因。

通过上述计算，BHR评分将序列比对结果转化为一个直观的数值，便于后续设定阈值进行筛选。在KEGG注释流程中，BHR常与BLAST得分一起使用，以提高直系同源基因识别的准确性。

理解BHR评分算法其实很容易，就是因为BBH方法因为是取Top1的交集，所以结果会很少。而BHR方法则是认为直系同源的注释的结果不一定为第一条，所以在这里通过一个评分来做阈值筛选完成可能不取第一条的双向比对注释，提升注释结果数量。

Database: ./UniProt_KO.fasta
           8,771,660 sequences; 3,406,738,582 total letters

Query= Synpcc7942_0001 DNA polymerase III, beta subunit
Length=390
                                                            Score     E
Sequences producing significant alignments:                 (Bits)  Value

K02338 Beta sliding clamp                                   786     0.0
K02338 Beta sliding clamp                                   751     0.0
K02338 Beta sliding clamp                                   485     1e-169
K02338 Beta sliding clamp                                   482     1e-168
K02338 Beta sliding clamp                                   480     7e-168
K02338 Beta sliding clamp                                   476     4e-166

现在我们遵从一个比对的顺序定义，就是forward方向为query vs subject，此为Rf；而reverse方向则是为subject vs query，此为Rr。那在这个比对方向的前提基础上，我们可以得到：

• Rf = Bits_score[q_s] / MaxBits_score[q_s]
• Rr = Bits_score[s_q] / MaxBits_score[s_q]
• Blast bits score > 60
• BHR(Rf*Rr) > 0.95

Blast Bits Score 是在 Blast raw score 换算过来的。上面的两个公式中，分子都是 query 和 subject 中的一个基因的Bit_score，分母是 某一个方向上的一个query比对结果中最大的 bit_score。对的，BHR是不是非常简单？

BHR评分的范围与筛选标准

在实际应用中，BHR评分主要用于筛选直系同源候选基因。KEGG注释工具（如KAAS）会先根据BLAST得分选取一定数量的候选同源基因，然后利用BHR进一步过滤，以确保选出的候选基因在两个方向上都是最佳匹配。具体筛选标准如下：

• 设定阈值：通常采用双向最佳命中（BBH）策略时，要求BHR高于一个较高的阈值（例如0.95）。这意味着候选基因对必须在彼此基因组中都接近唯一最佳匹配。BHR>0.95保证了两个基因互为最佳匹配的程度极高，从而大大提高它们是直系同源基因的可信度。
• 双向最佳命中（BBH）：BBH方法是BHR最严格的应用，要求两个基因互为对方的最佳匹配。在BBH策略下，只有当A在B基因组中是最佳匹配且B在A基因组中也是最佳匹配时，才认为它们是直系同源基因。这实际上等同于BHR=1的情况（或非常接近1，考虑到BLAST得分可能有微小差异）。BBH方法非常保守，能最大程度避免将旁系同源误判为直系同源，但可能漏掉一些真实的直系同源基因（例如当存在旁系同源干扰时）。
• 单向最佳命中（SBH）：SBH策略相对宽松，只要求查询基因在参考基因组中找到最佳匹配即可，不严格要求反向匹配。在SBH策略下，通常设定一个较低的BHR阈值（例如只要求正向相对得分Rf>0.95）。这意味着只要查询基因在参考基因组中有一个非常接近的最佳匹配，就将其作为候选，而不管该参考基因在查询基因组中是否也是最佳。SBH方法能提高注释的覆盖率（更多基因获得注释），但可能引入一些假阳性（将旁系同源误认为直系同源）。

Yuki Moriya, Masumi Itoh, Shujiro Okuda, Akiyasu C. Yoshizawa, Minoru Kanehisa, KAAS: an automatic genome annotation and pathway reconstruction server, Nucleic Acids Research, Volume 35, Issue suppl_2, 1 July 2007, Pages W182–W185, https://doi.org/10.1093/nar/gkm321

gkm321.PDF

通过上述标准，BHR评分帮助KEGG注释流程在准确性和覆盖度之间取得平衡。在BBH模式下，高BHR阈值确保注释结果高度可靠；在SBH模式下，适当放宽BHR阈值则可以注释更多基因，适用于需要快速初步注释的场景（如宏基因组数据量巨大时）。

BHR的代码实现

实现上面所描述的BHR评分的计算代码，大家可以阅读GCModeller之中的这个源代码文件：BBH/Algorithm/BHR.vb。

按照KEGG官方文档中的定义，我们可以构建出以下所示的基因同源性注释等级来描述双向比对的结果：

''' <summary>
''' The blast result alignment identify levels
''' </summary>
Public Enum Levels
    ''' <summary>
    ''' A和B之间没有同源性
    ''' </summary>
    NA
    ''' <summary>
    ''' A vs B 以及 B vs A 都是得分最高的(BHR = 1)
    ''' </summary>
    BBH
    ''' <summary>
    ''' A vs B 以及 B vs A 都不是得分最高的,
    ''' 但是BHR得分是最高的 (BHR >= threshold)
    ''' </summary>
    PartialBBH
    ''' <summary>
    ''' 只有一个方向的比对结果是最高的
    ''' </summary>
    SBH
End Enum

那，现在我们在已经解析出blastp的输出结果的基础上，针对某一个方向的blastp比对得到的query基础结果之上，可以有相应的Hit rate计算函数：

''' <summary>
''' Calculate R score for hit
''' ```
''' R = bits / max(bits)
''' ```
''' </summary>
''' <param name="threshold">
''' bit score filter threshold, Those hits with bit scores less than 60 are removed.
''' </param>
''' <returns></returns>
<Extension>
Public Iterator Function HitRate(query As Query, threshold As Double) As IEnumerable(Of NamedValue(Of Double))
    Dim bits = (From hit As SubjectHit
                In query.SubjectHits
                Where hit.Score.Score >= threshold
                Select New NamedValue(Of Double)(hit.Name, hit.Score.Score)).ToArray

    If Not bits.Any Then
        Return
    End If

    Dim maxBits As Double = Aggregate score As NamedValue(Of Double)
                            In bits
                            Into Max(score.Value)

    For Each hit As NamedValue(Of Double) In bits
        Yield New NamedValue(Of Double) With {
            .Name = hit.Name,
            .Value = hit.Value / maxBits
        }
    Next
End Function

从上面的代码中也可以看得出来，计算Rf或者Rr得分其实就是对query中的每一个hits的得分做一次[0,1]正则化，将数据缩放到0到1之间。那再计算出Rf以及Rr之后呢，将两个乘起来，取得分最高的结果即可，下面的代码会将最高得分的匹配结果取出来：

''' <summary>
''' 在这里计算出最高的BHR比对结果
''' </summary>
''' <param name="q"></param>
''' <param name="r"></param>
''' <returns></returns>
<Extension>
Private Function GetTopBHR(q As NamedCollection(Of NamedValue(Of Double)),
                           r As Dictionary(Of String, Double)) As Map(Of (q$, r$), Double)

    Return q.Where(Function(hit) r.ContainsKey(hit.Name)) _
        .Select(Function(hit)
                    Dim bhr_score = hit.Value * r(hit.Name)
                    Dim align = (q.name, hit.Name)

                    Return New Map(Of (String, String), Double) With {
                        .Key = align,
                        .Maps = bhr_score
                    }
                End Function) _
        .OrderByDescending(Function(bhr) bhr.Maps) _
        .FirstOrDefault
End Function

最后，我们就可以基于上面的两个评价函数，按照KEGG官方帮助文档中的方法，从两个基因组数据集之间的双向比对结果中组装出直系同源结果：

<Extension>
Private Function EvaluateBHR(Rf As NamedCollection(Of NamedValue(Of Double)), 
                             Rr As NamedCollection(Of NamedValue(Of Double)), 
                             threshold#) As BiDirectionalBesthit
    If Rr = 0 Then
        Return New BiDirectionalBesthit With {
            .QueryName = Rf.name,
            .length = Rf.description,
            .HitName = HITS_NOT_FOUND,
            .level = Levels.NA
        }
    End If

    Dim topBHR As Map(Of (q$, r$), Double) = Rr _
        .ToDictionary(Function(hit) hit.Name,
                        Function(hit)
                            Return hit.Value
                        End Function) _
        .DoCall(Function(scores)
                    Return Rf.GetTopBHR(scores)
                End Function)

    If topBHR.Maps >= threshold Then
        ' is a BBH
        Return New BiDirectionalBesthit With {
            .level = If(topBHR.Maps = 1.0, Levels.BBH, Levels.PartialBBH),
            .QueryName = Rf.name,
            .length = Rf.description,
            .HitName = topBHR.Key.r
        }
    Else
        Dim maxR = Rf.OrderByDescending(Function(hit) hit.Value).First

        If maxR.Value >= threshold Then
            ' is an acceptable sbh hit
            Return New BiDirectionalBesthit With {
                .level = Levels.SBH,
                .QueryName = Rf.name,
                .length = Rf.description,
                .HitName = maxR.Name
            }
        Else
            ' no hit
            Return New BiDirectionalBesthit With {
                .QueryName = Rf.name,
                .length = Rf.description,
                .HitName = HITS_NOT_FOUND,
                .level = Levels.NA
            }
        End If
    End If
End Function

上面的EvaluateBHR这个函数是最终的决策逻辑，其会根据bhr评分实现不同等级的同源性注释：

1. If topBHR.Maps >= threshold: 首先检查最高BHR得分是否超过阈值。
2. If topBHR.Maps = 1.0: 如果BHR等于1，则判定为BBH。
3. Else: 如果BHR在阈值和1之间，则判定为PartialBBH。
4. Else (当 BHR < threshold): 这个函数会退而求其次，检查前向R值（Rf）中最高的那个是否超过阈值。
   • If maxR.Value >= threshold: 如果是，则判定为SBH。
   • Else: 否则，判定为NA，即没有找到任何有效的同源性结果。

• Author
• Recent Posts

谢桂纲

高级数据科学家 at 苏州帕诺米克

Working on Engineered bacteria CAD design on its genome from scratch. Writing scientific computing software for Tianhe & Sunway TaihuLight supercomputer. Do scientific computing programming in R/R# language, he is also the programming language designer of the R# language on the .NET runtime.

Latest posts by 谢桂纲 (see all)

• KEGG KO编号得分计算 - 2025年12月20日
• 【虚拟细胞】转录事件建模 - 2025年12月20日
• 布隆过滤器在宏基因组测序reads数据分类处理中的作用 - 2025年12月18日