2023
2022
2021
2020
2018
2017
2016
2015
Задания по материалу дополнительной лекции
1. Исследование зависимости списка находок от параметров BLAST
Проведите исследование зависимости числа "хороших" находок (пусть это будут находки с E-value < 0,001) от параметров. Можно исследовать параметры:
- -word_size
- -threshold
- -window_size
- -xdrop_ungap
- -xdrop_gap и -xdrop_gap_final (кстати, не могу нигде найти точный смысл этих двух параметров, точнее, в чём разница между ними! С.А.С.)
- -comp_based_stats
- -matrix
- -gapopen и -gapextend
Примерный дизайн исследования:
- Выберите какой-нибудь большой и полный протеом (например, мыши) и скачайте его из Uniprot в fasta-формате
- Проиндексируйте его под BLAST программой makeblastdb
Выберите какой-нибудь протеом поменьше (например, Haemophilus influenzae или какой-нибудь микоплазмы) и скачайте его тоже
- Запустите blastp по созданному банку, указав в качестве query маленький протеом (да, BLAST может принимать в качестве запроса много последовательностей сразу и иногда это очень удобно!), с порогом на E-value в 0.001 и с выходным форматом (-outfmt) 6 или 7 (табличным). Для импорта в электронные таблицы удобнее 6, а для обработки скриптом более или менее всё равно.
- Определите список белков из маленького протеома, для которого нашлись хорошие гомологи в большом.
- Теперь повторяйте два предыдущих шага, варьируя выбранный вами параметр, и обработайте результаты.
Неплохо бы также фиксировать время работы blastp в зависимости от значения параметра (в bash есть команда time, которая позволяет определять время работы любой программы: просто в командной строке перед своей командой напишите "time ")
2. Исследование корректности вычисления E-value
Проиндексируйте под BLAST какой-нибудь протеом. Теперь запустите по нему blastp много (минимум 100, а лучше 1000) раз, давая в качестве запроса случайные последовательности. Для разных порогов на E-value (скажем, 0.25, 0.5, 1, 2, 4 и 8) посчитайте среднее число находок с E-value ниже этого порога. Сравните с самими порогами и сделайте выводы.
Смысл исследования в том, что принятый в BLAST метод вычисления E-value основан на предположении, что банк тоже состоит из случайных последовательностей с независимым появлением букв (иногда такие последовательности называются бернуллиевскими). Между тем реальные банки не бернуллиевские и к тому же обладают внутренней структурой (гомологичные белки). Вопрос: насколько это искажает реальное матожидание числа случайных находок?
Советы:
- Индексация банка: программа makeblastdb (запустите с опцией -halp, чтобы увидеть, как задавать параметры). Название банка (-out) удобнее сделать коротким.
- Для получения случайных последовательностей лучше всего взять белок средней величины и программой shuffleseq перемешать его последовательность нужное (параметр -shuffle) число раз, а потом разнести по отдельным файлам программой seqretsplit.
Чтобы запустить blastp много раз, можно написать bash-скрипт (если умеете), а можно вызвать blastp в цикле из Python, функция os.system выполняет свой аргумент (строку) в системе. Есть проблема: os.system запускает параллельный процесс, поэтому нельзя следующей строкой открывать выходной файл, он ещё не будет готов! Нужно сделать все выходные файлы, потом все обработать. Эту проблему можно обойти и другим способом: освоить библиотеку subprocess (оттуда понадобятся Popen, PIPE и метод communicate). Полезно знать, что blastp выдаст результат на stdout, если не задать выходной файл.
Всё можно делать на kodomo. Для пользователей Windows: если всё же захотите установить BLAST на свой компьютер, не скачивайте последнюю версию (2.10.0), у неё под Windows возникают проблемы с makeblastdb. Установите предыдущую (2.9.0), вот ссылка на установщик: https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.9.0/ncbi-blast-2.9.0+-win64.exe