Genere todas las permutaciones con un máximo de d desajustes

Aug 20 2020

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 Counterpara encontrar k-mer más frecuentes. Pero, esto solo funciona si d = 1. ¿Cómo generalizar esto para cualquiera d <= k?

Respuestas

2 SayanDey Aug 20 2020 at 13:24

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