ВВЕДЕНИЕ
Задача распознавания регуляторных сигналов является одной из классических задач вычислительной биологии. После того, как началось проведение массовых экспериментов по анализу экспрессии генов (например, [68, 65,]; обзоры в [4, 7]), с одной стороны, и бьши опубликованы полные геномы многих бактерий и некоторых эукариот, что сделало возможным сравнительный анализ процессов регуляции (в частности, [22, 55]; обзор в [21,30]).
Распознавание сайтов регуляции транскрипции (операторов) является сложной вычислительной проблемой, которая далека от решения. В большинстве случаев, недостаточный размер выборки и низкая степень консервативности последовательностей не позволяет создать надежное распознающее правило. В то же время, точные (например, переборные) методы занимают столь большое время, что практически не осуществимы и не пригодны для решения реальных задач.
Исходные понятия
Молекула дезоксирибонуклеиновой кислоты (ДНК) состоит из двух длинных комплементарных полимерных цепей, закрученных в правильную двойную спираль. Каждая цепь — это полинуклеотид, т.е. регулярный полимер, в котором остатки сахара двух соседних нукпеотидов связаны при помощи фосфатных групп. Каждый нуклеотид состоит из остатка дезоксирибозы, фосфатной группы и пуринового или пиримидинового основания. В состав ДНК входят два пурина — аденин (А) и гуанин (G) и два пиримидина - тимин (Т) и цитозин (С). Вдоль цепи последовательность пуриновых и пиримидиновых остатков информативна. Пуриновые и пиримидиновые основания — это плоские плохо растворимые в воде молекулы, которые стремятся расположиться друг над другом перпендикулярно оси спирали. Две цепи удерживаются вместе при помощи водородных связей между парами оснований. Аденин всегда
образует пару с тимином, а гуанин с цитозином. Строгая специфичность спаривания обусловливает комплементарность последовательностей оснований в двух цепях, закрученных в двойную спираль. Сахарофосфатный скелет у двух цепей направлен в противоположные стороны, поэтому вторая цепи читается в противоположном направлении [86].
Белки - это цепи аминокислот, последовательно соединенных пептидными связями. Последовательность нуклеотидов в участке ДНК, кодирующем белок, соответствует последовательности аминокислот в этом белке. При синтезе белка определенные участки ДНК (называемые кодирующими участками или генами) копируются в виде другого полинуклеотида - рибонуклеиновой кислоты, мРНК, отличающийся от ДНК как по химическому составу, так и по выполняемой функции. Подобно ДНК, мРНК образована в виде линейной последовательности нуклеотидов, но имеет два химических отличия: сахарофосфатный скелет вместо дезоксирбозы содержит сахар рибозу, и вместо основания тимина (Т) в ней содержится близкородственное основание урацила (U). мРНК сохраняет информационное содержание того участка ДНК, копией которого она является, а также сохраняет способность к спариванию комплементарных оснований, поскольку U спаривается с А таким же образом, как и Т с А. Синтез молекул мРНК называется транскрипцией ДНК (рис. 1) [81]. Синтез молекул мРНК катализируется ферментами, которые назьшаются РНК-полимеразами. Молекулы мРНК одноцепочечные, относительно небольшой длины, соответствуют одному или нескольким белкам или их цепям. Последовательность нуклеотидов в молекуле мРНК
ДНК
i
Транскрипция
мРНК
Трансляция Б*пок<
Рисунок 1. Процесс синтеза белка.
о
считывается в соотвествии с генетическим кодом и этот процесс называется трансляцией.
Потребность клетки в некоторых белках значительно изменяется во времени, поэтому имеются механизмы регуляции, обеспечивающие изменение уровня синтеза белков в соответствие с потребностью в них. В частности, специальная группа белков контролирует синтез мРНК, регулируя таким образом концентрацию соответствующих ферментов. Такая регуляция может быть как положительной (тогда сам регулирующий белок называется активатором), так и отрицательной (репрессором). Аминокислотная последовательность белка-регулятора, как и любого другого белка, также кодируется в ДНК; определяющий ее ген называется геном-регулятором. Регуляторы специфичны, то есть каждый из них влияет на синтез какого-либо одного или нескольких определенных белков. В работе исследуется случай прокариотов хотя предлагаемый нами алгоритм применим и для гораздо более широкого класса задач. В простейших (прокариотических организмах) эта специфичность достигается специфичностью связывания молекулы белка-регулятора с определенными некодирующими участками молекул ДНК, расположенными непосредственно перед участком, кодирующим мРНК регулируемого набора ферментов. Специфические участки нуклеотидной последовательности, с которыми связываются регуляторные белки, называются сайтами. Много более длинный участок в ДНК, включающий эти сайты и расположенный перед кодирующим участком, называется лидерной областью или апстримом.
Общепринятое и подтвержденное предположение состоит в том, что все сайты связывания одного белка достаточно сходны между собой. Это предположение позволяет поставить задачу поиска набора сайтов связывания одного белка-регулятора в исходном наборе родственных (относительно него) лидерных областей как задачу нахождения набора наиболее сходных фрагментов в этой выборке лидерных областей. Сам такой набор называется сигналом, а слова, входящие в него, естественно
называются сайтами (или иногда - потенциальными сайтами). Обычно структура и некоторые численные характеристики (например, длина) искомого сигнала заранее фиксируются или подбираются в ходе вычислений. Сигналу приписывается некоторое качество, которое тем выше, чем более попарно похожи друг на друга входящие в него сайты. Возможны разные точные определения качества сигнала. Саму задачу нахождения оптимального (наилучшего) по качеству сигнала в данном исходном наборе (выборке) родственных регуляторных областей называют задачей поиска оптимального сигнала. Она имеет некоторую связь с задачей множественного локального выравнивания, но, конечно, никак не сводится к ней. Существует подход, использующий метод поиска сигнала для построения выравнивания нескольких последовательностей [77].
Оценки качества сигнала и способы его описания
Существует несколько способов оценки качества полученного сигнала. Один из них - использование матрицы позиционных весов, элементы которой вычисляются по формуле [84]:
(1),
где а е {А, С, Т, G), C(i, а) — количество появлений нуклеотида а в позиции /, п — количество последовательностей. Коэффициенты А и Ва подбираются таким образом,
' 1 '
чтобы 2^jW (i,a)'ра = 0, а — Y^ У~]\?2 (i,a)-pa =1, где ра- фоновая
«=1 4 ae{A,C,T,G) t=l
вероятность нуклеотида а, 0.1 — общее число позиций в сайте. Фоновые вероятности ра букв исходного алфавита {А, С, Т, G) определяются как частоты вхождений букв в полный геном рассматриваемого организма или в исходную выборку регуляторных областей геномов; иногда в этом качестве используются априорные частоты, как-то характеризующие исходный генетический материал. Используются также более
10
сложные, например, Марковские, статистические модели для нерегуляторного генетического материала [72, 46, 47]. Матрица (1) использовалась в диссертационной работе при исследовании метаболизма глицерол-3-фосфта (см. главу 4).
А А Т Т G А
A G G Т С С
A G G А Т G
A G G С G Т
1 2 3 4 S O
А 4 1 0 1 0 1
С 0 0 0 1 1 1
G 0 3 3 0 2 1
Т 0 0 1 2 1 1
consensus: A G G T G N Рисунок 2. Матрица выравнивания для 4 слов длины б.
Другой способ состоит в том, что сигнал описывается матрицей выравнивания, каждый элемент иа/ которой показывает число появлений каждой буквы а (из того же алфавита) в /-ой позиции сигнала (см. рис.2). По ней строится вероятностная позиционная матрица сигнала:
п Л- п
(2)
Ъ (na,i ae{A,C.T,G)
Значения поправок
обычно выбираются так, чтобы
выполнялось
сс= -JU* гДе п ' число последовательностей в выравнивании (п
= 4 на рис. 2), а сами эти поправки были пропорциональны вероятностям ра появления букв в том материале, где ищется регуляторный сигнал [43]. Заметим, что
/(а, г) = 1 в любом столбце (позиции) /.
ae{A,T,C,G)
11
С помощью этой матрицы вычисляется информационное содержание данного набора сайтов по формуле:
г _ V""* V"^ f(n i\ 1„ / У**'V cx\
i=l a€{A,T,C,G) Pa
Величина информационного содержания (3) иногда используется как характеристика качества найденного сигнала, а вероятностная позиционная матрица -как решающее правило для поиска новых сайтов в исходных и новых регуляторных областях. Другой часто употребляемой характеристикой качества, применяемой как для представления (1), так и для (2), является качество соответствия сайтов, порождающих статистическую модель, этой модели [1,2,3,25,43,65,72,47].
Также для оценки сигнала в диссертационной работе используются сумма попарного сходства всех сайтов, входящих в сигнал, и среднее этого сходства.
F(x>y) - функция, отражающая степень сходства для двух слов д; и у длины /, в данном случае, количество совпадающих букв в них.
S -сигнал длины /, sly...,sk(kпоследовательностей.
к Качество сайта Sj q(sj) = ^F(s,,Sj).
Среднее качество сайта s, в сигнале: p(s,) =------q{s,). Если в найденном сигнале
к-\
все слова одинаковы, то эта величина равна /. Качество сигнала S - Q(S) =
«-I
1 * Среднее качество сигнала S - P(S) =—У\р(я,)
Лучшим считается сигнал с наибольшим значением P(S), а все сайты сигнала имеют качество p(Sj).
12
Иногда удобно представлять сигнал в виде консенсуса - слова, состоящего из наиболее часто встечающихся букв в каждой позиции сигнала. Иногда два, реже три, нуклеотида встречаются одинаково часто, тогда в такие позиции ставится буква отличная от А, Т, С, G, означающая, что здесь может стоять один из нескольких нуклеотидов. Например, w означает, что в данной позиции может стоять А или Т. [11]
Основные алгоритмы поиска регуляторных сигналов
Известные подходы и алгоритмы выделения потенциальных регуляторных сигналов в исходном наборе (предполагаемых) регуляторных областей условно делятся на две группы: оптимизационные и комбинаторные. Некоторые алгоритмы сочетают в себе черты обеих групп. Оптимизационные алгоритмы основаны на некоторой характеристике качества сигнала (например, на его информационном содержании). Далее производится построение цепочки сигналов, так чтобы их качество (иногда говорят: значение функционала) постепенно возрастало. Таким образом, процедура сводится к поиску экстремума некоторого функционала в пространстве всех допустимых сигналов. Таковы алгоритмы максимизации ожидания МЕМЕ [1], стохастические и жадные алгоритмы Gibbs sampler [23] и ряд других (например, имитация теплового отжига и DMS [36]).
Комбинаторные алгоритмы работают также с пространством сигналов, однако в этом случае цель состоит в построении специального слова (консенсуса), представленного во всех или во многих последовательностях из исходной выборки в том смысле, что искомые сайты отклонялись от него (и в этом смысле друг от друга) наименьшим образом (т.е. здесь также присутствует некоторая функция качества — какая-то мера компактности полученного набора сайтов, например его диаметр). К их
13
числу относятся CONSENSUS [28], PROJECTION [8], WINNOWER, SP-STAR [60], MITRA [13] и другие (например, Conslnd и Matlnd, ITB, WORDUP [59]).
Теперь рассмотрим подробнее несколько наиболее известных описанных в литературе алгоритмов распознавания регуляторных сигналов.
Метод максимизации ожидания (ЕМ). Пусть g(j, к) - вероятность того, что искомый сигнал начинается в у позиции в к последовательности, a f(a, i) - вероятность появления буквы а в колонке / для каждой из букв {А, С, G, Т} и 1 <г <1. В процессе работы алгоритма происходит переопределение g и/до тех пор, пока/не будет мало меняться. Это переопределение происходит с помощью правила байесовской вероятности.
v ' ' Р(В) К J
В данной формуле в качестве А берется событие, состоящее в том, что сигнал в / последовательности начинается в j позиции, а в качестве В берется объединение двух событий:
1) Совокупность значений Да, i) для всех / и ее,
2) Что / последовательность равна данной последовательности 5,.
Р(В\А) подсчитывается как произведение значений функции/по всем буквам слова длины /, начинающегося в последовательности St cj-ой позиции.
Р(А) является начальной вероятностью того, что сигнал начнется в /-ой последовательности с у-ой позиции, и она равна 1/(т-1+1), где т — длина последовательности.
Р(В) вычисляется по формуле полной вероятности с перебором всех возможных начал сигнала в /-ой последовательности.
По формуле (4) переопределяется значения g(j, к), которые позволяют нам определить новые значениями; i).
14
Ниже приведен псевдокод основного цикла алгоритма:
ЕМ (выборка, /){
Выбираем начальную точку/(случайно или определяется пользователем) do{
переопределяем g в соответствии с/
переопределяем/в соответствии с g } until (изменения fАлгоритм ЕМ находит модель сигнала.
Алгоритм МЕМЕ представляет собой расширение метода максимизации ожидания для поиска нескольких сигналов в одной выборке последовательностей. В качестве начальных точек берутся все слова длины / из выборки и запускается одна итерация алгоритма ЕМ. Каждый запуск ЕМ определяет одну вероятностную модель сигнала/и g, после чего выбирается лучшая и уже из начальной точки определенной этой моделью запускается алгоритма ЕМ до его сходимости. Получившийся сигнал запоминается, из выборки удаляются все вхождения этого сигнала, и запускается следующая итерация алгоритма МЕМЕ.
МЕМЕ (выборка, /, количество итераций к){ for i = 1 to к {
for каждого слова длины / из каждой последовательности выборки { запускаем 1 итерацию ЕМ с начальной точки определенной из этого слова
}
выбираем лучшую модель сигнала
запускаем ЕМ из начальной точки определенной этим сигналом
запоминаем получившийся сигнал
удаляем этот сигнал из выборки
Эвристическим аналогом ЕМ является довольно широко используемый алгоритм Gibbs sampling. Метод Gibbs Sampling был разработан Геманом и Геманом [23] для восстановления изображений. Впервые был применён для решения задачи МЛВ
15
Лоренсом и соавторами [43]. Различные модификации этого метода часто применяются для поиска слабых сигналов [43,37,72,47,29]. В самой простой версии мы ищем лучший консервативный неразрывный сигнал длины / в виде вероятностной позиционной матрицы (2). Мы считаем, что сигнал встречается во всех последовательностях.
Алгоритм итеративен. Сначала случайно выбираем одно слово длины / из каждой последовательности. Эти слова формируют начальное множество вхождений сигнала. Обозначим позицию слова в г-ой последовательность как (7.
Итерационный шаг:
- берём одну последовательность /. Обычно, их выбирают по очереди, хотя возможны и другие варианты, например, случайный выбор. Существенно, что у всех последовательностей шансы быть выбранными равны.
Вычисляем вероятностную позиционную матрицу из выбранных слов всех последовательностей, кроме /. Обозначим эту матрицу как Р. Возьмем каждое слово длины / из /-ой последовательности и вычислим вес этого слова, относительно матрицы Р. Вес вычисляется как отношение вероятностей данного слова быть случайно порождённым из позиционной модели вероятностной позиционной матрицы и из фоновой модели.
- Разыграем следующее Or случайно из всех слов в / длины / используя распределение вероятностей, определенное весами (вероятность — это вес, нормированный на 1)
Заменяем О, на Or.
Повторять итерационный шаг, пока ситуация не станет стабильной.
Основная идея алгоритма CONSENSUS заключается в том, что в матрицу выравнивания (описанную выше) по одному добавляют слова длины / из последовательностей, еще не включенных в матрицу. После каждого добавления
16
выбирается набор слов с наибольшим информационным содержанием. Алгоритм работает до тех пор, пока в матрице выравнивания не будет содержаться по одному слову из каждой последовательности или то количество, которое определяет пользователь. В качестве начального слова берутся по очереди все слова длины / из всех последовательностей или только из части входных данных, определенной пользователем. На рисунке 3 приведен пример поиска сигнала длины 4 в трех последовательностях длины 5.
Существует группа алгоритмов, решающих задачу поиска так называемого (/, d)-сигнала, т.е. сигнала длины /, который входит в каждую последовательность не более чем с J заменами. Сложность нахождения такого сигнала заключается в том, что два разных вхождения одного сигнала могут различаться в 2d позициях.
Алгоритм WINNOWER строит многодольный граф, вершинами которого являются слова длины /, долями - слова из одной последовательности, а ребрами связаны те вершины, которые отличаются друг от друга не более чем на 2d замен. Многодольный граф - это граф, в котором ребрами соединяются только вершины из разных долей. Сигналу соответствует клика в таком графе. Клика - это такой подграф, в котором любые две вершины соединены ребром. При поиске клики алгоритм сокращает число ребер в графе путем удаления тех, которые не могут быть частью никакой клики заданного размера.
Алгоритм SP-STAR отличается от WINNOWER тем, что строит взвешенный многодольный граф, в котором каждому ребру приписывается вес того, насколько отличаются слова соединенные этим ребром. Теперь задача заключается в поиске клики минимального веса, которая будет соответствовать искомому сигналу.
17
|