Поиск строки, допускающей одно несоответствие в любом месте строки - PullRequest
23 голосов
/ 10 марта 2010

Я работаю с последовательностями ДНК длиной 25 (см. Примеры ниже). У меня есть список из 230 000, и мне нужно искать каждую последовательность во всем геноме (токсоплазма, паразит gondii). Я не уверен, насколько велик геном, но намного длиннее, чем 230000 последовательностей.

Мне нужно искать каждую из моих последовательностей по 25 символов, например (AGCCTCCCATGATTGAACAGATCAT).

Геном отформатирован в виде непрерывной строки, т.е.

Мне все равно, где и сколько раз он найден, только вне зависимости от того, есть он или нет.
Это просто, я думаю -

str.find(AGCCTCCCATGATTGAACAGATCAT)

Но мне также нужно найти точное совпадение, определенное как неправильное (несоответствующее) в любом месте, но только в одном месте, и записать местоположение в последовательности. Я не уверен, как это сделать. Единственное, о чем я могу думать, это использовать подстановочный знак и выполнять поиск с подстановочным знаком в каждой позиции. То есть, поиск 25 раз.

Например,
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT

Точное совпадение с несоответствием в позиции 13.

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

Существуют программы, которые делают это - находят совпадения и частичные совпадения - но я ищу тип частичного совпадения, который невозможно обнаружить в этих приложениях.

Вот аналогичный пост для perl, хотя они сравнивают только последовательности и не ищут непрерывную строку:

Похожие записи

Ответы [ 13 ]

23 голосов
/ 11 марта 2010

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

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

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

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

Некоторые работы по параллельному битовому поиску Левенштейна с пропуском были выполнены Наварро и Раффином.ot (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
20 голосов
/ 11 июня 2012

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']
6 голосов
/ 11 марта 2010

Я гуглил «геном паразита токсоплазмы гонди», чтобы найти некоторые из этих файлов генома онлайн. Я нашел то, что, по моему мнению, было близко, файл с названием «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 символов.)

Таким образом, после первого прохода pyparsing для извлечения последовательностей генов, с которыми нужно сравнить, этому основанному на RE поисковику потребовалось еще около 1-1 / 2 минуты, чтобы отсканировать все нетекстовые последовательности, чтобы найти все одинаковые 1 -согласовать записи, которые сделал решение pyparsing.

3 голосов
/ 10 марта 2010

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

2 голосов
/ 10 марта 2010
>>> 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))))
1 голос
/ 07 февраля 2012

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

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

1 голос
/ 11 марта 2010

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

Преимущество этого метода в том, что вы просматриваете все последовательности одновременно. Это сэкономит вам много сравнений. Например, когда вы начинаете с верхнего узла и спускаетесь по ветви «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 означает, что не было несоответствия. Во всяком случае, я надеюсь, что я был достаточно ясен, чтобы вы получили картину.

1 голос
/ 10 марта 2010

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

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

0 голосов
/ 24 февраля 2018

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

б == д * +1027 *

массив ([True, False, True])

сумма (б == г) * * одна тысяча тридцать семь

2

0 голосов
/ 30 января 2016

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

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

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

...