生信自學筆記（六）Needleman-Wunsch算法的python實作

Needleman-Wunsch算法是一種全局聯配算法，原理說起來不算難，就是首先把兩條目标DNA序列的最前面加一個空格代表空位，然後畫一個大小為（一條DNA序列長度*另一條DNA序列的長度）的矩陣，矩陣橫縱坐标分别代表對應序列的堿基或是空位，這樣就可以依據打分規則把這個矩陣填充起來。

而矩陣的初始化方法是這樣的，首先填完第一行和第一列，也就是簡單的其中一個是空位的配對情況，然後依次填第二行第二列、第三行第三列。每次填新的分數時，首先要選取目标位置上側、左上、左側的分數，然後計算從這幾個方向趨近分别會再得多少分，分别和原來的分數相加得到一個估計分數，最後填上最大的估計分數或者是最少的懲罰值。

得分矩陣畫好之後，就可以回溯聯配了。打分我們是從左上到右下，回溯就是從右下到左上。每一次都要選擇這樣的位置，即，從這個位置走到目前位置，按照打分規則恰好可以獲得目前位置現在的分數，這就是十分典型的深搜問題了。此外，由于回溯時可以有多個位置能夠達到相同的分數，是以最終的聯配結果也是有多個的，當然它們的得分也是一緻的。

樣例用的是教科書上的:CTGTATC vs CTATAATCCC

計分規則：堿基比對不罰分，錯配罰一分，引入空位罰三分。

計算出的得分矩陣如下

[[ 0.  3.  6.  9. 12. 15. 18. 21. 24. 27. 30.]
 [ 3.  0.  3.  6.  9. 12. 15. 18. 21. 24. 27.]
 [ 6.  3.  0.  3.  6.  9. 12. 15. 18. 21. 24.]
 [ 9.  6.  3.  1.  4.  7. 10. 13. 16. 19. 22.]
 [12.  9.  6.  4.  1.  4.  7. 10. 13. 16. 19.]
 [15. 12.  9.  6.  4.  1.  4.  7. 10. 13. 16.]
 [18. 15. 12.  9.  6.  4.  2.  4.  7. 10. 13.]
 [21. 18. 15. 12.  9.  7.  5.  3.  4.  7. 10.]]
           

DNA比對結果

源碼:

import numpy as np

# 對任意兩個字元進行比較
def compare(m, n, match, gap, n_match):
    if m == n:
        return match
    elif " " == m or n == " ":
        return gap
    else:
        return n_match
#建構打分矩陣
def mark(seq1, seq2, match, n_match, gap):
    "用來生成最終的打分矩陣,傳入形參分别為兩個序列和比對罰分、錯配罰分、引入空位罰分"
    a = len(seq1)
    b = len(seq2)
    # 初始化打分矩陣
    S = np.zeros((a + 1, b + 1))
    # 為打分矩陣的第一行和第一列指派
    S[0, :] = np.fromfunction(lambda x, y: gap * (x + y), (1, b + 1), dtype=int)
    S[:, 0] = np.fromfunction(lambda x, y: gap * (x + y), (1,a + 1), dtype=int)
    # 在字元串最前添加前導空格
    seq1 = " " + seq1[:]
    seq2 = " " + seq2[:]
    # 三目運算符要這樣寫
    for r in range(1, b+1 if a>b else a+1):
        # up = S[r-1,r-1:]
        # left = S[r-1:,r-1]
        for c in range(r, b+1):
            # 分别獲得上方元素、對角元素、左側元素的分數, 上方和左側補了空格之後恒為3
            last_score = [S[r - 1, c], S[r - 1, c - 1], S[r, c - 1]]
            change_score = [gap,compare(seq1[r],seq2[c], match, gap, n_match),gap]
            #兩個數組中的元素相加要這麼寫
            new_score = np.add(last_score, change_score)
            # print(last_score)
            # print(change_score)
            # print(new_score)
            # print("r = ",r," c = ", c)
            S[r,c] = min(new_score)
            # print(S[r,c])
        for c in range(r+1, a+1):
            # 分别獲得上方元素、對角元素、左側元素的分數
            last_score = [S[c-1,r], S[c-1, r-1], S[c, r - 1]]
            change_score = [gap,compare(seq1[c],seq2[r], match, gap, n_match),gap]
            new_score = np.add(last_score, change_score)
            S[c,r] = min(new_score)
    return S
#回溯
def traceback(seq1, seq2, match, n_match, gap, S, current_x, current_y, S1, S2, t1, t2):
    if current_x == 1 and current_y == 1:
        S1.append(seq1[1] + t1[:])
        S2.append(seq2[1] + t2[:])
        return
    if S[current_x,current_y]-S[current_x-1,current_y] == gap:
        traceback(seq1, seq2, match, n_match, gap, S, current_x-1, current_y, S1, S2, seq1[current_x]+t1[:], "-"+t2[:])
    if S[current_x,current_y]-S[current_x,current_y-1] == gap:
        traceback(seq1, seq2, match, n_match, gap, S, current_x, current_y-1, S1, S2, "-" + t1[:], seq2[current_y] + t2[:])
    if S[current_x,current_y]-S[current_x-1,current_y-1] == compare(seq1[current_x], seq2[current_y], match, gap, n_match):
        traceback(seq1, seq2, match, n_match, gap, S, current_x-1, current_y-1, S1, S2, seq1[current_x] + t1[:], seq2[current_y] + t2[:])
           

調用

F = mark("CTGTATC","CTATAATCCC",0,1,3)
S1 = []
S2 = []
t1 = ""
t2 = ""
traceback(" CTGTATC", " CTATAATCCC", 0, 1, 3, F, F.shape[0]-1, F.shape[1]-1, S1, S2, t1, t2)
print(S1)
print(S2)
           

學了自然語言處理之後，換了個思路寫，好像更加簡單一些。。。不過缺點是不能得到所有可行解了

# -*- coding: UTF-8 -*-
import numpy as np
def dynamic_edit(str1, str2):
    len_str1 = len(str1) + 1
    len_str2 = len(str2) + 1
    newstr1 = newstr2 = ""
    matrix = np.zeros([len_str1, len_str2])
    matrix[0] = np.arange(len_str2)
    matrix[:, 0] = np.arange(len_str1)
    for i in range(1, len_str1):
        for j in range(1, len_str2):
            matrix[i][j] = min(matrix[i - 1][j] + 1, matrix[i][j - 1] + 1,
                               matrix[i - 1][j - 1] + 2*(int)(str1[i - 1] != str2[j - 1]))
    i = len_str1-1
    j = len_str2-1
    while i != 0 or j != 0:
        if (matrix[i][j] == matrix[i - 1][j - 1] and str1[i - 1] == str2[j - 1]) or matrix[i][j] == matrix[i - 1][j - 1] + 2:
            newstr1 = str1[i - 1].__add__(newstr1)
            newstr2 = str2[j - 1].__add__(newstr2)
            i = i-1
            j = j-1
        elif matrix[i][j] == matrix[i - 1][j] + 1:
            newstr1 = str1[i - 1].__add__(newstr1)
            newstr2 = "*".__add__(newstr2)
            i = i-1
        elif matrix[i][j] == matrix[i][j - 1] + 1:
            newstr1 = "*".__add__(newstr1)
            newstr2 = str2[j - 1].__add__(newstr2)
            j = j-1
    print(matrix)
    print(newstr1)
    print(newstr2)


str1 = "EXECUTION"
str2 = "INTENTION"
dynamic_edit(str1, str2)