Générer toutes les permutations avec au plus d discordances

Question

Je résolvais un problème de correspondance de motif avec une distance de frappe jusqu'à d pour une séquence d'ADN. Regex m'a sauvé là-bas. Mais maintenant, j'ai rencontré un problème différent. Étant donné une longue séquence d'ADN, je dois trouver le k-mer le plus souvent dépareillé avec au plus d discordances. Ici, k-mer fait référence à une sous-séquence de longueur k.

Remarque: la séquence d'ADN ne peut être représentée qu'en utilisant quatre lettres: {A, C, G, T}

Par exemple,

DNA sequence= "AGGC"
k = 3
d = 1

Ici, seuls deux k-mers sont possibles: "AGG", "GGC"

Maintenant, je peux les permuter individuellement avec 1 incompatibilités en exécutant le code suivant pour "AGG" et "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" donnera:

['TGG', 'ATG', 'AGG', 'GGG', 'AGT', 'CGG', 'AGC', 'AGA', 'ACG', 'AAG']

Et, "GCC" donnera:

['GCC', 'GAC', 'GGT', 'GGA', 'AGC', 'GTC', 'TGC', 'CGC', 'GGG', 'GGC']

Ensuite, je peux utiliser Counterpour trouver les k-mer les plus fréquents. Mais cela ne fonctionne que si d = 1. Comment généraliser cela pour tout d <= k?

Answer 1

C'est une méthode coûteuse en calcul. Mais oui devrait aller chercher désiré. Ce que j'ai fait ici, c'est de calculer toutes les discordances avec hamming dist 1. puis calculer la nouvelle discordance avec ham dist 1 à partir de la discordance précédente et recurse jusqu'à 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)]))

trouvé un meilleur code de performance presque 10 à 20 fois plus rapide

%%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)]))

Temps CPU: utilisateur 1 min 32 s, sys: 1,81 s, total: 1 min 34 s Temps mural: 1 min 34 s