有没有一种函数可以根据对齐参数计算对齐序列的分数？

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，但我感觉已经有人实现了这种库。

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

提前感谢。

- Jessada Thutkawkorapin

---

为什么不直接用Biopython blastall来进行Blast呢？ - joaquin

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。肯定有更漂亮的方法来完成所有这些操作，但这应该是一个不错的开始。

- david w

Answer 2

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))

因为使用了yield而非return，所以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)))

- eyquem

---

虽然我没有测试过，但是这个答案的代码返回的结果相同，而且似乎比被接受的答案快得多。 - thposs