Программы пакета 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).

Назад