有没有一个函数可以计算给定对齐参数的对齐序列的得分？

Question

我试着给已经对齐的序列打分，比方说

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

带给定参数

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

我看了biopython的食谱，但我能得到的是替换矩阵blogsum62，但我觉得一定有人已经实现了这种库。

所以，谁能建议任何库或最短的代码，可以解决我的问题？

谢谢

Answer 1

杰萨达。

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项的独断。

我更倾向于用缺少的项目来完成，因为它代表的是只有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()是一个生成器函数，因为有 屈服 而不是 返回.

编辑：更新了Python 3的代码。

from Bio.SubsMat.MatrixInfo import blosum62 as blosum

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

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e, gap = True):
    for A,B in zip(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)))