Busque una cadena que permita una falta de coincidencia en cualquier ubicación de la cadena

Estoy trabajando con secuencias de ADN de longitud 25 (ver ejemplos a continuación). Tengo una lista de 230,000 y necesito buscar cada secuencia en todo el genoma (parásito toxoplasma gondii). No estoy seguro de cuán grande es el genoma, pero es mucho más largo que 230,000 secuencias.

Necesito buscar cada una de mis secuencias de 25 caracteres, por ejemplo, (AGCCTCCCATGATTGAACAGATCAT).

El genoma está formateado como una cadena continua, es decir (CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGGGGGCCGGGGGGGTTG …

No me importa dónde ni cuántas veces se encuentra, solo si lo es o no.
Esto es simple, creo.

str.find(AGCCTCCCATGATTGAACAGATCAT) 

Pero también tengo que encontrar una coincidencia cercana definida como incorrecta (no coincidente) en cualquier ubicación, pero solo en una ubicación, y registrar la ubicación en la secuencia. No estoy seguro de cómo hacer esto. Lo único que se me ocurre es usar un comodín y realizar la búsqueda con un comodín en cada posición. Es decir, buscar 25 veces.

Por ejemplo,
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT

Un partido cercano con un desajuste en la posición 13.

La velocidad no es un gran problema porque solo lo hago 3 veces, aunque sería bueno si fuera rápido.

Hay progtwigs que hacen esto, encontrar coincidencias y coincidencias parciales, pero estoy buscando un tipo de coincidencia parcial que no sea detectable con estas aplicaciones.

Aquí hay una publicación similar para perl, aunque solo están comparando secuencias y no buscando una cadena continua:

Publicación relacionada

Antes de seguir leyendo , ¿has mirado biopython ?

Parece que desea encontrar coincidencias aproximadas con un error de sustitución y cero errores de inserción / eliminación, es decir, una distancia de Hamming de 1.

Si tiene una función de concordancia de distancia de Hamming (ver, por ejemplo, el enlace provisto por Ignacio), puede usarla así para hacer una búsqueda de la primera concordancia:

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

pero esto sería bastante lento, porque (1) la función de distancia de Hamming continuaría moliendo después del segundo error de sustitución (2) después de la falla, avanza el cursor en uno en lugar de saltar hacia adelante según lo que vio (como un Boyer) La busqueda de moore lo hace).

Puedes superar (1) con una función como esta:

 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 

Nota: intencionalmente no es Pythonic, es Cic, porque necesitarías usar C (quizás a través de Cython) para obtener una velocidad razonable.

Navarro y Raffinot (google “Navarro Raffinot nrgrep”) han realizado algunos trabajos en las búsquedas de Levenshtein aproximadas de bits paralelos y esto podría adaptarse a las búsquedas de Hamming. Tenga en cuenta que los métodos de bit paralelo tienen limitaciones en la longitud de la cadena de consulta y el tamaño del alfabeto, pero los suyos son 25 y 4 respectivamente, por lo que no hay problemas. Actualización: omitir probablemente no sea de mucha ayuda con un tamaño de alfabeto de 4.

Cuando busques en Google para la búsqueda a distancia de Hamming, notarás muchas cosas sobre su implementación en hardware y no mucho en software. Este es un gran indicio de que tal vez cualquier algoritmo que surja debería implementarse en C o en algún otro lenguaje comstackdo.

Actualización: Código de trabajo para un método de bit paralelo.

También proporcioné un método simplista para ayudar con la verificación de la corrección, y empaqueté una variación del código de Paul para algunas comparaciones. Tenga en cuenta que el uso de re.finditer () proporciona resultados no superpuestos, y esto puede hacer que una coincidencia de distancia-1 impida una coincidencia exacta; ver mi último caso de prueba

El método de bit paralelo tiene estas características: comportamiento lineal garantizado O (N) donde N es la longitud del texto. Tenga en cuenta que el método ingenuo es O (NM) al igual que el método de expresiones regulares (M es la longitud del patrón). Un método estilo Boyer-Moore sería el peor de los casos O (NM) y el esperado O (N). Además, el método de bit paralelo puede usarse fácilmente cuando la entrada tiene que ser almacenada en búfer: se puede alimentar un byte o un megabyte a la vez; Sin mirar hacia adelante, no hay problemas con los límites del búfer. La gran ventaja: la velocidad de ese código de byte de entrada simple cuando se codifica en C.

Desventajas: la longitud del patrón se limita efectivamente al número de bits en un registro rápido, por ejemplo, 32 o 64. En este caso, la longitud del patrón es 25; No hay problema. Utiliza memoria extra (S_table) proporcional al número de caracteres distintos en el patrón. En este caso, el “tamaño del alfabeto” es solo 4; No hay problema.

Detalles de este informe técnico . El algoritmo que hay para la búsqueda aproximada en la distancia de Levenshtein. Para convertir a usar la distancia de Hamming, simplemente (!) Quité los fragmentos de la statement 2.1 que manejan la inserción y eliminación. Notará muchas referencias a “R” con un superíndice “d”. “d” es la distancia. Solo necesitamos 0 y 1. Estas “R” se convierten en las variables R0 y R1 en el código siguiente.

 # 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 

La biblioteca de expresiones regulares de Python admite la coincidencia de expresiones regulares difusas. Una ventaja sobre TRE es que permite encontrar todas las coincidencias de expresiones regulares en el texto (también admite coincidencias superpuestas).

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

Busqué en Google el “genoma del parásito de toxoplasma gondii” para encontrar algunos de estos archivos del genoma en línea. Encontré lo que creo que estaba cerca, un archivo titulado “TgondiiGenomic_ToxoDB-6.0.fasta” en http://toxodb.org , de unos 158Mb de tamaño. Utilicé la siguiente expresión pyparsing para extraer las secuencias de genes, tomó poco menos de 2 minutos:

 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) 

(¡Sorpresa! ¡Algunas de las secuencias de genes incluyen series de ‘N’s! ¡¿Qué diablos es eso ?!)

Luego escribí esta clase como una subclase de la clase Token de reproducción, para hacer coincidencias cercanas:

 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) 

Para cada coincidencia, esto devolverá una tupla que contiene la cadena real que coincidió y una lista de las ubicaciones que no coinciden. Las coincidencias exactas, por supuesto, devolverían una lista vacía para el segundo valor. (Me gusta esta clase, creo que la agregaré a la próxima versión de pyparsing).

Luego ejecuté este código para buscar coincidencias de “hasta 2 discrepancias” en todas las secuencias leídas del archivo .fasta (recuerde que genedata es una secuencia de grupos ParseResults, cada uno de los cuales contiene un id, una longitud entera y una secuencia de secuencia):

 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 "" print "at location", startLoc print print 

Tomé la secuencia de búsqueda al azar de uno de los bits del gen, para asegurarme de que podía encontrar una coincidencia exacta, y solo por curiosidad para ver cuántos desajustes de 1 y 2 elementos había.

Esto tomó un poco de tiempo para correr. Después de 45 minutos, tuve esta salida, enumerando cada id y la longitud del gen, y cualquier coincidencia parcial encontrada:

 scf_1104442825154 (964) ------------------------ scf_1104442822828 (942) ------------------------ scf_1104442824510 (987) ------------------------ scf_1104442823180 (1065) ------------------------ ... 

Me estaba desanimando, por no ver ninguna coincidencia hasta que:

 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 

Y finalmente mi coincidencia exacta en:

 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  at location 848 MATCH: ATCATCGAATGGAATCTAATGGAAT FOUND: ATCGTCGAATGGAGTCTAATGGAAT * * at location 969 

Así que, aunque esto no estableció ningún récord de velocidad, terminé el trabajo y también encontré 2 partidos, en caso de que puedan ser de interés.

A modo de comparación, aquí hay una versión basada en RE, que solo encuentra coincidencias de 1-coincidencia:

 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 

(Al principio, intenté buscar la fuente de archivos FASTA sin procesar, pero me desconcertó el porqué de tan pocas coincidencias en comparación con la versión de pyparsing. Luego me di cuenta de que algunas de las coincidencias deben cruzar los saltos de línea, ya que la salida del archivo fasta está ajustada a n caracteres.)

Entonces, después de la primera fase de creación de archivos para extraer las secuencias de genes con las que coincidir, este buscador basado en RE luego tomó otros 1-1 / 2 minutos para escanear todas las secuencias envueltas sin texto, para encontrar todas las mismas entradas de 1 no coincidencia que hizo la solución pyparsing.

Puede encontrar las diversas rutinas en Python-Levenshtein de algún uso.

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

Este detiene la primera partida, por lo que puede ser un poco más rápido cuando hay varias coincidencias

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

para un genoma de 10.000.000 de longitud, está buscando alrededor de 2,5 días para que un solo hilo escanee 230,000 secuencias, por lo que puede dividir la tarea en unos pocos núcleos / CPU.

Siempre puedes comenzar a implementar un algoritmo más eficiente mientras este se está ejecutando 🙂

Si desea buscar elementos sueltos o agregados, cambie la expresión regular a este

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

Esto sugiere el problema más largo de subsecuencias comunes . El problema con la similitud de cadenas aquí es que necesita probar contra una cadena continua de 230000 secuencias; así que si estás comparando una de tus 25 secuencias con la cadena continua obtendrás una similitud muy baja.

Si calcula la subsecuencia común más larga entre sus 25 secuencias y la cadena continua, sabrá si está en la cadena si las longitudes son las mismas.

Podrías hacer un trío de todas las diferentes secuencias que deseas emparejar. Ahora es la parte difícil de hacer una primera función de búsqueda profunda en el trie que permita a lo sumo una falta de coincidencia.

La ventaja de este método es que está buscando a través de todas las secuencias a la vez. Esto te ahorrará muchas comparaciones. Por ejemplo, cuando comienzas en el nodo superior y bajas por la twig ‘A’, acabas de guardarte muchos miles de comparaciones porque se han combinado instantáneamente con todas las secuencias que comienzan con ‘A’. Para un argumento diferente, considere una porción del genoma que coincida exactamente con una secuencia dada. Si tiene una secuencia diferente en su lista de secuencias que difiere solo en el último símbolo, entonces el uso de un trie le ha guardado 23 operaciones de comparación.

Aquí hay una forma de implementar esto. Convierta ‘A’, ‘C’, T ‘, G’ a 0,1,2,3 o una variante de eso. Luego usa tuplas de tuplas como tu estructura para tu trie. En cada nodo, el primer elemento de su matriz se corresponde con ‘A’, el segundo con ‘C’ y así sucesivamente. Si ‘A’ es una twig de este nodo, entonces hay otra tupla de 4 elementos como el primer elemento de la tupla de este nodo. Si no hay una twig ‘A’, establezca el primer elemento en 0. En la parte inferior de la lista se encuentran los nodos que tienen el ID de esa secuencia para que se pueda incluir en la lista de coincidencias.

Aquí están las funciones de búsqueda recursiva que permiten una falta de coincidencia para este tipo de 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) 

Aquí, empiezo la búsqueda con algo como

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

y addtomatches es algo similar a

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

donde igual a -1 significa que no hubo un desajuste. De todos modos, espero que haya sido lo suficientemente claro para que se haga una idea.

Probé algunas de las soluciones, pero como ya se ha escrito, son lentas cuando se trata de una gran cantidad de secuencias (cadenas).

Se me ocurrió utilizar Bowtie y mapear la subcadena de interés (soi) con un archivo de referencia que contiene las cadenas en formato FASTA. Puede proporcionar el número de desajustes permitidos (0..3) y volver a obtener las cadenas a las que se asignó el soi dados los desajustes permitidos. Esto funciona bien y bastante rápido.

Podría usar la capacidad incorporada de Pythons para realizar la búsqueda con la coincidencia de expresiones regulares.

re módulo en python http://docs.python.org/library/re.html

cebador de expresiones regulares http://www.regular-expressions.info/

Supongo que esto puede llegar un poco tarde, pero hay una herramienta llamada PatMaN que hace exactamente lo que quiere: http://bioinf.eva.mpg.de/patman/

Puede utilizar la biblioteca de coincidencias de expresiones regulares TRE, para “coincidencia aproximada”. También tiene enlaces para Python, Perl y 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] 

que dará salida

 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/

Pensé que el código de abajo es simple y conveniente.

 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 

Puede insertar fácilmente los genomas y segmentos que desea verificar y también configurar el valor de la falta de coincidencia.

Esto es bastante antiguo, pero quizás esta simple solución podría funcionar. Recorre la secuencia tomando rebanadas de 25 caracteres. convertir el sector en una matriz numpy. Compare con la cadena de 25 caracteres (también como una matriz numpy). Sume la respuesta y, si la respuesta es 24, imprima la posición en el bucle y la discrepancia.

Las siguientes líneas lo muestran funcionando.

importar numpy como np

a = [‘A’, ‘B’, ‘C’]

b = np.array (a)

segundo

array ([‘A’, ‘B’, ‘C’], dtype = ‘

c = [‘A’, ‘D’, ‘C’]

d = np.array (c)

b == d

array ([Verdadero, Falso, Verdadero])

sum (b == d)

2