这篇博文以简化的 Needleman-Wunsch 比对算法为例详解比对算法的原理及代码实现。 Needleman-Wunsch 算法是最著名的全局比对算法，在此基础上形成的 Smith-Waterman 算法是最著名的局部比对算法，虽然两者都非常巧妙，但基本原理和代码相差不大，这里仅讲解 Needleman-Wunsch 算法 。 1. Needleman-Wunsch 算法的原理 1970 年 Saul B. Needleman and Christian D. Wunsch 首次将动态规化的思路应用到生物信息学领域，形成了大名鼎鼎的 Needleman-Wunsch 算法，该算法在当前生物信息学领域得到广泛应用，是全局比对算法中最重要的算法。 闲言少述，直接进入正题。以下两条序列进行全局比对为例， Needleman-Wunsch 算法的具体步骤如下： sequence1 GCATGCU sequence2 GATTACA 1.1 初始化矩阵 首先建立一个空的矩阵，矩阵上的行名为 sequence1 的碱基，矩阵的列名为 sequence2 的碱基。因为需要初始值，所有数据区第一行和第一列依次为 -1, -2, -3, -4... 这组递减数列，相连两项之间的差值为 gap 罚分。 这里为了方便起见，采用最一般的罚分，即 match 得 1 分， mismatch 和 gap 罚 1 分（即得 -1 分，以下叙述均用得分表示，若为罚分则得分为负）。 图 1. 初始空矩阵 注： match=1, mismatch=-1, gap=-1 1.2 三个方向移动到当前位点时的综合得分 依次从左上往右下计算出每个位点的得分，计算时先算出从左，从上及从左上角移动到当前位点时的得分，这个得分值为： 不同方向移动综合得分 = 移动前位点的得分 + 移动过程的得分 移动前位点的得分为移动前位点方框中的值，移动过程的得分按 1.1 中的得分约定计算如下： Ø 从上往下和从左往右移动时都会引入 gap ，前者是在横向这条序列上引入 gap ，后者是在纵向这条序列上引入 gap ，因此都会得 -1 分； Ø 从左上往右下方向移动时，如果当前位点横向和纵向对应碱基一致，表明为 match ，得 1 分；如果当前位点横向和纵向对应碱基不一致，表明为 mismatch ，得 -1 分。 具体示意图如下： 图 2. 三个方向移动时综合得分 值得一提的是， Smith-Waterman 算法 仅在此基础上加入一个 0 值，让上述三个方向上的值与 0 ，共四个值比较大小。而且在最开始初始化矩阵时将初始行初始列的值都设为 0 。这样整个矩阵中的值没有负值。因此可以从任何位置开始，后面回溯时从矩阵中最大的值开始回溯，这样就可以达到局部比对的效果，真的是太精妙了！ 1.3 当前位点得分 从三个方向（从上到下，从左到右，从左上到右下）移动到当前位点的综合得分的最大值，即为当前位点得分。 那么由上图 2 中可得三个方向移动到当前位点时的最大得分为 max(1, -2, -2) = 1 即当前位点得分为 1 （图中表格第三行第三列值为 1 ） 按照以上原则，将矩阵上每个位点都计算出来，填满整个表格。 强列建议大家手动算一次 ，实际计算会让思路更加清晰，这对后面写代码实现非常有帮助！ 我手动计算结果如下，与 wiki 百科 上的一致。 图 3. 整个矩阵的结果图 1.4 回溯获取最佳比对结果 沿右下角向左上角回溯，每个位点依然有三个位置，左上，左边和上边，如果最大值出现在上面，则横向这条序列引入一个 gap (-) ，纵向这条序列取该处碱基；如果最大值出现在左边，则纵向这条序列引入一个 gap (-) ，横向这条序列取该处碱基 ; 如果最大值出现在左上角，则不引入 gap ，纵向和横向均取该处碱基。这样获取到两段序列，再反转过来（因为序列是从后往前回溯的）即为最终结果。 图 4. 回溯找最佳路径 2. 用 Python 实现全局比对 方便起见，这里仅用了原生python实现算法，具体代码如下，有兴趣的朋友也可以使用cython或python数据分析包，或者其它语言实现，并欢迎贴到讨论区，谢谢。 --------------------- 可下载代码源文件 global_alignment.py --------------------- #!/usr/bin/env python import sys __author__ = 'luria@sohu.com' __date__ = '2018.05.04' __version__ = 'v1.0' def main(self, subject, query): match, mismatch, gap = 1, -1, -1 # if you wanna to use other score matrix, # you could do code reactoring as a practice pos_dict = {(i, j): 0 for i in xrange(len(subject)+1) for j in xrange(len(query)+1)} for k in pos_dict: if not k and not k : pos_dict = 0 elif not k and k : pos_dict = k * gap elif not k and k : pos_dict = k * gap # print pos_dict # this step must be follow first loop for i in xrange(1, len(subject)+1): for j in xrange(1, len(query)+1): up2down = pos_dict + gap left2right = pos_dict + gap if subject == query : topleft2bottomright = pos_dict + match else: topleft2bottomright = pos_dict + mismatch pos_dict = max(up2down, left2right, topleft2bottomright) # print matrix ''' for i in xrange(len(subject)+1): temp = )) print \\t.join(temp) ''' out_subject, out_query = '', x, y = len(subject), len(query) while 1: if not x and not y: break direct_dict = { up : pos_dict , left : pos_dict , top_left : pos_dict } #print direct_dict order = sorted(direct_dict.iteritems(), key=lambda o:o , reverse=True) # only get one best path at this program, # you can get all best path if you like if order == up: out_subject += - out_query += query x -= 1 elif order == 'left': out_query += - out_subject += subject y -= 1 else: out_subject += subject out_query += query x -= 1 y -= 1 print out_subject print .join( == out_query else for i in xrange(len(out_subject))] ) print out_query print if __name__ == '__main__': if len(sys.argv) == 1: sys.exit( global_align.py subject query) main(*sys.argv) 参考材料： https://en.wikipedia.org/wiki/Smith-Waterman_algorithm https://en.wikipedia.org/wiki/Needleman-Wunsch_algorithm
标签: 序列比对
