在生物信息学研究中，序列比对是一个非常基础的问题，在很多研究中都会用到。本文对其中涉及到的全局比对（Needleman-Wunsch算法）、局部比对算法（Smith-Waterman 算法）、Blast、PSI-Blast、PSSM、HMM的原理进行总结和介绍。

Needleman-Wunsch算法——基于动态规划的序列全局比对

Needleman-Wunsch算法是序列全局比对的核心算法，具体实现如下。

打分表

对于不完全匹配的序列，会存在空位或者不匹配的位点，打分表是表示一种氨基酸（或核苷酸）变为另外一种氨基酸（核苷酸）的打分，对于空位则要罚分，最终综合所有位点的打分结果，获得两个序列的匹配分值，分值越高的表示两个序列相似度越高。

对于蛋白质序列，常用的打分表有PAM250和BLOSUM62，其中PAM250是基于高同源序列构建，BLOSUM62基于远程同源序列构建，因此寻找远程同源序列一般用BLOSUM62。

计算步骤

为便于计算和表示，以核苷酸序列为例。

1. 定义打分表

空位罚分为 -5。

2. 待比对序列的表示

对于长度为3的两个序列“AAG”和“AGC”，设计一个5*5的表用于序列比对。

3. 打分规则

对于表格中的空格， 依据下列打分规则进行填充：

其中，$s(x_i,y_j)$表示打分表中对应氨基酸替代情况的打分。

(1): $f(i-1,j-1)+s(x_i,y_j)$ 表示 $x_i$ 比对到 $y_j$

(2): $f(i-1,j)+d$ 表示 $x_i$ 比对到缺失“-”

(3): $f(i,j-1)+d$ 表示 $y_j$ 比对到缺失“-”

4. 打分结果

按照上一步的打分结果，整个比对的打分结果为：

5. 回溯获得比对结果

根据上述打分结果进行回溯，即可得到比对结果：

这里得到两个得分为“-6”的全局比对结果，分别为：

1
2
AAG-
-AGC

1
2
AAG-
A-GC

这两个即最好的全局比对，而得到为“-6”。

Smith-Waterman 算法——基于动态规划的序列局部比对

Needleman-Wunsch算法是对两个目标序列进行全局比对，但是局部比对也非常重要，因此后来在其基础上，Smith和Waterman提出了Smith-Waterman 算法。

1. 打分规则的修改

Smith-Waterman 算法对Needleman-Wunsch算法的打分规则做了改进，使之适用于局部序列比对，具体为：

2. 打分结果

相当于将其中负值的打分情况终止了，重新开始比对，获得了局部比对。

3. 回溯结果

由此获得了两个局部比对：

1
2
AG
AG
1
2
A
A

BLAST 算法——启发式近似算法

BLAST是速度比Smith-Waterman快50倍，不如Smith-Waterman算法精确的序列局部比对算法。具体流程如下：

1. 序列预处理

移除待比对序列的低复杂度区域或重复序列。

2. 构建待查询的K-元 词表

对于待比对序列，从序列第一个元素开始，对每个位置开始的连续K个字母作为一个词（word），将其key和位置存入hash表。

对于数据库中的序列做同样的“分词”工作。

对于DNA序列，通常选择长度为11的核苷酸链，对于蛋白质序列，通常选择3个氨基酸作为一个词。

3. 比对打分

根据打分矩阵为所有待比对序列的“word”和数据库中的“word”的比对打分，找出所有打分高于T的比对，作为“种子”（seed）。

4. Seed优化

减少seed的数目可以减少下一步extend的消耗，可以将距离相近的seed合并为一个seed。

5. 基于“种子”扩展（extend）

对于每一对选择出来的种子，将其向两边延伸，使其在尽可能长的距离得到尽可能多的分数。具体方法是逐渐向两边扩展，规定常数D，在扩展到分数为HighestScore-D时停止，如果最终得分大于得分阈值S，则将其设为高分区域HSP（high-scoring segment pair）。

6. 合并相邻的或距离较近的HSP
7. Smith-Waterman局部比对结果

对于每部分HSP使用Smith-Watermans算法进行局部比对，为每部分打分，作为最终结果。

PSI-BLAST 算法

1. 初始序列比对

用初始的一条蛋白质序列作为query，用BLAST基于BLOSUM62在蛋白质数据库中搜索，获得初始结果

2. 构建PSSM

基于初始的多序列比对（MSA，multiple sequence alignments）结果构建初始profile，即PSSM（position-specific score matrice）

3. 基于PSSM的再搜索

用PSSM作为新的打分矩阵，搜索数据库，发现较远的相似性序列，并给出E-value

4. 基于PSSM的迭代搜索

根据新的MSA，构建新的PSSM，再迭代搜索过程，直到没有符合要求(E-value threshold,比如0.005,有时可放宽至0.01)的新序列，或者达到指定的迭代次数，一般为3。

PSSM——position-specific score matrice

PSSM是从MSA中构建的，其典型计算过程为：

1. 计算MSA中每个位点不同氨基酸的出现次数

2. 基于氨基酸次数计算频率

第1列：$f{A,1} = \frac{0}{5} = 0, f{G,1} = \frac{5}{5} = 1, …$

第2列：$f{A,2} = \frac{0}{5} = 0, f{H,2} = \frac{5}{5} = 1, …$

$… … $

第15列：$f{A,15} = \frac{2}{5} = 0.4, f{C,15} = \frac{1}{5} = 0.2, …$

为了解决有些氨基酸类型出现为0的情况，给分子、分母增加pseudo-counts：

第1列：$f{A,1}^{‘} = \frac{0+1}{5+20} = 0.04, f{G,1}^{‘} = \frac{5+1}{5+20} = 0.24, …$

第2列：$f{A,2}^{‘} = \frac{0+1}{5+20} = 0.04, f{H,2}^{‘} = \frac{5+1}{5+20} = 0.24, …$

$… … $

第15列：$f{A,15}^{‘} = \frac{2+1}{5+20} = 0.12, f{C,15}^{‘} = \frac{1+1}{5+20} = 0.08, …$

3. 计算对数似然

$q_{i}$为相应氨基酸类型在随机序列中出现的期望频率，比如可以取平均值，对于氨基酸为$\frac{1}{20}=0.05$，对于核苷酸为：$\frac{1}{4}=0.25$。

说明

以上为典型的PSSM计算步骤，但是在具体实践上，不同的软件可能会对细节进行修改，从而得到的PSSM并不一样，如PSI-BLAST会采用BLOSUM62处理pseudo-counts问题。

另外，PSSM不处理MSA中位点插入和删除的情况。

基于HMM的序列比对

由于PSSM不能处理位点插入和删除的情况，因此后来开发了基于HMM（Hidden Markov Models）的序列比对，采用HMM对所有情况进行建模，这里对其基本原理进行介绍。

下面是对于一个序列比对的示例：

如果采用PSSM或其他类似的Profile方式，都无法对位点的插入和删除进行建模，而HMM则可以。这里不对HMM的原理进行介绍，进对其在序列比对中的应用方式进行简要介绍。

1. 初始的HMM

2. 增加对插入的处理

3. 增加对删除的处理

4. 完整的用于序列比对的HMM

5. HMM的表现形式

最终，对于一个氨基酸位点的HMM为包含20种氨基酸和7种转移状态（M->M、M->I、M->D、I->M、I->I、D->M、D->D ）概率的数组。

其中M、I、D的含义分别为，M：match，I：insertion，D:deletion。

参考：

  1. 生物信息学：导论与方法(北京大学)
  2. BLAST 算法
  3. 生物序列局部比对之Blast算法
  4. PSI-BLAST寻找同源蛋白
  5. An introduction to Patterns, Profiles, HMMs and PSI-BLAST
  6. Profile HMMs for Sequence Alignment

最后更新： 2019年04月19日 21:37

原始链接： http://andersjing.com/2018/11/25/2018-11-25-sequence-alignment/

赏
× 请打赏~
打赏二维码