Прежде чем продолжить, смотрели ли вы на биопайтон?
Похоже, что вы хотите найти приблизительные совпадения с одной ошибкой замены и нулевыми ошибками вставки / удаления, то есть с расстоянием Хэмминга, равным 1.
Если у вас есть функция сопоставления расстояния Хэмминга (см., Например, ссылку, предоставленную Игнасио), вы можете использовать ее следующим образом для поиска первого совпадения:
any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))
но это было бы довольно медленно, потому что (1) функция расстояния Хэмминга будет продолжать измельчать после второй ошибки подстановки (2) после неудачи, она перемещает курсор на единицу, а не пропускает вперед, основываясь на том, что она увидела (как функция Бойера- Поиск Мура).
Вы можете преодолеть (1) с помощью такой функции:
def Hamming_check_0_or_1(genome, posn, sequence):
errors = 0
for i in xrange(25):
if genome[posn+i] != sequence[i]:
errors += 1
if errors >= 2:
return errors
return errors
Примечание: это намеренно не Pythonic, а Cic, потому что вам нужно будет использовать C (возможно, через Cython), чтобы получить разумную скорость.
Некоторая работа над побитно-параллельным приближенным поиском Левенштейна с пропуском была проделана Наварро и Раффино (Google "Navarro Raffinot nrgrep"), и это может быть адаптировано для поиска Хэмминга. Обратите внимание, что у бит-параллельных методов есть ограничения по длине строки запроса и размеру алфавита, но у вас 25 и 4 соответственно, так что проблем нет. Обновление: пропуск, вероятно, не очень помогает с размером алфавита 4.
Когда вы загуглите поиск по расстоянию Хэмминга, вы заметите много вещей о его аппаратной реализации, а не о программном обеспечении. Это большой намек на то, что, возможно, любой алгоритм, который вы придумаете, должен быть реализован на C или другом скомпилированном языке.
Обновление: Рабочий код для параллельного метода
Я также предоставил упрощенный метод, помогающий с проверкой правильности, и я упаковал вариант кода Paul re для некоторых сравнений. Обратите внимание, что использование re.finditer () дает неперекрывающиеся результаты, и это может привести к тому, что совпадение с расстоянием 1 затеняет точное совпадение; см. мой последний тестовый пример.
Бит-параллельный метод имеет следующие особенности: гарантированное линейное поведение O (N), где N - длина текста. Обратите внимание, что наивный метод - O (NM), как и метод регулярных выражений (M - длина шаблона). Метод Бойера-Мура будет наихудшим случаем O (NM) и ожидаемым O (N). Также бит-параллельный метод можно легко использовать, когда ввод должен быть буферизован: он может подаваться по байту или мегабайту за раз; нет предвидения, нет проблем с границами буфера. Большое преимущество: скорость этого простого байтового кода ввода при кодировании на C.
Недостатки: длина шаблона фактически ограничена количеством бит в быстром регистре, например. 32 или 64. В этом случае длина шаблона равна 25; нет проблем. Он использует дополнительную память (S_table), пропорциональную количеству различных символов в шаблоне. В данном случае «размер алфавита» всего 4; нет проблем.
Подробности см. В этом техническом отчете. Приведен алгоритм приближенного поиска по расстоянию Левенштейна. Чтобы преобразовать в использование расстояния Хэмминга, я просто (!) Удалил части инструкции 2.1, которые обрабатывают вставку и удаление. Вы заметите много ссылок на букву «R» с надстрочным индексом «d». "d" - расстояние. Нам нужны только 0 и 1. Эти "R" становятся переменными R0 и R1 в приведенном ниже коде.
# coding: ascii
from collections import defaultdict
import re
_DEBUG = 0
# "Fast Text Searching with Errors" by Sun Wu and Udi Manber
# TR 91-11, Dept of Computer Science, University of Arizona, June 1991.
# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.20.8854
def WM_approx_Ham1_search(pattern, text):
"""Generate (Hamming_dist, start_offset)
for matches with distance 0 or 1"""
m = len(pattern)
S_table = defaultdict(int)
for i, c in enumerate(pattern):
S_table[c] |= 1 << i
R0 = 0
R1 = 0
mask = 1 << (m - 1)
for j, c in enumerate(text):
S = S_table[c]
shR0 = (R0 << 1) | 1
R0 = shR0 & S
R1 = ((R1 << 1) | 1) & S | shR0
if _DEBUG:
print "j= %2d msk=%s S=%s R0=%s R1=%s" \
% tuple([j] + map(bitstr, [mask, S, R0, R1]))
if R0 & mask: # exact match
yield 0, j - m + 1
elif R1 & mask: # match with one substitution
yield 1, j - m + 1
if _DEBUG:
def bitstr(num, mlen=8):
wstr = ""
for i in xrange(mlen):
if num & 1:
wstr = "1" + wstr
else:
wstr = "0" + wstr
num >>= 1
return wstr
def Ham_dist(s1, s2):
"""Calculate Hamming distance between 2 sequences."""
assert len(s1) == len(s2)
return sum(c1 != c2 for c1, c2 in zip(s1, s2))
def long_check(pattern, text):
"""Naively and understandably generate (Hamming_dist, start_offset)
for matches with distance 0 or 1"""
m = len(pattern)
for i in xrange(len(text) - m + 1):
d = Ham_dist(pattern, text[i:i+m])
if d < 2:
yield d, i
def Paul_McGuire_regex(pattern, text):
searchSeqREStr = (
'('
+ pattern
+ ')|('
+ ')|('.join(
pattern[:i]
+ "[ACTGN]".replace(c,'')
+ pattern[i+1:]
for i,c in enumerate(pattern)
)
+ ')'
)
searchSeqRE = re.compile(searchSeqREStr)
for match in searchSeqRE.finditer(text):
locn = match.start()
dist = int(bool(match.lastindex - 1))
yield dist, locn
if __name__ == "__main__":
genome1 = "TTTACGTAAACTAAACTGTAA"
# 01234567890123456789012345
# 1 2
tests = [
(genome1, "ACGT ATGT ACTA ATCG TTTT ATTA TTTA"),
("T" * 10, "TTTT"),
("ACGTCGTAAAA", "TCGT"), # partial match can shadow an exact match
]
nfailed = 0
for genome, patterns in tests:
print "genome:", genome
for pattern in patterns.split():
print pattern
a1 = list(WM_approx_Ham1_search(pattern, genome))
a2 = list(long_check(pattern, genome))
a3 = list(Paul_McGuire_regex(pattern, genome))
print a1
print a2
print a3
print a1 == a2, a2 == a3
nfailed += (a1 != a2 or a2 != a3)
print "***", nfailed
person
John Machin
schedule
10.03.2010
