สร้าง Permutations ทั้งหมดโดยไม่ตรงกันมากที่สุด d

Question

ฉันกำลังแก้ปัญหาสำหรับการจับคู่รูปแบบที่มีระยะการตอกถึง d สำหรับลำดับดีเอ็นเอ Regex ช่วยฉันไว้ที่นั่น แต่ตอนนี้ฉันพบปัญหาที่แตกต่างออกไป จากลำดับดีเอ็นเอที่ยาวฉันพบว่า k-mer ที่ไม่ตรงกันบ่อยที่สุดโดยที่ d ไม่ตรงกันเกือบทั้งหมด ในที่นี้ k-mer หมายถึงลำดับย่อยของความยาว k

หมายเหตุ:ลำดับดีเอ็นเอสามารถแสดงได้โดยใช้ตัวอักษรสี่ตัวเท่านั้น: {A, C, G, T}

ตัวอย่างเช่น,

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

ที่นี่สามารถทำได้เพียงสอง k-mers: "AGG", "GGC"

ตอนนี้ฉันสามารถอนุญาตทีละรายการโดยมี 1 รายการที่ไม่ตรงกันโดยเรียกใช้รหัสต่อไปนี้สำหรับ "AGG" และ "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" จะให้:

['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 วินาที