生物信息学序列比对算法——动态规划

前言

序列比对是生物信息学中非常重要的一个概念，对分析生物数据具有不可或缺的作用。目前绝大多数的序列比对工具均包含了基于动态规划的序列比对的步骤。序列比对问题类似于最长公共子序列问题 （Longest Common Subsequence problem, LCS) ，进而衍生出了全局比对算法尼德曼-翁施算法 （Needleman-Wunsch Algorithm） 与局部比对算法史密斯-沃特曼算法 （Smith-Waterman algorithm） 。

一、LCS问题

1. 子序列

如果给出一条序列 X = { x 1 , x 2 , . . . , x m } X = \{x_1,x_2,...,x_m\} X={ x1​,x2​,...,xm​}，并且有另外一条序列 Z = { z 1 , z 2 , . . . , z k } Z=\{z_1,z_2,...,z_k\} Z={ z1​,z2​,...,zk​}, 存在严格递增的下标序列 { i 1 , i 2 , . . . , i k } \{i_1,i_2,...,i_k\} { i1​,i2​,...,ik​}，并且对于$j=1,2,...,k$有$z_j=x_{i_j}$，那么称序列$Z$为序列$X$的子序列。如序列 Z = { B , C , D , B } Z = \{B,C,D,B\} Z={ B,C,D,B}为序列 X = { A , B , C , B , D , A , B } X=\{A,B,C,B,D,A,B\} X={ A,B,C,B,D,A,B}的子序列，并且存在对应的递增下标序列 { 2 , 3 , 5 , 7 } \{2,3,5,7\} { 2,3,5,7}。

2. 公共子序列

给出两条序列$X$和$Y$，如果序列$Z$是序列$X$和$Y$的共同的子序列，则称$Z$为序列$X$和$Y$的公共子序列。
设 Z = { z 1 , z 2 , . . . , z k } Z=\{z_1,z_2,...,z_k\} Z={ z1​,z2​,...,zk​}是序列 X = { x 1 , x 2 , . . . , x m } X=\{x_1,x_2,...,x_m\} X={ x1​,x2​,...,xm​}与序列 Y = { y 1 , y 2 , . . . , y n } Y=\{y_1,y_2,...,y_n\} Y={ y1​,y2​,...,yn​}的最长公共子序列，则有：
（1）如果$x_m=y_n$，则有$z_k=x_m=y_n$，并且$Z_{k-1}$是序列$X_{m-1}$与$Y_{n-1}$的最长公共子序列；
（2）如果$x_m\ne y_n$，且$z_k\ne x_m$，则序列$Z$是序列$X_{m-1}$和序列$Y$的最长公共子序列；
（3）如果$x_m\ne y_n$，且$z_k\ne y_n$，则序列$Z$是序列$X$和序列$Y_{n-1}$的最长公共子序列；

由此可见，两个序列的最长公共子序列包含两个序列前缀的最长公共子序列。因此，最长公共子序列问题具有最优子结构性质。
根据最长公共子序列问题的最优子结构性质，可以建立子问题最优值的递归关系。使用$c[i][j]$表示序列 X i = { x 1 , x 2 , . . . , x i } X_i=\{x_1,x_2,...,x_i\} Xi​={ x1​,x2​,...,xi​}与序列 Y j = { y 1 , y 2 , . . . , y j } Y_j=\{y_1,y_2,...,y_j\} Yj​={ y1​,y2​,...,yj​}的最长公共子序列的长度。当$i=0$或者$j=0$时，显然最长公共子序列是空的，此时$c[i][j]=0$。进而可以根据最优子结构特性构建递归关系，如下所示：
c [ i ] [ j ] = { 0 , i = 0 ∨ j = 0 c [ i − 1 ] [ j − 1 ] + 1 , i , j > 0 ∧ x i = y j max ⁡ { c [ i ] [ j − 1 ] , c [ i − 1 ] [ j ] } , i , j > 0 ∧ x i ≠ y j \begin{equation} c[i][j]= \begin{cases} 0, & i=0 \vee j=0 \\ c[i-1][j-1]+1, & i,j > 0 \land x_i=y_j \\ \max\{c[i][j-1], c[i-1][j]\}, & i,j>0 \land x_i\neq y_j \end{cases} \end{equation} c[i][j]=⎩ ⎨ ⎧​0,c[i−1][j−1]+1,max{ c[i][j−1],c[i−1][j]},​i=0∨j=0i,j>0∧xi​=yj​i,j>0∧xi​=yj​​​​
如果需要构造最长公共子序列的话，我们需要一个辅助矩阵$b$来记录构造信息。每一条记录$b[i][j]$指示到达$c[i][j]$值的子问题。
在构造最长公共子序列时，我们首先从$b[m][n]$开始，沿着箭头指示的方向搜索矩阵$b$。
（1）当$b[i][j]=“\nwarrow ”$时（这意味着元素$x_i=y_j$ LCS的一部分），表示序列$X_i$与序列$Y_j$中的公共子序列为序列$X_{i-1}$与序列$Y_{j-1}$的公共子序列在尾部增加元素$x_i$或者$y_j$；
（2）当$b[i][j]=“\uparrow ”$时，这表明此时序列$X_i$与$Y_j$的最长公共子序列和序列$X_{i-1}$与序列$Y_j$的最长公共子序列相同；
（3）当$b[i][j]=“\leftarrow ”$时，这表明此时序列$X_i$与$Y_j$的最长公共子序列和序列$X_i$与序列$Y_{j-1}$的最长公共子序列相同。
构造最长公共子序列的矩阵$c$的时间复杂度为$O(mn)$，构造最长公共子序列的时间复杂度为$O(m+n)$。

二、Needleman Wunsch

Needeman-Wunsch也采用动态规划算法，与LCS算法类似，但构造矩阵$c$不同。LCS中长度的概念不再用于构造矩阵$c$，确认代之的是使用最优比对结果（使用得分表示）用于构造递归关系。与此同时，会引入罚分的概念。
为了与LCS算法表示区别，使用函数$P(s_t,t_j)$表示打分函数，使用矩阵$d$表示得分情况。在这里，我们设置如下的打分函数（单位阵）：
$$\begin{array}{c}P(s_t,t_j)=\{\begin{array}{ll}1,& s_i=t_j\\ 0,& s_i\ne t_j\\ -1,& s_i=gap\vee t_j=gap\end{array}\end{array}$$
$$\begin{array}{c}d[i][j]=\max \{\begin{array}{}\end{array}\end{array}$$