Байесовский подход в филогенетике позволяет получить наиболее вероятное филогенетическое дерево при заданных исходных данных, последовательностях ДНК или белков рассматриваемых организмов и эволюционной модели замен. Для снижения вычислительной сложности алгоритма расчёт апостериорной вероятности реализуется различными алгоритмами, использующими метод Монте-Карло для марковских цепей. Главными преимуществами байесовского подхода по сравнению с методами максимального правдоподобия и максимальной экономии является вычислительная эффективность, способность работать со сложными моделями эволюции, а также то, что, в отличие от методов, указывающих на единственное наилучшее по заданному критерию дерево, он позволяет выбрать несколько вариантов филогенетического дерева с наибольшим значением апостериорной вероятности.
История и основы метода
Байесовский подход является развитием вероятностного метода, разработанного английским математиком и священником Томасом Байесом на основе теоремы Байеса. Этот метод был опубликован в 1763 году, через два года после его смерти. Позднее современную формулировку теоремы вывел Пьер-Симон Лаплас .
В 1953 году Николас Метрополис ввёл методы Монте-Карло для марковских цепей (MCMC, Markov chain Monte Carlo). Преимущества в скорости вычислений и возможность интеграции с методами MCMC позволили байесовскому подходу стать одним из самых популярных методов статистического вывода. Байесовский подход имеет множество применений в молекулярной филогенетике и систематике. По сравнению с другими методами построения филогенетических деревьев (метод максимальной экономии, метод максимального правдоподобия), он позволяет учитывать филогенетическую неопределенность, использовать априорную информацию и сложные модели эволюции, для которых традиционные методы имеют вычислительные ограничения.
Применение байесовского подхода в филогенетике состоит в следующем. Всё множество допустимых филогенетических деревьев описывается дискретными параметрами (топология деревьев) и непрерывными параметрами (длины ветвей деревьев и параметры эволюционной модели замен). Чтобы вычислить значение апостериорной плотности распределения вероятностей для дерева c топологией и параметрами , при заданных исходных данных , применяется формула Байеса , где — условная плотность распределения вероятностей исходных данных . Знаменатель в этой формуле вычисляется по формуле полной вероятности в виде суммы по интегралов от произведения по , где — априорная плотность распределения для деревьев. Явные аналитическое расчеты по этой формуле не всегда возможны, а численные — требуют большого количества вычислений при поиске максимума функции по . Применение метода статистических испытаний (который также называется методом Монте-Карло) на цепях Маркова позволяет получить приближенные значения апостериорных вероятностей и уменьшить вычислительную сложность алгоритма поиска наиболее вероятного дерева по критерию максимума апостериорной вероятности.
В методах MCMC апостериорная плотность вычисляется за счет имитации работы цепи Маркова, состояниями которой являются филогенетические деревья. Расчет апостериорной плотности выполняется как частота посещения этих состояний в установившемся режиме. Наиболее вероятное дерево определяется по максимальной частоте того состояния, которое чаще всего посещается, или нескольких из них наиболее часто посещаемых. Методы MCMC можно описать двумя этапами: на первом применяется стохастический механизм для получения нового состояния цепи Маркова; на втором выполняется расчет вероятности перехода в это состояние и разыгрывается случайное событие смены состояния. Эта процедура повторяется тысячи или миллионы раз. Доля времени, в течение которого одиночное дерево посещается в процессе работы цепи Маркова, является достаточно точной аппроксимацией для его апостериорной вероятности. К числу наиболее часто применяемых алгоритмов, использующихся в методах MCMC, относятся алгоритм Метрополиса — Гастингса, алгоритм Метрополиса в сочетании с MCMC (MC³) и алгоритм LOCAL Ларгета и Симона.
Алгоритм Метрополиса — Гастингса Алгоритм Метрополиса — Гастингса является одним из наиболее распространенных методов MCMC и представляет собой модифицированную Гастингсом версию алгоритма Метрополиса. Алгоритм Меторополиса — Гастингса строит случайную реализацию цепи Маркова, состояниями которой являются филогенетические деревья. При имитации изменения состояния на каждом шаге выполняется переход от одного дерева к другому за счет изменения топологии или параметров эволюционной модели по определённому правилу. Алгоритм состоит из следующих шагов: выбирается стартовое дерево со случайной топологией и параметрами модели и принимается как текущее; строится дерево с новой топологией или новыми параметрами модели; вычисляется отношение вероятностей (или функций плотности вероятности) нового дерева и старого дерева : (под подразумевается условная вероятность или плотность распределения при заданных исходных данных ); если , то в качестве текущего дерева принимается новое дерево ; если , то выбор дерева происходит с вероятностью (для этого генерируется случайное равномерно распределенное число в интервале , и если это число меньше , то в качестве текущего дерева принимается , иначе остается ;
далее процесс повторяется с шага 2) n раз до тех пор, пока не достигнет равновесного распределения.
В оригинальном алгоритме Метрополиса предполагается, что вероятности переходов от дерева к дереву и обратно равны. Если это условие не выполняется, то применяются поправки Гастингса, состоящие в следующем: вероятность перехода вычисляется по формуле , где — совместная функция распределения.
Алгоритм Метрополиса в сочетании с МСМС Алгоритм Метрополиса в сочетании с MCMC (Metropolis-coupled MCMC, MC³), известный также как алгоритм параллельного отжига, является модифицированной версией алгоритма Метрополиса — Гастингса для марковских цепей со сложным и многомодальным распределением вероятностей состояний. Для этих случаев алгоритмы эвристического поиска деревьев критериями MP (метод максимальной экономии, maximum parsimony), ML (метод максимального правдоподобия) и ME (метод минимальной эволюции), а также МСМС могут выйти на локальный максимум, что приведёт к неверной аппроксимации апостериорной плотности распределения вероятностей. Алгоритм MC³ за счёт смешивания марковских цепей с разной температурой позволяет правильно аппроксимировать распределение апостериорных вероятностей и избегать попадания в локальные оптимумы.
Алгоритм параллельно запускает цепей, по итераций в каждой цепи с разными стационарными распределениями , , где первое распределение с целевой плотностью называется холодной цепью, а другие цепи с распределениями , называются разогретыми. Плотности распределений разогретых цепей имеют вид:
где — температурный фактор.
Возведение плотности в степень при имеет эффект уплощения распределения по аналогии с нагреванием металла. В таком распределении легче перемещаться между пиками, разделёнными долинами, чем в первоначальном распределении. После каждой итерации алгоритм предписывает выполнить обмен состояний между двумя случайно выбранными цепями с помощью предложенного Метрополисом шага. Обмен между состояниями и происходит с вероятностью:
где — текущее состояние в цепи с номером , .
Эвристически, разогретые цепи будут посещать локальные пики довольно легко, и обмен состояниями между цепями позволит холодной цепи иногда перепрыгивать через долины. Если слишком мало, обмен состояниями будет редко выполняться, поэтому в алгоритме используются несколько цепей с разными температурными факторами для улучшения смешивания.
Для получения стационарного распределения вероятностей используются только состояния из холодной цепи, а состояния из разогретых цепей отбрасываются.
Алгоритмы GLOBAL и LOCAL Для генерации нового состояния марковской цепи существуют различные вероятностные способы модификации деревьев, например, бисекция с последующим переприсоединением, обмен ветвей, замена на ближнее соседнее дерево. Алгоритмы LOCAL и GLOBAL предлагают другой способ построения нового дерева по текущему за счёт изменения топологии и длин ветвей. Это приводит к значительному сокращению вычислений для больших деревьев по сравнению с алгоритмами бутстрэпа для методов максимального правдоподобия и максимальной экономии.
Общая идея заключается в том, что дерево представляется в виде следующих параметров: топология дерева и длина его ветвей, а также параметры модели замен. При изменении состояний марковской цепи выполняются последовательные шаги, в которых отдельно либо меняется топология дерева и длина его ветвей, либо меняются только параметры модели замен. Решение о переходе к новому дереву в качестве текущего состояния цепи Маркова принимается так же, как в алгоритме Метрополиса — Гастингса, но значение пороговой вероятности вычисляется с использованием параметров модифицированного дерева.
В алгоритме GLOBAL, представленном Мау, Ньютоном и Ларгетом в 1999 году, все длины ветвей дерева изменяются на небольшую величину в каждом цикле. Алгоритм Ларгета и Симона LOCAL предполагает модификацию дерева в небольшой окрестности случайно выбранной внутренней ветви дерева.
Построение нового дерева в алгоритме LOCAL при модификации топологии и длин ветвей выполняется по следующему правилу: равновероятно выбирается произвольное внутреннее ребро дерева с вершинами и . Вследствие того, что филогенетическое дерево должно быть бинарным, а ребро внутреннее, у каждой из вершин обязательно есть две смежные. Смежные вершины для произвольным образом обозначаются буквами и , а смежные вершины для — буквами и . Далее для вершин и равновероятно выбирается по одной смежной, например, и , и рассматривается путь между вершинами и , состоящий из трёх рёбер. Длины этих рёбер модифицируются пропорционально путём умножения на случайное число по правилу , где — старая длина пути, — новая длина пути, — это равномерно распредёленная случайная величина на отрезке , а — это положительный настраиваемый параметр. Следующий шаг модификации дерева состоит в отсоединении одной из вершин, или , выбираемых равновероятно, и присоединения её в случайно выбранной по равномерному закону точке на пути от вершины до вершины вместе с её дочерней ветвью. При такой модификации возможно изменение топологии дерева, если порядок следования вершин и вдоль пути от к изменился, иначе — топология дерева не изменяется. Поправка Гастингса равна квадрату отношения длин нового и старого пути: .
При модификации параметров модели в алгоритме рассматриваются два варианта: в первом варианте, когда один параметр ограничен набором значений , новое значение параметра вычисляется прибавлением равномерно распределённой случайной величины из интервала . Если новое значение выходит за пределы допустимого диапазона , то остаток отражается внутрь этого отрезка. Поправка Гастингса принимается равной 1. Второй вариант составляет случай, когда модифицируется множество параметров, сумма которых равна константе. В этом случае новое множество значений этих параметров выбирается из распределения Дирихле, центрированного по текущим значениям параметров. Поправка Гастингса вычисляется как отношение плотностей Дирихле с новыми и старыми параметрами.
Критика и обсуждение Значения бутстрэп против апостериорных вероятностей. Было показано, что значения бутстрэпа, вычисленные методами максимальной экономии или правдоподобия, обычно ниже, чем апостериорные вероятности, полученные байесовским методом. Это приводит ряду вопросов, например: Приводят ли апостериорные вероятности к завышению вероятности результата? Являются ли значения бутстрэпа более устойчивыми по сравнению с апостериорными вероятностями?