Программа getorf.

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

Задание 1. Работа с программой getorf пакета EMBOSS

1. Воспользовавшись сервисом системы SRS, получили файл d89965.entret с записью D89965 банка EMBL. Его также можно было получить с помощью команды entret:

entret embl:D89965 -auto

2. Пользуясь бактериальным кодом, получим из этой записи набор трансляций всех открытых рамок (ORF) данной последовательности длиной более 30 нуклеотидов.
Открытой рамкой будем считать последовательность триплетов, начинающуюся со старт-кодона и заканчивающуюся стоп-кодоном.
Для этого в командную строку вводим следующую команду:

getorf d89965.entret d89965.orf -find 1

Минимальный размер рамки в 30 нуклеотидов задан по умолчанию.
На выходе получаем файл d89965.orf, содержащий 9 найденных открытых рамок.

3. При сравнении этих рамок с кодирующими последовательностями из записи D89965 банка EMBL, можно отметить, что пятая открытая рамка (163 - 432) практически полностью соответствует CDS (163 - 435), приведенной в рассматриваемой записи EMBL:

>D89965_5 [163 - 432] Rattus norvegicus mRNA for RSS, complete cds.
MALMHFQFTFKQFEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHM
AVTAYAYYSCHELTPWLRIQSTNPVQKYGA

FT   CDS             163..435
FT                   /product="RSS"
FT                   /note="Rat Stomach Serotonin receptor-related gene"
FT                   /db_xref="GOA:P0A7B8"
FT                   /db_xref="InterPro:IPR001353"
FT                   /db_xref="InterPro:IPR022281"
FT                   /db_xref="PDB:1E94"
FT                   /db_xref="PDB:1G4A"
FT                   /db_xref="PDB:1G4B"
FT                   /db_xref="PDB:1HQY"
FT                   /db_xref="PDB:1HT1"
FT                   /db_xref="PDB:1HT2"
FT                   /db_xref="PDB:1NED"
FT                   /db_xref="UniProtKB/Swiss-Prot:P0A7B8"
FT                   /protein_id="BAA14040.1"
FT                   /translation="MALMHFQFTFKQFEQRKSIRSTARKARDDFVVVQTADLFHVAFHY
FT                   GIAQRGLTITSDDHMAVTAYAYYSCHELTPWLRIQSTNPVQKYGA"

4. Из поля FT мы узнаём, что данная запись EMBL ссылается на запись Swiss-Prot с кодом доступа P0A7B8:

FT                   /db_xref="UniProtKB/Swiss-Prot:P0A7B8"

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

seqret sw:p0a7b8

5. С помощью команды blastp посмотрим какой из ранее найденных открытых рамок соответствует эта последовательность:

blastp -query hslv_ecoli.fasta -subject d89965.orf -out blastp.out

На выходе получаем файл blastp.out, из которого узнаём, что последовательность записи P0A7B8 соответствует девятой найденной нами открытой рамке (294 - 1):

Subject= D89965_9 [294 - 1] (REVERSE SENSE) Rattus norvegicus mRNA for RSS,
complete cds.

Length=98

Score = 200 bits (509), Expect = 2e-57, Method: Compositional matrix adjust.
Identities = 98/98 (100%), Positives = 98/98 (100%), Gaps = 0/98 (0%)

Query 28   MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR 87
                MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR
Sbjct 1    MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR 60

Query 88   MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS 125
                 MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS
Sbjct 61    MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS 98

При этом последовательность значительно длиннее как справа, так и слева. Это объясняется тем, что белок из Swiss-Prot принадлежит E.coli, которая используется в исследованиях для его синтеза с данного гена Rattus norvegicus. Вполне возможно, что полипептид, синтезированный в бактерии, в ходе модификаций соединяется с некоторыми последовательностями в начале и в конце, чтобы не причинить бактерии вреда.

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

С помощью программы blastn определим сколько гомологов каждой из тРНК (взятой из файла trna_bacsu.fasta с последовательностями всех тРНК, проаннотированных в полном геноме Bacillus subtilis BSn5) найдётся в геноме родственной бактерии Geobacillus thermodenitrificans:

blastn -query trna_bacsu.fasta -db gt -out trna_bacsu_gt_blastn.txt -outfmt 6 -evalue 0.01 -task blastn

Чтобы получить столбец со значениями числа находок для каждой последовательности в файле trna_bacsu_gt_blastn.txt, разделим задание на несколько пунктов:

1. С помощью команды grep составим список названий входных последовательностей:

grep ">" trna_bacsu.fasta > trnas.txt

2. Список названий последовательностей перенесём в Excel

3. Используя функцию "СЦЕПИТЬ" в Excel, создаём строки для скрипта вида:

grep -c 't2****' trna_bacsu_gt_blastn.txt >> trna_in_gt.txt

4. Делаем скрипт исполняемым:

chmod +x script_trna_in_gt.scr

5. Запускаем скрипт:

./script_trna_in_gt.scr

После чего получаем столбец со значениями числа находок для каждой последовательности.

6. Переносим полученные значения в Excel (столбец "BLASTN default").

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

Повторим поиск гомологов ещё раз, но на этот раз изменив определённые параметры:

1. Изменим параметры весовой матрицы → -reward 5 и -penalty -4. Из предложенных программой изменений значений параметров -gapopen и -gapextent выберем равные 10 и 6 соответственно.

blastn -query trna_bacsu.fasta -db gt -out trna_in_gt_1.txt -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -task blastn

Получим столбец со значениями.

2. Теперь, помимо изменений в 1 пункте, поменяем значение параметра -word_size на минимальное, то есть на 4.

blastn -query trna_bacsu.fasta -db gt -out trna_in_gt_2.txt -evalue 0.01 -outfmt 6 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -word_size 4 -task blastn

Получим ещё один столбец со значениями.

3. В третий раз выполним поиск при минимальном значении -word_size и при значениях параметров вычисления веса выравнивания по умолчанию.

blastn -query trna_bacsu.fasta -db gt -out trna_in_gt_3.txt -evalue 0.01 -outfmt 6 -word_size 4 -task blastn

Снова получим столбец со значениями.

Используя команду time, записываемую перед исполняемой командой, проследим за соотношением времени работы программы при изменении параметров:

reward	penalty	gapopen	gapextend	word_size	time: real	time: user	time: sys
1	-3	5	2	11	0m0.409s	0m0.216s	0m0.180s
5	-4	10	6	11	0m0.501s	0m0.328s	0m0.160s
5	-4	10	6	4	0m39.547s	0m39.370s	0m0.172s
1	-3	5	2	4	0m29.237s	0m29.042s	0m0.184s

Как мы видим, изменение длины слова (параметра word_size) со стандартного (11) на 4, в разы увеличило время работы blastn.

В Excel можно посмотреть окончательную таблицу результатов работы по заданиям 2 и 3.

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

Из предыдущих результатов выбрали tRNA BSn5_t21000, которая встречается:
1 - со стандартными параметрами → 4 раза
2 - с измененными штрафами → 7 раз
3 - с измененными штрафами и минимальной длиной слова → 9 раз
4 - со стандартными параметрами и минимальной длиной слова → 6 раз

>BSn5_t21000 tRNA-Gln
tgggctatatccaagcggtaaggcaacggattttgactccgtcatgcgttggttcgaatc
cagctagcccagtca

Нужный нам гомолог находится при изменённых параметрах (2-4), но не находится при стандартных параметрах:

BSn5_t21000 CP000557 75.76 66 15 1 3 67 554234 554299 7e-06 44.6

Это логично, так как изменяя параметры мы "ослабляем" требования к выдаче (смягчаем критерии поиска), следовательно получаем больше находок.

Вырежем нужный участок (гомологичный выбранной tRNA), командой seqret -sask, в результате получим:

>CP000557 CP000557.1 Geobacillus thermodenitrificans NG80-2, complete genome.
ggctatggcgaagtggttaacgcaccagattgtggctctggcatgcgtgggttcgattcc
cactag

Ниже приведены результаты выравнивания программой needle со стандартными параметрами:

########################################
# Program: needle
# Rundate: Tue 30 Oct 2012 20:21:51
# Commandline: needle
#    -asequence 21000nat.fasta
#    -bsequence 21000.fasta
# Align_format: srspair
# Report_file: bsn5_t21000.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: BSn5_t21000
# 2: CP000557
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 76
# Identity:      50/76 (65.8%)
# Similarity:    50/76 (65.8%)
# Gaps:          11/76 (14.5%)
# Score: 180.0
# 
#
#=======================================

BSn5_t21000        1 tgggctatatccaagcgg-taaggcaacggattttgactccgtcatgcgt     49
                       ||||||..|.|||.|| |||.|||.|.||||.||.|||.|.|||||||
CP000557           1 --ggctatggcgaagtggttaacgcaccagattgtggctctggcatgcgt     48

BSn5_t21000       50 tggttcgaatccagctagcccagtca     75
                     .|||||||.|||..||||        
CP000557          49 gggttcgattcccactag--------     66


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

Как мы видим идентичность последовательностей - 65,8%, так что они далеко не однозначные гомологи. Гомологичным был отмечен участок меньшей длины, чем тРНК. В середине последовательностей гэпов практически нет, встречаются достаточно продолжительные одинаковые участки, например, соответствующий акцепторному стеблю (вначале). При уменьшении длины слова смягчаются критерии поиска, и эта последовательность находится (даже при неизменных остальных параметрах); при стандартных же параметрах - нет.

Информация из поля FT генома бактерии (CP000557):

FT   gene            554232..554304
FT                   /locus_tag="GTNG_t048"
FT   tRNA            554232..554304
FT                   /locus_tag="GTNG_t048"
FT                   /product="tRNA-His"

Мы убеждаемся, что это действительно тРНК, но гистидиновая, а не глутаминовая. Плюс ко всему её длина действительно несколько больше, чем было найдено blastn.