3. E-value 的计算原理：得分矩阵、空位罚分、Karlin-Altschul 统计

好，这一章我们聊聊E-value到底是怎么算出来的。很多教程直接扔给你一个公式，然后告诉你阈值设多少。我个人觉得，不理解计算原理，阈值设得再漂亮也是瞎蒙。

E-value的计算，说白了就三块基石：得分矩阵、空位罚分和Karlin-Altschul统计。前两个决定了你比对的原始得分有多高，第三个决定了这个得分在随机情况下出现的概率有多大。

3.1 得分矩阵：给碱基或氨基酸匹配打分

你想想看，两条序列比对，凭什么说它们像？就是看对应的位置上，字符是匹配还是错配。得分矩阵就是定义了每种匹配或错配的分数。

对于核酸BLAST，常用的就是匹配+1分、错配-1分或-2分。但蛋白质比对就复杂了。不同的氨基酸替换，生物学意义完全不同。

我早期做项目时，有一次用BLASTP搜一个代谢酶的同源序列。默认用的是BLOSUM62矩阵，结果一个关键功能位点怎么都对不上。后来换成BLOSUM80，一下就找到了。说白了，不同矩阵适用于不同进化距离的序列，千万别死磕一个。

常见的得分矩阵有这么几种：

得分矩阵中的值，你可以直接理解为一种权重。比如BLOSUM62矩阵，对于Ile→Leu的替换，得分是+2（鼓励）；而对于Ile→Asp，得分是-3（重罚）。这些值都是基于大量真实比对的统计得到的。

3.2 空位罚分：插入缺失必须付出代价

进化过程中，序列会发生插入或缺失（indel）。比对时如果允许出现空位，能更好地对齐功能区域。但空位不能乱开，否则任何两条序列都能强行对齐，那就毫无意义了。

空位罚分通常采用仿射罚分模型，意思就是：

• 开启空位罚分（Gap Open Penalty）：较高，通常5~15分。一旦在序列中插入一个空位，先扣掉这个分数。
• 延伸空位罚分（Gap Extension Penalty）：较低，通常1~3分。空位长度每增加一个字符，再扣这个分数。

为什么这样设计？生物学上，单个indel事件往往产生一段连续的缺失，而不是一个字符一个字符地分开缺失。仿射罚分鼓励较少的、连续的indel，而不是一大堆零散的小空位。

我曾经在一个大型宏基因组项目中，用默认的空位罚分跑BLASTX。结果某些低复杂度区域被强行打开十几个空位，E-value虽然很低，但实际是假阳性。后来我把开启罚分从11提高到了15，那些假阳性瞬间消失。嗯，这里要注意，空位罚分太严会漏掉真同源，太松又会引入噪声，需要根据你的序列特征微调。

3.3 Karlin-Altschul统计：把原始得分变成概率

好，到目前为止，你有了得分矩阵和空位罚分，BLAST引擎可以给出一条比对的原始得分S。但原始得分不能直接拿来设阈值。为什么？因为得分高低依赖于序列长度和搜索空间。

想象一下，你在一本1000万字的书里找一句话，和在一封100字的短信里找同一句话，难度完全不一样。Karlin-Altschul统计就是用来消除这种长度依赖的。

它的核心思想是：对于随机序列（无真实同源性），原始得分S的分布服从一种叫做极值分布的特殊形态。公式长这样（不用背，理解就行）：

E = K × m × n × e^(-λ × S)

这里：

• S：原始比对得分（由得分矩阵和空位罚分计算得到）
• m、n：两条序列的长度（数据库维度越大，E-value越容易升高）
• λ（lambda）：归一化因子，主要依赖于得分矩阵和空位罚分参数
• K：缩放常数，同样依赖于打分系统

说白了，λ和K就是Karlin-Altschul统计给你的两个“魔法参数”，它们把不同打分系统下的S值，拉到同一个尺度上，再结合序列长度，算出这个比对在随机情况下发生的期望次数（E-value）。

我的习惯是：用E-value = 0.01做初筛，然后用E-value = 1e-5做严格过滤。为什么是这两个数？因为0.01允许一定的假阳性，适合探索性分析；而1e-5基本只有显著同源才能达到，适合发表论文或做实验验证用。

3.4 一张图看清计算全流程

下面这张图，把得分矩阵、空位罚分和Karlin-Altschul统计的关系串了起来。我建议你多看两遍，比光看公式管用得多。

从图中你能看到，左边的得分矩阵和右边的空位罚分共同决定原始得分S。然后把S代入Karlin-Altschul公式，结合数据库的大小（m、n）和那两个魔法参数（λ、K），最终算出E-value。