Standalone BLAST

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

Определим, закодированы ли белки похожие на URHG2_BACSU в геноме термофильной бактерии Bacillus licheniformis. Для этого создадим локальную базу данных на основе генома Bacillus licheniformis:

makeblastdb -in bl_genome.fasta -dbtype nucl

Для поиска гомологов белка используем программу TBLASTN c применением некоторых опций:

tblastn -query URHG2_BACSU.fasta -db bl_genome.fasta -out homologues.out -evalue 0.001

Таким образом, получаем файл homologues.out, в котором представлены данные о двух гомологах. Информация о лучшем из них показана в табл.1.

Табл.1. Результаты поиска гомолов белка URHG2_BACSU в геноме бактерии Bacillus licheniformis.

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

Запустим программу BLASTN, указав в качестве последовательностей для поиска файл trna_bacsu.fasta, а в качестве банка - геном бактерии Bacillus licheniformis bl_genome.fasta, отформатированный в предыдущем задании.

 blastn -task blastn -query trna_bacsu.fasta -db bl_genome.fasta -out blastn.out -outfmt 7 -evalue 0.01

Далее создаем колонку из названий входных последовательностей командой:

grep ">" trna_bacsu.fasta >grep.fasta

Создаем скрипт из команд, выдающих число находок для каждой последовательности. Для этого напишем простую программу на языке программирования Python: grep.py. Получаем скрипт: script.sh. Результаты работы скрипта приведены в Excel-файле: trna.xlsx.

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

Повторим проделанные выше операции еще по два раза, но изменив некоторые параметры программы, а именно весовую матрицу(-reward 5, -penalty -4, -gapopen 25 и -gapextend 10). В первом случае запускаем команду следующего вида:

 blastn -task blastn -query trna_bacsu.fasta -db bl_genome.fasta -out blastn1.out -outfmt 7 -evalue 0.01 -reward 5 -penalty -4 -gapopen 25 -gapextend 10

Во втором случае, оставив те же значения параметров -reward, -penalty, -gapopen и -gapextend, поменяем также значение параметра -word_size на минимально возможное, то есть на 4.

blastn -task blastn -query trna_bacsu.fasta -db bl_genome.fasta -out blastn1.out -outfmt 7 -evalue 0.01  -word_size 4 -reward 5 -penalty -4 -gapopen 25 -gapextend 10

Для выходных файлов применяем тот же скрипт, что и в предыдущем задании. Результаты представлены в Exel-файле: trna1.xlsx.

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

При изменении веса матрицы количество гамологов некоторых тРНК несколько увеличивается, как и при минимальной длине слова.

Для дальнейшего анализа выберем пару tRNA B.subtilis и гомологичный участок, найденный в геноме бактерии Bacillus licheniformis. К сожалению, среди трех выходов программы BLASTN нет пары, которая находится при одном наборе параметров и не находится при другом, поэтому выбираем любую понравившуюся, например BSn5_t20976 tRNA-Met и гомологичный участок AE017333 с координатами 74750...74824. Вырежем гомологичный участок бактерии Bacillus licheniformis в файл ae017333.fasta, а исходную последовательность тРНК B.subtilis в файл trna_bsub.fasta.Выровняем две последовательности программой needle и получим файл ae017333.needle.

# Aligned_sequences: 2
# 1: AE017333
# 2: BSn5_t20976
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 77
# Identity:      75/77 (97.4%)
# Similarity:    75/77 (97.4%)
# Gaps:           2/77 ( 2.6%)
# Score: 375.0
# 
#
#=======================================

AE017333           1 ggcggtgtagctcagctggctagagcgtacggttcatacccgtgaggtcg     50
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
BSn5_t20976        1 ggcggtgtagctcagctggctagagcgtacggttcatacccgtgaggtcg     50

AE017333          51 ggggttcgatcccctccgccgctac--     75
                     |||||||||||||||||||||||||  
BSn5_t20976       51 ggggttcgatcccctccgccgctacca     77

Выравнивание очень хорошее, гомологичный участок действительно соответствует tRNA-Met.

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

Для того, чтобы узнать время работы программы при изменении параметров, при запуске в командной строке перед всей командой допишем слово time. Тогда получим, что при стандартных параметрах время работы программы 0m0.573s, при изменении параметров весовой матрицы - 0m0.658s, а при изменении длины слова - 0m46.592s. Очевидно, что при добавлении параметров, программа затрачивает больше времени на выполнение команды.