Модель замен

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

Существует два основных подхода к созданию модели:эмпирический и параметрический (Lio and Goldman, 1998)[1].

Для описания эволюции нуклеотидных последовательностей наиболее часто используются параметрические модели из семейства GTR (General Time-Reversible models). К этому семейству относятся все стандартные модели нуклеотидных замен (JC69, K80 или K2P, F81, HKY85, TN93, GTR), реализованные в стандартных пакетах для филогенетического анализа, таких как MEGA, PAUP, PHYLIP и PHYML. Эти модели различаются по числу оцениваемых параметров и, соответственно, по реалистичности предпосылок.

В самом общем виде модель GTR (Tavaré 1986[2]) не делает никаких предпосылок относительно частот нуклеотидов и вероятности замены между ними. Она требует оценки четырёх параметров для равновесных частот нуклеотидов ${\displaystyle \Pi =(\pi _{1},\pi _{2},\pi _{3},\pi _{4})}$ и шести параметров скорости замен между ними:

${\displaystyle {\begin{aligned}\alpha =r(T\rightarrow C)=r(C\rightarrow T)\\\beta =r(T\rightarrow A)=r(A\rightarrow T)\\\gamma =r(T\rightarrow G)=r(G\rightarrow T)\\\delta =r(C\rightarrow A)=r(A\rightarrow C)\\\epsilon =r(C\rightarrow G)=r(G\rightarrow C)\\\eta =r(A\rightarrow G)=r(G\rightarrow A)\end{aligned}}}$

Итоговая матрица скорости замен в такой модели будет иметь вид:

${\displaystyle Q={\begin{pmatrix}{-(\alpha \pi _{C}+\beta \pi _{A}+\gamma \pi _{G})}&{\alpha \pi _{C}}&{\beta \pi _{A}}&{\gamma \pi _{G}}\\{\alpha \pi _{T}}&{-(\alpha \pi _{T}+\delta \pi _{A}+\epsilon \pi _{G})}&{\delta \pi _{A}}&{\epsilon \pi _{G}}\\{\beta \pi _{T}}&{\delta \pi _{C}}&{-(\beta \pi _{T}+\delta \pi _{C}+\eta \pi _{G})}&{\eta \pi _{G}}\\{\gamma \pi _{T}}&{\epsilon \pi _{C}}&{\eta \pi _{A}}&{-(\gamma \pi _{T}+\epsilon \pi _{C}+\eta \pi _{A})}\end{pmatrix}}}$

Напротив, наиболее простая модель JC (Jukes and Cantor 1969)[3], делает наибольшее число упрощающих предпосылок и требует оценки всего одного параметра. Модель предполагает равные частоты всех нуклеотидов ${\displaystyle \left(\pi _{A}=\pi _{G}=\pi _{C}=\pi _{T}={1 \over 4}\right)}$ и одинаковую скорость замен между всеми нуклеотидами ${\displaystyle \mu }$ (единственный параметр модели). Матрица скорости замен в такой модели принимает вид:

${\displaystyle P={\begin{pmatrix}{{1 \over 4}+{3 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}\\\\{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}+{3 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}\\\\{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}+{3 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}\\\\{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}-{1 \over 4}e^{-t\mu }}&{{1 \over 4}+{3 \over 4}e^{-t\mu }}\end{pmatrix}}}$

Модель JC плоха тем, что не учитывает некоторых важных свойств эволюции ДНК, таких как разные частоты транзиций (замен одного пуринового основания на другое, либо одного пиримидинового основания на другое.) и трансверсий (замен пуринового основания на пиримидиновое, и наоборот), а также неравенство частот нуклеотидов в реальных нуклеотидных последовательностях. Недостатком модели GTR является большое число параметров, затрудняющее оценивание модели. Поэтому на практике часто применяются промежуточные модели, ослабляющие часть предпосылок JC. Например, модель K80 (Kimura 2-parameter 1980)[4] предполагает равные частоты нуклеотидов, но разные скорости транзиций и трансверсий.

Расстояния между кодирующими последовательностями измеряются более точно, если использовать подходы, основанные на сравнении не нуклеотидных, а аминокислотных последовательностей. Изменение аминокислоты в белке происходит вследствие случайной замены нуклеотидов. Однако вероятность того, что произошедшая замена закрепится, зависит от сходства пары аминокислот, между которыми происходит замена. Сходные по свойствам (гидрофобность, размер, заряд и др) аминокислоты чаще замещаются между собой, так как это не влияет на функцию белка. (Это верно для большинства белков, в соответствии с теорией нейтральности и по причине преобладания отрицательного отбора. Для антигенных детерминант ситуация может быть иной). Кроме того, вероятность замещения зависит от частоты встречаемости данных аминокислот в природе, и от количества нуклеотидных замен (1, 2 или 3), которые отличают пару аминокислот.

В отличие от моделей замен ДНК (нуклеотидной модели замен), наиболее распространенные аминокислотные модели замен разработаны эмпирически. Исходя из экспериментальных наблюдений, для каждой аминокислоты можно оценить вероятность её сохранения или замещения каждой другой аминокислотой в различных группах гомологичных белков. Соответственно, каждому событию (сохранение аминокислоты или её замещение другой аминокислотой) можно присвоить определённые значения (положительные или штрафы), в зависимости от вероятности этих событий. Эти значения можно представить в виде таблицы (матрицы) и использовать при расчете аминокислотных дистанций. При таком подходе менее вероятные (менее частые) замещения будут учитываться как значительно больше увеличивающие дистанцию между последовательностями, чем более вероятные замещения. Получаемые матрицы замен симметричны, то есть вероятность замены L->S равна вероятности S->L.

В моделях замен аминокислот не разделяются параметры частоты встречаемости аминокислот и вероятности их замещения друг на друга. В качестве модели замен используются матрицы, в которых эти параметры совмещены (в силу их эмпирической природы). Матрица замен аминокислот — это таблица (матрица) размером 20x20, каждая ячейка которой описывает вероятность замещения аминокислоты из i-й строки на аминокислоту из j-го столбца за определённый отрезок времени. Направление замещения не играет роли (матрицы замен симметричны). Проще всего предположить, что замещение аминокислоты на другую запрещено. Тогда вероятность того, что аминокислота не изменилась, равна 1. Матрица замен для этой модели выглядит следующим образом:

${\displaystyle {\begin{bmatrix}1&0&\cdots &0&0\\0&1&&0&0\\\vdots &&\ddots &&\vdots \\0&0&&1&0\\0&0&\cdots &0&1\end{bmatrix}}}$

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

Матрицы BLOSUM

Рис. 1. Матрица BLOSUM62

Использованный Дэйхофф методологический подход сравнения близкородственных последовательностей не всегда оптимален при анализе последовательностей, связанных дальним родством. Эта проблема решается в матрицах класса BLOSUM (Block substitution matrices, Henikoff and Henikoff, 1992)[11]. Разработка матриц BLOSUM основывается на анализе только консервативных участков (блоков) не близкородственных последовательностей. Эти участки расцениваются как функционально значимые. Используются матрицы BLOSUM62, BLOSUM50, BLOSUM30 и др. (большее число соответствует большей схожести последовательностей). Матрица BLOSUMn вычислена на основе выравнивания последовательностей, сходных не более, чем на n%. (Например, матрица BLOSUM62 основана на выравнивании белков, сходных не более, чем на 62 %). Таким образом, в отличие от матриц PAM, здесь не используется экстраполяция, и матрицы, предназначенные для анализа сильно дивергировавших последовательностей, основаны на анализе блоков, имеющих тот же уровень дивергенции.

Оценки в матрице BLOSUM представляют собой логарифм отношения вероятности того, что две аминокислоты в выравнивании родственны к вероятности того, что они не родственны, а сопоставлены в выравнивании по случайным причинам:

${\displaystyle S_{ij}=\left({\frac {1}{\lambda }}\right)\log {\left({\frac {p_{ij}}{q_{i}*q_{j}}}\right)}}$

Здесь ${\displaystyle p_{ij}}$ — вероятность того, что аминокислоты ${\displaystyle i}$ и ${\displaystyle j}$ родственны и заместили друг друга в процессе эволюции, а ${\displaystyle q_{i}}$ и ${\displaystyle q_{j}}$ — вероятности встречи соответствующих аминокислот ${\displaystyle i}$ и ${\displaystyle j}$ в случайной белковой последовательности. Коэффициент ${\displaystyle \lambda }$ добавлен для приведения значений матрицы к целочисленным.

Матрицы BLOSUM62 используются (по умолчанию) для оценки качества выравнивания белковых последовательностей в таких программах, как BLAST.

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

Цель подхода состоит в том, чтобы проранжировать деревья или в соответствии с их вероятностью Р(данные|lдерево) или (используя Байесовский подход) по их постериорной вероятности P(дерево|данные). Одной из предпосылок этого подхода является выбор модели замен, в соответствии с которой изменяются последовательности на ветках дерева.

Пусть P(x*|T,t.) - это вероятность набора данных при условии дерева. Данные x* - это набор из n последовательностей x^j , j=1...n. T - это дерево с n листьями и последовательностью j на листе j, а t. - это длины ребер дерева. Далее редположим, мы можем определить P(x|y,t), вероятность того, что предковая последовательность y преобразовалась в последовательность x за период времени t. Тогда вероятность дерева Т с предковыми последовательностями, приписанными к вершинам, может быть получена путём перемножения вероятностей для всех ветвей дерева, например:

P(x¹,...,x⁵|T,t.)= P(x¹|x⁴,t₁)P(x²|x⁴,t₂)P(x³|x⁵,t₃)P(x⁴|x⁵,t₄)P(x⁵), где P(x⁵) - это вероятность того, что последовательность x⁵ находится в корне дерева.

Для задания каждой конкретной P(x|y,t), вероятности того, что последовательность х образовалась из последовательности y за время t, используется модель замен. Модели могут применяться как для нуклеотидных, так и для аминокислотных последовательностей. В простейшем случае предполагается, что изменения в каждом сайте происходят независимо, а делеции и вставки не происходят. Такая модель, дополненная вероятностями делеций и вставок, позволяет более реалистично описать марковский процесс, в результате которого последовательность символов (нуклеотидов или аминокислот) изменяется с течением времени.

Обозначим за P(b|a) вероятность того, что символ a был заменен символом b за период времени длиной t. Тогда для двух выровненных последовательностей без гэпов х и у,  ${\displaystyle P(x|y,t)=\prod _{u}P(x_{u}|y_{u},t)}$, где u - позиция в выравнивании.

Для алфавита размером К вероятности замен P(b|a,t) для всех пар символов а и b образуют матрицу S(t) размером ${\displaystyle {K\times K}}$ :

${\displaystyle S(t)={\begin{pmatrix}P(A_{1}|A_{1},t)&\ldots &P(A_{k}|A_{1},t)\\\vdots &\ddots &\vdots \\P(A_{1}|A_{k},t)&\ldots &P(A_{k}|A_{k},t)\end{pmatrix}}\,.}$

Конкретный вид матрицы S(t) рассчитывается на основе выбранной модели замен. Так, в модели JC эта матрица принимает вид:

${\displaystyle S(t)={\begin{pmatrix}r_{t}&s_{t}&s_{t}&s_{t}\\s_{t}&r_{t}&s_{t}&s_{t}\\s_{t}&s_{t}&r_{t}&s_{t}\\s_{t}&s_{t}&s_{t}&r_{t}\end{pmatrix}}\,.}$

При этом r_t=s_t=1/4 при ${\displaystyle t\rightarrow \infty }$

В модели К80 матрица условных вероятностей замен S(t) примет вид:

${\displaystyle S(t)={\begin{pmatrix}r_{t}&s_{t}&u_{t}&s_{t}\\s_{t}&r_{t}&s_{t}&u_{t}\\u_{t}&s_{t}&r_{t}&s_{t}\\s_{t}&u_{t}&s_{t}&r_{t}\end{pmatrix}}\,.}$

Аналогичным образом могут использоваться эмпирические модели замен для белковых последовательностей.

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

Методы, которые используются для выбора наилучшей модели замен:

1)Иерархический тест отношения правдоподобия Тест отношения правдоподобия применяется только для попарного сравнения “вложенных” моделей (Felsenstein, 1981; Huelsenbeck and Crandall, 1997) Для вычисления LRT нужно сравнить оценки правдоподобия двух сравниваемых моделей L1 и L2: LR = 2*(lnL1-lnL2) Распределение LRT статистики близко к распределению хи-квадрат. Поэтому для определения статистической значимости различий между двумя моделями нужно определить число степеней свободы. Оно соответствует числу дополнительных параметров, вводимых в более сложную модель для того, чтобы она совпала с простой. Например, нужно настроить 4 параметра в HKY85, чтобы она совпала с GTR. Этой информации достаточно, чтобы определить критическое значение тестовой статистики из стандартных таблиц.

2)Информационный критерий Акаике (Hurvich and Tsai 1989, Sugiura 1978).

3)Байесовский информационный критерий (Schwarz 1978).

Программы, с помощью которых можно подобрать наилучшую модель: jModelTest, PartFinder, MEGA, TreeFinder.

Р.Дурбин, Ш.Эдди, А.Крогг, Г.Митчисон (2006). "Анализ биологических последовательностей". М.-Ижевск: НИЦ "Регуляторная и хаотическая динамика", Институт компьютерных исследований - 480с. В.В.Лукашов (2009). "Молекулярная эволюция и филогенетический анализ". М.-Бином. Лаборатория знаний - 228с.