5

我有大量包含序列（'ATCGTGTGCATCAGTTTCGA...'）的记录，最多 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 距离？这可能需要太多时间。

谢谢。

UPDATE #Naive 实现在查找完全匹配后生成所有子字符串。

def find_tag(pattern,text,errors):       
    m = len(pattern)
    i=0
    min_distance=errors+1
    while i<=len(text)-m:
        distance = Levenshtein.distance(text[i:i+m],pattern)
        print text[i:i+m],distance #to see all matches.
        if distance<=errors:
            if distance<min_distance:
                match=[i,distance]
                min_distance=distance
        i+=1
    return match

#Real example. In this case just looking for one pattern, but we have about 50.
import re, Levenshtein

text = 'GACTAGCACTGTAGGGATAACAATTTCACACAGGTGGACAATTACATTGAAAATCACAGATTGGTCACACACACATTGGACATACATAGAAACACACACACATACATTAGATACGAACATAGAAACACACATTAGACGCGTACATAGACACAAACACATTGACAGGCAGTTCAGATGATGACGCCCGACTGATACTCGCGTAGTCGTGGGAGGCAAGGCACACAGGGGATAGG' #Example of a record
pattern = 'TGCACTGTAGGGATAACAAT' #distance 1
errors = 2 #max errors allowed

match = re.search(pattern,text)

if match:
    print [match.start(),0] #First we look for exact match
else:
    find_tag(pattern,text,errors)
4

1 回答 1

8

假设允许的最大 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

现在：

>>> list(find_near_matches('bde', 'abcdefg', 0))
[]
>>> list(find_near_matches('bde', 'abcdefg', 1))
[Match(start=1, end=5, dist=1), Match(start=3, end=5, dist=1)]
>>> list(find_near_matches('cde', 'abcdefg', 0))
[Match(start=2, end=5, dist=0)]
>>> list(find_near_matches('cde', 'abcdefg', 1))
[Match(start=2, end=5, dist=0)]
>>> match = _[0]
>>> 'abcdefg'[match.start:match.end]
'cde'

编辑：

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

现在，试试这个find_near_matches_with_ngrams()功能！它应该在您的用例中表现得特别好。

于 2013-11-01T13:29:49.473 回答