Поиск строки, допускающей одно несоответствие в любом месте строки
Я работаю с последовательностями ДНК длиной 25 (см. Примеры ниже). У меня есть список из 230000, и мне нужно искать каждую последовательность во всем геноме (токсоплазма, паразит gondii). Я не уверен, насколько велик геном, но намного длиннее, чем 230000 последовательностей.
Мне нужно искать каждую из моих последовательностей по 25 символов, например (AGCCTCCCATGATTGAACAGATCAT).
Геном отформатирован как непрерывная строка, т.е. (CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT....)
Мне все равно, где и сколько раз он найден, только вне зависимости от того, есть он или нет.
Это просто, я думаю -
str.find(AGCCTCCCATGATTGAACAGATCAT)
Но мне также нужно найти точное совпадение, определенное как неправильное (несоответствующее) в любом месте, но только в одном месте, и записать местоположение в последовательности. Я не уверен, как это сделать. Единственное, о чем я могу думать, это использовать подстановочный знак и выполнять поиск с подстановочным знаком в каждой позиции. То есть искать 25 раз.
Например,
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT
Близкий матч с несовпадением в позиции 13.
Скорость не является большой проблемой, потому что я делаю это только 3 раза, хотя было бы неплохо, если бы она была быстрой.
Существуют программы, которые делают это - находят совпадения и частичные совпадения - но я ищу тип частичного совпадения, который нельзя обнаружить в этих приложениях.
Вот аналогичный пост для Perl, хотя они сравнивают только последовательности и не ищут непрерывную строку:
14 ответов
Прежде чем читать дальше, вы смотрели на биопионе?
Похоже, что вы хотите найти приблизительные совпадения с одной ошибкой замещения и нулевыми ошибками вставки / удаления, т.е. расстоянием Хэмминга, равным 1.
Если у вас есть функция сопоставления расстояний Хэмминга (см., Например, ссылку, предоставленную Ignacio), вы можете использовать ее следующим образом для поиска первого совпадения:
any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))
но это было бы довольно медленно, потому что (1) функция расстояния Хемминга продолжала бы размалывать после 2-й ошибки замещения (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.
Когда вы воспользуетесь поиском расстояния Хэмминга в Google, вы заметите множество вещей о его внедрении в аппаратном обеспечении, а не в программном обеспечении. Это большой намек на то, что, возможно, любой алгоритм, который вы придумаете, должен быть реализован на C или другом компилируемом языке.
Обновление: рабочий код для параллельного метода
Я также предоставил упрощенный метод, помогающий с проверкой правильности, и для некоторых сравнений я упаковал вариант ре-кода Пола. Обратите внимание, что использование re.finditer() дает непересекающиеся результаты, и это может привести к тому, что совпадение на расстоянии 1 затеняет точное совпадение; посмотрите мой последний контрольный пример.
Битопараллельный метод имеет следующие особенности: гарантированное линейное поведение O(N), где N - длина текста. Обратите внимание, что наивным методом является O(NM), как и метод регулярных выражений (M - длина шаблона). Метод в стиле Бойера-Мура будет иметь наихудший случай O(NM) и ожидаемый O(N). Также параллельный битовый метод может быть легко использован, когда вход должен быть буферизован: он может быть запитан байтом или мегабайтом за раз; нет забегая вперед, нет проблем с границами буфера. Большое преимущество: скорость этого простого байтового кода на входе при кодировании на C.
Недостатки: длина шаблона эффективно ограничена количеством битов в быстром регистре, например, 32 или 64. В этом случае длина шаблона составляет 25; нет проблем. Он использует дополнительную память (S_table), пропорциональную количеству различных символов в шаблоне. В этом случае "размер алфавита" составляет всего 4; нет проблем.
Подробности из этого технического отчета. Алгоритм существует для приближенного поиска на расстоянии Левенштейна. Чтобы перейти к использованию расстояния Хэмминга, я просто (!) Удалил части утверждения 2.1, которые обрабатывают вставку и удаление. Вы заметите много ссылок на "R" с верхним индексом "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
Библиотека регулярных выражений Python поддерживает нечеткое сопоставление регулярных выражений. Одним из преимуществ TRE является то, что он позволяет находить все совпадения регулярного выражения в тексте (также поддерживает перекрывающиеся совпадения).
import regex
m=regex.findall("AA", "CAG")
>>> []
m=regex.findall("(AA){e<=1}", "CAAG") # means allow up to 1 error
m
>>> ['CA', 'AG']
Я гуглил "геном паразита токсоплазмы гонди", чтобы найти некоторые из этих файлов генома онлайн. Я нашел то, что, по-моему, было близко, файл с названием "TgondiiGenomic_ToxoDB-6.0.fasta" на http://toxodb.org/, размером около 158 МБ. Я использовал следующее выражение pyparsing для извлечения последовательностей генов, это заняло чуть менее 2 минут:
fname = "TgondiiGenomic_ToxoDB-6.0.fasta"
fastasrc = open(fname).read() # yes! just read the whole dang 158Mb!
"""
Sample header:
>gb|scf_1104442823584 | organism=Toxoplasma_gondii_VEG | version=2008-07-23 | length=1448
"""
integer = Word(nums).setParseAction(lambda t:int(t[0]))
genebit = Group(">gb|" + Word(printables)("id") + SkipTo("length=") +
"length=" + integer("genelen") + LineEnd() +
Combine(OneOrMore(Word("ACGTN")),adjacent=False)("gene"))
# read gene data from .fasta file - takes just under a couple of minutes
genedata = OneOrMore(genebit).parseString(fastasrc)
(Сюрприз! Некоторые из последовательностей генов включают серии "N"! Что за чертовщина?!)
Затем я написал этот класс как подкласс класса pyparsing Token для выполнения близких совпадений:
class CloseMatch(Token):
def __init__(self, seq, maxMismatches=1):
super(CloseMatch,self).__init__()
self.name = seq
self.sequence = seq
self.maxMismatches = maxMismatches
self.errmsg = "Expected " + self.sequence
self.mayIndexError = False
self.mayReturnEmpty = False
def parseImpl( self, instring, loc, doActions=True ):
start = loc
instrlen = len(instring)
maxloc = start + len(self.sequence)
if maxloc <= instrlen:
seq = self.sequence
seqloc = 0
mismatches = []
throwException = False
done = False
while loc < maxloc and not done:
if instring[loc] != seq[seqloc]:
mismatches.append(seqloc)
if len(mismatches) > self.maxMismatches:
throwException = True
done = True
loc += 1
seqloc += 1
else:
throwException = True
if throwException:
exc = self.myException
exc.loc = loc
exc.pstr = instring
raise exc
return loc, (instring[start:loc],mismatches)
Для каждого совпадения будет возвращаться кортеж, содержащий фактическую строку, которая была сопоставлена, и список местоположений несоответствия. Точные совпадения, конечно, вернут пустой список для второго значения. (Мне нравится этот класс, я думаю, я добавлю его в следующий выпуск pyparsing.)
Затем я запустил этот код для поиска совпадений "до-2-несоответствия" во всех последовательностях, считанных из файла.fasta (напомним, что genedata - это последовательность групп ParseResults, каждая из которых содержит идентификатор, целую длину и строка последовательности):
searchseq = CloseMatch("ATCATCGAATGGAATCTAATGGAAT", 2)
for g in genedata:
print "%s (%d)" % (g.id, g.genelen)
print "-"*24
for t,startLoc,endLoc in searchseq.scanString(g.gene):
matched, mismatches = t[0]
print "MATCH:", searchseq.sequence
print "FOUND:", matched
if mismatches:
print " ", ''.join(' ' if i not in mismatches else '*'
for i,c in enumerate(searchseq.sequence))
else:
print "<exact match>"
print "at location", startLoc
print
print
Я взял случайную последовательность поиска из одного из битов гена, чтобы быть уверенным, что смог найти точное совпадение, и просто из любопытства посмотрел, сколько было 1- и 2-элементных несовпадений.
Это заняло немного времени, чтобы бежать. Через 45 минут у меня был этот вывод, перечисляющий каждый идентификатор и длину гена, и любые найденные частичные совпадения:
scf_1104442825154 (964)
------------------------
scf_1104442822828 (942)
------------------------
scf_1104442824510 (987)
------------------------
scf_1104442823180 (1065)
------------------------
...
Я был обескуражен, чтобы не видеть никаких матчей, пока:
scf_1104442823952 (1188)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAACGGAATCGAATGGAAT
* *
at location 33
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 175
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 474
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 617
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATAGAAT
* *
at location 718
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGATTCGAATGGAAT
* *
at location 896
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGTAT
* *
at location 945
И наконец мое точное соответствие:
scf_1104442823584 (1448)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGACTCGAATGGAAT
* *
at location 177
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCAAATGGAAT
*
at location 203
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
* *
at location 350
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAA
* *
at location 523
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
* *
at location 822
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCTAATGGAAT
<exact match>
at location 848
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCGTCGAATGGAGTCTAATGGAAT
* *
at location 969
Так что, хотя это не устанавливало никаких рекордов скорости, я выполнил свою работу и нашел несколько 2-х матчей, на случай, если они могут быть интересны.
Для сравнения приведена версия на основе RE, которая находит только совпадения с 1 несоответствием:
import re
seqStr = "ATCATCGAATGGAATCTAATGGAAT"
searchSeqREStr = seqStr + '|' + \
'|'.join(seqStr[:i]+"[ACTGN]".replace(c,'') +seqStr[i+1:]
for i,c in enumerate(seqStr))
searchSeqRE = re.compile(searchSeqREStr)
for g in genedata:
print "%s (%d)" % (g.id, g.genelen)
print "-"*24
for match in searchSeqRE.finditer(g.gene):
print "MATCH:", seqStr
print "FOUND:", match.group(0)
print "at location", match.start()
print
print
(Сначала я попытался найти сам исходный файл необработанного файла FASTA, но был озадачен, почему так мало совпадений по сравнению с версией pyparsing. Затем я понял, что некоторые совпадения должны пересекать разрывы строк, поскольку вывод файла fasta переносится в n персонажи.)
Таким образом, после первого прохода синтаксического анализа для извлечения последовательностей генов для сравнения этому основанному на RE поисковику потребовалось еще около 1-1/2 минуты для сканирования всех нетекстовых последовательностей, чтобы найти все одинаковые записи с 1 несоответствием что сделал решение pyparsing.
>>> import re
>>> seq="AGCCTCCCATGATTGAACAGATCAT"
>>> genome = "CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT..."
>>> seq_re=re.compile('|'.join(seq[:i]+'.'+seq[i+1:] for i in range(len(seq))))
>>> seq_re.findall(genome) # list of matches
[]
>>> seq_re.search(genome) # None if not found, otherwise a match object
Это останавливает первое совпадение, поэтому может быть немного быстрее, когда есть несколько совпадений
>>> print "found" if any(seq_re.finditer(genome)) else "not found"
not found
>>> print "found" if seq_re.search(genome) else "not found"
not found
>>> seq="CAT"
>>> seq_re=re.compile('|'.join(seq[:i]+'.'+seq[i+1:] for i in range(len(seq))))
>>> print "found" if seq_re.search(genome) else "not found"
found
для генома длиной 10 000 000 вы ищете около 2,5 дней для одного потока, чтобы отсканировать 230 000 последовательностей, поэтому вы можете разделить задачу на несколько ядер / процессоров.
Вы всегда можете начать реализовывать более эффективный алгоритм, пока он работает:)
Если вы хотите искать отдельные удаленные или добавленные элементы, измените регулярное выражение на это
>>> seq_re=re.compile('|'.join(seq[:i]+'.{0,2}'+seq[i+1:] for i in range(len(seq))))
Вы можете найти различные подпрограммы в Python-Levenshtein для некоторого использования.
Это указывает на самую длинную общую проблему подпоследовательности. Проблема схожести строк здесь заключается в том, что вам нужно проверить непрерывную строку из 230000 последовательностей; так что если вы сравниваете одну из своих 25 последовательностей с непрерывной строкой, вы получите очень низкое сходство.
Если вы вычислите самую длинную общую подпоследовательность между своими 25 последовательностями и непрерывной строкой, вы узнаете, находится ли она в строке, если длины одинаковы.
Я попробовал некоторые из решений, но, как уже было написано, они работают медленно при работе с большим количеством последовательностей (строк).
Я придумал использовать bowtie и сопоставить интересующую подстроку (soi) с эталонным файлом, который содержит строки в формате FASTA. Вы можете указать количество допустимых несоответствий (0..3) и получить обратно строки, на которые сопоставлены данные с учетом допустимых несоответствий. Это работает хорошо и довольно быстро.
Вы могли бы сделать три из всех различных последовательностей, которые вы хотите сопоставить. Теперь самое сложное - создать функцию поиска в глубину, которая допускает самое большее одно несоответствие.
Преимущество этого метода в том, что вы просматриваете все последовательности одновременно. Это сэкономит вам много сравнений. Например, когда вы начинаете с верхнего узла и спускаетесь по ветви "A", вы просто спасли себе тысячи сопоставлений, потому что просто мгновенно сопоставили все последовательности, начинающиеся с "A". Для другого аргумента рассмотрим фрагмент генома, который точно соответствует данной последовательности. Если у вас есть другая последовательность в вашем списке последовательностей, которая отличается только последним символом, то использование дерева только сэкономило вам 23 операции сравнения.
Вот один из способов реализации этого. Преобразовать 'A','C',T',G' в 0,1,2,3 или вариант этого. Затем используйте кортежи кортежей в качестве структуры для вашего дерева. В каждом узле первый элемент в вашем массиве соответствует "A", второй - "C" и так далее. Если "A" является ветвью этого узла, то есть еще один кортеж из 4 элементов в качестве первого элемента кортежа этого узла. Если ветки "A" нет, то установите для первого элемента значение 0. В нижней части дерева находятся узлы, которые имеют идентификатор этой последовательности, чтобы его можно было включить в список совпадений.
Вот рекурсивные функции поиска, допускающие одно несоответствие для этого вида дерева:
def searchnomismatch(node,genome,i):
if i == 25:
addtomatches(node)
else:
for x in range(4):
if node[x]:
if x == genome[i]:
searchnomismatch(node[x],genome,i+1)
else:
searchmismatch(node[x],genome,i+1,i)
def searchmismatch(node,genome,i,where):
if i == 25:
addtomatches(node,where)
else:
if node[genome[i]]:
searchmismatch(node[genome[i]],genome,i+1,where)
Здесь я начинаю поиск с чем-то вроде
searchnomismatch(trie,genome[ind:ind+25],0)
и addtomatches это что-то похожее на
def addtomatches(id,where=-1):
matches.append(id,where)
где равно -1 означает, что не было несоответствия. Во всяком случае, я надеюсь, что я был достаточно ясен, чтобы вы получили картину.
Я тоже столкнулся с этой проблемой, но предложилregexpackage ведет себя не так, как ожидалось (с буквой «e» в регулярном выражении), особенно при вставке/удалении в запросе и последовательности поиска. Например, не удалось найти совпадение:regex.search("xxgyy" + "{e<=1}", "xxggyy").
Наконец-то я нашел пакет, который работает! find_near_matches()вfuzzysearchупаковка! https://github.com/taleinat/fuzzysearch
Он использует расстояние Левенштейна (max_l_dist).
>>> from fuzzysearch import find_near_matches
# search for 'PATTERN' with a maximum Levenshtein Distance of 1
>>> find_near_matches('PATTERN', '---PATERN---', max_l_dist=1)
[Match(start=3, end=9, dist=1, matched="PATERN")]
если вы хотите получить возвращаемое значение:
>>> x= find_near_matches('PATTERN', '---PATERN---', max_l_dist=1)
## get the start position
>>> x[0].start
3
## get the matched string
>>> x[0].matched
'PATERN'
Вы можете использовать встроенную возможность Pythons для поиска по регулярному выражению.
Повторный модуль в Python http://docs.python.org/library/re.html
учебник по регулярным выражениям http://www.regular-expressions.info/
Я думал, что код ниже прост и удобен.
in_pattern = "";
in_genome = "";
in_mistake = d;
out_result = ""
kmer = len(in_pattern)
def FindMistake(v):
mistake = 0
for i in range(0, kmer, 1):
if (v[i]!=in_pattern[i]):
mistake+=1
if mistake>in_mistake:
return False
return True
for i in xrange(len(in_genome)-kmer+1):
v = in_genome[i:i+kmer]
if FindMistake(v):
out_result+= str(i) + " "
print out_result
Вы можете легко вставить геномы и сегменты, которые вы хотите проверить, а также установить значение несоответствия.
Вы можете использовать библиотеку соответствия TRE для регулярного соответствия. Он также имеет привязки для Python, Perl и Haskell.
import tre
pt = tre.compile("Don(ald)?( Ervin)? Knuth", tre.EXTENDED)
data = """
In addition to fundamental contributions in several branches of
theoretical computer science, Donnald Erwin Kuth is the creator of
the TeX computer typesetting system, the related METAFONT font
definition language and rendering system, and the Computer Modern
family of typefaces.
"""
fz = tre.Fuzzyness(maxerr = 3)
print fz
m = pt.search(data, fz)
if m:
print m.groups()
print m[0]
который будет выводить
tre.Fuzzyness(delcost=1,inscost=1,maxcost=2147483647,subcost=1, maxdel=2147483647,maxerr=3,maxins=2147483647,maxsub=2147483647)
((95, 113), (99, 108), (102, 108))
Donnald Erwin Kuth
Это довольно старое, но, возможно, это простое решение может сработать. цикл через последовательность, берущую 25 знаков характера. преобразовать ломтик в массив NumPy. Сравните со строкой 25char (также в виде пустого массива). Суммируйте ответ, и если ответ 24, распечатайте позицию в цикле и несоответствие.
следующие несколько строк показывают, что это работает
импортировать numpy как np
a = ['A', 'B', 'C']
b = np.array (a)
б
массив (['A', 'B', 'C'], dtype = '
c = ['A', 'D', 'C']
d = np.array (c)
б == д
массив ([ True, False, True])
Сумма (б == д)
2
Я думаю, что это может быть немного поздно, но есть инструмент под названием PatMaN, который делает именно то, что вы хотите: http://bioinf.eva.mpg.de/patman/