Правдоподобие и баейсов подход – как это работает?

Тагир Самигуллин 2 октября 2014

Модели эволюции нуклеотидных последовательностей

JC K2P GTR Процесс замещения одного основания другим…

Модели семейства GTR

JC K2P K3P SYM F81 HKY85 TN93 GTR Частоты оснований равны Частоты оснований не равны Одна скорость замен 2 3 6 скорости замен (транзиции и трансверсии) скорости замен (трансверсии и 2 типа транзиций) скоростей замен + + + + + + + + + + + + + + + +

Подразумевается, что: •Эволюция последовательностей – случайный процесс •Частоты оснований – постоянные •Замены происходят независимо друг от друга •Вероятности замен не меняются со временем (гомогенный эволюционный процесс)

Модели для аминокислотных последовательностей

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

Эмпирические матрицы : Дейхов (Dayhoff)

Dayhoff et al. 1978:

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

значения умножены на 10 000

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

Эмпирические матрицы : другие и BLOSUM

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

Правдоподобие – вероятность данных для выбранной модели

Модель

= дерево и

модель эволюции признаков Модель эволюции признаков

= состав … Одна последовательность, один нуклеотид

А

Правдоподобие = ?

Для модели “

100% А

” = 1 Для модели “

100% С

” = 0 Для модели “

30% А

” = 0.3

Одна последовательность, два нуклеотида

АC

Правдоподобие = ?

Для модели “

4 равновероятных нуклеотида

“ : ¼ x ¼ = 1/16 Для модели “

40% A, 10% C

“ : 0.4 x 0.1 = 0.04

Модель эволюции признаков

= состав и процесс

Данные:

Две последовательности по одному нуклеотиду

А C A

(

Модель: 0.25) C( 0.25) A↔C = 0.4

Правдоподобие ветви между последовательностями: 0.25 x 0.4 = 0.1

A C G T

состав процесс

Модель: A C G A C G T Правдоподобие ветви между последовательностями : T Данные:

Две последовательности по 4 нуклеотида

ССАT CCGT

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

Likelihood

И, наконец, правдоподобие простейшего дерева:

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

Модели эволюции нуклеотидных последовательностей еще раз…

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

Ключевое понятие – апостериорная вероятность

априорная вероятность дерева Pr(T|D) = Pr(T) Pr(D|T) Pr(D) апостериорная вероятность топологии Pr (T,D) = Pr (D,T) Pr(T) Pr(D|T) = Pr(D) Pr(T|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 шагов имеем:

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

После 10 000 шагов :

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

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

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

MCRobot Lewis, 2006

Короткие ветви и байесовская филогения

100 наборов, 1000 оснований, модель К2Р

Поддержка

неверных узлов

Количество

верных узлов Байесовский метод может присвоить коротким ветвям очень высокие значения апостериорной вероятности (коротким – это от 1,3-1,4 ожидаемых замен). Для парсимонии, например, для поддержки 95% требуется минимум 3 ожидаемых замены. Количество верно разрешенных узлов для байесовского метода выше, чем для парсимонии, но и очень короткие 50%,

неверные

ветви получают поддержку выше, чем дает метод максимальной экономии. Эта поддержка в некоторых случаях превышает то есть в 50% консенсусном дереве могут появиться неверные короткие ветви!

Alfaro et al., 2003

Конвергенция МСМС

…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% будут иметь такую же высокую поддержку.

“

Essentially, all models are wrong, but some are useful

.

” George Box <=тенденция или y= a x+ b все-нюансность=> y= a x 6 + b x 5 + c x 4 + d x 3 + e x 2 + f x+ g “

With four parameters I can fit an elephant, and

with five I can make him wiggle his trunk.

John von Neumann ” соответствующим повышением правдоподобия, в противном случае выбирается более простая модель.

?

Критерии для выбора модели

AIC = Akaike Information criterion AIC = -2 lnL +

2k,

где

k = число параметров модели

преимущество при

Δ

AIC > 10 сильное > 4 < 2 слабое никакое AICс= AIC с поправкой на малые наборы данных (

n/k<40

)

n = длина выравнивания

BIC = Bayesian Information criterion BIC = -2 lnL +

k

ln(

n

) Модель с меньшим BIC предпочтительнее Bayes factor , преимущество : 3 - 20 = есть 20 - 150 = сильное 150+ = очень сильное DT (Decision Theory)

Сравнение топологий

Если топология дерева не совпадает с имеющейся гипотезой, значит ли это, что данные отвергают гипотезу?

AU (Approximately Unbiased) test даст ответ CONSEL, PAUP, …?

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

Практический вывод: наилучшая реконструкция seq 1 seq 2 seq 3 2 наименьшим числом замен. seq 6

эвристическим путем

.

А

3 6

С

7

Т С А

6 3

А

максимально экономным (МР tree). Поиск такой топологии идет 2

Т С

7

Словарик Гомоплазия : Синапоморфия Аутапоморфия – мутация, унаследованная потомками от предкового таксона – мутация, характерная для таксона – независимое появление одной и той же мутации у разных таксонов Симплезиоморфия – унаследованное потомками от предка «древнее» состояние признака

seq 1 seq 2 seq 3 seq 4 seq 5 seq 6 seq 7

О признаках в методе МР

Признаки: •постоянные (инвариабельные) •изменчивые (вариабельные).

Последние делятся на информативные и неинформативные. T A T seq2 seq3 seq4 T seq1 G G A seq2 seq3 seq4 A seq1 A G A G seq4 seq2 seq1 seq3 T Т А seq2 seq4 seq3 T seq1

дерево 2.1

дерево 2.1 и дерево 2.2

требуют одинакового числа замен

дерево 2.2

дерево 1.1

дерево 1.2

требует меньшего числа замен, т. е. более экономно, чем дерево 1.2

Информативные признаки позволяют предпочесть одну топологию дерева другой.

Оценка длины дерева

Алгоритм Санкова Позволяет придать разным заменам разный вес: Tv :Ts = 4 : 1 Алгоритм Фитча 3 2 {GA} 1 {A} 1 С A Топология оптимальна?

Для первой клетки: [1;12;2;12] + [8;4;9;6] = 5

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

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