Standalone BLAST

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

Зная аминокислотную последовательность белка TENA_BACSU из Bacillus subtilis, нужно определить, закодированы ли похожие белки в геноме другого организма. С этой целью можно воспользоваться программами пакета BLAST.

Создадим индексные файлы пакета BLAST+ для поиска по геному бактерии Streptococcus agalactiae с помощью следующей команды:


  makeblastdb -in sa_genome.fasta -out sa -dbtype nucl

Затем воспользуемся программой TBLASTN для поиска гомологов белка TENA_BACSU по геному бактерии Streptococcus agalactiae. Для этого выполним соответствующую команду с необходимыми параметрами:


  tblastn -query tena_bacsu.fasta -db sa -out tena_sa.txt -evalue 0.001

По результатам поиска, сохранённым в файле tena_sa.txt, заполним таблицу.

Поиск гомологов белка TENA_BACSU в геноме бактерии Streptococcus agalactiae

Число находок с E-value < 0,001 1
E-value лучшей находки 6e-15
Название последовательности с лучшей находкой Streptococcus agalactiae NEM316 complete genome, segment 5 (AL766847)
Координаты лучшей находки 97660 - 98277
Доля последовательности белка TENA_BACSU, вошедшая в выравнивание с лучшей находкой 199/236 ≈ 0.84

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

Поиск был проведен по файлу, содержащему все тРНК из генома Bacillus subtilis. Был использован алгоритм BLASTN. Команда с помощью которой были получены первые результаты: 
 blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out trna.out -evalue 0.01 -outfmt 7
Затем с помощью команды grep из исходного файла был получен список названий тРНК, а из конечного - количество хитов на каждую тРНК. Проделанные выше операции повторили еще 2 раза, но с измененными параметрами. Была изменена весовая матрица. В первом случае были изменены параметры reward и penalty, равные 5 и -4 соответственно. После этого были определены параметры gapopen и gapextend, которые принимают значения 25 и 10 соответственно. Команда, с помощью которой были получены результаты: 
blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out trna.out -reward 5 -penalty -4 -gapopen 25 -gapextend 10 -evalue 0.01 -outfmt 7
Затем помимо этих двух параметров был задан параметр word_size со значением 7 (минимальное значение для работы программы). 
blastn -task blastn -query trna_bacsu.fasta -db sa_genome.fasta -out trna.out -reward 5 -penalty -4 -gapopen 25 -gapextend 10 -word_size 7 -evalue 0.01 -outfmt 7
Для выходных файлов команд с измененными параметрами был запущен скрипт, который использовался с параметрами по умолчанию. Таблица с количеством хитов на каждый из 3 запусков BLASTN приведена в файле.

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

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

Была выбрана треониновая тРНК: 

Это треониновая  тРНК из Streptococcus agalactiae
>AL766845 AL766845.1 Streptococcus agalactiae NEM316 complete genome, segment 3
gccggcttagctcagttggtagagcatctgatttgtaatcagagggtcgcgtgttcaagt
catgtagccggca
Это треониновая тРНК из Bacillus subtilis
>BSn5_t20996 tRNA-Thr
gccggtgtagctcaattggtagagcgcctgacttgtaagcagtggattgggggttcaagt
cctcttgccggca
Их выравнивание
########################################
# Program: needle
# Rundate: Tue  3 Dec 2013 21:50:07
# Commandline: needle
#    [-asequence] 5.fasta
#    [-bsequence] 6b.fasta
#    -outfile olo.needle
# Align_format: srspair
# Report_file: olo.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: AL766845
# 2: BSn5_t20996
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 73
# Identity:      58/73 (79.5%)
# Similarity:    58/73 (79.5%)
# Gaps:           0/73 ( 0.0%)
# Score: 230.0
# 
#
#=======================================

AL766845           1 gccggcttagctcagttggtagagcatctgatttgtaatcagagggtcgc     50
                     |||||..|||||||.||||||||||..||||.||||||.|||.||.|.|.
BSn5_t20996        1 gccggtgtagctcaattggtagagcgcctgacttgtaagcagtggattgg     50

AL766845          51 gtgttcaagtcatgtagccggca     73
                     |.|||||||||.|.|.|||||||
BSn5_t20996       51 gggttcaagtcctcttgccggca     73


#---------------------------------------
#---------------------------------------

Их выравнивание вроде бы очень хорошее. Наверное потому, что это тРНК, и без этой консервативной структуры она не смогла бы выполнять свои функции. И, кстати, в нотации EMBL в Streptococcus agalactiae это тоже треониновая тРНК. Почему она нашлась только в поиске с более коротким словом (см. параметры) - видимо это связано с какими-то отклонениями.

© Nosikova Kate, 2012