Учебная страница курса биоинформатики,
год поступления 2017
Указания по EMBOSS и Standalone BLAST
1. Команды EMBOSS
wossname ищет программу пакета EMBOSS по ключевому слову из её описания.
Если уже знаете название программы, то программа tfm выдаёт её подробное описание. Например:
tfm getorf
выдаст описание программы getorf.
К 13 ноября необходимо будет выучить предназначение программ, перечисленных в начале задания.
2. Скрипт
Есть два варианта: (1) написать скрипт на bash, который вызывает программы и ваши же скрипты на python; (2) написать всё, в том числе вызов программ, на python.
Как вызывать программу из скрипта на python
Есть два способа:
os.system.(проще, но не рекомендован специалистами). Например, вам нужен список описаний тех белков курицы (из SwissProt), которые содержат в своём составе слово "Trypsin". Если такой список вам нужен один раз, то выполняете на kodomo команду:
infoseq sw:*_chick | grep Trypsin > mylist.txt
и в файле myslist.txt получаете искомое. А если слово не обязательно "Trypsin", а организм не обязательно с мнемоникой "chick" и вы хотите написать скрипт, который ищет описания в SwissProt по заданным мнемонике и слову из описания? Тогда, например, так:
import os from sys import argv mnem = argv[1] word = argv[2] command = "infoseq sw:*_" + mnem + " | grep " + word + " > mylist.txt" os.system(command)
Если такому скрипту дать в командной строке мнемонику и слово, он сделает файл result.txt с нужным содержимиым. Если же вам нужен не просто файл, а какие-то результаты обработки данных, то нужно продолжить скрипт: открыть файл "result.txt" и извлекать информацию из его строк.
subprocess.Popen. Скрипт, делающий то же, в этом варианте:
from subprocess import Popen, PIPE from sys import argv mnem = argv[1] word = argv[2] command1 = ["infoseq", "sw:*_" + mnem] process1 = Popen(command1, stdout = PIPE, stderr = PIPE) (out1, err1) = process1.communicate() command2 = ["grep", word] process2 = Popen(command2, stdin = PIPE, stdout = PIPE, stderr = PIPE) (out2, err2) = process2.communicate(input = out1) lines = out2.split('\n')
Здесь наоборот: то, что вам нужно, попало в список lines, и с этой информацией можно дальше работать (например, записать в файл "mylist.txt", тогда результат будет точно таким же, как у первого скрипта).
3. Табличная выдача BLAST
Программы blastp, blastn, blastx, tblastn и tblastx имеют опцию -outfmt, которую можно использовать, если выходной файл нужно не разглядывать, а обрабатывать скриптами. Самые простые из 12 предлагаемых форматов — табличные, то есть -outfmt 6 "tabular" и -outfmt 7 "tabular with comment lines". В этих форматах выдаётся только информация о выравниваниях, ни списка находок, ни самих выравниваний нет. Но для задачи о числе находок с e-value лучше заданного как раз такой формат очень удобен: e-value находится в предпоследнем столбце таблицы.
Не забывайте, что файл query может содержать не одну, а много последовательностей в формате fasta. Для всех этих последовательностей будут выданы находки.
Ещё одна особенность программ blast, которую можно использовать, состоит в том, что если -query явно не указать, то последовательность для поиска берётся со стандартного входа (stdin). Аналогично, если не указать выходной файл (параметр -out), то результат будет выдан на стандартный выход (stdout). Поэтому можно писать конвейеры вида:
seqret file.fasta:name stdout -auto | blastn -db mybase -outfmt 6 | grep somename