Анализ ортологов в геномах трех бактерий Ralstonia с помощью OrthoMCL

By | 10 July 2015

Один из наиболее мощных инструментов для кластеризации ортологов – это программа OrthoMCL. Несмотря на свою популярность, это также одна из наиболее трудно устанавливаемых и настраиваемых программ в биоинформатике. Мы постарались максимально упростить работу с этим инструментом в рамках InsideDNA и описываем ниже полный пайплайн для анализа ортологов в OrthoMCL / MCL .

Использование сайта InsideDNA бесплатно и открыто для всех!

1. Входные данные

Мы будем обрабатывать геномы для трех видов бактерий рода Ralstonia полученных с ftp NCBI. Все три генома можно скачать в одном общем архиве вот здесь (792Kb). Исходные данные для OrthoMCL могу быть представлены как в виде fasta файлов с аминокислотными последовательностями, так и нуклеотидными. Кроме того, это могут быть как аннотированные данные, так и просто результаты секвенирования RNA-seq или single reads/paired end данных Illumina. Чаще всего анализируются данные по нескольким геномам.

2. Сборка всех нужных инструментов в проект

Залогиньтесь в приложение InsideDNA и перейдите во вкладку Tools. На первом этапе мы создадим проект Bacterial_orthomcl

 1

и добавим туда все необходимые для анализа инструменты (7):  mcl,  orthomclAdjustFasta, orthomclFilterFasta, orthomclBlastParser, orthomclLoadBlast, blastall,  formatdb.

2

3. Запуск orthomclAdjustFasta

Первый этап  – это подготовка входных fasta-файлов. На этом этапе мы форматируем строку названия последовательности. Например, в нашем случае мы приводим заголовки каждого из трех входных fasta файлов из вида:

к виду:

>RE|NP_942643.1

Где RE – это краткое название вида (например, Ralstonia_eutropha) Процедуру коррекции заголовков нужно повторить для всех входных fasta файлов (например, для всех видов бактерий, в которых мы хотим найти ортологи). Мы используем вот такие параметры для вида Ralstonia_solanacearum:

 3

В нашем случае, поле с уникальным номер последовательности находится на четвертой позиции – но в случае других данных, позиция может отличаться. Нам нужно повторить процедуру для всех трех геномов

 4

В результате мы получим три откорректированных файла

 5

  Вот так выглядит fasta-файл до коррекции

 6

И после

 7

Переместим данные в отдельную папку adjusted

 8

4. Запуск orthomclFilterFasta

На следующем шаге отфильтруем все последовательности во всех геномах по длине и числу стоп-кодонов (минимальная длина протеина – 20 АА и максимальный процент стоп-кодонов – 20%) и объеденим их в один общий файл. 9

Именно для этой команды мы перемещали все откорректированные геномы в отдельную папку.

10

На выходе мы получим единый файл со всеми данными из всех трех геномов который называется goodProteins.fasta. Также мы получим файл poorProteins.fasta с отфильтрованными (плохими) последовательностями.

 11

5. Запуск formatdb

Создадим новую папку db

 12

И переместим файлы goodProteins.fasta и poorProteins.fasta в эту папку

 13

Теперь запустим форматирование данных для reciprocal blast с помощью команды formatdb. Эта команда трансформирует фаста-файл в специальный формат базы данных BLAST.

 14

На выходе мы получаем три дополнительных файла – phr, pin, psd

 15

5. Запуск Blastall

Теперь запустим одну из наиболее важных команд, которая и выполняет reciprocal blast (one-to-one blast – бласт один к одиному) и проверяет парное сходство между всеми последовательностями во всех геномах. Откроем окно параметров Blastall. Во-первых, нам нужно указать нашу пользовательскую базу BLAST – ту самую, что мы форматировали с formatdb. Приставка этой базы – goodProteins.fasta Во-вторых, мы должны передать серверу все входные файлы относящиеся к базе – phr, pin, psd, fasta – выбираем их в третьем поле сверху В-третьих, мы задаем т.н. query file – это собственно говоря тот же самый goodProteins.fasta – именно потому, что мы шлем запрос с данными к базе данных на основе этого же файла процедура называется reciprocal blast.

 16

В-четвертых, мы задаем выходной файл – blastallvall_prot.csv В-пятых, наши данные – это аминокислотные последовательности – следовательно мы используем алгоритм Blastp В-шестых, мы можем указать e-value – пороговое значение по которому будут определяться парные последовательности. В-седьмых, важно не забыть указать формат выходного файла – нам нужен тип 8 – tabulated format. Именно он может быть корректно прочитан программой orthoMCLBlastParser на следующем этапе.

 17

Для blastall лучше указать очередь помощнее (8 ядер и 30 Gb RAM)

 18

В результате обработки программы будет создан blastallvall_prot.csv с содержимым вида:

 19

6. Запуск orthoMCLBlastParser

После того, как мы получили файл blastallvall_prot.csv мы должны его обработать программой OrthoMCLBlastParser. Для этой программы нам потребуется как файл blastallvall_prot.csv, так и директория adjusted с тремя откорректированными фаста-файлами с геномами бактерий. Назовем выходной файл similarSequences.txt.

220

Это программа в результате создаст файл similarSequences.txt в папке blast_parser со следующим содержимым

221

7. Запуск orthoMCLLoadBlast

Теперь мы готовы к финальному шагу с использованием OrthoMCLLoadBlast – на вход нам надо предоставить файл similarSequences.txt. На вход этой программе нам надо передать всего лишь один файл, а также указать название лог-файла. Остальные выходные файлы будут автоматически созданы в той же директории, где расположен файл similarSequences.txt (то есть в нашем случае –  в папке blast_parser)

222

В результате работы этой программы будут созданы следующие выходные файлы и папки

223

Внутри папки pairs находятся следующие файлы

224

Нас же интересует текстовый файл mclInput – который и требуется для финально кластеризации ортологов с помощью программы MCL

225

8. Запуск MCL

Финальным шагом анализа ортологов в OrthoMCL является запуск программы MCL. Для кластеризации нам необходимо указать созданный mclInput , а также выбрать пороговое значение для кластеризации на графе – параметр -I, который мы ставим равным 1.5. Важно понимать, что с подбором данного параметра нужно экспериментировать и его оптимальное значение будет определяться тем, насколько эволюционно дивергентными являеются виды.

226

Результатом обрабоки данного файла является еще один файл с разделителями табуляция, где каждая строка представляет собой список генов, относящихся к общей ортологической группе

227

Всего, в случае с тремя геномами бактерий Ralstonia мы получили 562 группы ортологов. Обычно т.н. conserved genes – то есть гены представляющие core genome – представлены однократно в каждом геноме и расположены в нижней части csv файла, в то время как псевдогены (к примеру, ретротранспозоны) будут расположеные в верхних строчках csv файла и будут представлены общирными группами с сотнями генов. И те и другие группы ортологов интересны – например, core genome genes могут быть использованы для реконсрукции консенсусного дерева, в то время как ретротранспозоны – для анализа истории дупликаций или горизонтального переноса генов.