有限制的散列序列

Question

我必須修改以下代碼以查找允許每三個鹼基（一個核苷酸）不匹配的所有100-mers（100個核苷酸的塊）。任何關於如何解決這個問題的邏輯將不勝感激。謝謝！

# length of hash key 
kmerlen = 30 

# hash table for finding hits 
lookup = defaultdict(list) 

# store sequence hashes in hash table 
print("hashing seq1...") 
for i in xrange(len(seq1) - kmerlen + 1): 
    key = seq1[i:i+kmerlen] 
    lookup[key].append(i) 

# look up hashes in hash table 
print("hashing seq2...") 
hits = [] 
for i in xrange(len(seq2) - kmerlen + 1): 
    key = seq2[i:i+kmerlen] 

    # store hits to hits list 
    for hit in lookup.get(key, []): 
     hits.append((i, hit)) 

# hits should be a list of tuples 
# [(index1_in_seq2, index1_in_seq1), 
# (index2_in_seq2, index2_in_seq1), 
# ...] 

Answer 1

我想你只需要一個step項添加到您的切片表達式：

kmerlen = 30 
kmerstep = 1  # new variable 

# ... 
for i in xrange(len(seq1) - kmerlen + 1): 
    key = seq1[i:i+kmerlen:kmerstep] # add step to slice 
    lookup[key].append(i) 

# ... 
for i in xrange(len(seq2) - kmerlen + 1): 
    key = seq2[i:i+kmerlen:kmerstep] # here too 

這會爲你的工作2-4工作，具有調節步長。任務5有點棘手，因爲每三個項目需要兩個匹配。我建議連接兩個切片爲：

key = seq1[i:i+kmerlen:3]+seq1[i+1:i+kmerlen:3] 

的項目是不是爲了，但如果你切其他序列以同樣的方式，他們應該完全對應。您也可能需要調整迴路range以獲得i。