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

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

   Имея аминокислотную последовательность белка DPS_ECOLI из Escherichia coli К-12, определим, закодированы ли похожие белки в геноме бактерии Xanthomonas campestris.
   Для этого вначале создадим в рабочей директории индексные файлы пакета BLAST для поиска по геному этой бактерии с помощью команды:
   

   formatdb -i xc_genome.fasta -p F -n xc

   После проведем поиск гомологов в программе TBLASTN с порогом E-value, равным 0,001:
   

   blastall -p tblastn -d xc -i dps.fasta -o homol.txt -e 0.001

   На основе полученного файла составим таблицу:
   

   Поиск гомологов белка DPS_ECOLI в геноме бактерии Xanthomonas campestris

   Число находок с Е-value<0,001		1
   Характеристика лучшей находки:
   	E-value находки	6 . 10 -10
   	Название геномной последовательности	AE012409
   	Координаты выравнивания в найденной последовательности	7268-7714




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

   C помощью программы seqret был получен файл с последовательностью участка генома, найденного в предыдущем упражнении. На сайте EBI был проведен поиск этой последовательности в банке "EMBL standard prokaryote". В результате была найдена запись со стопроцентным совпадением последовательностей. Ее AC - AE008922, это последовательность генома бактерии Xanthomonas campestris pv. campestris str. ATCC 33913. Найденный в предыдущем задании участок последовательности соответствует участкам последовательности с 3506767-го номера нуклеотидной пары до 3507213-й (цепь +).
   С помощью программы entret была получена запись EMBL последовательности (AE008922). Об участке 3506767-3507213 в поле FT имеется следующая информация:
   

   FT   gene            3506674..3507225
   FT                   /locus_tag="XCC2946"
   FT   CDS             3506674..3507225
   FT                   /codon_start=1
   FT                   /transl_table=11
   FT                   /locus_tag="XCC2946"
   FT                   /product="DNA-binding related protein"
   FT                   /note="identified by sequence similarity; putative; ORF
   FT                   located using Blastx/Glimmer/Genemark"
   FT                   /db_xref="GOA:Q8P6M1"
   FT                   /db_xref="HSSP:1JI5"
   FT                   /db_xref="InterPro:IPR012347"
   FT                   /db_xref="UniProtKB/TrEMBL:Q8P6M1"
   FT                   /protein_id="AAM42218.1"
   FT                   /translation="MSKKKNAANKPAAPTDALPVAAPSAPHIDIGIKDADRKQISDGLA
   FT                   RYMADAFTLYLKTHNFHWNVTGSMFNSLHTMFETQYTEQWGALDEVAERIRALGYNAPG
   FT                   SYREFVALTSIPEEPGLSDSADWREMVRQLVSGNEAVCRTARKVLGTADDAGDDPTVDL
   FT                   LTQRLQTHEKYAWMLRSLLQ"
   

   Как видно, кодирующая последовательность имеет координаты 3506674-3507225 и соответствует записи Q8P6M1 банка UniProt. Этот белок так же, как и DPS_ECOLI, принадлежит к семейству dps и без сомнений является его гомологом. Таким образом, BLASTN дал положительный результат и справился с поставленной задачей.



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

   Была взята запись EMBL X69337, в ней последовательность, кодирующая белок DPS_ECOLI, имеет координаты CDS 296-799. Программой seqret эта последовательность была вырезана в отдельный файл.
   Программой BLASTN былы найдены гомологи белка DPS_ECOLI в геноме бактерии Xanthomonas campestris (но на этот раз по нуклеотидной последовательности). Для этого была выполнена команда:
   

   blastall -p blastn -d xc -i parta.fasta -o homolog.txt

   В полученном файле содержится информация о выравниваниях последовательностей. Лучшая находка имеет E-value 0.032, причем это уже не та находка, что в первом упражнении. Длина выравнивания составляет всего 18 нуклеотидов. Остальные находки имеют E-value больше 2, что является совсем плохим результатом. В результате, можно предположить, что в геноме бактерии Xanthomonas campestris не имеется белков, гомологичных DPS_ECOLI. Находка, найденная в первом упражнении, не нашлась в этом, скорее всего, потому, что в первом задании имеющаяся аминокислотная последовательность может иметь очень много соответствующих нуклеотидных последовательностей (вследствие вырожденности генетического кода), поэтому одна находка случайным образом нашлась (не будучи гомологом). В третьем же задании условия были более жесткие (так как имелись две нуклеотидные последовательности), поэтому ни одной достоверной находки не нашлось.



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

   Вначале был получен файл с последовательностью записи D89965. С помощью команды "tfm getorf" изучил описание программы getorf. Использование бактериального кода задается значением параметра "-table 11", длина открытых рамок данной последовательности более 30 нуклеотидов задается значением параметра "-minsize 30", соответствие начала открытых рамок старт-кодону задается значением параметра "-find 1". Итак, полученная в результате команда:

   getorf -table 11 -minsize 30 -find 1 -sequence d89965.entret

   Получен файл с открытыми рамками, удовлетворяющими данным условиям. Согласно этому файлу 5-я рамка соответствует приведенной в записи CDS, а 13-я - записи Swiss-Prot (P0A7B8).



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

   Была дана последовательность всех тРНК, проаннотированных в полном геноме E. coli K12. С помощью программ BlastN и MegaBlast требовалось определить, сколько гомологов каждой из тРНК находится в геноме родственной бактерии Xanthomonas campestris. Для того, чтобы определить это число с помощью программы BlastN была введена следующая команда:

   blastall -p blastn -d xc -i trna_ecoli.fasta -m 8 -o trnahom.txt

   Был получен файл trnahom.txt со списком находок в виде таблицы. Для того, чтобы узнать теперь количество находок для каждой последовательности, в Excel был импортирован список названий входных последовательностей. А далее был создан скрипт, подсчитывающий количество находок для каждой последовательности. С помощью скрипта был получен файл blastn.txt. Результат был оформлен в таблице Excel (на странице names - имена последовательностей и число находок для каждой, на странице Лист1 - строчки скрипта для каждой последовательности).
   
   Чтобы получить теперь число находок с E-value, меньшим чем 0.001, была введена следующая команда:

   blastall -p blastn -d xc -i trna_ecoli.fasta -m 8 -e 0.001 -o trnahom2.txt

   В результате команды был получен файл с находками (в виде таблицы). Далее были проведены те же операции, что и при поиске находок без порога на E-value (соответствующий скрипт был сохранен в файле blastn2.scr, а число находок для каждой последовательности - в файле blastn2.txt). Результат - в той же таблице Excel (на тех же листах, что и при предыдущем поиске).
   
   Получим число гомологов с помощью программмы MegaBlast. Для этого выполним следующую команду:

   megablast -d xc -i trna_ecoli.fasta -m 8 -o mega.txt

   Полученный файл - mega.txt. Проведем те же операции. Скрипт сохранен в файле blastn3.scr, количество находок - в файле blastn3.txt.
   
   Для получения числа гомологов с помощью discontigous megablast введем команду:

   megablast -d xc -i trna_ecoli.fasta -m 8 -D 2 -t 18 -W 11 -N 1 -o mega2.txt

   В этой команде параметр -D задает тип выдачи (значение "2" задает стандартную выдачу blast), параметр -t задает длину слов из тРНК, которые будут искаться в геноме бактерии (может принимать значения "16", "18", "21"), параметр -W задает длину слов из генома бактерий, по которым ведется поиск последовательности (может принимать значения "11" и "12"), параметр -N задает тип разрывов в матрице (может принимать значения "0", "1" и "2"). Был получен файл mega2.txt. Скрипт для подсчета находок сохранен в файле blastn4.scr, файл с числом находок - в файле blastn4.txt.
   
   Для анализа была выбрана последовательность тРНК aspV. BlastN обнаружил находки для этой последовательности с E-value, меньшим 0.001, в отличие от MegaBlast, не обнаружившего для этой тРНК ни одной находки (MegaBlast не обнаружил ни одной находки, потому что эта программа ищет слова длиной в 28 букв в геноме бактерий, а BlastN - слова в 11 букв; как видно из выравнивания, совпадений в 28 букв подряд в нем нет, но есть совпадения в 11 букв, именно поэтому BlastN обнаружил находки, а MegaBlast - нет). Выравнивание BlastN последовательности aspU с одной из последовательностей генома Xanthomonas campestris сохранено в файле trnahomol.txt. Оно выглядит так:

   >AE012195 AE008922 Xanthomonas campestris pv. campestris str. ATCC
               33913,  section 103 of 460 of the complete genome.
             Length = 10083
   
    Score = 89.7 bits (45), Expect = 3e-19
    Identities = 69/77 (89%)
    Strand = Plus / Plus
   
                                                                           
   Query: 1    ggagcggtagttcagtcggttagaatacctgcctgtcacgcagggggtcgcgggttcgag 60
               |||||||||||||||  ||||||||| |  ||||||||||| || |||||||||||||||
   Sbjct: 9722 ggagcggtagttcagctggttagaatgctggcctgtcacgccggaggtcgcgggttcgag 9781
   
                                
   Query: 61   tcccgtccgttccgcca 77
               ||||||||| |||||||
   Sbjct: 9782 tcccgtccgctccgcca 9798
   

   Теперь выравним две последовательности программой needle. Для этого вначале вырежем соответствующие фрагменты командой seqret. Получаем файлы aspV.fasta и aspVxc.fasta. Теперь направим их программе needle командой:

   needle aspV.fasta aspVxc.fasta aspV.needle

   Полученное выравнивание сохранено в файле aspV.needle. Выравнивание выглядит так:

   # Aligned_sequences: 2
   # 1: aspV
   # 2: AE012195
   # Matrix: EDNAFULL
   # Gap_penalty: 10.0
   # Extend_penalty: 0.5
   #
   # Length: 77
   # Identity:      69/77 (89.6%)
   # Similarity:    69/77 (89.6%)
   # Gaps:           0/77 ( 0.0%)
   # Score: 313.0
   # 
   #
   #=======================================
   
   aspV               1 ggagcggtagttcagtcggttagaatacctgcctgtcacgcagggggtcg     50
                        |||||||||||||||..|||||||||.|..|||||||||||.||.|||||
   AE012195           1 ggagcggtagttcagctggttagaatgctggcctgtcacgccggaggtcg     50
   
   aspV              51 cgggttcgagtcccgtccgttccgcca     77
                        |||||||||||||||||||.|||||||
   AE012195          51 cgggttcgagtcccgtccgctccgcca     77
   

   Как видно из выравнивания, процент идентичности составляет 89,6%, что является очень высоким процентом. Процент гэпов, естественно, равен нулю. Такой высокий процент идентичности связан с консервативностью тРНК, так как она выполняет в клетках одни и те же незаменимые функции (для которых очень важна консервативность последовательности).
   Гомологичный участок проаннотирован в записи EMBL, описывающей геном бактерии, следующим образом:

   FT   gene            1135403..1135479
   FT                   /locus_tag="XCC0980"
   FT   tRNA            1135403..1135479
   FT                   /locus_tag="XCC0980"
   FT                   /product="tRNA-Asp"
   FT                   /note="Found by tRNAscan"
   

   Чтобы найти этот участок, пришлось зайти на сайт EBI и найти запись последовательности AE012195 (удаленную 27-го февраля 2009 года), а в ней следующие значения поля FT:

   FT   tRNA            9722..9798
   FT                   /gene="XCC0980"
   FT                   /product="tRNA-Asp"
   FT                   /note="Found by tRNAscan"
   

   Как видно, номера нуклеотидов в последовательности в точности совпадают с номерами в выравнивании. Ну а дальше по AC гена (XCC0980) была найдена эта последовательность в полной записи генома бактерии Xanthomonas campestris. Как и следовало ожидать, эта последовательность в ней соответствует той же самой тРНК (аспарагиновой), что и исходная тРНК из E. coli.



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

   Проведем поиск гомологов тРНК E. coli в геноме бактерии Xanthomonas campestris c помощью программы fasta35. Так как программа задумана как интерактивная, вначале создадим скрипт с ответами на вопросы, задаваемые программой (каждый ответ с новой строки). Запишем в нем (с помощью Excel) названия файлов, в которые попадут данные о гомологах каждой тРНК, после каждого названия поставим три пробела (нас удовлетворяют ответы, предлагаемые программой по умолчанию). В итоге получим скрипт. Теперь запустим программу fasta35 следующей командой:
   

   fasta35 trna_ecoli.fasta xc_genome.fasta 6 < fas.txt
   

   На выходе имеем файлы с набором гомологов для каждой тРНК.
   Теперь создадим скрипт, который сможет подсчитать число гомологов с E-value < 0.001 для каждого тРНК. Легче всего создать два скрипта. Первый подсчитает число гомологов с E-value < 0.0001 (записанный в виде [(число)e-(число)] ), а второй подсчитает число гомологов с E-value, лежащим от 0.0001 до 0.001 (записанный в виде [0.000(число)] ). Строчки для этих скриптов можно найти на странице Лист1 таблицы Excel. В результате работы скриптов были получены соответственно файлы fas1.txt и fas2.txt. В таблице Excel на странице names последние три столбца соответствуют числу находок программой Fasta с E-value < 0.0001; с E-value < 0.001, но > 0.0001; а так же сумма этих результатов (то есть число находок с E-value < 0.001). Колонка с конечным результатом окрашена в желтый цвет.
   Как видно, Fasta для каждого тРНК E. coli нашла в геноме бактерии Xanthomonas campestris меньше гомологов, чем BlastN. Более того, для всех гомологов, найденных обеими программами, E-value одной и той же находки при поиске в BlastN ниже, чем при поиске программой Fasta (например, гомолог тРНК valV AE012195 имеет при поиске BlastN E-value, равное 3e-19, а при поиске программой Fasta - 2.2e-11).