Needleman-Wunsch实现的对齐在cogent和skbio中不同

Question

在skbio中的实现结果和在pycogent中的实现结果比较起来，感觉有点奇怪。

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的结果，可以在调用global_pairwise_align_nucleotide时传入penalize_terminal_gaps=True，像这样：

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的结果不完全相同，但这只是一个小的实现差异：两个对齐结果的得分是一样的（差异在于--是对齐到第一个AT还是第二个AT，在ATAT重复中），而两个实现选择打破这个平局的方式不同。

如果你回到你发布的对齐结果（scikit-bio的默认结果），你会发现返回的对齐结果非常好——实际上，如果不惩罚末尾的空缺，它就是得分最优的对齐（根据定义，因为它返回的就是得分最优的对齐）。不过，你说得对，这个结果看起来有点奇怪，因为在这个特定情况下，scikit-bio返回的对齐结果可能不是生物学上最相关的对齐。如果你知道你的序列都是从同一个位置开始并在同一个位置结束，你可以传入penalize_terminal_gaps=True。但是，你的例子是个简单的示例，而在大多数真实序列的情况下，我认为使用默认参数返回的对齐结果会更符合生物学的相关性。