python - Needleman-Wunsch 实现在 cogent 和 skbio 中给出了不同的对齐方式

Question

与在 pycogent 中实现的结果相比，skbio 中的实现给出了一个奇怪的结果。

from cogent.align.algorithm import nw_align as nw_align_cogent
from skbio.alignment import global_pairwise_align_nucleotide as nw_align_scikit

seq_1 = 'ATCGATCGATCG'
seq_2 = 'ATCGATATCGATCG'

print "Sequences: "
print "     %s" % seq_1
print "     %s" % seq_2
print

alignment = nw_align_scikit(seq_1, seq_2)
al_1, al_2 = [alignment.get_seq(_id).__str__() for _id in alignment.ids()]

print "    nw alignment using scikit:"
print "        %s" % al_1
print "        %s" % al_2
print

al_1, al_2 = nw_align_cogent(seq_1, seq_2)

print "    nw alignment using cogent:"
print "        %s" % al_1
print "        %s" % al_2
print

输出看起来像这样:

nw alignment using scikit:
    ------ATCGATCGATCG
    ATCGATATCGATCG----

nw alignment using cogent:
    ATCGAT--CGATCG
    ATCGATATCGATCG

Answer 1

这是一个很好的问题，并且与 scikit-bio 和 PyCogent 中比对评分​​的差异有关。默认情况下，在 scikit-bio 中，终端间隙不会受到惩罚，因为这会导致一些非常奇怪的对齐。这issue was briefly discussed here和插图here (见笔记本的最后一个单元格)。

如果您想实现更类似于 PyCogent 中的解决方案，您可以将 penalize_terminal_gaps=True 传递给 global_pairwise_align_nucleotide，如下所示:

alignment = nw_align_scikit(seq_1, seq_2, penalize_terminal_gaps=True)
al_1, al_2 = [alignment.get_seq(_id).__str__() for _id in alignment.ids()]

print "    nw alignment using scikit:"
print "        %s" % al_1
print "        %s" % al_2

输出:

nw alignment using scikit:
        ATCG--ATCGATCG
        ATCGATATCGATCG

您会注意到对齐方式仍然与您从 PyCogent 获得的对齐方式不同，但这是一个微小的实现差异:两个结果对齐方式具有相同的分数(不同之处在于 -- 与 ATAT 重复中的第一个 AT 或第二个 AT 对齐)，这两种实现方式在如何选择上有所不同他们打破了那种关系。

如果您返回到您发布的对齐方式(scikit-bio 的默认设置)，您会注意到返回的对齐方式非常好 - 事实上，如果不惩罚终端间隙(根据定义，因为最佳评分对齐是它返回的)。然而，你是对的，这很奇怪，因为 scikit-bio 在这种特定情况下返回的比对可能不是最具生物学相关性的比对。如果您知道您的序列都在相同位置开始并在相同位置结束，您可以传递 penalize_terminal_gaps=True。但是，您的示例只是一个玩具示例，在大多数情况下具有真实序列，我认为最具生物学相关性的比对将使用默认参数返回。