2. BLAST算法原理：基本功，躲不开的
BLAST，全称 Basic Local Alignment Search Tool。这工具在病原菌基因组注释里，基本是天天见。我个人最早接触BLAST还是读研时跑一个16S rRNA鉴定，当时的师兄跟我说：“这玩意就是你序列的身份证查询系统”。后来干了这么多年，我慢慢理解，它远不止查个身份证那么简单。
2.1 序列比对基础：全局 vs 局部
分析BLAST之前，先搞清楚序列比对的两个大方向。
- 全局比对（Global Alignment）：把两条序列从头到尾全对齐。适合亲缘关系很近、长度差不多的序列。比如同一个毒力基因的完整序列比对。
- 局部比对（Local Alignment）：只找出两条序列中高度相似的局部片段。适合那些序列长度差很多、或只有部分保守区域的场景。
BLAST用的是局部比对。为什么？因为病原菌基因组里，大段大段的插入缺失太常见了。比如两个菌株的CRISPR区域，可能只有中间几个重复单元是类似的，两端完全不同。用全局比对效果很差，但局部比对能精准抓住那个保守核心。我在项目中遇到过多次这种场景——用全局比对硬跑，结果一堆gap，反而把真的同源区域给淹没了。
局部比对的数学基础，就是所谓的Smith-Waterman算法（动态规划）。它保证能找到最优局部比对，但问题是太慢。BLAST没直接用这个算法，而是做了个巧妙近似。
2.2 BLAST的基本搜索策略：三步走
BLAST比Smith-Waterman快在哪？说白了，它用了一个“先找线索，再展开验证”的策略。我们拆解来看：
- 种子词过滤（Seeding）
先把你的查询序列切成很多小片段（默认长度：蛋白3个氨基酸，DNA 11个碱基）。这些片段叫“种子词”（word）。然后去数据库里快速查这些种子词有没有精确匹配。这一步是速度关键，因为精确匹配可以用哈希表或者BWT索引，秒级完成。 - 种子词延伸（Extension）
找到种子词匹配后，在匹配位置向两边延伸，用不带gap的局部比对不断打分。当分值低于一个阈值时，停住。注意，这里最初是不允许gap的，这是BLAST早期的速度来源。后来的BLAST 2.0加入了gap扩展，但核心逻辑仍然是无gap先做一轮筛选。 - 高分对筛选与E值评估
延伸出来的序列片段称为HSP（High-scoring Segment Pair）。每对HSP都有一个原始分值。最后用统计学模型算出一个E值，用来评估这个比对结果是不是“运气好撞上的”。
# 以一个最简单的BLAST命令为例（人眼可读的步骤逻辑）
# 查询序列：ATCGATCGATCG
# 数据库序列：ATCGATCGATCG（完美匹配）
# 步骤1：种子词（默认11-mer）
# 从查询序列中提取所有11-mer
# ATCGATCGATC
# TCGATCGATCG
# 子词...（每个位置一个）
# 步骤2：在数据库中索引这些11-mer
# 命中第二个序列的匹配位置（第1位和第2位）
# 步骤3：向两侧扩展
# 原始得分 + 匹配得分 - 不匹配罚分
# 最终HSP分值 = 24（假设匹配+2，不匹配-3）
# 步骤4：计算E值
# E = K * m * n * e^(-λ * S)
# 实际公式更复杂，但基本就是：
# 数据库规模越大、期望随机出现的次数越多，E值越大。
2.3 得分矩阵：用什么尺子量相似度？
BLAST在延伸和打分阶段，需要一个替换矩阵（Substitution Matrix），也叫得分矩阵。它定义了不同残基替换的得分和罚分。对于DNA序列，通常用简单的匹配/不匹配模型（例如+1/-2）。但对蛋白序列，事情就复杂多了。
蛋白序列用BLOSUM或PAM系列矩阵。比如BLOSUM62，意思是基于62%序列一致性的保守区构建的矩阵。值越高，越偏向保守替换。我个人习惯用BLOSUM62，因为它在敏感性和特异性之间比较平衡。但如果你比对的是亲缘关系很远的物种（比如病原菌和真核宿主），可以试试BLOSUM45；如果同源性很高（比如亚种之间），BLOSUM80可能更好。
| 得分矩阵 | 适用场景 | 矩阵大小 | 我个人建议 |
|---|---|---|---|
| BLOSUM62 | 大多数通用蛋白比对 | 24×24 | 新手首选，默认就是它 |
| BLOSUM80 | 同源性高的序列（＞80%一致） | 24×24 | 比对高度保守毒力基因时用 |
| BLOSUM45 | 远源序列（＜30%一致） | 24×24 | 跨属病原菌基因组比对推荐 |
| PAM250 | 远源关系，但不如BLOSUM45常用 | 24×24 | 备选，一般不优先 |
矩阵选择影响E值和比对结果。选错了矩阵，就像用错误的尺子量东西——量出来的结果看起来是数字，但物理意义完全不对。
2.4 E值：那个让你心里有底的数字
BLAST输出结果里，你最该关注的其实是E值（Expect value）。
E值的定义很简单：在随机情况下，你期望看到至少一次同样分值的匹配的次数。注意，不是P值（P值是指本次匹配是随机产生的概率）。但实际解读时可以近似理解：E值越小，比对结果越不可能来自随机。
- E值 = 0：几乎不可能是随机产生的。你找到了真实的同源区域。
- E值 = 1e-50：50个零一个点后跟50个零，极小的概率。可以放心用。
- E值 = 0.01：仍有1%的概率是随机匹配。对于病原菌基因注释，通常不建议用这么高的E值。
- E值 = 5：每一百次比对应出现500次随机匹配。不靠谱。
E值的计算公式涉及数据库大小（m）和查询序列长度（n）。同样的原始比对分值，在更大的数据库里E值会变大。这就是为什么同一个比对，在nr库比在专用毒力因子库里的E值要大一些——数据库大了，随机匹配的机会就多。
# E值的经典公式（Karlin-Altschul statistic）：
# E = K * m * n * e^(-λ * S)
# 其中：
# S = 原始比对分值
# λ = 根据得分矩阵和序列组成计算出的常数
# K = 也是常数，通常≈0.1
# m = 数据库总长度（氨基酸或碱基数量）
# n = 查询序列长度
# 实际运行时BLAST会内部计算这些参数，不需要你手动算。
# 但了解这个公式，能帮你理解为什么E值会受数据库大小影响。
有些人在跑BLAST时只看百分比相似度（%identity），觉得这样更直观。但我觉得这是个大坑。因为高相似度（比如99%）也可能是基于很短的对齐长度，统计意义不大。E值或者位分值（bit score）才更能体现比对的可信度。