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

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

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

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

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

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

Я искал в Google "геном паразита toxoplasma gondii", чтобы найти некоторые из этих файлов генома в Интернете. Я нашел то, что мне показалось очень близким, файл под названием «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 для выполнения близких совпадений:

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 символов.)

Итак, после первого прохода pyparsing для извлечения последовательностей генов для сопоставления, этот основанный на RE Затем поисковику потребовалось еще около 1-1 / 2 минуты, чтобы просканировать все последовательности без обтекания текстом, чтобы найти все те же записи с одним несоответствием, что и решение 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 дня для одного потока для сканирования 230000 последовательностей, поэтому вы можете хотите разделить задачу на несколько ядер / процессор.

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

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

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

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

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

Вот один из способов реализации этого. Преобразуйте 'A', 'C', T ', G' в 0,1,2,3 или его вариант. Затем используйте кортежи кортежей в качестве структуры для вашего дерева. В каждом узле первый элемент в вашем массиве соответствует букве «А», второй - букве «С» и так далее. Если «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)

, а дополнительные совпадения похожи на

def addtomatches(id,where=-1):
    matches.append(id,where)

, где значение -1 означает, что не было Несоответствие. В любом случае, я надеюсь, что я был достаточно ясен, чтобы вы поняли.

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

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