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


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


Получим файл с записью D89965 банка EMBL. Для этого отыщем его на сайте SRS или используем команду 
entret embl:D89965
С помощью tfm getorf получим все необходимые данные о программе. Так, чтобы получить набор трансляций всех открытых рамок данной последовательности D89965, нужно:
Запустим getorf с указанными параметрами на последовательности из записи D89965. 
                                                            
 getorf embl:D89965 -outseq D89965_30_ORF.fasta -minsize 30 -find 1
Или при логичном требовании длины белка в 60 и более а.о. 
 
 getorf embl:D89965 -outseq D89965_ORF.fasta -minsize 180 -find 1
В результате получим файлы: D89965_30_ORF.fasta и D89965_ORF.fasta.

Среди найденных открытых рамок одна полностью (если не считать разночтения со стоп-кодоном) соответствует единственной, приведённой в поле FT кодирующей последовательности (CDS). Это рамка №5 (№1 из 3 длинных), [163 - 432] Rattus norvegicus mRNA for RSS, complete cds.

Создадим файл с последовательностью записи Swiss-Prot P0A7B8_full.fasta, на которую ссылается данная запись EMBL.

С помощью blastp 
 blastp -query D89965_30_ORF.fasta -subject P0A7B8_full.fasta -out ORF_1.txt
выясним, что эта последовательность соответствует (не считая самых коротких и слабо идентичных рамок) только последней ORF (Query= D89965_9 [294 - 1]) с идентичностью 100% (файл выравнивания).

Итак, очевидны явные несоответствия данных. Во-первых, исходная м-РНК принадлежит крысе, а ссылка в описании белка приводит нас к белку E. Coli.
Во-вторых, описанная в файле открытая рамка считывания не совпадает с рамкой считывания, сответствующей белку, на который ссылаются авторы. При этом, выравнивание этого, изначально описанного гена показывает, что такого белка не описано вовсе (по крайней мере пока).

Разрешить возникшие несоответствия нам поможет допущение, что при секвенировании была получена м-РНК кишечной палочки, которую случайно выделили из экстракта крысиных клеток.
В таком случае все становится более объяснимым, жаль только, что ликвидировать несостыковки, похоже, никто не собирается.



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

Используя последовательности всех тРНК, проаннотированных в полном геноме Bacillus subtilis BSn5, определим, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии L.monocytogenes.

Для этого:
запустим программу blastn. 
blastn -query trna_bacsu.fasta -task blastn -db lm -out trna_lm.txt -evalue 0.01 -outfmt 7
Подсчитаем количество находок именно для одной конкретной последовательности тРНК командой grep. 
grep -c BSn5_t20894.e trna_lm.txt
С помощью 
grep ">" trna_bacsu.fasta > file.txt
получим названия всех тРНК, затем немного обработаем их для Excel, где сцепим данные в команды для скрипта командой 
 =СЦЕПИТЬ("grep -c '";A4;E4;"' trna_lm.txt >> count.txt")
После отладки скрипта с помощью указаний получим файл Excel_ex2.xls, где на листе 2 находится итоговая таблица отчета.



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


Повторим предыдущие действия ещё два раза с изменёнными параметрами программы, каждый раз сохраняя результаты в новый файл.

В первый раз изменим весовую матрицу, установив -reward 5 и -penalty -4. При этом придется также применить -gapopen 10 и -gapextend 6.

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

Результат занесем в отчетный файл Excel_ex2.xls в столбцах №2 и №3.

Строки с командами blastn для всех 4 (с дополнительным) столбцов отчета: 
                                                                                     

blastn -query trna_bacsu.fasta -task blastn -db lm -out trna_lm.txt -evalue 0.01 -outfmt 7

blastn -query trna_bacsu.fasta -task blastn -db lm -out trna_lm_1.txt -evalue 0.01 -outfmt 7 -reward 5 -penalty -4 -gapopen 10 -gapextend 6
                                                                                                                                                                                                                         
blastn -query trna_bacsu.fasta -task blastn -db lm -out trna_lm_2.txt -evalue 0.01 -outfmt 7 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -word_size 4

time blastn -query trna_bacsu.fasta -task blastn -db lm -out trna_lm_3.txt -evalue 0.01 -outfmt 7 -word_size 4
Дополнительные подсчеты

Прогоним Blastn еще раз по нашим последовательностям со значениями матрицы по умолчанию, но минимальным размером слов.

Также подсчитаем время выполнения комманд.

Как и ожидалось, получим, что очень значительно на время работы программы влияет -word_size, при изменении которого с 11 (по умолчанию) на 4 (минимум) увеличивает время выравниваний с 0.506 секунд до 20.775 секунд.
Дополнительные изменения матрицы значимы при малом размере слов, где время выравнивания увеличивается таким образом до 27.470 секунд.

Заметим, что все выравнивания проводились для сравнительно небольшого числа последовательностей (83 т-РНК/ 1 геном).



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

Если говорить в общем, как меняется число найденных гомологов при изменении параметров расчёта веса выравнивания и изменении длины слова, то можно заметить следующее:
В одном из полученных при выполнении заданий 2 и 3 выходных файлов BLASTN выберем пару из tRNA B.subtilis и найденного в геноме L.monocytogenes гомологичного участка.

Желательно было выбрать такую находку, которая находится программой BLASTN при одном наборе параметров и не находится при другом, и постараться объяснить причину этого.

Сравним две "пары", одна из которых (верхняя) находится всегда, а вторая - с меньшим размером слов и "лучшей" матрицей.
Можно заметить, что обе они выравнивают одни и те же молекулы, считываются в обратном направлении [41863 -41792] но процент идентичности второй гораздо хуже, то есть это скорее всего просто другой вариант выравнивания двух последовательностей, но с худшим весом. 

BSn5_t20894	embl|AL591983|AL591983	93.06	72	4	1	1	71	41863	41792	1e-15
BSn5_t20894	embl|AL591983|AL591983	69.12	68	21	0	2	69	41694	41627	6e-05
Вырежем гомологичный участок (AL591983 [41627 -41694], complement) в отдельный файл командой 
seqret lm_genome.fasta -sask
Исходную последовательность т-РНК также выделим в отдельный файл.
Выровняем две последовательности программой needle. 
                                                                              
needle -asequence lm_compl.fasta -bsequence trna_gln.fasta -outfile align_final.txt
Характеристики выравнивания оказались очень даже похожи на выдачу blastn с минимальным размером слов. 

# Aligned_sequences: 2                                                          
# 1: AL591983                                                                   
# 2: BSn5_t20894                                                                
# Matrix: EDNAFULL                                                              
# Gap_penalty: 10.0                                                             
# Extend_penalty: 0.5                                                           
#                                                                               
# Length: 71                                                                    
# Identity:      47/71 (66.2%)                                                  
# Similarity:    47/71 (66.2%)                                                  
# Gaps:           3/71 ( 4.2%)                                                  
# Score: 151.0                                                                  
#                                                                               
#                                                                               
#=======================================                                        
                                                                                
AL591983           1 -gcgccatagccaagtggtaaggcagaggtctgcaaaacctttatcaccg     49  
                      |.||.|||||||||.|||||||||..||.||...|..||.|.|||...|         
BSn5_t20894        1 tgggctatagccaagcggtaaggcaatggactttgactccgtgatcgttg     50  
                                                                                
AL591983          50 gttcaaatccggttggcgc--     68                               
                     ||||.|||||.|.|.||.|                                        
BSn5_t20894       51 gttcgaatccagctagcccag     71
Однако сперва посмотрим аннотацию гомологичного участка в поле FT записи EMBL генома бактерии L.monocytogenes. 

FT   tRNA            complement(41622..41695)
FT                   /product="transfert RNA-Cys"
Мы видим, что это, действительно, т-РНК, причем границы гена почти совпадают, но это другая т-РНК, не RNA-Gln, а RNA-Cys.
Однако это абсолютно не должно нас пугать.

Обратимся к выравниванию needle выше.

Можно выделить участки с высокой консервативностью (скорее всего это стебли) и достаточно вариабельные участки - петли.

Таким образом, очень вероятно, что эти т-РНК являются гомологами, пусть даже в ходе эволюции участок связывания с аминокислотой поменялся и стал переносить другую кислоту.
Нам остается только порадоваться, что мы все же сумели детектировать эту гомологию с помощью blastn, что с параметрами по умолчанию было бы невозможно в принципе. И все же для этого пришлось достаточно повозиться...




На страницу 3 семестра

© Aleshin Vasily