Standalone BLAST
Поиск в геноме участков, кодирующих белки, похожих на YQGN_BACSU
Определим, закодированы ли похожие на YQGN_BACSU белки в геноме из другого организма (S.agalactiae), не пользуясь аннотацией генома.
Программа makeblastdb принимает файл с последовательностями в fasta-формате (параметр "-in") и создаёт локальную базу данных.
Параметр "-dbtype" указывает на тип последовательности – в случае нуклеотидной последовательности нужно указать "-dbtype nucl".
Создадим с помощью этой комнады базу, с которой будем работать:
sa_genome.fasta – полный геном бактерии Streptococcus agalactiae. Дальнейший поиск будем осуществлять программой tblastn, так как мы ищем гомолог белка в формальной трансляции нуклеотидного банка.
Получим результат программы в файле tblastn.out
Таблица 1. Поиск гомологов белка YQGN_BACSU в геноме S.agalactiae
Число находок с E-value < 0,001 |
1 |
E-value лучшей находки |
2e-27 |
Название последовательности с лучшей находкой |
AL766845 Streptococcus agalactiae NEM316 complete |
Координаты лучшей находки (от-до) |
76831 - 77343 |
Доля последовательности белка YQGN_BACSU, |
38% |
Поиск гомологов некодирующих последовательностей программой BLASTN
В файле trna_bacsu.fasta лежат последовательности всех тРНК, проаннотированных в полном геноме Bacillus subtilis BSn5.
Нам нужно определить, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии (Streptococcus agalactiae).
Запустим программу blastn, указав в качестве последовательностей для поиска файл trna_bacsu.fasta, в качестве банка – геном бактерии,
отформатированный в предыдущем задании.
Установим табличный формат выдачи (опция "-outfmt 6" или "-outfmt 7"). Порог на E-value установим равным 1e-2.
Получим файл - blast_standart.out.
Поиск гомологов при изменённых параметрах программы BLASTN
Повторим предыдущее задание ещё два раза с изменёнными параметрами программы. В первый раз изменим весовую матрицу, то есть параметры -reward и -penalty. Установим -reward 5 и -penalty -4. При этом программа откажется работать, пока мы не поменяем также параметры -gapopen и -gapextend (-gapopen 10 -gapextend 6).
Получим файл - blast_2.out.
Во второй раз, оставив те же (изменённые по сравнению со значениями по умолчанию) значения параметров -reward, -penalty, -gapopen и -gapextend, поменяем также значение параметра -word_size на минимально возможное (-word_size 4).
Получим файл - blast_3.out.
Теперь переведем всю полученную информацию после работы программ в удобный для просмотра вид. С помощью команд UNIX из полученных
файлов blast_standart.out, blast_2.out, blast_3.out вынесем в таблицу Excel названия последовательностей из файла trna_bacsu.fasta
и количество для каждого из них числа находок (гомологов каждой из тРНК в геноме бактерии Streptococcus agalactiae).
Файл с результатами можно посмотреть здесь (trna.xls).
Анализ результатов
Как видно из файла trna.xls, вцелом, количество найденных предположительно гомологичных последовательностей с переходом от стандартных настроек к настройкам с -reward 5 -penalty -4 -gapopen 10 -gapextend 6 и -word_size 4 довольно сильно увеличивается. Это можно объяснимо тем, что при изменении параметров мы ослабляем параметры, по которым программа определяет гомологов.
В выходном файле blast_3.out выберем какую-нибудь пару: tRNA B.subtilis и гомологичный участок, найденный в геноме бактерии Streptococcus agalactiae, которая находится при данном наборе параметров и не находится при других (BSn5_t20892 tRNA-Asn - [28016:28076]). Вырежем гомологичный участок в отдельный файл командой seqret -sask. Выделим исходную последовательность также в отдельный файл. Выровняем две последовательности программой needle.
Read and write (return) sequences
Begin at position [start]: 28016
End at position [end]: 28076
Reverse strand [N]: N
output sequence(s) [al766843.fasta]:
needle BSn5_t20892.fasta al766843.fasta -out task_4_out.fasta
Привожу в отчете выравнивание программой needle.
########################################
# Program: needle
# Rundate: Sat 14 Dec 2013 16:57:25
# Commandline: needle
# [-asequence] BSn5_t20892.fasta
# [-bsequence] al766843.fasta
# -outfile task_4_out.fasta
# Align_format: srspair
# Report_file: task_4_out.fasta
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: BSn5_t20892
# 2: AL766843
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 74
# Identity: 44/74 (59.5%)
# Similarity: 44/74 (59.5%)
# Gaps: 14/74 (18.9%)
# Score: 146.0
#
#
#=======================================
BSn5_t20892 1 gctctagtagcacagc-ggatagtgcagcagtttcctaaactgcaggtcg 49
|||||.|||| ||.|||.||..|.|.|||.||..|.|.||||||
AL766843 1 ------gtagctcagctggctagagcgtccggttcatacccgggaggtcg 44
BSn5_t20892 50 ggagttcgaatctctcctagagcg 73
||.||||||..|.||||
AL766843 45 ggggttcgatcccctcc------- 61
Identity выравнивания 59.5%, что не так много. Это говорит о том, что последовательности отдаленно гомологичны.
Время работы программы BLASTN
Время работы в зависимости от настроек алгоритма поиска отличается. Полученные данные представлены в Табл.2
Таблица 2. Время работы программы BLASTN
Настройки алгоритма |
Время |
Стандартные |
real 0m0.231s |
-reward 5 -penalty -4 -gapopen 10 -gapextend 6 |
real 0m0.276s |
-reward 5 -penalty -4 -gapopen 10 -gapextend 6 -word_size 4 |
real 0m11.979s |
Из таблицы 2 видно, что меньше всего времени программа затрачивает на поиск со стандартными настройками. Переопределение весов не сильно усложняет работу, а вот изменение длины слова приводит к затратам значительно большего времени.