允许字符串中任意位置一个不匹配的字符串搜索
我正在处理长度为25的DNA序列（见下面的例子）。我有一个包含230,000个序列的列表，需要在整个基因组中查找每个序列（毒浆虫基因组）。我不确定基因组有多大，但肯定比230,000个序列要长得多。
我需要查找每个25个字符的序列，比如说，(AGCCTCCCATGATTGAACAGATCAT)。
基因组的格式是一个连续的字符串，也就是像这样：(CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT....)
我不在乎它出现的位置或次数，只想知道它是否存在。
我觉得这很简单 --
str.find(AGCCTCCCATGATTGAACAGATCAT)
但我还想找到一个近似匹配，定义为在任何位置有一个错误（不匹配），并记录在序列中的位置。我不太确定该怎么做。我能想到的唯一方法是使用通配符，在每个位置进行搜索，也就是说，要搜索25次。
例如，
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT
一个在第13个位置有不匹配的近似匹配。
速度不是大问题，因为我只需要做3次，虽然如果能快一点那就更好了。
有一些程序可以做到这一点——找到匹配和部分匹配——但我在寻找一种特定类型的部分匹配，这些应用程序无法发现。
这里有一个类似的帖子，关于perl，虽然他们只是比较序列，而不是在一个连续字符串中搜索：
14 个回答
我在网上搜索“弓形虫基因组”，想找到一些相关的基因组文件。我找到了一份我认为比较接近的文件，名字叫“TgondiiGenomic_ToxoDB-6.0.fasta”，大小大约是158MB，地址在http://toxodb.org。我用下面的pyparsing表达式提取基因序列，花了不到2分钟：
fname = "TgondiiGenomic_ToxoDB-6.0.fasta"
fastasrc = open(fname).read() # yes! just read the whole dang 158Mb!
"""
Sample header:
>gb|scf_1104442823584 | organism=Toxoplasma_gondii_VEG | version=2008-07-23 | length=1448
"""
integer = Word(nums).setParseAction(lambda t:int(t[0]))
genebit = Group(">gb|" + Word(printables)("id") + SkipTo("length=") +
"length=" + integer("genelen") + LineEnd() +
Combine(OneOrMore(Word("ACGTN")),adjacent=False)("gene"))
# read gene data from .fasta file - takes just under a couple of minutes
genedata = OneOrMore(genebit).parseString(fastasrc)
（惊讶！有些基因序列中竟然有一串'N'！这到底是怎么回事？！）
然后我写了一个类，作为pyparsing的Token类的子类，用于进行近似匹配：
class CloseMatch(Token):
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
return loc, (instring[start:loc],mismatches)
每次匹配时，这个类会返回一个元组，里面包含实际匹配的字符串和一个不匹配位置的列表。完全匹配的情况下，第二个值会返回一个空列表。（我很喜欢这个类，觉得可以把它加到下一个pyparsing的版本里。）
接着，我运行了这段代码，搜索从.fasta文件中读取的所有序列中的“最多2个不匹配”的匹配（记得genedata是一个ParseResults组的序列，每个组包含一个id、一个整数长度和一个序列字符串）：
searchseq = CloseMatch("ATCATCGAATGGAATCTAATGGAAT", 2)
for g in genedata:
print "%s (%d)" % (g.id, g.genelen)
print "-"*24
for t,startLoc,endLoc in searchseq.scanString(g.gene):
matched, mismatches = t[0]
print "MATCH:", searchseq.sequence
print "FOUND:", matched
if mismatches:
print " ", ''.join(' ' if i not in mismatches else '*'
for i,c in enumerate(searchseq.sequence))
else:
print "<exact match>"
print "at location", startLoc
print
print
我随机从某个基因片段中选了一个搜索序列，以确保能找到一个完全匹配的，同时也好奇有多少个1个和2个元素的不匹配。
这段代码运行花了一点时间。经过45分钟，我得到了这个输出，列出了每个id和基因长度，以及找到的任何部分匹配：
scf_1104442825154 (964)
------------------------
scf_1104442822828 (942)
------------------------
scf_1104442824510 (987)
------------------------
scf_1104442823180 (1065)
------------------------
...
我有点沮丧，因为一直没有看到匹配，直到：
scf_1104442823952 (1188)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAACGGAATCGAATGGAAT
* *
at location 33
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 175
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 474
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
*
at location 617
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATAGAAT
* *
at location 718
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGATTCGAATGGAAT
* *
at location 896
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGTAT
* *
at location 945
最后我找到了我的完全匹配在：
scf_1104442823584 (1448)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGACTCGAATGGAAT
* *
at location 177
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCAAATGGAAT
*
at location 203
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
* *
at location 350
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAA
* *
at location 523
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
* *
at location 822
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCTAATGGAAT
<exact match>
at location 848
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCGTCGAATGGAGTCTAATGGAAT
* *
at location 969
虽然这没有创造任何速度记录，但我完成了任务，还找到了几个2个不匹配的情况，可能会有用。
为了比较，这里有一个基于正则表达式的版本，只找到1个不匹配的匹配：
import re
seqStr = "ATCATCGAATGGAATCTAATGGAAT"
searchSeqREStr = seqStr + '|' + \
'|'.join(seqStr[:i]+"[ACTGN]".replace(c,'') +seqStr[i+1:]
for i,c in enumerate(seqStr))
searchSeqRE = re.compile(searchSeqREStr)
for g in genedata:
print "%s (%d)" % (g.id, g.genelen)
print "-"*24
for match in searchSeqRE.finditer(g.gene):
print "MATCH:", seqStr
print "FOUND:", match.group(0)
print "at location", match.start()
print
print
（起初，我试着直接搜索原始的FASTA文件，但对比pyparsing版本，发现匹配的数量少得多，感到困惑。后来我意识到，有些匹配可能跨越了换行符，因为fasta文件的输出在n个字符处换行。）
所以在第一次用pyparsing提取基因序列以便匹配之后，这个基于正则表达式的搜索器又花了大约1分半钟扫描所有未换行的序列，找到了与pyparsing解决方案相同的所有1个不匹配的条目。
在继续阅读之前，你有没有查看过biopython？
看起来你想找到大约匹配的字符串，允许一个字符的替换错误，而不允许插入或删除错误，也就是汉明距离为1。
如果你有一个计算汉明距离的匹配函数（可以参考Ignacio提供的链接），你可以像这样使用它来搜索第一个匹配：
any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))
不过这样做会比较慢，因为（1）汉明距离函数在第二个替换错误后仍会继续计算；（2）在失败后，它只会向前移动一个字符，而不是根据之前的结果跳过一些字符（就像Boyer-Moore搜索那样）。
你可以用这样的函数来解决（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
注意：这故意不是Python风格的，而是C风格，因为你需要用C（可能通过Cython）来获得合理的速度。
Navarro和Raffinot做了一些关于位并行近似Levenshtein搜索的跳过方法的工作（可以在网上搜索“Navarro Raffinot nrgrep”），这可以适应汉明搜索。需要注意的是，位并行方法在查询字符串的长度和字母表大小上有一些限制，但你的情况是25和4，所以没有问题。更新：对于字母表大小为4的情况，跳过可能帮助不大。
当你在网上搜索汉明距离搜索时，你会发现很多关于在硬件中实现它的内容，而软件方面的资料不多。这是一个很大的提示，说明你可能需要在C或其他编译语言中实现你想出的算法。
更新：位并行方法的工作代码
我还提供了一种简单的方法来帮助检查正确性，并且我打包了保罗的re代码的一个变体以进行比较。注意，使用re.finditer()会返回不重叠的结果，这可能导致距离为1的匹配覆盖精确匹配；请查看我的最后一个测试案例。
位并行方法具有以下特点：保证线性行为O(N)，其中N是文本长度。注意，简单方法是O(NM)，正则表达式方法也是如此（M是模式长度）。Boyer-Moore风格的方法在最坏情况下是O(NM)，预期是O(N)。此外，当输入需要缓冲时，位并行方法也很容易使用：它可以一次处理一个字节或一个兆字节；没有前瞻，也没有缓冲边界的问题。最大的优势是：当用C编写时，这种简单的每字节输入代码速度非常快。
缺点是：模式长度实际上受限于快速寄存器中的位数，例如32或64。在这种情况下，模式长度是25；没有问题。它使用额外的内存（S_table），与模式中不同字符的数量成正比。在这种情况下，“字母表大小”只有4；没有问题。
详细信息请参考这份技术报告。那里的算法是用于近似搜索的Levenshtein距离。要转换为使用汉明距离，我只是简单地去掉了处理插入和删除的第2.1条中的部分内容。你会注意到很多提到“R”的地方，带有“d”的上标。“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