Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2017

Указания по EMBOSS и Standalone BLAST

1. Команды EMBOSS

wossname ищет программу пакета EMBOSS по ключевому слову из её описания.

Если уже знаете название программы, то программа tfm выдаёт её подробное описание. Например:

tfm getorf

выдаст описание программы getorf.

К 13 ноября необходимо будет выучить предназначение программ, перечисленных в начале задания.

2. Скрипт

Есть два варианта: (1) написать скрипт на bash, который вызывает программы и ваши же скрипты на python; (2) написать всё, в том числе вызов программ, на python.

Как вызывать программу из скрипта на python

Есть два способа:

  1. 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" и извлекать информацию из его строк.

  1. 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