Байесовский подход в филогенетике

Байесовский подход в филогенетике позволяет получить наиболее вероятное филогенетическое дерево при заданных исходных данных, последовательностях ДНК или белков рассматриваемых организмов и эволюционной модели замен. Для снижения вычислительной сложности алгоритма расчёт апостериорной вероятности реализуется различными алгоритмами, использующими метод Монте-Карло для марковских цепей. Главными преимуществами байесовского подхода по сравнению с методами максимального правдоподобия и максимальной экономии является вычислительная эффективность, способность работать со сложными моделями эволюции, а также то, что, в отличие от методов, указывающих на единственное наилучшее по заданному критерию дерево, он позволяет выбрать несколько вариантов филогенетического дерева с наибольшим значением апостериорной вероятности.

История и основы метода

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

Метафора, иллюстрирующая шаги метода MCMC

Байесовский подход является развитием вероятностного метода, разработанного английским математиком и священником Томасом Байесом на основе теоремы Байеса. Этот метод был опубликован в 1763 году, через два года после его смерти. Позднее современную формулировку теоремы вывел Пьер-Симон Лаплас .

В 1953 году Николас Метрополис ввёл методы Монте-Карло для марковских цепей (MCMC, Markov chain Monte Carlo). Преимущества в скорости вычислений и возможность интеграции с методами MCMC позволили байесовскому подходу стать одним из самых популярных методов статистического вывода. Байесовский подход имеет множество применений в молекулярной филогенетике и систематике. По сравнению с другими методами построения филогенетических деревьев (метод максимальной экономии, метод максимального правдоподобия), он позволяет учитывать филогенетическую неопределенность, использовать априорную информацию и сложные модели эволюции, для которых традиционные методы имеют вычислительные ограничения.

Применение байесовского подхода в филогенетике состоит в следующем. Всё множество допустимых филогенетических деревьев $T$ описывается дискретными параметрами $\tau$ (топология деревьев) и непрерывными параметрами $\alpha$ (длины ветвей деревьев и параметры эволюционной модели замен). Чтобы вычислить значение апостериорной плотности распределения вероятностей $f(T\mid B)$ для дерева $T$ c топологией ${\tau }_0$ и параметрами ${\alpha }_0$, при заданных исходных данных $B$, применяется формула Байеса $f(T({\tau }_0,{\alpha }_0)\mid B)=\frac{f(B\mid T)f(T)}{P(B)}$, где $f(B\mid T)$ — условная плотность распределения вероятностей исходных данных $B$. Знаменатель в этой формуле вычисляется по формуле полной вероятности в виде суммы по $\tau$ интегралов от произведения $f(T(\tau ,\alpha ))f(B\mid T(\tau ,\alpha ))$ по $\alpha$, где $f(T(\tau ,\alpha ))$ — априорная плотность распределения для деревьев. Явные аналитическое расчеты по этой формуле не всегда возможны, а численные — требуют большого количества вычислений при поиске максимума функции $f(T\mid B)$ по $T$. Применение метода статистических испытаний (который также называется методом Монте-Карло) на цепях Маркова позволяет получить приближенные значения апостериорных вероятностей и уменьшить вычислительную сложность алгоритма поиска наиболее вероятного дерева по критерию максимума апостериорной вероятности.

В методах MCMC апостериорная плотность вычисляется за счет имитации работы цепи Маркова, состояниями которой являются филогенетические деревья. Расчет апостериорной плотности выполняется как частота посещения этих состояний в установившемся режиме. Наиболее вероятное дерево определяется по максимальной частоте того состояния, которое чаще всего посещается, или нескольких из них наиболее часто посещаемых. Методы MCMC можно описать двумя этапами: на первом применяется стохастический механизм для получения нового состояния цепи Маркова; на втором выполняется расчет вероятности перехода в это состояние и разыгрывается случайное событие смены состояния. Эта процедура повторяется тысячи или миллионы раз. Доля времени, в течение которого одиночное дерево посещается в процессе работы цепи Маркова, является достаточно точной аппроксимацией для его апостериорной вероятности. К числу наиболее часто применяемых алгоритмов, использующихся в методах MCMC, относятся алгоритм Метрополиса — Гастингса, алгоритм Метрополиса в сочетании с MCMC (MC³) и алгоритм LOCAL Ларгета и Симона.

Алгоритм Метрополиса — Гастингса

Алгоритм Метрополиcа — Гастингса является одним из наиболее распространенных методов MCMC и представляет собой модифицированную Гастингсом версию алгоритма Метрополиса. Алгоритм Меторополиса — Гастингса строит случайную реализацию цепи Маркова, состояниями которой являются филогенетические деревья. При имитации изменения состояния на каждом шаге выполняется переход от одного дерева к другому за счет изменения топологии или параметров эволюционной модели по определённому правилу. Алгоритм состоит из следующих шагов:

1. выбирается стартовое дерево $Ti$ со случайной топологией и параметрами модели и принимается как текущее;
2. строится дерево $Tj$ с новой топологией или новыми параметрами модели;
3. вычисляется отношение вероятностей (или функций плотности вероятности) нового дерева $Tj$ и старого дерева $Ti$: $R=\frac{f(Tj)}{f(Ti)}$

(под $f(T)$ подразумевается условная вероятность или плотность распределения при заданных исходных данных $B$);

• если $R⩾1$, то в качестве текущего дерева принимается новое дерево $Tj$;
• если $R<1$, то выбор дерева $Tj$ происходит с вероятностью $R$ (для этого генерируется случайное равномерно распределенное число в интервале $(0,1)$, и если это число меньше $R$, то в качестве текущего дерева принимается $Tj$, иначе остается $Ti$;

1. далее процесс повторяется с шага 2) n раз до тех пор, пока не достигнет равновесного распределения.

В оригинальном алгоритме Метрополиса предполагается, что вероятности переходов от дерева $Tj$ к дереву $Ti$ и обратно равны. Если это условие не выполняется, то применяются поправки Гастингса, состоящие в следующем: вероятность перехода вычисляется по формуле $\frac{f(Tj)}{f(Ti)}\frac{Q(Ti|Tj)}{Q(Tj|Ti)}$, где $Q$ — совместная функция распределения.

Алгоритм Метрополиса в сочетании с МСМС

Алгоритм Метрополиса в сочетании с MCMC (Metropolis-coupled MCMC, MC³), известный также как алгоритм параллельного отжига, является модифицированной версией алгоритма Метрополиса — Гастингса для марковских цепей со сложным и многомодальным распределением вероятностей состояний. Для этих случаев алгоритмы эвристического поиска деревьев критериями MP (метод максимальной экономии, maximum parsimony), ML (метод максимального правдоподобия) и ME (метод минимальной эволюции), а также МСМС могут выйти на локальный максимум, что приведёт к неверной аппроксимации апостериорной плотности распределения вероятностей. Алгоритм MC³ за счёт смешивания марковских цепей с разной температурой позволяет правильно аппроксимировать распределение апостериорных вероятностей и избегать попадания в локальные оптимумы.

Алгоритм параллельно запускает $m$ цепей, по $n$ итераций в каждой цепи с разными стационарными распределениями $f_j(.)$, $j=1,2,\dots ,m$, где первое распределение с целевой плотностью $f_1=f$ называется холодной цепью, а другие цепи с распределениями $f_j$, $j=2,3,\dots ,m$ называются разогретыми. Плотности распределений разогретых цепей имеют вид:

$f_j(T)=f(T{)}^{1/[1+\lambda (j-1)]},\lambda >0,$ где $\lambda$ — температурный фактор.

Возведение плотности $f(.)$ в степень $1/T$ при $T>1$ имеет эффект уплощения распределения по аналогии с нагреванием металла. В таком распределении легче перемещаться между пиками, разделёнными долинами, чем в первоначальном распределении. После каждой итерации алгоритм предписывает выполнить обмен состояний между двумя случайно выбранными цепями с помощью предложенного Метрополисом шага. Обмен между состояниями $i$ и $j$ происходит с вероятностью:

$\alpha =\frac{f_i(T^{(j)})f_j(T^{(i)})}{f_i(T^{(i)})f_j(T^{(j)})},$ где $f^{(j)}$ — текущее состояние в цепи с номером $j$, $j=1,2,\dots ,m$.

Эвристически, разогретые цепи будут посещать локальные пики довольно легко, и обмен состояниями между цепями позволит холодной цепи иногда перепрыгивать через долины. Если $f_i(T)/f_j(T)$ слишком мало, обмен состояниями будет редко выполняться, поэтому в алгоритме используются несколько цепей с разными температурными факторами для улучшения смешивания.

Для получения стационарного распределения вероятностей используются только состояния из холодной цепи, а состояния из разогретых цепей отбрасываются.

Алгоритмы GLOBAL и LOCAL

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

Общая идея заключается в том, что дерево представляется в виде следующих параметров: топология дерева и длина его ветвей, а также параметры модели замен. При изменении состояний марковской цепи выполняются последовательные шаги, в которых отдельно либо меняется топология дерева и длина его ветвей, либо меняются только параметры модели замен. Решение о переходе к новому дереву в качестве текущего состояния цепи Маркова принимается так же, как в алгоритме Метрополиса — Гастингса, но значение пороговой вероятности вычисляется с использованием параметров модифицированного дерева.

В алгоритме GLOBAL, представленном Мау, Ньютоном и Ларгетом в 1999 году, все длины ветвей дерева изменяются на небольшую величину в каждом цикле. Алгоритм Ларгета и Симона LOCAL предполагает модификацию дерева в небольшой окрестности случайно выбранной внутренней ветви дерева.

Построение нового дерева в алгоритме LOCAL при модификации топологии и длин ветвей выполняется по следующему правилу: равновероятно выбирается произвольное внутреннее ребро дерева с вершинами $u$ и $v$. Вследствие того, что филогенетическое дерево должно быть бинарным, а ребро внутреннее, у каждой из вершин обязательно есть две смежные. Смежные вершины для $u$ произвольным образом обозначаются буквами $a$ и $b$, а смежные вершины для $v$ — буквами $c$ и $d$. Далее для вершин $u$ и $v$ равновероятно выбирается по одной смежной, например, $a$ и $c$, и рассматривается путь между вершинами $a$ и $c$, состоящий из трёх рёбер. Длины этих рёбер модифицируются пропорционально путём умножения на случайное число по правилу $m^{\star }=m\exp (\lambda (U_1-0.5))$, где $m$ — старая длина пути, $m^{\star }$ — новая длина пути, $U_1$ — это равномерно распредёленная случайная величина на отрезке $(0,1)$, а $\lambda$ — это положительный настраиваемый параметр. Следующий шаг модификации дерева состоит в отсоединении одной из вершин, $u$ или $v$, выбираемых равновероятно, и присоединения её в случайно выбранной по равномерному закону точке на пути от вершины $a$ до вершины $c$ вместе с её дочерней ветвью. При такой модификации возможно изменение топологии дерева, если порядок следования вершин $u$ и $v$ вдоль пути от $a$ к $c$ изменился, иначе — топология дерева не изменяется. Поправка Гастингса равна квадрату отношения длин нового и старого пути: $(m^{\star }/m{)}^2$.

При модификации параметров модели в алгоритме рассматриваются два варианта: в первом варианте, когда один параметр ограничен набором значений $[0,M]$, новое значение параметра вычисляется прибавлением равномерно распределённой случайной величины из интервала $(-\delta ,\delta )$. Если новое значение выходит за пределы допустимого диапазона $[0,M]$, то остаток отражается внутрь этого отрезка. Поправка Гастингса принимается равной 1. Второй вариант составляет случай, когда модифицируется множество параметров, сумма которых равна константе. В этом случае новое множество значений этих параметров выбирается из распределения Дирихле, центрированного по текущим значениям параметров. Поправка Гастингса вычисляется как отношение плотностей Дирихле с новыми и старыми параметрами.

Критика и обсуждение

• Значения бутстрэп против апостериорных вероятностей. Было показано, что значения бутстрэпа, вычисленные методами максимальной экономии или правдоподобия, обычно ниже, чем апостериорные вероятности, полученные байесовским методом. Это приводит ряду вопросов, например: Приводят ли апостериорные вероятности к завышению вероятности результата? Являются ли значения бутстрэпа более устойчивыми по сравнению с апостериорными вероятностями?

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

Программа MRBAYES

MrBayes Шаблон:Wayback — это бесплатная программа, осуществляющая байесовский анализ филогении. Первоначально написана Джоном Хюльсенбеком и Фредериком Ронкустом в 2001 году. Когда байесовские методы приобрели популярность, многие молекулярные филогенетики стали выбирать MrBayes. Программа использует стандартный алгоритм MCMC и алгоритм Метрополиса, связанный с MCMC.

MrBayes использует МСМС для приближённого вычисления апостериорных вероятностей деревьев. Пользователь может поменять предположения о модели замен, априорных вероятностей и деталей МС анализа. Также программа позволяет удалять и добавлять таксоны и символы для анализа. В программе можно использовать широкий спектр моделей замен -- от стандартной модели подстановок DNA 4х4, также называемой JC69, в которой считается, что частоты оснований равны и все замены нуклеотидов происходят с равной вероятностью, до наиболее общей модели GTR, в которой различаются и частоты оснований, и вероятности замен. Также программа включает несколько 20х20 моделей замен аминокислот, кодоновые и дублетные модели замены ДНК. Программа предлагает различные методы для ослабления предположения о равных скоростях замен в нуклеотидных позициях. MrBayes также может выводить наследственные состояния, содержащие неопределённость филогенетического дерева и параметров модели.

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

MrBayes 3.2 — это новая версия MrBayes, выпущенная в 2012 году. Новая версия позволяет пользователю запускать несколько анализов параллельно. Также она обеспечивает более быстрое вычисление вероятностей и даёт возможность использования ресурсов графического процессора для выполнения этих вычислений. Версия 3.2 предоставляет больше опций для выходных данных, совместимых с программой FigTree и другими программами для просмотра деревьев.

Список программ, использующих байесовский подход

Название программы	Описание	Метод	Авторы	Ссылка на сайт
Armadillo Workflow Platform	Программа, предназначенная для филогенетического и общего биоинформатического анализа	Вывод филогенетических деревьев с использованием методов ML, MP, байесовскийого подхода и др.	E. Lord, M. Leclercq, A. Boc, A.B. Diallo, V. Makarenkov	https://web.archive.org/web/20161024081942/http://www.bioinfo.uqam.ca/armadillo/.
Bali-Phy	Одновременное получение выравнивания и дерева на основе байесовскийого похода	Байесовский вывод выравниваний и филогенетических деревьев	M.A. Suchard, B. D. Redelings	http://www.bali-phy.org Шаблон:Wayback
BATWING	Вывод деревьев методом Байеса с созданием внутренних узлов	Байесовский анализ, демографическая история, метод расщепления популяций	I. J. Wilson, D. Weale, D.Balding	http://heidi.chnebu.ch/doku.php?id=batwing Шаблон:Wayback
Bayes Phylogenies	Вывод деревьев методом Байеса с использованием методов Монте-Карло для марковских цепей и Метрополиса в сочетании с MCMC	Байесовский анализ, множественные, смешанные модели (с автоматическим разбиением)	M. Pagel, A. Meade	http://www.evolution.rdg.ac.uk/BayesPhy.html Шаблон:Wayback
PhyloBayes / PhyloBayes MPI	MCMC семплер для филогенетических реконструкций.	MCMC, вероятностная модель CAT, учитывающая сайт-специфичные нуклеотиды или аминокислоты	N. Lartillot, N. Rodrigue, D. Stubbs, J. Richer	https://web.archive.org/web/20181218053945/http://www.phylobayes.org/
BEAST	Анализ молекулярных последовательностей с помощью MCMC (Bayesian Evolutionary Analysis Sampling Trees)	Байесовский анализ, релаксированные молекулярные часы, демографическая история	A. J. Drummond, A. Rambaut & M. A. Suchard	http://beast.bio.ed.ac.uk Шаблон:Wayback
BUCKy	Байесовское соответствие филогенетических деревьев для генов	Байесовское соответствие с использованием модифицированного жадного консенсуса для неукоренённых квартетов	C. Ané, B. Larget, D.A. Baum, S.D. Smith, A. Rokas, B. Larget, S.K. Kotha, C.N. Dewey, C. Ané	http://www.stat.wisc.edu/~ane/bucky/ Шаблон:Wayback
Geneious (MrBayes plugin)	Средства для исследования геномов и протеомов	Neighbor-joining, UPGMA, плагины MrBayes, PHYML, RAxML, FastTree, GARLi, PAUP*	A. J. Drummond,M.Suchard,V.Lefort и др.	http://www.geneious.com Шаблон:Wayback
TOPALi	Филогенетический вывод	Выбор филогенетической модели, байесовскийий анализ и оценка филогенетических деревьев методом максимального правдоподобия, определение сайтов, находящихся под позитивным отбором, анализ положения точек рекомбинации	I.Milne, D.Lindner и др.	http://www.topali.org Шаблон:Wayback

Приложения

Байесовский подход широко используется молекулярными филогенетиками для различных приложений:

• Вывод филогений
• Вывод и оценка неопределенности филогений
• Датирование с помощью молекулярного анализа
• Вывод черт предковых форм
• Вывод ареалов предковых форм
• Модель динамики диверсификации и вымирания видов
• Объяснение закономерностей распространения патогенных организмов