利用Needleman–Wunsch算法進行DNA序列全局比對

生物信息學原理作業第二彈：利用Needleman–Wunsch算法進行DNA序列全局比對。

具體原理：https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm。

利用Needleman–Wunsch算法進行DNA序列全局比對

轉載請保留出處！

貼上python代碼：

# -*- coding: utf-8 -*-
"""
Created on Sat Nov 25 18:20:01 2017

@author: zxzhu
后需修改：
1.加命令行參數
2.給出多種比對結果
"""

import numpy as np
import pandas as pd
sequence1 = 'AACGTACTCA'
sequence2 = 'TCGTACTCA'
s1 = ''
s2 = ''
gap = -4
score_matrix = pd.read_excel('score.xlsx')        #score matrix
best_matrix = np.empty(shape= (len(sequence2)+1,len(sequence1)+1),dtype = int)

def get_match_score(s1,s2):
    score = score_matrix[s1][s2]
    return score

for i in range(len(sequence2)+1):
    for j in range(len(sequence1)+1):
        if i == 0:
            best_matrix[i][j] = gap * j
        
        elif j == 0:
            best_matrix[i][j] = gap *i
        else:
            match = get_match_score(sequence2[i-1],sequence1[j-1])
            gap1_score = best_matrix[i-1][j]+gap
            gap2_score = best_matrix[i][j-1]+gap
            match_score = best_matrix[i-1][j-1]+match
            best_matrix[i][j] = max(gap1_score,gap2_score,match_score)
print(best_matrix)
i,j = len(sequence2),len(sequence1)
while(i>0 or j>0):
    match = get_match_score(sequence2[i-1],sequence1[j-1])
    if i>0 and j>0 and best_matrix[i][j] == best_matrix[i-1][j-1]+match:
        s1 += sequence1[j-1]
        s2 += sequence2[i-1]
        i-=1;j-=1
    elif i>0 and best_matrix[i,j] == best_matrix[i-1,j]+gap:
        s1+='-'
        s2+=sequence2[i-1]
        i-=1
    else:
        s1+=sequence1[j-1]
        s2+='-'
        j-=1
print(s1[::-1]+'\n'+s2[::-1])

后面會加入命令行。

多種結果這里只取了一種，這個問題有待解決。

如果有其他的方法我會及時添加。