第4章 BLAST参数详解：e-value、word size、matrix、gap penalties、filtering
很多人用BLAST，上来就是一把梭，输入序列、点击运行、看结果。但你有没有想过，为什么有时候比对出上百条结果，有时候一条都没有？为什么有些比对明显是错的，但BLAST还给了高评分？
说白了，BLAST的底层逻辑就藏在这几个参数里。不懂它们，你看到的图谱就是天书。懂了你就能让BLAST听你的话。
4.1 e-value：我该信你几分？
e-value，全称是“期望值”。这个概念很多人都搞混了。它不是p-value，不是百分比，更不是置信度。它描述的是：在随机情况下，能匹配到这个得分的次数。
举个例子。你设定e-value为1，意味着理论上允许出现1次随机匹配。你设定0.001，那就是一千次随机搜索中才可能出现一次。这个值越小，结果越可靠。
e-value越小，结果越有生物学意义。
e-value越大，假阳性风险越高。
我刚开始做宏基因组分析时，有一次跑BLASTx，e-value我设了10。结果返回了上千条匹配，我还挺高兴，觉得覆盖率高。后来仔细一查，大量结果都是随机匹配的转座子片段，气得我重跑了一周。从那以后，我习惯分场景设定：
| 场景 | 建议 e-value | 理由 |
|---|---|---|
| 严格同源搜索（如基因注释） | 1e-10 以下 | 减少假阳性，确保同源性 |
| 中等严格（如物种鉴定） | 1e-5 ~ 1e-10 | 平衡灵敏度和特异性 |
| 宽松搜索（如发现远源同源） | 0.001 ~ 1 | 捕获更多序列，但注意验证 |
| 极度宽松（如查找污染序列） | 10 ~ 100 | 宁可错杀，不可放过 |
4.2 Word size：搜索粒度的艺术
Word size（词长）是BLAST搜索的基本单元。它决定了算法在比对时的“最小匹配种子”长度。你把word size设成3，只要两条序列有3个连续字符一致，就触发延伸比对。设成11，需要11个字符连续一致才触发。
我个人习惯是这样区分的：短序列（比如30bp的引物）用短word size，长序列（比如完整基因、全长基因组）用长word size。
- 短word size (2~7)：敏感度高，能捕获远源同源，但速度慢、假阳性多。适用于蛋白比对或短序列搜索。
- 长word size (15~30)：速度快、假阳性少，但容易漏掉远源匹配。适用于高度相似的同源物搜索。
我这里有个真实案例。之前做病毒序列分型，用默认的11查找，结果总有一条序列比对不上任何数据库。我把它单独拿出来调低word size到7，一下就匹配到了相关的病毒蛋白。为什么？因为我那条序列突变了几个位点，正好打断了11个连续匹配。word size一降，就能跳过去了。
4.3 Matrix：谁像谁，我说了算
矩阵（Matrix）是BLAST给你的“密码本”。它决定了匹配给多少分，错配扣多少分，不同类型的替换怎么算分。
常用矩阵有三类：
| 矩阵名称 | 适用场景 | 特点 |
|---|---|---|
| BLOSUM62 | 标准蛋白比对，中度进化距离 | 默认，平衡灵敏度和特异性 |
| BLOSUM80 | 高度相似的蛋白 | 对保守区域赋分更高 |
| BLOSUM45 | 远源同源搜索 | 容忍更多替换 |
| PAM250 | 远源进化关系分析 | 基于点突变模型 |
我记得有一次做古老蛋白家族分析，用BLOSUM62跑出一堆低分匹配，怎么看都不像。后来换成BLOSUM45，得分上升了、匹配也多了。说白了，矩阵选不对，好结果可能就被埋了。
对于核酸比对，一般用+1/-3的简单计分，也可以用IUB码矩阵来处理简并碱基。不过核酸比对本身评分差异不大，重点还是矩阵和gap参数。
4.4 Gap penalties：空位的代价
插入和删除，在比对里表现为空位（gap）。BLAST通过gap open penalty和gap extend penalty来扣分。
gap open penalty：创建空位的初始惩罚。通常较大，防止随意插入。
gap extend penalty：空位每多一个字符扣的分。通常较小，因为一旦创建了空位，伸长它的代价相对低。
我踩过一个坑，在做多序列比对优化时，我设gap open=15，gap extend=5，结果出来的比对全是短空位，看着挺漂亮。但下游构建进化树时，分支关系全都乱了。仔细一看，那些短空位是强制压缩出来的，真实进化事件被拆碎了。后来我改成gap open=10，gap extend=2，空位分布才合理。
经验值参考：
- 核酸比对：gap open 5~7，gap extend 1~2
- 蛋白比对（严）：gap open 11~14，gap extend 2~4
- 蛋白比对（松）：gap open 8~10，gap extend 1~2
4.5 Filtering：别让重复元素欺骗了你
基因组里有大量重复序列、低复杂度区域（比如AAAAAAA、GCGCGCGCGC、或poly-Q）。这些区域在比对时会产生大量高分匹配，但毫无生物学意义。
Filtering就是做这件事的。它用算法（比如SEG、DUST、RepeatMasker）把这些区域标记成“XXX”，不让它们参与比对。你可以设置过滤条件，或者关掉过滤。
我曾经遇到过这样的情况：一个转录组项目中，我用BLASTn比对功能注释，结果有一半结果都指向人基因组的ALU元件。我很纳闷，后来发现是过滤没开。打开低复杂度过滤后，真正的功能匹配才浮出水面。
| 过滤模式 | 适用场景 | 注意事项 |
|---|---|---|
| 低复杂度过滤 | 基因组序列、重复区域多的序列 | 默认开启，不要轻易关 |
| 重复序列屏蔽 | 整个基因组分析 | 会移除已知重复元件 |
| 不过滤 | 短序列、探针设计 | 需要人工判断重复匹配 |
但过滤也不是万能的。如果你故意想找重复元件或转座子，那必须关掉过滤。同样的，短序列（比如60bp以下）可能过滤后就没剩多少有效信息了，这时候关掉或许更好。
4.6 把这些参数串起来
参数不是孤立的。它们是一个整体。低矩阵+宽松gap+短word size，意味着超高灵敏度，但结果会很冗长。高矩阵+严gap+长word size，结果精确但可能遗漏远源同源。
我建议你记住一个原则：先明确目标，再选参数调校。
# 一个典型的严格BLASTp参数设定
blastp -query myseq.fa -db nr \
-evalue 1e-10 \
-word_size 3 \
-matrix "BLOSUM62" \
-gapopen 11 \
-gapextend 1 \
-seg yes \
-out result.txt
看到这个配置了吗？这是做蛋白严格同源搜索的常用组合。如果你做的是宽松扫描，可以把e-value调高到1e-3，word_size降到2，gapopen降到9。