python - 如何从 DNA 序列中获得比对分数？

Question

我对 Biopython 的 pairwise2 函数有点熟悉，但我注意到它在序列中添加了破折号以获得最佳的对齐分数。例如，

for a in pairwise2.align.globalxx("ACCGT", "ACG"):
  print(format_alignment(*a))

会产生这个结果：

ACCGT
|||||
A-CG-
Score=3
<BLANKLINE>
ACCGT
|||||
AC-G-
Score=3
<BLANKLINE>

即使第二个序列中的前 2 个字符（A 和 C）与第一个序列对齐。有没有办法找到对齐碱基对的数量而不是对齐碱基对的最高数量（例如：ACTGAA 序列相对于 GCCGTA 序列的得分为 3）？

score 0 · Accepted Answer

所以你只想计算两个序列（相同长度）中的相同碱基而不做任何比对？

像这样：

seq1 = 'ACTGAA'
seq2 = 'GCCGTA'

score = 0

for a, b in zip(seq1, seq2):
    if a == b:
        score +=1

print(score)

以更 Pythonic 的方式：

seq1 = 'ACTGAA'
seq2 = 'GCCGTA'

score = sum([1 for a, b in zip(seq1, seq2) if a == b])
print(score)

请注意，这个分数的倒数（不同碱基的数量）将是汉明距离。虽然您可以通过强制非常高的差距罚分来强制 Biopythonspairwise2返回您想要的结果，但上面显示的解决方案似乎更简单。

# I don't recommend this
pairwise2.align.globalxs(seq1, seq2, -1000, -1000)

score 0 · Accepted Answer

如果您只是试图阻止函数添加间隙，则可以增加间隙惩罚。对齐采用参数来设置匹配分数、不匹配罚分、创建空位罚分和扩展空位罚分：

pairwise2.align.globalxx("ACCGT", "ACG", 2, -1, -1, -0.5)

2 回答 2