Monthly Archives: February 2014

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

Среди сравнительных методов одним из наиболее старых является метод реконструкции предковых дискретных состояний (ancestral character state reconstruction) с использованием филогенетических деревьев.  Исторически, в основе метода  лежали две математических модели:

  • Модель бережливости (parsimony, Maddison et al., 1984)
  • Модель на основе непрерывного марковского процесса (continuous-time Markov process, Pagel , 1994)

В настоящее время оценка параметров моделей на основе непрерывного марковского процесса реализована с помощью:

  • Метода максимальнго правдоподобия (maximum likelihood, Pagel, 1999)
  • Метода на основе теоремы Байеса и Монте-Карло с цепями Маркова (Bayesian MCMC, Pagel, 2004)

Кратко вспомним что такое модели на основе непрерывного марковского процесса (НМП).

НМП – это такой процесс, в котором состояния Марковской цепи (пространство состояний цепи) принимают значения конечного или счетного множества и для которой время проведенное в каждом из этих состояний распределено экспоненциально. Переводя на чуть более простой язык: модель НМП подразумевает, что число конечных состояний цепи (например, дискретных биологических признаков) – счетно (то есть мы можем назвать сколько есть состояний есть у нашего объекта заранее), а вероятность  пребывания в каждом из состояний по времени имеет экспоненциальное распределение.

pr40_1

Экспоненциальное распределение (waiting time)

Марковское свойство данного процесса заключается в том, что любое следующее состояние процесса полностью определеяется его текущим состоянием – таким образом, процесс этот является процессом без памяти.

Математически НМП задается:

  1. пространством состояний цепи S (states set)
  2. матрицей мгновенных переходных вероятностей Q (transition rate matrix)
  3. начальное вероятностное распределение признаков цепи Маркова (initial probability distribution for a state space)

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

Попробуем перевести вышесказанное на биологический язык. Проще всего этого сделать, если рассмотреть пример элементарного филогенетического дерева из двух ветвей. Пусть дерево будет ультраметричное и длина каждой ветви составляет 5 единиц (например, миллионов лет). На концах дерева зададим два дискретных состояния бинарного признака [0;1] – эти два состояния формируют пространство состояний цепи S. Пусть в реальности 0 кодирует белый окрас, а 1 – серый окрас шкурки мыши.

Простейший пример для предковой реконструкции

Простейший пример для предковой реконструкции

Наша задача – определить вероятности признаков 0 или 1 для основания дерева (то есть какой был цвет шкурки в корне дерева).

Для того, чтобы посчитать вероятность состояний 0 или 1 для корня дерева, нам необходимо знать вероятности переходов  P01 и P10  с  течением времени  – именно это задается матрицей переходных вероятностей Q и временем t. Для непрерывного марковского процесса  вероятность состояния расчитывается по формуле (всю математику можно найти в уравнении Колмогорова):

pr40_3

Где t – время, которое длится случайный процесс  (определяется длиной ветвей дерева, в нашем случае равно 5), а Q – это матрица мгновенных переходов (instantaneous rate matrix):

pr40_4

q01 и q10 – это параметры мгновенного перехода из состояния 0 в 1, и обратно. Обратите внимание, что согласно свойствам Марковских цепей сумма величин по строкам должна равняться нулю. Таким образом, легко сделать вывод, что -q01 это просто q01 с отрицательным знаком. Перемножение на время t и потенциирование (возведение в степень) матрицы Q дают на выходе вероятности переходов P01(t) и P10(t) . Обратите внимание, что P00(t) и P11(t) равны (1 – P01(t)) и (1 – P10(t)) соответственно:

pr40_5

Важно: в нашем примере с мышками параметры матрицы мгновенных переходов Q неизвестны (!) и должны быть оценены с помощью метода максимального правдоподобия или других численных методов.  Чтобы исключить шаг из рассмотрения, для простоты мы просто предположим что q01 = 0.4 и q10 = 0.9. Теперь мы легко можем подсчитать вероятности переходов P01(t) и P10(t),  подставив параметры q01 = 0.4 и q10 = 0.9 и t=5 в формулу выше.

На следующем шаге мы можем рассчитать полную вероятность P(t) для наших данных d при заданных параметрах (m= [q01 = 0.4, q10=0.9,  t=5]) по формуле:

pr40_6

Где w(0) и w(1) – начальное вероятностное распределение признаков цепи – то есть начальные вероятности состояний анцестрального признака (a) – initial probability distribution for a state space. В нашем случае, предположим, что и w(0) и w(1) равны 0.5, ведь оснований полагать различные начальные вероятности для корневого признака у нас нет.

Сделаем новое допущение – представим, что нас интересует не вероятности как таковые, а правдоподобие модели.  Решим (согласно Edwards, 1972), что что правдоподобие пропорционально (с точностью до константы) вероятности, тогда:

pr40_7

Кроме того, теперь мы можем посчитать и условную вероятность того, что в основании дерева была белая или серая мышка по формуле:

pr40_8

Где i={0;1}. Сравнив  два полученных праводоподобия, мы можем сделать вы вывод о том, какое анцестральное состояние было более вероятным.

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

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

pr40_9

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

Расшифровка терминов часто используемых при сборке генома

Начнем с термина использованного в заголовке статьи.

Сборка генома (genome assembly) – процесс создания генома из большого числа коротких нуклеотидных последовательностей (ридов) длинной от 50 нк (нуклеотидов) до нескольких тысяч. Обычно сборка генома включает в себя ряд этапов:

  1. Первичная обработка и чистка данных
  2. Сборка ридов в контиги (contings)
  3. Сборка контигов в скаффолды (scaffolds)
  4. Секвенирование и закрытие дыр (гэпов – от англ. gaps)

Каждый из перечисленных этапов осуществляется многократно и часто в несколько итераций в зависимости от результатов следующего или текущего шага.

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

Рассмотрим пример рида полученного методом секвенирования paired-end и структуру блоков, которые его составляют.

Размер вставки и размер фрагмента

Размер вставки и размер фрагмента

Размер вставки (insert size) – часть нуклеотидной последовательности между двумя парными адаптарами (paired-end adaptors).

Размер фрагмента (fragment size) –  часть нуклеотидной последовательности с учетом длины адапторов на концах.

Внутреннее парное расстояние (mate inner distance) – внутреннее расстояние между концами двух ридов.

Таким образом, важно помнить, что размер вставки обычно включает в себя длины ридов и расстояние между ними. Это расстояние зависит от метода секвенирования и его нужно уточнять у компании или коллег, которые делают лабораторную часть.

Глубина секвенирования (sequencing depth) – при получении сырых геномных данных, каждый участок генома оказывается (чаще всего) пройденным несколько раз. Количество раз, которое был отсеквенирован каждый нуклеотид, называется глубиной покрытия и обычно описывается как 10Х, 20Х, 50Х и т.д. Глубина покрытия позволяет выявить возможные ошибки считывания нуклеотидов на машины и определить истинно полиморфные позиции в геноме.

Покрытие (coverage) – в настоящее время понятие покрытие используется как минимум в трех значения

  1. Теоретическая глубина секвенирования – рассчитывается по формуле:

    (общее число полученных ридов * длина рида ) / теоретическая длина генома

  2. Теоретический или эмпирический размах покрытия собранного генома, рассчитывается как:

    размер собранного генома / теоретическая длина генома

  3. Эмпирическая глубина покрытия, рассчитывает как:

    (общее число полученных ридов * длина рида ) / размер собранного генома

Граф де Брейна (de Bruijn graph) – математический алгоритм, лежащий в основе большинства программ-сборщиков ридов в контиги (этап 2 выше). В связи с тем, что интернете (в том числе в википедии) есть много статей детально описывающих и математическую, и техническую компоненту алгоритма, мы не будем останавливаться на принципах работы графа де Брейна в данной статье. Однако отметим, что наиболее важным с прикладной точки зрения параметром алгоритма является т.н. размер кмера (от англи. k-mer size, читается как кей-мер).

Кмер (kmer) –это фрагмент нуклеотидной последовательности фиксированной, часто небольшой, длины (К). Обычно длина должна быть делима на 4. Типичный размер кмера – 24, 36, 48, 56, 96 и т.д. Кмер используется для нахождения перекрывающихся участков между отдельным ридами. Логичным является следствие, что чем длина кмера меньше, тем более вероятно нахождение последовательности кмера во многих ридах.

Контиг (contig) – это консенсусная непрерывная нуклеотидная последовательность составленная из отдельных перекрывающихся ридов. Контиги являются первым исходныи результатом обработки сырых ридов с помощью программы-сборщика генома.

Скаффолд (scaffold) – это последовательность контигов расположенных и ориентированных в порядке их расположения на (теоретической) хромосоме. Для получения скаффолдов необходимы как сами контиги, так и вспомогательная информация (например, референсный геном, или данные paired-end) которая позволит расположить континги относительно друг друга.