【生信】第二章 序列比对——Needleman-Wunsch全局比对
主要为基因组测序比对相关知识,部分内容作笔记自查使用。如有错误或遗漏还请海涵,可评论或邮箱联系。
最后修改时间:2020-04-03 07:53:55 星期五
Needleman-Wunsch全局比对
一、符号约定
符号 | 说明 |
---|---|
| 两个序列 |
| 在{x_i} 比对到{y_i} 情况下,序列x和y目前的最好比对得分(括号里的下标,代表当前已经比对完毕的位置) |
| 在{x_i} 比对到gap的情况下,序列x和y目前的最好比对得分(括号里的下标,代表当前已经比对完毕的位置) |
| 在{y_i} 比对到gap的情况下,序列x和y目前的最好比对得分(括号里的下标,代表当前已经比对完毕的位置) |
| A替换成B的分数 |
| 开启一个gap的罚分(正值) |
| 延伸一个gap的罚分(正值) |
二、简单模型
2.1 模型假设
- 符合全局比对要求,即x全串比对到y全串上
- 先不考虑开启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 模型假设
- 符合全局比对要求,即x全串比对到y全串上
- 引入开启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
的最右下角出发,按照正向路径相反的方向回溯,找到最优比对方法。