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

Posted

技术标签:

【中文标题】Needleman-Wunsch 实现在 cogent 和 skbio 中给出了不同的对齐方式【英文标题】:Needleman-Wunsch implementation gives different alignments in cogent and in skbio 【发布时间】:2014-08-15 19:21:59 【问题描述】:

与您在 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

【问题讨论】:

【参考方案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 获得的对齐方式不同,但这是一个小的实现差异:两个结果对齐方式具有相同的分数(差异在于 -- 是否对齐到第一个ATATAT 重复中的第二个AT),并且这两个实现在如何打破这种关系方面做出了不同的选择。

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

【讨论】:

以上是关于Needleman-Wunsch 实现在 cogent 和 skbio 中给出了不同的对齐方式的主要内容,如果未能解决你的问题,请参考以下文章

blast及其格式输出简介

200高分问下操作系统专家几个问题

Dapper 啥时候“实现”? Like Entity 在 ToList() 实现

@Service注解是标注在实现类上的的接口中添加注解还是在实现类impl

是否可以在接口中实现本机方法?

尽管没有在表格视图中实现尾随滑动动作,但删除按钮会在滑动时自动实现