python 比对匹配_用Python从头实现Needleman-Wunsch序列比对算法
本文将介绍生物信息学中一个基本的序列比对算法,Needleman-Wunsch比对算法(以下简称NW算法),并展示如何编写一个简单的Python程序来实现该算法。全文分为上下两个部分:上篇介绍并阐释NW算法的原理,下篇展示如何用Python实现。如果你对该算法不陌生并且仅对实现方法感兴趣,那么可以直接跳到下篇。
——————————————————— 上篇 ————————————————————
序列比对是生物信息学领域的经典问题,出现了很多算法,也衍生出了一大批生物信息学工具,有些已经逐渐淡出历史舞台,有些依然被频繁使用,比如MUSCLE,BLAST,BWA,Bowtie2等。
两序列间的比对(“一对一”)是序列比对话题中最基本的问题,也是本文将要探讨的主题,更复杂的比对问题都可以由此衍生而来,比如,“多对多”:多重序列比对;“一对多”:目标序列对数据库的比对;“多对一”:大量二代测序短序列比对到人类基因组上。
两序列比对,从解决问题的方向上来说,可下分为全局比对(Global Alignment)和局部比对(Local Alignment)两种。前者试图将两条序列从头到尾整体对齐,而后者则只是针对两条序列中高度相似的局部区域进行比对。这两种方式都各有特点,局部比对的意义在于,在实际研究中,我们往往更关注序列间高度保守的区域,因为这暗示着重要的生物学功能。同样的两条序列,采用全局比对和局部比对都可以,但得到的结果可能差别很大:同样的序列,不同的比对方式。上:全局比对;下:局部比对 ##图片来自网络##
NW算法就是一种经典的全局比对算法,由Saul B. Needleman和Christian D. Wunsch于1970年首次提出,当时是为了解决蛋白序列的比对问题,至今依然在使用。值得一提的是,同时期涌现了许多重要的比对算法,比如Smith-Waterman算法。但NW算法是最早的、也是最经典的一个,所以你可能会在很多生物信息学专著中见到它,比如:Richard Durbin 教授的这本 Biological Sequence AnalysisDavid Mount教授的 Bioinformatics: Sequence and Genome Analysis
下面开始介绍NW算法。
第一步,定一个打分策略。这是基础, 定义了两个碱基在各种比对情况下的分值(包括gap),我们这里将会用到一个简单的:
若是两个碱基一样,即完美匹配,即 +1 分
若是两个碱基不一样,即错配,则 -1 分
若是在任一条链上开一个gap,则 -1 分
我们用 θ(a, b)函数来表示碱基 a 和 b 之间的打分。
第二步,计算分值矩阵。分值矩阵的概念是NW算法的重点,矩阵由格子组成,一个比对方案就是一条从左上角格子移动到右下角格子的路径,一个矩阵的路径包含了两条序列间所有可能的比对情况,我们要做的就是找出这其中的最优的路径,这其中的方法就是根据分值矩阵来搜索。这里我们假设想要比对的序列分别为seq1和seq2:
>seq1
TCATC
>seq2
TCATGGC
在计算矩阵的时候,需要在两条序列之前加上一个 “-”,表示gap。左上角的格子为起点,记为0。接下来看矩阵里剩下的格子,每个格子都由别的格子推算而来,用 C(i, j) 表示第 i 行,第 j 列的格子的得分,计算方法:
其实对应的就是从“上”,”对角线“,”左“ 三个方向分别计算分值,取最大值。第一行和第一列的意义有点特殊,表示序列的起始阶段比对到gap上,因此这里的打分就是开 1 个gap,扣1 分,gap越多,扣的越多。
接下来就是剩下的格子,这就需要比较三个方向计算出的最大值了。以下图为例,要计算T对T这个格子的分值,首先T对T自身是一个完美匹配,+1,从对角线而来就是0+1,得到1;从上面来就是-1-1,得-2;从左边来就是-1-1,得-2。显然取最大值1,同时我们知道了这个T-T位置的1是来自对角线的。取最大值1
剩下的所有格子都可以此类推,每到一个格子干两件事:计算分值,记录路径:红色标示可走的路径
最终计算到了右下角,这里值得注意的是从最左上角的根出发,整个路径其实是一颗树,只有一条能够到达最右下角(严格来说并不是唯一的,这里先不考虑),因此引出接下来的第三步:
第三步,矩阵生成完之后,开始回溯(traceback)。我们需要知道到底最右下角的值是经历怎样的路径计算而来,这路径就是我们需要的比对结果。能走到右下角的路径就是全局比对解
——————————————————— 下篇 ————————————————————
关于Python代码部分,我先做几点说明,考虑到本文定位比较基础,所以:1,不尝试优化代码执行效率,只展示算法逻辑;2,仅使用Python内置的数据类型;3,采用简单的函数式编程,不采用OOP。下面就开始撸代码:
首先,搭个架子:
def main():
pass
return
if _name == '__main__':
main()
我们定义一个 theta() 函数来实现任意两个碱基之间的打分:
def theta(a, b):
if a == '-' or b == '-' or a != b: # gap or mismatch
return -1
elif a == b: # match
return 1
接下来,算打分矩阵。为了方便,这里使用双层字典来构造矩阵,并且同时构造两个矩阵,一个是score_mat,用来记录每个格子的分值;一个是trace_mat,用来记录相应格子里的值是由那个方向计算而来(以后回溯会用到),我们约定:0表示从对角线来,也就是左上角,1表示从左边来,2表示从上面来:
def make_score_matrix(seq1, seq2):
"""return score matrix and map(each score from which direction)0: diagnosis1: up2: left"""
seq1 = '-' + seq1
seq2 = '-' + seq2
score_mat = {}
trace_mat = {}
for i,p in enumerate(seq1):
score_mat[i] = {}
trace_mat[i] = {}
for j,q in enumerate(seq2):
if i == 0: # first row, gap in seq1
score_mat[i][j] = -j
trace_mat[i][j] = 1
continue
if j == 0: # first column, gap in seq2
score_mat[i][j] = -i
trace_mat[i][j] = 2
continue
ul = score_mat[i-1][j-1] + theta(p, q) # from up-left, mark 0
l = score_mat[i][j-1] + theta('-', q) # from left, mark 1, gap in seq1
u = score_mat[i-1][j] + theta(p, '-') # from up, mark 2, gap in seq2
picked = max([ul,l,u])
score_mat[i][j] = picked
trace_mat[i][j] = [ul, l, u].index(picked) # record which direction
return score_mat, trace_mat
写一个辅助函数用来打印矩阵,看看效果(主要为了自己学习和调试用)。这个函数打印 score_mat 和 trace_mat 都可以:
def print_m(seq1, seq2, m):
"""print score matrix or trace matrix"""
seq1 = '-' + seq1; seq2 = '-' + seq2
print()
print(' '.join(['%3s' % i for i in ' '+seq2]))
for i, p in enumerate(seq1):
line = [p] + [m[i][j] for j in range(len(seq2))]
print(' '.join(['%3s' % i for i in line]))
print()
return
运行一下试试(截图的显示效果比较理想):
接下来就是准备回溯,寻找”比对解“了。
但立即就遇到一个现实问题:即便知道了如何寻找比对解,知道了路径,以什么形式返回呢,也就是如何记录比对解呢?这里我采用办法是用一段字符串记录比对解,即 path_code。
path_code是一个字符串,从左至右记录了打分矩阵中路径从左上角起点到右下角终点是如何移动的,0表示沿对角线走,1表示往下走(在seq2上开gap),2表示往右走(在seq1上开gap)。
不忙着寻找path_code,我们先来看看如何根据 path_code 打印比对结果:
def pretty_print_align(seq1, seq2, path_code):
'''return pair alignment result string frompath code: 0 for match, 1 for gap in seq1, 2 for gap in seq2'''
align1 = ''
middle = ''
align2 = ''
for p in path_code:
if p == '0':
align1 = align1 + seq1[0]
align2 = align2 + seq2[0]
if seq1[0] == seq2[0]:
middle = middle + '|'
else:
middle = middle + ' '
seq1 = seq1[1:]
seq2 = seq2[1:]
elif p == '1':
align1 = align1 + '-'
align2 = align2 + seq2[0]
middle = middle + ' '
seq2 = seq2[1:]
elif p == '2':
align1 = align1 + seq1[0]
align2 = align2 + '-'
middle = middle + ' '
seq1 = seq1[1:]
print('Alignment:\n\n' + align1 + '\n' + middle + '\n' + align2 + '\n')
return
运行效果:
那么接下来就只剩最后一个最难,也是最核心的问题了,寻找“比对解”,也就是回溯(trace back)。
这里必须要介绍动态规划算法了。总体上来说,动态规划算法就是把一个大问题分解为一个个可重复的小问题,然后一个一个解决小问题,也就解决了整个大问题。
这里有两个关键点:1,整个问题可以拆分为若干个相似的小问题,能够循环;2, 循环有终点。我们的回溯问题就是一个动态规划问题:从最右下角出发,每到一个格子,查找他的源格子,移动到源格子,再查找源格子的源格子,移动到源格子的源格子,以此类推,直到回到矩阵最左上角的根。记录这一路走来的路径就是答案了,也就是path_code。
def traceback(seq1, seq2, trace_mat):
'''find one optimal traceback path from trace matrix, return path code-!- CAUTIOUS: if multiple equally possible path exits, only return one of them -!-'''
seq1, seq2 = '-' + seq1, '-' + seq2
i, j = len(seq1) - 1, len(seq2) - 1
path_code = ''
while i > 0 or j > 0:
direction = trace_mat[i][j]
if direction == 0: # from up-left direction
i = i-1
j = j-1
path_code = '0' + path_code
elif direction == 1: # from left
j = j-1
path_code = '1' + path_code
elif direction == 2: # from up
i = i-1
path_code = '2' + path_code
return path_code
至此整个问题就都解决了,完整可运行的脚本:
import sys
def theta(a, b):
if a == '-' or b == '-' or a != b: # gap or mismatch
return -1
elif a == b: # match
return 1
def make_score_matrix(seq1, seq2):
"""return score matrix and map(each score from which direction)0: diagnosis1: up2: left"""
seq1 = '-' + seq1
seq2 = '-' + seq2
score_mat = {}
trace_mat = {}
for i,p in enumerate(seq1):
score_mat[i] = {}
trace_mat[i] = {}
for j,q in enumerate(seq2):
if i == 0: # first row, gap in seq1
score_mat[i][j] = -j
trace_mat[i][j] = 1
continue
if j == 0: # first column, gap in seq2
score_mat[i][j] = -i
trace_mat[i][j] = 2
continue
ul = score_mat[i-1][j-1] + theta(p, q) # from up-left, mark 0
l = score_mat[i][j-1] + theta('-', q) # from left, mark 1, gap in seq1
u = score_mat[i-1][j] + theta(p, '-') # from up, mark 2, gap in seq2
picked = max([ul,l,u])
score_mat[i][j] = picked
trace_mat[i][j] = [ul, l, u].index(picked) # record which direction
return score_mat, trace_mat
def traceback(seq1, seq2, trace_mat):
'''find one optimal traceback path from trace matrix, return path code-!- CAUTIOUS: if multiple equally possible path exits, only return one of them -!-'''
seq1, seq2 = '-' + seq1, '-' + seq2
i, j = len(seq1) - 1, len(seq2) - 1
path_code = ''
while i > 0 or j > 0:
direction = trace_mat[i][j]
if direction == 0: # from up-left direction
i = i-1
j = j-1
path_code = '0' + path_code
elif direction == 1: # from left
j = j-1
path_code = '1' + path_code
elif direction == 2: # from up
i = i-1
path_code = '2' + path_code
return path_code
def print_m(seq1, seq2, m):
"""print score matrix or trace matrix"""
seq1 = '-' + seq1; seq2 = '-' + seq2
print()
print(' '.join(['%3s' % i for i in ' '+seq2]))
for i, p in enumerate(seq1):
line = [p] + [m[i][j] for j in range(len(seq2))]
print(' '.join(['%3s' % i for i in line]))
print()
return
def pretty_print_align(seq1, seq2, path_code):
'''return pair alignment result string frompath code: 0 for match, 1 for gap in seq1, 2 for gap in seq2'''
align1 = ''
middle = ''
align2 = ''
for p in path_code:
if p == '0':
align1 = align1 + seq1[0]
align2 = align2 + seq2[0]
if seq1[0] == seq2[0]:
middle = middle + '|'
else:
middle = middle + ' '
seq1 = seq1[1:]
seq2 = seq2[1:]
elif p == '1':
align1 = align1 + '-'
align2 = align2 + seq2[0]
middle = middle + ' '
seq2 = seq2[1:]
elif p == '2':
align1 = align1 + seq1[0]
align2 = align2 + '-'
middle = middle + ' '
seq1 = seq1[1:]
print('Alignment:\n\n' + align1 + '\n' + middle + '\n' + align2 + '\n')
return
def usage():
print('Usage:\n\tpython nwAligner.py seq1 seq2\n')
return
def main():
try:
seq1, seq2 = map(str.upper, sys.argv[1:3])
except:
seq1, seq2 = 'TCATC','TCATGGC'
usage()
print('--------Demo:-------\n')
print('1:%s' % seq1)
print('2:%s' % seq2)
score_mat, trace_mat = make_score_matrix(seq1, seq2)
#print_m(seq1, seq2, score_mat)
#print_m(seq1, seq2, trace_mat)
path_code = traceback(seq1, seq2, trace_mat)
pretty_print_align(seq1, seq2, path_code)
#print(' '+path_code)
if __name__ == '__main__':
main()
运行效果:
如果不加任何参数,运行帮助信息:
常规使用,直接接上两条序列:
代码开源在:https://github.com/Xiongyanshi/needleman_wunsch_pairAlignergithub.com
希望能对大家有所帮助,如有不当之处,还望大家留言评论或私信我。
———————————————— 20191120 终 ——————————————————
python 比对匹配_用Python从头实现Needleman-Wunsch序列比对算法相关推荐
- python标准词匹配_用 Python 自动化办公能做到哪些有趣或有用的事情?
我想介绍一下我是如何从每天工作8小时,进化成每天工作10分钟的. 不涉及太多的技术细节,毕竟知乎是一个分(现)享(编)知(故)识(事)的地方 0.先自我介绍一下: 我不是程序员,大学学的也不是IT专业 ...
- 原 python实现模糊匹配_使用python中的fuzzywuzzy库进行模糊匹配实例
fuzzywuzzy库是Python中的模糊匹配库,它依据 Levenshtein Distance 算法 计算两个序列之间的差异. Levenshtein Distance 算法,又叫 Edit D ...
- python正则表达式模糊匹配_用python正则表达式编译模糊正则表达式
当我发现python regex模块可以进行模糊匹配时,我感到非常高兴,因为它似乎是解决我许多问题的简单方法. 但是现在我遇到了一个问题,我没有从文档中找到任何答案. 如何使用新的模糊性值功能将字符串 ...
- python 地址模糊匹配_使用python处理selenium中的xpath定位元素的模糊匹配问题
# 用contains,寻找页面中style属性值包含有sp.gif这个关键字的所有div元素,其中@后面可以跟该元素任意的属性名. self.driver.find_element_by_xpath ...
- python最长匹配_二分图最大匹配:匈牙利算法的python实现
二分图匹配是很常见的算法问题,一般用匈牙利算法解决二分图最大匹配问题,但是目前网上绝大多数都是C/C++实现版本,没有python版本,于是就用python实现了一下深度优先的匈牙利算法,本文使用的是 ...
- python字符串去掉空行_从python中的字符串中删除空格
python字符串去掉空行 如何在python中删除字符串中的空格 (How to remove whitespaces in a string in python) str.lstrip()str. ...
- python delimiter分隔符用法_使用Python文件读写,自定义分隔符(custom delimiter)
众所周知,python文件读取文件的时候所支持的newlines(即换行符),是指定的.这一点不管是从python的doucuments上还是在python的源码中(作者是参考了python的io版本 ...
- python大牛 关东升_《Python从小白到大牛》第4章 Python语法基础
本章主要为大家介绍Python的一些语法,其中包括标识符.关键字.常量.变量.表达式.语句.注释.模块和包等内容. 标识符和关键字 任何一种计算机语言都离不开标识符和关键字,因此下面将详细介绍Pyth ...
- python之禅 中文_《Python之禅》中对于Python编程过程中的一些建议
<Python之禅>中对于Python编程过程中的一些建议 来源:中文源码网 浏览: 次 日期:2018年9月2日 [下载文档: <Python之禅>中对于Pyt ...
- python列表添加元组_【Python爬虫】列表、元组、集合练习
列表: pop() 函数用于移除列表中的一个元素(默认最后一个元素),并且返回该元素的值. list.append(obj) 在列表末尾添加新的对象 list.count(obj) 统计某个元素在列表 ...
最新文章
- \r与\n有何差别,编码的时候应该怎样使用
- 使用微信开发者工具创建小程序项目
- SAP RETAIL 通过自动补货功能触发的采购申请有些啥特殊的地方?
- Zonbu-售价 99 美元的袖珍电脑
- 西工大matlab计算机实验题,西工大信号系统上机实验一实验二
- Eclipse使用SVN操作说明
- 实现树状结构_钢结构设计 | “生命之树”景观案例赏析
- MSYS2+MinGW32 编译 QEMU需做的准备工作
- linux命令音乐视频合并,Linux下基于命令行的音乐播放器 (1)
- jQuery ui widget和jQuery plugin的实现原理简单比较
- 总结陈丹琦博士论文(一):NEURAL READING COMPREHENSION AND BEYOND
- ESD笔记(四)_击穿电压规律
- CC8编译报错:error #10099-D 解决方案(已解决)
- java hl7v3_HL7标准V3开发框架中个模型的关系
- 高通android q 通过backtrace,使用addr2ine工具,定位crash问题记录
- 大数据专业学校课程安排 (仅供参考)
- 【Cmake实战:番外】库、动态库和静态库(.dll,.so,.lib,.a)
- 从银行角度看二代征信
- 支持向量机原理之线性SVM与非线性SVM
- execute 等返回值
热门文章
- 【修真院web小课堂】如何理解html结构的语义化
- word表格内文字行间距调整方法
- Node-RED使用指南:7:配置与设定总结:其他配置
- 主板螺丝是机箱配还是主板配_MATX主板配什么机箱好?曜越Tt启航者A3装机记
- mysql 统计当个用户从当前时间连续登录天数,以及多用户某时间段,最长连续登录天数查询
- 安卓修改Airpods的双击功能,改“播放暂停”————下一首
- 秀!搭建一个永久运行的个人服务器!
- 发布订阅模式(一):tiny-emitter
- 加密的pdf文件如何解密?
- Spark从SQL的解析、执行与调优到Sparksql的解析的史上最全介绍