BLAST核心步骤：词索引、种子匹配、延伸、打分
BLAST 的核心理念，说白了就是“用空间换时间”。它不是把两条序列逐碱基从头比到尾，而是先建索引、找种子，再精准扩展。我最早接触这个思路时，觉得它有点像查字典——先找部首，再看偏旁，最后读释义。
下面这四个步骤，你理解了它们，也就理解了 BLAST 为什么快，为什么准，又为什么会有漏检。
第一步：词索引（Word Indexing）——构建“单词本”
我先问一个问题：给你 10 万条蛋白质序列，让你快速找到和一条查询序列相似的片段，你会怎么干？
逐条对齐肯定太慢了。BLAST 的做法是：先给查询序列做一套“单词碎片”。
具体来说，对于蛋白质序列，BLAST 把查询序列切成一个个长度为 3 的短肽片段（称为 word）。比如：
查询序列（蛋白质）: MEPALV...
切成 word（窗口大小 = 3）:
MEP
EPA
PAL
ALV
...然后，BLAST 会为这些 word 建一个哈希表（说白了就是一张“单词→位置”的对照表）。
关键参数：
- w（word size）：蛋白质默认 3，核苷酸默认 11。
- word 越短，索引越密，敏感度越高，但内存占用也越大。
- 我习惯在探索性分析中把 w 调小 1，来捕获更多远源同源。
嗯，这里要注意：不是所有 word 都会被保留。BLAST 只会留下那些分数高于某个阈值的 word。这个阈值用 T（neighborhood word threshold） 来控制。
个人经验：
曾经我在做宏基因组注释时，发现默认的 T 参数总是漏掉一些低丰度物种的片段。后来我把 T 从默认的 11 降到了 9，召回率提升了 12%。代价就是索引建的时间多了一倍，但结果质量明显好了。
第二步：种子匹配（Seed Matching）——在索引里“捞”候选
索引建好了，接下来干什么？
很简单：拿着查询序列的 word，去数据库里找完全一样的 word。
比如我们有一个 word 叫“PAL”。数据库索引里有这同一个三肽，出现在第 2001、4503、7890 位。好，这三个位置就叫做种子（seed）。
查询 word: PAL
数据库命中:
db|seqA| 位置 2001 -> 命中
db|seqB| 位置 4503 -> 命中
db|seqC| 位置 7890 -> 命中这个过程是 O(1) 的——直接从哈希表里查。你想想看，如果每次比对都要从头对齐，那就是 O(m × n)。这里差距有多大？一个序列有几万个位置，差距是 4 到 5 个数量级。
| 算法 | 时间复杂度 | 典型查询时间（1000 aa 序列） |
|---|---|---|
| Smith-Waterman | O(m × n) | 几秒到几十秒 |
| BLAST（索引+种子） | O(m × hit_rate) | 毫秒级 |
避坑指南：
我曾经在构建宿主-病原互作数据库时，误以为种子匹配的数量越多越好。结果内存被撑爆了。后来我明白了：种子匹配只是“候选者池”，不是越多越好，而是要精准。通过控制 word 长度和阈值，可以在捕获率和计算量之间找到平衡。
第三步：延伸（Extension）——从种子向两侧生长
种子只是一个点。真正的比对结果要从这一点向两侧扩展。
BLAST 的做法是：以种子对齐位置为起点，同时向左和向右延伸比对。每延伸一个残基，计算一次分数。只要分数在下降不超过一个阈值（通常用 X（dropoff阈值）），就继续延伸。
种子位点：PAL (查询) ↔ PAL (数据库)
向左延伸：
E → E 分数 +1
M → M 分数 +2
向右延伸：
V → V 分数 +1
R → K 分数 +0 （BLOSUM62 打分）
...为什么要有 dropoff 机制？因为有些延伸方向分数越走越低，比如出现空位或者保守性差的区域。如果一直延下去，时间和噪音都受不了。设定 X 值，就是在告诉算法：“分数掉太多就停，别浪费资源。”
我个人习惯的参数调整：
- X（dropoff）：蛋白质默认 7 ~ 20。我一般先跑默认，再看 GAP 区域是否太多。
- 如果目标是找非常保守的域，我会把 X 设小一点（比如 5），结果会更集中。
记得有一次我处理一个含大量重复区域的序列，延伸阶段几乎把整个序列都连起来了。后来调大了 X 并增加了过滤重复的步骤，才摆脱了那个假阳性。嗯，重复序列真的是 BLAST 的大敌。
第四步：打分（Scoring）——质量如何？
延伸到某个长度后，就要给这个局部比对打个分。
BLAST 用的是 局部比对分数，基于 BLOSUM 或 PAM 矩阵。具体来说：
- 每匹配一个残基对，从打分矩阵（如 BLOSUM62）查到一个分值。
- 累加所有分值，得到 raw score（原始分）。
- 再用 raw score 和序列长度、数据库大小，算出 E-value（期望值）。
BLOSUM62 打分示例（部分）：
查询: M E P A L
数据库: M E P A V
得分: +5 +1 +7 +2 +1 = 16 (raw score)
E-value: 0.003 （说明随机出现的概率只有千分之三）打分矩阵到底选哪个？
我一般遵循一条简单规则：
- 如果两条序列的相似度 ≥ 50%（高相似），用 BLOSUM80。它能更好地区分高度保守的位点。
- 如果相似度只有 20%～30%（低相似，远源），用 BLOSUM45。它阈值更低，能扒出更弱的信号。
- 中间情况（40%～50%），BLOSUM62 是通吃型。
四个步骤的完整流程
好，我们把四个步骤串起来看一张图。
图里可以看到，整个流程是线性的。但实际 BLAST 做了很多并行：所有种子同时延伸，互不干扰。这也是它能跑得飞快的原因之一。
额外的优化：两轮命中策略（Two-Hit）
早期的 BLAST 版本要求一个 seed 就能触发延伸。但后来发现，孤立的一个种子容易产生假阳性。于是有了 two-hit 方法：
- 要求两条对角线上有两个种子，距离在一定范围内（比如 40 个残基以内），才启动延伸。
- 这样能大幅减少随机匹配带来的假阳性，尤其在基因组规模的数据集上效果显著。
实用建议：
如果你在用 BLASTN（核苷酸 vs 核苷酸）、BLASTP（蛋白质 vs 蛋白质）、PSI-BLAST（迭代搜索）等不同模式时感觉速度慢，可以先检查种子密度：
# 查看 BLAST 输出中的统计信息（以 BLASTP 为例）
# 注意下面两行（通常在输出尾部）
Query: lcl|12345
Length = 300
# Hits 里的 Number of words: 9850
# 如果 word 数太多（远超过序列残基数 × k），说明索引太密。
# 试着把 word size 增大 1，同时调高 T。一个小技巧：
个人习惯在正式跑全库搜索前，先抽一小部分序列做参数扫描。用一个小 shell 脚本跑不同 w、T、X 组合，对比 E-value 分布和运行时间。5 分钟的预调试能省下后面几小时的重跑。
总结
这四个步骤环环相扣：
- 词索引——把序列切成小块，建好索引。
- 种子匹配——用哈希表快速定位完全匹配。
- 延伸——从种子处向两侧扩展对齐区域。
- 打分——用 BLOSUM/PAM 矩阵评估比对质量，算 E-value 显著性。
你可能会问：为什么用 word 而不用整条序列？原因很简单：局部同源性往往只需要一个短的核心片段就能锚定。就像两个人在茫茫人海中认亲，不用看完全身，看一眼脸就行了。BLAST 就是那个“看一眼脸”的算法，而 Smith-Waterman 更像是“全身扫描”。
明白这四步，你对 BLAST 的理解就从“黑盒子”变成了“白盒子”。后面我们可以接着聊参数怎么选、结果怎么解读、不同 BLAST 变体（blastp、blastn、blastx……）之间的细微差别。