Содержание к диссертации
Введение
ГЛАВА 1. Алгоритм поиска выделения регуляторных сигналов белок-днкового взаимодействия 21
ГЛАВА 2. Тестирование программы 27
2.1. На искусственных выборках 27
2.2. На природных выборках 30
ГЛАВА 3. Применение программы для поиска потенциальных сигналов связывания транскрипционных факторов в ортологичных рядах генов организмов из групп ентеробактерии и бациллы/клостридии 36
ГЛАВА 4. Применение программьі для исследования регуляции метаболизма глицерол-3-фосфата 53
Выводы 64
Работы автора по теме диссертации 65
Список литературы 66
- На искусственных выборках
- На природных выборках
- Применение программы для поиска потенциальных сигналов связывания транскрипционных факторов в ортологичных рядах генов организмов из групп ентеробактерии и бациллы/клостридии
- Применение программьі для исследования регуляции метаболизма глицерол-3-фосфата
Введение к работе
Задача распознавания регуляторных сигналов является одной из классических задач вычислительной биологии. После того, как началось проведение массовых экспериментов по анализу экспрессии генов (например, [68, 65,]; обзоры в [4, 7]), с одной стороны, и были опубликованы полные геномы многих бактерий и некоторых эукариот, что сделало возможным сравнительный анализ процессов регуляции (в частности, [22, 55]; обзор в [21,30]).
Распознавание сайтов регуляции транскрипции (операторов) является сложной вычислительной проблемой, которая далека от решения. В большинстве случаев, недостаточный размер выборки и низкая степень консервативности последовательностей не позволяет создать надежное распознающее правило. В то же время, точные (например, переборные) методы занимают столь большое время, что практически не осуществимы и не пригодны для решения реальных задач.
Исходные понятия
Молекула дезоксирибонуклеиновой кислоты (ДНК) состоит из двух длинных комплементарных полимерных цепей, закрученных в правильную двойную спираль. Каждая цепь — это полинуклеотид, т.е. регулярный полимер, в котором остатки сахара двух соседних нуклеотидов связаны при помощи фосфатных групп. Каждый нуклеотид состоит из остатка дезоксирибозы, фосфатной группы и пуринового или пиримидинового основания. В состав ДНК входят два пурина - аденин (А) и гуанин (G) и два пиримидина - тимин (Т) и цитозин (С). Вдоль цепи последовательность пуриновьгх и пиримидиновых остатков информативна. Пуриновые и пиримидиновые основания - это плоские плохо растворимые в воде молекулы, которые стремятся расположиться друг над другом перпендикулярно оси спирали. Две цепи удерживаются вместе при помощи водородных связей между парами оснований. Аденин всегда
образует пару с тимином, а гуанин с цитозином. Строгая специфичность спаривания обусловливает комплементарность последовательностей оснований в двух цепях, закрученных в двойную спираль. Сахарофосфатный скелет у двух цепей направлен в противоположные стороны, поэтому вторая цепи читается в противоположном направлении [86].
Белки - это цепи аминокислот, последовательно соединенных пептидными связями. Последовательность нуклеотидов в участке ДНК, кодирующем белок, соответствует последовательности аминокислот в этом белке. При синтезе белка определенные участки ДНК (называемые кодирующими участками или генами) копируются в виде другого полинуклеотида - рибонуклеиновой кислоты, мРНК, отличающийся от ДНК как по химическому составу, так и по выполняемой функции. Подобно ДНК, мРНК образована в виде линейной последовательности нуклеотидов, но имеет два химических отличия: сахарофосфатный скелет вместо дезоксирбозы содержит сахар рибозу, и вместо основания тимина (Т) в ней содержится близкородственное основание урацила (U). мРНК сохраняет информационное содержание того участка ДНК, копией которого она является, а также сохраняет способность к спариванию комплементарных оснований, поскольку U спаривается с А таким же образом, как и Т с А. Синтез молекул мРНК называется транскрипцией ДНК (рис. 1) [81]. Синтез молекул мРНК катализируется ферментами, которые называются РНК-полимеразсши. Молекулы мРНК одноцепочечные, относительно небольшой длины, соответствуют одному или нескольким белкам или их цепям. Последовательность нуклеотидов в молекуле мРНК
ДНК
I Транскрипция
мРНК ^^^м
I Трансляция
Булок —в<ш———
Рисунок 1. Процесс синтеза белка.
считывается в соотвествии с генетическим кодом и этот процесс называется трансляцией.
Потребность клетки в некоторых белках значительно изменяется во времени, поэтому имеются механизмы регуляции, обеспечивающие изменение уровня синтеза белков в соответствие с потребностью в них. В частности, специальная группа белков контролирует синтез мРНК, регулируя таким образом концентрацию соответствующих ферментов. Такая регуляция может быть как положительной (тогда сам регулирующий белок называется активатором), так и отрицательной (репрессором). Аминокислотная последовательность белка-регулятора, как и любого другого белка, также кодируется в ДНК; определяющий ее ген называется геном-регулятором. Регуляторы специфичны, то есть каждый из них влияет на синтез какого-либо одного или нескольких определенных белков. В работе исследуется случай прокариотов хотя предлагаемый нами алгоритм применим и для гораздо более широкого класса задач. В простейших (прокариотических организмах) эта специфичность достигается специфичностью связывания молекулы белка-регулятора с определенными некодирующими участками молекул ДНК, расположенными непосредственно перед участком, кодирующим мРНК регулируемого набора ферментов. Специфические участки нуклеотидной последовательности, с которыми связываются регуляторные белки, называются сайтами. Много более длинный участок в ДНК, включающий эти сайты и расположенный перед кодирующим участком, называется лидерной областью или апстримом.
Общепринятое и подтвержденное предположение состоит в том, что все сайты связывания одного белка достаточно сходны между собой. Это предположение позволяет поставить задачу поиска набора сайтов связывания одного белка-регулятора в исходном наборе родственных (относительно него) лидерных областей как задачу нахождения набора наиболее сходных фрагментов в этой выборке лидерных областей. Сам такой набор называется сигналом, а слова, входящие в него, естественно
называются сайтами (или иногда - потенциальными сайтами). Обычно структура и некоторые численные характеристики (например, длина) искомого сигнала заранее фиксируются или подбираются в ходе вычислений. Сигналу приписывается некоторое качество, которое тем выше, чем более попарно похожи друг на друга входящие в него сайты. Возможны разные точные определения качества сигнала. Саму задачу нахождения оптимального (наилучшего) по качеству сигнала в данном исходном наборе (выборке) родственных регуляторных областей называют задачей поиска оптимального сигнала. Она имеет некоторую связь с задачей множественного локального выравнивания, но, конечно, никак не сводится к ней. Существует подход, использующий метод поиска сигнала для построения выравнивания нескольких последовательностей [77].
Оценки качества сигнала и способы его описания
Существует несколько способов оценки качества полученного сигнала. Один из них - использование матрицы позиционных весов, элементы которой вычисляются по формуле [84]:
W(i,a)=A>
flnihi_«l
0).
где а є {А, С, Т, G}, С (і, а) - количество появлений нуклеотида а в позиции і, п —
количество последовательностей. Коэффициенты А и Ва подбираются таким образом,
' 1 '
юновая
чтобы J2w(.^a)'Pa=0' а Т S Х^2(*>)'Р« =1> где Р«~ Фс
»=1 4 a{A,C,T,G} t=l
вероятность нуклеотида а, 0.1 — общее число позиций в сайте. Фоновые вероятности ра букв исходного алфавита {А, С, Т, G} определяются как частоты вхождений букв в полный геном рассматриваемого организма или в исходную выборку регуляторных областей геномов; иногда в этом качестве используются априорные частоты, как-то характеризующие исходный генетический материал. Используются также более
сложные, например, Марковские, статистические модели для нерегуляторного генетического материала [72, 46, 47]. Матрица (1) использовалась в диссертационной работе при исследовании метаболизма глицерол-3-фосфта (см. главу 4).
А А Т Т G А
A G G Т С С
A G G А Т G
A G G С G Т
consensus: A G G Т G N Рисунок 2. Матрица выравнивания для 4 слов длины б.
Другой способ состоит в том, что сигнал описывается матрицей выравнивания, каждый элемент nai которой показывает число появлений каждой буквы а (из того же алфавита) в /-ой позиции сигнала (см. рис. 2). По ней строится вероятностная позиционная матрица сигнала:
п . +с, а,г) =
(2)
ae{A,C,T,G)
Значения поправок
обычно выбираются так, чтобы
выполнялось ^2 са~ -J где п " числ0 последовательностей в выравнивании (п
ae{A,C,T,G)
= 4 на рис. 2), а сами эти поправки были пропорциональны вероятностям ра появления букв в том материале, где ищется регуляторный сигнал [43]. Заметим, что
J~] f(a, г) = 1 в любом столбце (позиции) /.
ae{A,T,C,G)
С помощью этой матрицы вычисляется информационное содержание данного набора сайтов по формуле:
Lq = ± Е /faO-in^M (3)
і=ї a{A,T,C,G) Ра
Величина информационного содержания (3) иногда используется как характеристика качества найденного сигнала, а вероятностная позиционная матрица -как решающее правило для поиска новых сайтов в исходных и новых регуляторных областях. Другой часто употребляемой характеристикой качества, применяемой как для представления (1), так и для (2), является качество соответствия сайтов, порождающих статистическую модель, этой модели [1,2,3,25,43, 65,72,47].
Также для оценки сигнала в диссертационной работе используются сумма попарного сходства всех сайтов, входящих в сигнал, и среднее этого сходства.
F(.x>y) - функция, отражающая степень сходства для двух слов х к у длины /, в данном случае, количество совпадающих букв в них.
S -сигнал длины /, sl,...,sk(k
Качество сайта Sj q(Sj) = ^F(sltSj).
/-і
Среднее качество сайта Sj в сигнале: p(Sj) = a(Sj). Если в найденном сигнале
/С — 1
все слова одинаковы, то эта величина равна /.
Качество сигнала S - Q(S) = ^q(s,)
«-і
1 * Среднее качество сигнала S - P(S) = -^^()
Лучшим считается сигнал с наибольшим значением P(S), а все сайты сигнала
имеют качество p(s,).
Иногда удобно представлять сигнал в виде консенсуса - слова, состоящего из наиболее часто встечающихся букв в каждой позиции сигнала. Иногда два, реже три, нуклеотида встречаются одинаково часто, тогда в такие позиции ставится буква отличная от А, Т, С, G, означающая, что здесь может стоять один из нескольких нуклеотидов. Например, w означает, что в данной позиции может стоять А или Т. [11]
Основные алгоритмы поиска регуляторных сигналов
Известные подходы и алгоритмы выделения потенциальных регуляторных сигналов в исходном наборе (предполагаемых) регуляторных областей условно делятся на две группы: оптимизационные и комбинаторные. Некоторые алгоритмы сочетают в себе черты обеих групп. Оптимизационные алгоритмы основаны на некоторой характеристике качества сигнала (например, на его информационном содержании). Далее производится построение цепочки сигналов, так чтобы их качество (иногда говорят: значение функционала) постепенно возрастало. Таким образом, процедура сводится к поиску экстремума некоторого функционала в пространстве всех допустимых сигналов. Таковы алгоритмы максимизации ожидания МЕМЕ [1], стохастические и жадные алгоритмы Gibbs sampler [23] и ряд других (например, имитация теплового отжига и DMS [36]).
Комбинаторные алгоритмы работают также с пространством сигналов, однако в этом случае цель состоит в построении специального слова {консенсуса), представленного во всех или во многих последовательностях из исходной выборки в том смысле, что искомые сайты отклонялись от него (и в этом смысле друг от друга) наименьшим образом (т.е. здесь также присутствует некоторая функция качества -какая-то мера компактности полученного набора сайтов, например его диаметр). К их
числу относятся CONSENSUS [28], PROJECTION [8], WINNOWER, SP-STAR [60], MITRA [13] и другие (например, Conslnd и Matlnd, ITB, WORDUP [59]).
Теперь рассмотрим подробнее несколько наиболее известных описанных в литературе алгоритмов распознавания регуляторных сигналов.
Метод максимизации ожидания (ЕМ). Пусть g(j, к) - вероятность того, что искомый сигнал начинается в у позиции в к последовательности, а Да і) - вероятность появления буквы а в колонке / для каждой из букв {А, С, G, Т} и 1 <1. В процессе работы алгоритма происходит переопределение g и/до тех пор, пока/не будет мало меняться. Это переопределение происходит с помощью правила байесовской вероятности.
=Р(ДрУ(Л) <4>
В данной формуле в качестве А берется событие, состоящее в том, что сигнал в / последовательности начинается в j позиции, а в качестве В берется объединение двух событий:
Совокупность значений Да і) для всех / и ос,
Что / последовательность равна данной последовательности St.
Р(В\А) подсчитывается как произведение значений функции/по всем буквам слова длины /, начинающегося в последовательности St с/ой позиции.
Р(А) является начальной вероятностью того, что сигнал начнется в /-ой последовательности с /ой позиции, и она равна 1/(т-1+1), где т — длина последовательности.
Р(В) вычисляется по формуле полной вероятности с перебором всех возможных начал сигнала в /-ой последовательности.
По формуле (4) переопределяется значения g(j, к), которые позволяют нам определить новые значения Да, і).
Ниже приведен псевдокод основного цикла алгоритма:
ЕМ (выборка, /){
Выбираем начальную точку/(случайно или определяется пользователем) do{
переопределяем g в соответствии с/
переопределяем/в соответствии с g } until (изменения/<г) return
}
Алгоритм ЕМ находит модель сигнала.
Алгоритм МЕМЕ представляет собой расширение метода максимизации ожидания для поиска нескольких сигналов в одной выборке последовательностей. В качестве начальных точек берутся все слова длины / из выборки и запускается одна итерация алгоритма ЕМ. Каждый запуск ЕМ определяет одну вероятностную модель сигнала/и g, после чего выбирается лучшая и уже из начальной точки определенной этой моделью запускается алгоритма ЕМ до его сходимости. Получившийся сигнал запоминается, из выборки удаляются все вхождения этого сигнала, и запускается следующая итерация алгоритма МЕМЕ.
МЕМЕ (выборка, /, количество итераций к){ for і = 1 to к {
for каждого слова длины / из каждой последовательности выборки { запускаем 1 итерацию ЕМ с начальной точки определенной из этого слова
}
выбираем лучшую модель сигнала
запускаем ЕМ из начальной точки определенной этим сигналом запоминаем получившийся сигнал удаляем этот сигнал из выборки } }
Эвристическим аналогом ЕМ является довольно широко используемый алгоритм Gibbs sampling. Метод Gibbs Sampling был разработан Геманом и Геманом [23] для восстановления изображений. Впервые был применён для решения задачи МЛВ
Лоренсом и соавторами [43]. Различные модификации этого метода часто применяются для поиска слабых сигналов [43,37,72,47,29]. В самой простой версии мы ищем лучший консервативный неразрывный сигнал длины / в виде вероятностной позиционной матрицы (2). Мы считаем, что сигнал встречается во всех последовательностях.
Алгоритм итеративен. Сначала случайно выбираем одно слово длины / из каждой последовательности. Эти слова формируют начальное множество вхождений сигнала. Обозначим позицию слова в г-ой последовательность как О..
Итерационный шаг:
берём одну последовательность /. Обычно, их выбирают по очереди, хотя возможны и другие варианты, например, случайный выбор. Существенно, что у всех последовательностей шансы быть выбранными равны.
Вычисляем вероятностную позиционную матрицу из выбранных слов всех последовательностей, кроме /. Обозначим эту матрицу как Р. Возьмем каждое слово длины / из /-ой последовательности и вычислим вес этого слова, относительно матрицы Р. Вес вычисляется как отношение вероятностей данного слова быть случайно порождённым из позиционной модели вероятностной позиционной матрицы и из фоновой модели.
Разыграем следующее Or случайно из всех слов в / длины / используя распределение вероятностей, определенное весами (вероятность — это вес, нормированный на 1)
Заменяем О і на Or.
Повторять итерационный шаг, пока ситуация не станет стабильной.
Основная идея алгоритма CONSENSUS заключается в том, что в матрицу выравнивания (описанную выше) по одному добавляют слова длины / из последовательностей, еще не включенных в матрицу. После каждого добавления
выбирается набор слов с наибольшим информационным содержанием. Алгоритм работает до тех пор, пока в матрице выравнивания не будет содержаться по одному слову из каждой последовательности или то количество, которое определяет пользователь. В качестве начального слова берутся по очереди все слова длины / из всех последовательностей или только из части входных данных, определенной пользователем. На рисунке 3 приведен пример поиска сигнала длины 4 в трех последовательностях длины 5.
Существует группа алгоритмов, решающих задачу поиска так называемого (/, d)-сигнала, т.е. сигнала длины /, который входит в каждую последовательность не более чем с d заменами. Сложность нахождения такого сигнала заключается в том, что два разных вхождения одного сигнала могут различаться в 2d позициях.
Алгоритм WINNOWER строит многодольный граф, вершинами которого являются слова длины /, долями - слова из одной последовательности, а ребрами связаны те вершины, которые отличаются друг от друга не более чем на 2d замен. Многодольный граф - это граф, в котором ребрами соединяются только вершины из разных долей. Сигналу соответствует клика в таком графе. Клика — это такой подграф, в котором любые две вершины соединены ребром. При поиске клики алгоритм сокращает число ребер в графе путем удаления тех, которые не могут быть частью никакой клики заданного размера.
Алгоритм SP-STAR отличается от WINNOWER тем, что строит взвешенный многодольный граф, в котором каждому ребру приписывается вес того, насколько отличаются слова соединенные этим ребром. Теперь задача заключается в поиске клики минимального веса, которая будет соответствовать искомому сигналу.
sequence 1
ACTGA
sequence 2 T AGCG
sequence 3 CTTGC
CYCLE 1
'seq — 5.5
CYCLE 2
А С T G T A G С
*seq
/,eq = 4.2
/seq = 4.2
«seq — 2.8
CYCLE 3
А С T G AGCG С T T G
А С T G AGCG T T G С
/.eq =3.2
/.eq = 2.1
Рисунок 3. Пример работы алгоритма CONSENSUS.
Априорная вероятность каждого основания 25%.
Цикл 1. Для простоты, на этом цикле рассматривается матрица выравнивания для первого слова первой последовательности. На самом деле, создаются матрицы выравнивания для всех слов из всех последовательностей, если пользователь не определил подвыборку последовательностей в качестве начальной (это существенно, если выборка очень большая). Для этого примера на этом цикле создается б матриц выравнивания.
Цикл 2. Добавляя по слову из оставшихся двух последовательностей, получается 4 матрицы. Из них выбирается матрица с самым большим информационным содержанием (она обведена жирной линией).
Цикл 3. Только две матрицы создаются на этом цикле путем добавления двух слов из третьей последовательности, еще не вошедшей в матрицу. Жирной линией обведена матрица с наибольшим информационным содержанием, но она не обязательно является лучшей, потому что мы рассмотрели работу алгоритма только для первого слова из первой последовательности. Тоже самое надо проделать для всех слов, и после этого выбрать лучшую матрицу. (Конец подписи к рисунку 3)
*
Алгоритм MITRA (Mismatch Tree Algorithm) использует префиксное дерево для разбиения всех возможных сигналов на непересекающиеся подмножества, имеющие общий префикс с точностью до d замен. Это разбиение делит проблему поиска сигнала на более мелкие задачи, в которой можно применить уже известные алгоритмы, например WINNOWER (см. пример на рис. 4).
В диссертационной работе предложен и тестирован новый алгоритм для выделения сигнала в выборке нуклеотидных последовательностей. Он является промежуточным с точки зрения приведенной выше классификации - производится оптимизация функции качества, однако она определяется через попарное сходство слов, а не как информационное содержание набора.
Дано 3 последовательности длины 6.
ААТАТС
ATCCGT
GCATAC
Ищем (4, 1)-сигнал.
Имеем 9 слов длины 4 из трех последовательностей длины 6:
ААТА 4) АТСС 7) GCAT
АТАТ 5) TCCG 8) САТА
ТАТС 6) CCGT 9) АТАС
На каждом шаге все слова разделяются на части. В узлах дерева находятся слова, совпадающие с префиксом сигнала, с точностью до 1 замены. Мы ищем по одному слову из каждой последовательности, поэтому в процессе работы отсекаются ветки, не содержащие необходимого количества слов.
В данном случае сигналом является последовательность АТАС. Для простоты рассмотрены не все возможные варианты, а только одна ветвь дерева, соответствующая сигналам, начинающимся с А, поэтому этот сигнал может быть не единственным.
Рисунок 4. Пример работы алгоритма MITRA.
На искусственных выборках
Далее итеративно повторяем такое разбиение «вглубь» графа G. А именно, каждую из двух полученных частей графа G снова разбиваем на две (в том же смысле) равные части. Относительно уже этих разбиений один конец ребра определен однозначно: это А для одной части и В для другой, а второй конец ребра выбирается произвольно. На рисунке 6 ребра второго уровня деления - это (1,3) и (2,4). И так далее: каждую появившуюся в этой процедуре не одновершинную часть Р разбиваем на две равные части Pi и Рг так, чтобы ребра новых частей выходили из концов ребра предыдущей части Р. Конечно, при этом каждое Р равно объединению его частей Pj и Р.?. Процедура разбиений прекращается, когда все части Р станут одновершинными; на самом деле, можно остановиться, когда эти части станут просто мелкими (из 1-3 вершин).
ЭтапЗ. Внешний цикл алгоритма состоит во взаимно однозначном приписывании каждой вершине графа G одной из исходных последовательностей. Такое приписывание назовем расстановкой последовательностей по вершинам графа и обозначим г. Каждая следующая расстановка выбирается таким образом, чтобы как можно больше пар последовательностей, не соединенных на предыдущих итерациях, сейчас соединились. Заметим, что из вершин ребра, полученного при первом делении, выходит больше всего ребер, поэтому на следующей итерации в эти вершины ставятся последовательности, из которых, на текущий момент, меньше всего выходило ребер.
Условием выхода из цикла является момент, когда любая пара последовательностей хотя бы раз будет приписана вершинам графа G так, что окажется соединенной в нем ребром. Это обеспечивает разумное количество итераций внешнего цикла при достаточном разнообразии обработанных расстановок (порядка п) и позволяет избежать полного перебора. При каждой расстановке находится своя система слов, которые обрабатываются на этапе 6.
Этап 4. Для текущей расстановки г происходит поиск системы наиболее похожих слов, который мы называем сборкой. Граф G помогает организовать этот поиск. Отметим, что г (А) - это последовательность, приписанная вершине А, в которой, как и во всех других приписанных последовательностях, ищется сигнал, и далее мы будем говорить просто А, понимая, что это некоторая последовательность приписанная вершине А, в соответствии с текущей расстановкой г. Процесс сборки является обратным к процессу деления графа. Теперь мы объединяем мелкие части и при этом находим систему наиболее похожих слов из последовательностей, находящихся в этих частях, выполняя индуктивный шаг, который состоит в следующем (рис. 7). Пусть для двух частей Р; и Р; с ребрами соответственно (А, С) и (В, D), полученных разбиением подграфа Р (в исходном графе G) с ребром (А, В), уже определены два набора каждый из лучших систем. В первом наборе находятся слова из последовательностей, находящихся в части Pi, а во втором наборе - из Р.?, назовем их набор над Pi и Р?, соответственно. Первый набор состоит из продолжений любых двух слов а и с из соответственно последовательностей А и С (где а и с такие, что значение F(a, с) больше некоторого фиксированного порога и) на множество Pi, т.е. любая пара (а, с) определяет лучшую систему слов из последовательностей части Р/. И аналогично, второй набор состоит из продолжений на все множество Рг слов Ъ и d, произвольно взятых в последовательностях В и D. Объединяя таким образом части, мы получим набор над всем графом G, т.е. искомую систему слов (сигнал).
Например, на рисунке 6 пусть Р/ = (1, 5, б, 3} и Рг = {2, 7, 4}; (А, С) = (1, 3) и (В, D) = (2, 4); Р равно объединению множеств Pi и Р (то, что в данном случае Р совпало с G, конечно, случайно); так всегда будет в конце этой индукции, но не на предшествующих ему шагах).
Этап 6. Немаловажное значение имеет обработка результатов. Имеется несколько подходов, один из них получение достаточно представительной статистики. А именно, каждому слову ставится в соответствие число, которое отражает меру того, что это слово входит в искомый сайт. Это число равно сумме качеств по всем полученным сигналам, которые включают данное слово. Здесь под качеством понимается качество не всего сигнала, а качество слова из сигнала (который его содержит) по отношению ко всему этому сигналу, т.е. сумма значений F(x,y), где х - упомянутое слово, а у пробегает все остальные слова этого сигнала. Таким образом, слова, входящие в истинный сигнал, будут помечены в исходных последовательностях числами, которые заметно больше чисел, которые помечают другие слова.
Еще один вариант состоит в том, чтобы рассматривать каждую систему как отдельный потенциальный сигнал или выбрать одну систему, в которой попарная схожесть слов наибольшая. В программе IRSA реализованы оба варианта и для каждой задачи выбирается более подходящий способ представления результатов. Наша практика показывает, что удобнее работать с одной наилучшей системой.
Алгоритм был реализован в виде компьютерной программы, которая использовалась в основном для поиска сигнала в наборе генетических последовательностей с характерной длиной от 100 до 300 нуклеотидов; в этом частном случае сигнал представлял собой систему отдельных слов с характерной длиной от 15 до 30. В этом варианте программа была написана на языке Object Pascal с помощью среды программирования Delphi.
Схема компьютерной реализации алгоритма представлена на рисунке 5. На этапе 1 входные последовательности разбиваются на подслова длины /, вычисляются и запоминаются все значения функции F(x,y) - попарная похожесть всех слов. Затем (2) строится вспомогательный граф G, который определяет порядок, в котором будут перебираться последовательности из выборки. На этапе 3 создается расстановка, т.е. вершинам графа приписываются входные последовательности до тех пор, пока каждая пара не будет хотя бы один раз соединена ребром. Для каждой расстановки происходит процедура сборки. И на последнем этапе происходит обработка результатов. Реализованы две возможности представления результатов, в виде статистики и в виде списка всех систем с пометкой самого лучшего (см. введение).
На природных выборках
Сначала рассматривались регуляторные области из четырех геномов Е. coli, Е. carotovora, Y. enterocolitica, К. pneumoniae, и по ним был получен базисный сигнал, включающий уже известные сайты Е. coli (табл. 12а) с консенсусом TGTTCGATAACGAACA (рис. 12а). По базисному сигналу как по обучающей выборке была построена матрица позиционных весов (табл. 13а) для поиска палиндромных сайтов длины 16. С помощью этой матрицы с порогом 4,1 были еще найдены сайты в дополнительных геномах Y. pestis, S. typhimurium, S. typhi (табл. 12a).
В регуляторных областях из трех геномов V. cholerae, V. vulnificus, V. fischeri был выделен базисный палиндромный сигнал длины 18 с консенсусом AATGCTCGATCGAGCATT (рис. 126). Базисный сигнал включает сайты в геноме V. cholerae перед ортологами генов glpA, glpD, glpT, в V.vulnificus - перед glpA, glpD, glpT, glpF, в V. fischeri - перед glpA, glpD, glpF. Найденные сайты и матрица позиционных весов представлены в табл. 126 и табл. 136 соответственно. При сканировании геномов с использованием этой матрицы, новых потенциальных сайтов обнаружено не было.
В регуляторных областях генов glpD, glpF из четырех геномов семейства P. aeruginosa, P. fluorescens, P. putida, P. syringae был найден палиндромный базисный сигнал wTTTTCGTATACGAAAAw длины 18 (рис. 12в), включающий сайты, ранее предсказанные в работе [67] у P. aeruginosa. По этому базисному сигналу была построена позиционная матрица (табл. 13в), с помощью которой были найдены новые потенциальные сайты связывания GlpR перед генами glpT в P. aeruginosa, P. syringae и P. fluorescens, а также еще один сайт в регуляторной области гена glpD в P. aeruginosa (табл. 12в).
В регуляторных областях гена glpD в геномах а-протеобактерий М. loti, S. meliloti, A. tumefaciens, В. melitensis, R. palustris и гена glpK в S. meliloti и еще одного ортолога гена glpD в A. tumefaciens (на рис. 11 и в табл. Юг он обозначен как glpDl) не удалось обнаружить палиндромного сигнала, но была найдена достаточно большая консервативная область (порядка 30 нуклеотидов), в которой, при подробном изучении, были замечены тандемные повторы. Тогда данный результат был доработан с помощью программы SignalX и, в случае гена glpD в A. tumefaciens и гена glpK в S. meliloti, удалось продлить сигнал. Таким образом, были найдены 3-4 тандемных повтора слова TTTCGTT (рис. 12г), идущих друг за другом через 3-5 нуклеотидов (табл. 12г), которые составили базисный сигнал. При исследовании р-протеобактерий с помощью матрицы позиционных весов (табл. 13г), построенной по этому базисному сигналу, аналогичные повторы были обнаружены перед генами glpD в геномах бактерий рода Burkholderia: В. fungorum, В. pseudomallei, В. cepacia.
Снижение порога при исследовании геномов из группы Enterobacteriaceae, приводит к сильному перепредсказанию, но при этом обнаруживаются относительно слабые сайты перед генами, входящими в ГЗФ регулон (табл. 12г). В то же время, даже при низком пороге не во всех областях футпринта, показанных в статьях [78,42, 80] удается обнаружить потенциальные сайты. Поскольку данные о трехмерной структуре регуляторов семейства DeoR отсутствуют, нет оснований полагать, что эти регуляторы во всех случаях образуют димеры, связывающиеся с палиндромными сайтами. При выравнивании областей перед геном glpD в организмах из группы Enterobacteriaceae: Е. coli, S. typhimurium, К. pneumoniae, Y. enterocolitica можно усмотреть как палиндромную симметрию, так и повторы (табл. 14). Более того, принимая во внимание, что в некоторых других группах организмов нам удалось выделить сигнал в виде тандемного повтора, можно было предположить, что структура сайта в группе Enterobacteriaceae тоже представляет собой тандемный повтор. В тоже время построить хорошую матрицу позиционных весов для этих регуляторных областей не удается и в предположении палиндромной симметрии, и для симметрии типа тандемного повтора.
Итак, нами были идентифицированы сигналы связывания GlpR в Vibrionaceae и Pseudomonadaceae (палиндром) и а-, р-протеобактериях (тандемный повтор). Ситуация в Enterobacteriaceae остается открытой. Выравнивание регуляторных областей и их анализ с помощью программ выделения сигналов позволяет предположить существование как палиндромного сигнала, так и сигнала в виде фазированного прямого повтора. Нельзя исключать того, что мономеры GlpR в этих бактериях могут связываться с регуляторными областями в произвольных направлениях, образуя кооперативные комплексы (ср. регуляторы NarL [24] и FUR [56]). Ситуация может проясниться после появления экспериментальных данных, секвенирования новых геномов из этой группы, или анализа других регуляторов из семейства DeoR.
На рисунке 11 представлены примеры оперонов в некоторых из рассмотренных геномов, регулируемых репрессором GlpR и имеющие соответствующий сайт. Закрашенными кружочками отмечены известные сайты, а незакрашенными — предсказанные, с их весами.
Применение программы для поиска потенциальных сигналов связывания транскрипционных факторов в ортологичных рядах генов организмов из групп ентеробактерии и бациллы/клостридии
Задача распознавания регуляторных сигналов является одной из классических задач вычислительной биологии. После того, как началось проведение массовых экспериментов по анализу экспрессии генов (например, [68, 65,]; обзоры в [4, 7]), с одной стороны, и были опубликованы полные геномы многих бактерий и некоторых эукариот, что сделало возможным сравнительный анализ процессов регуляции (в частности, [22, 55]; обзор в [21,30]).
Распознавание сайтов регуляции транскрипции (операторов) является сложной вычислительной проблемой, которая далека от решения. В большинстве случаев, недостаточный размер выборки и низкая степень консервативности последовательностей не позволяет создать надежное распознающее правило. В то же время, точные (например, переборные) методы занимают столь большое время, что практически не осуществимы и не пригодны для решения реальных задач.
Молекула дезоксирибонуклеиновой кислоты (ДНК) состоит из двух длинных комплементарных полимерных цепей, закрученных в правильную двойную спираль. Каждая цепь — это полинуклеотид, т.е. регулярный полимер, в котором остатки сахара двух соседних нуклеотидов связаны при помощи фосфатных групп. Каждый нуклеотид состоит из остатка дезоксирибозы, фосфатной группы и пуринового или пиримидинового основания. В состав ДНК входят два пурина - аденин (А) и гуанин (G) и два пиримидина - тимин (Т) и цитозин (С). Вдоль цепи последовательность пуриновьгх и пиримидиновых остатков информативна. Пуриновые и пиримидиновые основания - это плоские плохо растворимые в воде молекулы, которые стремятся расположиться друг над другом перпендикулярно оси спирали. Две цепи удерживаются вместе при помощи водородных связей между парами оснований. Аденин всегда образует пару с тимином, а гуанин с цитозином. Строгая специфичность спаривания обусловливает комплементарность последовательностей оснований в двух цепях, закрученных в двойную спираль. Сахарофосфатный скелет у двух цепей направлен в противоположные стороны, поэтому вторая цепи читается в противоположном направлении [86].
Белки - это цепи аминокислот, последовательно соединенных пептидными связями. Последовательность нуклеотидов в участке ДНК, кодирующем белок, соответствует последовательности аминокислот в этом белке. При синтезе белка определенные участки ДНК (называемые кодирующими участками или генами) копируются в виде другого полинуклеотида - рибонуклеиновой кислоты, мРНК, отличающийся от ДНК как по химическому составу, так и по выполняемой функции. Подобно ДНК, мРНК образована в виде линейной последовательности нуклеотидов, но имеет два химических отличия: сахарофосфатный скелет вместо дезоксирбозы содержит сахар рибозу, и вместо основания тимина (Т) в ней содержится близкородственное основание урацила (U). мРНК сохраняет информационное содержание того участка ДНК, копией которого она является, а также сохраняет способность к спариванию комплементарных оснований, поскольку U спаривается с А таким же образом, как и Т с А. Синтез молекул мРНК называется транскрипцией ДНК (рис. 1) [81]. Синтез молекул мРНК катализируется ферментами, которые называются РНК-полимеразсши. Молекулы мРНК одноцепочечные, относительно небольшой длины, соответствуют одному или нескольким белкам или их цепям. Последовательность нуклеотидов в молекуле мРНК считывается в соотвествии с генетическим кодом и этот процесс называется трансляцией.
Потребность клетки в некоторых белках значительно изменяется во времени, поэтому имеются механизмы регуляции, обеспечивающие изменение уровня синтеза белков в соответствие с потребностью в них. В частности, специальная группа белков контролирует синтез мРНК, регулируя таким образом концентрацию соответствующих ферментов. Такая регуляция может быть как положительной (тогда сам регулирующий белок называется активатором), так и отрицательной (репрессором). Аминокислотная последовательность белка-регулятора, как и любого другого белка, также кодируется в ДНК; определяющий ее ген называется геном-регулятором. Регуляторы специфичны, то есть каждый из них влияет на синтез какого-либо одного или нескольких определенных белков. В работе исследуется случай прокариотов хотя предлагаемый нами алгоритм применим и для гораздо более широкого класса задач. В простейших (прокариотических организмах) эта специфичность достигается специфичностью связывания молекулы белка-регулятора с определенными некодирующими участками молекул ДНК, расположенными непосредственно перед участком, кодирующим мРНК регулируемого набора ферментов. Специфические участки нуклеотидной последовательности, с которыми связываются регуляторные белки, называются сайтами. Много более длинный участок в ДНК, включающий эти сайты и расположенный перед кодирующим участком, называется лидерной областью или апстримом.
Общепринятое и подтвержденное предположение состоит в том, что все сайты связывания одного белка достаточно сходны между собой. Это предположение позволяет поставить задачу поиска набора сайтов связывания одного белка-регулятора в исходном наборе родственных (относительно него) лидерных областей как задачу нахождения набора наиболее сходных фрагментов в этой выборке лидерных областей. Сам такой набор называется сигналом, а слова, входящие в него, естественно называются сайтами (или иногда - потенциальными сайтами). Обычно структура и некоторые численные характеристики (например, длина) искомого сигнала заранее фиксируются или подбираются в ходе вычислений. Сигналу приписывается некоторое качество, которое тем выше, чем более попарно похожи друг на друга входящие в него сайты. Возможны разные точные определения качества сигнала. Саму задачу нахождения оптимального (наилучшего) по качеству сигнала в данном исходном наборе (выборке) родственных регуляторных областей называют задачей поиска оптимального сигнала. Она имеет некоторую связь с задачей множественного локального выравнивания, но, конечно, никак не сводится к ней. Существует подход, использующий метод поиска сигнала для построения выравнивания нескольких последовательностей [77].
Применение программьі для исследования регуляции метаболизма глицерол-3-фосфата
Для определения характеристик алгоритма в контролируемых условиях было проведено его тестирование на искусственных выборках. Выборки были созданы двумя разными способами для различных целей.
При первом способе исходная выборка сначала содержала искомый сигнал и затем он ослаблялся путем добавления к выборке новых последовательностей, уже не содержащих сигнала («мусорных последовательностей»), а кроме того - и путем «порчи» сайтов самого исходного сигнала. А именно, генерировались выборки из 10 бернуллиевских последовательностей каждая длиной 200 в четырехбуквенном алфавите {А, С, Т, G} и в каждую последовательность сначала подставлялось одно и тоже слово длины 16. Затем в каждом из вхождений этого слова случайным образом «портилось» несколько букв (имитация ослабления сигнала), а также добавлялись новые бернуллиевские последовательности, не содержащие сигнала (мусорные последовательности - имитация загрязнения выборки). Такой искусственный сайт считался найденным, если полученный в результате работы нашей программы сайт перекрывался с ним не менее, чем на половину его длины. Результаты таковы: сайты длиной 16 устойчиво находились при внесении в исходный сигнал до 3 независимых ошибок (в каждый из его сайтов для выборки из этих 10 последовательностей), а также — когда число мусорных последовательностей не превышало число всех 10 исходных последовательностей (табл. 1). При ошибках в 4 позициях исходного сигнала результат зависел от чистоты выборки: приемлемые результаты получались, когда число мусорных последовательностей не превосходило чисел 3-4. А именно, в большинстве испытаний искомые сайты правильно определялись практически во всех исходных последовательностях. При дальнейшем загрязнении выборки некоторые сайты в сигнале могли не обнаружиться, а доля таких результатов, естественно, повышалась с увеличением числа мусорных последовательностей. При ошибках в 5 позициях сайтов исходного сигнала менее пяти из этих сайтов в каждом сигнале обнаруживалось (табл. 1).
При втором способе ослабление сигнала достигалось за счет увеличения длин исходных последовательностей. Такое тестирование аналогично тому, которое выбрано в [60] для демонстрации качества представленных там алгоритмов. Там оно применялось к 8 исходным выборкам, каждая из которых содержит по 20 последовательностей длины п, а п меняется от 100 до 1000, с заранее имеющимися в них сигналами длиной 15 с 4 бернуллиевскими независимыми заменами в каждом сайте каждого сигнала. В [60] такие сигналы названы (15,4)-сигналами. В работе [60] ее авторы на этих 8 выборках тестировали ряд типовых программ и ряд их собственных программ для поиска оптимального сигнала с целью их сравнения между собой. Это были программы CONSENSUS, Gibbs sampler, МЕМЕ, WINNOWER, SP-STAR. В табл. 2 для всех этих программ приведены результаты из [60] и к ним добавлен результат такого же тестирования и нашей программы. В табл. 2 на пересечении строки и столбца приводится средний (по всем сайтам и всем выборкам) коэффициент нахождения сайта, где последний определен в [60] следующим образом. Если для данной последовательности обозначить К множество позиций исходного сайта и обозначить Р множество позиций сайта, предсказанного каким-то одним из перечисленных алгоритмов, то коэффициент нахождения сайта равен числу общих позиций у К и Р, деленному на число позиций в объединении множеств К и Р. Это тестирование позволяет сравнить эффективность нашей программы с другими наиболее употребительными программами. Из табл. 2 видно, что наша программа IRSA на выборках, предложенных в [60] для тестирования всех таких алгоритмов, находится на втором месте после алгоритмов WINNOWER и SP-STAR, предложенных самими авторами работы [60]. Отметим, что выборки, предложенные в [60], специально ориентированы на поиск именно (15, 4)-сигналов. Что касается второго места нашего алгоритма, то заметим, что для алгоритмов, занявших первое место на этих фиксированных выборках, известна только экспоненциальная верхняя оценка числа их шагов, а для нашего алгоритма нами получена полиномиальная верхняя оценка с низкими степенями вида пг-тъ-1г, где и — число последовательностей, т — длина последовательности, / - длина искомого сигнала.
В данном параграфе приводятся результаты тестирования нашей программы IRSA для поиска регуляторных сайтов на природных выборках регуляторных областей, которые постепенно портились. А именно, в качестве 3 исходных выборок были взяты регуляторные области перед генами (бактерии Escherichia coli), которые регулируются соответственно тремя белками-регуляторами PurR (пуриновый регулон), ArgR (аргининовый регулон), CRP (регулон катаболитной репрессии). Для каждой из трех выборок сигнал постепенно портился путем удаления из выборки по одному наилучшему2 из имеющихся в ней биологических сайтов. Таким образом, могли появляться мусорные последовательности и уменьшалось число сайтов в сигнале. Сайты удалялись таким образом до тех пор, пока их в общей сложности оставалось не менее 3 и пока среднее попарное сходство всех остающихся сайтов строго превышало число 1. Наш алгоритм IRSA искал сигнал с сайтами той же длины, что и у сайтов рассматриваемого биологического сигнала. Последовательности регуляторних сайтов были взяты из базы данных dpinteract [63], а фрагменты последовательностей, содержащие эти сайты, были извлечены из полного генома Е. coli при помощи программы GenomeExplorer [85].
Результаты тестирования оценивались с помощью двух функций S/ и &. Первая из них определяется как доля найденных биологических сайтов (в %) к общему числу таких сайтов, где биологический сайт считается найденным, если алгоритмически полученный сайт пересекается с ним не менее, чем на половину их общей длины. Вторая функция определяется как доля всех найденных сайтов (в %) к числу всех выданных алгоритмом сайтов. Перейдем к описанию результатов.
Пуриновий регулон. Здесь на вход алгоритма подавалась выборка регуляторных областей генов, регулируемых пуриновым репрессором PurR. Она состояла из 19 последовательностей каждая длиной 200 нуклеотидов и содержала в общей сложности 21 сайт длиной по 16 нуклеотидов. Две последовательности содержали по два сайта, остальные - по одному. Результаты таковы (табл. 3): даже если выборка более чем наполовину состояла из мусорных последовательностей, то больше половины остающихся сайтов опознавалось правильно в том смысле, что найденный нашей программой сайт и биологический сайт (одинаковой длины) совпадали не менее, чем на половину их длины. Когда в одной последовательности содержится два сайта, то после удаления первого из них второй находился правильно. Первые ошибки появляются при удалении 8 последовательно наилучших из этих сайтов.
