基于动态规划进行双序列全局比对

说明

核酸序列打分算法脚本,基于动态规划进行双序列全局比对,得到两条DNA序列的相似度并打分,但程序还有一些问题,在匹配长序列的时候还不够完善.(20年初写的放在简书上,这次copy到博客园,欢迎访问我的简书:https://www.jianshu.com/p/6bedaa9f6de1)

环境

Linux、Python3.6

实例

 command

 result

以下为代码

# -*- coding: utf-8 -*-
"""
Created on Wed Feb 19 13:42:52 2020
@author: moDD
Email:312198221@qq.com
"""
import pandas as pd
import re,argparse
x_seq = 'ATCGATGTGTAGTATATATCGATCAGTTGA'
y_seq = 'ATCGATGTCTAAGTATAT'
def parameter():
    '''设置三个参数'''
    parser = argparse.ArgumentParser(prog=" Pairwise Sequence Alignment ", usage='python3.6 ./Pairwise_Sequence_Alignment.py -seq_a ATCGATGTGTAGTATATATCGATCAGTTGA -seq_b ATCGATGTCTAAGTATAT  -o ./output/result.txt',
                                     description="description:此脚本基于动态规划进行双序列全局比对,输入数据为a序列和b序列,输出为文本文件,得到两条DNA序列的相似度", epilog="Tip:此脚本使用python3.6编写完成, 请尽量使用python3.6版本执行")
    parser.add_argument("-seq_a", dest='seq_a',type=str, help="first sequence. for example:ATCGATGTGTAGTATATATCGATCAGTTGA")
    parser.add_argument("-seq_b", dest='seq_b',type=str, help="second  sequence. for example:ATCGATGTCTAAGTATAT")
    parser.add_argument("-o", dest='outfilename',type=str, help="The name of result. for example:result.txt")
    
    (para, args) = parser.parse_known_args()
    try:
        x_seq= para.seq_a
        y_seq= para.seq_b
        if  len(y_seq) > len(x_seq):   #确保x为长序列 y为短序列
            x_seq= para.seq_b
            y_seq= para.seq_a
        out_file_name  = para.outfilename
    except:
        print('Missing parameters or Parameter error! Please check parameters!')
        raise ValueError
    #没有设置这些参数的外部输入 ,如有需要直接添加即可
    gap = -5
    wrong = -5
    right = 2
    base_pair = -2
    return (x_seq,y_seq,out_file_name,right,base_pair,wrong,gap)

def Generating_scoring_matrix(right,base_pair,wrong,gap):
    '''创建分数数据框'''
    scoring_matrix =  pd.DataFrame({
                            '-':[0,gap,gap,gap,gap],
                            'A':[gap,right,base_pair,wrong,wrong],
                            'T':[gap,base_pair,right,wrong,wrong],
                            'C':[gap,wrong,wrong,right,base_pair],
                            'G':[gap,wrong,wrong,base_pair,right]
                            },
                            index = ['-','A','T','C','G']
                            )
    return scoring_matrix
    
def cutText(text, sec):
    '''切割字符串为多段  每段长度为sec'''
    return [text[i:i+sec] for i in range(0,len(text),sec)]

def Adjust_output_format_and_output(align_a,Middle_row,align_b,out_file_name):
    '''切割字符串为固定长度 并保存为指定文件名的文件'''
    with open(out_file_name , 'w') as f:
        index = 1
        for (row_1,row_2,row_3) in zip(cutText(align_a,50),cutText(Middle_row,50),cutText(align_b,50)):
            end_len_row_1 = len(re.sub('-',"",row_1)) #去除减号 得到长度 加在字符串末尾
            end_len_row_3 = len(re.sub('-',"",row_3)) #同上
            element = str('Query' + '\t' + str(index) + '\t' + row_1  + '\t' + str(end_len_row_1) +'\n'+
                          '     ' + '\t' + ' '        + '\t' +  row_2 + '\n'+
                          'sbjct' + '\t' + str(index) + '\t' +  row_3 + '\t' + str(end_len_row_3) +'\n\n')
            f.write(element)  #写入
            index += 1
            
def compute_result_matrix(x_seq, y_seq, scoring_matrix):
    '''得到一个高为length_x+1,宽为length_y+1 的数据框.  即(length_x+1) * (length_y+1) '''
    length_x = len(x_seq)
    length_y = len(y_seq)
    result_matrix = [[0 for i in range(length_y + 1)] for j in range(length_x + 1)]
    result_matrix = pd.DataFrame(result_matrix)
    
    #根据动态规划算法 , 首先,生成第0列的数据 依次为 0 -5 -10 -15 
    for x_index in range(1, length_x+1):
        result_matrix[0][x_index] = result_matrix[0][x_index-1] + scoring_matrix[x_seq[x_index-1]]['-']  #数据框 列index在前面 行index在后面
    #之后,生成第0行的数据 依次为 0 -5 -10 -15 
    for y_index in range(1, length_y+1):
        result_matrix[y_index][0] = result_matrix[y_index-1][0] + scoring_matrix[y_seq[y_index-1]]['-']
    #最后从数据框的左上角开始,向右下角逐步计算每一个值
    for x_index in range(1,length_x+1):
        for y_index in range(1, length_y+1):
            #取以下三者的最大值 这三个数分别为: 1,此位置左上角的值 + 得分矩阵中两个字符对应的值
            #                               2,此位置上面的值 + 得分矩阵中的gap
            #                               2,此位置左边的值 + 得分矩阵中的gap
            result_matrix.iloc[x_index,y_index]=max(result_matrix.iloc[x_index-1,y_index-1]+ scoring_matrix.loc[y_seq[y_index-1],x_seq[x_index-1]],
                                               result_matrix.iloc[x_index,y_index-1]  + scoring_matrix.loc['-',x_seq[x_index-1]],  #x序列对应y的空值
                                               result_matrix.iloc[x_index-1,y_index]  + scoring_matrix.loc[y_seq[y_index-1],'-']   #y序列对应x的空值
                                               )
    return (result_matrix)

def compute_global_alignment(x_seq, y_seq, scoring_matrix, result_matrix):
    '''将 矩阵数据逆推回序列数据'''
    
    #确定累积得分最大值是在数据框的最后一列还是最后一行,用于确定累积得分最大值所在的索引 如[17,18]
    length_x = len(x_seq)
    length_y = len(y_seq)
    terminal_max = max(max(result_matrix[length_y]),  #最后一列最大值
                   max(result_matrix.loc[length_x,:]) #最后一行最大值
                )
    if terminal_max in list(result_matrix[length_y]):
        the_value_x_index = list(result_matrix[length_y]==terminal_max).index(True)
        the_value_x_y_index = [the_value_x_index , length_y]
        x_index=the_value_x_y_index[0]
        y_index=the_value_x_y_index[1]
    else:
        the_value_y_index = list(result_matrix.loc[length_x,:]==terminal_max).index(True)
        the_value_x_y_index = [length_x , the_value_y_index]
        x_index=the_value_x_y_index[0]
        y_index=the_value_x_y_index[1]
        
    #取此位置以后的两端序列, 开始向前依次添加ATCG或者'-'
    section_x_seq = x_seq[x_index:]
    section_y_seq = y_seq[y_index:]
    #因为从右下角到左上角依次检索,所以先给字符串反转,然后再尾部添加. 这一过程相当与从尾部向头部依次添加
    section_x_seq = section_x_seq[::-1]
    section_y_seq = section_y_seq[::-1]
    #此过程为从后往前把序列补齐 
    while x_index and y_index:
        #如果是1,相同的碱基
        #     2,AT GC互补 , 
        #     3,AG TC错配          这三种情况之一则分别添加原序列的原位置的碱基
        if result_matrix.iloc[x_index,y_index] == result_matrix.iloc[x_index-1,y_index-1] + scoring_matrix[x_seq[x_index-1]][y_seq[y_index-1]]:
            section_x_seq += x_seq[x_index-1]#;print(1)
            section_y_seq += y_seq[y_index-1]
            x_index -= 1
            y_index -= 1
        #否则 , 分别添加原序列的原位置的碱基和'-'
        else:
            if result_matrix.iloc[x_index,y_index] == result_matrix.iloc[x_index-1,y_index] + scoring_matrix[x_seq[x_index-1]]['-']:
                section_x_seq += x_seq[x_index-1]#;print(1)
                section_y_seq += '-'
                x_index -= 1
            else:
                section_x_seq += '-'#;print(1)
                section_y_seq += y_seq[y_index-1]
                y_index -= 1
    #如果x_index或者y_index 其中一个未归零,另个为零, 则直接给未归零的序列补'-'
    while x_index:
        section_x_seq += x_seq[x_index-1]#;print(1)
        section_y_seq += '-'
        x_index -= 1
    while y_index:
        section_x_seq += '-'#;print(1)
        section_y_seq += y_seq[y_index-1]
        y_index -= 1
    #把倒转的序列再转回来
    result_x_seq = section_x_seq[::-1]
    result_y_seq = section_y_seq[::-1]
        
    # 使section_x_seq 和section_y_seq为同一长度 , 短序列补值'-'
    length_x = len(result_x_seq)
    length_y = len(result_y_seq)
    if length_x < length_y:
        result_x_seq += '-' * (length_y - length_x)#;print(1)
    else:
        result_y_seq += '-' * (length_x - length_y)#;print(1)
        
    #依据补值完成的两列数据和得分矩阵 , 计算总得分
    Total_score = sum([scoring_matrix[result_x_seq[x_index]][result_y_seq[x_index]] for x_index in range(len(result_x_seq))])
    ###################################################################################
    #得到输出结果的中间行 例如 '|||||||||| |||| ||||||'
    Middle_row=''
    for (x_element,y_element) in zip(result_x_seq,result_y_seq):
        if x_element==y_element:
            Middle_row += '|'
        else:
            Middle_row += ' '
    return Total_score, result_x_seq, result_y_seq,Middle_row

################################################################################
if __name__ == '__main__':
        (x_seq,y_seq,out_file_name,right,base_pair,wrong,gap) = parameter() #得到所有参数
        scoring_matrix = Generating_scoring_matrix(right,base_pair,wrong,gap) #生成得分矩阵
        result_matrix = compute_result_matrix(x_seq=x_seq, y_seq=y_seq, scoring_matrix=scoring_matrix) #生成序列得分矩阵
        score, result_x_seq, result_y_seq,Middle_row = compute_global_alignment(x_seq=x_seq, y_seq=y_seq, scoring_matrix=scoring_matrix, result_matrix=result_matrix) #将矩阵转化为结果序列
        Adjust_output_format_and_output(result_x_seq,Middle_row,result_y_seq,out_file_name) #整理数据,写入数据