轻松理解BLAST搜索原理

BLAST: Basic Local Alignment Search Tool 基于局部比对搜索工具 （基于北京大学的生物信息学课程总结）

序列比对一般使用动态规划算法 (Simith-waterman 算法) 来实现， 对于两两比对(pairwise-alignment)而言，大O时间为Cmn (c为常量，表示单位操作时间， m， n分别为序列长度)。 如果使用动态规划算法比对数据库中的每一条序列，将会非常耗时。 故要采取其他的方法，BLAST是使用最为广泛的一种方法， BLAST也包含动态规划算法的步骤， 但其整个步骤是非常不同的。

大致来讲，BLAST采取seeding-and-extending方法。 首先找到query序列和subject序列间匹配的seed， 然后沿种子向两边延伸，产生High Scoring Segement Pairs (HSPs)， 然后在这些特定区域运用Simith-waterman 算法，最后评估比对的可信度。

Seeding过程: 沿query序列，以一定窗口(word length：w， 蛋白质序列默认为3个氨基酸，DNA序列默认为11bp核苷酸) 进行滑动，将query序列变成许多连续的”seed words”。

Seed搜索过程: 如果我们使用过stand-alone BLAST， 就会知道，在使用BLAST工具之前，需要为蛋白或核酸数据库建立索引 (makeblastdb 命令)， 待索引建立后， 所有在第一步中生成的seed可以根据索引快速定位到相似的seed和拥有该seed的序列。在BLAST的原始版本中,所有T得分大于18的hit会被用于后续的extending. T得分的计算基于核苷酸或氨基酸得分矩阵,如BLOSUM62矩阵。

Seed hit cluster: 然后我们就可以为hit到的sequence和query序列绘制对角线图， 好的比对一般是沿对角线平行的。 我们可以去掉其他比较凌乱的分布。

seed延伸过程: 沿每个找到的hit向序列两边延伸，此处使用Simith-waterman 算法，直到比对的得分低于某一阈值。并输出得分较高的比对, 一般默认输出50个比对。

可以采取一些方法来加快BLAST算法。 如屏蔽低复杂度(low-complexity)区域。如CACACACACACACACACA， KLKLKLKLKLKLKL。 低复杂度区在seed搜索时会产生很多假阳性。但比较麻烦的是我们不能通过定义一个列表来说明哪些是低复杂度区。只能基于query序列来定义哪些是低复杂度区。 这是通过K值来定义的， 小于某值则认为是低复杂度区。 L: 窗口长度 (默认蛋白质为3， DNA为11) N: 4 (DNA)或20 (蛋白质) ni表示每个碱基或氨基酸的个数。 如微卫星序列 CACACACACACACACA (假设窗口长度为6)， 其k值为0。36。

为了增加敏感性， 除了seed word序列自身，BLAST还使用他们高度相似的邻居种子(neighbourhood word)进行种子搜索。 如假设seed word是DKT， 基于氨基酸替代矩阵(如BLOSUM62)， DRT的得分为13 (6 +2+5)， 可以视为邻居种子。 一般得分大于11可以作为邻居种子。

对于一个给定的氨基酸，我们有1/20的机会有一个随机匹配。 所以对于一个长度为L的序列，随机产生一模一样的序列的概率是(1/20)L。 假设有一个长度为6个氨基酸的序列， 随机产生一模一样序列的概率是(1/20)6。 如果我们搜索swiss-prot数据库(有540958条序列共192206270氨基酸) 故可能产生一模一样序列的条数为: (1/20)6*192206270 = 3。

常用e-value来评估一个匹配是由于随机产生的可能性大小。 一个得分为S的序列比对，在一个数据库中随机产生一样得分的比对的期望数目即为e-value。 如果E=10，表示有10个有着这个得分的比对是由于随机产生的。 e-value的计算公式如下: m表示query 长度， n 表示数据库大小， S为匹配的得分， k 和λ 是用于normalization的，用来平衡不同矩阵和搜索空间对比对的影响。

当e-value小于0。1时， 和P值(概率值)线性关系非常好。 如下图， 当P等于0。05时， e-value等于0。0513，所以常用0。05作为e-value的cutoff。

BLAST是一种启发式算法，并不能保证找到最优的序列比对，但可以保证在一定时间内找到得分比较高的比对。 BLAST比纯粹的动态规划算法要快1000-10000倍，但灵敏度不如动态规划算法，特别是亲缘关系比较远的序列。