BLAST基础：序列比对算法、E-value与比对分数含义

BLAST这个工具，我想只要做过序列分析的同行，应该没有不知道的。但说实话，很多人只是把它当成一个"搜相似序列的黑盒子"——输入一条FASTA，点提交，等结果。我自己刚入行时也这样，直到有一次项目里跑出来的结果怎么看怎么不对劲，才逼着自己把底层的逻辑啃了一遍。

这一章，我就把BLAST最核心的三个东西讲透：它是怎么比对的、E-value到底在算什么、bit score又代表什么。

BLAST的核心思想：以快制胜

BLAST的全称是Basic Local Alignment Search Tool。关键在"Basic"和"Local"这两个词上。它不像Smith-Waterman那样做全局动态规划，而是只找局部高度相似的片段。

BLAST用的策略叫"seed-and-extend"。

1. 建索引：把查询序列切成固定长度的短片段（通常是3个氨基酸或11个核苷酸），这叫"seed"或"word"。
2. 找命中：在目标数据库里快速查找这些seed的精确匹配或近似匹配。
3. 延伸比对：从命中位置向两边延伸，直到得分不再增加。
4. 评估显著性：用统计模型计算这个比对的E-value。

你可能会问，为什么不用动态规划？我的项目经验告诉我，数据库动不动就是几十G的序列，动态规划在这个量级下根本跑不动。BLAST牺牲了精确性，但换来了实际可用的速度。

E-value：不是你以为的那个p-value

学习BLAST时最容易被坑的就是E-value。很多人初看觉得它跟p-value差不多，都是越小越显著。但两者其实差别很大。

E-value的定义：在同一数据库大小的随机搜索中，期望出现的比对数目。

举个例子。E-value = 0.01意味着：如果你拿一条随机序列去搜索这个数据库，预期每100次搜索会有1次碰到同样得分的比对。你想想看，如果数据库越大，随机命中的机会越大，E-value自然会膨胀。

E-value与p-value的关系（当E-value很小时成立）：

E = 1 - exp(-p) ≈ p

当E-value非常小（< 0.001）时，两者数值几乎相等。但E-value依赖数据库大小，而p-value不依赖。

Bit score：归一化的比对质量

E-value虽然直观，但它跟数据库大小绑定。如果你要在不同数据库之间比较比对的显著性，E-value就不太适合了。

Bit score (S') 就是来解决这个问题的。它把原始比对得分做了归一化处理，消除了数据库大小和评分矩阵参数的影响。

公式长这样：

S' = (λ * S - ln K) / ln 2

其中S是原始比对得分，λ和K是数据库相关的统计参数。说白了，bit score就是一个"标准化"后的分数，单位是bits。

我个人判断的标准很简单：

• Bit score > 200：强同源信号，基本可以闭眼认
• 50 - 200：中等置信度，需要结合功能域注释判断
• < 50：很可能是弱匹配，或者只是共享了一个小motif

有意思的是，bit score跟序列长度有一个大致关系。我统计过不少项目数据，大概每100个氨基酸的保守区域能贡献约50 bits的得分。如果一条全长蛋白只有100个残基，跑出80的bit score，那基本意味着整条序列都在保守区域内。

BLAST与Smith-Waterman的对比

我记得第一次做比对算法选型时，组长问我："为什么不用Smith-Waterman？它更准。" 这是个好问题，我把两者放在一起对比过：

我的做法是：先用BLAST做粗筛，把E-value < 1e-5的候选序列挑出来；如果这些候选后面要用来构建系统发育树或做结构建模，再用Smith-Waterman跑一遍做精校。这样既兼顾了速度，又不牺牲关键路径的精度。

命令行实战：跑一次BLASTP

光讲理论不过瘾，我们来看一个实际的命令行例子。假设你有一条蛋白质序列，想搜一个本地数据库：

blastp -query query.faa -db swissprot -out result.txt -outfmt 6 -evalue 1e-5 -num_threads 8

参数说明：

• -outfmt 6：输出为表格格式，便于下游分析
• -evalue 1e-5：阈值设得不算太严
• -num_threads 8：多线程加速

输出结果前几行大致长这样：

seq1    sp|P12345|GENE_HUMAN    85.2    352    52    1    1e-112    412.3
seq2    sp|Q98765|KIN_YEAST     62.5    280    88    2    3e-45     188.7

每列依次是：查询ID、目标ID、一致性百分比、比对长度、错配数、空位数、E-value、bit score。我个人习惯先看bit score，再看E-value。如果bit score超过200而E-value低于1e-50，基本可以放心往下推功能了。

BLAST的算法流程

下面这张SVG图把BLAST的完整流程说清楚了。我建议你对照着上面的步骤一起看。

从图中你能看到，BLAST并不是从零开始做一遍全序列动态规划。它先通过索引锁定候选区域，然后只在候选位置上做精细比对。这就像查字典——先看偏旁部首缩小范围，再逐字核对，而不是从第一页翻到最后一页。

如何在实际项目中用好E-value和Bit score

理论知识讲完了，我分享两个实际中的判断经验。

经验一：不要只看E-value

有一次我做细菌全基因组注释，一条序列的E-value是1e-30，看起来不错。但我扫了一眼bit score，只有45。问题出在哪里？因为模板序列太短，只有30多个氨基酸，统计上容易产生低E-value。后来我把方案调整成：

• 先看bit score是否超过50
• 再检查比对覆盖度是否超过70%
• 最后才看E-value是否低于1e-5

经验二：跨数据库比较时只用bit score

假设你在Swiss-Prot里搜到一个bit score 120的匹配，又在TrEMBL里搜到一个E-value 1e-80的匹配。E-value不能直接比较，因为两个数据库大小不同。但bit score是归一化的，120就是120。这种情况下，我优先相信bit score更高的那个，不管数据库多大。