Реализация 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, вы можете передать 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, но это незначительная разница в реализации: два итоговых выравнивания имеют одинаковую оценку (разница заключается в том, -- выравнивается с первым AT или второй AT в ATAT повторяю), и две реализации по-разному выбирают способ разрыва этой связи.

Если вы вернетесь к выровненному выравниванию (по умолчанию от scikit-bio), то вы заметите, что возвращаемое выравнивание очень хорошее - на самом деле, это оптимальное выигрышное выравнивание, если не штрафовать зазоры клемм (по определению, потому что оптимальное выравнивание выигрыша - это то, что он возвращает). Тем не менее, вы правы, что это странно, потому что выравнивание, которое возвращает scikit-bio в этом конкретном случае, вероятно, не является наиболее биологически значимым выравниванием. Если вы знаете, что все ваши последовательности начинаются в одной и той же позиции и заканчиваются в одной и той же позиции, вы можете передать penalize_terminal_gaps=True, Тем не менее, ваш пример - игрушка, и в большинстве случаев с реальными последовательностями, я думаю, наиболее биологически релевантное выравнивание будет возвращено с параметрами по умолчанию.