Getorf program. In search of homologs of non-coding sequences

Задание 1.

Работа с программой getorf пакета EMBOSS

Через систему SRS найдем запись D89965 банка EMBL.
Выполним команду

tfm getorf

чтобы понять, какое опции нужно использовать, чтобы получить набор трансляций всех открытых рамок данной последовательности длиной более 30 нуклеотидов, считая открытой рамкой последовательность триплетов, начинающуюся со старт-кодона и заканчивающуюся стоп-кодоном, при использовании стандартного кода:

getorf D89965.entret -table 0 -minsize 30 -find 1

-minsize - minimum nucleotide size of ORF to report; 30 - by default
-find - what to consider as an ORF; values (0 - by default):
- 0 - Translation of regions between STOP codons;
- 1 - Translation of regions between START and STOP codons;
- 2 - Nucleic sequences between STOP codons;
- 3 - Nucleic sequences between START and STOP codons;
- 4 - Nucleotides flanking START codons;
- 5 - Nucleotides flanking initial STOP codons;
- 6 - Nucleotides flanking ending STOP codons.
-table - code to use; values (0 - by default)
- 0 - Standard;
- 1 - Standard (with alternative initiation codons);
- 2 - Vertebrate Mitochondrial;
- 3 - Yeast Mitochondrial;
- 4 - Mold, Protozoan, Coelenterate Mitochondrial and Mycoplasma/Spiroplasma;
- 5 - Invertebrate Mitochondrial;
- 6 - Ciliate Macronuclear and Dasycladacean;
- 9 - Echinoderm Mitochondrial;
- 10 - Euplotid Nuclear;
- 11 - Bacterial;
- 12 - Alternative Yeast Nuclear;
- 13 - Ascidian Mitochondrial;
- 14 - Flatworm Mitochondrial;
- 15 - Blepharisma Macronuclear;
- 16 - Chlorophycean Mitochondrial;
- 21 - Trematode Mitochondrial;
- 22 - Scenedesmus obliquus;
- 23 - Thraustochytrium Mitochondrial;

После запуска getorf с указанными параметрами на последовательности из записи D89965 нашлась ORF, полностью соответсвующая, приведенной в записи в поле FT кодирующей последовательности (CDS):

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"

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

Данная запись EMBL ссылается на запись P0A7B8 в UniProtKB/Swiss-Prot:

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

Командой

entret sw:P0A7B8

создадим файл с записью Swiss-Prot (hslv_ecoli.entret) и вырежем оттуда последовательность в файл hslv_ecoli.fasta. (Можно просто скачать последовательность со страницы записи, найденной через SRS). Видим, что эта последовательность принадлежит кишечной палочке и очень удивляемся! Внимательно просматриваем запись SwissProt и находим предупреждение:

caution

Что означает, что все это - ошибка эксперимента.

Заметим (и не удивимся теперь), что девятая найденная рамка соответствует последовательности белка, но только частично:

>sp|P0A7B8|HSLV_ECOLI ATP-dependent protease subunit HslV;
MTTIVSVRRNGHVVIAGDGQATLGNTVMKGNVKKVRRLYNDKVIAGFAGGTADAFTL
FELFERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPE
NDLIAIGSGGPYAQAAARALLENTELSAREIAEKALDIAGDICIYTNHFHTIEELSYKA

>D89965_9 [294 - 1] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds.
MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR
MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS

Этот соответсвующий участок (28—125 а.о.) начинается с метионина (т.е. подходит как старт-кодон), а то, почему рамка заканчивается раньше, становится понятно, если вспомнить, что мы странным образом соотносим геномы кишечной палочки и крысы.

Задание 2.

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

Задача: определить, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии.

Так как в предыдущем задании мне была дана для сравнения бактерия Listeria monocytogenes (lm_genome.fasta), уже есть все файлы, отформатированные программой makeblastdb, т.е. три файла с расширениями (для нуклеотидной базы) nhr, nin и nsq: lm.nhr, lm.nin, lm.nsq.

В файле /P/y11/Term_3/Block_2/trna_bacsu.fasta лежат последовательности всех тРНК, проаннотированных в полном геноме Bacillus subtilis BSn5. Для удобства, скопируем все эти файлы в рабочую директорию.

Для последовательности trna_bacsu.fasta в качестве банка будем использовать отформатированный при выполнении задания 6 геном бактерии Listeria monocytogenes.

-query - input file name;
-db - BLAST database name;
-out - output file name;
-task - task to execute, values: 'blastn', 'blastn-short', 'dc-megablast','megablast' 'vecscreen', by default = `megablast';

-outfmt - alignment view options:

     0 = pairwise,
     1 = query-anchored showing identities,
     2 = query-anchored no identities,
     3 = flat query-anchored, show identities,
     4 = flat query-anchored, no identities,
     5 = XML Blast output,
     6 = tabular,
     7 = tabular with comment lines,
     8 = Text ASN.1,
     9 = Binary ASN.1,
    10 = Comma-separated values,
    11 = BLAST archive format (ASN.1)

-evalue - expectation value (E) threshold for saving hits; default = `10'.

Командой

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

выполняем поиск. Получили таблицу с выравниваниями lmtrna.txt.

Команда

grep -i '# Q' lmtrna.txt | sed -e 's/# Query: //' -e "s/ .*//" > list.txt

создаст список из названий входных последовательностей. grep ищет и выводит строчки, которые содержат '# Q', а sed -e выполняет скрипт для убирания ненужного (всего, кроме названия последовательности). Скрипт из команды s/// как раз заменит ненужное на пусто, как будто вырежет:

s/регулярное выражение/на что заменить/

Второе регулярное выражение содержит пробел, а после него конструкция .* — это любой символ (.), повторенный сколько угодно раз (*), что означает, что заменится все, что после пробела.

Теперь имеем list.txt список. Импортируем его в Excel.

Создадим скрипт из команд, выдающих число находок для каждой последовательности, как в примере. Команда

 grep -c 'BSn5_t20966' lmtrna.txt

выдаст на stdout число строк, содержащих слово "BSn5_t20966". Но т.к. слов в списке много, можно написать скприт wordcount.scr вида:

grep -c 'BSn5_t20894' lmtrna.txt >> wordcount.txt
grep -c 'BSn5_t20966' lmtrna.txt >> wordcount.txt
grep -c 'BSn5_t20968' lmtrna.txt >> wordcount.txt
...

Далее, выполним команду

noreturn wordcount.scr wc.scr

а затем сделаем файл со скриптом исполняемым:

chmod +x wс.scr

и запустим

 ./wс.scr

В результате в рабочей директории образуется файл wordcount.txt, в котором будет содержаться колонка чисел.

Результат работы скрипта импортируем в Excel. Получили отчетный Excel-файле (trna.xls) с двумя колонками:"Names" с названиями последовательнотей и "BLASTN default" с числами находок.

Задание 3.

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

Повторим предыдущее задание с изменёнными параметрами программы:

-gapopen 
   Cost to open a gap
-gapextend 
   Cost to extend a gap
-penalty 
   Penalty for a nucleotide mismatch
-reward =0>
   Reward for a nucleotide match
-word_size =4>
   Word size for wordfinder algorithm (length of best perfect match)

0. Параметры по умолчанию

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

Время: 0m0.405s

grep -i '# Q' lmtrna.txt | sed -e 's/# Query: //' -e "s/ .*//" > list.txt
noreturn wordcount.scr wc.scr
chmod +x wс.scr
./wс.scr

1. -gapopen 8 -gapextend 6 -reward 5 -penalty -4

time blastn -query trna_bacsu.fasta -task blastn -db lm -out lmtrna1.txt -gapopen 8
 -gapextend 6 -reward 5 -penalty -4 -outfmt 7 -evalue 0.01

Время: 0m0.493s

grep -i '# Q' lmtrna1.txt | sed -e 's/# Query: //' -e "s/ .*//" > list1.txt
noreturn wordcount1.scr wc1.scr
chmod +x wс1.scr
./wс1.scr

2. -gapopen 8 -gapextend 6 -reward 5 -penalty -4 -word_size 4

time blastn -query trna_bacsu.fasta -task blastn -db lm -out lmtrna2.txt -gapopen 8 
-gapextend 6 -reward 5 -penalty -4 -outfmt 7 -evalue 0.01 -word_size 4

Время: 0m27.426s

grep -i '# Q' lmtrna2.txt | sed -e 's/# Query: //' -e "s/ .*//" > list2.txt
noreturn wordcount2.scr wc2.scr
chmod +x wс2.scr
./wс2.scr

3.-word_size 4

time blastn -query trna_bacsu.fasta -task blastn -db lm -out lmtrna3.txt
 -outfmt 7 -evalue 0.01 -word_size 4

Время: 0m20.844s

grep -i '# Q' lmtrna3.txt | sed -e 's/# Query: //' -e "s/ .*//" > list3.txt
noreturn wordcount3.scr wc3.scr
chmod +x wс3.scr
./wс3.scr

Все фалы лежат в рабочей директории ~/Term3/Practice7/.
Окончательный вариант таблицы.

Задание 4.

Анализ результатов

После уменьшения значения параметра -word_size программа прорабатывает больше возможных выравниваний и процесс значительно замедляется, однако и находится больше выравниваний.Так как выравниваний при изменении стандартных параметров на -gapopen 8 -gapextend 6 -reward 5 -penalty -4 больше, можно понять, что они немного свободнее.

В качестве примера для анализа возьмем выравнивание BSn5_t20972 tRNA-Ile с последовательностью из записи embl AL591977.

Индентичность: 76.27%
Длина выравнивания: 59
Несовпадения: 13
Открытие gap: 1
В tRNA 5...63
В геноме 46660...46717
Evalue: 2e-04
Вес: 38.2

Выравнивание с помощью blastn:

> embl|AL591977|AL591977 Listeria monocytogenes strain EGD, complete 
genome, segment 5/12 
Length=270050

 Score = 38.2 bits (159),  Expect = 2e-04
 Identities = 45/59 (77%), Gaps = 1/59 (1%)
 Strand=Plus/Plus

Query  5      CTGTAGCTCAGCTGGTTAGAGCGCACGCCTGATAAGCGTGAGGTCGGTGGTTCGAGTCC  63
              | ||||||||| |||| |||||   || ||| ||| ||   |||||  |||||||||||
Sbjct  46660  CAGTAGCTCAGTTGGT-AGAGCAATCGGCTGTTAACCGATCGGTCGCAGGTTCGAGTCC  46717

Теперь вырежем последовательности в отдельные файлы (lm.fasta и t20972.fasta) и построим выравнивание с помощью needle:

needle lm.fasta t20972.fasta al.needle -auto

#=======================================
#
# Aligned_sequences: 2
# 1: AL591977
# 2: BSn5_t20972
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 77
# Identity:      45/77 (58.4%)
# Similarity:    45/77 (58.4%)
# Gaps:          19/77 (24.7%)
# Score: 163.0
# 
#
#=======================================

AL591977           1 ----cagtagctcagttgg-tagagcaatcggctgttaaccgatcggtcg     45
                         |.|||||||||.||| ||||||...||.|||.|||.||...|||||
BSn5_t20972        1 gggcctgtagctcagctggttagagcgcacgcctgataagcgtgaggtcg     50

AL591977          46 caggttcgagtcc--------------     58
                     ..|||||||||||              
BSn5_t20972       51 gtggttcgagtccactcaggcccacca     77


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

Получили файл с выраниванием al.needle.

Основные участки, которые примерно соответствуют стеблям тРНК, обнаруживают приличную степень сходства. Если прогнать этот кусок генома через mfold, получим весьма милую картинку:

caution

Гомологичный участок (из генома) проаннотирован в в поле FT записи EMBL, описывающей геном бактерии, так что blastn действительно выровнял нам тРНК с тРНК! Ура! Только наша исходная тРНК носила изолейцин, а у этой бактерии она носит аспарагин.

FT   tRNA            46656..46728
FT                   /product="transfert RNA-Asn"
FT                   /note="tRNAscan-SE vs 1.3 result - Cove score = 85.46"

Наверх