我想找到所有的模式 1) 字符串中出现次数最多的 2) 最多有 d-不匹配。
对于这个给定的任务，我实现了一个函数，该函数计算给定模式在 d 不匹配的字符串中出现的次数。该算法的思想基于使用字符串子模式的位掩码和给定模式的位掩码的卷积。它产生正确的结果。这是该算法的代码:
def create_bit_mask(letter, text):
buf_array=[]
for c in text:
if c==letter:
buf_array.append(1)
else:
buf_array.append(0)
return buf_array
def convolution(bit_mask1, bit_mask2):
return sum(b*q for b,q in zip(bit_mask1, bit_mask2))
def number_of_occurances_with_at_most_hamming_distance(genome,pattern,hamming_distance):
alphabet=["A","C","G","T"]
matches=0
matrix_of_bit_arrays_for_pattern=[]
matrix_of_bit_arrays_for_genome=[]
buf_output=0
buf=0
positions=""
for symbol in alphabet:
matrix_of_bit_arrays_for_pattern.append(create_bit_mask(symbol,pattern))
matrix_of_bit_arrays_for_genome.append(create_bit_mask(symbol, genome))
for i in xrange(len(genome)-len(pattern)+1):
buf_debug=[]
buf=sum(convolution(bit_mask_pattern,bit_mask_genome[i:i+len(pattern)]) for bit_mask_pattern, bit_mask_genome in zip(matrix_of_bit_arrays_for_pattern,matrix_of_bit_arrays_for_genome))
hamming=len(pattern)-buf
if hamming<=hamming_distance:
buf_output+=1
#print "current window: ", genome[i:i+len(pattern)], "pattern :", pattern,"number of mismatches : ", hamming, " @ position : ",i
return buf_output
鉴于上面的功能，任务的解决方案应该相当简单，即
def fuzzy_frequency():
genome="CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGCGGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCCGGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCGCCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTACCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACACACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGGGCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"
length=10
hamming_distance=2
for i in range(len(genome)-length):
#print genome[i:i+10], " number of occurances: ", number_of_occurances_with_at_most_hamming_distance(genome, genome[i:i+10],hamming_distance)
print (genome[i:i+length],number_of_occurances_with_at_most_hamming_distance(genome, genome[i:i+length],hamming_distance))
但是，上面的代码没有找到以下子字符串:
GCACACAGAC
你能帮我打一下吗，为什么？我不想让你发布正确的代码，而是给我一个想法，我的错误在哪里(我认为错误可能在第二个函数中)。
附注我确实意识到我必须在 Stepic 在线类(class)上解决以下任务，但由于没有来自学习小组的在线社区的反馈，我已在 StackOverflow 上发布了我的代码。
最佳答案
我在 pyparsing 列表上提出了类似的基因组解析问题，我想出了这个 CloseMatch 解析器类。它将大部分字符串遍历和测试代码封装在 pyparsing 自己的字符串解析框架中，但这仍然可以让您深入了解自己的代码:
genome = "CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGCGGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCCGGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCGCCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTACCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACACACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGGGCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"
length=10
hamming_distance=2
from pyparsing import Token, ParseException
# following from pyparsing.wikispaces.com Examples page
class CloseMatch(Token):
"""A special subclass of Token that does *close* matches. For each
close match of the given string, a tuple is returned giving the
found close match, and a list of mismatch positions."""
def __init__(self, seq, maxMismatches=1):
super(CloseMatch,self).__init__()
self.name = seq
self.sequence = seq
self.maxMismatches = maxMismatches
self.errmsg = "Expected " + self.sequence
self.mayIndexError = False
self.mayReturnEmpty = False
def parseImpl( self, instring, loc, doActions=True ):
start = loc
instrlen = len(instring)
maxloc = start + len(self.sequence)
if maxloc <= instrlen:
seq = self.sequence
seqloc = 0
mismatches = []
throwException = False
done = False
while loc < maxloc and not done:
if instring[loc] != seq[seqloc]:
mismatches.append(seqloc)
if len(mismatches) > self.maxMismatches:
throwException = True
done = True
loc += 1
seqloc += 1
else:
throwException = True
if throwException:
#~ exc = self.myException
#~ exc.loc = loc
#~ exc.pstr = instring
#~ raise exc
raise ParseException(instring, loc, self.errmsg)
return loc, (instring[start:loc],mismatches)
# first walk genome, get all unique N-character patterns
patterns = set()
for i in range(len(genome)-length):
patterns.add(genome[i:i+length])
print len(patterns)
# use pyparsing's CloseMatch to find close matches - each match
# returns the substring and the list of mismatch locations
matches = {}
for p in sorted(patterns):
matcher = CloseMatch(p, hamming_distance)
matches[p] = list(matcher.scanString(genome, overlap=True))
# Now list out all patterns and number of close matches - for the most
# commonly matched pattern, dump out all matches, where they occurred and
# an annotated match showing the mismatch locations
first = True
for p in sorted(matches, key=lambda m: -len(matches[m])):
if first:
first = False
for matchdata in matches[p]:
matchvalue, start, end = matchdata
substring,mismatches = matchvalue[0]
print ' ', substring, 'at', start
if mismatches:
print ' ', ''.join('^' if i in mismatches else ' ' for i in range(length))
else:
print ' ', "***EXACT***"
print
print p, len(matches[p])
给定基因组中有 254 个独特的 10 个字符模式。
以下是最常匹配模式的输出:
CGGCACACAC at 12
^^
GCACACACAG at 14
^ ^
GGGTACACAC at 126
^ ^
GGGCGCACAC at 138
^ ^
GCGCACACAC at 140
***EXACT***
GCACACACAG at 142
^ ^
CGGCACACAC at 213
^^
GCACACACAG at 215
^ ^
GCCCACACAC at 227
^
GCGCACACAC at 253
***EXACT***
GCACACACAC at 255
^
ACACACACAC at 257
^ ^
GCGCACAGCC at 272
^^
CCGCCCACAC at 280
^ ^
GCCCACACAC at 282
^
CCACACACAC at 284
^ ^
GGGCGCACAC at 316
^ ^
GCGCACACAC at 318
***EXACT***
GCACACACAC at 320
^
GCGCACAGCC at 351
^^
关于python - 字符串中的模糊模式搜索 : the most frequent patterns with d-mismatches，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/19775978/