基因序列比较

基因序列比较

设计算法，计算两给定基因序列的相似程度。

人类基因由4种核苷酸，分别用字母ACTG表示。要求编写一个程序，按以下规则比较两个基因序列并确定它们的相似程度。即给出两个基因序列AGTGATG和GTTAG，它们有多相似呢？测量两个基因相似度的一种方法称为对齐。使用对齐方法可以在基因的适当位置加入空格，让两个基因的长度相等，然后根据基因的分值矩阵计算分数。

基因分数表：

Score	A	C	G	T	-
A	5	-1	-2	-1	-3
C	-1	5	-3	-2	-4
G	-2	-3	5	-2	-2
T	-1	-2	-2	5	-1
-	-3	-4	-2	-1	*

 

例：比较AGTGATG与GTTAG

第一种对齐方案为：

首先可以给AGTGATG插入一个空格得：AGTGAT-G

         GTTAG插入3个空格即得：  -GT--TAG

上面的匹配分值为:-3+5+5+(-2)+(-3)+5+(-3)+5=9.

第二种对齐方案为：

AGTGATG

-GTTA-G

得到的分值为：(-3)+5+5+(-2)+ 5+(-1)+5=14.

当然还有其它对齐方式，但以上对齐方式是最优的，所以两个基因的相似度就为14。

问题分析与解决思路

设两条序列分别为X序列和Y序列，长度分别为m和n。

定义一个表结构来抽象该问题。比较两个序列的相似度，可以理解为寻找一条由坐标（0,0）到坐标（m-1，n-1），并且只能向下、向右或者向右下延伸的路径，使得这条路径得分最高。

如下图是例子中，第二中对齐方案（图中“↓”意味着X序列的该核苷酸与Y序列的“-”匹配，“→”则相反）：

抽象后的序列相似度问题：

Y	G	T	T	A	G
A	↓
G	↘
T		↘
G			↘
A				↘
T					↓
G					↘

 

如果用暴力搜索方法求解该问题，就要穷举X的所有子序列和Y的所有子序列，并对它们进行匹配，记录所有匹配情况下，两条序列的得分。

根据分析，可以得出基因序列比较问题具有最优子结构性质。规模为m，n的问题的解，可以由规模为m-1，n-1的子问题，规模为m-1，n的子问题，和规模为m，n-1的子问题中，分别加上规模为规模为m，n时最后一个字符的得分结果中，最大的那一个。用公式表示为：

s1 = result[i-1][j-1] + get_score(X[i], Y[j])

s2 = result[i-1][j] + get_score(X[i], '-')

s3 = result[i][j-1] + get_score('-', Y[j])

result[i][j] = max(s1,s2,s3)

模型建立与算法描述

记两条序列分别为X序列和Y序列，长度分别为m和n。

我们将上述分析过程和解决的思路进一步归纳为以下步骤：

（1）初始化结果数组result的第一行及第一列。由于当X的规模m等于1时，意味着出了一个匹配的位置之外，Y序列的其他元素需要和空格匹配，所以该情况下，子问题由result[i][j] = max(s1,s2,s3)改变为result[i][0] = result[i-1][0] + get_score(X[i], '-')。当Y的规模n等于1时同理。

（2）从i等于2和j等于2开始，根据result[i][j] = max(s1,s2,s3)，由底向上地计算结果，同时也可以利用另一个表用于记录路径。

（3）返回结果表，其中result[m-1][n-1]为该问题的结果。

复杂度分析

本算法在M和N前面插入一个空字符串时，新建了两个数组来存放新的序列，所以copy操作花费了m+n。

初始化时有m+n次赋值操作，加上运算时(m×n)次赋值操作，一共(m×n)+m+n次赋值，所以时间复杂度为θ(m×n)。

运算时创建了一个额外的二维数组——结果表，所以空间复杂度为(m×n)+(m+n)×2，所以空间复杂度为θ(m×n)。

运行结果分析

以下为比较AGTGATG与GTTAG的运算结果：

根据图4.2，可以发现当i或j等于1时，无法继续向前回溯了，需要另外处理。图4.3标注了从path[m-1][n-1]开始回溯时，实际上的起点和终点。

得分表

路径记录情况

 

 匹配结果及得分标题

代码： 

from copy import deepcopy
import numpy as np
def get_score(x, y):
    index = {
        'A':0,
        'C':1,
        'G':2,
        'T':3,
        '-':4,
    }
    score_table = [
        [5, -1, -2, -1, -3],
        [-1, 5, -3, -2, -4],
        [-2, -3, 5, -2, -2],
        [-1, -2, -2, 5, -1],
        [-3, -4, -2, -1, float('-inf')]
    ]
    return score_table[index[x]][index[y]]
def ACTG_sim(M, N):
    X = deepcopy(M)
    Y = deepcopy(N)
    X.insert(0,'')
    Y.insert(0,'')
    m = len(X)
    n = len(Y)
    path = [['' for x in range(n)] for y in range(m)]
    path = np.array(path, dtype=str)
    result = np.zeros((m,n)) #存值
    #注意初始化
    for i in range(1,m):
        result[i][0] = result[i-1][0] + get_score(X[i], '-')
    for j in range(1,n):
        result[0][j] = result[0][j-1] + get_score('-', Y[j])
    for i in range(1, m):
        for j in range(1, n):
            s1 = result[i-1][j-1] + get_score(X[i], Y[j])
            s2 = result[i-1][j] + get_score(X[i], '-')
            s3 = result[i][j-1] + get_score('-', Y[j])
            max_score = max(s1,s2,s3)
            result[i][j] = max_score
            #第一种情况
            if s1 == max_score:
                path[i][j] = '`'
            #第二种情况
            if s2 == max_score:
                path[i][j] = '|'
            #第三种情况
            if s3 == max_score:
                path[i][j] = '-'
    return path, result
#回溯法输出结果
def print_sim(path, X, Y, i, j):
    if i == 0 or j == 0:
        return
    if path[i, j] == '`':
        if i == 1 and j != 1:
            path[i, j - 1] = '-'
            print_sim(path, X, Y, i, j - 1)
        elif j == 1 and i != 1:
            path[i - 1, j] = '|'
            print_sim(path, X, Y, i - 1, j)
        else:
            print_sim(path, X, Y, i - 1, j - 1)
        print(X[i-1], Y[j-1], '\t', get_score(X[i-1], Y[j-1]))
    if path[i, j] == '|':
        print_sim(path, X, Y, i - 1, j)
        print(X[i-1], '-', '\t', get_score(X[i-1], '-'))
    if path[i, j] == '-':
        print_sim(path, X, Y, i, j - 1)
        print('-', Y[j-1], '\t', get_score('-', Y[j-1]))
if __name__ == '__main__':
    # M = 'AGTGATG'
    # N = 'GTTAG'
    M = input('序列M：')
    N = input('序列N：')
    import re
    if re.match('[^ACTG]', M) or re.match('[^ACTG]', N):
        print('input error')
        exit()
    M = list(map(str, M))
    N = list(map(str, N))
    path, memory = ACTG_sim(M, N)
    print(path)
    print(memory)
    m = len(M)
    n = len(N)
    print_sim(path, M, N, m, n)
    print('得分为：%d' % memory[m][n])