Hasilkan semua Permutasi dengan paling banyak d Mismatch
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 Counter
untuk menemukan k-mer yang paling sering. Tapi, ini hanya berhasil jika d = 1
. Bagaimana cara menggeneralisasi ini untuk semua d <= k
?
Jawaban
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