Я гуглил «геном паразита токсоплазмы гонди», чтобы найти некоторые из этих файлов генома онлайн. Я нашел то, что, по моему мнению, было близко, файл с названием «TgondiiGenomic_ToxoDB-6.0.fasta» размером http://toxodb.org, размером около 158 МБ. Я использовал следующее выражение 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.)
Затем я запустил этот код для поиска совпадений «до-2-несоответствия» во всех последовательностях, считанных из файла .fasta (напомним, что genedata - это последовательность групп ParseResults, каждая из которых содержит идентификатор, целое число длина и строка последовательности):
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 минут у меня был этот вывод, перечисляющий каждый идентификатор и длину гена, и любые найденные частичные совпадения:
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 / 2 минуты, чтобы отсканировать все нетекстовые последовательности, чтобы найти все одинаковые 1 -согласовать записи, которые сделал решение pyparsing.