Задание 1. Работа с getorf
Необходимые параметры для getorf: -minsize 30 — стоит по умолчанию; -table 0 — стоит по умолчанию; -find 1 — указываем (выдача трансляций участков между старт- и стоп-кодонами). Результат работы на записи embl:D89965:
>D89965_1 [66 - 155] Rattus norvegicus mRNA for RSS, complete cds. MQFHPRLPAVLQVCAACDRYASLLPAQRRL >D89965_2 [56 - 169] Rattus norvegicus mRNA for RSS, complete cds. MISDAVSSATASSASSLRSMRSVRQSFASSTAALTRWP >D89965_3 [245 - 316] Rattus norvegicus mRNA for RSS, complete cds. MTLSLYRRRTFFTLPFITVLPNVA >D89965_4 [332 - 379] Rattus norvegicus mRNA for RSS, complete cds. MTTWPLRRTLTIVVTS >D89965_5 [163 - 432] Rattus norvegicus mRNA for RSS, complete cds. MALMHFQFTFKQFEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHM AVTAYAYYSCHELTPWLRIQSTNPVQKYGA >D89965_6 [433 - 350] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds. MPRTFVRGLYSVFVTKGSARDNYSKRTP >D89965_7 [341 - 297] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds. MWSSLVMVRPRWAIP >D89965_8 [218 - 3] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds. MLLRCSNCLNVNWKCIRAIWSKPPLSWQKTGVPIACCANLKHCWQSRMKLHRLSSPVTVT WCSQKTILLLSA >D89965_9 [294 - 1] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds. MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS
Найденная рамка считывания №5 соответствует CDS в поле FT оригинальной последовательности (за исключением стоп-кодона, который в записи EMBL включен в последовательность, а в выдаче getorf — нет).
Однако в записи из банка UniProt, на которую ссылается эта запись из EMBL, мы находим частичное соответствие девятой рамке считывания: она в точности соответствует аминокислотным остаткам 28—125 (всего их 176). А если поискать бластом вообще всю последовательность D89965, то окажется, что это вообще отрывок из генома кишечной палочки, и крысы тут совершенно ни при чем. Об этом упоминается и в записи UniProt: хотя изначально предполагалось, что эта последовательность из крысы, учитывая схожесть последовательностей, вероятно, это просто загрязнение от самой кишечной палочки, в которой, видимо, пытались экспрессировать нужный крысиный ген. Дальнейшее копание с помощью бласта показывает, что эта последовательность встречается во всевозможных штаммах кишечной палочки, но никогда — у крыс. Резюме: экспериментальная ошибка.
Задание 2. Поиск гомологов некодирующих последовательностей с помощью blastn
Запустим blastn с нужными параметрами на последовательностях тРНК; получим результат. Далее сначала вытащим список последовательностей:
grep -i '# Q' blastn_0.txt | sed -e 's/# Query: //' -e "s/ .*//" > list.txt
Потом для каждой найдем количество вхождений:
while read s; do grep -i $s blastn_0.txt | grep -v '# Query' | wc -l >> list0.txt; done < list.txt
Задание 3. Поиск гомологов при измененных параметрах blastn
Прогоним BLAST с разными параметрами. Процесс при работе по умолчанию занимает 0.4 секунды.
- blastn_1.txt —
-reward 5 -penalty -4 -gapopen 8 -gapextend 6(0.5 с) - blastn_2.txt —
-reward 5 -penalty -4 -gapopen 8 -gapextend 6 -word_size 4(39 c) - blastn_3.txt —
-word_size 4(29 с)
Получим аналогичным образом файлы list1, list2, list3, подклеим их к таблице: окончательный вариант.
Анализ результатов
При уменьшении значения параметра word_size программе приходится прорабатывать больше потенциальных выравниваний и процесс значительно замедляется.
Заметно, что параметры 5/–4/8/6 немного более свободные, чем те, что по умолчанию, а изменение word_size на минимальное значение (4) приводит к заметно большему количеству выравниваний.
Для сравнения выберем вторую последовательность (BSn5_t20966 tRNA-Ile). При наборе параметров № 2 (см. выше) находится 10 выравниваний, а при наборе № 3 (изменение word_size) — 11. Это выравнивание имеет такие параметры: покрытие оригинальной последовательности — 67.57 (неплохо, и для сравнения подойдет); длина выравнивания — 74, несовпадений 24, гэпов 0 (ага, значит, почти точно это настоящее выравнивание), позиции выравнивания в оригинальной последовательности 4..77, а в геноме — 2168909..2168836 (), e-value = 6e-04.
Теперь вырежем оригинальную и гомологичную последовательности в отдельные файлы и построим выравнивание с помощью needle (чтобы избежать гэпов, которые в данном случае нежелательны и которых BLAST не находил, я поставил параметр gapopen 40, поскольку отдельного параметра nogap или подобного не нашел):
BSn5_t20966 1 gggcctgtagctcagctggttagagcgcacgcctgataagcgtgaggtcggtggttcgagtccactcaggcccacca 77
||.|||||||||.|||.||||||....||||..|||||....||||||..||||||.||....|.||..|.|||
CP000557 1 ---ccagtagctcagttggatagagcagcggccttctaagccgtcggtcgggagttcgaatctctcctgggacgcca 74
(Для сравнения выравнивание со всеми параметрами по умолчанию, погоды не делает, по-моему.)
Основные участки, которые соответствуют стеблям тРНК, показывают (на глаз) хорошую степень сходства. Попробуем что-нибудь более точное, чем прикидка, хотя бы mfold:
В соответствии с этим теперь раскрасим выравнивание по стеблям и посмотрим, что получится (в mfold антикодоновый стебель вместе с антикодоном немножко съехал, но править это я не буду, мы же только прикидываем):
BSn5_t20966 1 gggcctg tagctc agctggttagagc gcacgcctgataagcgtg aggtcggtgg ttcgagtccact cagg cccacca 77||.| |||||| ||.|||.|||||| ....||||..|||||... .||||||..| |||||.||.... |.|| ..|.||| CP000557 1 ---ccag tagctc agttggatagagc agcggccttctaagccgt cggtcgggag ttcgaatctctc ctgg gacgcca 74
Итак, это действительно тРНК. Осталось проверить, проаннотирована ли она соответствующим образом в геноме другой бактерии. Проаннотирована! Это tRNA-Arg. Интересно, что остатки, соответствующие антикодонам, как раз не совпадают (выделены жирным).
Осталось сделать вывод, почему это выравнивание находится при параметре word_size = 4, но не по умолчанию. Чтобы не гадать на кофейной гуще, попробуем постепенно увеличивать этот параметр и посмотрим, когда выравнивание пропадет. Опыт показывает, что выравнивание находится при word_size = 9, но не находится при значении 10. Если взглянуть на выравнивание, видно, что наибольший полностью совпадающий участок имеет длину как раз 9. Описание этого параметра в справке на сайте NCBI подтверждает вывод: этот параметр определяет минимальную длину "затравки" для выравнивания.