文章阅读目录大纲
在前面写了一篇文章来介绍我们可以如何通过KEGG的BHR评分来注释直系同源。在今天重新翻看了下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），这正是信息论中度量信息的基本单位。现在我们再次来审视上面的这个计算公式：
大致的，我们可以将这这个计算函数划分为三大部分。
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>
''' <param name="empiricalProbability">经验概率 p0</param>
''' <returns>得分最高的KO分配候选，如果没有则返回Nothing</returns>
Public Function AssignBestKO(allBHRHitsForQuery As IEnumerable(Of BestHit),
koGeneCounts As Dictionary(Of String, Integer),
bhrThreshold As Double,
empiricalProbability As Double) As KOAssignmentCandidate
' 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
Let groupHits = hit_group.Group.ToArray
Select groupHits.KOCandidates(ko, koGeneCounts, bhrThreshold, empiricalProbability)
).ToArray
' 3. 选择得分最高的KO
If Not candidates.Any() Then
Return Nothing
End If
Return candidates.OrderByDescending(Function(c) c.AssignmentScore).First()
End Function
<Extension>
Private Function KOCandidates(groupHits As BestHit(), ko As String,
koGeneCounts As Dictionary(Of String, Integer),
bhrThreshold As Double,
empiricalProbability As Double) As KOAssignmentCandidate
' 获取KO的总基因数 x
If Not koGeneCounts.ContainsKey(ko) Then Return Nothing
Dim x As Integer = koGeneCounts(ko)
' 找到该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 BestHit
In groupHits
Where h.SBHScore >= bhrThreshold
Into Count
' 计算S_KO
Dim score As Double = CalculateAssignmentScore(Sh, m, ni, x, N, empiricalProbability)
Return New KOAssignmentCandidate With {
.QueryName = topHit.QueryName,
.KO = ko,
.AssignmentScore = score,
.Sh = Sh,
.N = N,
.X = x
}
End Function
- KEGG KO编号得分计算 - 2025年12月20日
- 【虚拟细胞】转录事件建模 - 2025年12月20日
- 布隆过滤器在宏基因组测序reads数据分类处理中的作用 - 2025年12月18日
No responses yet