Hasilkan semua Permutasi dengan paling banyak d Mismatch

Aug 20 2020

Saya sedang memecahkan masalah pencocokan pola dengan jarak hamming hingga d untuk urutan DNA. Regex menyelamatkan saya di sana. Tetapi sekarang saya mengalami masalah yang berbeda. Mengingat urutan DNA yang panjang, saya menemukan k-mer yang paling sering tidak cocok dengan paling banyak d ketidakcocokan. Di sini, k-mer mengacu pada sub-urutan panjang k.

Catatan: Urutan DNA dapat direpresentasikan hanya dengan menggunakan empat huruf: {A, C, G, T}

Sebagai contoh,

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

Di sini, hanya mungkin dua k-mers: "AGG", "GGC"

Sekarang, saya dapat mengubah satu per satu dengan 1 ketidakcocokan dengan menjalankan kode berikut untuk "AGG" dan "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" akan memberikan:

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

Dan, "GCC" akan memberikan:

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

Kemudian saya dapat menggunakan Counteruntuk menemukan k-mer yang paling sering. Tapi, ini hanya berhasil jika d = 1. Bagaimana cara menggeneralisasi ini untuk semua d <= k?

Jawaban

2 SayanDey Aug 20 2020 at 13:24

Ini adalah metode yang mahal secara komputasi. Tapi ya harus mengambil yang diinginkan. Apa yang saya lakukan di sini adalah menghitung semua ketidakcocokan dengan hamming dist 1. kemudian menghitung ketidakcocokan baru dengan ham dist 1 dari sebelumnya mismatch dan recurse sampai 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)]))

menemukan kode kinerja yang lebih baik hampir 10-20X lebih cepat

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

Waktu CPU: pengguna 1 menit 32 detik, sys: 1.81 detik, total: 1 menit 34 detik Waktu dinding: 1 menit 34 detik