2

我想找到 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

你能打我吗，为什么？我不希望您发布正确的代码，而是给我一个想法，我的错误在哪里（我认为错误可能在第二个函数中）。

PS 我确实意识到我必须在 Stepic 在线课程上解决以下任务，但是没有来自学习组的在线社团的反馈，我已经在 StackOverflow 上发布了我的代码。

4

2 回答 2

1

Aaaaand，这是一个基于 re-based 的解决方案：

genome = "CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGCGGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCCGGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCGCCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTACCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACACACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGGGCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"
length=10
hamming_distance=2


import re
from itertools import product

charset = "ACGT"

# first walk genome, get all unique N-character patterns
patterns = set()
for i in range(len(genome)-length):
    patterns.add(genome[i:i+length])


# for each pattern, create re's with all possible alternates
allmatches = {}
for p in sorted(patterns):
    options = [[c,"["+charset+"]"] for c in p]
    proditer = product(*options)

    matches = set()
    for pr in proditer:
        # count how many elements in this product start with '['
        # only want up to hamming_distance number of alternatives
        numalts = sum(prpart.startswith('[') for prpart in pr)
        if numalts > hamming_distance:
            continue

        compiled_pattRE = re.compile(''.join(pr))
        for match in filter(None, (compiled_pattRE.match(genome,i) 
                                        for i in range(len(genome)-length+1))):
            matches.add((match.start(), match.group(0)))

    allmatches[p] = matches

for p,matches in sorted(allmatches.items(), key=lambda m: -len(m[1])):
    print p, len(matches)
    for m in sorted(matches):
        print m
    print
    break

印刷：

GCGCACACAC 20
(12, 'CGGCACACAC')
(14, 'GCACACACAG')
(126, 'GGGTACACAC')
(138, 'GGGCGCACAC')
(140, 'GCGCACACAC')
(142, 'GCACACACAG')
(213, 'CGGCACACAC')
(215, 'GCACACACAG')
(227, 'GCCCACACAC')
(253, 'GCGCACACAC')
(255, 'GCACACACAC')
(257, 'ACACACACAC')
(272, 'GCGCACAGCC')
(280, 'CCGCCCACAC')
(282, 'GCCCACACAC')
(284, 'CCACACACAC')
(316, 'GGGCGCACAC')
(318, 'GCGCACACAC')
(320, 'GCACACACAC')
(351, 'GCGCACAGCC')
于 2013-11-05T02:45:23.670 回答
1

我在 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
         ^^
于 2013-11-04T20:29:02.837 回答