Monthly Archives: August 2013

Детектирование позитивного отбора в геномных последовательностях. Теория.Detecting positive selection in genetic sequences. Theory

Одним из наиболее интересных объектов приложения филогенетических реконструкций в эру геномики является поиск нуклеотидных позиций и ветвей филогенетического дерева находящихся (или находившихся) под позитивным отбором. Сразу оговоримся, что в данном случае термин позитивный отбор (positive selection) не является полной аналогией термина «положительный отбор» принятого в эволюционной биологии. Под позитивным отбор подразумевается следующее.

При сравнении двух нуклеотидных белок-кодирующих последовательностей может быть расчитана частота замен нуклеотидов приводящих к несмысловым (кодирование той же аминокислоты, synonymous substitution) и смысловым  (кодирование отличающейся аминокислоты, non-synonymous substitution) заменам аминокислот. Число теоретически возможных смысловых замен для данной последовательности (суммируется для всех кодонов) обозначается как Мs, а число теоретически возможных несмысловых как Мa (также сумма всех кодонов). Далее, подсчитывается общее число несмысловых замен для двух сравниваемых последовательностей и обозначается как Na. Так же расчитывается число смысловых замен для двух сравниваемых последовательностей и обозначается как Ns. Теперь мы можем расчитать скорость несмысловых замен Ka как отношение числа реализованных несмысловых замен (то, что произошло) к теоретически возможному числу несмысловых замен. Аналогично мы расчитываем скорость смысловых замен Ks как отношение числа реализованных (видимых) смысловых замен к теоретически возможному числу смысловых замен. Отношение этих двух величин Ka/Ks (или ωdN/dS– читается как «диэн-диэс» или омега) и явлется основной для оценки типа отбора действующего на две нуклеотидных последовательности. Если отношение  < 1 , то считается, что последовательности находятся под т.н. стабилизирующим отбором, если значение числа примерно равно 1, то последовательность эволюционирует нейтрально, если же значение > 1, то считается, что последовательность находится под положительным (позитивным) отбором. Получение этой эмпирической закономерности (<1, ~1, >1) следудет из основ популяционной генетики и не является предметом данной статьи.

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

Оговоримся сразу, что приведенный выше расчет является наиболее примитивным. Так, этот метод не учитывает такие важные факторы, как систематическая неравномерность в заменах транзиций и трансверсий (transition/transversion rate bias), частота использования кодонов (codon-usage bias), насыщение заменами с течением времени и т.д.

Начиная с 1994 года были разработаны более совершенные модели, позволяющие учитывать перечисленные выше факторы. Существует четыре класса моделей dN/dS:

  1. Расчет dN/dS для пары выровненных последовательностей (whole alignment model)
  2. Расчет dN/dS для отдельных нуклеотидной позиции (site model)
  3. Расчет dN/dSдля отдельных ветвей филогенетического дерева (branch model)
  4. Расчет dN/dS  для отдельных нуклеотидных позиций с учетом топологии филогенетического дерева (branch-site model).

Попробуем разобраться с механикой современных моделей, используемых для расчета dN/dS.

Для моделирования эволюции кодонов используются модели на основе непрерывных цепей Маркова (НЦМ). НЦМ это математическая модель, которая принимает конечное число состояний (states), и где время, проведенное в каждом состоянии должно быть неотрицательным и иметь экспоненциальное распределение. Непрерывная цепь Маркова низших порядков хороша тем, что условное распределение последующего состояния цепи зависит только от текущего состояния и не зависит от всех предыдущих состояний. Вероятности переходов между состояниями задаются т.н. Q-матрицей переходов.

Переведем на более понятный биологический язык написанное выше. Упрощенный пример цепи Маркова для двух нуклеотидов приведен на рисунке (а). В данном случае, цепь состоит из двух состояний – А и Т. С течением времени между ними могут осуществляться переходы (нуклеотидные замены), определяемые вероятностью с. Кроме того, есть некоторая вероятность, что нуклеотид Т останется в прежнем состоянии, и вероятность этого равна а. Аналогично, есть вероятность b, что нуклеотид А также не изменится. Матрица переходов Q будет задана тогда формой представленной на рисунке (b). Разумеется, это очень упрощенная схема, и в реальности матрица переходов должна была бы включать не два, а четыре нуклеотида.

Схема (a) и (b)

Схема (a) и (b)

Усложним модель.

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

Схема (с) и формула 1)

Схема (с) и формула 1)

Самым важным элементом формулы 1 является матрица Q. Она подобна той матрице, что мы рассмотрели выше для двух нуклеотидов, и элементы этой матрицы состоят из следующих параметров: (1) стационарные частоты нуклеотида π в нуклеотидной последовательности; (2) параметр каппа k определяющий скорость трансверсий и транзиций;  (3) параметр омега ω – отношение несмысловых замен к смысловым. Кроме того, в формулу 1 входит параметр t, являющиеся параметром определяющим длины ветвей дерева (в данном случае t1 и t2). Таким образом, на основе матрицы переходов Q и длин ветвей дерева t мы можем посчитать вероятность перехода из некоторого предкового состояния (обозначено знаком вопроса на схеме) в текущее. Так как мы не знаем, какой был нуклеотид в основании дерева, нам нужно просуммировать вероятности переходов из всех возможных предковых состояний в наблюдаемые текущие  – формула 2 на схема (d).

pr34_3

Таким образом, мы можем перебрать возможные параметры t1, t2, π, ω и k и найти те из них, которые дают максимальную вероятность наблдения нуклеотидов T и А. Обратите внимание, что если речь идет лишь о нуклеотидных последовательностях, а не о кодонах, то параметр ω не имеет никакого смысла и не может быть определен. В связи с этим еще раз немного усложним модель.

Для моделей с кодонами (а именно такая модель позволяет оценить наличие позитивного отбора) мы считаем, что процесс изменения кодонов во времени может быть описан НЦМ, где состояниями цепи Маркова являются 61 смысловых кодона (три несмысловых кодона не учитываются, потому что мы рассматриваетм только белок-кодирующие последовательности). Раз мы моделируем переходы между кодонами, то часть из этих переходов становится смысловым. Кроме того, мы знаем какие нуклеотидные замены являются трансверсиями, а какие транзициями. В связи с этим для каждого элемента матрицы Q мы вводим дополнительное «дробление» на категории учитывающие является переход смысловым или нет, с трансверсией или с транзицией. Обратите внимание, на формулу схемы (с) – для каждого возможного перехода стационарная частота нуклеотида (а в данном случае мы будем говорить уже о кодонах) π разбивается на пять возможных классов. Это связано именно с тем, что часть замен кодонов будут смысловыми, плюс мы учитываем приводит ли к замене трансверсия или транзиция.

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

Feature image source: Positive selection acting on orthologs. Map of the amino acids involved in the binding of ligand [13] on the SAL1 3D structure are in blue, and positively selected sites in pink (PDB: 1GM6).