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

By | 23 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

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