对齐匹配序列答案

【问题标题】：Aligning matching sequences对齐匹配序列
【发布时间】：2013-09-09 23:53:10
【问题描述】：

我正在尝试根据前 5 个和后 5 个字符对齐匹配序列。所以，例如：

AAATGCEGAIRPVOGJKKK
KKKTGKAFKEJWKLJFFFF
FFFKEJFWKFJWEKFJIGK

将对齐和连接创建：

AAATGCEGAIRPVOGJKKKTGKAFKEJWKLJFFFFKEJFWKFJWEKFJIGK

请注意，映射区域不应重复。我实际上有超过 3 行，并且它们没有排序，因此我试图创建一个循环以将所有它们对齐在一起。我不确定解决这个问题的最佳方法。

【问题讨论】：

你总是从第一行开始吗？如果不是，你如何选择从哪个序列开始？
如果匹配的可能性不止一种，应该怎么办？如果没有怎么办？
另外，你如何处理不匹配其他任何内容的序列。或者，如果您有一个与其他 10 个序列匹配的序列怎么办？有没有办法确定您的顺序偏好？
en.wikipedia.org/wiki/De_Bruijn_graph
我认为re 如果您正在使用蛮力，那将是矫枉过正。使用startswith 和endswith 字符串方法会更有效。

标签： python bioinformatics

【解决方案1】：

我认为你可以通过制作一个前缀字典来解决这个问题。一旦你有了它，从一个任意序列开始，发现它的后缀在前缀列表中。然后转到具有该前缀的序列，构建一个链。

这里有一些代码：

def make_circular_overlapping_sequence(sequences, min_overlap=1, max_overlap=5):
    # start by mapping prefixes to full sequences
    prefixes = {}
    for seq in sequences:
        for length in range(min_overlap, max_overlap+1):
            prefixes[seq[:length]] = seq

    # pick arbitrary a start sequence
    start = current = sequences[0]

    # build a chain of sequences with overlapping suffixes and prefixes
    chain = [start]
    while True:
        # try longest suffixes first
        for length in range(max_overlap, min_overlap-1, -1):
            suffix = current[-1-length:]
            if suffix in prefixes:
                current = prefixes[suffix]
                if current == start: # looped around, so we're done
                    return "".join(chain)
                chain.append(current[length+1:]) # don't duplicate the prefix
                break
        else: # for loop ended without breaking
            raise ValueError("No match found for sequence {!r}"
                             .format(current))

测试输出：

>>> sequences = '''AAATGCEGAIRPVOGJKKK
KKKTGKAFKEJWKLJFFFF
FFFKEJFWKFJWEKFJIGK
GKXYZ1234AAAT'''.split()
>>> make_circular_overlapping_sequence(sequences)
'AAATGCEGAIRPVOGJKKKTGKAFKEJWKLJFFFFKEJFWKFJWEKFJIGKXYZ1234AAAT'

注意事项：

此代码适用于明确的循环序列，如果存在不明确的匹配项（例如，"...ABC" 后面可以跟 "ABC..." 或 "BC..."），则可能无法正常工作。目前它总是选择最长的重叠。
如果链在没有循环的情况下结束，代码会引发异常，但如果存在不包含起始元素的短循环（例如["A...B", "B..C", "C..B"]），它可能会一直运行。
该算法不保证所有输入序列都包含在输出中。它只是找到一个循环并停止。
第一个和最后一个值的匹配部分不会被剪掉。如果你不想这样，我建议从返回值中切出最后的后缀：return "".join(chain)[:-length]

【讨论】：