Поиск в геноме участков, кодирующих белки, похожие на заданный
С помощью Standalone BLAST программы tblastn в геноме бактерии Geobacillus thermodenitrificans проведён поиск гомологов белка YojM_BacSu из Bacillus subtilis. Нашёлся 1 гомолог, результаты представлены в таблице 1.
| Число находок с E-value < 0.001 | 1 |
| E-value лучшей находки | 1*10-35 |
| Название последовательности с лучшей находкой | CP000557.1 Geobacillus thermodenitrificans NG80-2, complete genome |
| Координаты лучшей находки | 2994751 : 2994272 |
| Доля последовательности белка, вошедшая в выравнивание с лучшей находкой | 196*3/3550319 = 0.0165% |
Поиск гомологов некодирующих последовательностей программой BLASTN
С помощью программы blastn проведён поиск гомологов всех тРНК Bacillus subtilis BSn5 в геноме бактерии Geobacillus thermodenitrificans. Результаты в файле trna.out 
. Далее было подсчитано количество гомологов для каждой тРНК (для этой и дальшейших задач использовались скрипты), и результат импортирован в Excel: trna.xls.
Поиск гомологов при изменённых параметрах программы BLASTN
Теперь проведён тот же поиск, но с изменёнными параметрами:
- Изменена весовая матрица. А именно следующие параметры приобрели новые значения:
-reward 5: бонус за match равен 5-penalty -4: штраф за mismatch равен 4-gapopen 10: штраф за открытие гэпа равен 10-gapextend 6: штраф за продление гэпа равен 6
Запрос имеет вид:blastn -task blastn -query /P/y12/term3/block3/tRNA/trna_bacsu.fasta -db gt.fasta -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -out trna_rp.out
Изменена весовая матрица (тем же образом) и длина слова для поиска выставлена минимальной - 4
Запрос имеет вид:blastn -task blastn -query /P/y12/term3/block3/tRNA/trna_bacsu.fasta -db gt.fasta -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -word_size 4 -out trna_rpws.out
Изменена только длина слова для поиска: на 4 нуклеотида
Запрос имеет вид:blastn -task blastn -query /P/y12/term3/block3/tRNA/trna_bacsu.fasta -db gt.fasta -evalue 0.01 -outfmt 6 -word_size 4 -out trna_ws.out
Для каждого случая подсчитано количество гомологов для каждой тРНК. Результаты внесены в файл trna.xls в колонки BLAST RP, BLAST RPWs и BLAST Ws соответсвенно для запросов 1, 2 и 3.
Анализ результатов
Описание различий
Примечательно, что каждый из изменнёных запросов в сумме находит больше гомологов, чем запрос с параметрами по умолчанию. Для того чтобы сравнить полученные результаты, была подсчитана разность между количеством всех найденных гомологов изменнёных запросов и запроса по умолчанию. Обозначим запрос с параметрами по умочание как запрос №0. Тогда:
- Δ(1,0) = 159
- Δ(2,0) = 1067
- Δ(3,0) = 293
Как видно из этих цифр, изменение длины слова приводит к большему числу найденных гомологов, чем изменение параметров матрицы. Особенно хорошо эта разница проявляется во втором запросе, когда изменение длины слова на порядок увеличивает число находимых гомологов (но вряд ли они все достоверные).
Сравнение находок
Рассмотрим находку для tRNA-Ser бактерии B.subtilis из запроса №3:
Query= BSn5_t21004 tRNA-Ser
Length=91
Score = 51.8 bits (56), Expect = 6e-08
Identities = 68/92 (74%), Gaps = 2/92 (2%)
Strand=Plus/Plus
Query 1 GGAGAAGTACTCAGGTGGCTGAAGAGG-CGCCCCTGCTAAGGGCGTAGGTCGCGTAAGCG 59
|||| ||||| || |||| ||||| || || | || || | | ||| || |||
Sbjct 554027 GGAGGAGTACCCAAGTGGTTGAAGGGGTCGGTCTTGAAAACCGAG-AGGCGTCGCAAGGC 554085
Query 60 GCGCGAGGGTTCAAATCCCTCCTTCTCCGCCA 91
||||| |||||| |||||||||| ||||||||
Sbjct 554086 GCGCGGGGGTTCGAATCCCTCCTCCTCCGCCA 554117
Она не находится в запросе №0. Это может быть объяснено тем, что в запросе №3 используется наименьшая допустимая длина слова для поиска (4 нуклеотида), т.е. blast будет пробовать выровнять каждые 4 нуклеотида одной последовательности против целой другой последовательности. Таким образом, у программы будет больший шанс найти совпадающий фрагмент, а значит и больше вероятность обнаружить гомолог.
Дополнительное выравнивание
Для дополнительного анализа были вырезаны последовательности tRNA-Ser и найденного гомолога и выровнены программой needle:
Length: 97
Identity: 68/97 (70.1%)
Similarity: 68/97 (70.1%)
Gaps: 12/97 (12.4%)
Score: 237.5
BSn5_t21004 1 ggagaagtactcaggtggctgaagagg-cgcccctg-----ctaagggcg 44
||||.|||||.||.||||.|||||.|| ||..|.|| |.|..||||
CP000557 1 ggaggagtacccaagtggttgaaggggtcggtcttgaaaaccgagaggcg 50
BSn5_t21004 45 taggtcgcgtaagcggcgcgagggttcaaatccctccttctccgcca 91
| ||.|||..|||||.||||||.||||||||||.||||||||
CP000557 51 t------cgcaaggcgcgcgggggttcgaatccctcctcctccgcca 91
Видно, что выравнивания blast и needle различаются. Сложно сказать почему, скорее всего это связано просто с разницей в алгоритмах, к тому же штраф needle за продление гэпа всего 0.5, из-за чего ему было выгодно менять mismatch на гэп. Если смотреть на это с биологической точки, то эта тРНК содержит вариационную петлю длиной около 20 нуклеотидов (нормальная длина тРНК ~70-75 нуклеотидов). Поэтому понятно, что начало и конец будут достаточно хорошо выравниваться (что мы и видим, особенно конец), а участок в середине - не очень (из-за этого, вероятно, пограммы и дают разные результаты касательно середины последовательности).
Проверка
В записи EMBL полного генома бактерии Geobacillus thermodenitrificans действительно содержится ген в позиции 554027..554117, который действительно кодирует сериновую тРНК:
FT gene 554027..554117 FT /locus_tag="GTNG_t046" FT tRNA 554027..554117 FT /locus_tag="GTNG_t046" FT /product="tRNA-Ser"Blast нашёл реальный гомолог.
Время работы программы BLAST
При запуске запросов blastn из заданий 2-3 проверялось время работы программ:
0: user 0m0.312s sys 0m0.252s 1: user 0m0.396s sys 0m0.248s 2: user 0m40.144s sys 0m0.288s 3: user 0m29.908s sys 0m0.220s
По логике, время работы программы напрямую зависит от количества находимых гомологов. Так и есть. Однако стоит заметить, что при поиске по слову длиной 4 время работы программы резко возрастает (с 0.3с до 29.9с), в то время как изменение параметров весовой матрицы тоже увеличивает время работы, но совсем незначительно (на ~0.1с). При этом изменение параметров матрицы при поиске по слову длиной 4 даёт прирост ко времени поиска в 10 секунд.