BLAST：基于 seed-and-extend 的局部搜索

快速概览

BLAST (Basic Local Alignment Search Tool) 是生物信息学中使用最广泛的工具之一。它通过牺牲一定的敏感性，利用启发式算法在巨大的数据库中极速寻找局部相似片段。

理解 BLAST 的三步走策略：Word 列表生成、扫描命中与局部扩展
掌握 Seed-and-Extend 思想如何将全局搜索压力转化为局部计算
理解 E-value (期望值) 在衡量比对显著性中的核心作用
区分不同任务场景下的 BLAST 变体（blastn, blastp, etc.）

所属板块 核心方法

索引、比对、组装与概率模型构成的核心算法主轴。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系，再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉，再进入本页。

1. 是什么

BLAST（Basic Local Alignment Search Tool）是由 Altschul 等人于 1990 年开发的序列相似性搜索工具。它不追求找到数学上最优的局部比对，而是通过一套启发式策略（Heuristic），在极短的时间内找到数据库中与查询序列（Query Sequence）高度相似的片段。

BLAST 的核心思想可以概括为：先找到短的精确匹配（Seed），再从这些 Seed 向两侧延伸（Extend），最终形成高分局部比对（HSP）。这一策略将搜索空间从 $O(mn)$ 降低到接近 $O(m + n)$ ，使得在大规模数据库中的搜索成为可能。

2. 要解决什么生物信息学问题

核心问题

给定一条查询序列 $Q$ 和一个大型数据库 $D$ （可能包含数十亿条序列），找到 $D$ 中所有与 $Q$ 具有显著局部相似性的序列，并评估这种相似性的统计显著性。

典型应用场景

基因功能注释：将新发现的基因序列与已知数据库比对，推断其可能的功能。
物种鉴定：通过 16S rRNA 或 ITS 等标记基因的比对，鉴定微生物物种。
同源基因查找：在不同物种中寻找某一基因的同源物（Ortholog / Paralog）。
结构域预测：通过蛋白质序列比对识别保守的结构域。
引物特异性检查：验证 PCR 引物是否会与非目标序列结合。
基因组注释：对新测序的基因组进行自动功能注释。

3. 为什么不能直接用 Smith-Waterman？

虽然 Smith-Waterman 算法能保证找到最优局部比对，但它的复杂度是 $O(mn)$ 。

挑战：当你在搜索拥有数千亿碱基的 GenBank 数据库时，对每一条序列运行动态规划在计算上是不可接受的。假设查询序列长度为 300 bp，数据库总长为 $10^{11}$ bp，使用 Smith-Waterman 需要 $3 \times 10^{13}$ 次基本操作，在单机上可能需要数天。
BLAST 的哲学：大多数数据库序列与查询序列完全不相关。我们只需要快速找到那些”看起来像”有联系的片段（Seeds），然后只对这些区域进行精细比对。这种策略将搜索时间从”天”级降低到”秒”级。

敏感性与速度的权衡

BLAST 牺牲了一定的敏感性（Sensitivity）——它可能漏掉一些真实的弱同源关系——但极大地提高了速度。研究表明，对于大多数实际应用，BLAST 能找到 > 95% 的 Smith-Waterman 会找到的结果。

4. BLAST 的启发式三步法

第一步：构建 Word 列表（Seed Generation）

将查询序列（Query）拆分成长度为 $w$ 的小片段（Words）。

核酸（blastn）：

通常 $w = 11$ （默认值），且要求精确匹配。
查询序列被分解为所有长度为 11 的子串，共 $m - 10$ 个 word。

蛋白质（blastp）：

通常 $w = 3$ 。由于氨基酸的替换矩阵允许相似氨基酸之间的替换，BLAST 会使用得分矩阵（如 BLOSUM62）生成一组与原始 word 相似且得分超过阈值 $T$ 的”邻近词”（Neighboring Words）。
对于每个长度为 3 的 word，计算所有可能的单点突变后的得分，保留得分 $\geq T$ 的突变体。
阈值 $T$ 是敏感性与速度的关键调节参数： $T$ 越高，邻近词越少，搜索越快但可能漏掉弱同源； $T$ 越低，邻近词越多，搜索越慢但更全面。

第二步：数据库扫描（Scanning for Hits）

利用预先构建好的索引（如哈希表），在数据库中极速定位这些 word 的出现位置。

每个匹配的位置被称为一个 Hit（命中）。
精确地说，当数据库中的某个长度为 $w$ 的子串与查询的某个 word（或其邻近词）匹配时，就产生了一个 Hit。
命中位置定义了后续扩展的起点。

第三步：局部扩展（Extension）

从命中的 Seed 开始，向左右两个方向延伸比对。

策略：使用 Smith-Waterman 式的动态规划，计算延伸后的得分。只要得分在持续增加，就继续扩展。
X-drop 机制：如果当前得分从历史最高值下降了超过预设的门槛 $X$ ，则停止扩展。X-drop 的引入防止了在低质量区域浪费计算资源。
输出：最终形成 HSP (High-scoring Segment Pair)。只有得分超过阈值 $S$ 的 HSP 才会被报告。

扩展过程的形式化描述

设 $H$ 为当前最高得分， $X$ 为 X-drop 阈值。扩展过程中，如果当前累计得分 $s$ 满足：

$s < H - X$

则立即停止该方向的扩展。默认的 X-drop 参数通常为 15（蛋白质）或 20（核酸）。

5. Worked Example：blastn 的三步法

假设查询序列为 $Q = \text{ATCGATCGATCG}$ ，数据库中有一条目标序列 $D = \text{...GGATCGATCCGA...}$ 。

第一步：Word 列表（ $w = 4$ ，为简化演示）

$Q$ 被分解为：

ATCG, TCGA, CGAT, GATC, ATCG, TCGA, CGAT, GATC, ATCG

第二步：数据库扫描

在 $D$ 中搜索这些 word。假设 $D$ 的子串 ATCG 在位置 4 处被找到。

命中： $Q$ 的第 1 个 word “ATCG” 与 $D$ 的位置 4-7 匹配。

第三步：局部扩展

从命中位置向两侧延伸：

查询: ...ATCGATCGATCG...
            |  |  |  |
目标: ...GGATCGATCCGA...

向右扩展：

位置 8：Q 的 A vs D 的 A，匹配（+2），累计 +10
位置 9：Q 的 T vs D 的 T，匹配（+2），累计 +12
位置 10：Q 的 C vs D 的 C，匹配（+2），累计 +14
位置 11：Q 的 G vs D 的 C，错配（-3），累计 +11
得分下降 3，如果 $X = 5$ ，当前最高 14， $11 \geq 14 - 5 = 9$ ，继续
位置 12：Q 的 …，继续或停止

向左扩展：

位置 3：Q 的 - vs D 的 G，Gap 或 Q 的 A vs D 的 G，错配（-3），累计 -1
如果 $X = 5$ ，当前最高 14， $-1 < 14 - 5 = 9$ 成立，停止向左扩展

最终 HSP：

查询: ATCGATCG
        |||||
目标: ATCGATCC

得分：7 个匹配（+14） + 1 个错配（-3） = +11。如果此得分超过阈值 $S$ ，则报告该 HSP。

6. 结果的衡量：E-value 是什么？

BLAST 最重要的输出指标是 E-value (Expect value)。

定义

$E = K \cdot m \cdot n \cdot e^{-\lambda S}$

其中：

$m$ 是查询序列的有效长度。
$n$ 是数据库的有效长度。
$S$ 是 HSP 的原始得分（Raw Score）。
$K$ 和 $\lambda$ 是 Karlin-Altschul 统计量，取决于所用的打分矩阵和 Gap 罚分。
$e^{-\lambda S}$ 是得分为 $S$ 或更高的 HSP 在随机序列中出现的概率。

生物学解读

E-value 的含义：在一个规模相同的随机数据库中，由于偶然机会产生当前得分 HSP 的期望次数。
E-value 越小越好：E-value 越接近 0，说明比对结果越显著，越不可能是随机发生的噪音。
常用阈值：
- $E < 10^{-5}$ ：通常被认为具有显著的生物学同源性。
- $E < 10^{-10}$ ：非常显著的匹配，可信度很高。
- $E < 10^{-100}$ ：极其显著，几乎可以确定是同源关系。
- $E > 1$ ：该匹配在随机情况下也可能出现，没有生物学意义。

Bit Score

BLAST 同时报告 Bit Score，它是对原始得分的标准化：

$S' = \frac{\lambda S - \ln K}{\ln 2}$

Bit Score 的优势在于它不依赖于数据库大小，因此不同数据库搜索之间的 Bit Score 可以直接比较。而 E-value 会随着数据库的增大而增大（更大的数据库更容易产生偶然的高分匹配）。

P-value 与 E-value 的关系

P-value 是至少找到一个得分 $\geq S$ 的 HSP 的概率，而 E-value 是期望次数。当 E-value 很小时，两者近似相等：

$P \approx 1 - e^{-E} \approx E \quad (\text{当 } E \ll 1)$

7. 常见的 BLAST 变体

工具	查询序列（Query）	数据库（DB）	应用场景
blastn	DNA	DNA	寻找同源基因、物种鉴定。
blastp	蛋白质	蛋白质	功能预测、结构域分析。
blastx	DNA (6 帧翻译)	蛋白质	寻找新基因、分析 EST 序列。
tblastn	蛋白质	DNA (6 帧翻译)	在基因组中寻找特定的蛋白编码区。
tblastx	DNA (6 帧翻译)	DNA (6 帧翻译)	跨物种寻找保守编码区（最慢）。

blastx 的特殊价值

blastx 将核酸查询序列在 6 个阅读框（3 个正向 + 3 个反向互补）中全部翻译为蛋白质，然后与蛋白质数据库比对。这在以下场景中特别有用：

分析 EST（Expressed Sequence Tag）数据时，EST 可能包含测序错误或部分编码区。
在新测序的基因组中预测蛋白编码基因。
当核酸序列的保守性太低（同义突变积累），但蛋白质水平仍有信号时。

PSI-BLAST：迭代搜索

PSI-BLAST（Position-Specific Iterated BLAST）是 blastp 的增强版本。它在第一轮搜索后，用所有显著匹配构建一个位置特异性得分矩阵（PSSM / Profile），然后在后续轮次中使用这个 Profile 替代标准替换矩阵进行搜索。这使得 PSI-BLAST 能够发现首次搜索中遗漏的远缘同源关系。

PSI-BLAST 的流程：

使用标准 blastp 搜索。
收集所有 $E < 0.005$ 的匹配，构建多重序列比对。
从比对中提取 PSSM。
使用 PSSM 重新搜索数据库。
重复步骤 2-4，直到没有新的显著匹配出现（通常 3-5 轮）。

8. 复杂度分析

步骤	时间复杂度	说明
Word 列表生成	$O(m)$	扫描查询序列一次
数据库扫描	$O(m + n)$	使用哈希索引快速查找
局部扩展	$O(k \cdot L)$	$k$ 为命中数， $L$ 为平均扩展长度
总体	$O(m + n + k \cdot L)$	远优于 Smith-Waterman 的 $O(mn)$

在实际运行中，BLAST 的速度主要取决于命中数 $k$ 。对于低复杂度序列（如富含 AT/GC 的区域）， $k$ 可能非常大，导致速度下降。BLAST 通过低复杂度过滤（如 SEG 和 DUST 算法）来缓解这一问题。

适用前提

查询序列不能太短。对于 blastn，查询序列通常应至少 50-100 bp；对于 blastp，至少 15-20 个氨基酸。
数据库需要预先建立索引。NCBI 的在线 BLAST 服务已经完成了这一步；本地运行时需要使用 makeblastdb 构建索引。
对于高度重复的序列（如微卫星、转座子），BLAST 可能产生大量假阳性命中。

9. 与真实工具的连接

NCBI 在线 BLAST

NCBI 提供的在线 BLAST 服务（https://blast.ncbi.nlm.nih.gov）是最常用的序列搜索平台。它提供了上述所有变体的网页界面，并支持对结果进行过滤、下载和可视化。

本地 BLAST+

NCBI 的 BLAST+ 命令行工具套件是本地运行的推荐选择。主要命令包括：

makeblastdb：构建本地数据库索引。
blastn、blastp、blastx、tblastn、tblastx：各变体的搜索命令。
psiblast：PSI-BLAST 迭代搜索。

本地运行的优势在于可以搜索自定义数据库（如特定物种的基因组），不受网络延迟限制，且可以精确控制所有参数。

DIAMOND

DIAMOND（Buchfink et al., 2014）是 blastp 的超高速替代工具，速度可达 BLAST+ 的 500-20000 倍，同时保持相近的敏感性。它使用双索引策略和更高效的内存管理，特别适合大规模蛋白质数据库搜索（如宏基因组学中的蛋白质注释）。

BLAT

BLAT（BLAST-Like Alignment Tool）由 UCSC 开发，专门针对同一物种内的高相似度序列搜索（如将 cDNA 比对到基因组）。它使用索引策略而非 seed-and-extend，在 > 95% 一致性时比 BLAST 快得多。