Standalone BLAST

1. Поиск в геноме участков, кодирующих белки, похожие на заданный

Определим, закодированы ли похожие на белок PDAA из Bacillus subtilis белки в геноме Streptococcus agalactiae организма, не пользуясь аннотацией генома. Для этого выберем подходящую для решения данной задачи программу из пакета BLAST+ и запустим ее локально (standalone). Наиболее подходящей является программа tblastn, так как она осуществляет поиск гомологов нуклеотидной последовательности в нуклеотидной базе данных.

Создадим для начала локальную базу данных на основе полного генома бактерии Streptococcus agalactiae (sa_genome.fasta):

makeblastdb -in sa_genome.fasta -dbtype nucl

И получим файл(pdaa_bacsu.fasta) с аминокислотной последовательностью моего белка (PDAA_BACSU):

seqret sw:pdaa_bacsu

Далее осуществим поиск гомологов:

tblastn -query pdaa_bacsu.fasta -db sa_genome.fasta -evalue 0.001 -out tblast.out

Поиск гомологов белка PDAA в геноме Streptococcus agalactiae
Число находок с E-value < 0,001 2
E-value лучшей находки 2e-31
Название последовательности с лучшей находкой (embl) AL766847
Координаты лучшей находки (от-до) 83757-84350
Доля последовательности вашего белка, вошедшая в выравнивание с лучшей 112/205=0.55

2. Поиск гомологов некодирующих последовательностей программой BLASTN

Определим, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии (Streptococcus agalactiae).Для этого запустим команду 1:

blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out 1blastn.out -evalue 0.01 -outfmt 6

С помощью grep и Excel узнаем сколько находок было для каждой тРНК, получаем файл trna.xls (см. comand1). Использованный скрипт (как писать тут).


3. Поиск гомологов при изменённых параметрах программы BLASTN

Повторим предыдущее задание при измененных параметрах.

Сначала изменим весовую матрицу (-reward и -penalty), теперь будет -reward 5 и -penalty -4. А так же поменяем -gapopen и -gapextend на 25 и 10 (можно было так же 10 и 6, 8 и 6), команда 2:

blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out 2blast.out -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 25 -gapextend 10

Получили файл 2blast.out.

Изменим еще и word_size на минимально возможное, т.е. на 4, команда 3:

blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out 3blast.out -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 25 -gapextend 10 -word_size 4

Получили файл 3blast.out.

Итоговый файл с названиями тРНК и числом находок, полученных от каждой команды (1-3) - здесь


4. Анализ результатов

При изменении весов подчета количество находок незначительно увеличивается. А при изменении длины слова на минимальное очень значительно увеличивается (что очень логично).

В одном из полученных при выполнении заданий 2 и 3 выходных файлов BLASTN выбрала пару: tRNA B.subtilis - BSn5_t20966 и гомологичный участок, найденный в геноме другой бактерии.

Получила что ее координаты 100237-100182 (=> с - цепи). Далее вырезала этот фрагмент из Streptococcus agalactiae. И выровним эту последовательность с последовательностями тРНК Bacillus subtilis программой needle:

needle 1.fasta trna_bacsu.fasta Needleman-Wunsch global alignment of two sequences Gap opening penalty [10.0]: Gap extension penalty [0.5]: Output alignment [al766856.needle]:

.

Получим:

# Aligned_sequences: 2
# 1: AL766856
# 2: BSn5_t20966
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 78
# Identity:      46/78 (59.0%)
# Similarity:    46/78 (59.0%)
# Gaps:          23/78 (29.5%)
# Score: 174.0
# 
#
#=======================================

AL766856           1 -------tagctcagttgg-tagtagcgcatgactgttaatcatgatgtc     42
                            ||||||||.||| ||| ||||||.|.|||.|||.|.|||.|||
BSn5_t20966        1 gggcctgtagctcagctggttag-agcgcacgcctgataagcgtgaggtc     49

AL766856          43 gcaggttcgagtcc--------------     56
                     |..|||||||||||              
BSn5_t20966       50 ggtggttcgagtccactcaggcccacca     77

Как видно из выравнивания фрагмент генома Streptococcus agalactiae не очень хорошо выравнивается с BSn5_t20966 - треониновой тРНК B.subtilis, хотя с этой тРНК лучше, чем со всеми остальными. Следовательно, можно сделать вывод, что поиск гомологов с минимальной длиной слова не очень хорош, т.к. он может запросто найти просто похожие последовательности, но они врятли будут гомологами.

BSn5_t20966 em|AL766856 78.57 56 12 0 8 63 100237 100182 6e-08 50.0

5*. Время работы программы BLAST

Заметим, что при прибавлении команд время работы увеличивается, но это и логично, тк программа постепенно прибавляет новые параметры, тем самым увеличивая время своей работы. Время команды 3 особенно отличается от времени ожидания результата от других команд, скорее всего это связано с тем, что длина 4 очень маленькая и результатов (хитов) много.


© Tishina Sofia, 2013