EMBOSS. Выравнивание геномов.


Упражнения EMBOSS

  1. (seqret) Несколько файлов в формате fasta собрать в единый файл.

    Были использованы файлы HSP7C_HUMAN.fasta, TERT_HUMAN.fasta и CISY_HUMAN.fasta из предыдущего практикума.
    В этом случае было удобно собрать файлы по маске.

    Использованная команда: seqret "*_HUMAN.fasta" HUMAN.fasta
    Выходной файл: HUMAN.fasta

  2. (seqretsplit) Один файл в формате fasta с несколькими последовательностями разделить на отдельные fasta файлы

    Был взят файл HUMAN.fasta из предыдущего задания.

    Использованная команда: seqretsplit HUMAN.fasta
    Выходные файлы: p11142.1.fasta, o14746.1.fasta, o75390.2.fasta

  3. (seqret) Из файла с хромосомой в формате .gb вырезать три кодирующих последовательности по указанным координатам "от", "до", "ориентация" и сохранить в одном fasta файле

    Я взяла последовательность полного генома Streptococcus suis GZ1. (файл CP000837.1) и выбрала отттуда три кодирующих послдовательности. Для этого сначала был создан список с USA интересующих нас участков с указанием координат и, где требовалось, ориентации "r". Затем выбранные фрагменты последовательности были объединены в один fasta-файл.

    Использованные команды:

    echo -e "gb::genbank:CP000837[1528:2664]\ngb::genbank:CP000837[3647:3850]\ngb::genbank:CP000837 [57196:58134:r]">cds.list

    seqret @cds.list fasta:cds.fasta

    Выходной файл: cds.fasta

  4. (transeq) Транслировать кодирующие последовательности, лежащие в одном fasta файле, в аминокислотные, используя указанную таблицу генетического кода. Результат - в одном fasta файле.

    Полученные в предыдущем задании кодирующие последовательности были транслированы в последовательности аминокислот с использованием стандартной таблицы генетического кода (параметр -table 0, стоит по умолчанию). Стоит заметить, что все они начинались стартовым кодоном, который был транслирован в метионин, и заканчивались стоп-кодоном, который в выходном файле обозначен "*". На выходе было получено столько же последовательностей, сколько было во входном файле.

    Использованная команда: transeq -table 0 cds.fasta proteins.fasta

    Выходной файл: proteins.fasta

  5. (transeq) Транслировать данную нуклеотидную последовательность в шести рамках.

    Чтобы транслировать последовательность в шести рамках, необходимо было добавить опцию -frame 6. В выходном файле оказалось 18 последовательностей (каждая из трех исходных последовательностей была транслирована в 6 рамках, из которых 3 прямых и 3 обратных).

    Использованная команда: transeq -frame 6 cds.fasta proteins6.fasta

    Выходной файл: proteins6.fasta

  6. (seqret) Перевести выравнивание из fasta формата в формат .msf

    Исходный файл - alignment.fasta. После перевода в msf формат получен файл alignment.msf, в котором, в отличие от исходного fasta файла, приведено выравнивание всех последовательностей вместе, а не по отдельности, а также содержится дополнительная информация о последовательностях.

    Использованная команда: seqret alignment.fasta msf::alignment.msf

    Выходной файл: alignment.msf

  7. (infoalign) Выдать в выходной поток число совпадающих букв между второй последовательностью выравнивания и всеми остальными (на выходе только имя последовательности и число).

    Infoalign выдает различную информацию о последовательностях во входном множественном выравнивании (USA, имя, длину, количество гэпов, совпадений и пр.) в сравнении с эталонной последовательностью. По умолчанию эталонной последовательностью является вычисленная консенсусная последовательность, но ее можно задать и вручную по имени или порядковому номеру в файле. Нам необходимо было провести сравнение со второй последовательностью, поэтому была использована опция -refseq 2. Чтобы получить на выходе только имя и информацию о количестве совпадений были использованы опции -only -name -idcount.

    Использованная команда: infoalign alignment.msf -refseq 2 -only -name -idcount stdout



  8. (featcopy) Перевести аннотации особенностей в записи формата .gb в табличный формат .gff

    Featcopy читает таблицы особенностей и переводит их в любой из поддерживаемых форматов. Исходный файл - sequence.gb.

    Использованная команда: featcopy sequence.gb sequence.gff

    Выходной файл: sequence.gff

  9. (extractfeat) Из данного файла в формате .gb получить fasta файл с кодирующими последовательностями; (*) добавить в описание каждой последовательности функцию белка (из поля product)

    Extractfeature предназначена для извлечения участков последовательности, аннотированных какой-дибо особенностью (-type). Участки с одной и той же особенностью могут быть извлечены как отдельные последовательности или сцеплены вместе, если используется опция -join. Нам необходимо было получить файл с кодирующими последовательностями, поэтому команда была запущена с опцией -type CDS. Чтобы добавить в описание функцию белка была также использована опция -describe. Исходный файл - sequence.gb

    Использованная команда: extractfeat sequence.gb -type CDS -describe product coding.fasta

    Выходной файл: coding.fasta



  10. (shuffle) Перемешать буквы в данной нуклеотидной последовательности; (*) проверить с помощью blastn сколько "достоверных" находок (с E-value < 0.1) найдется в нуклеотидном банке данных (запустите с порогом E = 10 - по умолчанию)

    Shuffle прочитывает одну или несколько последовательностей и случайным образом перемешивает их. Число перемешиваний может быть задано. Чтобы результат выводился в файл, необходимо использовать опцию -о. Для работы был взят файл cds.fasta.

    Использованная команда: shuffle -o shuffled.fasta cds.fasta

    Выходной файл: shuffled.fasta

    Для выходного файла был запущен blastn с параметрами по умолчанию. "Достоверных" находок (с E-value < 0.1) не нашлось ни одной. Выдачу blastn можно увидеть на рисунке.



  11. (cusp)Найдите частоты кодонов в данных кодирующих последовательностях.

    Cusp высчитывает таблицу использования кодонов для одной или более нуклеотидных последовательностей. В данной таблице для каждого кодона указано: последовательность кодона, закодированная аминокислота, частота данного кодона по отношению к остальным, кодирующим ту же аминокислоту (fraction), частота кодона в данной последовательности (frequency), число кодонов данного типа в последовательности.
    Для работы был взят файл cds.fasta.

    Использованная команда: cusp cds.fasta frequency.cusp

    Выходной файл: frequency.cusp

  12. (compseq) Найдите частоты динуклеотидов в данной нуклеотидной последовательности и сравните их с ожидаемыми

    Compseq производит подсчет состава слов указанной длины во входной последовательности. В выходной файл записывается само слово, сколько раз оно встречается, его наблюдаемая частота встречаемости (Obs Frequency), предполагаемая частота (Exp Frequency), а также соотношение этих частот (Obs/Exp Frequency). Чтобы найти наблюдаемые и ожидаемые частоты динуклеотидов необходимо запустить данную команду для слов длины 2 (опция -word 2). Для более точного подсчета ожидаемых частот можно использовать опцию -calcfreq, при активировании которой данный подсчет ведется на основании частот встречаемости отдельных оснований в данной последовательности (если опция выключена, то частоты всех оснований считаются равными, что не является верным). Исходный файл - cds.fasta.

    Использованная команда: compseq -word 2 -calcfreq cds.fasta cds.compseq

    Выходной файл: cds.compseq

    Практически для всех динуклетидов наблюдаемая частота встречаемости не особо отличается от ожидаемой, соотношение Obs/Exp Frequency для всех, кроме двух динуклеотидов, отклоняется от единицы не более, чем на 0,2. Для динуклеотидов CG и TA оно составляет 0.7695578 и 0.7000311 соответственно. Это значит, что данные динуклеотиды встречаются в последовательности несколько реже, чем было предсказано.

  13. (tranalign)Выровняйте кодирующие последовательности соответственно выравниванию белков - их продуктов

    Tranalign принимает на вход набор из невыровненных нуклеотидных последовательностей и соответсвующий им набор выровненных транслированных с них белковых последовательностей. В выходной файл записывается выравнивание нуклеотидных последовательностей. Каждая нуклеотидная последовательность транслируется во всех трех прямых рамках по указанному генетическому коду и трансляции сравниваются с соответсвующими во входном выравнивании. Важно, чтобы соответствующие друг друг последовательности во входных файлах располагались в одном порядке. Для работы были взяты следующие файлы: nuc.fasta - невыровненные нуклеотидные последовательности, pep.fasta - выровненные белковые последовательности.

    Использованная команда: tranalign nuc.fasta pep.fasta nucalign.fasta

    Выходной файл: nucalign.fasta

Выравнивание геномов


Построение карты локального сходства

Для выполнения задания были взяты бактерии Yersinia pestis Nepal516 и Yersinia pestis Antiqua. Это патогенные бактерии, вызывающие чуму.

С помощью blast2seq было построено выравнивание и получена карта локального сходства для этих двух штаммов.

Карта локального сходства Yersinia pestis Nepal516 (OX) и Yersinia pestis Antiqua (OY)

Анализ полученнной карты

Данная карта (dot matrix view) показывает участки сходства на основании результатов работы BLAST. Выравнивания показаны в виде непрерывных линий. Совпадения прямой цепи (plus strand) отображены линиями, идущими из левого нижнего в правый верхний угол, совпадения комплементарной цепи (minus strand) - из левого верхнего в правый нижний. Количество линий соответсвует количеству найденных BLAST выравниваний.
Последовательность на оси OY (Query) имеет длину 4702289, полседовательность на оси OX (Reference) - 4534590. Минимальная цена деления на осях составляет 50000.

Cходство (Identity %) между гомологичными участками в данном выравнивании - 99% (было найдено как среднее сходство по нескольким наиболее длинным выравниваниям).

Синтеничные области
Синтеничные области - участки геномов, состоящие из ортологичных областей с сохранением их порядка на хромосоме для сравниваемых геномов. Наиболее крупные синтеничные области для геномов Yersinia pestis Nepal516 (OX) и Yersinia pestis Antiqua (OY) подчеркнуты зеленым.
Наиболее крупные эволюционные события
Рис. 1
Красным выделены инверсии.
Рис. 2
Желтым выделены вставки в Query
Рис. 3
В выделенной голубым области вероятно произошла транслокация. Наблюдается перестановка порядка участков - AB в одном геноме, BA - в другом.

Построение нуклеотидного пангенома

Чтобы построить нуклеотидный пангеном, я использовала программу NPG explorer на сервере kodomo.

В домашней директории была создана специальная папка ann_karpukhina@kodomo:~/term3/block2/pr9/npg1.
В ней я создала файл genomes.tsv, содержащий адреса последовательностей геномных ДНК для выбранных мной штаммов бактерий: Yersinia pestis Java9, Yersinia pestis Nairobi и Yersinia pestis Shasta.

Затем я запустила команду npge Prepare, чтобы скачать и переименовать геномные ДНК. В результате были получены файлы: genomes-raw.fasta с исходными последовательностями, genomes-renamed.fasta с переименованными последовательностями, features.embl с особенностями (features) выбранных организмов и genes/features.bs. с блоками генов. Файлы genomes-raw.fasta и features.embl содержат необработанную входную информацию и на последующих этапах не используются, поэтому их можно удалить.

С помощью команды npge Examine была проведена оценка сходства геномов и создан файл identity_recommended.txt, на основании которого я определила рекомендуемый для моих последовательностей параметр MIN_IDENTITY - 0.900.

Далее я создала файл с параметрами командой npge -g npge.conf и изменила в нем необходимые пункты (WORKERS = 1 для kodomo; MIN_IDENTITY менять не пришлось, так как он совпал со значением по умолчанию).

Наконец была вызвана команада npge MakePangenome, которая создала файл pangenome.bs - необходимый нам пангеном в формате BlockSet. Для дальнейшего анализа использовалась команда npge PostProcessing, которая создала множество файлов с аналитической информацией о пангеноме.