рабочий BLAST (bl2seq), не создавая файлы последовательности
У меня есть сценарий, который выполняет запросы BLAST (bl2seq)
Сценарий работает как это:
- Получите последовательность a, упорядочьте b
- запишите последовательность в filea
- запишите последовательность b в fileb
- команда выполнения 'bl2seq-i filea-j fileb-n blastn'
- будьте произведены от STDOUT, синтаксического анализа
- повторите 20 миллионов раз
Программа bl2seq не поддерживает передачу по каналу. Там какой-либо путь состоит в том, чтобы сделать это и постараться не писать/читать в жесткий диск?
Я использую Python BTW.
5 ответов
Откуда вы знаете, что bl2seq не поддерживает pipeing? Кстати, pipes - это особенность ОС, а не программы. Если ваша программа bl2seq что-то выводит, будь то в STDOUT или в файл, вы должны быть в состоянии разобрать вывод. Посмотрите в файле справки программы bl2seq опции для вывода в файл, например, -o. Тогда вы сможете разобрать файл.
Также, поскольку вы используете Python, альтернативой может быть модуль BioPython.
В зависимости от того, на какой ОС вы работаете, вы можете использовать что-то вроде замены процесса bash . Я не уверен, как вы это настроили в Python, но в основном вы используете именованный канал (или именованный дескриптор файла). Это не сработает, если bl2seq попытается выполнить поиск в файлах, но это должно сработать, если он просто читает их последовательно.
Это программа bl2seq из BioPerl? Если да, то не похоже, что в ней можно что-то сделать. Однако вы можете написать свой собственный хак, используя Bio::Tools::Run::AnalysisFactory::Pise, что является рекомендуемым способом. Однако вам придется делать это на Perl.
Если это другой bl2seq, то не обращайте внимания на это сообщение. В любом случае, вам, вероятно, следует предоставить больше деталей.
Вау. Я разобрался.
Ответ - использовать модуль подпроцесса и каналы Python!
РЕДАКТИРОВАТЬ: забыл упомянуть, что я использую blast2, который поддерживает трубопровод.
(это часть класса)
def _query(self):
from subprocess import Popen, PIPE, STDOUT
pipe = Popen([BLAST,
'-p', 'blastn',
'-d', self.database,
'-m', '8'],
stdin=PIPE,
stdout=PIPE)
pipe.stdin.write('%s\n' % self.sequence)
print pipe.communicate()[0]
где self.database - это строка, содержащая имя файла базы данных, т.е. 'nt.fa' self.sequence - это строка, содержащая последовательность запроса
Это выводит результат на экран, но вы можете легко его проанализировать. Нет медленного дискового ввода-вывода. Нет медленного синтаксического анализа XML. Я собираюсь написать для этого модуль и выложить его на github.
Кроме того, я еще не дошел до этого, но думаю, что вы можете выполнять несколько запросов, так что базу данных blast не нужно читать и загружать в оперативную память для каждого запроса.
Я вызываю blast2 с помощью сценария R:
....
system("mkfifo seq1")
system("mkfifo seq2")
system("echo sequence1 > seq1"), wait = FALSE)
system("echo sequence2 > seq2"), wait = FALSE)
system("blast2 -p blastp -i seq1 -j seq2 -m 8", intern = TRUE)
....
Это в 2 раза медленнее (!) По сравнению с записью и чтением с жесткого диска!