Needleman–Wunsch 算法的代码实现

# -*- coding: utf-8 -*-
"""
    :Author: huangsh
    :Date: 19-7-28 下午19:17
    :Description: 使用bidu Needleman–Wunsch 算法来计算两条序列的最大相似得分
    如果您对此算法不熟悉,可以去看看我写的一篇拙文:https://www.jianshu.com/p/002bbebcaaef
"""

from collections import namedtuple

F = namedtuple('F', ('score', 'pointer'))

## 初始化二维矩阵, # 生成 x行,y列的二维矩阵,初始化第0行,0列的元素
def init_array(x, y):
    array = [[0] * (y) for _ in range(x)]
    array[0][0] = F(0, None)
    for j in range(1, y):
        array[0][j] = F((-5)*j, [0, j-1])
    for i in range(1, x):
        array[i][0] = F((-5)*i, [i-1, 0])
    return array

## 一行一行的计算矩阵中的每个各自中的最优结果。当前格子中的最优结果由它的三个来源推出
def compute(array, seq1, seq2):
    row, col = len(seq2), len(seq1)
    for i in range(1, row+1):
        for j in range(1, col+1):
            if seq1[j-1] == seq2[i-1]:  # 这里简化了得分矩阵,完全匹配得10分,不完全得5分,有gap减5分
                s = 10
            else:
                s = 5
            lu = [array[i-1][j-1].score+s, [i-1, j-1]] # idx 0:最大得分,idx 1:来源坐标
            left = [array[i-1][j].score-5, [i-1, j]]
            up = [array[i][j-1].score-5, [i, j-1]]
            max_choice = max([lu,left, up], key=lambda x: x[0])
            score= max_choice[0]
            pointer = max_choice[1]
            array[i][j] = F(score, pointer)  # 在当前保存最大得分,和来源坐标,方便回溯。
    return array

## 回溯。从(m,n)一直回溯到(0,0)
def backtrack(array, seq1, seq2):
    s1 = []
    s2 = []
    row, col = len(seq2), len(seq1)
    while array[row][col].score != 0:
        i, j = array[row][col].pointer # pointer 指向来源方的坐标
        if i+1 == row and j+1 == col: # 左上方
            s1.append(seq1[col-1])
            s2.append(seq2[row-1])
            row, col = i, j
        elif row == i+1 and col == j: # 来源:上方
            s1.append("-")
            s2.append(seq2[i])
            row, col = i, j
        elif row == i and col == j+1: # 左方
            s1.append(seq1[j])
            s2.append("-")
            row, col = i, j
    s1 = ''.join(s1[::-1])  #因为是从最后往前回溯的,需要将逆转一下list
    s2 = ''.join(s2[::-1])
    return s1, s2

def main(seq1, seq2):
    x, y = len(seq2)+1 , len(seq1)+1 # x是矩阵行数,y是矩阵列数

    array = init_array(x, y)
    array = compute(array, seq1, seq2)
    s1, s2 = backtrack(array, seq1, seq2)
    max_score = array[x-1][y-1].score

    print("最大得分:", max_score)
    print(s1)
    print(s2)

if __name__ == '__main__':
    seq1 = "ATCGCGCAACTGCGCGC"
    seq2 = "ACGCGCACTGCGGC"
    main(seq1, seq2)

原文地址:https://www.cnblogs.com/huanping/p/11273391.html

时间: 2024-10-03 10:41:06

Needleman–Wunsch 算法的代码实现的相关文章

文本比较算法Ⅱ——Needleman/Wunsch算法

在"文本比较算法Ⅰ--LD算法"中介绍了基于编辑距离的文本比较算法--LD算法. 本文介绍基于最长公共子串的文本比较算法--Needleman/Wunsch算法. 还是以实例说明:字符串A=kitten,字符串B=sitting 那他们的最长公共子串为ittn(注:最长公共子串不需要连续出现,但一定是出现的顺序一致),最长公共子串长度为4. 定义: LCS(A,B)表示字符串A和字符串B的最长公共子串的长度.很显然,LSC(A,B)=0表示两个字符串没有公共部分. Rev(A)表示反转

利用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

文本比较算法:Needleman/Wunsch算法

本文介绍基于最长公共子序列的文本比较算法--Needleman/Wunsch算法.还是以实例说明:字符串A=kitten,字符串B=sitting那他们的最长公共子序列为ittn(注:最长公共子序列不需要连续出现,但一定是出现的顺序一致),最长公共子序列长度为4. 和LD算法类似,Needleman/Wunsch算法用的都是动态规划的思想,两者十分相似. 举例说明:A=GGATCGA,B=GAATTCAGTTA,计算LCS(A,B). 第一步:初始化动态转移矩阵 Needleman/Wunsch

Python实现各种排序算法的代码示例总结

Python实现各种排序算法的代码示例总结 作者:Donald Knuth 字体:[增加 减小] 类型:转载 时间:2015-12-11我要评论 这篇文章主要介绍了Python实现各种排序算法的代码示例总结,其实Python是非常好的算法入门学习时的配套高级语言,需要的朋友可以参考下 在Python实践中,我们往往遇到排序问题,比如在对搜索结果打分的排序(没有排序就没有Google等搜索引擎的存在),当然,这样的例子数不胜数.<数据结构>也会花大量篇幅讲解排序.之前一段时间,由于需要,我复习了

排序算法总结---代码+性能

// data_sort_alg.cpp : 定义控制台应用程序的入口点. // #include "stdafx.h" #include "sort_alg.h" #include <iostream> #include <vector> void show(std::vector<int> &a) { std::vector<int>::iterator it=a.begin(); while(it!=a.

机器学习算法的代码实现之第四章节:回归之梯度上升法

二种类别的点在平面上分布,我想找到一条直线,将平面划为两半边,每一边的点类别尽可能的统一,如何找到效果最佳的分界线,这就是最佳拟合问题,也叫作回归问题. 这次,代码很少.logRegres.py # coding:utf-8 from numpy import * #=============================================================================== # 数据集 #=============================

计算机视觉算法与代码集锦

计算机视觉算法与代码集锦 计算机视觉是结合了传统摄影测量,现代计算机信息技术.人工智能等多学科的一个大学科,是一片开垦不足的大陆,路很远,但很多人都在跋涉! 本文转自CSDN(地址http://blog.csdn.net/whucv/article/details/7907391),是一篇很好的算法与代码总结文档,转载在此供大家学习参考. 原文如下: UIUC的Jia-Bin Huang同学收集了很多计算机视觉方面的代码,链接如下: https://netfiles.uiuc.edu/jbhua

谱聚类算法及其代码(Spectral Clustering)

简介 文章将介绍谱聚类(spectral clustering)的基本算法,以及在matlab下的代码实现.介绍内容将包括: 从图分割角度直观理解谱聚类 谱聚类算法步骤 数据以及实现代码 本文将不会涉及细节化的证明和推导,如有兴趣可参考july大神的文章从拉普拉斯矩阵说到谱聚类. 对谱聚类的理解 这一节将从图分割聚类的角度直观理解谱聚类.不过,因为本人是从事社交媒体分析的,将从一种社会关系网络的角度来介绍网络图分割成多个子图的概念. 图的分割 首先将社会关系网络看成是一个整体,每一个个体(use

php短网址算法实例代码分享

php实现的短网址算法,理论上支持1,073,741,824个短网址. 每个网址用6个字符代替,(6^32) 最多可以拥有1,073,741,824个短网址.当然,你还可以记录更详细的信息,如访问记录,创建时间等.如果真不够用了,还可以删掉很久不用的. function shorturl($input) { $base32 = array ( 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p