1.Работа с программой getorf пакета EMBOSSБыла запущена программа getorf, чтобы получить набор трансляций всех открытых рамок данной последовательности - длиной более 30 нуклеотидов, - считая открытой рамкой последовательность триплетов,начинающуюся со старт-кодона и заканчивающуюся стоп-кодоном, - при использовании бактериального кода. Сначала был создан файл D89965.entret с записью банка EMBL. После этого была выполнена следующая команда:
getorf -sequence d89965.entret -minsize 30 -table 11 -find 1 -outseq d89965.orf
Параметры :
-sequence d89965.entret = файл с последовательностью, поданной на вход
-find = принимает параметры от 0 до 6.
0 транслирует рамку между стоп кодонами.
1 транслирует рамку между старт и стоп кодоном.
2 ищет нуклеотидную рамку между стоп кодонами.
3 ищет нуклеотидную рамку между старт и стоп кодоном.
4 нуклеотидная рамка фланкирующая к старт кодону
5 нуклеотидная рамка фланкирующая к инициаторному стоп кодону
6 нуклеотидная рамка фланкирующая к завершающему стоп кодону
-table = принимает параметра от 0 до 23.Коды для разных групп организмов (11 - для бактерий)
-minsize = минимальная длина рамки(в нуклеотидах)
-outseq d89965.orf = файл d89965.orf с результатом
На выходе получили файл d89965.orf В поле CDS файла d89965.entret указаны координаты участка [163-435]. Не нашлось такой рамки считывания, которая бы идеально описывала участок поля CDS, но наиболее похожей на нужную нам оказалась пятая [19 - 432] рамка. Чтобы сравнить полученные рамки с записью их в Swiss-Prot, получим нужный файл hslv_ecoli.entret: entret sw:p0a7b8Результаты для поиска пятой рамки приведены ниже:
>dbj|D89965.1| Geo Rattus norvegicus mRNA for RSS, complete cds
Length=448
Score = 289 bits (739), Expect = 2e-76, Method: Compositional matrix adjust.
Identities = 137/138 (99%), Positives = 138/138 (100%), Gaps = 0/138 (0%)
Frame = +1
Query 1 MVFWLHHVTVTGDDKRCSFIRDCQQCFKFAQHAIGTPVFCQLNGGFDQMALMHFQFTFKQ 60
+VFWLHHVTVTGDDKRCSFIRDCQQCFKFAQHAIGTPVFCQLNGGFDQMALMHFQFTFKQ
Sbjct 19 IVFWLHHVTVTGDDKRCSFIRDCQQCFKFAQHAIGTPVFCQLNGGFDQMALMHFQFTFKQ 198
Query 61 FEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHMAVTAYAYYSCHE 120
FEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHMAVTAYAYYSCHE
Sbjct 199 FEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHMAVTAYAYYSCHE 378
Query 121 LTPWLRIQSTNPVQKYGA 138
LTPWLRIQSTNPVQKYGA
Sbjct 379 LTPWLRIQSTNPVQKYGA 432
С помощью программы BLASTP видим, что записи hslv_ecoli.entret
банка Swiss-Prot соответствует только последняя рамка [375-1] с кодирующей последовательностью, расположенной на комплементарной цепи.
>sp|P0A7B8.2|HSLV_ECOLI RecName: Full=ATP-dependent protease hslV; AltName: Full=Heat
shock protein hslV
Score = 253 bits (647), Expect = 2e-67, Method: Compositional matrix adjust.
Identities = 125/125 (100%), Positives = 125/125 (100%), Gaps = 0/125 (0%)
Query 1 MTTIVSVRRNGHVVIAGDGQATLGNTVMKGNVKKVRRLYNDKVIAGFAGGTADAFTLFEL 60
MTTIVSVRRNGHVVIAGDGQATLGNTVMKGNVKKVRRLYNDKVIAGFAGGTADAFTLFEL
Sbjct 1 MTTIVSVRRNGHVVIAGDGQATLGNTVMKGNVKKVRRLYNDKVIAGFAGGTADAFTLFEL 60
Query 61 FERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPENDL 120
FERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPENDL
Sbjct 61 FERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPENDL 120
Query 121 IAIGS 125
IAIGS
Sbjct 121 IAIGS 125
2.Поиск некодирующих последовательностей программой BLASTN.1. Необходимо определить, сколько гомологов каждой из тРНК E.coli находит программа BLASTN в геноме родственной бактерии (Pasteurella multocida).
Для этого запускаем программу blast по банку из генома бактерии (без указания порога E-value)
2. Теперь нужно просмотреть выходной файл и придумать, как для данной последовательности из trna_ecoli.fasta запустить grep так,
чтобы на выходе получилось число - количество находок именно для данной последовательности.
grep "valV" trna.txt -c Результат: 32 Параметры: "valV" - выбранная последовательность, trna.txt - файл для поиска, -c - показывать только число находок. Далее я создала колонку из названий входных последовательностей командой grep ">" trna_ecoli.fasta > names Импортируем её в Excel. Далее необходимо создать скрипт из команд, выдающих число находок для каждой последовательности. Файл scr.scr был написан в редакторе Far и сохранен в соответствующем формате. Сделаем файл исполняемым: chmod +x script.scr Запуск файла осуществляется следующей командой:./script.scr Получаем файл trnacount.txt со столбцом цифр. Импортируем этот результат в Exel.Повторим поиск, на этот раз указав порог на E-value, равный 0.001. blastall -p blastn -d pm -i trna_ecoli.fasta -o trna_ev.txt -m 9 -e 0.001 Создадим скрипт: Запустим его: ./script2.scr Получаем файл trnacount2.txt со столбцом цифр. Добавляем в отчётную таблицу соответствующий столбец. Получили файл |
3.Поиск некодирующих последовательностей программой megablast.
1. Поиск гомологичных тРНК при помощи программы megablast.
megablast -d pm -i trna_ecoli.fasta -o megablast.txt -m 9
Полученный файл: megablast.txt
2. Поиск гомолочичных тРНК при помощи программы discontigous megablast.
megablast -d pm -m 8 -N 1 -W 11 -t 16 -i trna_ecoli.fasta -o discontigous.txt
Полученный файл: discontig.txt
Описание значений параметров:
Полученный отчетный Excel-файл: trna.xls
Следует отметить, что в параметах всех команд участвовало m9, т.е. скрип считывал еще и названия заголовков, поэтому в выходных файлах всех находок было ровно на 1 больше. Вручную это было исправлено.
4.Анализ результатов
В одном из полученных при выполнении заданий 2 и 3 выходных файлов BLAST была выбрана пара из tRNA E.coli и найденного в геноме другой бактерии гомологичного участка glyW - AE006138 (10908-10983). Эта пара находится программой BLASTN (e-value = 6e-31), но не находится программой megablast. С помощью программы seqret гомологичный участок был вырезан в отдельный файл:
seqret -sask
Reads and writes (returns) sequences
Input (gapped) sequence(s): pm_genome.fasta:Ae006138
Begin at position [start]: 10908
End at position [end]: 10983
Reverse strand [N]:
output sequence(s) [ae006138.fasta]:
Исходную последовательность также вырезали в отдельный файл.
Далее эти две последовательности были выровнены программой needle командой: needle glyW.fasta ae06138.fasta aln.needle
Вот, что получилось:
########################################
# Program: needle
# Rundate: Wed 18 Nov 2009 22:07:12
# Commandline: needle
# [-asequence] glyW.fasta
# [-bsequence] ae006138.fasta
# [-outfile] aln2.needle
# Align_format: srspair
# Report_file: aln2.needle
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: glyW
# 2: AE006138
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 76
# Identity: 73/76 (96.1%)
# Similarity: 73/76 (96.1%)
# Gaps: 0/76 ( 0.0%)
# Score: 353.0
#
#
#=======================================
glyW 1 gcgggaatagctcagttggtagagcacgaccttgccaaggtcggggtcgc 50
|||||||||||||||||||||||||||.|||||||||||||.||||||||
AE006138 1 gcgggaatagctcagttggtagagcacaaccttgccaaggttggggtcgc 50
glyW 51 gagttcgagtctcgtttcccgctcca 76
|||||||||.||||||||||||||||
AE006138 51 gagttcgagcctcgtttcccgctcca 76
#---------------------------------------
#---------------------------------------
Мы наблюдаем практически полное совпадение.
Не стоит забывать, что рограмма BLASTN ищет по словам длины 11, а megablast - 28. Если посмотреть на выравнивание, то видно, что в нем нет 28 совпавших нуклеотидов, расположенных подряд. Только в первой строчке выравнивания совпадает 27 нуклеотидов, в остальных и того меньше. Поэтому megablast и не нашел этот участок.
Как проаннотирован гомологичный участок в записи EMBL, описывающей геном бактерии?
FT gene 10905..11369 FT /gene="ccmH_1" FT /locus_tag="PM0012" FT CDS 10905..11369 FT /codon_start=1 FT /transl_table=11 FT /gene="ccmH_1" FT /locus_tag="PM0012" FT /product="CcmH" FT /db_xref="InterPro:IPR005616" FT /db_xref="UniProtKB/TrEMBL:Q9CPM5" FT /protein_id="AAK02096.1" FT /translation="MKKLTALFLLCLSFSSLAAIEGVQFSSTQQEKDYHALTQELRCPQ FT CQNNNIADSNATIAVDMRHKVLELLQEGKSKQDVVNFMVERYGHFVTYDPPLTVATVSL FT WVIPALFVILGFRLLFRRHTKQVVTAQASEPRLSDEQKQRLQRLLQEKKE"Можно сдеоать вывод, что данный гомологичный участок входит в состав гена ccmH_1.