font - subplot python title
Busque una cadena que permita una falta de coincidencia en cualquier ubicación de la cadena (13)
Busqué en Google "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 , con un tamaño de aproximadamente 158Mb. Usé la siguiente expresión pyparsing para extraer las secuencias de genes, tardó un 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 no coincidentes. 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 faltas de coincidencia" 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 "<exact match>"
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
<exact match>
at location 848
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCGTCGAATGGAGTCTAATGGAAT
* *
at location 969
Entonces, 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.)
Así que 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 discrepancia que hizo la solución pyparsing.
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 de 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 (CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGGGGGCCGTGGGAGGCTT ...
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 programas que hacen esto: encontrar coincidencias y coincidencias parciales, pero estoy buscando un tipo de coincidencia parcial que no se pueda detectar con estas aplicaciones.
Aquí hay una publicación similar para perl, aunque solo están comparando secuencias y no buscando una cadena continua:
Esto es bastante antiguo, pero quizás esta simple solución podría funcionar. recorrer la secuencia tomando 25 cortes de 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])
suma (b == d)
2
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á comparando una de sus 25 secuencias con la cadena continua, obtendrá 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.
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'']
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.
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/
Podrías hacer un trie 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 comienza en el nodo superior y baja por la rama ''A'', acaba de ahorrarse muchos miles de comparaciones, ya que se ha 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 usar 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 rama de este nodo, entonces hay otra tupla de 4 elementos como el primer elemento de la tupla de este nodo. Si no hay una rama ''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í hay 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ó usar bowtie y mapear la subcadena de interés (soi) contra 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.
Puede encontrar las diversas rutinas en Python-Levenshtein de algún uso.
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
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/
Antes de biopython 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 en función de 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 búsquedas Levenshtein aproximadas de bits y saltos, 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 la 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 compilado.
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 expresión regular (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 declaración 2.1 que manejan la inserción y la 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
>>> 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))))