字符串匹配的基本算法以及动态规划算法简析

新的一期算法课作业是DNA匹配。DNA可以看作一串字符串，于是可以转化为字符串的匹配问题。字符串的匹配并不容易，因为这是一个计算量及其巨大的工程。比如“abc”和“ac”的匹配，最佳匹配应该是“abc”和“a-c”，可是为了计算这么短的字符串，就必须要考虑到很多种情况。

于是引入动态规划算法。Dynamic Programming 是一个对于重复计算优化的算法。比如一个问题，分拆成几个小问题，然后再分成几个小问题，如果小问题的计算是重复的，就可以直接代入即可。比如计算匹配的矩阵的时候：

def compute_alignment_matrix(seq_x, seq_y, scoring_matrix, global_flag):

    length_x = len(seq_x)
    length_y = len(seq_y)
    alignment_matrix = [[0 for dummy_col in range(length_y + 1)] for dummy_row in range(length_x + 1)]
    
    for idx in range(1, length_x+1):
        alignment_matrix[idx][0] = alignment_matrix[idx-1][0] + scoring_matrix[seq_x[idx-1]]['-']
    for idx in range(1, length_y+1):
        alignment_matrix[0][idx] = alignment_matrix[0][idx-1] + scoring_matrix['-'][seq_y[idx-1]]
    for idx_x in range(1, length_x+1):
        for idx_y in range(1, length_y+1):
            alignment_matrix[idx_x][idx_y] = max(
                alignment_matrix[idx_x-1][idx_y-1] + scoring_matrix[seq_x[idx_x-1]][seq_y[idx_y-1]],
                alignment_matrix[idx_x-1][idx_y] + scoring_matrix[seq_x[idx_x-1]]['-'],
                alignment_matrix[idx_x][idx_y-1] + scoring_matrix['-'][seq_y[idx_y-1]],
            )
                
    return alignment_matrix

scoring_matrix是得分矩阵，比如A与C匹配得分为－5，A与A匹配得分10，诸如此类。

这段代码主要是用枚举的方式，把所有可能的组合都计算一遍，得出的分数放到矩阵相应的位置上。前面计算的数值可以给后面的计算用。

注意看其中这样一段：

for idx in range(1, length_x+1):
    alignment_matrix[idx][0] = alignment_matrix[idx-1][0] + scoring_matrix[seq_x[idx-1]]['-']

比如计算(1, 0)位置上的数值，其实就是(0, 0) 加上当前位置的得分。然后(2, 0) 就是(1, 0)加上当前位置的得分。这样计算量就大大减轻了。因为所有的子问题都是已经解决的问题。y方向上的计算与x方向是一致的，中间区域则是超各个方向计算，然后取其中的最大值。

这个过程有点像选择图片区域的时候，从左上角开始滑，网格慢慢覆盖到右下角。

得到的是一个得分矩阵，我们从中选择最大值，其实就是最佳匹配，我们把上述过程推回去，就可以分别得到两个匹配的字符串。

def compute_global_alignment(seq_x, seq_y, scoring_matrix, alignment_matrix):
 
    score_optimal = max([max(item) for item in alignment_matrix])
    for item in alignment_matrix:
        if score_optimal in item:
            idx_y = item.index(score_optimal)
            idx_x = alignment_matrix.index(item)
    
    #initial the align_x and align_y, reserve it 
    align_x = seq_x[idx_x:]
    align_y = seq_y[idx_y:]
    align_x = align_x[::-1]
    align_y = align_y[::-1]

    while idx_x and idx_y:
        if alignment_matrix[idx_x][idx_y] == alignment_matrix[idx_x-1][idx_y-1] + scoring_matrix[seq_x[idx_x-1]][seq_y[idx_y-1]]:
            align_x += seq_x[idx_x-1]
            align_y += seq_y[idx_y-1]
            idx_x -= 1
            idx_y -= 1
        else:
            if alignment_matrix[idx_x][idx_y] == alignment_matrix[idx_x-1][idx_y] + scoring_matrix[seq_x[idx_x-1]]['-']:
                align_x += seq_x[idx_x-1]
                align_y += '-'
                idx_x -= 1
                #import pdb; pdb.set_trace()
            else:
                align_x += '-'
                align_y += seq_y[idx_y-1]
                idx_y -= 1
    while idx_x:
        align_x += seq_x[idx_x-1]
        align_y += '-'
        idx_x -= 1
    while idx_y:
        align_x += '-'
        align_y += seq_y[idx_y-1]
        idx_y -= 1
        
    #reserve the string
    align_x = align_x[::-1]
    align_y = align_y[::-1]
        
    # align_x and align_y have same length
    length_x = len(align_x)
    length_y = len(align_y)
    if length_x < length_y:
        align_x += '-' * (length_y - length_x)
    else:
        align_y += '-' * (length_x - length_y)
        
    #The score need to recalculate, cause there may be '-' more
    score_updated = sum([scoring_matrix[align_x[idx]][align_y[idx]] for idx in range(len(align_x))])

    return score_updated, align_x, align_y

值得注意的是，匹配出来的字符串是反序的，因此我在这个模块中两次颠倒字符串。这个模块看起来很长，其中核心就是那个while循环，把整个过程返回去，得到匹配的字符串。

在while之前，先求出矩阵的最大值和最大值的Index，然后把不进入while循环的字符先取出来。

不进入while循环的字符是匹配的最后，全部保持原有字符的部分。

while循环之后，是把还有不足的头部（生成时的尾部）加入"-"。

最后因为加入了"-"，因此得分需要重新计算。这样全局匹配就完成了。

与全局匹配相对应的是局部匹配，指的是在局部形成最大值。比如"tttabcttt"与"ac"的匹配，最佳应该是

tttabcttt
   a-c

局部匹配与全局匹配的原理基本相同，只不过计算得分矩阵的时候，最小值为0，以及去掉首位的 '-' 。具体还是看源代码吧。

Dynamic Programming的应用范围不像迭代那么通用，但是效率的提升是很明显。字符串的匹配，计算相似度，拼写检查都是可以用的。我会在这个系列的博客的最后一篇构想一个计算人与人相似度的算法。

源代码请参照：Github