Создав индексные файлы для базы данных по геному Geobacillus thermodenitrificans, проведём tblastn по последовательности белка yxiM из Bacillus subtilis. Результаты приведены в таблице 1.
Таблица 1.
| Число находок с E-value < 0,001 | 1 |
| E-value лучшей находки | 1e-24 |
| Название последовательности с лучшей находкой | CP000557.1 Geobacillus thermodenitrificans NG80-2, complete genome |
| Координаты лучшей находки (от-до) | 2012702..2013325 |
| Доля последовательности белка, вошедшая в выравнивание с лучшей находкой | (372-179)/382 * 100% = 50,5% |
По созданной ранее базе данных проведём поиск тРНК Bacillus subtilis с помощью программы blastn. Для этого используем команду
blastn -task -blastn -outfmt 7 -query trna_bacsu.fasta -db gt_genome.fasta -out trna.out -evalue 0.01
Теперь нужно получить таблицу с числом находок для каждой тРНК. Для этого был использован следующий скрипт. Сведём в одну таблицу в Excel столбец с названиями последовательностей тРНК и столбец с полученными числами находок (итоговый файл).
Теперь проведём поиск с тем же запросом, по той же базе данных, но с изменёнными параметрами.
В первом случае изменим весовую матрицу, запрос имеет вид
blastn -task -blastn -outfmt 7 -query trna_bacsu.fasta -db gt_genome.fasta -out trna.out -evalue 0.01 -penalty -4 reward 5 -gapopen 10 -gapextend 6
Во втором случае изменим весовую матрицу и размер слова на минимальный, запрос имеет вид
blastn -task -blastn -outfmt 7 -query trna_bacsu.fasta -db gt_genome.fasta -out trna.out -evalue 0.01 -penalty -4 reward 5 -gapopen 10 -gapextend 6 -word_size 4
В третьем случае изменим только размер слова
blastn -task -blastn -outfmt 7 -query trna_bacsu.fasta -db gt_genome.fasta -out trna.out -evalue 0.01 -word_size 4
В результате получим три файла с выравниваниями, используя почти тот же скрипт (можно просмотреть по ссылкам для первого, для второго и для третьего вариантов) подсчитаем число находок для каждой тРНК и добавим три столбца в xlsx-файл.
Заметим, что при изменении параметров количество находок изменятеся, происходит увеличение числа гомологов. Наиболее сильно этот эффект проявляется при изменении и длины слова, и весовой матрицы. Если же сравнивать между собой результаты поиска только с изменённой матрицей и только с изменённой длиной слова, то выигрывает длина слова. Можно сказать, что для при отборе гомологов длина слова играет большую роль.
Выберем гомолог, который есть в выдаче по второму запросу (изменена только весовая матрица) и отсутствует в выдаче по третьему (изменены весовая матрица и длина слова), и вырежем его программой seqret. Сохраним также в отдельный файл последовательность BSn5_t20966 tRNA-Ile, по которой он был найден. Выравняем программой needle. Полученное выравнивание:
CP000557 1 gattccgtagctcagctgg-gagagcgccac-cttga-cagggtggaggt 47
|...|.||||||||||||| .|||||| ||| |.||| .||.|| |||||
BSn5_t20966 1 gggcctgtagctcagctggttagagcg-cacgcctgataagcgt-gaggt 48
CP000557 48 cgctggttcgagcccagtcggaatca--- 73
||.|||||||||.|||.||.|...||
BSn5_t20966 49 cggtggttcgagtccactcaggcccacca 77
Идентичность выравнивания 70.9%, тогда как у выравнивания BLAST'ом - 75.38%. Это связано с более высоким штрафом за продолжение гэпа (6 вместо 0.5). В середине последовательности выровнены хуже, это связано с расположением там вариабельной петли, тогда как на концах выравнивание лучше.
Открыв запись об аннотированном геноме, например, на сайте NCBI, можно найти место, соответствующее координатам этой тРНК (819739-819811). На этом участке аннотирована валиновая тРНК, значит, BLAST нашёл верный гомолог.