如何从 DNA 序列中获得比对分数?

How to get an alignment score from DNA sequences?

我对 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)?

如果您只是想防止函数添加间隙,您可以增加间隙惩罚。对齐采用参数来设置匹配分数、非匹配惩罚、创建空位惩罚和扩展空位惩罚:

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

所以您只想计算两个序列(相同长度)中相同的碱基而不进行任何比对?

像这样:

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)

请注意,这个分数的倒数(不同碱基的数量)就是汉明距离。虽然您可以通过强制非常高的间隙惩罚来强制 Biopythons pairwise2 达到 return 您想要的结果,但上面显示的解决方案似乎更简单。

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