BLAST学习笔记 | Larry’s Space

type

status

date

slug

summary

BLAST——一种启发式算法

为了实现数据库搜索这个目标，我们最朴素的想法便是将序列和数据库中的所有序列通过Smith-Waterman算法进行Local Alignment，但是随着数据库数据量的爆炸式增长，这种朴素算法实际上变得不可行。以比较理想的情况来讲，两条短序列进行基于Smith-Waterman算法的Local Alignment需要0.2s，那么即使成百上千条短序列对比都会使得时间成本难以接受。

这时，美国国家生物信息中心（NCBI）的Eugene Myers、Stephen Altschul、Warren Gish、David J. Lipman及Webb Miller博士开发了一套启发式算法用来解决这个问题，即大名鼎鼎的BLAST。

BLAST的基本算法原理

想象人脑在对比两条相似序列的策略，大部分人在拿到两条序列后，都会先将最相似的部分对齐，之后再做调整。这种方法虽然得到的可能不是算法上的最优解，但也能说明很多问题了。BLAST的想法类似于上面这种场景。

BLAST的工作原理可以分为下面的一些步骤：

1.预处理库中已有数据

首先，我们需要引入一个数据结构中的概念：哈希表。哈希表是一种数据结构，用于实现高效的键值对（Key-Value）的存储和检索。它通过一个哈希函数，将键（Key）映射到表中的一个位置（索引），从而实现快速的数据访问。他有快速查找的特性，理想情况下，哈希表的查找、插入和删除操作的时间复杂度都是O(1)。

BLAST在预处理阶段会为数据库中的所有可能的word（并且会扩展到与其相似的算法）（通常来说，对于蛋白质是3-mer，对于核苷酸是11-mer）构建一个哈希表。每个word及其在数据库中的位置（序列ID和起始位置）会被存储在哈希表中。下面是一个示例：

哈希键（Key）	值（Value）
ATGCGTACGTA	DatabaseGene1: pos1, DatabaseGene2: pos11, DatabaseGene3: pos3
TGCGTACGTAG	DatabaseGene1: pos2, DatabaseGene2: pos12, DatabaseGene3: pos4
GCGTACGTAGC	DatabaseGene1: pos3, DatabaseGene2: pos13, DatabaseGene3: pos5

2.种子查找

与上一步类似BLAST会从需要查询序列中提取出特定所有可能的word（通常来说，对于蛋白质是3-mer，对于核苷酸是11-mer）。这些片段会在数据库中寻找相应的匹配。这一步骤的目的是快速锁定潜在的相似区域，而不需要对整个序列进行全局比对。

3.扩展

扩展(Extension)是BLAST工作流程中的第二个主要步骤。在种子查找（Word Seeding）阶段，BLAST已经锁定了一些潜在的匹配区域。扩展步骤的目标是将这些短暂的匹配区域扩展成更长的比对，以提高比对的准确性和覆盖度。

扩展的过程包括：

双向扩展：

左右两端同时扩展匹配区域，尽可能延伸到整个匹配片段。

允许一定的错配和缺口：

在扩展过程中，BLAST允许比对中存在一定数量的错配（mismatches）和缺口（gaps），但会根据设定的评分系统对这些变化进行惩罚。

评分与终止：

扩展过程会持续，直到比对得分不再增加或达到预设的阈值（如E-value阈值）。

生成局部对齐：

最终的扩展结果将形成一个局部对齐区域，显示查询序列与数据库序列之间的相似部分。

4.评分与过滤

评分机制

Raw Score（原始得分）

定义：

Raw Score 是BLAST在比对过程中直接计算的分数，基于匹配、错配和缺口的数量和性质。

类似于SP-Score（之前MSA算法提到的）。

与序列长度有关，不利于批判比对好坏。

Bit Score（位分数）

定义：

Bit Score 是一个标准化的比对得分，消除了数据库大小和评分矩阵的影响，使不同比对结果之间具有可比性。

计算公式：

λ（Lambda） 和 K 是与数据库大小和搜索空间相关的常数，由BLAST通过统计方法估计。

S 是Raw Score。

高Bit Score：比对结果更强，表示序列之间更高的相似性。

低Bit Score：比对结果较弱，表示序列之间相似性较低。

过滤机制

低复杂性区域过滤（Low Complexity Filtering）

去除序列中的重复或简单区域，减少由于这些区域引起的非特异性匹配，从而降低假阳性结果。

Dust（核酸）或 SEG（蛋白质）算法，这些算法识别并掩盖序列中的低复杂性区域：

Dust 主要用于核酸序列，识别重复或简单的碱基序列。

SEG 主要用于蛋白质序列，识别重复或简单的氨基酸序列。

软掩码（Soft Masking）与硬掩码（Hard Masking）：

软掩码：通过将低复杂性区域的碱基或氨基酸转换为小写字母来掩盖，不影响比对结果的显示。

硬掩码：通过完全移除或替换低复杂性区域的碱基或氨基酸来掩盖。

高复杂性区域强调（High Complexity Filtering）

类似于低复杂性区域过滤，强调序列中的高复杂性区域，进一步减少噪音。

5.统计评估（Statistical Evaluation）

E-value（期望值）介绍

定义：

E-value 表示在一个随机数据库中，预期会出现与当前比对得分相同或更高得分的比对次数。

它是评估比对结果显著性的关键指标。

计算公式：

K 和 λ 是与数据库大小和搜索空间相关的常数。

m 是查询序列的长度。

n 是数据库中序列的总长度。

S 是Raw Score。

E-value越低：比对结果越显著，意味着序列相似性更有可能具有生物学意义。

E-value越高：比对结果可能是随机发生的，相似性不一定具有生物学意义。

阈值设定：

常用阈值：1e-5（0.00001）

这意味着在随机数据库中，预期只有1e-5次这样的比对会出现。只有E-value小于设定阈值的比对结果才被认为是显著的。

统计学原理

BLAST假设比对得分（Score）遵循Gumbel极值分布，这是一种用于描述一组随机变量中的极值（最大值或最小值）的概率分布。在大量独立随机事件中，极端值（如最大得分）往往符合Gumbel分布。具体来说，BLAST使用Gumbel分布来建模在随机序列比对中出现的最高得分。

Gumbel分布的概率密度函数（PDF）和累积分布函数（CDF）分别为：

其中：

S：比对得分。

λ（Lambda）：尺度参数，反映得分分布的陡峭程度。

μ（Mu）：位置参数，通常与数据库大小和查询序列长度相关。

计算整个搜索空间中有 ≥1 个比对区域得分≥ S 的期望值得到：

令，得到：

NCBI BLAST 在线工具

NCBI（美国国家生物技术信息中心）提供了一套功能强大的在线BLAST（Basic Local Alignment Search Tool）工具，允许研究人员在无需本地安装BLAST软件的情况下，快速进行序列比对和分析。

NCBI BLAST的类型

此外，还有megablast用于高相似性核酸序列的比对，适用于非常相似的序列（如物种内部的基因比对）。discontiguous megablast适用于中等相似性的核酸序列比对，比megablast更具灵活性。

比对参数的设置

在“Algorithm parameters”部分，可以根据需求调整比对参数，常用的参数包括：

Max target sequences：每个查询序列返回的最大比对结果数。默认通常为100，但可以根据需要调整，如设置为10。

Expect threshold (E-value)：E-value阈值，用于过滤不显著的比对结果。常用阈值为1e-5或更低。

Word size：启动词的长度。较小的word size增加比对灵敏度，较大的word size提高比对速度。

Match/Mismatch Scores：核酸比对的得分系统。默认设置通常适用大多数情况。

Gap costs：缺口的开启和延伸罚分。默认设置通常适用，但在特定需求下可调整。

Low complexity filter：低复杂性区域过滤。开启可减少由于简单重复序列引起的非特异性匹配。

Database coverage：选择是否覆盖整个数据库，适用于大规模比对需求。

高级参数：

Subject filters：进一步过滤数据库中的序列，如掩盖低复杂性区域。

Organism：限制比对到特定物种的序列，提高比对效率。