Поиск строки, допускающей одно несоответствие в любом месте строки

Ответ 1

Прежде чем читать, просмотрели ли вы biopython?

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

Если у вас есть функция сопоставления расстояния Хэмминга (см., например, ссылку, предоставленную Ignacio), вы можете использовать ее так, чтобы выполнить поиск для первого совпадения:

any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))

но это будет довольно медленным, потому что (1) функция расстояния Хэмминга будет продолжать измельчать после второй ошибки замены (2) после отказа, она продвигает курсор на один, а не пропускает вперед в зависимости от того, что он видел (например, поиск Boyer-Moore делает).

Вы можете преодолеть (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 или другом компилированном языке.

Обновление: Рабочий код для бит-параллельного метода

Я также предоставил упрощенный метод для помощи в проверке правильности, и я собрал вариант кода Paul re для некоторых сравнений. Обратите внимание, что использование re.finditer() обеспечивает неперекрывающиеся результаты, и это может привести к тому, что совпадение расстояния 1 будет тень точно совпадающим; см. мой последний тестовый пример.

Метод бит-параллель имеет следующие функции: гарантированное линейное поведение O (N), где N - длина текста. Примечание. Наивный метод - это O (NM), как и метод регулярных выражений (M - длина шаблона). Метод Boyer-Moore был бы наихудшим случаем O (NM) и ожидаемым O (N). Также бит-параллельный метод может быть легко использован, когда вход должен быть буферизирован: он может быть загружен байтом или мегабайтом за раз; нет перспектив, никаких проблем с границами буфера. Большое преимущество: скорость этого простого кода с байтовым кодом при кодировании в C.

Нижние стороны: длина шаблона фактически ограничена количеством бит в быстром регистре, например. 32 или 64. В этом случае длина шаблона равна 25; без проблем. Он использует дополнительную память (S_table), пропорциональную количеству различных символов в шаблоне. В этом случае "размер алфавита" равен всего 4; нет проблем.

Детали из этот технический отчет. Алгоритм существует для приближенного поиска расстояния Левенштейна. Чтобы преобразовать в использование расстояния Хэмминга, я просто (!) Удалил фрагменты инструкции 2.1, которые обрабатывают вставку и удаление. Вы увидите много ссылок на "R" с надписью "d". "d" - расстояние. Нам нужны только 0 и 1. Эти "R" s становятся переменными 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

Ответ 2

Python regex библиотека поддерживает нечеткое регулярное выражение. Одним из преимуществ над TRE является то, что он позволяет находить все совпадения регулярного выражения в тексте (также поддерживает совпадающие совпадения).

import regex
m=regex.findall("AA", "CAG")
>>> []
m=regex.findall("(AA){e<=1}", "CAAG") # means allow up to 1 error
m
>>> ['CA', 'AG']

Ответ 3

Я googled для "генома генома токсоплазмы gondii", чтобы найти некоторые из этих файлов генома в Интернете. Я нашел то, что, по моему мнению, было близко, файл под названием "TgondiiGenomic_ToxoDB-6.0.fasta" на http://toxodb.org размером около 158 МБ. Для извлечения последовательностей генов я использовал следующую экспрессию пиперазина, потребовалось менее 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"! Что это такое?!)

Затем я написал этот класс как подкласс класса Token класса pyparsing для выполнения близких совпадений:

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

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

Для сравнения, здесь приведена версия на основе 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 -mismatch, которые выполнили решение pyparsing.

Ответ 4

Вы можете найти различные процедуры в Python-Levenshtein для использования.

Ответ 5

>>> 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 последовательностей, поэтому вы можете разбить задачу на несколько ядер/процессор.

Вы всегда можете начать реализацию более эффективного алгоритма, пока он работает:)

Если вы захотите найти одиночные отброшенные или добавленные элементы, измените regexp на это

>>> seq_re=re.compile('|'.join(seq[:i]+'.{0,2}'+seq[i+1:] for i in range(len(seq))))

Ответ 6

Эти подсказки самая длинная общая проблема подпоследовательности. Проблема с сходством строк здесь заключается в том, что вам нужно протестировать непрерывную строку из 230000 последовательностей; поэтому, если вы сравниваете одну из ваших 25 последовательностей с непрерывной строкой, вы получите очень низкое сходство.

Если вы вычислите самую длинную общую подпоследовательность между вашими 25 последовательностями и непрерывной строкой, вы узнаете, есть ли она в строке, если длины одинаковы.

Ответ 7

Вы можете сделать trie из всех разных последовательностей, которые вы хотите сопоставить. Теперь это сложная часть создания первой функции поиска по глубине, которая позволяет не более одного несоответствия.

Преимущество этого метода заключается в том, что вы просматриваете все последовательности сразу. Это сэкономит вам много сравнений. Например, когда вы начинаете вверху node и спускаетесь по ветке "А", вы только что сохранили себе много тысяч сравнений, потому что сразу же согласовали все последовательности, начинающиеся с "А". Для другого аргумента рассмотрим фрагмент генома, который точно соответствует заданной последовательности. Если у вас есть другая последовательность в списке последовательностей, которая отличается только в последнем символе, то использование trie только что спасло вам 23 сравнения.

Вот один из способов реализации этого. Преобразуйте 'A', 'C', T ', G' в 0,1,2,3 или вариант этого. Затем используйте кортежи кортежей в качестве вашей структуры для вашей игры. В каждом node первый элемент в вашем массиве соответствует "A" , второй - "C" и т.д. Если "A" является ветвью этого node, то в качестве первого элемента этого кортежа node имеется еще один набор из 4 элементов. Если ветка "A" не существует, установите первый элемент равным 0. В нижней части trie находятся узлы, у которых есть идентификатор этой последовательности, чтобы он мог быть помещен в список совпадений.

Здесь представлены функции рекурсивного поиска, допускающие одно несоответствие для такого рода trie:

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 означает, что не было несоответствия. Во всяком случае, я надеюсь, что я был достаточно ясен, чтобы вы получили картину.

Ответ 8

Я пробовал некоторые из решений, но, как уже было написано, они медленны при работе с большим количеством последовательностей (строк).

Я придумал использовать bowtie и сопоставляю подстроку, представляющую интерес (soi), в отношении ссылочного файла, который содержит строки в формате FASTA, Вы можете указать количество разрешенных несоответствий (0..3), и вы вернете строки, на которые нанесено soi, с учетом разрешенных несоответствий. Это работает хорошо и довольно быстро.

Ответ 9

Вы можете использовать встроенную функцию Pythons для выполнения поиска с использованием регулярного выражения.

re модуль в python http://docs.python.org/library/re.html

праймер регулярного выражения http://www.regular-expressions.info/

Ответ 10

Я думаю, это может немного запоздать, но есть инструмент под названием PatMaN, который делает именно то, что вы хотите: http://bioinf.eva.mpg.de/patman/

Ответ 11

Вы можете использовать библиотеку соответствий регулярных выражений 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

http://en.wikipedia.org/wiki/TRE_%28computing%29

http://laurikari.net/tre/

Ответ 12

Я думал, что код ниже прост и удобен.

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

Вы можете легко вставить геномы и сегменты, которые хотите проверить, а также настроить значение несоответствия.