Wyszukiwanie ciągu pozwalającego na jedną niedopasowanie w dowolnym miejscu łańcucha

Biblioteka Pythona regex obsługuje dopasowywanie rozmytych wyrażeń regularnych. Jedną z zalet w stosunku do TRE jest to, że pozwala znaleźć wszystkie dopasowania wyrażenia regularnego w tekście (obsługuje również nakładające się dopasowania).

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

Wyszukałem w googlach "Toxoplasma gondii parasite genome", aby znaleźć niektóre z tych plików genomu w Internecie. Znalazłem to, co myślę, że było blisko, plik zatytułowany " TgondiiGenomic_ToxoDB-6.0.fasta " at http://toxodb.org , rozmiar około 158Mb. Użyłem następującej ekspresji pyparsingu do wyodrębnienia sekwencji genów, zajęło to niecałe 2 minuty:

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)

(niespodzianka! niektóre sekwencje genów obejmują Run ' N ' S! O co tu chodzi?!)

Potem napisałem tę klasę jako podklasę z klasy Pyparsing Token, do wykonywania zbliżeń:

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)

Dla każdego dopasowania zwróci krotkę zawierającą rzeczywisty ciąg, który został dopasowany, oraz listę lokalizacji niedopasowania. Dokładne dopasowanie zwróci oczywiście pustą listę dla drugiej wartości. (Lubię tę klasę, myślę, że dodam ją do następnego wydania pyparsingu.)

Następnie uruchomiłem ten kod, aby wyszukać dopasowania" do-2-niedopasowania " we wszystkich sekwencjach odczytanych z .FASTA file (przypomnijmy, że genedata jest sekwencją grup ParseResults, z których każda zawiera id, długość całkowitą i ciąg sekwencji):

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

Wziąłem sekwencję wyszukiwania losowo z jednego z bitów genu, aby być pewnym, że mogę znaleźć dokładne dopasowanie, i tak z ciekawości, aby zobaczyć, ile było 1-i 2-elementowych niedopasowań.

scf_1104442825154 (964)
------------------------

scf_1104442822828 (942)
------------------------

scf_1104442824510 (987)
------------------------

scf_1104442823180 (1065)
------------------------
...

I was getting zniechęcony, Nie zobaczyć żadnych meczów do:

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

I wreszcie mój dokładny mecz na:

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

Więc chociaż to nie ustanowiło żadnych rekordów prędkości, wykonałem robotę i znalazłem kilka 2-pasujących, na wypadek, gdyby mogły być interesujące.

Dla porównania, tutaj jest wersja RE-based, który znajduje 1-niedopasowania pasuje tylko:

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

(na początku próbowałem przeszukiwać surowe źródło pliku FASTA, ale byłem zaskoczony, dlaczego tak mało pasuje w porównaniu do wersji pyparsing. Wtedy Ja zdałem sobie sprawę, że niektóre dopasowania muszą przekraczać podziały linii, ponieważ wyjście pliku fasta jest zawinięte na n znaków.)

Więc po pierwszym przejściu pyparsingu, aby wyodrębnić sekwencje genów, aby dopasować się, ten ponownie oparty Wyszukiwacz zajął około 1-1 / 2 minut, aby zeskanować wszystkie sekwencje nie-texttwrapped, aby znaleźć wszystkie te same wpisy 1-niedopasowania, które rozwiązanie pyparsing zrobił.

Możesz znaleźć różne procedury w Python-Levenshtein jakiegoś zastosowania.

>>> 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

Ten zatrzymuje pierwszy mecz, więc może być nieco szybszy, gdy jest wiele meczów

>>> 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

Dla genomu o długości 10,000,000 patrzymy na około 2,5 dnia na pojedynczy wątek do skanowania 230,000 sekwencji, więc możesz chcieć podzielić zadanie na kilka rdzeni/procesorów.

Zawsze możesz zacząć implementować wydajniejszy algorytm, gdy ten jest uruchomiony:)

Jeśli chcesz wyszukać pojedyncze upuszczone lub dodane elementy Zmień Wyrażenie regularne na this

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

To podpowiedź najdłuższy wspólny problem . Problem z podobieństwem ciągów polega na tym, że musisz przetestować ciąg ciągły 230000 sekwencji; więc jeśli porównujesz jedną z Twoich 25 sekwencji do ciągów ciągłych, uzyskasz bardzo niskie podobieństwo.

Jeśli obliczysz najdłuższą ciągłość wspólną pomiędzy 25 sekwencjami a ciągiem ciągłym, będziesz wiedział, czy jest w ciągu, jeśli długości są takie same.

Możesz zrobić trie ze wszystkich różnych sekwencji, które chcesz dopasować. Teraz jest trudną częścią tworzenia funkcji wyszukiwania głębi w dół trie, która pozwala co najwyżej na jedną niedopasowanie.

Oto jeden ze sposobów realizacji tego. Przelicz "A", "C", T", G " na 0,1,2,3 lub jego wariant. Następnie użyj krotek krotek jako struktury dla trie. Na każdym węzeł, pierwszy element w tablicy odpowiada 'A', drugi ' C ' i tak dalej. Jeśli' A ' jest gałęzią tego węzła, to istnieje kolejna krotka z 4 elementów jako pierwsza element krotki tego węzła. Jeśli nie ma gałęzi 'A', to ustaw pierwszą pozycję na 0. Na dole trie znajdują się węzły, które mają id tej sekwencji, dzięki czemu można ją umieścić na liście dopasowań.

Oto rekurencyjne funkcje wyszukiwania pozwalające na jedną niezgodność dla tego rodzaju 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)

Tutaj zaczynam poszukiwania od czegoś w rodzaju

searchnomismatch(trie,genome[ind:ind+25],0)

I addtomatches jest czymś podobnym do

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

Gdzie równe -1 oznacza, że nie było niedopasowania. W każdym razie, mam nadzieję, że wyraziłem się jasno, żebyś zrozumiał.

Próbowałem niektórych rozwiązań, ale jak już napisano, są one powolne, gdy mamy do czynienia z dużą ilością sekwencji (ciągów).

Wymyśliłem użycie bowtie i mapowanie podłańcucha zainteresowania (soi) na plik referencyjny, który zawiera ciągi znaków w formacie FASTA. Możesz podać liczbę dozwolonych niedopasowań (0..3) i dostajesz z powrotem ciągi, do których soi mapowane biorąc pod uwagę dozwolone niedopasowania. To działa dobrze i dość szybko.

Możesz użyć Pythonów wbudowanych w możliwości wyszukiwania z dopasowaniem wyrażeń regularnych.

Moduł Re w Pythonie http://docs.python.org/library/re.html

Podkład wyrażenia regularnego http://www.regular-expressions.info/

Myślę, że to może przyjść trochę późno, ale jest narzędzie o nazwie PatMaN, które robi dokładnie to, co chcesz: http://bioinf.eva.mpg.de/patman/

Możesz użyć regex matching library TRE, dla "approximate matching". Posiada również powiązania dla Pythona, Perla i Haskella.

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]

Który wyświetli

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/

Myślałem, że poniższy kod jest prosty i wygodny.

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

Możesz łatwo wstawić genomy i segmenty, które chcesz sprawdzić, a także ustawić wartość niedopasowania.

To jest dość stare, ale być może to proste rozwiązanie mogłoby zadziałać. pętli przez sekwencję biorąc 25charakter plastry. Konwertuj plasterek na tablicę numpy. Porównaj z łańcuchem 25char (również jako tablica numpy). Sum odpowiedź i jeśli odpowiedź jest 24 wydrukuj pozycję w pętli i niedopasowanie.

Te kolejne wiersze pokazują, że działa

Import numpy jako np

A = ["A", "B", "C"]

B = np.array (a)

B

Array (["A", "B", "C"], dtype="

C = ["A", "D", "C"]

D = np.array (c)

B= = d

Array ([True, False, True])

Sum (b==d)

2