Программы пакета BLAST для работы с нуклеотидными последовательностями

1. Поиск в геноме участков, кодирующих белки, похожие на заданный

Задан белок OTC1_ECOLI из Escherichia coli K-12. Имеем аминокислотную последовательность этого белка: otc1_ecoli.fasta.
Для поиска в геноме Xanthomonas campestris участков, кодирующих сходные с заданным белки, будем использовать программу TBLASTN. Создадим в рабочей директории индексные файлы BLAST для нужного генома:
formatdb -i xc_genome.fasta -n xc -p F
Порог E-value устанавливаем равным 0,001. В командную строку вводим следующее:
blastall -p tblastn -d xc -i otc1_ecoli.fasta -o homologs.txt -e 0.001
В результате получаем файл с результатами поиска: homologs.txt. По данным из этого файла заполним таблицу:

Поиск гомологов белка OTC1_ECOLI в геноме Xanthomonas campestris

Число находок с Е-value<0,001 2
Характеристика лучшей находки:
   E-value находки 5e-30
Название геномной последовательности AE012332 AE008922 Xanthomonas campestris
pv. campestris str. ATCC 33913,
section 240 of 460 of the complete genome.
Координаты выравнивания(-ий) в найденной последовательности 7069 - 6071

2. Нахождение записи EMBL по последовательности с помощью программы BLASTN

Создадим с помощью программы seqret файл с последовательностью лучшей находки из прошлого задания: sequence.fasta.
Затем проведем поиск этой последовательности по банку "EMBL standard prokaryote" на сайте EBI (http://www.ebi.ac.uk/Tools/). Из списка результатов выберем первую последовательность (AC записи в EMBL: CP000050).
Выравнивание первой находки с последовательностью, по которой проводился поиск, находится здесь. Видно, что нам необходим участок последовательности с координатами 2263153 - 2264151.
Теперь при помощи entret создадим соответствующий файл EMBL (cp000050.entret), а в нем в поле FT найдем указанный участок последовательности: 
FT   gene            2263144..2264163
FT                   /locus_tag="XC_1869"
FT                   /note="XC1869"
FT   CDS             2263144..2264163
FT                   /codon_start=1
FT                   /transl_table=11
FT                   /locus_tag="XC_1869"
FT                   /product="ornithine carbamoyltransferase"
FT                   /db_xref="GOA:Q4UVJ1"
FT                   /db_xref="InterPro:IPR006131"
FT                   /db_xref="UniProtKB/TrEMBL:Q4UVJ1"
FT                   /protein_id="AAY48932.1"
FT                   /translation="MSLKHFLNTQDWSRAELDALLTQAALFKRNKLGSELKGKSIALVF
FT                   FNPSMRTRTSFELGAFQLGGHAVVLQPGKDAWPIEFNLGTVMDGDTEEHIAEVARVLGR
FT                   YVDLIGVRAFPKFVDWSKDREDQVLKSFAKYSPVPVINMETITHPCQELAHALALQEHF
FT                   GTPDLRGKKYVLTWTYHPKPLNTAVANSALTIATRMGMDVTLLCPTPDYILDERYMDWA
FT                   AQNVAESGGSLQVSHDIDSAYAGADVVYAKSWGALPFFGNWEPEKPIRDQYQHFIVDER
FT                   KMALTNNGVFSHCLPLRRNVKATDAVMDSPNCIAIDEAENRLHVQKAIMAALVGQSRP"
Как видим, искомый участок является частью аннотированной кодирующей последовательности, координаты которой 2263144..2264163. Эта CDS соответствует записи Q4UVJ1 в UniProt и кодирует орнитинкарбамоилтрансферазу.

3. Поиск гомологов с помощью программы BLASTN

Сначала получаем последовательность из генома E. coli, кодирующую исследуемый белок OTC1_ECOLI. Для этого при помощи entret получим запись EMBL c кодом доступа X00210 (см. задание по EMBL, пункт 4):
entret embl:X00210 -auto
В результате получен файл x00210.entret. Отсюда узнаем координаты CDS - это позиции 81..1085. Затем из записи EMBL вырезаем участок последовательности с этими координатами:
seqret "embl:X00210[81:1085]"
Получен файл с нужной нам последовательностью: x00210.fasta.

Далее при помощи BLASTN будем искать гомологи этого гена в геноме Xanthomonas campestris:
blastall -p blastn -d xc -i x00210.fasta -o homologs2.txt
Здесь предел значения E-value указывать не будем. В результате найдено 17 последовательностей с неприлично высокими значениями E-value: у первых двух находок 1.0, у остальных 4.0 (см. в файле homologs2.txt). Рассмотрим первую из находок с E-value = 1.0:
>AE012538 AE008922 Xanthomonas campestris pv. campestris str. ATCC
            33913,  section 446 of 460 of the complete genome.
          Length = 10676

 Score = 32.2 bits (16), Expect = 1.0
 Identities = 16/16 (100%)
 Strand = Plus / Plus

                            
Query: 280  gcccgcgtgcttggtc 295
            ||||||||||||||||
Sbjct: 6703 gcccgcgtgcttggtc 6718
Лучшая находка BLASTN не совпадает с лучшей находкой TBLASTN из п.1, они находятся в разных сегментах генома (здесь это 446-й, в п.1 был 240-й). Как видим, BLASTN находит практически идентичные, но очень короткие участки последовательностей. Сходства между более-менее удаленными друг от друга гомологами эта программа не видит: в списке результатов нет двух последовательностей, найденных в пункте 1 программой TBLASTN. В отличие от TBLASTN, здесь выравниваются не аминокислотные, а нуклеотидные последовательности, и не учитывается то, что разные триплеты могут кодировать одну и ту же аминокислоту. Кроме того, при выравнивании аминокислотных последовательностей учитывается также сходство аминокислот в случае неидентичности, а здесь выравниваются только одинаковые нуклеотиды.

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

При помощи программы entret получим файл с записью D89965 банка EMBL: d89965.entret.
entret embl:D89965 -auto
Затем, пользуясь бактериальныи кодом, получим из этой записи набор трансляций всех открытых рамок (ORF) данной последовательности длиной более 30 нуклеотидов. Открытой рамкой будем считать последовательность триплетов, начинающуюся со старт-кодона и заканчивающуюся стоп-кодоном. Для этого в командную строку вводим следующую команду:
getorf -table 11 -find 1
Минимальный размер рамки в 30 нуклеотидов задан по умолчанию.
Выходной файл: d89965.orf.

Данная запись EMBL ссылается на запись Swiss-Prot с кодом доступа P0A7B8. Ей соответствует 13-я рамка. Кроме того, 5-я рамка соответствует CDS, приведенной в рассматриваемой записи EMBL.

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

В файле trna_ecoli.fasta лежат последовательности всех тРНК, проаннотированных в полном геноме E.coli K12. Задача — определить, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии Xanthomonas campestris.

Вначале при помощи blastn найдем сходные последовательности в заданном геноме: 
blastall -p blastn -d xc -i trna_ecoli.fasta -m 8 -o rna.txt
Выходной файл: rna.txt.

Далее при помощи grep найдем число находок для каждой последовательности тРНК из E.coli. Для этого сначала нужно получить список этих тРНК командой 
grep ">" trna_ecoli.fasta >> names.txt
(список лежит в файле names.txt),
затем для каждой последовательности из списка ввести команду 
grep -c 'valV' rna.txt >> counter.txt
Из списка таких команд для всех тРНК, пользуясь функцией Excel CONCATENATE ("СЦЕПИТЬ"), составляем скрипт , при помощи которого получим файл с числами находок: counter.txt.

Теперь проведем еще один поиск с пороговым E-value = 0.001. Результаты поиска находятся в файле rna2.txt. Аналогично при помощи скрипта script2.scr получим для каждой тРНК E.coli число находок с заданным пороговым E-value (в файле counter2.txt).

Таблица Excel с результатами обоих поисков: trna.xls

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

Необходимо повторить предыдущее задание, используя вместо BLASTN megablast. Для этого нужно ввести команду:
megablast -d xc -i trna_ecoli.fasta -m 8 > mb.txt
Выходной файл: mb.txt

Discontigous megablast: 
megablast -d xc -i trna_ecoli.fasta -m 8 -W 11 -t 16 -N 1 > mb_d.txt
Выходной файл: mb_d.txt
Для discontigous megablast дополнительно указываем следующие значения опций:
-t - длина слов из тРНК, которые будут искаться в геноме (может равняться 16, 18 либо 21, мы выберем 16);
-W - длина слов из генома бактерии, по которым идет поиск последовательности (может равняться 11 или 12, выбираем 11);
-N - задает тип разрывов в матрице (может принимать значения 0, 1, 2, мы выберем 1).
Далее при помощи скриптов script3.scr и script4.scr получим из выдач megablast списки чисел находок и добавим эти данные в таблицу trna.xls.

Сравнение результатов megablast и blastall
Возьмем тРНК leuZ, которую нашел BLASTN и не нашел megablast. BLASTN выровнял с этой тРНК (позиции 1-45) участок генома Xenthomonas campestris AE012315 (позиции 5511-5555, т.е. это участок прямой цепи):
>AE012315 AE008922 Xanthomonas campestris pv. campestris str. ATCC
            33913,  section 223 of 460 of the complete genome.
          Length = 11231

 Score = 42.1 bits (21), Expect = 8e-05
 Identities = 39/45 (86%)
 Strand = Plus / Plus

                                                         
Query: 1    gcccggatggtggaatcggtagacacaagggatttaaaatccctc 45
            |||||| ||| | ||| ||||||| ||||||| ||||||||||||
Sbjct: 5511 gcccgggtggcgaaattggtagacgcaagggacttaaaatccctc 5555
При помощи seqret вырежем эти участки в отдельные файлы:leuZ.fasta и AE012315.fasta и выровняем их с помощью needle:
needle leuZ.fasta AE012315.fasta leuZ.needle
Получилось выравнивание:
########################################
# Program: needle
# Rundate: Tue 29 Dec 2009 13:44:49
# Commandline: needle
#    [-asequence] leuZ.fasta
#    [-bsequence] AE012315.fasta
#    [-outfile] leuZ.needle
# Align_format: srspair
# Report_file: leuZ.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: leuZ
# 2: AE012315
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 45
# Identity:      39/45 (86.7%)
# Similarity:    39/45 (86.7%)
# Gaps:           0/45 ( 0.0%)
# Score: 171.0
#
#
#=======================================

leuZ               1 gcccggatggtggaatcggtagacacaagggatttaaaatccctc     45
                     ||||||.|||.|.|||.|||||||.|||||||.||||||||||||
AE012315           1 gcccgggtggcgaaattggtagacgcaagggacttaaaatccctc     45
Видим, что гэпов нет вообще и процент идентичности довольно высок (86,7%). При этом E-value меньше 0,001. Это говорит о том, что это выравнивание биологически осмысленно. Эти гомологичные участки не были найдены megablast'ом, т.к. в выравнивании нет полностью совпадающих участков длины 28. Таким образом, megablast нельзя применять для поиска таких коротких последовательностей.
Нужный нам участок генома Xanthomonas campestris - позиции 2463969 - 2464013. Он проаннотирован так:
FT   gene            2463969..2464052
FT                   /locus_tag="XCC2091"
FT   tRNA            2463969..2464052
FT                   /locus_tag="XCC2091"
FT                   /product="tRNA-Leu"
FT                   /note="Found by tRNAscan"


Назад