Содержание к диссертации
Введение
ГЛАВА 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]. Синтез молекул мРНК катализируется ферментами, которые называются РНК-полгтеразами. Молекулы мРНК одноцепочечные, относительно небольшой длины, соответствуют одному или нескольким белкам или их цепям. Последовательность нуклеотидов в молекуле мРНК
Трянскрктщия
Трансляция
мРНК
Б*лок>
Рисунок 1. Процесс синтеза белка.
считывается в соотвествии с генетическим кодом и этот процесс называется трансляцией.
Потребность клетки в некоторых белках значительно изменяется во времени, поэтому имеются механизмы регуляции, обеспечивающие изменение уровня синтеза белков в соответствие с потребностью в них. В частности, специальная группа белков контролирует синтез мРНК, регулируя таким образом концентрацию соответствующих ферментов. Такая регуляция может быть как положительной (тогда сам регулирующий белок называется активатором), так и отрицательной (репрессором). Аминокислотная последовательность белка-регулятора, как и любого другого белка, также кодируется в ДНК; определяющий ее ген называется геном-регулятором. Регуляторы специфичны, то есть каждый из них влияет на синтез какого-либо одного или нескольких определенных белков. В работе исследуется случай прокариотов хотя предлагаемый нами алгоритм применим и для гораздо более широкого класса задач. В простейших (прокариотических организмах) эта специфичность достигается специфичностью связывания молекулы белка-регулятора с определенными некодирующими участками молекул ДНК, расположенными непосредственно перед участком, кодирующим мРНК регулируемого набора ферментов. Специфические участки нуклеотидной последовательности, с которыми связываются регуляторные белки, называются сайтами. Много более длинный участок в ДНК, включающий эти сайты и расположенный перед кодирующим участком, называется лидерной областью или апстримом.
Общепринятое и подтвержденное предположение состоит в том, что все сайты
связывания одного белка достаточно сходны между собой. Это предположение
позволяет поставить задачу поиска набора сайтов связывания одного белка-регулятора в
исходном наборе родственных (относительно него) лидерных областей как задачу
нахождения набора наиболее сходных фрагментов в этой выборке лидерных областей.
ц Сам такой набор называется сигналом, а слова, входящие в него, естественно
называются сайтами (или иногда - потенциальными сайтами). Обычно структура и некоторые численные характеристики (например, длина) искомого сигнала заранее фиксируются или подбираются в ходе вычислений. Сигналу приписывается некоторое качество, которое тем выше, чем более попарно похожи друг на друга входящие в него сайты. Возможны разные точные определения качества сигнала. Саму задачу нахождения оптимального (наилучшего) по качеству сигнала в данном исходном наборе (выборке) родственных регуляторных областей называют задачей поиска оптимального сигнала. Она имеет некоторую связь с задачей множественного локального выравнивания, но, конечно, никак не сводится к ней. Существует подход, использующий метод поиска сигнала для построения выравнивания нескольких последовательностей [77].
Оценки качества сигнала и способы его описания
Существует несколько способов оценки качества полученного сигнала. Один из них - использование матрицы позиционных весов, элементы которой вычисляются по формуле [84]:
W(i>a)=A-{ln?ih2l~B0
(1),
где а є {А, С, Т, G}, С(г, а) - количество появлений нуклеотида а в позиции i,n-количество последовательностей. Коэффициенты Л и Ва подбираются таким образом,
зоновая
чтобы J2w(l>ayP« ~> а Т 2 Z)^20'Q)'P« -1» где Р«" Фс
t=l 4 a&{A,C,T,G} t=l
вероятность нуклеотида о, а / — общее число позиций в сайте. Фоновые вероятности ра букв исходного алфавита {Л, С, Т, 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 слов длины 6.
Другой способ состоит в том, что сигнал описывается матрицей выравнивания, каждый элемент иа( которой показывает число появлений каждой буквы а (из того же алфавита) в /-ой позиции сигнала (см. рис. 2). По ней строится вероятностная позиционная матрица сигнала:
f(a,i) =
%,і+с
(2)
Значения поправок са обычно выбираются так, чтобы
выполнялось У2 са = -у/п , где п - число последовательностей в выравнивании («
аЄ{АС,ї\С}
= 4 на рис. 2), а сами эти поправки были пропорциональны вероятностям ра появления букв в том материале, где ищется регуляторный сигнал [43]. Заметим, что
У^ f(a, і) — 1 в любом столбце (позиции) /.
at{A,T,C,G)
С помощью этой матрицы вычисляется информационное содержание данного набора сайтов по формуле:
i=l at\A,T,C,G} Ра
Величина информационного содержания (3) иногда используется как характеристика качества найденного сигнала, а вероятностная позиционная матрица -как решающее правило для поиска новых сайтов в исходных и новых регуляторных областях. Другой часто употребляемой характеристикой качества, применяемой как для представления (1), так и для (2), является качество соответствия сайтов, порождающих статистическую модель, этой модели [1,2, 3,25,43, 65,72,47].
Также для оценки сигнала в диссертационной работе используются сумма попарного сходства всех сайтов, входящих в сигнал, и среднее этого сходства.
F(.x*y) функция, отражающая степень сходства для двух слов х и у длины /, в данном случае, количество совпадающих букв в них.
S -сигнал длины /, s s^k^n) - сайты, входящие в сигнал S, п - количество
последовательностей.
Качество сайта s} q(Sj) = ^ F(sl}Sj).
i-i
Среднее качество сайта Sj в сигнале: p(sj) = q(Sj). Если в найденном сигнале
все слова одинаковы, то эта величина равна /.
*
Качество сигнала S - Q(S) = ^g(st)
(-і
Среднее качество сигнала S - P(S) =—yjTp(sl) Лучшим считается сигнал с наибольшим значением Р(5), а все сайты сигнала имеют качество p(s,).
Иногда удобно представлять сигнал в виде консенсуса - слова, состоящего из наиболее часто встечающихся букв в каждой позиции сигнала. Иногда два, реже три, нуклеотида встречаются одинаково часто, тогда в такие позиции ставится буква отличная от А, Т, С, G, означающая, что здесь может стоять один из нескольких нуклеотидов, Например, w означает, что в данной позиции может стоять А или Т. [11]
Основные алгоритмы поиска регуляторных сигналов
Известные подходы и алгоритмы выделения потенциальных регуляторных сигналов в исходном наборе (предполагаемых) регуляторных областей условно делятся на две группы: оптимизационные и комбинаторные. Некоторые алгоритмы сочетают в себе черты обеих групп. Оптимизационные алгоритмы основаны на некоторой характеристике качества сигнала (например, на его информационном содержании). Далее производится построение цепочки сигналов, так чтобы их качество (иногда говорят: значение функционала) постепенно возрастало. Таким образом, процедура сводится к поиску экстремума некоторого функционала в пространстве всех допустимых сигналов. Таковы алгоритмы максимизации ожидания МЕМЕ [1], стохастические и жадные алгоритмы Gibos sampler [23] и ряд других (например, имитация теплового отжига и DMS [36]).
Комбинаторные алгоритмы работают также с пространством сигналов, однако в этом случае цель состоит в построении специального слова (консенсуса), представленного во всех или во многих последовательностях из исходной выборки в том смысле, что искомые сайты отклонялись от него (и в этом смысле друг от друга) наименьшим образом (т.е. здесь также присутствует некоторая функция качества -какая-то мера компактности полученного набора сайтов, например его диаметр). К их
числу относятся CONSENSUS [28], PROJECTION [8], WINNOWER, SP-STAR [60], MITRA [13] и другие (например, Conslnd и Matlnd, ГГВ, WORDUP [59]).
Теперь рассмотрим подробнее несколько наиболее известных описанных в литературе алгоритмов распознавания регуляторньгх сигналов.
Метод максимизации оэюидания (ЕМ), Пусть g(j, к) - вероятность того, что искомый сигнал начинается в у позиции в к последовательности, ъ/(а, І) - вероятность появления буквы а в колонке /' для каждой из букв {А, С, G, Т} и 1 <і <І. В процессе работы алгоритма происходит переопределение g и/до тех пор, пока/не будет мало меняться. Это переопределение происходит с помощью правила байесовской вероятности.
v ' } Р(В) К '
В данной формуле в качестве А берется событие, состоящее в том, что сигнал в / последовательности начинается в j позиции, а в качестве В берется объединение двух событий:
Совокупность значений/( і) для всех і и or,
Что і последовательность равна данной последовательности S/.
Р(В\А) подсчитывается как произведение значений функции/по всем буквам слова длины /, начинающегося в последовательности St су-ой позиции.
Р(А) является начальной вероятностью того, что сигнал начнется в /-ой последовательности с у-ой позиции, и она равна 1/(т-1+1), где т — длина последовательности.
Р(В) вычисляется по формуле полной вероятности с перебором всех возможных начал сигнала в /-ой последовательности.
По формуле (4) переопределяется значения g(f, k)t которые позволяют нам определить новые значениями; і).
Ниже приведен псевдокод основного цикла алгоритма:
ЕМ (выборка, /){
Выбираем начальную точку/(случайно или определяется пользователем) do{
переопределяем g в соответствии с/
переопределяем/в соответствии с g } until (изменения/^є) return
}
Алгоритм EM находит модель сигнала.
Алгоритм МЕМЕ представляет собой расширение метода максимизации ожидания для поиска нескольких сигналов в одной выборке последовательностей. В качестве начальных точек берутся все слова длины / из выборки и запускается одна итерация алгоритма ЕМ. Каждый запуск ЕМ определяет одну вероятностную модель сигнала/и g, после чего выбирается лучшая и уже из начальной точки определенной этой моделью запускается алгоритма ЕМ до его сходимости. Получившийся сигнал запоминается, из выборки удаляются все вхождения этого сигнала, и запускается следующая итерация алгоритма МЕМЕ,
МЕМЕ (выборка, /, количество итераций к){ for І = 1 to к {
for каждого слова длины / из каждой последовательности выборки { запускаем 1 итерацию ЕМ с начальной точки определенной из этого слова
}
выбираем лучшую модель сигнала
запускаем ЕМ из начальной точки определенной этим сигналом запоминаем получившийся сигнал удаляем этот сигнал из выборки } }
Эвристическим аналогом ЕМ является довольно широко используемый алгоритм Gibbs sampling. Метод Gibbs Sampling был разработан Геманом и Геманом [23] для восстановления изображений. Впервые был применён для решения задачи МЛВ
Лоренсом и соавторами [43]. Различные модификации этого метода часто применяются для поиска слабых сигналов [43,37,72,47,29]. В самой простой версии мы ищем лучший консервативный неразрывный сигнал длины / в виде вероятностной позиционной матрицы (2). Мы считаем, что сигнал встречается во всех последовательностях.
Алгоритм итеративен. Сначала случайно выбираем одно слово длины / из каждой последовательности. Эти слова формируют начальное множество вхождений сигнала. Обозначим позицию слова в і-ой последовательность как О..
Итерационный шаг:
- берём одну последовательность /. Обычно, их выбирают по очереди, хотя
возможны и другие варианты, например, случайный выбор. Существенно, что у
всех последовательностей шансы быть выбранными равны.
Вычисляем вероятностную позиционную матрицу из выбранных слов всех последовательностей, кроме /. Обозначим эту матрицу как Р. Возьмем каждое слово длины / из /-ой последовательности и вычислим вес этого слова, относительно матрицы Р. Вес вычисляется как отношение вероятностей данного слова быть случайно порождённым из позиционной модели вероятностной позиционной матрицы и из фоновой модели.
- Разыграем следующее Ot- случайно из всех слов в і длины / используя
распределение вероятностей, определенное весами (вероятность — это вес,
нормированный на 1)
Заменяем О і на О,-.
Повторять итерационный шаг, пока ситуация не станет стабильной.
Основная идея алгоритма CONSENSUS заключается в том, что в матрицу выравнивания (описанную выше) по одному добавляют слова длины / из последовательностей, еще не включенных в матрицу. После каждого добавления
выбирается набор слов с наибольшим информационным содержанием. Алгоритм работает до тех пор, пока в матрице выравнивания не будет содержаться по одному слову из каждой последовательности или то количество, которое определяет пользователь, В качестве начального слова берутся по очереди все слова длины / из всех последовательностей или только из части входных данных, определенной пользователем. На рисунке 3 приведен пример поиска сигнала длины 4 в трех последовательностях длины 5.
Существует группа алгоритмов, решающих задачу поиска так называемого (/, d)-сигнала, т.е. сигнала длины /, который входит в каждую последовательность не более чем с d заменами. Сложность нахождения такого сигнала заключается в том, что два разных вхождения одного сигнала могут различаться в 2d позициях.
Алгоритм WINNOWER строит многодольный граф, вершинами которого являются слова длины /, долями - слова из одной последовательности, а ребрами связаны те вершины, которые отличаются друг от друга не более чем на 2d замен. Многодольный граф — это граф, в котором ребрами соединяются только вершины из разных долей. Сигналу соответствует клика в таком графе. Клика — это такой подграф, в котором любые две вершины соединены ребром. При поиске клики алгоритм сокращает число ребер в графе путем удаления тех, которые не могут быть частью никакой клики заданного размера.
Алгоритм SP-STAR отличается от WINNOWER тем, что строит взвешенный многодольный граф, в котором каждому ребру приписывается вес того, насколько отличаются слова соединенные этим ребром. Теперь задача заключается в поиске клики минимального веса, которая будет соответствовать искомому сигналу.
sequence 1 AC TG A
sequence 2 T AGCG
sequence 3 CTTGC
CYCLE 1
Aeq — 5.5
CYCLE 2
А С Т G Т Т G С
CYCLE 3
J.eq = 2.8
Л*, = 4.2
А С Т G AGCG С Т Т G
/
А С Т G AGCG
Т Т G С
*&eq
hea — 2.8
*«eq — >.*
/^, = 2.1
Рисунок 3. Пример работы алгоритма CONSENSUS.
Априорная вероятность каждого основания 25%.
Цикл 1. Для простоты, на этом цикле рассматривается матрица выравнивания для первого слова первой последовательности. На самом деле, создаются матрицы выравнивания для всех слов из всех последовательностей, если пользователь не определил подвыборку последовательностей в качестве начальной (это существенно, если выборка очень большая). Для этого примера на этом цикле создается 6 матриц выравнивания.
Цикл 2. Добавляя по слову из оставшихся двух последовательностей, получается 4 матрицы. Из них выбирается матрица с самым большим информационным содержанием (она обведена жирной линией).
Цикл 3. Только две матрицы создаются на этом цикле путем добавления двух слов из третьей последовательности, еще не вошедшей в матрицу. Жирной линией обведена матрица с наибольшим информационным содержанием, но она не обязательно является лучшей, потому что мы рассмотрели работу алгоритма только для первого слова из первой последовательности. Тоже самое надо проделать для всех слов, и после этого выбрать лучшую матрицу. (Конец подписи к рисунку 3)
Алгоритм МІТКА (Mismatch Tree Algorithm) использует префиксное дерево для разбиения всех возможных сигналов на непересекающиеся подмножества, имеющие общий префикс с точностью до (/замен. Это разбиение делит проблему поиска сигнала на более мелкие задачи, в которой можно применить уже известные алгоритмы, например WINNOWER (см. пример на рис. 4).
В диссертационной работе предложен и тестирован новый алгоритм для выделения сигнала в выборке нуклеотидных последовательностей. Он является промежуточным с точки зрения приведенной выше классификации - производится оптимизация функции качества, однако она определяется через попарное сходство слов, а не как информационное содержание набора.
Дано 3 последовательности длины 6.
ААТАТС
ATCCGT
GCATAC
Ищем (4, 1)-сигнал.
Имеем 9 слов длины 4 из трех последовательностей длины 6:
ААТА 4) АТСС 7) GCAT
АТАТ 5) TCCG 8) САТА
ТАТС 6) CCGT 9) АТАС
На каждом шаге все слова разделяются на части. В узлах дерева находятся слова, совпадающие с префиксом сигнала, с точностью до 1 замены. Мы ищем по одному слову из каждой последовательности, поэтому в процессе работы отсекаются ветки, не содержащие необходимого количества слов.
В данном случае сигналом является последовательность АТАС. Для простоты рассмотрены не все возможные варианты, а только одна ветвь дерева, соответствующая сигналам, начинающимся с А, поэтому этот сигнал может быть не единственным.
Рисунок 4. Пример работы алгоритма MITRA.
Алгоритм поиска выделения регул яторных сигналов белокового взаимодействия
Этап 2. Образуем вспомогательный граф G, который остается фиксированным в процессе работы алгоритма и задает порядок перебора исходных последовательностей. Граф G состоит из п вершин и всех ребер, которые возникают в процессе выполнения следующей процедуры, см. пример на рисунке б, в котором п=7. На первом шаге все вершины графа G разбиваются на две равные (с точностью до единицы, если п нечетное) части и между этими частями произвольным образом проводится ребро (А, В) (на рисунке 6 (А, В) =(1,2)), Далее итеративно повторяем такое разбиение «вглубь» графа G. А именно, каждую из двух полученных частей графа G снова разбиваем на две (в том же смысле) равные части. Относительно уже этих разбиений один конец ребра определен однозначно: это А для одной части и В для другой, а второй конец ребра выбирается произвольно. На рисунке 6 ребра второго уровня деления - это (1,3) и (2,4). И так далее: каждую появившуюся в этой процедуре не одновершинную часть Р разбиваем на две равные части Pi и Рг так, чтобы ребра новых частей выходили из концов ребра предыдущей части Р. Конечно, при этом каждое Р равно объединению его частей Pi и Р2. Процедура разбиений прекращается, когда все части Р станут одновершинными; на самом деле, можно остановиться, когда эти части станут просто мелкими (из 1-3 вершин).
ЭтапЗ. Внешний цикл алгоритма состоит во взаимно однозначном приписывании каждой вершине графа G одной из исходных последовательностей. Такое приписывание назовем расстановкой последовательностей по вершинам графа и обозначим г. Каждая следующая расстановка выбирается таким образом, чтобы как можно больше пар последовательностей, не соединенных на предыдущих итерациях, сейчас соединились. Заметим, что из вершин ребра, полученного при первом делении, выходит больше всего ребер, поэтому на следующей итерации в эти вершины ставятся последовательности, из которых, на текущий момент, меньше всего выходило ребер.
Условием выхода из цикла является момент, когда любая пара последовательностей хотя бы раз будет приписана вершинам графа G так, что окажется соединенной в нем ребром. Это обеспечивает разумное количество итераций внешнего цикла при достаточном разнообразии обработанных расстановок (порядка и) и позволяет избежать полного перебора. При каждой расстановке находится своя система слов, которые обрабатываются на этапе 6.
Этап 4. Для текущей расстановки г происходит поиск системы наиболее похожих слов, который мы называем сборкой. Граф G помогает организовать этот поиск. Отметим, что г (А) - это последовательность, приписанная вершине А, в которой, как и во всех других приписанных последовательностях, ищется сигнал, и далее мы будем говорить просто А, понимая, что это некоторая последовательность приписанная вершине А в соответствии с текущей расстановкой г. Процесс сборки является обратным к процессу деления графа. Теперь мы объединяем мелкие части и при этом находим систему наиболее похожих слов из последовательностей, находящихся в этих частях, выполняя индуктивный шаг, который состоит в следующем (рис. 7). Пусть для двух частей Pi и Рг с ребрами соответственно (А, С) и (В, D)t полученных разбиением подграфа Р (в исходном графе G) с ребром (А, В), уже определены два набора каждый из лучших систем. В первом наборе находятся слова из последовательностей, находящихся в части Pj, а во втором наборе - из / , назовем их набор над Р; и Рг, соответственно. Первый набор состоит из продолжений любых двух слов а и с из соответственно последовательностей А и С (где а и с такие, что значение F(a, с) больше некоторого фиксированного порога и) на множество Pi, т.е. любая пара (а, с) определяет лучшую систему слов из последовательностей части Р/. И аналогично, второй набор состоит из продолжений на все множество Рз слов бис/, произвольно взятых в последовательностях В и D. Объединяя таким образом части, мы получим набор над всем графом G, т.е. искомую систему слов (сигнал).
Например, на рисунке 6 пусть Pi = {1, 5, б, 3} и Рг = {2, 7, 4}; (А, С) = (I, 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 создается расстановка, т.е. вершинам графа приписываются входные последовательности до тех пор, пока каждая пара не будет хотя бы один раз соединена ребром. Для каждой расстановки происходит процедура сборки. И на последнем этапе происходит обработка результатов. Реализованы две возможности представления результатов, в виде статистики и в виде списка всех систем с пометкой самого лучшего (см. введение).
На искусственных выборках
Для определения характеристик алгоритма в контролируемых условиях было проведено его тестирование на искусственных выборках. Выборки были созданы двумя разными способами для различных целей.
При первом способе исходная выборка сначала содержала искомый сигнал и затем он ослаблялся путем добавления к выборке новых последовательностей, уже не содержащих сигнала («мусорных последовательностей»), а кроме того - и путем «порчи» сайтов самого исходного сигнала. А именно, генерировались выборки из 10 бернуллиевских последовательностей каждая длиной 200 в четырехбуквенном алфавите (А, С, Т, G} и в каждую последовательность сначала подставлялось одно и тоже слово длины 16. Затем в каждом из вхождений этого слова случайным образом «портилось» несколько букв (имитация ослабления сигнала), а также добавлялись новые бернуллиевские последовательности, не содержащие сигнала (мусорные последовательности - имитация загрязнения выборки). Такой искусственный сайт считался найденным, если полученный в результате работы нашей программы сайт перекрывался с ним не менее, чем на половину его длины. Результаты таковы: сайты длиной 16 устойчиво находились при внесении в исходный сигнал до 3 независимых ошибок (в каждый из его сайтов для выборки из этих 10 последовательностей), а также — когда число мусорных последовательностей не превышало число всех 10 исходных последовательностей (табл. 1). При ошибках в 4 позициях исходного сигнала результат зависел от чистоты выборки: приемлемые результаты получались, когда число мусорных последовательностей не превосходило чисел 3-4. А именно, в большинстве испытаний искомые сайты правильно определялись практически во всех исходных последовательностях. При дальнейшем загрязнении выборки некоторые сайты в сигнале могли не обнаружиться, а доля таких результатов, естественно, повышалась с увеличением числа мусорных последовательностей. При ошибках в 5 позициях сайтов исходного сигнала менее пяти из этих сайтов в каждом сигнале обнаруживалось (табл. 1).
Таблица 1. Результаты тестирования выборок из 10 исходных последовательностей. Номер строки указывает на число добавленных мусорных последовательностей (от 0 до 10). В заголовке столбца указано количество измененных букв во вхождениях исходного слова (от 0 до 5). На пересечении строки и столбца приведено число найденных исходных сайтов, где каждый знак соответствует отдельному независимому испытанию и знак Х=Ш. В первом столбце по 1 испытанию для 4-х случаев, во втором — по 10 испытаний для одного случая, в третьем - по 4 испытания для одного случая). В скобках указано среднее число найденных исходных сайтов в % для соответствующей серии испытаний.
При втором способе ослабление сигнала достигалось за счет увеличения длин исходных последовательностей. Такое тестирование аналогично тому, которое выбрано в [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)-сигналов. Что касается второго места нашего алгоритма, то заметим, что для алгоритмов, занявших первое место на этих фиксированных выборках, известна только экспоненциальная верхняя оценка числа их шагов, а для нашего алгоритма нами получена полиномиальная верхняя оценка с низкими степенями вида п2 п3-13, где п — число последовательностей, т — длина последовательности, / - длина искомого сигнала.
В данном параграфе приводятся результаты тестирования нашей программы IRSA для поиска регуляторных сайтов на природных выборках регуляториых областей, которые постепенно портились. А именно, в качестве 3 исходных выборок были взяты регуляторные области перед генами (бактерии Escherichia соїг), которые регулируются соответственно тремя белками-регуляторами PurR (пуриновий регулон), ArgR (аргининовый регулон), CRP (регулон катаболитной репрессии). Для каждой из трех выборок сигнал постепенно портился путем удаления из выборки по одному наилучшему2 из имеющихся в ней биологических сайтов. Таким образом, могли появляться мусорные последовательности и уменьшалось число сайтов в сигнале. Сайты удалялись таким образом до тех пор, пока их в общей сложности оставалось не менее 3 и пока среднее попарное сходство всех остающихся сайтов строго превышало число 1. Наш алгоритм IRSA искал сигнал с сайтами той же длины, что и у сайтов рассматриваемого биологического сигнала.
Последовательности регуляторных сайтов были взяты из базы данных dpinteract [63], а фрагменты последовательностей, содержащие эти сайты, были извлечены из полного генома Е. coli при помощи программы GenomeExplorer [85].
Результаты тестирования оценивались с помощью двух функций .S/и /,. Первая из них определяется как доля найденных биологических сайтов (в %) к общему числу таких сайтов, где биологический сайт считается найденным, если алгоритмически полученный сайт пересекается с ним не менее, чем на половину их общей длины. Вторая функция определяется как доля всех найденных сайтов (в %) к числу всех выданных алгоритмом сайтов. Перейдем к описанию результатов.
Применение программы для поиска потенциальных сигналов связывания транскрипционных факторов в ортологичных рядах генов организмов из групп ентеробактерии и бациллы/клостридии
Сравнительный подход также применяется для анализа регулонов, отвечающих за метаболизм ферментов необходимых для жизнедеятельности организма, В данном случае изучался GlpR-регулон, отвечающий за метаболизм глицерола и глицерол 3-фосфата (ГЗФ) в геномах а-, (Ї- и у-протеобактерий.
Белок GlpR, принадлежащий к семейству регуляторов DeoR, является репрессором и контролирует экспрессию генов метаболизма глицерола и ГЗФ. GlpR-регулон хорошо изучен в Escherichia coli [78,42,80] и частично охарактеризован в Pseudomonas aeruginosa [67].
Глицерол поступает извне в цитоплазму путем облегченной диффузии (см. рис. 9), обеспечиваемой продуктом гена glpF, а ГЗФ активно транспортируется продуктом гена glpT. Внутриклеточный глицерол фосфорилируется глицеролкиназой (glpK), давая ГЗФ. ГЗФ затем может быть превращен в дигидроксиацетонфосфат под действием одной из двух имеющихся у Е. coli ГЗФ дегидрогеназ: аэробной (glpD) или анаэробной (glpA). Кроме того, к GlpR регулону Е. coli относится ген glpQ, кодирующий периплазматическую глицерофосфодиэстеразу, гидролизующую глицерофосфодиэфиры с высвобождением ГЗФ, гены glpB и gipC, кодирующие дополнительные структурные компоненты анаэробной ГЗФ дегидрогеназы, а также гены glpE, glpG и glpX, функции которых не ясны. Вышеназванные гены собраны в три локуса на хромосоме Е. coli: glpTQlglpABC, glpEGRIglpD и glpFKX (/ разделяет опероны, ориентированные в разные стороны). Близкие гомологи GlpR были обнаружены во многих геномах а-, р- и у-протеобактерий и построено дерево (рис. 10). Цель этой главы - поиск сайтов связывания белка GlpR. Для этого нами проведен дополнительный анализ гомологии GlpR-регулируемых генов и определена их оперонная структура в ряде геномов (рис. 11). В тех случаях, когда оперонная структура данного фрагмента ДНК неизвестна, гены относились к одному потенциальному оперону, если они имели одинаковое направления считывания, а расстояние между ними не превышало 100 п.н. За начало такого потенциального оперона принимался ген, перед которым был найден GlpR-сайт. Обозначения генов соответствуют их ортологам в Е. coll Были рассмотрены следующие геномы. у-протеобактерии: Escherichia coli [5], Salmonella typhi [58], S. typhimurium [53], Klebsiella pneumoniae [31], Erwinia carotovora, Yersinia pestis [34], Y. enterocolitica [34], Vibrio cholerae [57], V. vulnificus [34], V.fischeri, Pasteurella multocida [26], P. haemolytica, Haemophilus influenzae [52], H. ducrey, H. somnus, Pseudomonas aeruginosa [15], P.fluorescens [32], P.putida [34], P. syringae [35], Actinobacillus actinomycetemcomitans [33], Xanthomonas axonopodis, X. campestris; Р-протеобактерии: Burkholderiafungorum, B. pseudomallei, B. cepacia; а-лротеобактерии: Bordetella parapertussis, Ralstonia eutropha, R. solanacearum [34], Mesorhizobium loti [34], Sinorhizobium meliloti [34], Rhizobium hguminosarum, Agrobacterium tumefaciens [34], Rhodopseudomonas palustris, Brucella melitensis [34], Rhodobacter sphaeroides. Для выравнивания последовательностей белков и построения филогенетического дерева нами использовались соответственно программы ClustalW [73] и Phylip [45], Основные таксономические группы, соответствующие ветвям дерева белков GlpR, были рассмотрены отдельно. Области перед генами GIpR-регулона составили файл (свой для каждой таксономической группы), к которому применялась программа IRSA, представленная в диссертационной работе. В результате работы программы находился сигнал - набор потенциальных сайтов связывания GlpR-penpeccopa в этих последовательностях. Таким образом полученный сигнал назовем базисным (в табл. 12 он отмечен жирным шрифтом). Уже по нему строилась матрица позиционных весов (1) (таким образом, он служил обучающей выборкой). Для ее построения использовалась программа SignalX [85], а для сканирования геномов - программа GenomeExpIorer [85]. Все матрицы позиционных весов применялись к геномам для поиска новых сайтов в областях —400.,.+50 п.н. относительно начала гена. Выборки, соответствующие базисному сигналу, для нескольких таксономических групп с предсказанными сайтами представлены в приложении 1. Полученные результаты более подробно мы приведем ниже по группам организмов. у-Протеобактерии. семейство Enterobacteriaceae. Сначала рассматривались регуляторные области из четырех геномов Е. coli, Е. carotovora, К enterocolitica, К. pneumoniae, и по ним был получен базисный сигнал, включающий уже известные сайты Е. соїі (табл. 12а) с консенсусом TGTTCGATAACGAACA (рис. 12а), По базисному сигналу как по обучающей выборке была построена матрица позиционных весов (табл. 13 а) для поиска палиндромных сайтов длины 16. С помощью этой матрицы с порогом 4,1 были еще найдены сайты в дополнительных геномах У. pestis, S. typhlmurium, S. typhi (табл. 12a).
Применение программы для исследования регуляции метаболизма глицерол-3-фосфата
Эвристическим аналогом ЕМ является довольно широко используемый алгоритм Gibbs sampling. Метод Gibbs Sampling был разработан Геманом и Геманом [23] для восстановления изображений. Впервые был применён для решения задачи МЛВ Лоренсом и соавторами [43]. Различные модификации этого метода часто применяются для поиска слабых сигналов [43,37,72,47,29]. В самой простой версии мы ищем лучший консервативный неразрывный сигнал длины / в виде вероятностной позиционной матрицы (2). Мы считаем, что сигнал встречается во всех последовательностях.
Алгоритм итеративен. Сначала случайно выбираем одно слово длины / из каждой последовательности. Эти слова формируют начальное множество вхождений сигнала. Обозначим позицию слова в і-ой последовательность как О.. Итерационный шаг: - берём одну последовательность /. Обычно, их выбирают по очереди, хотя возможны и другие варианты, например, случайный выбор. Существенно, что у всех последовательностей шансы быть выбранными равны. Вычисляем вероятностную позиционную матрицу из выбранных слов всех последовательностей, кроме /. Обозначим эту матрицу как Р. Возьмем каждое слово длины / из /-ой последовательности и вычислим вес этого слова, относительно матрицы Р. Вес вычисляется как отношение вероятностей данного слова быть случайно порождённым из позиционной модели вероятностной позиционной матрицы и из фоновой модели. - Разыграем следующее Ot- случайно из всех слов в і длины / используя распределение вероятностей, определенное весами (вероятность — это вес, нормированный на 1) Заменяем О І на О,-. Повторять итерационный шаг, пока ситуация не станет стабильной. Основная идея алгоритма CONSENSUS заключается в том, что в матрицу выравнивания (описанную выше) по одному добавляют слова длины / из последовательностей, еще не включенных в матрицу. После каждого добавления выбирается набор слов с наибольшим информационным содержанием. Алгоритм работает до тех пор, пока в матрице выравнивания не будет содержаться по одному слову из каждой последовательности или то количество, которое определяет пользователь, В качестве начального слова берутся по очереди все слова длины / из всех последовательностей или только из части входных данных, определенной пользователем. На рисунке 3 приведен пример поиска сигнала длины 4 в трех последовательностях длины 5.
Существует группа алгоритмов, решающих задачу поиска так называемого (/, d)-сигнала, т.е. сигнала длины /, который входит в каждую последовательность не более чем с d заменами. Сложность нахождения такого сигнала заключается в том, что два разных вхождения одного сигнала могут различаться в 2d позициях.
Алгоритм WINNOWER строит многодольный граф, вершинами которого являются слова длины /, долями - слова из одной последовательности, а ребрами связаны те вершины, которые отличаются друг от друга не более чем на 2d замен. Многодольный граф — это граф, в котором ребрами соединяются только вершины из разных долей. Сигналу соответствует клика в таком графе. Клика — это такой подграф, в котором любые две вершины соединены ребром. При поиске клики алгоритм сокращает число ребер в графе путем удаления тех, которые не могут быть частью никакой клики заданного размера.
Алгоритм SP-STAR отличается от WINNOWER тем, что строит взвешенный многодольный граф, в котором каждому ребру приписывается вес того, насколько отличаются слова соединенные этим ребром. Теперь задача заключается в поиске клики минимального веса, которая будет соответствовать искомому сигналу.
Алгоритм МІТКА (Mismatch Tree Algorithm) использует префиксное дерево для разбиения всех возможных сигналов на непересекающиеся подмножества, имеющие общий префикс с точностью до (/замен. Это разбиение делит проблему поиска сигнала на более мелкие задачи, в которой можно применить уже известные алгоритмы, например WINNOWER (см. пример на рис. 4).
В диссертационной работе предложен и тестирован новый алгоритм для выделения сигнала в выборке нуклеотидных последовательностей. Он является промежуточным с точки зрения приведенной выше классификации - производится оптимизация функции качества, однако она определяется через попарное сходство слов, а не как информационное содержание набора.
