我在谷歌上搜索“弓形虫基因组”，以在线查找一些基因组文件。我找到了一个看起来很接近的文件，标题为“TgondiiGenomic_ToxoDB-6.0.fasta”，大小约为158MB，位于http://toxodb.org。我使用以下pyparsing表达式提取基因序列，仅用了不到2分钟：
fname = "TgondiiGenomic_ToxoDB-6.0.fasta"
fastasrc = open(fname).read()
"""
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"))
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步骤中提取基因序列进行匹配后，这个基于RE的搜索器花了大约1-1/2分钟的时间扫描所有未经文本包装的序列，找到了与pyparsing解决方案相同的所有1个不匹配的条目。
shR0中的sh代表什么？ - RandomblueWM_approx_Ham1_search函数是否能够在文本长度小于100000，模式串长度小于文本长度且字母表基数为27的情况下顺利运行？ - Randomblue