Needleman-Wunsch 算法及实现

Needleman-Wunsch (NW) 算法常用于基因序列匹配以及单词匹配，是一种实现全局最优匹配的动态规划算法。其最差情况下的时间复杂度是 $O(mn)$ , 空间复杂度是 $O(mn)$ ，其中 $m$ 和 $n$ 是两个序列的长度。

不妨设有两个序列：
GCATGCU
GATTACA

NW算法的具体步骤如下：

构造一个如下形式的表格。
设计得分矩阵。例如：第 $i$ 行第 $j$ 列字母匹配+1，不匹配-1，插入和删除（字母与空白对比）操作-1。
将以上表格的第二行第二列的初始得分设为0，通过公式 $D(i,j)=\max\{D(i-1，j)+score(delete/insert), D(i,j-1)+score(delete/insert), D(i-1,j-1)+score(match/mismatch)\}$ 计算填充整个表格，并记录得分的方向（即通过哪个操作得到的）。特别的，当限制每个网格的最低分为0时，此时NW算法就变成了寻找局部最优匹配的Smith-Waterman算法
从表格的右下角开始，根据之前记录的路径逐级返回，并通过对应的操作得到最优匹配的序列。如下图所示：

得到的最优匹配为：
GCATGCU：GCATG-CU / GCAT-GCU / GCA-TGCU
GATTACA ：G-ATTACA / G-ATTACA / G-ATTACA

Python 实现

# -*- coding: utf-8 -*-
import itertools
import copy
seq1 = '*GCATGCU'
seq2 = '*GATTACA'

#定义得分矩阵
score_dic = {'match': 1, 'mismatch': -1, 'gap': -1}
# initial score_matrix
score_matrix = [[0 for column in range(len(seq1))] for row in range(len(seq2))]
trace_back = [[[]for column in range(len(seq1))] for row in range(len(seq2))]
for i in range(len(score_matrix[0])):
    score_matrix[0][i] = i * -1
    if i > 0:
        trace_back[0][i].append('left')
for i in range(len(score_matrix)):
    score_matrix[i][0] = i * -1
    if i > 0:
        trace_back[i][0].append('up')
trace_back[0][0].append('done')
# fill the table
for i in range(1, len(score_matrix)):
    for j in range(1, len(score_matrix[0])):
        if seq1[j] == seq2[i]:
            char_score = score_dic['match']
        else:
            char_score = score_dic['mismatch']
        top_score = score_matrix[i - 1][j] + score_dic['gap']
        left_score = score_matrix[i][j - 1] + score_dic['gap']
        diag_score = score_matrix[i - 1][j - 1] + char_score
        score = max(top_score, left_score, diag_score)
        score_matrix[i][j] = score
        if top_score == score:
            trace_back[i][j].append('up')
        if left_score == score:
            trace_back[i][j].append('left')
        if diag_score == score:
            trace_back[i][j].append('diag')

# 根据结果计算最优匹配的序列
# pointer = [seq2_index, seq1_index]
pointer = [len(score_matrix) - 1, len(score_matrix[0]) - 1]
align_seq1 = []
align_seq2 = []
arrow = trace_back[pointer[0]][pointer[1]]


def seq_letter_finder(current_arrow, current_pointer):
    if current_arrow == 'diag':
        letter = [seq1[current_pointer[1]], seq2[current_pointer[0]]]
        next_pointer = [current_pointer[0] - 1, current_pointer[1] - 1]
        next_arrow = trace_back[next_pointer[0]][next_pointer[1]]
        return letter, next_arrow, next_pointer
    elif current_arrow == 'left':
        letter = [seq1[current_pointer[1]], '-']
        next_pointer = [current_pointer[0], current_pointer[1] - 1]
        next_arrow = trace_back[next_pointer[0]][next_pointer[1]]
        return letter, next_arrow, next_pointer
    else:
        letter = ['-', seq2[current_pointer[0]]]
        next_pointer = [current_pointer[0] - 1, current_pointer[1]]
        next_arrow = trace_back[next_pointer[0]][next_pointer[1]]
        return letter, next_arrow, next_pointer


def align_seq_finder(rec_arrow, rec_pointer, rec_ls):
    if rec_arrow[0] == 'done':
        rec_ls = [rec_ls[0][::-1], rec_ls[1][::-1]]
        return rec_ls
    else:
        if len(rec_arrow) == 1:
            letter, rec_arrow, rec_pointer = seq_letter_finder(rec_arrow[0], rec_pointer)
            rec_ls[0] += letter[0]
            rec_ls[1] += letter[1]
            return align_seq_finder(rec_arrow, rec_pointer, rec_ls)

        elif len(rec_arrow) == 2:
            arrow1 = copy.deepcopy(rec_arrow[0])
            pointer1 = copy.deepcopy(rec_pointer)
            ls1 = copy.deepcopy(rec_ls)
            arrow2 = copy.deepcopy(rec_arrow[1])
            pointer2 = copy.deepcopy(rec_pointer)
            ls2 = copy.deepcopy(rec_ls)
            letter1, arrow1, pointer1 = seq_letter_finder(arrow1, pointer1)
            letter2, arrow2, pointer2 = seq_letter_finder(arrow2, pointer2)
            ls1[0] += letter1[0]
            ls1[1] += letter1[1]
            ls2[0] += letter2[0]
            ls2[1] += letter2[1]
            return list(itertools.chain(align_seq_finder(arrow1, pointer1, ls1),
                                        align_seq_finder(arrow2, pointer2, ls2)))
        else:
            arrow1 = copy.deepcopy(rec_arrow[0])
            pointer1 = copy.deepcopy(rec_pointer)
            pointer2 = copy.deepcopy(rec_pointer)
            pointer3 = copy.deepcopy(rec_pointer)
            ls1 = copy.deepcopy(rec_ls)
            ls2 = copy.deepcopy(rec_ls)
            ls3 = copy.deepcopy(rec_ls)
            letter, arrow1, pointer1 = seq_letter_finder(arrow1, pointer1)
            ls1[0] += letter[0]
            ls1[1] += letter[1]
            arrow2 = rec_arrow[1]
            letter, arrow2, pointer2 = seq_letter_finder(arrow2, pointer2)
            ls2[0] += letter[0]
            ls2[1] += letter[1]
            arrow3 = rec_arrow[2]
            letter, arrow3, pointer3 = seq_letter_finder(arrow3, pointer3)
            ls3[0] += letter[0]
            ls3[1] += letter[1]
            return list(itertools.chain(align_seq_finder(arrow1, pointer1, ls1),
                                        align_seq_finder(arrow2, pointer2, ls2),
                                        align_seq_finder(arrow3, pointer3, ls3)))

#打印结果
ls = align_seq_finder(arrow, pointer, ['', ''])
for i in range(0, len(ls), 2):
    print(ls[i], ls[i+1])

运行结果：

GCATG-CU G-ATTACA
GCAT-GCU G-ATTACA
GCA-TGCU G-ATTACA