第二章 BLAST算法详解：Seeding阶段、Extending阶段、E-value计算与统计显著性
如果你用过BLAST，很可能只是输个序列、点个按钮、等结果。但后台那几步到底怎么跑的？我当年第一次看源码时，说实话有点懵。后来折腾过几轮大规模比对优化，才算摸清门道。
BLAST的核心可以拆成三块：Seeding（种子生成）→Extending（延伸）→E-value评估。咱们一章一章啃明白。
Seeding阶段——从海量序列里快速锚点
说白了，Seeding就是穷举法加一点巧思。把查询序列切成固定长度的小片段，用这些片段去数据库里碰瓷。碰上的就叫“种子”，后面才值得延伸。
种子长度由 -word_size 控制。默认值：
| 程序 | 默认 word_size | 适用场景 |
|---|---|---|
| blastn | 11 | 核酸精确比对 |
| blastp | 3 | 蛋白质中度相似 |
| blastx | 3 | 翻译核酸 vs 蛋白 |
word_size越小，灵敏度越高，但计算量暴涨。我做过一个全长转录组注释的项目，刚开始用默认11，很多跨物种同源序列根本对不上。后来降到7，灵敏度上来了，但运行时间翻了三倍。嗯，这就是典型的trade-off。
避坑：我曾经为了追一个低相似度同源区域，把 word_size 降到6，结果一天还没跑完。后来一看，数据库是重复序列富集的——种子太多，假阳性全堆在延伸阶段。建议先扫一遍重复序列，或者用 -lcase_masking 软屏蔽。
Extending阶段——双向延伸与打分矩阵
找到种子后，BLAST会从seed位置向左右两边尝试延伸——不加空位，每次加一对碱基/氨基酸，根据打分矩阵算分。分数低于门限就停。这一步叫“无空位延伸”。
打分矩阵才是灵魂。蛋白质比对通常用BLOSUM62，核酸用+1/-2简单计分。我个人习惯先用BLOSUM62做一个粗筛，再用PAM250做精细边界调优。为什么呢？BLOSUM62对中度保守区域识别好，PAM250更擅长远处同源。
技巧：调整 -X（dropoff值）可以控制延伸的激进程度。默认值15可能漏掉长插入片段。我在比较真菌漆酶家族时把 -X 改到30，召回率明显提升。
无空位延伸结束后，BLAST会尝试在两端加空位——这一步叫“有空位延伸”。空位开放罚分（-gapopen）和空位延伸罚分（-gapextend）直接影响比对质量。我建议的起步参数：
blastp -query query.fa -db nr -outfmt 7 -gapopen 11 -gapextend 1 -evalue 1e-5gapopen 不要小于8，否则会出现大量莫名其妙的长空位。我踩过这个坑。
E-value计算与统计显著性——你match到的到底靠不靠谱
E-value应该是被误读最多的BLAST输出。很多人以为“E-value小于0.05就显著”，其实不完全对。
E-value的数学含义是：在随机背景下，期望出现等于或高于当前比对分数的最少条数。公式核心：
E-value = K * m * n * e^(-lambda * S)其中K、lambda是数据库相关常数，m是查询长度，n是数据库总长度，S是原始比对分。你看，E-value跟数据库大小直接相关——数据库越大，相同分数对应的E-value可能越大。
举个例子：你在一个小型病原菌库比对到一段序列，E-value=1e-10；换到NCBI nr库（大10倍），同样分数可能只剩1e-7。所以跨库比较时，别只盯着E-value绝对值。
面试常考：E-value = 0.05 意味着什么？意味着你这条查询序列，在随机序列库里跑100次，有5次可能撞到同样好的分数。并不是“有95%的概率是同源”。
我自己的经验阈值：
- E-value ≤ 1e-50：高度保守，直接正交验证
- 1e-10 ~ 1e-50：可能同源，需要看覆盖率
- 1e-3 ~ 1e-10：弱信号，结合结构域分析
- > 1e-3：基本不可信，除非有额外证据
曾经有一次，我比对某寄生物基因组到宿主数据库，E-value全在1e-5到1e-3之间。差点就以为是污染序列，后来发现是水平基因转移的真实案例——幸亏我同时检查了alignment长度（>100aa）和SCP（得分/长度比），没武断过滤。
最后提一下 -culling_limit 和 -hsp_max 的配合。默认BLAST对一条查询最多输出500个HSP（高分片段）。E-value低的片段可能被大量高分数击穿。我建议在注释管线里加上：
blastp -query query.fa -db ref_prot -out xml -culling_limit 1 -hsp_max 50这样每个查询只返回最优的一个HSP，避免冗余延伸把你的终端刷屏。