第二章 序列比对——Needleman-Wunsch全局比对

【生信】第二章 序列比对——Needleman-Wunsch全局比对

---

主要为基因组测序比对相关知识，部分内容作笔记自查使用。如有错误或遗漏还请海涵，可评论或邮箱联系。
最后修改时间：2020-04-03 07:53:55 星期五

---

Needleman-Wunsch全局比对

一、符号约定

符号	说明
xandy	两个序列
M\left( {i,j} \right)	在{x_i}比对到{y_i}情况下，序列x和y目前的最好比对得分（括号里的下标，代表当前已经比对完毕的位置）
X\left( {i,j} \right)	在{x_i}比对到gap的情况下，序列x和y目前的最好比对得分（括号里的下标，代表当前已经比对完毕的位置）
Y\left( {i,j} \right)	在{y_i}比对到gap的情况下，序列x和y目前的最好比对得分（括号里的下标，代表当前已经比对完毕的位置）
s\left( {A,B} \right)	A替换成B的分数
d	开启一个gap的罚分（正值）
e	延伸一个gap的罚分（正值）

---

二、简单模型

2.1 模型假设

1. 符合全局比对要求，即x全串比对到y全串上
2. 先不考虑开启gap和延伸gap的区别，将gap统一认定罚d分

---

2.2 模型构建

      设在{x_{1 \ldots i}}和{y_{1 \ldots j}}完全比对上时，最好的得分为F\left( {i,j} \right)，其中{x_{1 \ldots i}}和{y_{1 \ldots j}}不一定是x和y全串。此时可能出现3种情况：
       1. {{x_i}}比对上{{y_j}}，即前面的串{x_{1 \ldots i - 1}}和{y_{1 \ldots i - 1}}完全比对上
       2. {{x_i}}比对上gap，即前面的串{x_{1 \ldots i - 1}}和{y_{1 \ldots j}}完全比对上
       3. {{y_j}}比对上gap，即前面的串{x_{1 \ldots i}}和{y_{1 \ldots j - 1}}完全比对上
      综上，有

F\left( {0,0} \right) = 0

      及

F\left( {i,j} \right) = \max \begin{dcases}{F\left( {i - 1,j - 1} \right) + s\left( {{x_i},{y_j}} \right)} \\ {F\left( {i - 1,j} \right) - d} \\ {F\left( {i,j - 1} \right) -d} \end{dcases}

      有了上述公式，我们就可以利用动态规划思想，计算F\left( {i,j} \right)，填满下面每一个格子

 

 

 

---

2.3 模型使用

      以串x:AAG和y:AGC为例。gap罚分d=5。使用碱基替换矩阵如下：

 

 

      计算F\left( {i,j} \right)，填满动态规划矩阵，如下：

 

 

      注意，该矩阵是含有箭头的，意味着记录正向路径，方便反向回溯。
      整个矩阵的最右下角，即为两条序列全串的最优比对结果。此时，需要按正向路径（由箭头记录）相反的方向回溯，得到比对方法。

 

 

      在此例中，F\left( {2,1} \right)对应两个来源，因此上图中有两条路径可以回溯，得到2种比对方法。在此种打分规则下，2种方法都可以得到-6分，因此一样好，如下：

 

 

---

三、完整模型

3.1 模型假设

1. 符合全局比对要求，即x全串比对到y全串上
2. 引入开启gap的罚分d，和延伸gap的罚分e

---

3.2 模型构建

以有限状态机为思路，构建M、X、Y三种状态的转换公式。

 

 

• 状态M：
  状态M\left( {i,j} \right)的前一格共3种情况。
         1. 前一格{{x_{i - 1}}}和{{y_{j - 1}}}匹配上
         2. 前一格{{x_{i - 1}}}匹配gap
         3. 前一格{{y_{j - 1}}}匹配gap
  综上，有

  M\left( {i,j} \right) = \max\begin{dcases}{M\left( {i - 1,j - 1} \right) + s\left( {{x_i},{y_j}} \right)} \\ {X\left( {i - 1,j - 1} \right) + s\left( {{x_i},{y_j}} \right)} \\ {Y\left( {i - 1,j - 1} \right) + s\left( {{x_i},{y_j}} \right)} \end{dcases}

---

• 状态X：
  状态X\left( {i,j} \right)的前一格共2种情况。
         1. 前一格{{x_{i - 1}}}和{{y_j}}匹配上
         2. 前一格{{x_{i - 1}}}匹配gap
  综上，有

  M\left( {i,j} \right) = \max\begin{dcases}{M\left( {i - 1,j} \right) - d} \\ {X\left( {i - 1,j} \right) - e} \end{dcases}

  
  事实上，还有第3种情况，即前一格{{y_j}}匹配gap，但这会导致连续出现{{y_j}}对gap后{{x_i}}对gap，显然不如直接将两步合并成M\left( {i,j} \right)。

---

• 状态Y：
  状态Y\left( {i,j} \right)的前一格共2种情况。
         1. 前一格{{x_i}}和{{y_{j - 1}}}匹配上
         2. 前一格{{y_{j - 1}}}匹配gap
  综上，有

  M\left( {i,j} \right) = \max\begin{dcases}{M\left( {i,j - 1} \right) - d} \\ {X\left( {i,j - 1} \right) - e} \end{dcases}

  
        事实上，还有第3种情况，即前一格{{x_i}}匹配gap，但这会导致连续出现{{x_i}}对gap后{{y_j}}对gap，显然不如直接将两步合并成M\left( {i,j} \right)。

---

      有了上述公式，我们就可以利用动态规划思想，填充3张表，分别记录M\left( {{x_i},{y_j}} \right)、X\left( {{x_i},{y_j}} \right)和Y\left( {{x_i},{y_j}} \right)，如下：

 

 

      记录正向路径，正向路径可能会在2张表之间跳跃。之后从表M的最右下角出发，按照正向路径相反的方向回溯，找到最优比对方法。