文章阅读目录大纲
在前面写了一篇文章来介绍我们可以如何通过KEGG的BHR评分来注释直系同源。在KEGG数据库的同源注释算法中，BHR的核心思想是“双向最佳命中”。它比简单的单向BLAST搜索（例如，只看你的基因A在数据库里的最佳匹配是基因B）更为严格和可靠。在基因注释中，这种方法可以有效减少因基因家族扩张、结构域保守等原因导致的假阳性注释，从而更准确地识别直系同源基因，而直系同源基因通常具有相同的功能。在今天重新翻看了下KAAS的帮助文档之后，发现KAAS系统中更新了下面的Assignment score计算公式：
We define a score for each ortholog group in order to assign the best fitting K numbers to the query gene:
where Sh is the highest score among all ortholog candidates in the ortholog group, m and n are the sequence lengths of the query and the target of BLAST, respectively, N is the number of organisms in ortholog group, x is the number of organisms in the original ortholog group from which this group is derived, and p is the ratio of the size of the original ortholog group versus the size of the entire GENES database. The second term is for the normalization of the first term by sequence lengths, and the third term is a weighting factor to consider the number of ortholog candidates that are found in the original.
通过分析上面的KO得分计算公式我们可以发现，KEGG的Assignment Score（S_KO）算法设计得非常精妙，它不仅仅是一个简单的相似性得分，而是融合了信息论和统计推断的复合型评分系统。在信息论层面，它是一个综合评分，平衡了证据强度、搜索空间复杂度和事件罕见度，最终选择信息量最大的解释。在生物学层面，它超越了对“最佳匹配”的依赖，通过考察与整个KO直系同源群的整体匹配模式，来识别最特异、最可靠的“家族归属”，从而极大地提升了基因功能注释的精确度和生物学意义。
下面我们将详细的从信息熵原理和生物学意义两个层面来理解KEGG官方的KAAS注释系统中的KO编号的注释评分算法原理。
一、 信息熵原理
信息熵的核心思想是量化“不确定性”或“惊喜度”。一个事件发生的概率越小，它发生时所包含的信息量就越大。S_KO 公式巧妙地运用了对数（log_2），这正是信息论中度量信息的基本单位。现在我们再次来审视上面的这个计算公式：
如上所述，一个查询基因可能与多个不同的KO组都存在显著的BHR关系。Assignment Score的目的就是为每个候选KO打分，从而选出“最合适”的那一个。我们可以将上面的计算公式大致的给这个分数的计算综合考虑了三个因素：
- 匹配的质量（BHR分数本身）。
- 匹配的特异性（一个基因匹配了多个KO，或者一个KO里有多个基因匹配了它，都会降低特异性）。
- 匹配的统计学意义（这种匹配模式是偶然发生的还是真实存在的？）。
如上所述，大致的，依照上面提到的三个因素，我们在这里将这个计算函数划分为三大部分。
1. 第一项：S_h (比特分数) - 核心证据
信息论视角：BLAST的比特分数 S_h 本身就是一个信息量度量。它是一个对数几率分数，表示这个比对结果比随机比对更有可能发生的程度。S_h 越高，代表“查询基因与这个KO中的某个基因高度相似”这个事件所包含的信息量越大，是判断同源性的核心证据。
2. 第二项：-log_2(mn) (序列长度归一化) - 对搜索空间的惩罚
从信息论的视角来看，m 和 n 分别是查询基因和参考基因的长度。mn 近似代表了所有可能比对方式的总数，即搜索空间的大小。log_2(mn) 度量的是这个搜索空间的复杂度。为什么是惩罚？：我们是在 S_h 的基础上减去 log_2(mn)。这可以理解为：一个高分比对如果是在一个巨大的搜索空间（非常长的序列）中得到的，它的“惊喜度”会打折扣。因为序列越长，随机出现高相似片段的可能性就越大。这个项有效地对序列长度进行了归一化，消除了长度偏向性，使得不同长度基因的得分更具可比性。它就像在问：“考虑到所有可能的排列组合，你得到的这个高分还那么令人惊讶吗？”
3. 第三项：-log_2(P(k >= N)) (统计特异性奖励) - 算法的精髓
这是整个公式最巧妙的部分，也是最容易误解的地方。它实际上是一个奖励项，而不是惩罚项。我们从信息论的视角来看，会发现：
- 括号内的部分 P(k >= N) = sum_{k=N}^{x} ... 是一个概率值。它计算的是：如果“一个基因随机通过BHR阈值”的概率是 p，那么在一个包含 x 个基因的KO中，我们观察到 N 个或更多基因通过BHR的概率是多少？这是一个典型的二项分布右尾检验。
- log_2(P(k >= N)) 是对这个概率的度量。概率越小，这个值就越负。
- 关键点：公式是 S_h - ... - log_2(P)。当 P(k >= N) 非常小（即事件非常罕见）时，log_2(P) 是一个很大的负数。减去一个很大的负数，等于加上一个很大的正数！
- 因此，-log_2(P(k >= N)) 这一项，是在奖励统计上的特异性。一个事件越罕见（P越小），它提供的信息量（惊喜度）就越大，对总分的贡献就越多。
到这里，我们对上面的公式中三部分内容的信息熵观点上做了总结，得到了下面的信息熵原理：S_KO 的计算过程，就像一个理性的信息决策者：
- 获取核心证据 (S_h)：找到最强的相似性信号。
- 扣除背景噪音 (-log_2(mn))：考虑搜索空间的大小，对信号进行“纯化”。
- 奖励“意外发现” (-log_2(P))：如果这个信号不仅强，而且在群体层面（整个KO）上表现出极不寻常的聚集性（N远大于期望值），就给予巨大的加分。
最终，S_KO 最高的KO，就是那个“最不可能由随机因素造成，最具有信息含量的候选者”。
二、 可能的生物学意义
我们将上述信息论原理翻译成生物学语言，就能理解其深刻的生物学洞见。这个算法旨在解决一个核心问题：如何区分真正的功能同源和偶然的结构相似。
1. 从“单个相似”到“群体一致”
传统的BBH或SBH方法，本质上是在找一个“最佳拍档”。但这有风险：一个基因可能因为包含一个保守的结构域（如ATP结合域）而与许多不同功能的基因都有高分比对。S_KO 的第三项 -log_2(P(k >= N)) 将视角从个体提升到了群体。它不再只关心“我的基因和KO里的哪个基因最像？”，而是关心“我的基因和这个KO里的基因群体，整体上看起来像不像一家人？”
2. 识别“家族特征”而非“个人外貌”
- 场景A：一个大的、功能松散的KO。假设KO X 包含1000个基因（x=1000），功能都是“激酶”。我们所查询的基因是一个激酶，它可能和其中5个基因（N=5）有不错的BHR得分。虽然 S_h 可能很高，但由于 x 很大，N 相对较小，观察到 N=5 这个事件的随机概率 P(k >= 5) 可能并不算特别小。因此，S_KO 不会得到很高的奖励。算法在告诉我们：“你匹配的是一个大家族，证据不够特异。”
- 场景B：一个小的、功能特异的KO。假设KO Y 只包含20个基因（x=20），功能是“丝氨酸/苏氨酸蛋白激酶X亚型”。我们的查询基因和其中12个基因（N=12）都有高BHR得分。在一个只有20个成员的家族里，一下找到12个亲戚，这在随机情况下是极小概率事件（P(k >= 12) 非常小）。因此，S_KO 会得到巨大的奖励。算法在告诉我们：“这几乎可以肯定是铁证如山，你属于这个小而精的家族！”
- 场景C：无意义的匹配。如果 N=0 或 N=1，那么 P(k >= N) 会很大，log_2(P) 接近0，奖励项消失。S_KO 就退化成了 S_h - log_2(mn)。
3. 提高注释的特异性和可靠性
通过这种“群体一致性”的检验，S_KO 算法能够：
- 避免过宽泛的注释：防止将一个基因简单地注释为一个大类（如“激酶”），而是努力将其定位到更具体的功能亚类。
- 增加注释的置信度：一个高分 S_KO 意味着该注释不仅在单个序列上证据充分，而且在整个功能群层面上也具有高度的统计一致性，结果更可靠。
- 处理旁系同源体：在复杂的基因家族中，旁系同源体之间既有相似性又有功能分化。S_KO 通过考察与整个KO群体的匹配模式，能更好地将查询基因分配到功能最相近的那个旁系同源群中。
代码复现
基于上面的公式描述，我们可以编写出如下所示的代码来实现对应的得分公式计算：
''' <summary>
''' 计算KEGG的KO分配得分 S_KO
''' </summary>
''' <param name="Sh">最高比特分数</param>
''' <param name="m">查询基因长度</param>
''' <param name="n">参考基因长度</param>
''' <param name="x">KO中的总基因数</param>
''' <param name="N">满足BHR阈值的基因数</param>
''' <param name="p">单个基因满足BHR阈值的经验概率 (p0)</param>
''' <returns>Assignment Score (S_KO)</returns>
''' <remarks>
''' 这个函数直接对应了 KEGG 帮助文档中的公式：
'''
''' We define a score for each ortholog group in order to assign the best fitting K numbers to the query gene:
'''
''' ```
''' S_KO = S_h - log_2(mn) - log_2( sum_{k=N}^{x} C_k * p^k * (1-p)^{x-k} )
''' ```
'''
''' where Sh is the highest score among all ortholog candidates in the ortholog group,
''' m and n are the sequence lengths of the query and the target of BLAST, respectively,
''' N is the number of organisms in ortholog group, x is the number of organisms in the
''' original ortholog group from which this group is derived, and p is the ratio of the
''' size of the original ortholog group versus the size of the entire GENES database.
''' The second term is for the normalization of the first term by sequence lengths, and
''' the third term is a weighting factor to consider the number of ortholog candidates
''' that are found in the original.
''' </remarks>
Public Function CalculateAssignmentScore(Sh As Double,
m As Integer,
ni As Integer,
x As Integer,
N As Integer,
p As Double) As Double
If Sh <= 0 OrElse
m <= 0 OrElse
ni <= 0 OrElse
x <= 0 OrElse
p <= 0 OrElse
p >= 1 Then
Return Double.NegativeInfinity ' 无效输入
End If
' 第一项: S_h
Dim term1 As Double = Sh
' 第二项: log2(mn)
Dim term2 As Double = Log2(m * ni)
' 第三项: log2( P(k >= N) ) = log2( sum_{k=N}^{x} C_k * p^k * (1-p)^(x-k) )
' 为了数值稳定性，我们直接在对数空间计算
' 这一项是一个二项分布的右尾概率的对数，计算复杂且容易数值不稳定。
' 因此，我们将其委托给一个专门的函数来处理。
Dim term3 As Double = Log2BinomialRightTailProbability(x, N, p)
' S_KO = S_h - log2(mn) - log2(...)
Return term1 - term2 - term3
End Function
''' <summary>
''' 计算二项分布右尾概率的对数 (以2为底)
''' 即计算 log2( sum_{k=N}^{x} C(x,k) * p^k * (1-p)^(x-k) )
''' </summary>
''' <param name="x">总试验次数 (对应KO中的总基因数)</param>
''' <param name="N">成功的最小次数 (对应满足BHR阈值的基因数)</param>
''' <param name="p">单次试验成功的概率 (对应经验概率 p0)</param>
''' <returns>右尾概率的以2为底的对数</returns>
Private Function Log2BinomialRightTailProbability(x As Integer, N As Integer, p As Double) As Double
If N > x Then Return Double.NegativeInfinity ' 概率为0
If N <= 0 Then Return 0 ' 概率为1, log2(1)=0
Dim logTerms As New List(Of Double)()
Dim logP As Double = Log2(p)
Dim Log1MinusP As Double = Log2(1 - p)
Dim logTerm As Double
' 步骤1: 在对数空间中计算求和中的每一项
For k As Integer = N To x
' log2(C_k * p^k * (1-p)^(x-k)) = log2(C_k) + k*log2(p) + (x-k)*log2(1-p)
logTerm = Log2BinomialCoefficient(x, k) + k * logP + (x - k) * Log1MinusP
logTerms.Add(logTerm)
Next
' 使用 Log-Sum-Exp 技巧来稳定地计算 log(sum(exp(terms)))
' 这是为了避免直接对极小数求和时发生数值下溢。
' log2(sum(2^term_i)) = max_term + log2(sum(2^(term_i - max_term)))
' 找到最大的 log(term_i)
Dim maxLogTerm As Double = logTerms.Max()
' 通过减去最大值来缩放所有项，使指数部分在0和1之间，防止下溢
Dim sumOfExps As Double = logTerms.Sum(Function(t) 2 ^ (t - maxLogTerm))
' 恢复缩放，得到最终的对数和
Return maxLogTerm + Log2(sumOfExps)
End Function
''' <summary>
''' 计算二项式系数的对数 (以2为底) log2(C(n, k))，
''' 使用 LogGamma 函数来避免计算大数阶乘导致的溢出。
''' </summary>
Private Function Log2BinomialCoefficient(n As Integer, k As Integer) As Double
If k < 0 OrElse k > n Then Return Double.NegativeInfinity
' log2(n! / (k! * (n-k)!)) = (ln(n!) - ln(k!) - ln((n-k)!)) / ln(2)
' 使用 LogGamma 函数计算 ln(n!), 因为 ln(n!) = LogGamma(n+1)
Dim logGammaResult As Double = MathGamma.lngamm(n + 1) - MathGamma.lngamm(k + 1) - MathGamma.lngamm(n - k + 1)
' 将自然对数转换为以2为底的对数
Return logGammaResult / Math.Log(2)
End Function那现在我们有了上面的核心公式计算函数之后，就可以按照上面的函数所计算出来的得分结果对我们的匹配结果做排序，取最高得分的KO匹配结果作为我们最终的KEGG ID的注释结果了：
''' <summary>
''' 为单个查询基因分配最佳KO。它会评估所有可能的KO，并选择得分最高的一个。
''' </summary>
''' <param name="allBHRHitsForQuery">该查询基因与所有参考基因的BHR计算结果列表</param>
''' <param name="koGeneCounts">一个字典，包含每个KO的总基因数</param>
''' <param name="bhrThreshold">BHR阈值</param>
''' <returns>得分最高的KO分配候选，如果没有则返回Nothing</returns>
'''
<Extension>
Public Function AssignBestKO(allBHRHitsForQuery As IEnumerable(Of BHRHit),
koGeneCounts As Dictionary(Of String, Integer),
bhrThreshold As Double) As KOAssignmentCandidate
Dim totalGenesInDatabase As Integer = koGeneCounts.Values.Sum()
' 1. 将该查询基因的所有比对结果，按其所属的KO进行分组
Dim hitsByKO = From hit In allBHRHitsForQuery
Where Not String.IsNullOrEmpty(hit.HitName) ' 确保hit有KO信息
Group hit By ko = hit.HitName Into Group
' 2. 为每个KO组计算其S_KO得分，生成一个候选列表
Dim candidates As KOAssignmentCandidate() = (
From hit_group
In hitsByKO
Let ko As String = hit_group.ko
Where koGeneCounts.ContainsKey(ko)
Let groupHits = hit_group.Group.ToArray
Let koSize As Integer = koGeneCounts(ko)
Let p As Double = koSize / totalGenesInDatabase
Select groupHits.KOCandidates(ko, x:=koSize, bhrThreshold, empiricalProbability:=p)
).ToArray
' 3. 选择得分最高的KO
If Not candidates.Any() Then
Return Nothing
End If
Return candidates.OrderByDescending(Function(c) c.AssignmentScore).First()
End Function
''' <summary>
'''
''' </summary>
''' <param name="groupHits"></param>
''' <param name="ko"></param>
''' <param name="x">获取KO的总基因数 x</param>
''' <param name="bhrThreshold"></param>
''' <param name="empiricalProbability"></param>
''' <returns></returns>
<Extension>
Private Function KOCandidates(groupHits As BestHit(), ko As String,
x As Integer,
bhrThreshold As Double,
empiricalProbability As Double) As KOAssignmentCandidate
' 找到该KO组中得分最高的比对，用于获取 S_h, m, n
Dim topHit = groupHits.OrderByDescending(Function(h) h.score).First()
Dim Sh As Double = topHit.score
Dim m As Integer = topHit.query_length
Dim ni As Integer = topHit.hit_length ' n
' 计算满足BHR阈值的基因数 N
Dim N As Integer = Aggregate h As BHRHit
In groupHits
Where h.BHRScore >= bhrThreshold
Into Count
' 计算S_KO
Dim score As Double = CalculateAssignmentScore(Sh, m, ni, x, N, empiricalProbability)
Return New KOAssignmentCandidate(topHit) With {
.QueryName = topHit.QueryName,
.KO = ko,
.AssignmentScore = score,
.Sh = Sh,
.N = N,
.X = x,
.empiricalProbability = empiricalProbability
}
End Function在上面的函数代码中，有两个参数需要我们做额外的统计计算工作：
KO 基因数量统计
koGeneCounts As Dictionary(Of String, Integer)
''' <summary>
''' 基于reverse比对计算出每一个ko中的基因总数
''' </summary>
''' <param name="reverse">query is the KO, use input target as subject reference</param>
''' <returns></returns>
<Extension>
Public Function KOgeneCounts(reverse As IEnumerable(Of BestHit)) As Dictionary(Of String, Integer)
Return reverse.GroupBy(Function(a) v228.FirstToken(a.QueryName)) _
.ToDictionary(Function(a) a.Key,
Function(a)
Return a.Count
End Function)
End FunctionKO分布的先验概率
empiricalProbability As Double
- p 的真正含义：背景概率
根据KEGG文档和算法的逻辑，p 代表的是一个随机基因恰好属于某个特定KO组的先验概率。这个概率的大小，直接反映了该KO组的普遍程度：
- 一个大的、普遍的KO组（如核糖体蛋白、RNA聚合酶等核心成分），其 p 值会相对较高。
- 一个小的、特异的KO组（如某个物种特有的代谢途径中的酶），其 p 值会非常低。
算法的第三项（统计学惩罚项）的目的，正是要评估“观察到 N 个或更多有效命中”这件事，在给定的背景概率 p 下，是“意料之中”还是“出乎意料”。这个p概率值会需要我们根据匹配结果为每一个KO编号动态计算出来，如果 p 值对所有KO都一样，这个评估就失去了意义。
- p 的正确计算方法
p 的计算公式非常明确：p = 某个特定KO组中的基因总数 / 本地KEGG参考数据库中的基因总数
幸运的是，我们在上面缩写出来的的代码中已经有了 koGeneCounts 字典，这正是我们需要的：我们只需要从匹配结果中提取出对应的KO编号，计算出这个ko编号的基因数量在我们的数据库中的分布频率值即可。如上面的代码所示，在查询中即可以完成相应的参数计算：
From hit_group
In hitsByKO
Let ko As String = hit_group.ko
Where koGeneCounts.ContainsKey(ko)
Let groupHits = hit_group.Group.ToArray
' 从匹配结果中提取出对应的KO编号，拿到这个编号的基因总数
Let koSize As Integer = koGeneCounts(ko)
' 计算出这个ko编号的基因数量在我们的数据库中的分布频率值
' 用作为先验概率
Let p As Double = koSize / totalGenesInDatabase
Select groupHits.KOCandidates(ko, x:=koSize, bhrThreshold, empiricalProbability:=p)
从截图上来看，我们针对某一个需要做注释的蛋白进行了BHR筛选，得到了两个在BHR结果值上一摸一样的注释词条的注释结果。假若光靠BHR评分，这个蛋白的注释结果很可能是一个随机筛选的结果。但是在添加了KAAS Rank排序得分之后，一个注释词条得到了1617.32分，另一个注释词条得到了1616.99分。这两个注释词条以很微弱的差异得到了排序。从而基于KAAAS排序评分让我们的注释结果少了这种非常有可能出现的随机筛选的情况。
结合BHR和Assignment Score为完整KO注释流程
在这里，如果我们想要将前面所提到的BHR匹配方法和在这里的KAAS Assignment Score评分方法组合在一起，构建出一个完整的KEGG数据库所使用的KO编号注释流程的话，我们可以在这里将BHR看作为是一个过滤器，用于从海量的BLAST结果中筛选出高可信度的、双向验证过的同源关系。然后将KAAS Assignment Score 用作为是一个决策器，用于当一个基因通过了BHR过滤器但仍有多个候选KO时，通过综合考虑匹配质量、特异性和统计学意义，来选出唯一的、最佳的KO注释。
在之前的BHR算法代码的基础上，我们作了进一步修改，用来支持BHR过滤操作：
Public Iterator Function BHRGroups(forward As NamedCollection(Of BestHit)(),
reverse As NamedCollection(Of BestHit)(),
Optional threshold# = 0.95,
Optional bitsCutoff As Double = 60) As IEnumerable(Of NamedCollection(Of BHRHit))
Dim Rf As TopHitRates() = forward.GetHitRates(bitsCutoff, group:=True).ToArray
Dim Rr As Dictionary(Of String, NamedCollection(Of NamedValue(Of Double))) = reverse _
.GetHitRates(bitsCutoff, group:=True) _
.ReverseAssembly _
.ToDictionary(Function(hit)
Return hit.name
End Function)
For Each q As TopHitRates In Rf
Dim groupRr As NamedCollection(Of NamedValue(Of Double)) = Rr.TryGetValue(q.queryName)
If Not groupRr.value.IsNullOrEmpty Then
' handling of the possible duplicated key name error at here
Dim RrIndex As Dictionary(Of String, NamedValue(Of Double)) = groupRr _
.GroupBy(Function(a) a.Name) _
.ToDictionary(Function(a) a.Key,
Function(a)
Return a.OrderByDescending(Function(i) i.Value).First
End Function)
Yield New NamedCollection(Of BHRHit)(q.queryName, q.MakeBHRGroup(Rr:=RrIndex, threshold))
End If
Next
End Function
<Extension>
Private Iterator Function MakeBHRGroup(Rf As TopHitRates,
Rr As Dictionary(Of String, NamedValue(Of Double)),
threshold#) As IEnumerable(Of BHRHit)
For Each fscore As NamedValue(Of Double) In Rf.hits
Dim forward As BestHit = Rf.htop(fscore.Name)
Dim rscore As NamedValue(Of Double) = Rr.TryGetValue(fscore.Name)
If rscore.Value = 0 Then
Continue For
End If
Dim bhr As Double = fscore.Value * rscore.Value
If bhr >= threshold Then
Yield New BHRHit With {
.QueryName = Rf.queryName,
.HitName = forward.HitName,
.description = forward.description,
.score = forward.score,
.evalue = forward.evalue,
.identities = forward.identities,
.positive = forward.positive,
.length_hit = forward.length_hit,
.length_query = forward.length_query,
.length_hsp = forward.length_hsp,
.BHRScore = bhr,
.hit_length = forward.hit_length,
.query_length = forward.query_length
}
End If
Next
End Function基于上面所给出来的BHR匹配过滤代码，我们筛选出来很多高匹配的BHR结果。基于这些结果，我们就可以用作为KAAS Assignment Score的计算输入了：
Dim bhrGroups As NamedCollection(Of BHRHit)() = BHR.BHRGroups(forwardHits, reverseHits, threshold, bitsCutoff:=score_cutoff).ToArray
Dim assignments As New List(Of KOAssignmentCandidate)
Dim geneCounts As Dictionary(Of String, Integer) = reversePool.KOgeneCounts
For Each bhr As IGrouping(Of String, BHRHit) In bhrGroups
Dim koSize As Integer = Aggregate koid As String
In (From a As BHRHit In bhr Select a.HitName Distinct)
Into Sum(geneCounts(koid))
Dim assign As KOAssignmentCandidate = bhr.AssignBestKO(geneCounts, threshold)
If Not assign Is Nothing Then
Call assignments.Add(assign)
End If
Next
Return assignments.ToArray注释流程测试
我们在这里采用EC Number编号注释工作来测试上面的工作流程，首先我们的EC Number参考数据库中的fasta序列的标题格式为下面所示(EC Number编号位于标题中的第一个单词的位置)：
>3.1.-.-|BioCAD00000675413
MLTGYRLMWMVVMFDLPVITKAERKAATGFRNALLDVGFQMSQFSVYLRFCTSQAQVDTL
CRQVEQALPAGGKVHIFQFTDKQYERAISFHGRSRQPAQKAPDQFDLF
>4.2.1.129|BioCAD00000693803
MAEQLVEAPAYARTLDRAVEYLLSCQKDEGYWWGPLLSNVTMEAEYVLLCHILDRVDRDR
MEKIRRYLLHEQREDGTWALYPGGPPDLDTTIEAYVALKYIGMSRDEEPMQKALRFIQSQ
GGIESSRVFTRMWLALVGEYPWEKVPMVPPEIMFLGKRMPLNIYEFGSWARATVVALSIV
MSRQPVFPLPERARVPELYETDVPPRRRGAKGGGGWIFDALDRALHGYQKLSVHPFRRAA然后基于blastp分别进行query vs subj. reference和subj. reference vs query这样子的正向和反向的比对。得到两个方向上的blastp搜索结果。在GCModeller中提供了相应的R#函数用来完成对应的blastp结果的解析，在这里我们通过下面的一段R#脚本代码完成这样子的forward和reversed方向上的结果解析：
require(GCModeller);
imports "annotation.workflow" from "seqtoolkit";
let export_sbh = function(file, save) {
file
|> read.blast(type = "prot", fastMode = TRUE)
|> blasthit.sbh(topBest = FALSE,keepsRawName = TRUE)
|> stream.flush(stream = open.stream(save, type ="sbh", ioRead=FALSE))
;
}
export_sbh(relative_work("enzyme-ecoli.txt"), relative_work("enzyme-ecoli.csv"));
export_sbh(relative_work("ecoli-enzyme.txt"), relative_work("ecoli-enzyme.csv"));得到如何格式的单向匹配搜索的结果表格：
在上面的所讲解的流程代码中，在GCModeller中导出了下面的一个api函数用于上面所示的KAAS注释流程：
''' <summary>
''' do KO number assign based on the bbh alignment result.
''' </summary>
''' <param name="forward"></param>
''' <param name="reverse"></param>
''' <param name="env"></param>
''' <returns></returns>
<ExportAPI("assign_ko")>
Public Function KOannotations(forward As pipeline, reverse As pipeline,
Optional threshold As Double = 0.95,
Optional score_cutoff As Double = 60,
Optional kaas_rank As Boolean = True,
Optional env As Environment = Nothing) As Object
If forward Is Nothing Then
Return RInternal.debug.stop("forward data stream is nothing!", env)
ElseIf reverse Is Nothing Then
Return RInternal.debug.stop("reverse data stream is nothing!", env)
ElseIf Not forward.elementType.raw Is GetType(BestHit) Then
Return RInternal.debug.stop($"forward is invalid data stream type: {forward.elementType.fullName}!", env)
ElseIf Not reverse.elementType.raw Is GetType(BestHit) Then
Return RInternal.debug.stop($"reverse is invalid data stream type: {reverse.elementType.fullName}!", env)
End If
Dim forwardHits = forward.populates(Of BestHit)(env) _
.Where(Function(a) a.score > 0) _
.GroupBy(Function(a) a.QueryName) _
.Select(Function(q) New NamedCollection(Of BestHit)(q.Key, q)) _
.ToArray
Dim reversePool = reverse.populates(Of BestHit)(env).ToArray
Dim reverseHits = reversePool _
.Where(Function(a) a.score > 0) _
.GroupBy(Function(a) If(kaas_rank, a.QueryName, v228.FirstToken(a.QueryName))) _
.Select(Function(q) New NamedCollection(Of BestHit)(q.Key, q)) _
.ToArray
If kaas_rank Then
Dim bhrGroups As NamedCollection(Of BHRHit)() = BHR.BHRGroups(forwardHits, reverseHits, threshold, bitsCutoff:=score_cutoff).ToArray
Dim assignments As New List(Of KOAssignmentCandidate)
Dim geneCounts As Dictionary(Of String, Integer) = reversePool.KOgeneCounts
For Each bhr As IGrouping(Of String, BHRHit) In bhrGroups
Dim koSize As Integer = Aggregate koid As String
In (From a As BHRHit In bhr Select a.HitName Distinct)
Into Sum(geneCounts(koid))
Dim assign As KOAssignmentCandidate = bhr.AssignBestKO(geneCounts, threshold)
If Not assign Is Nothing Then
Call assignments.Add(assign)
End If
Next
Return assignments.ToArray
Else
Dim bhrResult = BHR.BHRResult(forwardHits, reverseHits, threshold, score_cutoff).ToArray
Return bhrResult
End If
End Function在上面的这个导出来的assign_ko
require(GCModeller);
imports "annotation.terms" from "seqtoolkit";
imports "annotation.workflow" from "seqtoolkit";
let forward = read.besthits(relative_work("ecoli-enzyme.csv"));
let reverse = read.besthits(relative_work("enzyme-ecoli.csv"));
let ec_numbers = annotation.terms::assign_ko(forward, reverse, kaas_rank=FALSE);
print(as.data.frame(ec_numbers));
write.csv(ec_numbers, file = relative_work("ec_number.csv"));
按照前面的文章中的描述，在这里我们按照不同的BHR评分范围将匹配结果划分为了不同的注释等级。在生成的结果数据框中，positive列为对应的BHR评分。
那现在，我们将上面测试的注释流程的R#脚本代码中的kaas_rank参数修改为TRUE，则整个流程将会走首先进行BHR评分过滤，然后再在BHR评分过滤的基础上完成KAAS Rank Assignment Score的计算，实现出一个更加灵活严谨的直系同源注释流程。在这里所示的BHR+KAAS Rank的流程将会生成下面所示的结果数据框：
在上面所生成KAAS Rank结果的数据框中，我们将KAAS Rank Assignment Score计算所需要的参数以及BHR筛选用到的得分，添加在了原先的Hits结果表格的末尾。
- CenterStar多序列比对算法 - 2026年1月8日
- 建立KEGG的KO序列数据库 - 2026年1月4日
- 环境微生物群落GEMs建模综述 - 2025年12月29日
No responses yet