NW全局比对算法原理及python实现 (考虑gap长度)

生信编程日常

发布于 2020-04-14 17:10:25

3.2K0

发布于 2020-04-14 17:10:25

文章被收录于专栏：生物信息学、python、R、linux生物信息学、python、R、linux

在序列比对的时候，有全局比对和局部比对两种方法，其中，Needleman-Wunsch比对算法是其中的一个很经典的全局比对算法。下面将用python从头实现，将考虑match，mismatch，gap和gap是否连续的因素。

先确定打分策略，先考虑match，mismatch和gap的分数。现将match定为1分，mismatch -1分，gap是-2分。

现有两个序列，一个是ATCG，一个是TCG，以下是初始分值矩阵：

first_score

第一列第一行表示一开始比对到gap上，因为计入gap的分数。

接着计算每行每列的分数，比如第三行第三列，可以得到三个分数，同一行列的为一个gap：

从三种路径中选择分数最高的作为这个位置的分数。接下来每个位置都要进行这样的计算，直到填满这个表格。

form

同时可以得到最佳路径：

path

以下用python实现：

比较两个碱基是否一致

# 比较两个碱基分数， gaps的分数考虑在cal_score
def diff(first, second):
    if first == second:
        return match 
    else:
        return mismatch

计算分数，并考虑多段gap的情况。一般来说，结果应该更倾向于仅出现一段长gap而不是多个短gap。所以将第一次出现的gap定位-2分，第二次及之后连续出现的gap定为一个更小的分数。即变量extend。

def cal_score(seq1, seq2):

    nrow = len(seq1)
    ncol = len(seq2)

    scores = np.zeros((nrow, ncol))
    path = np.zeros((nrow, ncol))
    # 初始化第一行第一列
    for index1 in range(nrow):
        scores[index1,0] = gap * index1
        path[index1,0] = 0 #上面

    for index2 in range(ncol):
        scores[0,index2] = gap * index2
        path[0, index2] = 2 #左侧


    for num1 in range(1, nrow):
        for num2 in range(1, ncol):

            # 得到上面，斜上和左侧的结果
            last_score = [scores[num1 - 1, num2], scores[num1 - 1, num2 - 1], scores[num1, num2 - 1]]
            change_score = diff(seq1[num1], seq2[num2])
            current_score = []

            if path[num1 - 1 , num2] == 0:
                current_score.append(scores[num1-1, num2] + extend) # 上面
            else:
                current_score.append(scores[num1-1, num2] + gap)

            current_score.append(scores[num1-1, num2 - 1] + change_score) # 斜上

            if path[num1 , num2 -1] == 2:
                current_score.append(scores[num1, num2-1] + extend) # 左侧
            else:
                current_score.append(scores[num1, num2-1] + gap)    

            current_index = current_score.index(max(current_score)) # 当前索引，不是0就gap

            scores[num1, num2] = max(current_score)

            path[num1, num2] = current_index


    return scores, path

回溯路径并输出结果

def cal_seq(scores, path):
    index1 = len(seq1) - 1
    index2 = len(seq2) -1

    top = ''
    middle = ''
    bottom = ''

    while True:
        if path[index1, index2] == 1:
            top += seq1[index1]
            bottom += seq2[index2]

            if seq1[index1] == seq2[index2]:
                middle += '|'
            else:
                middle += ' '

            index1 -= 1
            index2 -= 1

        elif path[index1, index2] == 0:
            top += seq1[index1]
            bottom += '-'
            middle += ' '

            index1 -= 1

        else:
            top += '-'
            bottom += seq2[index2]
            middle += ' '

            index2 -= 1 

        
        top_num = len(top) - top.count('-')
        bottom_num = len(bottom) - bottom.count('-')
        if top_num == max(len(seq1), len(seq2))-1 or bottom_num == max(len(seq1), len(seq2))-1:
            break

    return top, middle, bottom

4.命令行解析与主函数

import argparse
import numpy as np

def getArgs(argv=None):
    parser = argparse.ArgumentParser()
    parser.add_argument("-seq1", help = "", required = True)
    parser.add_argument("-seq2", help = "", required = True)
    args = parser.parse_args()
    return args

if __name__ == '__main__':
    args = getArgs()
    seq1 = "_" + args.seq1
    seq2 = "_" + args.seq2


    match=5  # 或者别的分数
    mismatch=-5
    gap=-10  
    extend = -0.5  

    scores,path = cal_score(seq1, seq2)
    top, middle, bottom = cal_seq(scores, path)

    print(top[::-1])
    print(middle[::-1])
    print(bottom[::-1])

调用：python extend_alignment.py -seq1 ATCGATGGTATATATCGATC -seq2 ATCGATGAGTATAT

输出：