在继续阅读之前 ,您是否看过biopython?
似乎您要查找一个替换错误,零插入/删除错误(即汉明距离为1)的近似匹配项。
如果您具有汉明距离匹配功能(例如,参见Ignacio提供的链接),则可以像这样使用它来搜索第一个匹配项:
any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))
但这会相当慢,因为(1)汉明距离函数会在发生第二次替换错误后继续打磨(2)在失败后会继续磨削,它会将光标前进一个,而不是根据所看到的内容向前跳过(例如,博耶-
摩尔搜索)。
您可以使用以下函数克服(1):
def Hamming_check_0_or_1(genome, posn, sequence): errors = 0 for i in xrange(25): if genome[posn+i] != sequence[i]: errors += 1 if errors >= 2: return errors return errors
注意:故意不是Pythonic,而是Cic,因为您需要使用C(可能通过Cython)来获得合理的速度。
Navarro和Raffinot(谷歌“ Navarro Raffinotnrgrep”)已经完成了一些有关跳过比特并行近似Levenshtein搜索的工作,并且可以将其应用于汉明搜索。请注意,位并行方法对查询字符串的长度和字母大小有限制,但是您的分别为25和4,因此在那里没有问题。更新:字母大小为4时,跳过的帮助可能不多。
当您在Google上搜索汉明距离搜索时,您会发现很多有关在硬件中实现它的内容,而在软件中却没有太多。这是一个很大的暗示,也许您想出的任何算法都应使用C或其他某种编译语言来实现。
更新: 位并行方法的工作代码
我还提供了一种简化方法来帮助进行正确性检查,并且打包了Paul
re的变体以进行比较。请注意,使用re.finditer()会产生不重叠的结果,并且这可能导致distance-1匹配遮盖完全匹配;看到我的最后一个测试案例。
位并行方法具有以下特征:保证线性行为O(N),其中N是文本长度。注意朴素的方法是O(NM),正则表达式的方法也是M(图案长度)。Boyer-
Moore型方法将是最坏的情况O(NM)和预期的O(N)。同样,当必须缓冲输入时,可以很容易地使用位并行方法:可以一次输入一个字节或一个兆字节;没有前瞻性,缓冲区边界没有问题。最大的优点是:用C编码时,每个输入字节的代码的速度很简单。
缺点:模式长度实际上限制为快速寄存器中的位数,例如32或64。在这种情况下,模式长度为25;模式长度为25。没问题。它使用的额外内存(S_table)与模式中不同字符的数量成比例。在这种情况下,“字母大小”仅为4;没问题。
本技术报告的详细信息。该算法用于使用Levenshtein距离进行近似搜索。为了转换为汉明距离,我简单地(!)删除了处理插入和删除的语句2.1。您会注意到很多引用带有“
d”上标的“ R”。“ d”是距离。我们只需要0和1。这些“ R”在下面的代码中成为R0和R1变量。
# coding: asciifrom collections import defaultdictimport re_DEBUG = 0# "Fast Text Searching with Errors" by Sun Wu and Udi Manber# TR 91-11, Dept of Computer Science, University of Arizona, June 1991.# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.20.8854def WM_approx_Ham1_search(pattern, text): """Generate (Hamming_dist, start_offset) for matches with distance 0 or 1""" m = len(pattern) S_table = defaultdict(int) for i, c in enumerate(pattern): S_table[c] |= 1 << i R0 = 0 R1 = 0 mask = 1 << (m - 1) for j, c in enumerate(text): S = S_table[c] shR0 = (R0 << 1) | 1 R0 = shR0 & S R1 = ((R1 << 1) | 1) & S | shR0 if _DEBUG: print "j= %2d msk=%s S=%s R0=%s R1=%s" % tuple([j] + map(bitstr, [mask, S, R0, R1])) if R0 & mask: # exact match yield 0, j - m + 1 elif R1 & mask: # match with one substitution yield 1, j - m + 1if _DEBUG: def bitstr(num, mlen=8): wstr = "" for i in xrange(mlen): if num & 1: wstr = "1" + wstr else: wstr = "0" + wstr num >>= 1 return wstrdef Ham_dist(s1, s2): """Calculate Hamming distance between 2 sequences.""" assert len(s1) == len(s2) return sum(c1 != c2 for c1, c2 in zip(s1, s2))def long_check(pattern, text): """Naively and understandably generate (Hamming_dist, start_offset) for matches with distance 0 or 1""" m = len(pattern) for i in xrange(len(text) - m + 1): d = Ham_dist(pattern, text[i:i+m]) if d < 2: yield d, idef Paul_McGuire_regex(pattern, text): searchSeqREStr = ( '(' + pattern + ')|(' + ')|('.join( pattern[:i] + "[ACTGN]".replace(c,'') + pattern[i+1:] for i,c in enumerate(pattern) ) + ')' ) searchSeqRE = re.compile(searchSeqREStr) for match in searchSeqRE.finditer(text): locn = match.start() dist = int(bool(match.lastindex - 1)) yield dist, locnif __name__ == "__main__": genome1 = "TTTACGTAAACTAAACTGTAA" # 01234567890123456789012345 # 1 2 tests = [ (genome1, "ACGT ATGT ACTA ATCG TTTT ATTA TTTA"), ("T" * 10, "TTTT"), ("ACGTCGTAAAA", "TCGT"), # partial match can shadow an exact match ] nfailed = 0 for genome, patterns in tests: print "genome:", genome for pattern in patterns.split(): print pattern a1 = list(WM_approx_Ham1_search(pattern, genome)) a2 = list(long_check(pattern, genome)) a3 = list(Paul_McGuire_regex(pattern, genome)) print a1 print a2 print a3 print a1 == a2, a2 == a3 nfailed += (a1 != a2 or a2 != a3) print "***", nfailed


