aibiology

Artificial intelligence in biology

0%

序列比对中的动态规划算法

Posted on Edited on

动态规划

动态规划(dynamic programming), 是一种求解多阶段决策过程最优化问题的算法。 核心思想就是: (1),将一个大问题拆解成一堆子问题，(2),求解出子问题的解，存储起来当再次需要求解此问题时候可以被重复利用，避免重复的计算过程。 能够使用DP必须满足下面三个特征:
(1), 最优子结构性质。一个问题的最优解包含其子问题的最优解。
(2), 重叠子问题性质。在求解子问题的过程中，会有大量的子问题的是重复的，可以将子问题的结果存储，后续直接查询，避免重复计算。 (3), 无后效性。子问题的解只与之前的阶段有关，与后面的阶段无关。

1 算法解释

对于序列比对问题，我们想知道两个序列之间的存在的差异，首先我们想了解全局比对的差异。 对于全局比对，前人开发出Needleman-Wunsch算法，该算法是基于动态规划的全局比对算法。

\[ F(i, j)=max \left\{ \begin{aligned} F(i-1, j-1) + s(x_i, y_j) \\ F(i-1, j) + s(x_i, -) \\ F(i, j-1) + s(-, y_j) \end{aligned} \right. \]

\(F(i-1, j-1) + s(x_i, y_j)\) 前半部分\(F(i-1, j-1)\)表示前面\(i-1\)和\(j-1\)比对计算完成的情况下， 比对两条序列的当前碱基\(x_i\)和\(y_j\), \(s(x_i, y_j)\)表示当前的替换分数。

\(F(i-1, j) + s(x_i, -)\) 前半部分\(F(i-1, j)\)表示前面\(i-1\)和\(j\)比对计算完成的情况下， 比对\(x\)序列的当前碱基\(x_i\)和\(y\)序列空白\(-\), \(s(x_i, -)\)表示\(x\)序列当前碱基和gap的替换分数。

\(F(i, j-1) + s(-, y_j)\) 前半部分\(F(i, j-1)\)表示前面\(i\)和\(j-1\)比对计算完成的情况下，
比对\(y\)序列的当前碱基\(y_j\)和\(x\)序列空白\(-\), \(s(-, y_j)\)表示\(y\)序列当前碱基和gap的替换分数。

2 罚分矩阵

常见的罚分矩阵有几种：

  • 等价矩阵：碱基相同为1，不同为0
  • 转换-颠换矩阵：嘌呤(AG)之间为转换，嘧啶(TC)之间为转换，嘧啶和嘌呤之间为颠换
  • BLAST矩阵：碱基相同为5，碱基不同为-4

NCBI Blast matrices

这里实现的代码采用自定的一个矩阵，具体如下：
对角线为自身比对上的分数0，即match。
嘌呤A/G之间比对分数为2，即mismatch。
嘧啶T/C之间比对分数为2，即mismatch。
嘧啶(T/C)和嘌呤(A/G)之间的比对分数为4, 即mismatch。
空白罚分为8，即gap。
Gap一般分2种:
Gap opening(只有出现一个就罚分)，
Gap extension(在一空位的扩大，拓展空位罚分，长度越长，罚分越多)

1
2
3
4
5
6
7
8
# 因为我们采用的是正数罚分矩阵，
# 所以最佳比对是罚分最小的
   A  T  G  C  -
A  0  4  2  4  8
T  4  0  4  2  8
G  2  4  0  4  8
C  4  2  4  0  8
-  8  8  8  8  8

3 全局比对python实现

1
2
3
4
5
6
7
alphabet = ['A', 'T', 'G', 'C']
score = [[0, 4, 2, 4, 8],
         [4, 0, 4, 2, 8],
         [2, 4, 0, 4, 8],
         [4, 2, 4, 0, 8],
         [8, 8, 8, 8, 8]]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import numpy as np 

def globalAlignment(x, y):
    """
    Needleman-Wunsch algorithm
    x : (str) sequence a
    y : (str) sequence b
    """
    # 创建一个距离矩阵，记录比对过程
    D = []
    # 初始化矩阵，全部填充0
    for i in range(len(x)+1):
        D.append([0] * (len(y)+1))
     
    
    # 初始化第一列，都填上空白的罚分，随着位置累加
    for i in range(1, len(x)+1):
        # 二维数组的第二索引为0
        D[i][0] = D[i-1][0] + score[alphabet.index(x[i-1])][-1]
        
    # 初始化第一行,都填上空白的罚分，随着位置累加
    for j in range(1, len(y)+1):
        # 二维数组的第一索引为0 
        D[0][j] = D[0][j-1] + score[-1][alphabet.index(y[j-1])]
    
    # 填充剩余的矩阵 
    for i in range(1, len(x)+1):
        for j in range(1, len(y)+1):
            # 水平距离，x(i)和y(j-1)比对上，yj位置gap罚分
            # D[i][j-1]未知，后一项罚分已知
            distHor = D[i][j-1] + score[-1][alphabet.index(y[j-1])]
            # 垂直距离，x(i-1)和y(j)比对上，xi位置gap罚分
            # D[i-1][j]未知，后一项罚分已知
            distVer = D[i-1][j] + score[alphabet.index(x[i-1])][-1]
            # 对角线距离，x(i-1)和y(j-1)比对上，xi和yj错配罚分
            # D[i-1][j-1]未知，后一项罚分已知
            distDiag = D[i-1][j-1] + score[alphabet.index(x[i-1])][alphabet.index(y[j-1])]
            # 取最小值的比对路线，记录比对值, 该值在后续上面的三个距离计算中重复使用，避免循环暴力递归
            # 这里才是动态规划算法的核心所在。
            D[i][j] = min(distHor, distVer, distDiag) # 即问题分解成上述三个子问题，罚分最小才是最佳比对,视罚分函数而定
    # 返回右下角的值, 最佳的比对值
    return D, D[-1][-1] 


1
2
3
4
5
6
7
8
x = 'TATGTCATGC'
y = 'AATGTCATGC'
D, value = globalAlignment(x,y)
import numpy as np 
d = np.array(D)
np.diagonal(d)
print(value)
print(D)

4 局部比对

在生物学序列比对中，由于内含子或者motif的因素的存在，我们常常需要算出局部最佳的比对情况， 类似最大子串问题，因此需要局部比对，前人实现了smith-waterman算法来解决这个问题。 和全局比对相比，局部比对在算法上只有一点不同。

\[ F(i, j)= max \left\{ \begin{aligned} F(i-1, j-1) + s(x_i, y_j) \\ F(i-1, j) + s(x_i, -) \\ F(i, j-1) + s(-, y_j) \\ 0 \end{aligned} \right. \]

5 局部比对python实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
import numpy

def costFunction(xc, yc):
    ''' 自定义一个罚分矩阵 
    2 to match, -6 to gap, -4 to mismatch 
    '''
    # 比对上match 
    if xc == yc:
        return 2
    # 空白 gap
    if xc == '-' or yc == '-': 
        return -6
    # 错配 mismatch 
    return -4

def localAlignment(x, y, s):
    '''
    Smith-Waterman algorithm
    Calculate local alignment values of sequences x and y using
    dynamic programming.  Return maximal local alignment value. 
    '''
    # 初始化二维数组记录子问题的解
    D = numpy.zeros((len(x)+1, len(y)+1), dtype=int)
    # 开始循环比对
    for i in range(1, len(x)+1):
        for j in range(1, len(y)+1):
            # 这里选最大值是因为match的罚分是正数, gap和mismatch的是负数
            D[i, j] = max(D[i-1, j-1] + s(x[i-1], y[j-1]), # 对角线，  前i-1和j-1比对完成, 计算i,j位置的罚分
                          D[i-1, j  ] + s(x[i-1], '-'),    # 垂直方向，前i-1和j比对完成，计算i和gap的罚分
                          D[i  , j-1] + s('-',    y[j-1]), # 水平方向，前i和j-1比对完成，计算gap和j的罚分
                          0)                               # 空值
    argmax = numpy.where(D == D.max()) # 最大值的索引
    return D, int(D[argmax])

x, y = 'GGTATGCTGGCGCTA', 'TATATGCGGCGTTT'
D, best = localAlignment(x, y, exampleCost)
print(D)
print("Best score=%d, in cell %s" % (best, numpy.unravel_index(numpy.argmax(D), D.shape)))