Практикум 7. EMBOSS, BLAST+
Этот практикум посвящён знакомству с пакетами EMBOSS и BLAST+ в применении к работе с нуклеотидными последовательностями.
Упражнения по пакету программ EMBOSS
Для семи из предложенных заданий были написаны команды EMBOSS, выполняющие их. Файл с командами можно просмотреть здесь (он также доступен по пути /home/students/y19/pork7007/term3/block2/pr9_emboss.txt).
Написание скрипта, подсчитывающего число случайных находок с заданным E-value
Этот скрипт выполняет первую из трёх задач, предложенных для практикума:
1. Проверить, сколько находок с E-value < 0.1 в среднем находит blastn для случайной последовательности данной длины в данном геноме бактерии.
Скачать скрипт можно по этой ссылке. Он написан на Python с использованием утилит EMBOSS и Standalone BLAST и получился достаточно громоздким исключительно благодаря моему желанию предусмотреть ввод неправильных аргументов. Скрипт работает следующим образом: аннотирует заданную длинную последовательность (бактериальный геном) как базу данных для поиска с помощью BLAST, затем некоторое количество раз генерирует случайную последовательность заданной длины, ищет её в базе и записывает число находок с E-value меньше некоторого порога, после чего подсчитывает среднее число таких находок среди всех случайных последовательностей.
Скрипт принимает на ввод следующие аргументы командной строки:
• длина случайных последовательностей
• файл или USA-адрес записи с последовательностью, выполняющей роль базы
• число генерируемых последовательностей (опционально, по умолчанию 100)
• порог на E-value (опционально, по умолчанию 0.1)
На выходе получается число (округлённое до 4 знаков после запятой) - среднее количество находок. Теоретически, оно должно быть близко по своему значению к установленному порогу на E-value (согласно определению E-value).
Ниже приведён пример работы скрипта:
pork7007@kodomo:~/term3/block2/pr7$ ./evaluecalc.py 100 seq.fasta 0.12
Здесь 100 - длина случайных последовательностей, seq.fasta - файл с последовательностью хромосомы цианобактерии Nostoc punctiforme (RefSeq AC: NC_010628), находящийся в текущей директории. При таких входных данных скрипт был запущен 3 раза, на выходе были получены значения 0.09, 0.08 и 0.12.
Просмотреть код
#! /usr/local/bin/python3
import os, subprocess
from sys import argv
# далее идёт обработка исключений, возникающих при некорректных аргументах (просто для красоты, чтобы скрипт не падал от отсутствия чего-то)
is_ok = True
try:
ran_seq_len = int(argv[1])
except IndexError:
is_ok = False
print('No argument: please input length of random sequences')
except ValueError:
is_ok = False
print('Invalid argument: please input length of random sequences')
try:
os.system('seqret {} subj.evalcalc -filter 2> /dev/null'.format(argv[2])) # пытаемся вытащить последовательность из файла / по USA-адресу
os.system('makeblastdb -in subj.evalcalc -out db.evalcalc -dbtype nucl -logfile /dev/null') # и сделать базу для blast
if not (os.path.isfile('db.evalcalc.nhr') and os.path.isfile('db.evalcalc.nin') and os.path.isfile('db.evalcalc.nsq')):
is_ok = False # отмечу, что все создаваемые в ходе работы временные файлы, для которых
print('Invalid argument: please input file or database entry to be used as subject') # это возможно, имеют расширение .evalcalc - чтобы случайно не
except IndexError: # перезаписать/удалить файл пользователя
is_ok = False
print('No argument: please input file or database entry to be used as subject')
try:
ran_seq_num = int(argv[3])
except IndexError:
ran_seq_num = 100
except ValueError:
is_ok = False
print('Invalid argument: please input number of random sequences')
try:
evalue_threshold = float(argv[4])
except IndexError:
evalue_threshold = 0.1
except ValueError:
is_ok = False
print('Invalid argument: please input E-value threshold')
# здесь начинается собственно подсчёт среднего числа находок с E-value не выше evalue_threshold
if is_ok:
res_num_list = []
res_sum = 0
for _ in range(ran_seq_num):
os.system('makenucseq -amount 1 -length {} -outseq query.evalcalc -filter'.format(ran_seq_len)) # генерируем случайную последовательность
os.system('blastn -task blastn -query query.evalcalc -db db.evalcalc -out blastout.evalcalc -word_size 4 -evalue {} -outfmt "6 evalue"'.format(evalue_threshold)) # ищем по базе, записываем в файл находки с нужным E-value
num = int(subprocess.check_output(['wc', '-l', 'blastout.evalcalc']).decode().split(' ')[0]) # считаем число находок
res_num_list.append(num) # и записываем в лист; повторяем цикл
for i in res_num_list:
res_sum += i
print(round(res_sum / len(res_num_list), 4)) # выводим результат - среднее число находок
os.system('rm db.evalcalc.*')
os.system('rm *.evalcalc')
Поиск гомологов белков в неаннотированном геноме
Для этой цели была использована предложенная сборка генома Amoeboaphelidium protococcarum, а также выбраны белки для поиска гомологов: метионин-тРНК лигаза (UniProt AC: P00958), бета-цепь тубулина (UniProt AC: P02557) и бета-субъединица АТФ-синтазы (UniProt AC: P00830) из Saccharomyces cerevisiae (strain ATCC 204508 / S288c). Поиск в UniProt проводился по следующим запросам:
"methionine trna ligase" organism:"saccharomyces cerevisiae"
"tubulin" organism:"saccharomyces cerevisiae"
"atp synthase" organism:"saccharomyces cerevisiae"
Последовательности белков были скачаны на kodomo при помощи seqret:
seqret sw:P00958 SYMC_YEAST.fasta
seqret sw:P02557 TBB_YEAST.fasta
seqret sw:P00830 ATPB_YEAST.fasta
Далее приведены ссылки на последовательности: SYMC_YEAST.fasta, TBB_YEAST.fasta, ATPB_YEAST.fasta
Сборка генома была проиндексирована как база данных BLAST:
makeblastdb -in "X5.fasta" -dbtype "nucl"
Затем для каждой последовательности белка был проведён поиск по индексированному геному при помощи tblastn:
tblastn -query "SYMC_YEAST.fasta" -db "X5.fasta" -out "SYMC_YEAST_search.txt" -outfmt "7 std qcovs"
tblastn -query "TBB_YEAST.fasta" -db "X5.fasta" -out "TBB_YEAST_search.txt" -outfmt "7 std qcovs"
tblastn -query "ATPB_YEAST.fasta" -db "X5.fasta" -out "ATPB_YEAST_search.txt" -outfmt "7 std qcovs"
Выдачу BLAST для каждого запуска можно просмотреть по ссылкам: SYMC_YEAST, TBB_YEAST, ATPB_YEAST
Легко заметить, что для каждого белка BLAST нашёл в геноме Amoeboaphelidium protococcarum последовательности, идентичные более чем на 50% (в случае TBB_YEAST и ATPB_YEAST - более чем на 80%), с покрытием запроса не менее 75% и E-value, близким к 0; это позволяет говорить о том, что у Amoeboaphelidium protococcarum присутствуют гомологи всех трёх рассматриваемых белков.