Получить позицию подпоследовательности с использованием Levenshtein-Distance

Question

У меня есть огромное количество записей, содержащих последовательности ('ATCGTGTGCATCAGTTTCGA ...'), до 500 символов. У меня также есть список меньших последовательностей, обычно 10-20 символов. Я хотел бы использовать расстояние Левенштейна, чтобы найти эти меньшие последовательности внутри записей, допускающие небольшие изменения или индексы (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).

Есть ли простой способ решить эту проблему, или мне нужно разбить большие последовательности на более мелкие, а затем запустить расстояние Левенштейна для всех из них? Это может занять слишком много времени.

Спасибо.

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) 

Answer 1

Предполагая, что максимальное расстояние Левенштейна позволило мала, это может быть сделано в один проход, сохраняя при этом список кандидатов на нечетких совпадений.

Вот пример реализации, который я только что обработал. Он не проверен, не документирован или не оптимизирован. Но, по крайней мере, это работает для простых примеров (см. Ниже). Я пытался избежать возврата нескольких совпадений из-за пропуска символов на краях подпоследовательности, но, как я уже сказал, я не проверил это полностью.

Если вам интересно, я буду рад очистить это, написать несколько тестов, сделать базовую оптимизацию и сделать это доступным как библиотека с открытым исходным кодом.

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' 

EDIT:

После этого вопроса, я пишу библиотеку Python для поиска почти соответствие суб-последовательности: fuzzysearch. Это все еще очень много работает.

На данный момент, попробуйте функцию find_near_matches_with_ngrams()! Он должен особенно хорошо проявляться в вашем случае использования.