问题描述
，您是否看过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: ascii
from collections import defaultdict
import 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.8854
def 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 + 1
if _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 wstr
def 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, i
def 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, locn
if __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
解决方法
我正在处理长度为25的DNA序列（请参见下面的示例）。我有一个230,000的列表，需要查找整个基因组中的每个序列（弓形虫寄生虫）。我不确定基因组有多大，但是比230,000个序列长得多。
我需要查找25个字符的每个序列，例如（AGCCTCCCATGATTGAACAGATCAT）。
基因组被格式化为连续字符串，即（CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT....）
我不在乎它在哪里或被发现多少次，无论它是否存在。
我认为这很简单-
str.find(AGCCTCCCATGATTGAACAGATCAT)
但是我也要找到一个在任何位置（只有一个位置）被定义为错误（不匹配）的紧密匹配项，并在序列中记录该位置。我不确定该怎么做。我唯一想到的就是使用通配符，并在每个位置使用通配符执行搜索。即搜索25次。
例如，
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT
位置13不匹配的近距离比赛。
速度不是一个大问题，因为我只做3次，不过如果速度快的话会很好。
有执行此操作的程序-查找匹配项和部分匹配项-但我正在寻找一种在这些应用程序中无法发现的部分匹配项类型。