Needleman-Wunsch实现给出了在实际和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中，终端间隙不会受到惩罚，因为这会导致一些非常奇怪的排列。这个这里简要地讨论了这个问题。和插图的这里 (见笔记本的最后一个单元格)。

如果希望实现与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。然而，你的是一个玩具的例子，在大多数情况下，与真实的序列，我认为最重要的生物相关对齐将返回与默认参数。