Содержание к диссертации
Введение
Глава 1. Биологические задачи: методологические различия физического и вычислительного модельных подходов 23
Модель. Понятие 23
Идеальность модели. Универсальность. Верифицируемость
Вычислительный метод как часть модели
Построение модели
Мотив как модель ДНК-белкового взаимодействия
Глава 2. Алгоритм поиска и определения длины и структуры консервативных участков ДНК, программная реализация 29
Представление мотива 30
Основной шаг поиска сильного мотива 31
Поиск длины и наилучшего сдвига множества найденных сайтов 33
Подробности организация алгоритма 36
Стадии работы программы , 36
Программная реализация алгоритма и важнейшие параметры поиска мотива 40
Интерфейс командной строки для запуска программы 41
Входи выход программы 41
Соглашение о представлении позиций 42
Ключи командной строки и метки конфигурационного файла , 42
Веб-сайт, реализующий внешний интерфейс к программе 56
Глава 3. Сравнение программы SeSiMCMC с аналогичными 58
Глава 4. Применение алгоритма для поиска участков связывания белков - регуляторов
дыхания АгсА и NarP на ДНК в бактерии Е. coli и описания соответствующих регулонов 60
Переключатель кислородного/бескислородного дыхания Phospho-ArcA 61
Регулятор азотного дыхания NarP 64
Выводы 76
Литература 77
- Идеальность модели. Универсальность. Верифицируемость
- Основной шаг поиска сильного мотива
- Программная реализация алгоритма и важнейшие параметры поиска мотива
- Сравнение программы SeSiMCMC с аналогичными
Введение к работе
Публикация почти завершённого человеческого генома — это удивительное достижение, но описание всех функциональных элементов, закодированных в человеческом и других геномах, остается неразрешённой задачей огромной сложности (Collins, Green et al. 2003). Директор Американского Национального Института Исследований генома человека Франсис Коллинз предполагает, что "... следующей стадией геномики будет каталогизация, описание и понимание всего множества функциональных элементов (включая и те, которые не кодируют белки), зашифрованных в человеческом и других геномах." (Collins, Green et al. 2003). Два из самых важных видов функциональных элементов генома — это гены, кодирующие транскрипционные факторы (ТФ, TF) и те участки ДНК (сайты связывания), к которым прикрепляются их молекулы. Взаимодействия транскрипционных факторов с ДНК регулируют в клетках множество важных процессов, например изменения характера развития или ответы на воздействия внешней среды, и дефекты в этих взаимодействиях могут приводить к развитию различных болезней, В последнее время, большой прогресс достигнут в сборе и анализе транскрипционных профилей мРНК многих типов тканей и клеток, в том числе тех типов, которые связаны с различными болезнями человека (Lockhart, Winzeler 2000); тем не менее, многое в цепях регуляции транскрипции, порождающих это профили, остаётся непонятным. Лучшее понимание структуры и функционирования транскрипционных факторов, соответствующих сайтов связывания на ДНК и их взаимодействий, могло бы привести к более широкому и при этом количественному описанию регуляторных каскадов в клетках, а также более глубокому понимание возможных функций индивидуальных генов, регулируемых вновь обнаруженными сайтами связывания ТФ (Gelfand 1999; Gelfand, Koonin et al. 2000; McGuire, Hughes et al. 2000; Bulyk 2003; Герасимова, Гельфанд et al. 2003; Liu, De Wulf 2004). (a)TACGAT TATAAT TATAAT GATACT TATGAT TATGTT (6)TATAAT (в)TATRNT (r)nos.l 2 3 4 5 6 A 0 6 0 3 4 0
С 0 0 1 0 1 0 G 1 0 0 3 0 0 T 5 0 5 0 1 6
Рисунок 1. Различные представления множества ССТФ: (а) список известных сайтов (б) консенсус (в) вырожденный консенсус (строка-шаблон) (г) матрица выравнивания
Лишь для немногих транскрипционных факторов хорошо описана их специфичность связывания. Сайты связывания ТФ (ССТФ, TFBS) обычно коротки (5-15 пар оснований, Ьр), и, как правило, представляют собой неточные (вырожденные) копии одного слова (рисунок 4 а, б, в); потенциальные сайты связывания могут очень часто встречаться в геномах. Вырожденность последовательности ССТФ - это следствие эволюционного отбора. Она полезна, поскольку даёт принципиальную возможность разным промоторам, регулируемым одним ТФ, проявлять разную активность, и тем самым с разной активностью транскрибировать разные гены, если это нужно для жизни клетки (Stormo 2000). Сайты связывания часто могут быть ориентированы в любую сторону. У прокариотов и у дрожжей, ССТФ находятся выше (раньше) точки начала регулируемого гена или серии вместе транскрибируемых генов (оперона), т.е в промоторной области. У высших эукариотов они могут быть выше гена, ниже его или находится внутри интрона, к тому же, они могут располагаться как вблизи, так и вдали от регулируемого гена.
Большая доля генома не кодирует белков, и искать ССТФ в некодирующей его части с помощью простых инструментов поиска участка последовательности, таких, как BLASTN или CLUSTALW (Cliften, Hillier et al. 2001) - это трудоёмкая и малоэффективная методика
Экспериментальные методы нахождения ССТФ.
Много информации о специфичности связывания ТФ было получено с помощью традиционных экспериментальных методик, таких как футпринтинг (обнаруживает участок ДНК, защищенный связанным белком - фактором), проверка связывания белка на нитроцеллюлозе, анализ уменьшения подвижности в геле (отслеживает изменения подвижности при связывании ДНК с белком), Саузерн блопинг ДНК и белка, репортёрные конструкции, а также точечный мутагенез. Тем не менее, все эти методы требуют больших временных затрат и не масштабируются на полные геномы. В последние годы, был разработан ряд высокопроизводительных технологий, позволяющих идентифицировать ССТФ in vivo и i« vitro. Один из таких in vitro методов - это рандомизация двух цепей исходной ДНК с последующей выборкой тех участков, которые хорошо связываются с белком (SELEX, systematic evolution of ligands by exponential evolution) (Oliphant, Brandl et al. 1989). Этот метод был также модифицирован в геномный SELEX, который использует геномные библиотеки как начальный материал для модификации и поиска (Gold, Brown et al. 1997). Затем, специфичность ДНК-связывающих белков стали определять непосредственно по связыванию белков с олигонуклеотидными матрицами двунаправленных ДНК (Bulyk, Gentalen et al. 1999; Bulyk, Huang etal.2001).
Так же, в мире разрабатывались высокопроизводительные методы измерений ДНК-ТФ взаимодействий in vivo. Самый широко применяемый из них в настоящее время (см. обзор Wyrick, Young
2002)- это иммунное осаждение с последующим анализом ДНК микрочипами CChlP-chip'), также известный как полногеномный анализ расположения сайтов связывания (Ren, Robert et al. 2000). Этот подход позволил описать функционирование ряда ТФ в дрожжах Saecharomyces cerevisiae (Reid, Iyer et al.
2000; Ren, Robert et al. 2000; Iyer, Horak et al. 2001; Lieb, Liu et al. 2001; Simon, Bamett et al. 2001; Lee, Rinaldi et al. 2002) и, позже, в клетках млекопитающих (Horak, Mahajan et al. 2002; Ren, Cam et al. 2002; Weinmann, Yan et al. 2002). Ещё эффективнее его сочетание с анализом данных экпрессии (Xu, Wang et al. 2004). Ещё один недавно появившийся метод основанный на применении ДНК-матриц для идентификации ССТФ in vivo, использует ТФ, сшитые с ДНК-адениметил-трансфераэой (Dam), которая в результате метилирует ДНК около сайтов связывания интегрального ТФ-Dam протеина (van Steensel, Henikoff2000; van Steensel, Delrow et al. 2001). Этим подходом были in vivo идентифицированы сайты в Drosophila (van Steensel, Delrow et al. 2001; van Steensel, Delrow et al. 2003) и vArabidopsis (Tompa, McCallum et al. 2002).
Поиск сайтов связывания транскрипционных факторов in silico.
Основой для алгоритмического поиска сайта связывания ТФ ранее неизвестного вида (структуры, консенсуса) является общепринятое и проверенное на практике предположение, что сайты связывания одного белка схожи между собой. Это предположение превращает задачу поиска набора сайтов, узнаваемых одним транскрипционным фактором, к задаче нахождения множества сходных фрагментов (мотива) в выборке (некодирующих) последовательностей ДНК. Если рассмотреть ту же задачу с точки зрения взаимного сходства последовательностей, содержащих сайты, она превращается в задачу множественного локального выравнивания (МЛВ) (Waterman 1986).
Описания мотива и оценки его качества.
Основа алгоритма поиска мотива (сигнала) — это представление этого мотива, или его вычислительная модель (этот аспект мы обсудим ниже, в главе 2). Кроме того, поиск наилучшего сигнала предполагает численное выражение качества сигнала. Это требование не столь тривиально, как кажется на первый взгляд. Непросто, например, сразу сказать, что лучше: слово длины 8, повторяющееся в некоторой последовательности 5 раз, или длины 10-3 раза.
Понятие "лучше" может содержать не только качество мотива как такового, но и его соответствие той выборке ДНК, по которой проводится поиск. Кроме того, численное представление качества, естественно, зависит от способа представления самого мотива. Рассмотрим подробнее некоторые из них.
Матрица позиционных весов,
Матрица позиционных весов (МПВ) для описания мотива приписывает каждому нуклеотиду (символу нуклеодитного алфавита) его "встречаемость", или вес, в каждой позиции мотива. Размер этой матрицы, таким образом, равен 4xs, где s - длина мотива. Математическое определение веса (элемента матрицы) может быть различным, но так или иначе, все представления МПВ монотонны друг по отношению к другу, и чем чаще встречается нуклеотид в данной позиции мотива, тем больше соответствующий вес. Простейший вариант МПВ - это матрица выравнивания С (і, а) (рисунок 4, г), считающая вхождения нуклеотида а в позиции і всех сайтов, формирующих мотив. Поскольку столбцы матрицы в простейшем случае независимы, соответствующие буква веса при подсчёте веса (качества) слова перемножаются. Поэтому удобно использовать не саму матрицу выравнивания, а её поэлементный логарифм. Кроме того, удобно работать с такой матрицей чтобы веса случайных слов, подсчитанные по ней, были распределены со средним 0 и дисперсией 1. Таким образом, получается (Миронов, Гельфанд 1999) следующая формула для веса слова: W(i,a)=K- ta^-*,) (0.1), где а є [А, С, Т, G}, С (і,а) - количество появлений нуклеотида а в позиции /. N-количество последовательностей. Коэффициенты KviBj подбираются таким образом, чтобы Z W(і,а) ра = 0, а ZW2 (і,а) ра = 1, где ра- фоновая вероятность нуклеотида а.
1=1 г=1
Также сигнал может описываться матрицей выравнивания, каждый элемент я, 0 которой показывает число появлений каждой буквы а (из набора А, Т, С, G) в /-ой позиции сигнала, (см. рис. 1). По ней строится вероятностная позиционная матрица сигнала: f(Ua)= > + С' (0.2)
1=1 Значения поправок с,, (pseudocounts) обычно берутся так, чтобы Zc/ = 4n , где N - количество последовательностей в выравнивании (4 на рис.1), а сами эти поправки были пропорциональны вероятностям ра появления букв в том материале, в котором мы ищем регуляторный сигнал (фоновым вероятностям) (см. например, Lawrence, Altschul et al. 1993).
Заметим, что Z f(i,&) = 1 в любом столбце (позиции) /. ек{Л,г,с,0}
По зтойматрице вычисляется информационное содержание данного набора сайтов по формуле Шеннона (Shannon, Weaver 1949):
4,, =Z I /(*,«)-ln/M
1=1 as{A,T,C,G) или его энтропийное расстояние (или условное информационное содержание) от фонового распределения частот по формуле Кульбака-Лейбера (Kullback 1997): '-,=t Z /0»-ln^^ (0.3). i-1 а<=[А,Т,С,0} Pa
Обычно, под информационным содержанием сигнала имеют в виду величину условного информационного содержания (0.3). Эта величина может эффективно использоваться как характеристика качества найденного сигнала (Favorov, Geltand et al. 2005). Другой часто употребляемой характеристикой качества, применяемой как для представления (0.1), так и для (0.2), является качество соответствия сайтов, порождающих статистическую модель, этой модели (Lawrence, Altschul et al. 1993; Bailey, Elkan 1994; Bailey, Elkan 1995; Bailey, Elkan 1995; Grundy, Bailey et al. 1996; Roth, Hughes et al. 1998; Liu, Brutlag et al. 2001; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005).
В качестве фоновых вероятностей букв можно брать частоту букв всех последовательностей организма (т.е. геномная частота нуклеотидных оснований) или частоты букв данной выборки последовательностей, или использовать заранее известные частоты, характеризующие данный генетический материал. Используются и более сложные, например, Марковские, статистические модели для нерегуляторного генетического материала (Liu, Brutlag et al. 2001; Thijs, Marchal et al. 2002).
Различные типы алгоритмов для поиска мотивов
Оба приведённые выше представления, (0.1) и (0.2), являются статистическими в том смысле, что в их основе лежат статистики встречаемости букв алфавита в позициях сайтов. Возможны и другие описания сигнала, например, грамматические (алгебраические), в которых сигнал выглядит как некоторое грамматическое выражение, которому должны соответствовать сайты. Такой тип описания представляет, например, (/,
В соответствие с типом представления сигнала, можно классифицировать алгоритмы поиска. Так, использующие статистическое представление (разного рода статистические модели порождения) мы относим к статистическим алгоритмам, а базирующиеся на алгебраической организации сигнала (строки, строки с шаблонами, графы соответствий) - к алгебраическим (грамматическим). Каждый возможный тип представления сигнала порождает свои алгоритмы для его поиска.
Задача поиска сигнала как множества наиболее сходных фрагментов - это всегда задача оптимизации. Оптимизируется и качество набора составляющих сигнал сайтов, и параметры текущей модели (шаблона), описывающей сигнал в представлении данного алгоритма. Параметры модели и набор слов связаны описанием модели, с одной стороны, и теми последовательностями ДНК, в которых происходит поиск, с другой. Целевая функция, которую непосредственно оптимизирует алгоритм, оптимизировать может быть как функцией от набора сайтов, так и от параметров модели. В первом случае алгоритм можно назвать комбинаторным, во втором - оптимизационным.
Примеры комбинаторных алгоритмов: Gibbs sampler (Lawrence, Altschul et al. 1993; Roth, Hughes et al. 1998; Liu, Brutlag et al. 2001; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005), Motif Sampler из (Rouchka 1997), MEME (Bailey, Elkan 1994; Bailey, Elkan 1995; Bailey, Elkan 1995; Grundy, Bailey et al. 1996), Conslnd и Matlnd (Freeh, Herrmann et al. 1993; Quandt, Freeh et al. 1995; Wolfertstetter, Freeh et al. 1996); ITB (Kielbasa, Korbel et al. 2001); WORDUP (Pesole, Prunella et al. 1992; Liuni, Prunella et al. 1993); CONSENSUS (Hertz, Stotmo 1999); WINNOWER, SP-STAR (Pevzner, Sze 2000), PROJECTION (Buhler, Tompa 2002); суффиксные и префиксные деревья (Jonassen 1997; Marsan, Sagot 2000); MITRA (Eskin, Pevzner 2002); а также другие алгоритмы (Fraenkel, Mandel et al. 1995; Ulyanov, Stormo 1995; Cho, Campbell et al. 1998; Rigoutsos, Floratos 1998; Rocke, Tompa 1998; van Helden, Andre et al. 1998; Tompa 1999; Wolfsberg, Gabrielian et al. 1999; Jensen, Knudsen 2000).
Оптимизационные алгоритмы - это CONSENSUS (Stormo, Hartzell 1989; Hertz, Hartzell et al. 1990; Hertz, Stormo 1999), алгоритмы максимизации ожидания (Lawrence, Reilly 1990; Cardon, Stormo 1992; Frishman, Mironov et al. 1999; Gelfatid, Koonin et al. 2000), DMS (Hu, Sandmeyer et al. 2000), Site Sampler из (Rouchka 1997).
Для многих алгоритмов целью является построить сигнал, представленный в каждой последовательности выборки с наименьшим отклонением, т.е. функцией качества является какая-то мера компактности полученного набора, например, его диаметр. К числу таких алгоритмов относятся Conslnd и Matlnd (Freeh, Herrmann et al, 1993; Quandt, Freeh et al. 1995); ITB (Kielbasa, Korbel et al. 2001); WORDUP (Pesole, Prunella et al. 1992; Liimi, Prunella et al. 1993); CONSENSUS (Hertz, Stormo 1999); WINNOWER, SP-STAR (Pevzner, Sze 2000), PROJECTION (Buhler, Tompa 2002); суффиксные и префиксные деревья (Jonassen 1997; Marsan, Sagot 2000); MITRA (Eskin, Pevzner 2002); а также другие алгоритмы (Fraenkel, Mandel et al. 1995; Ulyanov, Stormo 1995; Cho, Campbell et al. 1998; Rigoutsos, Floratos 1998; Rocke, Tompa 1998; van Helden, Andre et al. 1998; Tompa 1999; Wolfsberg, Gabrielian et al. 1999; Jensen, Knudsen 2000).
Другие алгоритмы оптимизируют вероятностные характеристики качества набора потенциальных сайтов (например, информационное содержание). Так работают жадные алгоритмы (Stormo, Hartzell 1989; Hertz, Hartzell et al. 1990), алгоритмы максимизации ожидания (Lawrence, Reilly 1990; Cardon, Stormo 1992; Frishman, Mironov et al. 1999; Gelfand, Koonin et al, 2000), DMS (Hu, Sandmeyer et al. 2000), MEME (Bailey, Elkan 1994; Bailey, Elkan 1995; Bailey, Elkan 1995; Grundy, Bailey et al. 1996); а также стохастические алгоритмы: имитация теплового отжига (Lukashin, Engelbrecht et al, 1992); Gibbs sampler (Lawrence, Altschul et al. 1993; Roth, Hughes et al. 1998; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005).
Алгоритмы поиска регуляторного сигнала, как и любые другие оптимизационные алгоритмы, бывают детерминированными (например, градиентные, или жадные алгоритмы) и стохастическими (эвристическими). Примеры жадных алгоритмов описаны в (Stormo, Hartzell 1989; Hert2, Hartzell et al. 1990; Bailey, Elkan 1994; Bailey, Elkan 1995; Bailey, Elkan 1995; Grundy, Bailey et al. 1996), а эвристических - в (Lawrence, Altschul et al. 1993; Roth, Hughes et al. 1998; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005). Первый тип даёт более проверяемые результаты, и, как правило, вычислительно выгоднее, но не может гарантировать глобальность найденного максимума. Второй тип сходится к глобальному максимуму, хотя бы в статистическом смысле.
Описание нескольких наиболее известных алгоритмов поиска мотивов
Теперь рассмотрим подробнее несколько наиболее известных описанных в литературе алгоритмов распознавания регуляторных сигналов,
Метод максимизации ожидания (ЕМ). Пусть gjk - вероятность того, что искомый сигнал начинается в у позиции в последовательности, a f(i, а) - вероятность буквы а в колонке і для каждой буквы из алфавита и 1
Выбираем начальную точку /(случайно или определяется пользователем) do{ переопределяем g в соответствии с/ переопределяем /в соответствии с g } until (изменения/) return
Описанный алгоритм ЕМ (Bailey, Elkan 1994) находит модель сигнала в виде ПВМ. Алгоритм МЕМЕ представляет собой расширение метода максимизации ожидания для поиска нескольких сигналов в одной выборке последовательностей. В качестве начальных точек берутся все слова длины / из выборки и запускается одна итерация алгоритма ЕМ. Каждый запуск ЕМ определяет одну вероятностную модель сигнала, после чего выбирается лучшая и уже из начальной точки определенной этой моделью запускается алгоритма ЕМ до сходимости. Получившийся сигнал запоминается, из выборки удаляются все вхождения этого сигнала, и запускается следующая итерация алгоритма МЕМЕ. МЕМЕ (выборка, /, количество итераций к){ for і = 1 to к { for каждого слова длины / из каждой последовательности выборки { запускаем 1 итерацию ЕМ с начальной точки определенной из этого слова выбираем лучшую модель сигнала запускаем ЕМ из начальной точки определенной этим сигналом запоминаем получившийся сигнал удаляем этот сигнал из выборки Gibbs sampling представляет довольно широко используемый эвристический аналог ЕМ. Метод Gibbs Sampling был разработан Геманом и Геманом (Geman, Geman 1984) для восстановления изображений. Впервые был применён для решения задачи МЛВ Лоренсом и соавторами (Lawrence, Altschul et al. 1993). Различные модификации этого метода часто применяются для поиска слабых сигналов (Lawrence, Altschul et al. 1993; Rouchka 1997; Hughes, Estep et al. 2000; Liu, Brutlag et al. 2001; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005). В самой простой версии мы ищем лучший консервативный неразрывный сигнал длины / в виде вероятносной позиционной матрицы (ВПМ) (0.2). Мы считаем, что сигнал встречается во всех последовательностях.
Алгоритм итеративен. Сначала случайно выбираем одно слово длины / из каждой последовательности. Эти слова формируют начальное множество вхождений сигнала. Обозначим позицию слова в і -ю последовательность как 0!.
Итерационный шаг: берём одну оследовательность L Обычно, их выбирают по очереди, хотя возможны и другие варианты, например, случайный выбор. Существенно, что у всех последовательностей шансы быть выбранными равны.
Вычисляем ВПМ из выбранных слов всех последовательностей, кроме /. Обозначим эту матрицу как Р. Возьмем каждое слово длины / из /-ой последовательности и вычислим вес этого слова, относительно матрицы Р. Вес вычисляется как отношение вероятностей данного слова быть случайно порождённым и позиционной модели ВПМ и из фоновой модели. - Разыграем следующее Or случайно из всех слов в і длины / используя распределение вероятностей, определенное весами (вероятность - зто вес, нормированный на 1) - Заменяем Оі на Or-
Повторять итерационный шаг, пока ситуация не станет стабильной.
Основная идея алгоритма CONSENSUS заключается в том, что в матрицу выравнивания (описанную выше) по одному добавляют слова длины / из последовательностей, еще не включенных в матрицу. После каждого добавления выбирается набор слов с наибольшим информационным содержанием. Алгоритм работает до тех пор, пока в матрице выравнивания не будет содержаться по одному слову из каждой последовательности или то количество, которое определяет пользователь. В качестве начального слова берутся по очереди все слова длины / из всех последовательностей или только из части входных данных, определенной пользователем. На рисунке 2 приведен пример поиска сигнала длины 4 в трех последовательностях длины 5.
Существует группа алгоритмов, решающих задачу поиска так называемого (/, с/)-сигнала, т.е. сигнала длины /, который входит в каждую последовательность не более чем с rf заменами. Сложность нахождения такого сигнала заключается в том, что два разных вхождения одного сигнала могут различаться в 2d позициях.
Алгоритм WINNOWER строит многодольный граф, вершинами которого являются слова длины /, долями - слова из одной последовательности, а ребрами связаны те вершины, которые отличаются друг от друга не более чем на 2d замен. Многодольный граф - это граф, в котором ребрами соединяются только вершины из разных долей. Сигналу соответствует клика в таком графе. Клика - это такой подграф, в котором любые две вершины соединены ребром. При поиске клики алгоритм сокращает число ребер в графе путем удаления тех, которые не могут быть частью никакой клики заданного размера.
Алгоритм SP-STAR отличается от WINNOWER тем, что строит взвешенный многодольный граф, в котором каждому ребру приписывается вес того, насколько отличаются слова соединенные этим ребром. Теперь задача заключается в поиске клики минимального веса, которая будет соответствовать искомому сигналу.
Алгоритм MITRA (Mismatch Tree Algorithm) использует префиксное дерево для разбиения всех возможных сигналов на непересекающиеся подмножества с общим префиксом с точностью до d замен. Это разбиение делит проблему поиска на более мелкие задачи, в которой можно применить уже известные алгоритмы, например WINNOWER.
Ключевая идея алгоритма PROJECTION заключается в том, чтобы разделить множество всех слов искомой длины из всех последовательностей на части, чтобы потом выбрать только те части, в которых могли бы содержаться вхождения сигнала, и в них применять уже известные алгоритмы, например МЕМЕ или Gibbs sample. Разделение происходит с помощью функции f(x), которая каждому слову х длины / сопоставляет слово длины к, полученное из х соединением оснований, находящихся в этих к фиксированных позициях. Таким образом, мы разделяем все слова длины / на части так, что в каждой слова совпадают в к одних и тех же позициях.
Идеальность модели. Универсальность. Верифицируемость
Строя модель интересующего нас предмета, мы можем ждать или не ждать от неё соответствия имени моделируемого предмета. Это, казалось бы, малосущественное, различие в подходах (в конце концов, от модели ждут адекватности поставленной задаче, а не семантической идеальности) существенным образом определяет парадигму моделирования. Мы будем говорить о более или менее идеальных моделях, имея в виду степень соответствия модели имени моделируемого объекта.
Построение идеальных в этом смысле моделей традиционно для естественных наук, близких к физике. Понятийный ряд этих наук постепенно складывался в процессе когнитивного диалога с миром, модели и понятия развивались вместе, между ними не проводилось различий. Такая модель выражает суть понятия как обобщения наблюдений. Ребёнок в какой-то момент понимает, что все эти все эти шумные маленькие - это погремушки, а наука понимает - что все эти качающиеся и длинные - это математические маятники. Такая модель универсальна, она адекватна всем задачам, в которых участвует моделируемый предмет. Вспомним, что предмет выделяется из мира понятийным фильтром. В нашем случае модель этому фильтру соответствует. Предмет как бы замещён одноименной с ним моделью. Кроме того, результат такого моделирования легко верифицируем: идеальность модели влечёт за собой адекватность результата её применения.
Другая, вычислительная, парадигма моделирования предполагает непосредственно строить адекватное численное описание свойств объекта в рамках поставленной задачи. При этом решённый или не решённый до этого вопрос об идеальности получающейся модели отходит на второй план. Например, моделирование простого числа решетом Эратосфена вполне идеально. Моделирование процесса диффузии параболическим уравнением адекватно большинству приложений и при этом существенно неидеально, в частности, не учитывает конечной скорости распространения фронта диффузии. Неидеальность лежит в основе модели, не учитывающей молекулярной природы вещества. Численное решение уравнения конечно-разностной сеткой делает скорость фронта конечной, модель становится ближе к идеальной: один гранулярный процесс моделируется другим.
Цель построения вычислительной модели - это адекватность описания, и число новых вводимых понятий минимально необходимо (действует принцип бритвы Оккама). Плата за такую понятийную простоту - это верификация применимости модели к конкретной задаче, которая требует привлечения дополнительного экспериментального материала.
Результат вычислительного моделирования — это числовое (пространственно-временное) описание моделируемого события. Аналитическое или алгоритмическое описание - это лишь часть модели. Другая её часть состоит из способов решения или оптимизации уравнений, выбора алгоритмов, их реализации и так далее - то есть представляет собой вычислительные методы. Эта вторая часть, так же, как и первая, влияет и на адекватность модели, и на её идеальность, как в приведённом выше примере про уравнение диффузии. Чем больше вычислительный метод схож с физической или понятийной природой моделируемого объекта, тем более идеальна, и, следовательно, более универсальна модель, включающая этот метод.
Простые модели, адекватные исследуемому событию часто достаточны, но в то же время применимость такой модели часто неочевидна. Другая крайность - идеальная модель - это перевод понятия на язык модели (математический, если это вычислительная модель). Объект появляется в восприятии как проявление этого понятия, и модель естественно его описывает. Но очень многие существенные в науке понятия не могут быть адекватно переведены на математический язык, например, понятие живой клетки. Для такого понятия невозможно привести идеальную модель, что не мешает всем остальным быть более или менее идеальными, в приведённом выше смысле. Разумным компромиссом, по-видимому, является, описание как можно более идеальными моделями максимально сложных компонентов, позволяющих такое описание, например, цитоскелета клетки (Фаворов, Волькенштейн 1991).
Следует заметить, что неидеальная модель необязательно плоха. Качество модели - это её адекватность задаче. Но неидеальная модель не универсальна и требует верификации своего применения.
Существующие модели ДНК-белкового связывания различаются представлением сайта связывания, и способом его поиска в выборке последовательностей. Модель, используемая семейством алгоритмов «Gibbs Sampler» (Lawrence, Altschul et al. 1993; Roth, Hughes et al. 1998; Liu, Bmtlag et al. 2001; Thijs, Marchal et al. 2002; Favorov, Gelfand et al. 2005), близка к идеальной в аспекте "стохастического узнавания".
Белок-регулятор, в соответствие с моделью Берга и вон Хиппеля (von Hippel, Berg 1986; Berg, von Hippel 1988), находится в динамическом равновесии трёх состояний: свободного, неспецифично связанного с нитью ДНК и специфично связанного с нитью ДНК. Переход из первого состояния происходит из-за того, что молекула белка в процессе трёхмерной диффузии встречает нить ДНК, Из второго белок отрывается от нити в свободное состоянии или находит в процессе одномерной диффузии вдоль нити сайт связывания и переходит в третье. Из третьего (специфически связного) молекула может оторваться или сместится по нити, переходя соответственно в первое или второе состояния. Понятно, что для белковой молекулы вероятность быть прикреплённой к своему сайту узнавания гораздо выше, чем к любому другому участку ДНК, или, иными словами, белок специфично связан со своим сайтом узнавания существенно большее время, чем неспецифично - с произвольным участком ДНК. При этом связывание молекул белка и ДНК может происходить на любом участке ДНК, а затем белок диффундирует вдоль нуклеотидной нити.
Используемый в SeSiMCMC при поиске мотива поочерёдный стохастический пересмотр присутствия и позиций сайта в каждой рассматриваемой последовательности по смыслу во многом сходен с описанной моделью трёх состояний, и в этом аспекте модель близка к идеальной.
С другой стороны, модель самого механизма узнавания белком участка ДНК существенно неидельна, поскольку механизм узнавания участка ДНК белком в подробностях неизвестен. Существующие модели исходят из предположения, что участки узнавания одного белка сходны между собой. При этом сами сайты связывания предполагаются реализациями (проекциями на последовательность) некоего мотива (идеала сайта). Модели различаются представлением мотива (примеры разных представлений представлены в обзорах Brejova, DiMarco et al. 2000; Bulyk 2003; Gelfand 2003), но все они используют понятие сайта как реализации мотива, пытаясь определить сайт -участок ДНК, узнаваемый белком. Представляемая работа (см. также Favorov, Gelfand et al. 2005) основана на этой же парадигме. Найденный мотив далеко не всегда несёт в себе искомый сайт узнавания, и результат применения требует верификации. В то же время применённая в данной работе и близкие ей модели, по-видимому, адекватно моделируют ДНК-белковое взамодействие, что подтверждается тем, что с их помощью были успешно идентифицированы сайты узнавания многих белков (Lawrence, Altschul et al. 1993; Герасимова, Гельфанд et al. 2003; Favorov, Gelfand et al. 2005).
Основной шаг поиска сильного мотива
Основная процедура поиска набора сходных сайтов следующая. Работа алгоритма начинается со случайно расположенных сайтов определенной длины, по одному на последовательность. Затем организуется цикл поочередных уточнений позиций сайтов. На каждом шаге выбирается одна (текущая) последовательность. Для однообразия мы рассматриваем отсутствие сайта последовательности как особую позицию сайта (нулевую). На каждом шаге мы подсчитываем нуклеотиды во всех позициях внутри сайта и для фона по всем последовательностям, кроме текущей. Мы оцениваем позиционные вероятности нуклеотидов внутри мотива по формулам (2.1), (2.3), (2.4) и фоновые вероятности по формуле (2.2).
Таким образом, объединяя априорные вероятности и правдоподобие обычным байесовским способом, мы получаем апостериорное распределение позиции сайта в текущей последовательности и разыгрываем новую позицию сайта (возможно, нулевую позицию) из этого распределения. Процесс последовательных итераций продолжается до тех пор, пока цепь последовательных множеств позиций сайтов не сойдется (то есть, изменения от шага к шагу не станут малыми). Приведенный алгоритм аналогичен описанному у Лоуренса в его работе 1993 года (Lawrence, Altschul et al. 1993), со следующими отличиями: мы рассматриваем возможность отсутствия сайта в последовательности стандартным байесовским способом при каждом уточнении позиции сайтов.
Поиск длины и наилучшего сдвига множества найденных сайтов. Фактически алгоритм оптимизирует взаимосогласованность множества позиций сайтов, что делает его очень чувствительным к изменению взаимного расположения сайтов, но слабо чувствительным к одновременным сдвигам всех сайтов как одного целого. Чтобы решить эту проблему, мы время от времени уточняем результат работы алгоритма после схождения основной цепи итераций и затем перезапускаем цепь. Каждое уточнение - это детерминированный поиск наилучшего решения среди всех возможных сдвигов всех сайтов как одного целого. На этом этапе наилучшее множество сайтов определяется по самому высокому информационному содержанию на позицию сайта (ИСП) в множестве вхождений сайтов. Это информационное содержание представляет собой сумму двух компонентов: структурного и пространственного. Оба они определяются как кульбаковские энтропийные расстояния. Структурный компонент - это расстояние между вероятностными моделями нуклеотидов внутри мотива (ППМ) и фоновым вероятностным распределением нуклеотидов.
Теперь счетчики c(i,r) и параметры модели q{i,r) и f{r) оцениваются по всем последовательностям. Формула (2.10) отличается от стандартной кульбаковской энтропии тем, что мы используем с (г, г) как множитель, a q(i,r) — как аргумент логарифма. Расстояния между оцененным распределением вероятностей символов в выравнивании q(i,r) (которые содержат псевдокаунты) и фоновым распределением f(r) рассчитывается исходя из наблюдаемых данных c(itr). В стандартной же кульбаковской энтропийной мере присутствуют два распределения и нет наблюдаемых данных. В этой ситуации q(i,r) участвовал бы в формуле и как множитель, и как аргумент логарифма. Заметим, что константа М в знаменателе уравнения (2.12), приведенного ниже, дает правильную нормализацию величины энтропийного расстояния.
Таким образом, на каждом уточняющем этапе мы находим мотив, для которого достигается максимальное значение ИСП на всей предыдущей цепи розыгрышей позиций сайтов, и затем варьируем длину мотива и абсолютное положение полного набора сайтов как целого, с тем, чтобы оптимизировать значение ИСП в соответствии с уравнением (2.12). Для разделенных мотивов на этом же этапе оценивается длина разделителя; для основания разделителя принимается фоновая вероятностная модель. Для каждой процедуры уточнения, то есть кооперативного сдвига позиций сайтов, как оптимальное принимается то минимальное значение длины разделителя, которое дает локальный максимум ИСП (2.12).
Такая процедура уточнения аналогична приведенной Лоуренсом (Lawrence, Altschul et al. 1993), но отличается от нее следующим образом. Вычисление информационного содержания по уравнению (2.12) содержит уточненную формулу для пространственного компонента. Каждый этап уточнения используется также для оценки оптимальной длины мотива и длины промежутка, если разрешены разделенные мотивы. Для каждой длины мотива длина промежутка принимается как минимальное значение, для которого достигается локальный максимум ИСП. Поскольку длина промежутка может оказаться нулевой, эта же процедура определяет, разделен мотив или нет,
Работа программы (рис. 2) состоит из двух стадий: отжига и поиска глобального максимума. На этапе отжига программа идет от случайного начального состояния к состоянию, находящемуся не очень далеко от оптимального. Уточнения на этом этапе не изменяют длины мотива. Мы считаем, что все последовательности содержат по одному сайту (вероятность отсутствия сайта в последовательности временно обнуляется). Все это необходимо для того, чтобы найти не очень плохую выборку сайтов, которая позволила бы продолжать более тонкую настройку на следующем этапе. Отжиг заканчивается, когда изменения набора сайтов становятся медленнее определенной скорости. Для этого предусмотрено два параметра: число последовательных полных циклов изменений позиций сайтов, приведших к малому изменению набора сайтов, и критический размер этого малого изменения.
Блок-схема работы алгоритма SeSiMCMC. Геометрический критерий близости наборов сайтов строится следующим образом. Сайты считаются совпавшими, если они перекрываются наполовину. Если они перекрываются меньше, чем на половину, такому наложению приписывается некий вес, меньший едишшицы. Если не перекрываются вообще - наложение нулевое. Затем мы можем сложить все попарные наложения сайтов разных наборов во всех последовательновтях, и поделить на число самих последовательностей. Этот аналог скалярного произведения и является мерой сходства наборов. Информационный критерий сходства - это информационное расстояние (условная информация) от первого набора до второго, отнесённое к информационному содержания первого набора.
Программная реализация алгоритма и важнейшие параметры поиска мотива
Программное обеспечение SeSiMCMC написано на языке C++ (gcc 3.x). Исполняемые файлы для FreeBSD и консоли Win32 доступны на сайте проекта http://bioinform.genctica.ru/SeSiMCMC. На данной странице также находится веб-интерфейс к этому программному обеспечению. Ниже мы дадим его описание с точки зрения пользователя. Здесь же остановимся подробнее на работе самой программы. Все параметры, кроме обязательной для ввода выборки последовательностей ДНК, изначально установлены в значении по умолчанию. Рассмотрим подробнее некоторые из них.
Несколько параметров позволяют пользователю выбирать пространственную геометрию мотива и устанавливать априорно ожидаемый диапазон длин и разумную стартовую длину мотива. Ожидаемая доля последовательностей, которые не содержат вхождения мотива, также может указывается параметром который называется называется motive absence prior. Дополнительно этот параметр выражает предпочтение пользователя о желаемом мотиве: предпочитается ли часто встречающийся слабый мотив или сильный, но редкий. Чем ниже значение параметра, тем большее внимание уделяется частоте встречаемости мотива. На стадии вычислений в алгоритме предполагается, что каждая последовательность содержит не больше одного вхождения мотива. При выдаче результатов программа сканирует все последовательности полученной на стадии вычислений позиционной вероятностной матрицей q{itr)jf{r) (2.1), (2.2) и извлекает все сайты, вес соответствия которых этой матрице выше, чем этот вес для наихудшего из сайтов, полученных в результате вычислений. Это послевычислительное сканирование опционально, и если программа не получила указаний проводить сканирование, то она выводит только наилучшее локальное выравнивание. Если послевычислительное сканирование извлекает слишком много сайтов, это значит, что наихудший из полученных в вычислении сайтов плохо соответствует мотиву. В этом случае разумно повторить вычисление, повысив априорную вероятность отсутствия сайта в последовательности. Также возможно указать программе установить порог извлекаемых на этапе сканирования сайтов равным весу соответствия і-того, а не первого наихудшего сайта. Как и другие инструменты (Thijs, Marchal et al. 2002), SeSiMCMC позволяет искать множественные мотивы, перезапуская вычисления на тех же входных данных с замаскированными сайтами, найденными ранее.
Программа читает выборку фрагментов последовательностей ДНК в формате FastA из стандартного ввода, выводит протокол работы в поток вывода стандартной ошибки и пишет результаты работы в поток стандартного вывода. Мы используем один символ, дополнительный по отношению к FastA стандарту. Это буква X (или х), которая означает "эта позиция не может находиться внутри мотива". Этот символ необходим, чтобы запрещать повторную идентификацию ранее найденного мотива при поиске последующих. Программа выводит протокол работы в стандартный поток ошибки, если этот вывод не подавляется ключами запуска программы. Результат пишется на стандартный выход в формате простого текста или html в зависимости от —html ключа запуска программы.
При выходе программа отмечает хорошие (с высоким информационным содержанием) столбцы мотива заглавными буквами. Хороший столбец имеет информационное содержание большее, чем любая трехбуквенная смесь из фонового распределения.
Все ключи командной строки могут быть опущены при запуске программы, при этом значения по умолчанию работают в большинстве случаев, и лишь иногда имеет смысл менять стандартное поведение программы. Если один и тот же параметров установлен ключом командной строки и меткой в конфигурационном файле, переключатель командной строки имеет приоритет. Если два ключа противоречат друг другу, выводится ошибка в файл диагностики и программа останавливает работу. Если две или более строки конфигурационного файла противоречат друг другу, программа работает в соответствии с последней.
Использует тестовый набор данных. Мы не читаем FastA-файл со стандартного ввода; вместо этого используем сгенерированную программой последовательность из 20 фрагментов ДНК длиной 120 со вставленным словом "atggccactt" в позициях: 84 81 78 75 absent 69 66 63 60(compl) 50 47 44 41 38 35 32 29 26 23 20. Absent означает, что пятая позиция не содержит сайта, а слово "compl" показывает, что восьмая последовательность содержит комплементарный сайт в 60-й позиции. Те же соглашения поддерживаются в тестовых выдачах. —config-file filename -f file Имя конфигурационного файла. Формат каждой строки конфигурационного файла следующий: метка = значение. Строка, начинающаяся с #, считается комментарием. Метки конфигурационного файла описываются ниже. Каждой метке соответствует ключ командной строки, который имеет приоритет перед меткой конфигурационного файла, если указаны оба. Последовательность меток в файле не имеет значения, если они не противоречат друг другу.
Сравнение программы SeSiMCMC с аналогичными
Программное обеспечение SeSiMCMC участвовало в сравнении с двенадцатью другими программами для предсказания мотивов, также доступными через Интернет (Tompa, Li et аІ. 2005) (см. http ://bio.cs.washtn gton.edu/assessnient/). Сравнение производилось на специально подготовленных выборках данных. Сайты связывания ТФ, их позиции и ориентации в последовательностях брались из базы данных TRANSFAC (http://transfac.gbf-braunschweig.de/TRANSFAC/) (Wingender, Dietze et al. 1996). Каждый ТФ был основой для одной выборки входных данных. Каждая выборка организовывалась из одного из трёх типов фоновых последовательностей, в который сайты из TRANSFAC встраивались в их позиции и в правильной ориентации. Эти три типа - это реальные промоторы, в которых расположены сайты, случайные промоторы из того же генома, и, наконец, случайные последовательности, порождённые Марковской моделью.
Сравнение коэффициентов корреляции идентифицированных и исходных наборов сайтов для SeSiMCMC и 12 других программ. Fly - сайты мухи; human - человека; mouse - мыши; yeast - дрожжей; all data - общие результаты. ГЛАВА 4. Применение алгоритма для поиска участков связывания белков - регуляторов дыхания АгсА и NarP на ДНК в бактерии Е. coli и описания соответствующих регулонов
Для проверки работы алгоритма, мы исследовали два транскрипционных фактора, сайты связывания которых общеизвестны как дивергентные и поэтому трудные для вычислительной идентификации, АгсА-Р и NarP, оба включённые в регуляцию дыхания. Многое указывает на то, что оба эти сайта имеют симметричную структуру. Эксперименты показали {Darwin, Tyson et al. 1997), что сайт NarP, вероятно, имеет палиндромную структуру. Сайт АгсА был описан как повторяющийся в нескольких копиях на одной нити ДНК. (Darwin, Tyson et al. 1997; McGuire, De Wulf et al. 1999; Liu, De Wulf 2004), что является сильным доводом в пользу предположения о структуре прямого повтора. Обе эти регуляторные системы жизненно важны для бактерий и хорошо экспериментально изучены (Darwin, Tyson et al. 1997; McGuire, De Wulf et al. 1999; Liu, De Wulf 2004).
В обоих случаях, мы сочетали вычислительную процедуру идентификации мотива с методами сравнительной геномики, анализируя сайты, найденные в регуляторних (upstream-ных) областях ортологичных генов из нескольких родственных геномов. Эта процедура в деталях описана в работе (Gelfand, Koonin et al. 2000). Параметры поиска мотива в обеих выборках, сами выборки и результаты поиска доступны как примеры запуска программы на сайте проекта SeSiMCMC, http;//bioinform. penetica.ru/SeSiMCMC. Анализ методами сравнительной геномики выполнялся с помощью программного продукта GenomeExplorer (Mironov, Vinokurova et al. 2000). Переключатель кислородного/бескислородного дыхания Phospho-ArcA.
Система АгсА регулирует экспрессию ряда генов в ответ на изменения аэробных/анаэробных условий окружающей среды (Iuchi, Weiner 1996). Эта система двухкомпонентна, она включает в себя мембранно-связанную сенсорную киназу АгсВ и белок-регулятор транскрипции (ТФ) АгсА. В анаэробных условиях АгсВ фосфорилирует сам себя и затем катализирует фосфорилирирование АгсА, стимулируя этим связывание АгсА с ДНК. В аэробных условиях трансфофорилирования АгсА не происходит. Фосфорилированный АгсА (Phospho-ArcA, ArcA-P), свою очередь, связывается с ДНК и подавляет транскрипцию нескольких оперонов, в том числе . icd, lid, git, glc, sdh и sodA, и активирует транскрипцию нескольких других, например, e.g. cyd and pjl. (Lynch, Lin 1996). Сейчас известно около дюжины АгсА - регулируемых оперонов, но в то же время есть данные, что что в этот тип регуляции как напрямую, так опосредованно вовлечены сотни генов (Liu, De Wulf 2004), что делает идентификацию сайта связывания крайне важной для понимания бактериального метаболизма.
В качестве входных данных были взят набор регуляторных участков генов, для которых различными методами было показано, что они регулируются системой АгсА (см. таблицу 1). SeSiMCMC запускалась для этой обучающей выборки последовательностей с различными параметрами для поиска (1) произвольного мотива, (2) мотива с возможным симметричным промежутком (спенсером) посередине, (3) прямого повтора, возможно, с симметричным спейсером и (4) палиндрома (встречного комплементарного повтора), возможно, с симметричным спейсером. Сайты искались на обеих нитях ДНК. Длина сайта могла быть между 6 и 22 основаниями. 2-1 1 0-5 CMrtTfirxoNcooiOT-cgco io
Наилучший найденный мотив имел структуру прямого повтора со спейсером (см. рис. 7, а). Этот мотив длиной в 15 оснований более консервативен и содержит больше информативных позиций, чем мотив АгсА, приведённый в работе (McGuire, De Wulf et al. 1999) (см. рис. 7, b). Мы считаем, что это уточнение произошло из-за предопределённой симметрии мотива при поиске. Когда идёт поиск мотива произвольной формы, более сильное плечи повторов выравниваются друг по другу, независимо от того, справа или слева они встречаются в повторе, и в результате информация из более слабого плеча игнорируется, поскольку оно выравнивается с неинформативными фланками. В случае поиска повтора, слабое плечо лучше отличается от неинформативной части, что позволяет более точно определить ядро мотива.
Высокая селективность полученного мотива позволила использовать его в исследовании регулона АгсА-Р в геномах четырёх гамма-протеобактерий: Escherichia coli, Yersiniapestis, Pasteurela multocida и Vibrio vulnificus методами сравнительной геномики. На основе всех сайтов, определённых программой как образующих мотив, был построен (Миронов, Гельфанд 1999) распознающий профиль (описанная выше формулой (0.1) разновидность ПВМ). Затем, регуляторные участки геномов этих четырёх бактерий сканировались этим профилем и в случае, если регулярная область гена в Exoli содержала участок, соответствующий профилю с вемом выше порога 4.25, а хотя бы два гена-ортолога из трех остальных геномов содержали в апстриме участок с ценой соответствия профилю выше порога 4.00, ген считался несущим консервативный бокс АгсА. Это распознающее правило позволило выявить ряд новых потенциально ArcA-регулируемых генов в гамма-протеобактериях, большинство их которых упоминаются в литературе как значимые для дыхательно-зависимой регуляции (Герасимова, Гельфанд et al. 2003) {см. таблицу 2).
