Для того, чтобы получить набор трансляций всех открытых рамок считывания последовательности D89965, воспользуемся следующей командой:
getorf D89965.fasta -minsize 30 -find 1 D89965.out
Приведенной в поле FT кодирующей последовательности CDS (163-435) соответсвует следующая открытая рамка:
>D89965.1_5 [163 - 432] Rattus norvegicus mRNA for RSS, complete cds. MALMHFQFTFKQFEQRKSIRSTARKARDDFVVVQTADLFHVAFHYGIAQRGLTITSDDHM AVTAYAYYSCHELTPWLRIQSTNPVQKYGA
Записи P0A7B8 в Swiss-Prot, на которую ссылается запись EMBL, лучше всего соответствует другая рамка:
>D89965.1_9 [294 - 1] (REVERSE SENSE) Rattus norvegicus mRNA for RSS, complete cds. MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFELFERKLEMHQGHLVKAAVELAKDWRTDR MLRKLEALLAVADETASLIITGNGDVVQPENDLIAIGS
Последовательности совпадают лишь частично, открытая рамка как бы "обрезает" реальные ген в его начале и конце из-за наличия там метионина и стоп-кодона.
1 MTTIVSVRRNGHVVIAGDGQATLGNTVMKGNVKKVRRLYNDKVIAGFAGGTADAFTLFEL 60 P0A7B8 HSLV_ECOLI
1 ---------------------------MKGNVKKVRRLYNDKVIAGFAGGTADAFTLFEL 33 D89965.1_9
*********************************
61 FERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPENDL 120 P0A7B8 HSLV_ECOLI
34 FERKLEMHQGHLVKAAVELAKDWRTDRMLRKLEALLAVADETASLIITGNGDVVQPENDL 93 D89965.1_9
************************************************************
121 IAIGSGGPYAQAAARALLENTELSAREIAEKALDIAGDICIYTNHFHTIEELSYKA 176 P0A7B8 HSLV_ECOLI
94 IAIGS--------------------------------------------------- 98 D89965.1_9
*****Причем если сравнить всю последовательность D89965 и данного белка, видно, что сходство еще больше.
Стоит отметить, что записи соответствуют разным организмам: из EMBL - Rattus norvegicus, а Swiss-Prot - E.coli. Однако в Uniprot есть замечание, что "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". Так как Swiss-Prot - курируемая база данных, ей мы доверяем больше, чем EMBL.
Для нахождения гомологов использовалась следующая команда:
time blastn -query trna_bacsu.fasta -db gt -out trna_gt_nucl.out -evalue 0.01 -task blastn -outfmt 6
Выходной файл - trna_gt_nucl.out.
Получим сначала файл со списком tRNA:
grep 'BSn5_t[0-9]{5}' trna_bacsu.fasta -o -E >> list.txtА потом с количеством их повторений в файле выхода:
while read l; do grep $l -c trna_gt_nucl.out >> list2.txt; done < list.txt
Использованные команды:
time blastn -query trna_bacsu.fasta -db gt -out trna_gt_nucl1.out -evalue 0.01 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -task blastn -outfmt 6
time blastn -query trna_bacsu.fasta -db gt -out trna_gt_nucl2.out -evalue 0.01 -word_size 4 -reward 5 -penalty -4 -gapopen 10 -gapextend 6 -task blastn -outfmt 6
time blastn -query trna_bacsu.fasta -db gt -out trna_gt_nucl3.out -evalue 0.01 -word_size 4 -task blastn -outfmt 6
Результаты
Видно (да и понятно с самого начала), что параметры -reward 5 -penalty -4 min -word_size являются самыми нестрогими, поиск занимает больше всего времени, результатов, соответственно, тоже больше всего. Значения же по умолчанию - самые строгие, результатов меньше всего. Естественно, все результаты, найденные по умолчанию, находятся при других использованных параметрах тоже.
Для выполнения этого задания была выбрана tRNA BSn5_t21000. Она встречается в результатах поиска со стандартными параметрами 4 раза, с измененными штрафами - 7, с измененныим штрафами и длиной "слова" - 9, стандартными + минимальная длина "слова" - 6.
>BSn5_t21000 tRNA-Gln tgggctatatccaagcggtaaggcaacggattttgactccgtcatgcgttggttcgaatc cagctagcccagtca Нужный гомолог: 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. Принципииально эти тРНК могут отличаться только антикодоном, поэтому естественно, что они были найдены как гомологичные.