4. BLAST基本原理
BLAST 这个名字，我读研那会儿第一次接触，觉得挺神秘。后来在项目中用它来对比转录组数据，跑了几千条序列之后，才真正理解它为什么能成为生物信息学的「瑞士军刀」。说白了，它的核心就三板斧：局部比对、评分矩阵 和 E 值。你把这三点啃透了，BLAST 就基本玩明白了。
4.1 为什么不选全局比对？
很多人一开始会困惑：为什么 BLAST 用局部比对，而不是咱们常用的 Needleman-Wunsch 或 Smith-Waterman？
我个人习惯这么解释：基因组上的同源区域往往是「碎片化」的。两个基因可能只有某个结构域高度保守，其它部分早就漂移得面目全非了。你要是强行做全局比对，那些不匹配的区域会把分数拉得很低，反而找不到真正的同源片段。
这就好比两个人长得很像，但你可能只注意到他们的眼睛一模一样，而不是从头到脚去对比每一个细节。BLAST 找的就是这种「眼睛」，而不是整张脸。
嗯，这里要注意：BLAST 不是直接拿原始序列去扫数据库。它会把查询序列拆成一个个小片段，叫「单词」（word）。比如默认单词长度是 11 个碱基或 3 个氨基酸。然后它用这些单词在数据库里快速索引查找。为什么能这么快？因为它建立了一张哈希表，单词直接映射到数据库位置，查表一次，O(1) 搞定。
# 示例：BLASTN 默认单词长度
# -word_size 参数可调整
blastn -query my_seq.fa -db nt -word_size 11你可能想问：单词长度设长点好，还是短点好？我在项目中测过：做跨物种比对时，核苷酸保守性低，设 7 比设 11 更容易找到远缘同源。但代价就是假阳性会增加，跑得也慢。
4.2 评分矩阵：从 PAM 到 BLOSUM
比对的结果好不好，很大程度取决于你用哪套评分规则。这就引出了评分矩阵——说白了就是一本「氨基酸相似度字典」。你用它给每一次替换、插入、缺失打分。
项目中我最早接触的是 PAM 矩阵。PAM 全称是「Point Accepted Mutation」，通过观察密切相关的蛋白序列统计出来。PAM 的编号越大，代表进化距离越远。比如 PAM250 比 PAM1 允许更多的替换。但我实际用下来，觉得 PAM 矩阵对短序列不太友好，统计样本也偏小。
下表列出了几个常用矩阵的对比，你可以在配置 BLAST 时通过 -matrix 参数指定：
| 矩阵名称 | 适用场景 | 相对保守度 | 典型参数 |
|---|---|---|---|
| BLOSUM62 | 通用蛋白比对 | 中等分割度 62% | 默认 |
| BLOSUM80 | 高度保守的近缘序列 | 较高 80% | -matrix BLOSUM80 |
| BLOSUM45 | 远缘序列（<30% 一致） | 较低 45% | -matrix BLOSUM45 |
| PAM250 | 早期进化模型 | 中等 | -matrix PAM250 |
| PAM30 | 短序列、快速扫描 | 较高 | -matrix PAM30 |
这里有句我踩过的坑：做细菌 16S rRNA 比对时，很多人默认用 BLOSUM62。结果跑出来一堆高分但不相关的 hits。为什么？因为 16S 是 DNA 序列，你却用了蛋白矩阵——氨基酸替换规则对核酸根本没用。记住：BLASTN 用匹配/错配分数，BLASTP 才用 BLOSUM。
4.3 E 值：你的比对有多可信？
我当年刚学会 BLAST 时，最迷惑的就是 E 值。老师讲多少遍「期望值」我都记不住。后来在项目中比对转录因子家族，看到同一个基因的 E 值从 1e-5 到 1e-150，才意识到它的威力。
E 值的统计学定义是：在随机情况下，预期出现一个与当前比对分数相同或更好的匹配的次数。说白了，它回答的是「这个比对结果是不是碰巧得到的？」。
- E 值非常小（比如 1e-50）：几乎不可能是巧合，匹配非常可靠。
- E 值接近 1 或更大：极有可能是噪声。我一般不会用它做进一步分析。
- E 值在 0.001 到 0.05 之间：属于灰色地带。可能是弱同源，也可能只是保守结构域的残差匹配。
你可能会问：E 值是不是越低越好？不完全是。E 值受数据库大小影响。假设你搜索一个超大基因组，数据量是 100 倍，E 值也会膨胀，因为「随机匹配」的机会变多了。我记得有一次比对人类 lncRNA，数据库太大，很多 genuine 的 hits E 值都在 0.01 附近，但通过比对覆盖度和一致性判定，它们确实是真的。
# 控制 E 值阈值
# -evalue 参数默认是 10
# 高严格搜索可以设到 1e-5
blastp -query protein.fa -db swissprot -evalue 1e-5 -outfmt 64.4 三大概念的协作流程
这三者的关系可以用一个流程来缕清。我刚才画了一张图，方便你理解它们如何联手工作：
从上图可以看到，局部比对把序列切成单词，用评分矩阵打分，再通过统计学模型算出 E 值，最后按 E 值排序输出结果。这个流程在 BLAST 算法里是闭环的——你改一个参数，比如单词长度或空位罚分，E 值就会跟着变。所以调 BLAST 时，你琢磨的是怎么让这三者搭配得当。
4.5 常见误区与避坑
最后，分享几个我在实战中遇到的坑，你注意绕过：
- 误区一：E 值 = 0 就代表绝对可靠。 E 值是一个浮点数，受限于浮点精度，显示为 0 的场合其实只是表示极小，不是绝对。如果数据库更新了，或者你缩小搜索范围，同样的比对可能 E 值就不再是 0 了。
- 误区二：只依赖 E 值不查覆盖度。 我遇到过一段 50 bp 的序列，E 值高得离谱，但覆盖度只有 15%。后来发现它只匹配了一个短保守区域，根本不是同源基因。所以只看 E 值不看覆盖度，容易误判。
- 误区三：DNA 搜索也用氨基酸矩阵。 这算是最常见的低级错误。BLASTN 默认用匹配/错配分数，如果你手动指定了 -matrix，参数会被忽略。但如果你用 TBLASTX（翻译后比对），这时矩阵就能起作用了——用之前确认一下程序的类型。
- 误区四：空位罚分设得太大。 我曾经做过一个启动子分析，为了减少 indels 把空位罚分设到 -20。结果大量真实的外显子跳跃被遗漏了。合适的做法是用 BLOSUM62 + 默认的空位罚分，除非你明确知道插入/缺失机制不常见。