【现学现卖】序列比对之算法

前期相关推送（序列比对结果中的一些值的意义）：

【现学现卖】序列比对之identity VS similarity

【现学现卖】序列比对之bit-score VS E-value

在bit-score和E-value的公式中，S为原始矩阵得分，其他常数与选择的算法相关。那么什么是序列比对的矩阵得分？序列比对都有哪些算法呢？

一、序列比对

1. 序列比对的理论基础

基础是进化学说，如果两个序列相似性高，则推测同源/有共同的进化祖先，它们是经过序列残基的一些列变化（替换、缺失、重组等）演化而来的。

比如下面两个序列：

这两条序列要实现更好的比对，就会引入空位。序列残基的替换导致单个残基不匹配（红色短线），残基的插入或缺失导致空位（红色长线）。

2. 序列比对的两种数学模型

主要模型是全局比对（Global alignment）和局部比对（Local alignment）。它们分别从整体和局部反映序列的特征，现实使用中，局部比对使用较多，生物序列往往不会全长相似，而是局部相似，所以局部比对有更高的灵敏度，结果更具生物学意义（当然如果本身就是选定的一段很短的目标序列，那么两种模型的结果应该差异不大；换句话说，全局比对适合比较长度相似的序列，而局部比对可以比对长度相差较大的序列）。

3.序列比对打分矩阵

（1）核苷酸序列

等价矩阵：相同核苷酸则赋值为1，不同为0。

BLAST矩阵：经验总结。

转移矩阵/转换-颠换矩阵：嘌呤-嘧啶的转换/颠换+经验总结。

（2）氨基酸序列（表太多，太大了，我就不在这里粘贴了）

氨基酸序列的这两个常用打分矩阵都属于替换矩阵，主要是考虑在进化过程中，不同氨基酸的替代对蛋白质功能和结构的影响不同，所以用简单的比对相同或者不同不足以描述两个氨基酸残基的关系。

PAM矩阵：Point accepted mutation matrices。基于进化的点突变模型，分析同源蛋白在进化中氨基酸变化的可能性。有PAM1-250，后面的数字越小表示亲缘关系越近，PAM1是similarity>85%的序列计算产生的。

BLOSUM矩阵：Blocks amino acid substitution matrices。基于蛋白质模块数据库，以序列片段为基础。有BLOSUM45，62，80等，数字越大，亲缘关系越近，BLOSUM80是identity>80%的序列计算得到的。

根据比对的两个序列的亲缘关系远近选择合适的矩阵，如果不清楚，一般选择PAM120或BLOSUM62。

二、序列两两比对

经典算法（动态规划算法）是Needleman-Wunsch算法（整体比对算法）和Smith-Waterman算法（局部比对算法）。两种算法均可以用于核苷酸和氨基酸序列，给定空位罚值和打分矩阵后，给出最高比对值的排列。

三、多序列比对

多序列比对的算法是基于渐进比对，在序列两两比对算法的基础上逐步优化的结果。

目前发展出来的程序有CLUSTALW（累进算法），MUSCLE（迭代算法），MAFFT等（上图选项，其中CLUSTALW最为人熟知，传说MUSCLE和MAFFT的精度和速度都优于CLUSTALW，由于我的数据还算简单，也不多，感觉没什么大差别。考虑到它们迭代可以纠正初始比对错误，所以理论上精度会提高）。

四、序列数据库搜索

数据库搜索可以说是序列相似性比对最有价值的应用，我最熟悉的是BLAST和antiSMASH。

BLAST算法

将查询序列分为短片段（局部比对），筛选库中具备这些片段的序列，然后将匹配的序列片段延伸（插入，gap等），根据矩阵计分排序，显示结果。

antiSMASH工具算法

它的算法是CASSIS（Cluster Assignment by Islands of Sites）。antiSMASH是对微生物次生代谢产物基因簇的序列搜索工具，由于生成这些化合物的基因都是成簇的，所以它的算法中重要的是找到基因簇。