Monthly Archives: December 2013

Сборка генома de novo с использованием программы Velvet Using Velvet for de novo genome assembly

Представим, что Вы получили свои первые paired-end данные от Illumina и хотите провести первичную сборку генома в т.н. контиги (contigs). Рассмотрим сборку на примере программы Velvet. Прежде чем приступать к сборке, убедитесь в следующих моментах:

  1. У вас есть два (или больше, но четное число) FASTQ файла (ов) с forward и reverse последовательностями (ридами, reads – обычно именно в таком виде их и предоставляют секвенирующие компании);
  2. Вы удалили адаптеры или убедились, что их удалила компания, производившая секвенирование;
  3. Вы оставили в файлах только те последовательности, которые прошли проверку качества самого секвенатора иллюмины. Это определяется кодом Y в строке-заголовке каждого рида (reads): @EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG Важно! Y обозначает то, что последовательность НЕ прошла контроль качества секвенирующего прибора, N обозначает, что все в порядке. Если в файлах присутствуют записи с кодом Y их лучше удалить.
  4. Вы произвели всю предварительную чистку данных и сделали несколько вариантов очистки (различные стратегии чистки данных могут привести к различным результатам сборки генома)
  5. Вы проверили, что в файлах нет т.н. последовательностей-сирот (последовательностей без пары)
  6. Сама программа Velvet установлена на достаточно быстром и мощном сервере (либо размер данных невелик, например, вы работаете с бактериальным геномо). Ззапускать сборку генома человека на домашнем компьютере совершенно бессмысленное занятие.

Velvet производит сборку в два шага. Сначала запускается программа velveth, которая подготавливает данные для второй программы velvetg. Velveth создает т.н. хэш-данных, на основе которого затем будет создан и проанализирован граф де Брейна. Если вы используете систему на основе UNIX обе программы могут быть запущены одна за другой с помощью символа &&. Наиболее важный параметр любого сборщика геномов на основе графа де Брейна это размер кмера (kmer size). Заранее знать, какой кмер будет оптимален для ваших данных невозможно, им может быть любое значение от 50 до 150 (в среднем, бывают и исключения!). Обычно принято запускать несколько сборок с постепенно возрастающим размером кмера и затем сравнивать результаты и выбирать тот, который дает оптимальную сборку (определяется по значению N50).

Для программы Velvet размер кмера должен быть четным.

Вот как может выглядеть примерная команда-сборщик на основе Velvet для Mac/UNIX: velveth96 <output_directory> <-fastq> <-shortPaired> <-separate> <trimmed_input_r1.fastq> <trimmed_input_r2.fastq> && velvetg96 <output_directory> -ins_length <600> -ins_length_sd <100> -scaffolding -exp_cov -cov_cutoff -unused_reads -min_contig_lgth <300> Что значат указанные выше параметры: Все что указано в <> скобках зависит от ваших собственных данных и следовательно должно быть вписано осмысленно, а не просто следуя приведенному выше примеру.

  • output_directory – папка, в которую будут сохранены результаты сборки
  • kmer – размер кмера
  • -fastq – формат входных данных
  • -shortPaired – используются paired-end данные секвенирования
  • -separate – данные forward и reverse расположенные в двух отдельных файлах
  • trimmed_input_R1{2}.fastq – файл с очищенными последовательностями
  • ins_length размер вставки (т.н. insert size), определяется тем, как готовилась библиотека для секвенирования. Иногда информацию можно найти в документации, которая поступает с сырыми данными. Вам нужна величина за вычетом длины адаптеров!
  • ins_length_sd  – вариация в размере вставки (insert size). Плюс-минус указанное число к размеру ins_length
  • exp_cov auto – включает опцию расчета покрытия
  • cov_cutoff auto – включает опцию того, исключать ли из анализа участки слишком высоким или слишком малым покрытием
  • unused_reads yes – включает опцию сохранения не использованных ридов в отдельный файл
  • min_contig_lgth   – минимальная длина контига для сохранения в выходной файл

В следующей статье мы детально рассмотрим каждый из подчеркнутых выше терминов и выясним как рассчитывать объем компьютерной памяти (RAM) необходимой для запуска сборщика. Более детальное описание каждого параметра можно найти на сайте SEQwiki

Введение в программу RIntroduction to R (basics)

Большая часть задач в филогенетике, эмолюционной биологии и биоинформатике решается в настоящее время в программе R. Это бесплатная программа, состоящая из множества отдельных пакетов (модулей, packages), доступ к которым осуществляется через обший интерфейс. Все пакеты разрабатываются различными людьми и обычно каждый пакет имеет какую-то свою специализацию. Например, есть пакеты (ggplot2), который позволяют создавать красивые графики; есть специализированные пакеты для оптимизации статистических функций (optim, nlm); есть пакеты в которых можно создавать карты и анализировать пространственные данные (raster, spatial). Для нас же наиболее интересными являются пакеты для работы с филогенетическими деревьями. Ключевым таким пакетом является пакет ape.

Поддержка R осуществляется на официальном сайте CRAN и именно оттуда можно скачать дистрибутив программы для любой операционной системы. Также, этот сайт является основным репозиторием отдельных пакетов.

Установить программу очень легко. Для Windows достаточно скачать последний дистрибутив программы и следовать инструкции по установке программы.

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

Рабочее окно R

Рабочее окно R

R не имеет графического интерфейса и это значит, что любые функции и запросы вызываются с помощью введенных в командную стороку команд. Командная строка начинается с символа больше > и подсвечена красным цветом (как на картинке выше).

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

Команда в R работает как кнопка или меню в программах с графическим интерфейсом. Если команда принимает на вход некоторые параметры – то это полный аналог меню, где есть какие-то дополнительные условия для реализации задачи, если команда ничего не принимает на вход (нет параметров) – то тогда она работает как кнопка.

Параметры команды обычно указываются в скобках после названия команды. Например,

>plot(data)

Здесь plot – это название команды (ее еще можно называть функция), а data – это передаваемый параметр (то есть данные). Команда Plot позволяет визуализировать данные на экране.  Например, если data – это набор значений от 0 до 10, то вызов данной команды приведет к следующему результату:

Результат выполнения команды plot

Результат выполнения команды plot

Пакеты в R можно рассматривать как аналоги дополнительных модулей к пакетам с графическим интерфейсом. Например, в пакете SPSS есть ряд модулей для анализа временных серий. В R тоже есть аналогичные пакеты для работы с временными сериями.

Установить дополнительный пакет очень просто. Для этого в командной строке нужно написать (не забудьте кавычки):

>install.packages(“название_пакета”)

После этого будет предложено выбрать географическое расположение репозитория пакетов. Для России нужно выбрать Russia. После того, как пакет установлен, его загрузка в рабочее окно R осуществляется с помощью команды:

>require(название_пакета)

Обратите внимание, что во втором случае название пакета дается без кавычек.

Рабочее окно R – это текущий рабочий сеанс. После закрытия программы все загруженные пакеты, все вызванные команды, все результаты работы не сохраняются и будут потеряны. Чтобы это избежать необходимо сохранять и сам проект (рабочий сеанс) и историю вызыванных команд.

Основы Байесовского моделированияBasics of Bayesian modelling

Чтобы провести анализ лог-файлов получающихся на выходе из программ BEAST/MrBayes/BEAST* (да и любых других моделей на основе Байесовского моделирования) разберемся в базовых принципах Байесовского моделирования.

Модель в биологии

Очень часто модель в биологическом анализе  – это de facto некоторое распределение или совокупность распределений, форма которого (которых) определена набором параметров. Например, распределение морфологических признаков в популяции некоторого вида очень часто можно описать нормальным распределением (как на картинке внизу – цвет шкурки кроликов). Форма нормального распределения полностью описывается двумя параметрами – средним и дисперсией. Зная эти параметры мы можем проводить некоторые сравнительные тесты – например, узнавать различаются ли средние значения признака или их дисперсии между двумя популяциями.

pr33_1

Пример распределения фенотипического признак

Очевидно, что для расчета среднего и дисперсии нормального распределения существуют относительно простые формулы, позволяющие вычислить эти два параметра (прим. – в статистике эти параметры также называются моментами и о них говорят как о моментах случайной величины, которая имеет нормальное распределение). Метод расчета формы распределения по таким формулам называется аналитическим. Он хорош тем, что очень точен и занимает мало вычислительного времени. Тем не менее, в реальности большое число биологических моделей описываются намного более комплексными распределениями или произведением распределений друг на друга и найти аналитическое решение и вычислить по нему значения параметров может быть довольно трудно.

Альтернативный метод нахождения параметров распределения случайной величины называется численным, вычислительным (numerical). В отличие от аналитического, численный метод не требует знания формул моментов распределений. Не в даваясь в сложные математические детали того, как это делается, рассмотрим общий принцип численного метода – этот принцип пригодится нам и для понимания методов на основе теоремы Байеса.

Численный метод

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

12.47; 9.98; 9.61; 9.27; 8.83 12.50; 9.49; 5.40 12.64; 5.62

Введем еще два предположения: (1) предположим, что в теоретической популяции вся совокупность размеров околоцветников описывается нормальным распределением и (2) у нас нет аналитических формул для расчета среднего и дисперсии (разумеется, в реальности посчитать эти параметры очень легко, тем не менее представим что у нас нету формул).

Как работает численный метод?

Численные методы используют  алгоритм  т.н. максимального правдоподобия (maximum likelihood). Принцип его работы довольно прост: на основе выборки данных и предполагаемой модели (нормальное распределение) алгоритм случайным образом предлагает значения параметров (среднего и дисперсии). Затем, эти два значения параметров, а также сами наблюдаемые данные, подставляются в формулу функции правдоподобия. Чаще всего функция правдоподобия это логарифм фунцкии плотности вероятности некоторого теоретического распределения. К примеру, нормальное распределение имеет функцию плотности вероятности

Нормальное распределение

Нормальное распределение

Ее логарифм – это и есть фунция правдоподобия

pr33_4

Likelihood

Что происходит когда мы подставляем значения среднего (μ) и дисперсии (σ2), а также сами данные (xi) в эту формулу? На выходе мы получаем одно число, так называемую правдоподобность (обратите внимание, это не тоже самое что вероятность!) Эта правдоподобность (в английской нотации likelihood) говорит нам насколько вероятно получить такой набор данных (десять измерений выше), при подставленных значениях среднего и дисперсии и при заданной модели.

Затем алгоритм с использованием специального математического аппарата (чаще всего используются частные производные и  Гессиан-функция, но зависит от конкретного метода численной оптимизации) предлагает новые два значения среднего и дисперсии (или одно из них), опять вычисляет правдоподобность и проверяется как изменилось ее значение – если оно возросло, алгоритм считает, что он движется в правильном направлении и предложенные новые параметры более точно описывают данные. Соответственно алгоритм принимает их как текущие. Если же значение правдоподобности уменьшилась, то параметры эти считаются неверными, и алгоритм продолжает поиск от тех параметров, что были до этого. Такой поиск идет до тех пор, пока не будет найден так называемый пик функции правдоподобия – такое ее состояние когда изменение параметров в любом направлении будет приводить к уменьшению значения прадоподобности.  Визуально можно представить себе этот алгоритм как слепого горного туриста, которому нужно забраться на пик горы. На каждом шаге,он проверяет поднялся ли он выше по склону или спустился за счет звукового альтиметра.

Таким образом, метод максимального правдоподобия (maximum likelihood) позволяет вычислить условную вероятность параметров относительно  данных и модели. Чаще всего это записывают в виде:

P (данные| параметры {, модель})

Или  более формально – pr33_5

Модель взята в фигурные скобки потому что очень часто ее пропускают в написании.  Самое главное, что нужно запомнить, это то, что метод максимального правдоподобия позволяет получить правдоподобие (вероятность) данных при заданных параметерах.

Теорема Байеса

Теперь перейдем к теореме Байеса.  Основное отличие моделей в основе которых лежит теорема Байеса заключается в том, что теорема Байеса, позволяет вычислить вероятность модели (гипотезы) при условии наличия некоторого набора наблюдений. Еще раз акценитируем внимание на главной особенности теоремы Байеса – в отличие от метода максимального правдоподия, который вычисляет правдоподобие данных относительно заданных параметров, эта формула вычисляет вероятность модели относительно данных. В математическом выражении это записывают следующим образом:

pr33_6

Теорема Байеса

Обратим внимание на левую часть. В ней у нас есть:

pr33_7 – что напоминает данная простая формула? Она напоминает то, что мы рассмотрели чуть выше, а именно формулу максимального правдопобия. Именно ей она и является! Это обычная формула правдоподобия данных относительно параметров

pr33_8 – это так называемое априорное распределение параметра модели ( если параметров много, то каждый параметр будет иметь свое априорное распределение)

pr33_9– это т.н. полная вероятность. Полная вероятность это произведение вероятностей данных относительно всех возможных значений параметров. То есть как если бы мы брали формулу максимального правдоподобия (как в числителе) и подставляли туда один за другим все возможные значения параметров, а затем перемножили их все между собой.

Сделаем небольшое заключение. Теорема Байеса позволяет вычислить вероятность модели относительно некоторых данных и параметров. В этом ее принципиальное отличие от метода максимального правдоподобия, который вычисляет правдоподобие параметров относительно данных и модели.