Genere todas las permutaciones con un máximo de d desajustes
Estaba resolviendo un problema de coincidencia de patrones con una distancia de Hamming hasta d para una secuencia de ADN. Regex me salvó allí. Pero ahora me he encontrado con un problema diferente. Dada una secuencia de ADN larga, tengo que encontrar los k-mer desajustados más frecuentes con como máximo d desajustes. Aquí, k-mer se refiere a una subsecuencia de longitud k.
Nota: la secuencia de ADN se puede representar usando solo cuatro letras: {A, C, G, T}
Por ejemplo,
DNA sequence= "AGGC"
k = 3
d = 1
Aquí, solo son posibles dos k-mers: "AGG", "GGC"
Ahora, puedo permutarlos individualmente con 1 discrepancias ejecutando el siguiente código para "AGG" y "GGC"
def permute_one_nucleotide(motif, alphabet={"A", "C", "G", "T"}):
import itertools
return list(
set(
itertools.chain.from_iterable(
[
[
motif[:pos] + nucleotide + motif[pos + 1 :]
for nucleotide in alphabet
]
for pos in range(len(motif))
]
)
)
)
"AGG" dará:
['TGG', 'ATG', 'AGG', 'GGG', 'AGT', 'CGG', 'AGC', 'AGA', 'ACG', 'AAG']
Y, "GCC" dará:
['GCC', 'GAC', 'GGT', 'GGA', 'AGC', 'GTC', 'TGC', 'CGC', 'GGG', 'GGC']
Entonces puedo usar Counter
para encontrar k-mer más frecuentes. Pero, esto solo funciona si d = 1
. ¿Cómo generalizar esto para cualquiera d <= k
?
Respuestas
Este es un método computacionalmente costoso. Pero sí, debería buscar lo deseado. Lo que hice aquí fue calcular todos los desajustes con hamming dist 1. luego calcular el nuevo desajuste con ham dist 1 a partir del anterior desajuste y recurrir hasta d.
import itertools
all_c=set('AGCT')
other = lambda x : list(all_c.difference(x))
def get_changed(sub, i):
return [sub[0:i]+c+sub[i+1:] for c in other(sub[i])]
def get_mismatch(d, setOfMmatch):
if d==0:
return setOfMmatch
newMmatches=[]
for sub in setOfMmatch:
newMmatches.extend(list(map(lambda x : ''.join(x), itertools.chain.from_iterable(([get_changed(sub, i) for i, c in enumerate(sub)])))))
setOfMmatch=setOfMmatch.union(newMmatches)
return get_mismatch(d-1, setOfMmatch)
dna='AGGC'
hamm_dist=1
length=3
list(itertools.chain.from_iterable([get_mismatch(hamm_dist, {dna[i:i+length]}) for i in range(len(dna)-length+1)]))
# without duplicates
# set(itertools.chain.from_iterable([get_mismatch(hamm_dist, {dna[i:i+length]}) for i in range(len(dna)-length+1)]))
encontró un código de mejor rendimiento casi 10-20 veces más rápido
%%time
import itertools, random
from cacheout import Cache
import time
all_c=set('AGCT')
get_other = lambda x : list(all_c.difference(x))
other={}
for c in all_c:
other[c]=get_other(c)
def get_changed(sub, i):
return [sub[0:i]+c+sub[i+1:] for c in other[sub[i]]]
mmatchHash=Cache(maxsize=256*256, ttl=0, timer=time.time, default=None)
def get_mismatch(d, setOfMmatch):
if d==0:
return setOfMmatch
newMmatches=[]
for sub in setOfMmatch:
newMmatches.extend(list(map(lambda x : ''.join(x), itertools.chain.from_iterable(([get_changed(sub, i) for i, c in enumerate(sub)])))))
setOfMmatch=setOfMmatch.union(newMmatches)
if not mmatchHash.get((d-1, str(setOfMmatch)), 0):
mmatchHash.set((d-1, str(setOfMmatch)), get_mismatch(d-1, setOfMmatch))
return mmatchHash.get((d-1, str(setOfMmatch)))
length_of_DNA=1000
dna=''.join(random.choices('AGCT', k=length_of_DNA))
hamm_dist=4
length=9
len(list(itertools.chain.from_iterable([get_mismatch(hamm_dist, {dna[i:i+length]}) for i in range(len(dna)-length+1)])))
# set(itertools.chain.from_iterable([get_mismatch(hamm_dist, {dna[i:i+length]}) for i in range(len(dna)-length+1)]))
Tiempos de CPU: usuario 1 min 32 s, sys: 1,81 s, total: 1 min 34 s Tiempo de pared: 1 min 34 s