Программа getorf

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

Создаю в своей директории файл с записью D89965 банка EMBL:

entret embl:d89965

Получен файл d89965.entret.

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

getorf d89965.entret d89965.orf -find 1

Получен файл d89965.orf.

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

seqret sw:p0a7b8

Файл d89965.orf содержит 9 открытых рамок. Сравниваю их с данными записи D89965 в поле FT кодирующей последовательности (CDS) на EMBL:

 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"

Рамка 5 почти полностью (последовательности идентичны, но рамка 5: 163-432, а в CDS указано: 163-435) соответствует приведенной в CDS последовательности:

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

Выясняю, что последовательность записи Swiss-Prot частично соответствует открытой рамке 9:

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

Об этом говорят данные, полученные в файле blastp.out: 

Query= HSLV_ECOLI P0A7B8 ATP-dependent protease subunit HslV (3.4.25.2)
(Heat shock protein HslV)

Length=176

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

Белок включает 176 остатков, в выравнивании участвует 98 остатков - не вошли начало и конец белка. При получении набора трансляций использовался стандартный код и открытой считалась рамка из более 30 нуклеотидов от старт- до стоп-кодона, что накладывает определенные ограничения на поиск. Видно, что некоторые рамки пересекаются друг с другом. Программа может не воспринять первый метионин из данной последовательности белка как старт-кодон из-за наличия стоп-кодона на растоянии меньше 30 нуклеотидов или отнести его к другой рамке считывания. Возможны другие причины, приводящие к неточностям: трансляционные модификации и различные вариации в генетическом коде.

В данном случае произошла ошибка при секвенировании, о чем говорится в предупреждении в поле CC: sequence is supposed to originate from rat but, based on sequence similarity, it seems that this is a case of bacterial contamination from E.coli. Вероятно, к тому же была отсеквенирована только часть последовательности, с потерянным концом и началом.

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

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

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

Time:

real 0m0.409s

user 0m0.216s

sys 0m0.180s

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

grep ">" trna_bacsu.fasta > trna.txt

Импортирую её в Excel.

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

grep -c 't20894' trna_bacsu_gt_blastn.txt >> trna_gt.txt

grep -c 't20966' trna_bacsu_gt_blastn.txt >> trna_gt.txt

...

Создаю скрипт trnagt.sh делаю его исполняемым, запускаю и импортирую результат работы в trna.xlsx.

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

Повторяю предыдущее задание. В первый раз изменяю весовую матрицу, то есть параметры -reward и -penalty, и параметры -gapopen и -gapextend:

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

Time:

real 0m0.501s

user 0m0.328s

sys 0m0.160s

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

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

Time:

real 0m39.547s

user 0m39.370s

sys 0m0.172s

В третий раз, использую минимальное значение word_size и параметры вычисления веса выравнивания, взятые по умолчанию:

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

Time:

real 0m29.237s

user 0m29.042s

sys 0m0.184s

Использовала команду time, записываемую перед исполняемой командой, чтобы узнать время работы программы при разных параметрах. Больше всего времени занимает выполнение команды сразу с пятью заданными параметрами, причем при выборе минимального значения word_size время поиска сразу увеличилось примерно на 29 секунд, т.е. возросло приблизительно в 70 раз, хотя поиск дал всего на 30% больше результатов. Добавление параметров -reward 5 -penalty -4 -gapopen 10 -gapextend 6 увеличило время поиска лишь на 0.1 секунды, причем гомологов нашлось также на 10% больше (было 570, стало 634).

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

В последнем полученном выходном файле BLASTN выбираю пару из tRNA BSn5_t20982 и найденного в геноме бактерии Geobacillus thermodenitrificans, CP000557 гомологичного участка 86461-86516. Данная пара находится программой BLASTN только при минимальном значении word_size, т.е. при предпоследнем и последнем поиске. Из выравнивания, приведенного ниже видно, что максимальная длина участка выравнивания без гэпов равна 9, поэтому, при значении параметра -word_size, равном 11, выравнивание не нашлось.

Вырезаю гомологичный участок в отдельный файл cp000557.fasta:

seqret CP000557.embl -sask

Выделяю исходную последовательность в отдельный файл t20982.fasta.

Выравниваю две последовательности программой needle:

needle t20982.fasta cp000557.fasta t20982.needle -auto

Получаю выравнивание:

BSn5_t20982        1 gcttccatagctcagcaggtagagcacttccatggtaaggaagaggtcag     50
                             |||||||||.||.|||||||||.|.|...|||.|||.||||.|
 CP000557           1 -------tagctcagctgggagagcacttgccttacaagcaaggggtcgg     43

 BSn5_t20982       51 cggttcgagcccgcttggaagct     73
                      ||||||||.||||
 CP000557          44 cggttcgatcccg----------     56

Вес выравнивания: 181.0, процент сходства и совпадения: 61.6%, гэпы: 23.3%, длина: 73.

Последовательности схожи, но длина гомологичного участка из CP000557 меньше, это приводит к добавлению гэпов в начале и в конце выравнивания и к сильному снижению процента сходства и совпадения.

Аннотация гомологичного участка в поле FT записи EMBL, описывающей геном бактерии:

       FT   gene            86454..86529
        FT                   /locus_tag="GTNG_t007"
        FT   tRNA            86454..86529
        FT                   /locus_tag="GTNG_t007"
        FT                   /product="tRNA-Val"

Т.е. участок, действительно, кодирует последовательностью тРНК.

© Eugenia Prokhorova 2011