최대 d 개의 불일치가있는 모든 순열 생성

Question

나는 DNA 염기 서열에 대해 d까지 해밍 거리와 패턴 매칭 문제를 해결하고 있었다. Regex는 저를 거기에서 구했습니다. 하지만 지금은 다른 문제에 직면했습니다. 긴 DNA 염기 서열을 감안할 때, 나는 최대 d 개의 불일치를 가진 가장 빈번한 불일치 k-mer를 찾아야합니다. 여기서 k-mer는 길이 k의 하위 시퀀스를 의미합니다.

참고 : DNA 염기 서열은 {A, C, G, T}의 네 글자로만 표현할 수 있습니다.

예를 들면

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

여기에서는 두 개의 k-mer 만 가능합니다 : "AGG", "GGC"

이제 "AGG"및 "GGC"에 대해 다음 코드를 실행하여 1 개의 불일치로 개별적으로 순화 할 수 있습니다.

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"는 다음을 제공합니다.

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

그리고 "GCC"는 다음을 제공합니다.

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

그런 다음 Counter가장 빈번한 k-mer를 찾는 데 사용할 수 있습니다 . 그러나 이것은 d = 1. 이것을 일반화하는 방법은 무엇 d <= k입니까?

Answer 1

이것은 계산 비용이 많이 드는 방법입니다. 그러나 네 원하는 것을 가져와야합니다. 내가 여기서 한 것은 hamming dist 1로 모든 불일치를 계산 한 다음 이전 불일치에서 ham dist 1로 새로운 불일치를 계산하고 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)]))

더 나은 성능 코드를 거의 10-20 배 더 빠르게 찾았습니다.

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

CPU 시간 : 사용자 1 분 32 초, 시스템 : 1.81 초, 총 : 1 분 34 초 벽 시간 : 1 분 34 초