Правдоподобие и баейсов подход – как это работает? Тагир Самигуллин 2 октября 2014

JC K2P GTR Модели эволюции нуклеотидных последовательностей Процесс замещения одного основания другим…

Модели семейства GTR JC K2P K3P SYM F81 HKY85 TN93 GTR Частоты оснований равны Частоты оснований не равны Одна скорость замен скорости замен + + (транзиции и трансверсии) 3 скорости замен + + (трансверсии и 2 типа транзиций) 6 скоростей замен + + Подразумевается, что: Эволюция последовательностей – случайный процесс Частоты оснований – постоянные Замены происходят независимо друг от друга Вероятности замен не меняются со временем (гомогенный эволюционный процесс)

Модели для аминокислотных последовательностей Матрица 20 х 20? => слишком много параметров для оптимизации, обычно недостаточно информации. К примеру, матрица скоростей модели GTR будет выглядеть так: Можно использовать математические модели (z.B. модель Пуассона, что эквивалентно модели JC для белков). Можно использовать модели эволюции кодонов (матрица 61 х 61 !!!) Чаще всего используются эмпирические матрицы

Эмпирические матрицы : Дейхов (Dayhoff) Dayhoff et al. 1978: Использованы последовательности ядерных водорастворимых белков, 72 белковых семейства, 1300 последовательности, 1572 замены. Поскольку сходство высоко (> 85%), сразу построено МР-дерево, реконструированы предковые последовательности, посчитаны замены: Таксон Предковая последовательность затем матрица посчитанных замен преобразована в матрицу вероятностей замен:

значения умножены на Значения в этой матрице справедливы для близкородственных белков (диагональ >> вне-диагональ) PAM1, что соответствует одной замене на 100 сайтов. Для отдаленных белков матрица преобразуется (возводится в степень): PAM10… PAM100 (D1) … PAM250…

На сегодняшний день предложены и другие матрицы, получены они либо с использованием похожего подхода (дистанционное дерево) на б Ольшем количестве данных (JTT-модель) либо с использованием похожего подхода (ML-дерево) на еще б Ольшем количестве данных (WAG-модель, mtREV-модель) либо непосредственно из парных выравниваний (BLOSUM) WWYIR CASILRKIYIYGPV GVSRLRTAYGGRKNRG WFYVR … CASILRHLYHRSPA … GVGSITKIYGGRKRNG WYYVR AAAVARHIYLRKTV GVGRLRKVHGSTKNRG WYFIR AASICRHLYIRSPA GIGSFEKIYGGRRRRG block 1block 2block 3 BLOSUM: 2000 блоков из выравнивания 500 семейств родственных белков разного уровня сходства, от 45 до 90% (серия BLOSUM45…90) Эмпирические матрицы : другие и BLOSUM

Правдоподобие – вероятность данных для выбранной модели Модель = дерево и модель эволюции признаков Модель эволюции признаков = состав … Одна последовательность, один нуклеотид А Правдоподобие = ? Для модели 100% А = 1 Для модели 100% С = 0 Для модели 30% А = 0.3 Одна последовательность, два нуклеотида АCАC Правдоподобие = ? Для модели 4 равновероятных нуклеотида : ¼ x ¼ = 1/16 Для модели 40% A, 10% C : 0.4 x 0.1 = 0.04

Модель эволюции признаков = состав и процесс Данные: Две последовательности по одному нуклеотиду АCАC Модель: A( 0.25) C( 0.25) AC = 0.4 Правдоподобие ветви между последовательностями: 0.25 x 0.4 = 0.1 состав процесс A C G T ACGTACGT Данные: Две последовательности по 4 нуклеотида ССАT CCGT Модель: Правдоподобие ветви между последовательностями :

Изменение длины ветви Likelihood

И, наконец, правдоподобие простейшего дерева: Данные Модель Дерево Первый столбец данных: Второй столбец данных: Третий столбец данных: Значение правдоподобия Четвертый столбец данных:

Модели эволюции нуклеотидных последовательностей еще раз… вероятности замен скорости замен JC K2P GTR … и нуклеотидный состав Процесс замещения одного основания другим…

Ключевое понятие – апостериорная вероятность Pr (T,D) = Pr (D,T) Pr(T) Pr(D|T) = Pr(D) Pr(T|D) Pr(T|D) = Pr(T) Pr(D|T) Pr(D) априорная вероятность дерева правдоподобие вероятность данных (маргинальная) апостериорная вероятность совместная вероятность маргинальная вероятность топологии длины ветвей

Распределение плотности апостериорной вероятности : Марковская цепь Монте-Карло (МСМС)

Начальная точка. Следущая точка выбирается случайно, переход на нее определяется следующим правилом: Если плотность её РР выше чем текущей позиции - шаг делается, если нет – делается с ненулевой вероятностью, которая пропорциональна отношению плотностей РР i+1 /РР i. Принципиально важна возможность перехода на более низкую позицию, иначе не удастся исследовать искомое распределение плотности РР!

Пусть соотношение РР i+1 / РР i равно 0.8 и это число мы сравниваем со случайным числом от 0 до 1. Если это число меньше 0.8, шаг принимается. Интервал значений меньше 0.8 шире, чем больше 0.8 ! Марковская цепь Монте-Карло (МСМС) Пусть соотношение РР i+1 / РР i равно 0.5. Интервалы равны, шансы 50/50. Если РРi +1 /РР i станет меньше 0.5, шаг будет чаще отвергаться, чем приниматься!

Марковская цепь Монте-Карло (МСМС) Первый шаг : принимается с вероятностью 1

Марковская цепь Монте-Карло (МСМС) Второй шаг : мог быть принят с вероятностью 0.144

Марковская цепь Монте-Карло (МСМС) Третий шаг : принимается с вероятностью 0.123

Марковская цепь Монте-Карло (МСМС) После 3 шагов имеем:

Марковская цепь Монте-Карло (МСМС) После шагов :

Марковская цепь Монте-Карло (МСМС) Конечный результат : Чем выше плотность РР в некотором интервале, тем чаще он посещается!

Марковская цепь Монте-Карло (МСМС) Lewis, 2006 MCRobot

Короткие ветви и байесовская филогения Alfaro et al., 2003 Поддержка неверных узлов Количество верных узлов 100 наборов, 1000 оснований, модель К2Р Байесовский метод может присвоить коротким ветвям очень высокие значения апостериорной вероятности (коротким – это от 1,3-1,4 ожидаемых замен). Для парсимонии, например, для поддержки 95% требуется минимум 3 ожидаемых замены. Количество верно разрешенных узлов для байесовского метода выше, чем для парсимонии, но и очень короткие неверные ветви получают поддержку выше, чем дает метод максимальной экономии. Эта поддержка в некоторых случаях превышает 50%, то есть в 50% консенсусном дереве могут появиться неверные короткие ветви!

…We then compare their [thirteen convergence diagnostics] performance in two simple models and conclude that all the methods can fail to detect the sorts of convergence failure they were designed to identify. Mary Kathryn Cowles and Bradley P. Carlin, 1996 How can we know that the chain we are sampling from has converged and mixes well? The disappointing answer is that it is impossible to know for certain. JOHN P. HUELSENBECK et al., 2002 Конвергенция МСМС

Схождение (конвергенция) Марковских цепей очень важно для получения корректного результата. Однако, даже отсутствие видимых проблем с конвергенцией не гарантирует, что цепи сошлись, и это главный недостаток метода. Главное преимущество – в разумные сроки можно получить результат в виде топологии с поддержкой ветвей! Nylander et al., 2008

Высокие значения бутстрапа (>85%) часто интерпретируют как высокую достоверность узла, что не совсем верно даже несмотря на статистическую природу бутстрапа. Строго говоря, бутстрап показывает, достаточно ли данных для поддержки узла, нет ли конфликта в данных. Даже полностью неверное дерево может иметь максимальную поддержку узлов! Значения бутстрап-поддержки некоторой группы зависят в первую очередь от количества признаков, поддерживающих группу, и уровня поддержки альтернативной группировки. Если трактовать бутстрап как показатель уровня достоверности, то BP 97% означает, что из неверных ветвей только 3% будут иметь такую же высокую поддержку. Интерпретация бутстрапа

With four parameters I can fit an elephant, and with five I can make him wiggle his trunk. John von Neumann y=ax+b y=ax 6 +bx 5 +cx 4 +dx 3 +ex 2 +fx+g Essentially, all models are wrong, but some are useful. George Box ? Очевидно, что более сложные модели лучше вписываются в данные и более правдоподобны. Однако усложнение модели должно быть оправдано соответствующим повышением правдоподобия, в противном случае выбирается более простая модель.

Критерии для выбора модели AIC = Akaike Information criterion AIC = -2 lnL +2k, где k = число параметров модели преимущество при Δ AIC> 10 сильное > 4 слабое < 2 никакое AICс= AIC с поправкой на малые наборы данных (n/k

Сравнение топологий Если топология дерева не совпадает с имеющейся гипотезой, значит ли это, что данные отвергают гипотезу? AU (Approximately Unbiased) test даст ответ CONSEL, PAUP, …?

6 3 С А 2 7 Т А С Принцип: В случае конфликта данных (наличие гомоплазий) выбирается гипотеза, которая поддерживается б Ольшим количеством синапоморфий. Следствие: Это приводит к уменьшению количества гомоплазий. Словарик : Синапоморфия – мутация, унаследованная потомками от предкового таксона Аутапоморфия – мутация, характерная для таксона Гомоплазия – независимое появление одной и той же мутации у разных таксонов Симплезиоморфия – унаследованное потомками от предка «древнее» состояние признака seq 1 seq 2 seq 3 seq 4 seq 5 seq 6 seq С А 6 7 Т Т Практический вывод: наилучшая реконструкция филогении (филогенетическое древо) та, которая объясняет наблюдаемые состояния признаков наименьшим числом замен. Дерево, для которого число замен является наименьшим, называют максимально экономным (МР tree). Поиск такой топологии идет эвристическим путем.

О признаках в методе МР Признаки: постоянные (инвариабельные) изменчивые (вариабельные). Последние делятся на информативные и неинформативные. seq 1 seq 2 seq 3 seq 4 seq 5 seq 6 seq 7 информативный неинформативный постоянный AGAG seq1seq3 seq2 seq4 GGAA seq1seq4seq2seq3 дерево 1.1 дерево 1.2 требует меньшего числа замен, т. е. более экономно, чем дерево 1.2 Информативные признаки позволяют предпочесть одну топологию дерева другой. TATT seq1seq4seq2seq3 TТАT seq1seq3seq2seq4 дерево 2.1 дерево 2.2 дерево 2.1 и дерево 2.2 требуют одинакового числа замен

Оценка длины дерева Алгоритм Санкова Позволяет придать разным заменам разный вес: Tv :Ts = 4 : 1 Для первой клетки: [1;12;2;12] + [8;4;9;6] = 5 Топология оптимальна? 2 1 ССGAAA {С}{С} {GA} {A} Алгоритм Фитча

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