生成最多d个不匹配的所有排列

如何解决生成最多d个不匹配的所有排列

我正在解决DNA序列汉明距离高达d的模式匹配问题。正则表达式将我保存在那里。但是现在我遇到了另一个问题。考虑到较长的DNA序列，我必须找到最频繁错配的k-mer，最多错配d个。这里，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','GTC','TGC','CGC','GGC']

然后，我可以使用Counter查找最频繁的k-mer。但是，这仅在d = 1下有效。如何对任何d <= k进行概括？

解决方法

这是计算上昂贵的方法。但是，是的，应该能够满足您的要求。我在这里所做的是计算与hamming dist 1的所有不匹配，然后从先前的不匹配和递归直到d来计算与ham dist 1的新不匹配。

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),c in enumerate(sub)])))))
    
    setOfMmatch=setOfMmatch.union(newMmatches)
    
    if not mmatchHash.get((d-1,str(setOfMmatch)),0):
        mmatchHash.set((d-1,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时间：用户1min 32s，sys：1.81 s，总计：1min 34s 挂墙时间：1分34秒