利用Needleman–Wunsch算法进行DNA序列全局比对
Posted orange1002
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了利用Needleman–Wunsch算法进行DNA序列全局比对相关的知识,希望对你有一定的参考价值。
生物信息学原理作业第二弹:利用Needleman–Wunsch算法进行DNA序列全局比对。
具体原理:https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm。
贴上python代码:
1 # -*- coding: utf-8 -*- 2 """ 3 Created on Sat Nov 25 18:20:01 2017 4 5 @author: zxzhu 6 后需修改: 7 1.加命令行参数 8 2.给出多种比对结果 9 """ 10 11 import numpy as np 12 import pandas as pd 13 sequence1 = ‘AACGTACTCA‘ 14 sequence2 = ‘TCGTACTCA‘ 15 s1 = ‘‘ 16 s2 = ‘‘ 17 gap = -4 18 score_matrix = pd.read_excel(‘score.xlsx‘) #score matrix 19 best_matrix = np.empty(shape= (len(sequence2)+1,len(sequence1)+1),dtype = int) 20 21 def get_match_score(s1,s2): 22 score = score_matrix[s1][s2] 23 return score 24 25 for i in range(len(sequence2)+1): 26 for j in range(len(sequence1)+1): 27 if i == 0: 28 best_matrix[i][j] = gap * j 29 30 elif j == 0: 31 best_matrix[i][j] = gap *i 32 else: 33 match = get_match_score(sequence2[i-1],sequence1[j-1]) 34 gap1_score = best_matrix[i-1][j]+gap 35 gap2_score = best_matrix[i][j-1]+gap 36 match_score = best_matrix[i-1][j-1]+match 37 best_matrix[i][j] = max(gap1_score,gap2_score,match_score) 38 print(best_matrix) 39 i,j = len(sequence2),len(sequence1) 40 while(i>0 or j>0): 41 match = get_match_score(sequence2[i-1],sequence1[j-1]) 42 if i>0 and j>0 and best_matrix[i][j] == best_matrix[i-1][j-1]+match: 43 s1 += sequence1[j-1] 44 s2 += sequence2[i-1] 45 i-=1;j-=1 46 elif i>0 and best_matrix[i,j] == best_matrix[i-1,j]+gap: 47 s1+=‘-‘ 48 s2+=sequence2[i-1] 49 i-=1 50 else: 51 s1+=sequence1[j-1] 52 s2+=‘-‘ 53 j-=1 54 print(s1[::-1]+‘\n‘+s2[::-1])
后面会加入命令行。
多种结果这里只取了一种,这个问题有待解决。
如果有其他的方法我会及时添加。
以上是关于利用Needleman–Wunsch算法进行DNA序列全局比对的主要内容,如果未能解决你的问题,请参考以下文章