蛋白质一级结构全局比对 Needleman-Wunsch 算法的 Python 实现

有关全局对比算法的基本思想可以看这篇文章，这里仅仅用 Python 浅浅实现一下。程序在最后。禁止转载。以下面这两个短肽为例说明一下程序的思路：

SIVK
SSIIVP

在这里我们设计了两个矩阵：matrix1为得分矩阵，用来计算和记录两序列各个配对的得分；matrix2为溯回矩阵，用于记录各配对得分的来源，以及最终对比结果的输出。

Initialize_Matrix部分初始化上面两个矩阵，初始化结果如下：

[[  0.   0.   0.   0.   0.   0.][  0.   0.  -3.  -6.  -9. -12.][  0.  -3.   0.   0.   0.   0.][  0.  -6.   0.   0.   0.   0.][  0.  -9.   0.   0.   0.   0.][  0. -12.   0.   0.   0.   0.][  0. -15.   0.   0.   0.   0.][  0. -18.   0.   0.   0.   0.]] [[0. 0. 0. 0. 0. 0.][0. 0. 2. 2. 2. 2.][0. 1. 0. 0. 0. 0.][0. 1. 0. 0. 0. 0.][0. 1. 0. 0. 0. 0.][0. 1. 0. 0. 0. 0.][0. 1. 0. 0. 0. 0.][0. 1. 0. 0. 0. 0.]]

Needleman_Wunsch部分中s代表得分(score)，d代表来源(direction)。这一部分需要接收对比序列seq1第j个残基和seq2第i个残基在矩阵中的位置(i,j)，其匹配得分w以及其上/左/左上方的最终得分(su,sl,sul)，根据这些输入及得分规则算出最终得分s和来源d。

blosum为配对分数的一个简易库（因为嫌麻烦所以只写了需要用到的，如果是核酸序列那就更简单了），并且可以查询得到配对分数。

check_d用来将matrix2来进行溯回以输出对比序列。

而历遍seq1和seq2的部分整合到了main()里面，由于matrix1索引为0的行(列)是留出来写序列的（后来发现好像numpy的matrix似乎写不进字符串），索引为1的行(列)是用来空对比的（或者说是gap对比），所以我们seq1第j个残基和seq2第i个残基 (j, i) 对应到matrix里面就成了 [i+2, j+2]，(后面的说明先咕一会

import numpy as npdef Initialize_Matrix(seq1,seq2,gap):matrix1 = np.zeros((len(seq2)+2,len(seq1)+2))matrix2 = np.zeros((len(seq2)+2,len(seq1)+2))matrix1[1, 1:] = np.linspace(0, (len(matrix1[1, 1:]) - 1) * gap, num=len(matrix1[1, 1:]), endpoint=True)matrix1[1:, 1] = np.linspace(0, (len(matrix1[1:, 1]) - 1) * gap, num=len(matrix1[1:, 1]), endpoint=True)matrix2[1, 2:] = 2matrix2[2:, 1] = 1return matrix1,matrix2def Needleman_Wunsch(gap,w,su,sl,sul):s = max(sul+w,su+gap,sl+gap)if sul+w == s:d = 3elif sl+gap == s:d = 2else:d = 1return [s,d]def blosum(res1,res2):blosum_bank = {'AS':1,'AI':-1,'AV':0,'AP':-1,'IS':-2,'II':4,'IV':3,'IP':-3,'VS':-2,'VV':4,'VP':-2,'KS':0,'KI':-3,'KV':-2,'KP':-1,'SP':-1,'SS':4}w = blosum_bank.get(res1+res2,blosum_bank.get(res2+res1))return wdef check_d(matrix2,seq1,seq2):i = -1j = -1alain_seq1 = ''alain_seq2 = ''while i >= -(len(seq2)+1) and j >= -(len(seq1)+1):#print(alain_seq1)if matrix2[i][j] == 3:alain_seq1 = seq1[j] + alain_seq1alain_seq2 = seq2[i] + alain_seq2i -= 1j -= 1elif matrix2[i][j] == 1:alain_seq1 = '-' + alain_seq1alain_seq2 = seq2[i] + alain_seq2i -= 1elif matrix2[i][j] == 2:alain_seq1 = seq1[j] + alain_seq1alain_seq2 = '-' + alain_seq2j -= 1elif matrix2[i][j] == 0:print('done')breakelse:print('wrong')breakreturn alain_seq1,alain_seq2def main():seq1 = 'SIVK'seq2 = 'SSIIVP'gap = -3matrix1,matrix2 = Initialize_Matrix(seq1,seq2,gap)#print(matrix1,'\n',matrix2)for j in range(len(seq1)+len(seq2)):for i in range(j+1):if j-i <= len(seq1)-1 and i <= len(seq2)-1:result = Needleman_Wunsch(gap,blosum(seq1[j-i],seq2[i]),matrix1[i+1][j-i+2],matrix1[i+2][j-i+1],matrix1[i+1][j-i+1])matrix1[i+2][j-i+2] = result[0]matrix2[i+2][j-i+2] = result[1]#print(j-i-1,i-1,seq1[j-i-1],seq2[i-1])#print(matrix1)else:print('wrong')continuealain_seq1, alain_seq2 = check_d(matrix2,seq1,seq2)print(matrix1,'\n\n',matrix2,'\n\n',alain_seq1,'\n',alain_seq2)main()

最终输出结果：

done[[  0.   0.   0.   0.   0.   0.][  0.   0.  -3.  -6.  -9. -12.][  0.  -3.   4.   1.  -2.  -5.][  0.  -6.   1.   2.  -1.  -2.][  0.  -9.  -2.   5.   5.   2.][  0. -12.  -5.   2.   8.   5.][  0. -15.  -8.  -1.   6.   6.][  0. -18. -11.  -4.   3.   5.]] [[0. 0. 0. 0. 0. 0.][0. 0. 2. 2. 2. 2.][0. 1. 3. 2. 2. 2.][0. 1. 3. 3. 3. 3.][0. 1. 1. 3. 3. 2.][0. 1. 1. 3. 3. 2.][0. 1. 1. 1. 3. 3.][0. 1. 1. 1. 1. 3.]] -S-IVK SSIIVP

但是这个方法还存在一些问题，就比如：如果在su,sl,sul计算得到的最总分数里存在两个相同的最大值，也就是存在多条回溯路径，那么该怎么实现多个输出呢，或者说设计一个优先级的判断……

等有兴趣再更

蛋白质一级结构全局比对 Needleman-Wunsch 算法的 Python 实现相关推荐

文本比较算法Ⅱ——Needleman/Wunsch算法
在"文本比较算法Ⅰ--LD算法"中介绍了基于编辑距离的文本比较算法--LD算法. 本文介绍基于最长公共子串的文本比较算法--Needleman/Wunsch算法. 还是以实例说明: ...
文本比较算法--Needleman/Wunsch算法
一.定义: 定义: LCS(A,B)表示字符串A和字符串B的最长公共子串的长度. 很显然,LCS(A,B)=0 表示两个字符串没有公共部分. 例如,字符串A=kitten,字符串B=sitting , ...
生物信息学算法之Python实现|Rosalind刷题笔记：002 中心法则：转录
我在生物信息学:全景一文中,阐述了生物信息学的应用领域非常广泛.但是有一点是很关键的,就是细胞内的生命活动都遵从中心法则,生物信息学很多时候就是在中心法则上做文章: 分子生物学中心法则:DNA --& ...
生物信息学算法之Python实现|Rosalind刷题笔记：003 中心法则：翻译
我在生物信息学:全景一文中,阐述了生物信息学的应用领域非常广泛.但是有一点是很关键的,就是细胞内的生命活动都遵从中心法则,生物信息学很多时候就是在中心法则上做文章: 分子生物学中心法则:DNA --& ...
Gillespie算法的Python简单实现（实例）
Gillespie算法的Python简单实现(实例) 文章目录 Gillespie算法的Python简单实现(实例) 前言一.Gillespie是什么? 二.题目三.代码 1.引入库 2.类定义 ...
基于朴素贝叶斯的垃圾分类算法（Python实现）
有代码和数据集的 https://blog.csdn.net/weixin_33734785/article/details/91428991 附有git库代码的 https://www.cnblog ...
手把手教你在多种无监督聚类算法实现Python（附代码）
来源: 机器之心本文约2704字,建议阅读6分钟. 本文简要介绍了多种无监督学习算法的 Python 实现,包括 K 均值聚类.层次聚类.t-SNE 聚类.DBSCAN 聚类. 无监督学习是一类用于 ...
八大排序算法的 Python 实现
八大排序算法的 Python 实现本文用Python实现了插入排序.希尔排序.冒泡排序.快速排序.直接选择排序.堆排序.归并排序.基数排序. 1.插入排序描述插入排序的基本操作就是将一个数据插入 ...
Fuzzy C Means 算法及其 Python 实现——写得很清楚，见原文
Fuzzy C Means 算法及其 Python 实现转自:http://note4code.com/2015/04/14/fuzzy-c-means-%E7%AE%97%E6%B3%95%E5% ...

蛋白质一级结构全局比对 Needleman-Wunsch 算法的 Python 实现

蛋白质一级结构全局比对 Needleman-Wunsch 算法的 Python 实现相关推荐

最新文章

热门文章