使用LevenshteinDistan获取子序列的位置

Question

我有大量包含序列的记录（'atcgttgcatcagttcga…'），最多500个字符。我也有一个小序列的列表，通常是10-20个字符。我想使用Levenshtein距离，以便在允许小变化或索引的记录中找到这些较小的序列（L_distance<；=2）。在

问题是我也想得到这些更小序列的起始位置，显然它只比较相同长度的序列。在

>>> import Levenshtein
>>> s1 = raw_input('first word: ')
first word: ATCGTAATACGATCGTACGACATCGCGGCCCTAGC
>>> s2 = raw_input('second word: ')
first word: TACGAT
>>> Levenshtein.distance(s1,s2)
29

在这个例子中，我想获得位置（7）和距离（在本例中是0）。在

有没有一个简单的方法来解决这个问题，或者我必须把较大的序列分解成更小的序列，然后对所有序列进行Levenshtein距离计算？那可能要花太多时间。在

谢谢。在

更新在寻找完全匹配后生成所有子字符串的朴素实现。

^{pr2}$

Answer 1

假设允许的最大Levenshtein距离很小，这可以在一次传递中完成，同时保留模糊匹配的候选列表。在

下面是我刚刚设计的一个实现示例。它没有经过彻底的测试、记录或优化。但至少它适用于简单的例子（见下文）。我试图避免由于跳过子序列边缘的字符而使它返回多个匹配项，但是正如我所说的，我还没有彻底测试过这一点。在

如果您感兴趣的话，我很乐意清理这个问题，编写一些测试，进行基本优化，并将其作为一个开源库提供。在

from collections import namedtuple

Candidate = namedtuple('Candidate', ['start', 'subseq_index', 'dist'])
Match = namedtuple('Match', ['start', 'end', 'dist'])

def find_near_matches(subsequence, sequence, max_l_dist=0):
    prev_char = None
    candidates = []
    for index, char in enumerate(sequence):
        for skip in range(min(max_l_dist+1, len(subsequence))):
            candidates.append(Candidate(index, skip, skip))
            if subsequence[skip] == prev_char:
                break
        new_candidates = []
        for cand in candidates:
            if char == subsequence[cand.subseq_index]:
                if cand.subseq_index + 1 == len(subsequence):
                    yield Match(cand.start, index + 1, cand.dist)
                else:
                    new_candidates.append(cand._replace(
                        subseq_index=cand.subseq_index + 1,
                    ))
            else:
                if cand.dist == max_l_dist or cand.subseq_index == 0:
                    continue
                # add a candidate skipping a sequence char
                new_candidates.append(cand._replace(dist=cand.dist + 1))
                # try skipping subsequence chars
                for n_skipped in range(1, max_l_dist - cand.dist + 1):
                    if cand.subseq_index + n_skipped == len(subsequence):
                        yield Match(cand.start, index + 1, cand.dist + n_skipped)
                        break
                    elif subsequence[cand.subseq_index + n_skipped] == char:
                        # add a candidate skipping n_skipped subsequence chars
                        new_candidates.append(cand._replace(
                            dist=cand.dist + n_skipped,
                            subseq_index=cand.subseq_index + n_skipped,
                        ))
                        break
        candidates = new_candidates
        prev_char = char

现在：

^{pr2}$

编辑：

在这个问题之后，我将编写一个Python库来搜索几乎匹配的子序列：^{}。这仍然是一项正在进行的工作。在

现在，试试^{}函数！它在您的用例中应该表现得特别好。在