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

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

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

Для поиска гомологов я пользовался программой tblastn пакета BLAST и командами:
    formatdb -i xc_genome.fasta -n index -p F
    blastall -p tblastn -d index -i P0ABH7.fasta -e 0.001 > P0ABH7.out
На выходе получился файл P0ABH7.out

Число находок с Е-value<0,001 2
Характеристика лучшей находки:  
   E-value находки e-160
Название геномной последовательности AE012440 AE008922 Xanthomonas campestris pv. campestris str. ATCC 33913, section 348 of 460 of the complete genome.
Координаты выравнивания(-ий) в найденной последовательности 10423-9152
Обе находки располагаются на одной геномной последовательности.

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

Найденная в предыдущем задании последовательность была вырезана в файл: 
    seqret -sask
    Reads and writes (returns) sequences
    Input (gapped) sequence(s): xc_genome.fasta:AE012440
         Begin at position [start]: 9152
           End at position [end]: 10423
            Reverse strand [N]: Y
    output sequence(s) [ae012440.fasta]:
Выходной файл - ae012440.
Далее был запущен поиск этой последовательности в банке "EMBL standard prokaryote".
У первой находки (AM920689) был выбран режим "Show Alignments". Получилось следущее выравнивание. В записи AM920689 последовательность имеет координаты 1090364 - 1091635. Далее была получена запись AM920689. Для этого я воспользовался командой: 
entret embl:AM920689 -auto
Информация о заданном участке в поле FT: 
FT   CDS             1090280..1091638
FT                   /transl_table=11
FT                   /gene="gltA"
FT                   /locus_tag="xcc-b100_0951"
FT                   /product="citrate (Si)-synthase"
FT                   /function="Citrate synthase"
FT                   /EC_number="2.3.3.1"
FT                   /db_xref="GOA:B0RPB4"
FT                   /db_xref="InterPro:IPR019810"
FT                   /db_xref="UniProtKB/TrEMBL:B0RPB4"
FT                   /protein_id="CAP50299.1"
FT                   /translation="MPWTSITCLPMWDQAPSTEGAYTVSDLDQVTLNAGDKSVVLPVLK
FT                   PTLGNDCVDISKLTKETGLFTYDSGFTATASCKSAITYIDGDNGVLLYRGYPIEQLAEK
FT                   SSFLEVSYLLMNGELPTADEFKKFDHEVTHHTMMHESLKNFLGGFRHDAHPMAMLAGSV
FT                   ASLSAFYHDTLDLNDPEQRRQAAIRLIAKVPTLAAAAYRYSIGWPIRYPRNNLNYVDRF
FT                   LHMMFEVPSEPLEINPVVAKALDLLFILHADHEQNASTSTVRLVGSTGANPYASVAAGI
FT                   TALWGPAHGGANEAVLKMLEEIGTADNVESAVAKAKDKNSSFRLMGFGHRVYKNFDPRA
FT                   KIIREMTHKVLGELGVNDPLLEVALKLEEAALKDDYFVQRKLYPNVDFYSGLIYKALNI
FT                   PVEMFTVMFAIARTAGWVSHWLEQQVDPEMKIGRPRQIYTGYDKRDYTDAGKR"
Координаты CDS: 1090280..1091638
Участок соответствует записи B0RPB4 UniProt

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

Выполняя прошлое задание, мы получили список записей EMBL с последовательностью, кодирующей белок CISY_ECOLI из генома E.coli. Из списка была выбрана запись V01501 EMBL. Эта запись была сохранена в отделный файл с помощью команды entret: 
entret embl:V01501 -auto
В записи содержались координаты CDS: 460..>755
Также в записи содержится информация о том, что участок GATGGTGAAA (с 732 по 744 нуклеотид) был изменён на GAATGGTGAAAAA.
Кодирующая последовательность была вырезана в отдельный файл: 
seqret "embl:V01501[460:755]"
Был получен файл с последовательностью. Далее программой BLASTN был проведён поиск гомологов этого гена в геноме Xanthomonas campestris: 
  blastall -p blastn -d index -i v01501.fasta > v01501.out
Был получен файл v01501.out
E-value лучшей находки 0.073, BLASTN не выдал ни одной находки с E-value < 0.001. Лучшая находка:
AE012258 AE008922 Xanthomonas campestris pv. campestris str. ATC... 34 0.073

  Аминокислотная последовательность Нуклеотидная последовательность
Число находок с Е-value<0.001 2 0
Характеристика лучших находок:
   E-value находки e-160 0.073
Длина находки 1272 17
Процент совпадений 60 100

Очевидно, поиск по аминокислотной последовательности намного лучше. При поиске по нуклеотидной последовательности находятся только короткие участки, причём такие же участки находятся и в негомологичных последовательностях. Это связано с тем, что гены кодируются всего четырьмя буквами. Кроме того, код ДНК вырожден, т.е. большая часть аминокислот представлена более чем одним кодоном, поэтому при поиске по нуклеотидной последовательности, многие гомологи "теряются".

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

Был создан файл D89965.entret с записью банка EMBL: 
entret embl:d89965 -auto
Была запущена программа getorf со следующими параметрами: Для этого была выполнена следущая команда: 
getorf -minsize 30 -find 1 -table 11
Получили файл d89965.orf 5 рамка соответствует CDS, а 13 - записи Swiss-Prot (P0A7B8).

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

Определим, сколько гомологов каждой из тРНК находит программа BLASTN в геноме родственной бактерии. Для этого запустим blastn с выводом в виде таблицы. 
blastall -p blastn -d index -i trna_ecoli.fasta -m 8 > t1.txt
Теперь повторим то же самое, но поставим порог E-value < 0.001 
blastall -p blastn -d index -i trna_ecoli.fasta -m 8 -e 0.001 > t2.txt
В результате получены 2 таблицы: t1.txt без ограничений и t2.txt с E-value < 0.001
Теперь напишем скрипт из команд, выдающих число находок для каждой последовательности.
  • Получим список всех тРНК
    grep ">" trna_ecoli.fasta
  • Копируем его в Excel и создаем таблицу, воспользовавшись формулой
    =СЦЕПИТЬ("grep -c '";B3;">' t1.txt >> list.txt")
    В получившейся таблице 2 столбца: один с названиями всех тРНК, второй - с текстом команд, которые будут использоваться в нашем скрипте.
    Ссылка на таблицу: trna1.xls
  • Копируем заполненную часть столбца B в редактор Far manager (перед выходом нажимаем <Shift+F2> и выбираем в меню "Unix format").
  • Теперь надо сделать файл со скриптом исполняемым: 

     chmod +x trna_linux.scr
    • Скрипт готов. Запускаем: 
     ./trna_linux.scr
    В результате работы скрипта был получен файл list.txt, содержащий только количества находок.
    Аналогично напишем второй скрипт и получим данные с ограничением порога e-value 0.001:
    скрипт и файл с количеством находок.

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

    Повторим предыдущий опыт с использованием программ megablast и discontigous megablast.
    Megablast: 
     megablast -d index -i trna_ecoli.fasta -m 8 > trna_megablast.fasta
    скрипт для Megablast
    Discontigous megablast: 
    megablast -d index -i trna_ecoli.fasta -m 8 -D 2 -t 18 -W 11 -N 1 > trna_discontigous.fasta
    скрипт для Discontigous megablast
    -D - определяет вид выходных файлов (2 - стандартный для программ пакета BLAST)
    -W 11 -t 16 -N 1 - определяют одну из форм запуска discontigous megablast с длиной слов 11, по некодирующей части генома.
    -o -d -i - тоже самое, что и для blastall (выходной файл, имя индексных файлов и входной файл соответственно)
    Результаты предыдущих заданий в виде таблицы

    7. Анализ результатов

    Рассмотрим тРНК valW. В выдаче BLASTN были находки из генома бактерии Xanthomonas campestris В выдаче megablast - нет.
    Выберем одну из гомологичных последовательностей бактерии Xanthomonas campestris и вырежем её в отдельный файл: 
    seqret -sask
Reads and writes (returns) sequences
Input (gapped) sequence(s): xc_genome.fasta:ae012274
     Begin at position [start]: 10086
       End at position [end]: 10106
        Reverse strand [N]: y
output sequence(s) [ae012274.fasta]: valW_xc.fasta
    Затем вырежем в отдельный файл тРНК бактерии. Построим выравнивание программой needle: 
    needle valW.fasta valW_xc.fasta valW.needle
    Характеристики выравнивания: 
    Length: 77
Identity:      20/77 (26.0%)
Similarity:    20/77 (26.0%)
Gaps:          56/77 (72.7%)
Score: 96.0
    Программа BLASTN ищет по словам длины 11, а megablast - 28. Поскольку в выравнивании нет 28 совпавших нуклеотидов, megablast не нашёл этот участок.
    Данный фрагмент не проаннотирован. В EMBL найти его не удалось.

    Назад