是否有一个函数可以在给定对齐参数的情况下计算对齐序列的分数？

Question

我尝试对已经对齐的序列进行评分.让我们说吧

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE'
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE'

给定参数

substitution matrix : blosum62
gap open penalty : -5
gap extension penalty : -1

我确实浏览了biopython cookbook,但我能得到的是替换矩阵blogsum62,但我觉得必须有人已经实现了这种类型的库.

那么有人可以建议任何可以解决我的问题的库或最短的代码吗？

Thx提前

Answer 1

Jessada,

Blosum62矩阵(注意拼写;)在Bio.SubsMat.MatrixInfo中,是一个字典,元组解析得分(所以('A', 'A')值得4分).它没有间隙,它只是矩阵的一个三角形(所以它可能是('T','A')而不是('A','T').在Biopython中有一些辅助函数,包括Bio.Pairwise中的一些,但这是我想出的答案:

from Bio.SubsMat import MatrixInfo

def score_match(pair, matrix):
    if pair not in matrix:
        return matrix[(tuple(reversed(pair)))]
    else:
        return matrix[pair]

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e):
    score = 0
    gap = False
    for i in range(len(seq1)):
        pair = (seq1[i], seq2[i])
        if not gap:
            if '-' in pair:
                gap = True
                score += gap_s
            else:
                score += score_match(pair, matrix)
        else:
            if '-' not in pair:
                gap = False
                score += score_match(pair, matrix)
            else:
                score += gap_e
    return score

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE'
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE'

blosum = MatrixInfo.blosum62

score_pairwise(seq1, seq2, blosum, -5, -1)

对齐返回82.这几乎是一种更加漂亮的方式来完成所有这些,但这应该是一个好的开始.

Answer 2

blosum62是一个276项的dictonary.

我更喜欢用缺少的项目来完成,因为它代表只有276转的迭代,而要分析的序列可能有超过276个元素.因此,如果您在函数score_match()的帮助下找到每对的得分,则此函数将必须对if pair not in matrix序列的每个元素执行测试 ,也就是说肯定远超过276次.

另一件需要花费大量时间的事情:每个都score += something创建一个新的整数并将名称分数绑定到这个新对象.每个绑定花费一定时间,该时间不是由生成器的整数流不存在的,其立即被添加到当前量.

from Bio.SubsMat.MatrixInfo import blosum62 as blosum
from itertools import izip

blosum.update(((b,a),val) for (a,b),val in blosum.items())

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e, gap = True):
    for A,B in izip(seq1, seq2):
        diag = ('-'==A) or ('-'==B)
        yield (gap_e if gap else gap_s) if diag else matrix[(A,B)]
        gap = diag

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE'
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE'

print sum(score_pairwise(seq1, seq2, blosum, -5, -1))

这个score_pairwise()是一个生成函数,因为有yield而不是return.