これらのゲノムファイルのいくつかをオンラインで見つけるために、「toxoplasmagondii寄生虫ゲノム」をグーグルで検索しました。http://toxodb.orgにある「TgondiiGenomic_ToxoDB-6.0.fasta」というタイトルのファイルで、サイズが約158Mbで、近いと思うものを見つけました。次の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)
一致するたびに、一致した実際の文字列と不一致の場所のリストを含むタプルが返されます。もちろん、完全一致は2番目の値の空のリストを返します。（私はこのクラスが好きです。次のリリースの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つからランダムに検索シーケンスを取得して、完全に一致するものを見つけられるようにしました。好奇心から、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つの一致も見つけました。
比較のために、これがREベースのバージョンで、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分半かかりました。 pyparsingソリューションが行ったこと。