是否有函数可以根据对齐参数计算对齐序列的分数？

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 是一个包含276个项目的字典。

我更倾向于补充缺少的项目，因为它只代表276个元素的一个迭代，而我们要分析的序列可能会有超过276个元素。因此，如果你使用函数 score_match() 来计算每对的分数，这个函数就得对序列中的每个元素进行测试 if pair not in matrix，这意味着肯定会比276次多得多。

还有一个耗时的地方是：每次 score += something 都会创建一个新的整数，并把名字 score 绑定到这个新对象上。每次绑定都需要一些时间，而如果用生成器直接生成整数并立即加到当前值上，就不会有这个时间消耗。

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。

编辑：
更新了适用于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)))

Answer 2

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分。虽然可能还有更好看的方法来实现这些，但这个应该是个不错的起点。