2. BLAST原理简介
嗯，我们正式开始聊BLAST的原理。说实话，很多教程一上来就甩公式，看得人头皮发麻。我个人觉得，搞懂BLAST，关键在于想明白一个问题：它凭什么能在海量序列里快速找到跟我手里这玩意儿相似的片段？
说白了，BLAST的全称是Basic Local Alignment Search Tool——基础局部比对搜索工具。注意“局部”这两个字，它不追求把两条序列从头比到尾，而是找那些局部相似的区域。我在项目中遇到过很多次，明明两条序列整体差异很大，但局部功能域保守得像亲兄弟，这时候BLAST就比那些全局比对工具好用多了。
核心思想：BLAST利用“短种子-双方向延伸”的策略，牺牲一点点灵敏度，换来指数级的速度提升。你想想看，如果每次都要做一遍Needleman-Wunsch，那比对一个基因可能要等一整天。
算法流程：三阶段走
BLAST的工作流程，大致可以拆成三个阶段。我习惯叫它们“切-找-扩”。
- 第一阶段：查询分割（Seed）
先把你的查询序列（query）切成一个个固定长度的小片段，这个长度叫
word size。蛋白比对的默认是3个氨基酸，DNA是11个碱基。举个例子，如果查询序列是
ACGTACGT，word size=4，那就会切出ACGT、CGTA、GTAC、TACG、ACGT。 - 第二阶段：种子查找（Search）
把这些种子片段，去目标数据库里快速匹配。用的是哈希表或者类似布隆过滤器的方法，看数据库里有没有完全一样或者高度相似的片段。
这一步快得离谱，因为它是精确匹配，不用算相似度。匹配上的位置叫“命中”（hit）。
- 第三阶段：延伸与扩展（Extend）
找到命中之后，BLAST不会停下来。它会向两端继续延伸比对，用动态规划（类似Smith-Waterman算法）算出一个高分比对片段（HSP，High-scoring Segment Pair）。
我踩过的坑：曾经我想提高对长片段重复序列的检出率，盲目把word size改得很小（比如DNA设成7）。结果呢？全是对齐垃圾，命中太多了，延伸阶段慢得像蜗牛。后来才发现，word size是个双刃剑——设小了，灵敏度高但噪音大、速度慢；设大了，速度快但可能漏掉真阳性。你得根据自己的数据特征去调。
评分系统：怎么才算“相似”？
找到候选区域之后，BLAST需要给每对区域打个分。这个分数决定了谁排在前面，谁被过滤掉。评分系统主要包含两个要素：替换矩阵和空位罚分。
替换矩阵（Substitution Matrix）
蛋白比对用的最多的是BLOSUM系列（比如BLOSUM62）和PAM系列。这些矩阵存的是：某个氨基酸突变成另一个氨基酸，在进化上“合理”的程度。
为什么不用简单的“匹配得1分，错配扣0分”？因为那太粗糙了。实际上，保守替换（比如亮氨酸Leu变成异亮氨酸Ile）应该拿高分，而激进替换（比如精氨酸Arg变成天冬氨酸Asp）应该被严重扣分。
| 矩阵 | 适用场景 | 我的经验 |
|---|---|---|
| BLOSUM62 | 通用蛋白比对，中等进化距离 | 默认不解释，95%的蛋白BLAST都用它 |
| BLOSUM80 | 高度相似的序列（比如同源物种） | 比对近缘物种时，用这个跑出来的结果更干净 |
| BLOSUM45 | 远缘序列，进化距离很远 | 容易漏掉真阳性，我一般用它做粗筛 |
| PAM250 | 远距离进化分析（古老蛋白家族） | 现在用得少了，被BLOSUM系列取代 |
空位罚分（Gap Penalties）
插入或删除一个碱基/氨基酸，BLAST需要惩罚。因为它破坏了序列的连续性。空位罚分由两个参数控制：
- 空位开放罚分（Gap Open Penalty, GOP）：打开一个空位就扣一次，通常设成11（蛋白）/5（DNA）。
- 空位延伸罚分（Gap Extension Penalty, GEP）：空位每多一个碱基，再扣一次，通常设成1（蛋白）/2（DNA）。
注意：空位开放罚分是一个“一次性”的成本。所以你想想看，如果GOP设得太高（比如15以上），BLAST就不愿意插入空位了，哪怕两边序列其实有很长的插入片段，它也会硬对齐，导致评分降低。反之，如果GOP设得太低（比如5以下），就会疯狂插入空位，比对结果看起来全是错位，噪声爆炸。我见过不少新手把空位罚分乱调，结果跑出来一堆假阳性。
统计显著性：E-value才是裁判
有了得分还不够。你得知道这个得分是“真牛”还是“运气好”。BLAST引入了一个叫期望值（E-value）的概念。
E-value表示：在随机情况下，拿到当前得分（或更高得分）的比对结果的期望个数。很简单的一个原则：E-value越小，比对结果越可靠。
- E-value < 1e-5（蛋白比对）：通常认为是显著的，可以接受。
- E-value < 1e-10：非常强的同源性证据。
- E-value > 0.01：基本是噪声，别拿它当证据。
举个例子，你跑了个BLAST，得到一条比对结果，得分为200分，E-value为1e-25。这意味着，如果你把查询序列和一库随机序列做比对，大概每10的25次方次随机实验里，才可能看到一次200分。这显然不太可能是凑巧的——所以这条比对结果很可能是真的同源片段。
实际案例：用核酸BLAST（blastn）做一次简单比对
咱们动手试一下，用命令行跑一个最简单的核酸BLAST。假设有一个查询文件query.fa和一个本地数据库my_db.fasta。
第一步，构建数据库索引：
makeblastdb -in my_db.fasta -dbtype nucl -parse_seqids -out my_db第二步，运行blastn：
blastn -query query.fa -db my_db -out results.txt -outfmt 6 -evalue 1e-5解释一下参数：
-outfmt 6：输出成表格格式，简练易读。-evalue 1e-5：只保留E-value小于1e-5的结果。
你会看到输出格式像这样：
query_1 ref_1 98.5 450 5 0 1 450 1 450 0.0 850
query_1 ref_2 75.0 300 30 2 50 349 100 399 1e-45 320每一列分别表示：查询ID、目标ID、相似度百分比、比对长度、错配数、空位数、查询起止、目标起止、E-value、得分。
嗯，看到这里，你对BLAST的算法骨架应该已经有感觉了。说白了，它就是用“切-找-扩”的思路，通过精心设计的评分系统来快速锁定相似区域。剩下的就是参数调优和结果解读的问题了。