Инструментальное средство поиска регуляторных мотивов в геномах

Инструментальное средство поиска регуляторных мотивов в геномах

В процессе жизнедеятельности клетки не все гены экспрессируются одновременно. Это достигается за счет регуляции. Понимание механизма регуляции экспрессии генов — важнейшая задача биологии. При изучении регуляции экспрессии на уровне транскрипции важно не только определить белки-регуляторы (транскрипционные факторы), но и участки их связывания с последовательностью ДНК. В настоящее время в открытом доступе находится большое количество секвенированных геномов и данных по экспрессии генов, что позволяет изучать регуляцию путем анализа последовательностей с помощью вычислительных методов. Задача поиска регуляторных мотивов в наборе последовательностей ДНК — классическая задача биоинформатики. К настоящему моменту создано огромное количество алгоритмов поиска мотивов, однако все они имеют свои ограничения, и не существует универсального алгоритма, который решает эту задачу.

Известно, что алгоритмы, комбинирующие различные методы, наиболее эффективны и универсальны. В данной работе мы представляем алгоритм поиска мотивов в последовательностях ДНК, совмещающий словарные техники и методики, использующие скрытые марковские модели (СММ).

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

Цель работы

Дополнить разработанный в лаборатории биоинформатики ФББ алгоритм поиска мотивов в нуклеотидных последовательностях возможностью идентификации мотивов сложной структуры с вариабельным спейсером; реализовать эту модификацию на языке программирования высокого уровня; разработать графический интерфейс и веб-ресурс для обеспечения открытого доступа к алгоритму и удобного просмотра результатов.

Задачи

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

a. Палиндром

b. Повтор

c. Инвертированный повтор.

2. Построить схемы скрытых марковских моделей (СММ) для поиска мотивов с более сложной структурой.

3. Реализовать алгоритм с помощью языка программирования Java.

4. Протестировать алгоритм на искусственных и реальных данных и сравнить с другими алгоритмами поиска мотивов.

Разработать графический интерфейс и веб-ресурс для обеспечения открытого доступа к алгоритму и удобного просмотра результатов.

1. Обзор литературы

1.1 Транскрипционные факторы

1.1.1 Общие сведения

Инициация транскрипции — сложный процесс, эффективность которого зависит от того, как устроена последовательность ДНК непосредственно вблизи начала транскрибируемой области (а у эукариот также и в более далеких участках генома — энхансерах и сайленсерах), а также от наличия или отсутствия различных белковых транскрипционных факторов.

Факторы транскрипции — белки, которые регулируют транскрипцию путем связывания со специфичными участками ДНК — сайтами связывания. Транскрипционные факторы выполняют свою функцию самостоятельно либо в комплексе с другими белками. Различают репрессорные и активирующие транскрипционные факторы, которые соответственно снижают или повышают константу связывания РНК-полимеразы с регуляторными последовательностями экспрессируемого гена.

Определяющая черта факторов транскрипции — наличие в их составе одного или более ДНК-связывающих доменов, которые взаимодействуют с характерными участками ДНК, расположенными в регуляторных областях генов.

Транскрипционные факторы бывают конститутивные (всегда активные в клетке) и активируемые (активны только в определенных условиях). Активируемые в свою очередь разделяют на тканеспецифические (участвуют в развитии организма) и сигнал-зависимые, или рецепторы (требуют внешнего сигнала для активации).

Для функционирования транскрипционных факторов чаще всего необходимо формирование гетеродимерного или гомодимерного комплекса. Например, гетеродимерные комплексы различных ядерных рецепторов с ретиноидным Х рецептором (RXR). Существуют также и гомодимерные комплексы RXR (рис. 1).

Рис. 1. Пример гетеродимера RXR/TR и гомодимера RXR

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

В зависимости от того, как части димера расположены друг относительно друга, сайт связывания такого димера с ДНК может представлять собой палиндром, прямой повтор или инвертированный повтор (Таблица 1) [5, 6].

Таблица 1. Типы структур сайтов связывания транскрипционных факторов

алгоритм поиск мотив марковский модель

Палиндром — это сайт, полностью соответствующий своему обратному комплементу (например, CACGTG). Сайты прямых и инвертированных повторов состоят из двух плеч, разделенных промежутком — спейсером. Длина спейсера часто постоянна, но иногда может варьироваться на несколько нуклеотидов. Плечи прямых повторов имеют совпадающую последовательность. Плечи инвертированных повторов обратно комплементарны. Если длина спейсера инвертированного повтора равна нулю, этот сайт можно назвать палиндромом.

1.1.2 Различные структуры сайтов связывания

Рассмотрим различные структуры сайтов связывания на примере рецепторов стероидных гормонов. Это внутриклеточные рецепторы, чаще всего локализованные в ядре и осуществляющие передачу сигнала от стероидных гормонов.

ДНК-связывающий домен стероидных рецепторов содержит аминокислоты, специфично связывающиеся с гормон-чувствительнным элементом (сайтом связывания) на последовательности ДНК. Этот участок состоит из 66−68 высоко консервативных аминокислот, из которых 8 цистеинов образуют 2 структуры типа цинковых пальцев, которые взаимодействуют с ДНК. Остальные аминокислоты ДНК-связывающего домена определяют специфичность связывания различных стероидных рецепторов (рис. 2).

Рис. 2. Структура ДНК-связывающего домена ядерных рецепторов. Синими кружками обозначены остатки цистеина, образующие координационные связи с цинком (Zn), оранжевыми — аминокислотные остатки, непосредственно контактирующие с нуклеотидами, зелеными — аминокислотные остатки, участвующие в димеризации рецепторов

Общая схема взаимодействия такова: два рецептора связываются с гормоном, а затем образуют гомодимер. Этот гомодимер связывается с гормон-чувствительным элементом. Затем в процесс транскрипции включаются другие транскрипционные факторы и РНК полимераза II, что стабилизирует преинициативный комплекс и запускает синтез мРНК (рис. 3).

Рис. 3. Общая схема механизма работы рецепторов стероидных гормонов. HRE — гормон-чувствительный элемент, pol II — РНК полимераза II

Сайт связывания чаще всего расположен в промоторной области или на расстоянии нескольких килобаз до TATA и CAAT боксов, которые находятся рядом с сайтом начала транскрипции. Предполагают, что в последнем случае позиционирование нуклеосомы может усиливать стимулирующее действие рецепторов на транскрипцию за счет образования петли (рис. 4).

Рис. 4. Участие нуклеосомы в образовании петли для усиления действия рецептора на процесс транскрипции. NR — ядерный рецептор, HRE — гормон-чувствительный элемент, TF — транскрипционный фактор, TFBS — его сайт связывания, Pol — РНК полимераза II

Гомодимеры рецепторов I типа связываются с сайтами, имеющими структуру типа палиндром или инвертированный повтор со спейсером длиной в 3 нуклеотида. Гомодимеры рецепторов II типа связываются с сайтами, имеющими структуру типа прямой повтор с вариабельным спейсером длины 0−5 нуклеотидов (рис. 5).

Рис. 5. Взаимодействие связанных с гормоном (черные треугольники) гомодимеров рецепторов стероидных гормонов с гормон-чувствительным элементом (HRE)

Размер спейсера между полусайтами гормон-чувствительных элементов определяет взаимодействие с ДНК димерных ядерных рецепторов. Чем больше длина спейсера, тем более специфичен гормон-чувствительный элемент (рис. 6).

Рис. 6. Зависимость специфичности сайтов связывания транскрипционных факторов от длины спейсера на примере различных гетеродимеров RXR

1.2 Способы представления регуляторных элементов

Наиболее распространенными способами представления последовательностей сайтов связывания белков с ДНК являются консенсус (регулярное выражение) и позиционная весовая матрица (PWM — position weight matrix, или PSSM — position-specific scoring matrix). Консенсус представляет собой общий вид последовательности сайта — слово, составленное из нуклеотидов, наиболее часто встречающихся в конкретных позициях сайта. Часто для учета вариаций в некоторых позициях консенсуса помимо основных четырех букв используют обозначения вырожденных нуклеотидов в соответствии с нормами IUPAC. Консенсусы хорошо подходят для описания сайтов связывания белков, которые связываются со строго консервативной последовательностью (например, белки системы рестрикции-модификации II-ого типа).

Однако консенсус не позволяет хорошо описать сайты в том случае, если последовательность сайта сильно варьируется. PWM, которые впервые были введены для характеристики сайтов инициации транскрипции и трансляции у E.coli [9, 10], значительно лучше подходят для описания сайтов связывания факторов транскрипции, так как способны количественно охарактеризовать частые и редкие вариации в последовательности сайтов, что невозможно в случае регулярных выражений.

PWM представляют собой матрицу L Ч 4 (L — длина сайта), каждый элемент которой отражает частоту встречаемости данного нуклеотида в данной позиции сайта. Вес, порождаемый матрицей при выравнивании с данным участком последовательности, обычно вычисляется как сумма элементов матрицы, соответствующих нуклеотидам, стоящим в каждой позиции рассматриваемого участка (рис. 7).

Рис. 7. Конструирование позиционной весовой матрицы (PWM) сайта связывания фактора транскрипции (TFBS). (а) Выравнивание десяти известных последовательностей TFBS. (b) Подсчет частот появления каждого нуклеотида в каждой позиции сайта (в данном случае величины не нормированы). Эта таблица обычно и называется позиционной весовой матрицей. (с) Для визуализации PWM часто используется диаграмма logo, на которой степень консервативности позиции показана высотой букв

Таким образом, PWM предоставляет достаточно полное описание участка ДНК, с которым способен связываться конкретный белок, и может быть применена при сканировании геномной последовательности для поиска сайтов, дающих достаточно хороший вес. Использование PWM позволяет достаточно эффективно предсказывать сайты связывания белков. Так, например, для 95% сайтов связывания тканеспецифического фактора печени HNF-1, найденных в последовательностях приматов из GenBank с использованием соответствующей PWM и отличающихся наиболее высоким весом, было экспериментально показано связывание с HNF-1 in vitro.

На настоящий момент существует две наиболее полные курируемые базы PWM сайтов связывания факторов транскрипции: TRANSFAC и JASPAR. JASPAR содержит значительно меньше данных, при том что каждому транскрипционному фактору соответствует только одна PWM, тогда как TRANSFAC содержит по несколько PWM для некоторых факторов. Кроме этого существует несколько баз данных, содержащих регуляторные области генов (SCPD [15], TRRD [16]), а также недавно созданная база данных UniPROBE [17], которая содержит сайты связывания транскрипционных факторов, полученные с помощью технологии белок-связывающих микрочипов (protein binding microarray, PBM).

Поиск сайтов связывания белков in silico — это только первый шаг к определению действительно функциональных сайтов. Регуляция генов сильно зависит также от структуры хроматина и ДНК-метилирования [19−21]. Большая часть хромосомной ДНК представляет собой компактно упакованный гетерохроматин и вследствие этого изолирована от взаимодействия с транскрипционными факторами. Метилирование ДНК тоже может препятствовать связыванию факторов с определенными участками ДНК, а также влиять на структуру хроматина. Поэтому многие потенциальные сайты, обнаруживаемые при полногеномном поиске без учета этих факторов, не являются функциональными in vivo, хотя они были бы способны связывать определенные транскрипционные факторы, будучи открытыми для взаимодействия.

Следует отметить, что, несмотря на все свои достоинства, PWM все-таки имеет несколько недостатков. Одним из них является то, что стандартная PWM не учитывает взаимное влияние соседних позиций сайта (мононуклеотидная модель). Однако наличие таких зависимостей было показано для некоторых факторов [22−24]. В таких случаях модели более высокого порядка (то есть учитывающие зависимость позиций сайта), например, динуклеотидные PWM, демонстрируют более аккуратное предсказание потенциальных сайтов [24−27].

К тому же, в некоторых случаях только половина (или даже меньше) позиций матрицы обладают достаточно высоким уровнем консервативности, в результате чего эффективность поиска с помощью такой матрицы падает. Иногда такая консервативность PWM отражает специфичность транскрипционного фактора, который сам по себе слабо взаимодействует с ДНК, а точность действия достигается только в контексте соседних сайтов связывания. Тем не менее, в большинстве случаев низкое качество PWM объясняется не свойствами транскрипционного фактора, а, скорее, недостаточно корректным составлением PWM.

Не очень высокое качество PWM может также объясняться малым числом сайтов, известных для данного фактора. В этом случае PWM может не отражать всех возможных вариаций сайтов, вследствие чего при поиске такой матрицей большое количество реальных сайтов не может быть найдено. В таких случаях имеет смысл объединять сайты связывания факторов одного семейства, так как последние часто имеют очень похожую структуру и способ связывания. Кроме того, при недостатке известных сайтов в некоторых случаях можно создавать модели регуляторных элементов, используя информацию о структуре ДНК-белковых взаимодействий. Методы, опирающиеся на информацию такого рода, пока не многочисленны (в основном из-за малого числа расшифрованных структур ДНК-белковых комплексов), однако в последнее время эта область активно развивается [29, 30]. Такие методы позволяют не только предсказать новые регуляторные мотивы [31], но и улучшить качество уже имеющихся PWM.

В случае, когда недостаток сайтов восполнить не удается или сайты слишком консервативны, при построении PWM используют искусственный прием «размывания» матрицы. Для этого часто используются псевдоотсчеты. Простейший вариант псевдоотсчетов — прибавить до нормировки к каждому счетчику нуклеотидов в позиции PWM какую-то величину. Величина псевдоотсчетов обычно выбирается так, чтобы их сумма была пропорциональна, где N — количество последовательностей в выравнивании.

Для оценки качества PWM часто используется энтропийное расстояние (или условное информационное содержание) от фонового распределения частот по формуле Кульбака-Лейбера [34]:

где I — информационное содержание, f (b, j) — наблюдаемая частота нуклеотида b в позиции j, p(b) — фоновая частота нуклеотида b.

1.3 Алгоритмы поиска мотивов

Исторически сложилось, что большинство существующих алгоритмов создано для поиска мотивов в регуляторных областях перед совместно регулируемыми генами из одного генома. Одновременное изменение экспрессии генов чаще всего вызвано совместной транскрипционной регуляцией. Таким образом, задачу поиска сайтов связывания транскрипционного фактора можно свести к задаче поиска мотива в наборе последовательностей ДНК. В случае прокариот связывание достигается в большей степени за счет аффинности сайта связывания и транскрипционного фактора, сайты связывания довольно длинные, и, как правило, перед геном присутствует один сайт. Поэтому для поиска таких мотивов чаще используются методики, способные искать один достаточно консервативный сайт, представленный в каждой последовательности. В случае эукариот сайты короткие и вырожденные, и связывание достигается в большей степени за счет большого количества сайтов в последовательности, нежели чем за счет аффинности. Поэтому, в случае эукариот поиск сильно осложняется: искомый мотив определяется как набор не очень консервативных сайтов, которые перепредставлены в исходных последовательностях.

Впоследствии стало известно, что у высших эукариот регуляторные сайты могут образовывать так называемые композиционные элементы (composite elements, CEs) [35], то есть небольшие группы сайтов, характеризующиеся определенным взаиморасположением. Биологические причины, ведущие к такому неслучайному расположению сайтов, понятны: транскрипционные факторы, связываясь с ДНК, также взаимодействуют между собой для достижения нужного влияния на уровень транскрипции [36, 37]. Другими словами, расположение регуляторных сайтов обусловлено трехмерной структурой белкового комплекса, вовлеченного в инициацию транскрипции. В самом простом случае СЕ — это пара сайтов связывания определенных транскрипционных факторов, совместно влияющих на экспрессию гена.

Массовое секвенирование геномов позволило использовать близкородственные геномы для анализа регуляции. Были разработаны алгоритмы, берущие на вход только промоторные участки ортологичных генов и использующие методы межвидового геномного сравнения, или филогенетического футпринтинга. Основная идея этого метода состоит в том, что функциональные элементы в последовательностях ДНК находятся под давлением отбора. Поэтому консервативные сайты в наборе регуляторных областей ортологичных генов скорее всего являются функциональными регуляторными элементами (рис. 8). Для определения таких элементов чаще всего строится множественное выравнивание промоторных областей ортологичных генов, а затем на нем выделяются консервативные участки.

Рис. 8. Применение сравнительной геномики к поиску регуляторных модулей. (а) выравнивание последовательностей далеких видов обнаруживает высоко консервативные некодирующие участки. Диаграммы демонстрируют высокую степень консервативности между последовательностями некодирующих областей перед геном Pax6 из геномов человека, мыши, крысы и рыбы Fugu. (b) Консервативность этого участка выше, чем ожидалось [39]

Позднее были созданы алгоритмы, комбинирующие два основных подхода для поиска мотивов в последовательностях ДНК, применяя их одновременно или по очереди. Авторы утверждают, что такие алгоритмы крайне эффективны, но их использование не всегда возможно из-за отсутствия данных.

Итак, с точки зрения исходных данных, алгоритмы поиска мотивов можно разделить на три основные группы:

1. Алгоритмы, использующие различные наборы последовательностей для поиска в них мотивов.

2. Алгоритмы, использующие методы сравнительной геномики для поиска мотивов в промоторных областях ортологичных генов из разных видов.

3. Алгоритмы, комбинирующие два подхода.

1.3.1 Алгоритмы поиска мотивов в наборе последовательностей

С алгоритмической точки зрения методы поиска мотивов в наборе последовательностей делят на:

1. переборные алгоритмы, основанные на словарных техниках

2. алгоритмы, использующие различные вероятностные модели.

Переборные алгоритмы обеспечивают нахождение глобального оптимального решения, но при этом на больших выборках работают довольно долгое время. К переборным алгоритмам относятся: Oligo-Analysis [40, 41], YMF [42−44], алгоритмы, использующие суффиксные деревья [45−48], и методы на основе графов [49, 50]

Алгоритмы, применяющие вероятностные модели, хороши тем, что находят приблизительное решение за реальное время. Это позволяет применять их к большим выборкам. Недостатком является то, что такие алгоритмы используют несколько параметров для поиска, которые необходимо подбирать. К сожалению, все вероятностные алгоритмы не гарантируют нахождения лучшего решения, так как используют различные формы локального поиска. К ним относятся: Consensus [51, 52], NestedMICA [53], алгоритмы, использующие метод максимизации ожидания (expectation maximization, EM) [54, 55], алгоритм Gibbs sampling [56, 57] и дополнения к нем.

Переборные алгоритмы, основанные на словарных техниках

Ван Хельден и др. разработали алгоритм поиска мотивов, названный Oligo-Analysis. Данный алгоритм ищет в последовательностях короткие перепредставленные слова — участки, частота встречаемости которых в начальных последовательностях выше соответствующих фоновых частот. Фоновые частоты были рассчитаны для каждого слова из всех последовательностей некодирующих участков геномов дрожжей. Несмотря на общую простоту, алгоритм показал высокую эффективность при поиске мотивов в регуляторных последовательностях дрожжей (Saccharomyces cerevisiae). К сожалению, данный алгоритм можно использовать только для поиска довольно коротких консервативных мотивов. Позднее, ван Хельден и др. усовершенствовали свой алгоритм, добавив в него возможность искать мотивы, состоящие из двух частей, разделенных спейсером. Так как спейсер может быть различным для одного мотива, длину промежутка можно варьировать от 0 до 16. Частота такого двойного мотива может быть вычислена как сумма частот двух плеч или же как общая частота двойного мотива. Основным недостатком алгоритма ван Хельдена является то, что в нем ищутся точные вхождения слов, то есть не учитывается вариабельность сайтов.

Томпа обратил внимание на эту проблему, и представил свой алгоритм, использующий словарную технику, для поиска коротких мотивов в последовательностях ДНК. В процессе его работы для каждого отрезка s длины k рассчитывается значение N_{s — количество вхождений слова s в исходные последовательности с допустимым количеством замен. Также рассчитывается значение N's, вычисленное для случайно сгенерированной последовательности той же длины. Мерой того, является ли s мотивом, считается разность Ns — N's.}

В дальнейших работах этот подход был усовершенствован. Пусть Х — отдельная случайная последовательность длины L. Фоновая частота каждого нуклеотида полагается равной 0.25, или же вычисляться по начальному набору данных. Предположим, что p_{s — это вероятность того, что Х содержит хотя бы одно слово s длины k или же любого его соседа (то есть слово, отличающееся в нескольких позициях). Если предположить, что в наборе из N случайных последовательностей длины L последовательности независимы, предполагаемое количество встреч слова s и его соседей в этом наборе есть, стандартное отклонение равно .}

Тогда

где — z-score или отклонение в стандартных единицах. Величина имеет стандартное нормальное распределение и позволяет сравнивать различные мотивы. Томпа предложил эффективный алгоритм оценки, использующий марковские модели.

Используя подобный подход, Синха и Томпа [43, 44] разработали алгоритм YMF (Yeast Motif Finder), в котором для расчета фонового распределения частот последовательности генерируются с помощью марковской модели. Для определения параметров модели используются все существующие последовательности ДНК дрожжей. Алгоритм возвращает мотивы с наибольшей величиной z-score. Авторы протестировали свой алгоритм на выборках из геномов дрожжей и показали его высокую эффективность.

Ванет и др. использовали суффиксные деревья для представления набора последовательностей при создании алгоритма для поиска единичных мотивов в полных геномах бактерий. Марсан и Сагот добавили в этот алгоритм поиск комбинаций мотивов. Представление набора последовательностей в виде суффиксного дерева давало огромное количество возможных решений, но, несмотря на это, методика оказалась эффективной.

Существуют и другие алгоритмы, использующие суффиксные деревья и их вариации, такие как Weeder и MITRA (Mismatch Tree Algorithm), созданные Павеси и др. и Эскиным и Певзнером соответственно, а также алгоритмы, использующие словарные техники совместно с графовыми методиками, такие как WINNOWER и cWINNOWER.

Вероятностные алгоритмы

Одним из первых вероятностных методов поиска сайтов связывания транскрипционных факторов стал вероятностный алгоритм Хертца и др. Он является жадным и ищет мотив, представленный в виде PWM, с наибольшим информационным содержанием. Предполагается, что каждая исходная последовательность содержит ровно один сайт. Позднее этот алгоритм был усовершенствован. В его последней версии (Consensus), разработанной Хертцем и Стормо [52], используется следующий метод. Строится PWM по одному случайному слову длины l. Далее по очереди из каждой последовательности выбирается слово, имеющее максимальный вес по PWM и добавляется к исходному слову. После каждого добавления выбирается набор слов с наибольшим информационным содержанием. По полученным словам PWM перестраивается.

Большинство вероятностных алгоритмов поиска мотивов используют эвристические методы, такие как метод максимизации ожидания и Gibbs sampling, а также дополнения к ним.

Метод максимизации ожидания

Одним из широко известных методов оценки параметров вероятностных моделей, позволяющих эффективно работать с большими объемами данных, является EM-алгоритм. Его название происходит от слов «expectation-maximization», что переводится как «ожидание-максимизация». Это связано с тем, что каждая итерация содержит два шага: вычисление математических ожиданий (expectation) и максимизацию (maximisation). Алгоритм основан на методике итеративного вычисления оценок максимального правдоподобия, предложенной в 1977 г.

EM-алгоритм впервые был применен для поиска мотивов Лоренцем и Рейлли. Их алгоритм — это дополнение к жадному алгоритму Хертца и др. Первоначально этот алгоритм был разработан для поиска белковых мотивов, но он также может использоваться и для поиска в последовательностях ДНК. Метод не требует никакого выравнивания сайтов в последовательностях, но изначально предполагает, что каждая из них включает только один сайт. Набор сайтов находится методом, описанным выше (см. обзор литературы, вероятностные модели). Неточность в расположении сайтов устраняется с помощью метода максимизации ожидания, работающего следующим образом. Пусть g_,jk — вероятность того, что искомый мотив начинается в j позиции в последовательности k, а f (i, ?) — вероятность символа ? в колонке i для каждого символа из алфавита и 1? i ?l. В процессе работы алгоритма происходит переопределение g и f до тех пор, пока f не будет мало изменяться. Это переопределение происходит с использованием байесовских методик.

Алгоритм MEME [55], созданный Бейли и Элканом, применяет стратегию EM для поиска мотивов. В алгоритме MEME используются три новаторские идеи для поиска мотивов. Во-первых, участки начальных последовательностей используются как отправные точки для алгоритма. Это позволяет повысить скорость работы алгоритма и вероятность нахождения глобально оптимального решения. Во-вторых, отменено требование о том, что в каждой последовательности должен встречаться в точности один сайт. В-третьих, за счет особенностей вероятностной модели, появилась возможность одновременно находить в одном наборе данных сразу несколько мотивов.

Метод Gibbs sampling

Gibbs sampling представляет собой довольно широко используемый стохастический аналог метода максимизации ожидания. Этот алгоритм ищет максимальное значение функции от нескольких переменных. Основная идея алгоритма заключается в следующем. На каждом шаге случайным образом выбирается одна переменная, и ее значение меняется при фиксированных значениях других переменных. Если это изменение приводит к возрастанию целевой функции, переменная принимает выбранное значение. Процесс повторяется до тех пор, пока значение функции не перестанет значимо меняться.

Метод Gibbs sampling был разработан Геманом и Геманом для восстановления изображений. Впервые в биоинформатике этот метод был применён для построения множественного выравнивания Лоренсом и др. Различные модификации этого метода часто применяются для поиска слабых мотивов. Gibbs sampling — это подход, использующий марковские цепи и метод Монте-Карло (MCMC). Для расчета вероятностей на данном шаге марковские цепи используют вероятности, полученные на предыдущем шаге. Суть метода Монте-Карло состоит в том, что выбор следующего шага осуществляется случайно с вероятностью, зависящей от текущего состояния. В самой простой версии метода Gibbs sampling мы ищем лучший консервативный неразрывный мотив длины l в виде PWM. Предполагается, что искомый сайт встречается во всех последовательностях.

Поиск осуществляется в несколько итераций. Сначала случайно выбирается по одному слову длины l из каждой последовательности. Эти слова формируют начальное множество вхождений мотива. Обозначим позицию слова в i-й последовательности через O_i.

Итерационный шаг:

Берём одну последовательность i. Обычно последовательности выбирают по очереди, хотя возможны и другие варианты, например, случайный выбор. Существенно, что у всех последовательностей шансы быть выбранными равны.

Строим PWM по выбранным словам из всех последовательностей, кроме i. Берем каждое слово длины l из i-ой последовательности и вычисляем Байесовскую вероятность того, что данное слово могло бы быть порождено PWM, а не фоном.

Разыграем следующее O_{i' случайно из всех слов в последовательности i длины l, используя полученное распределение вероятностей.}

Заменяем O_{i на Oi'.}

Итерационный шаг повторяется до тех пор, пока набор слов не станет неизменным.

Дополнения к методу Gibbs sampling

Рот и др. создали на основе метода Gibbs sampling алгоритм поиска мотивов AlignACE (Aligns Nucleic Acid Conserved Elements). Основные отличия от оригинального алгоритма Gibbs sampling состоят в следующем. Во-первых, фоновые частоты нуклеотидов фиксированы и соответствуют частотам в геноме (например, 62% A+T в случае дрожжей). Во-вторых, на каждом шаге алгоритм ищет мотив по двум цепям одновременно, и перекрывающиеся сайты исключены, даже если они находятся на разных цепях. В-третьих, в AlignACE используется метод MAP (maximum a posteriori log-likelihood) для оценки качества полученных в результате мотивов. Это мера того, что мотив встретился в последовательности не случайно. Важной особенностью метода MAP является то, что в нем учитывается не только распределение частот нуклеотидов в рассматриваемых геномах, но и некоторые другие особенности (например, А-богатые участки в ДНК дрожжей). В результате, мотивы, порожденные такими особенностями генома, исключаются из конечных результатов.

Позднее, Хьюгес и др. использовали AlignACE для поиска мотивов в наборе функционально важных генов дрожжей. Вместо MAP для оценки найденных мотивов в данном случае использовался усовершенствованный метод. Учитывались особенности конкретных начальных данных и выделялись те мотивы, которые более вероятно являются реальными сайтами, чем случайными мусорными последовательностями.

Фиджс и др. разработали алгоритм MotifSampler, использующий алгоритм Gibbs sampling со следующими изменениями. Во-первых, наличие только одного сайта в последовательности более не является обязательным. Во-вторых, в алгоритме используются марковские цепи высокого порядка для построения фоновой модели. Алгоритм применяли для поиска регуляторных мотивов в геномах бактерий и некоторых растений.

На основе алгоритма Gibbs sampling, Лью и др. разработали алгоритм BioProspector, обрабатывающий промоторные области перед совместно регулируемыми генами. Основные отличия от оригинального Gibbs sampling состоят в следующем. Во-первых, в нем используются марковские цепи от нулевого до третьего порядка для построения фонового распределения. Параметры для них задаются пользователем или вычисляются по исходным последовательностями. Во-вторых, алгоритм позволяет искать двойные мотивы, разделенные спейсером, и палиндромные мотивы. Алгоритм использовали для поиска сайтов связывания как в прокариотах, так и эукариотах (дрожжах).

Шида разработал алгоритм поиска мотивов GibbsSt, использующий метод имитации теплового отжига (simulated annealing [63]) совместно с алгоритмом Gibbs sampling. Позже стало известно, что этот метод гораздо лучше решает проблемы, связанные с нахождением локально лучшего решения. В биоинформатке метод имитации теплового отжига в основном применяется для улучшения методов поиска в пространстве решений [65, 66]. В алгоритме GibbsST метод имитации теплового отжига используется для улучшения работы алгоритма Gibbs sampling.

Другие подходы

Хью и др. использовали комбинированный подход для создания алгоритма поиска мотивов EMD. Алгоритм основан на кластеризации. В нем используется комбинация предсказаний, полученных из множества пробегов одного или более различных алгоритмов поиска: AlignACE, Bioprospector, MDScan [71], MEME и MotifSampler. Алгоритм в 22.4% случаев показал более высокий результат, нежели все компоненты алгоритма отдельно. EMD показал наибольшую эффективность в случае поиска в коротких последовательностях. В случае поиска в длинных последовательностях, он всегда более или по крайней мере также эффективен, как отдельные элементы алгоритма.

Каплан и др. создали алгоритм, использующий помимо последовательностей ДНК информацию о структуре ДНК-связывающих доменов известных транскрипционных факторов. По ним предсказываются возможные сайты связывания, которые ищутся в последовательностях.

Лью и др. разработали алгоритм, основанный на нейронных сетях, для поиска мотивов в последовательностях ДНК и белковых последовательностях. Сеть содержит несколько уровней. Предсказание мотивов происходит поступательно: на верхнем уровне последовательность разбивается на небольшие участки, а на нижнем эти участки классифицируются на мотивные и фоновые. При этом полученные данные сохраняются и используются для уточнения результатов в следующих тестах. Основное преимущество такого алгоритма в том, что он хорошо работает с длинными последовательностями

Кингсфорд и др. разработали алгоритм для поиска мотивов в последовательностях ДНК, который ищет набор подпоследовательностей определенного размера таким образом, чтобы сумма попарных расстояний между ними была минимальна. Для этого используется целочисленное линейное программирование (ILP). Преимуществом данного алгоритма является то, что он работает относительно небольшое время на выборках любой величины. Тестирование на последовательностях из E.coli показало эффективность алгоритма, сопоставимую с эффективностью некоторых методов, основанных на алгоритме Gibbs sampling.

Ле и др. создали генетический алгоритм HIGEDA, использующий в начальной стадии алгоритм EM для поиска лучших параметров модели мотива. Помимо этого, HIGEDA может искать мотивы не только с мутациями, но и с инсерциями и делециями.

1.3.2 Алгоритмы, основанные на методе филогенетического футпринтинга

Основное преимущество филогенетического футпринтинга по сравнению с подходом, использующим совместно регулируемые гены, состоит в том, что определить ортологичные гены часто бывает проще, чем совместно регулируемые. На сегодняшний день в открытом доступе находится большое количество аннотированных геномов, в том числе близкородственных, что позволяет применять технику филогенетического футпринтинга для поиска мотивов. Для определения регуляторных элементов в последовательностях чаще всего строится множественное выравнивание промоторных областей ортологичных генов, а затем на нем выделяются особо консервативные участки. Построение множественного выравнивания осуществляется при помощи таких алгоритмов, как CLUSTAL W.

К сожалению, было показано [76−78], что алгоритмы, использующие филогенетический футпринтинг не всегда применимы. Если сравниваемые виды слишком близки друг другу в смысле эволюционного расстояния (например, различные штаммы одного вида) выравнивание последовательностей очевидно, но не информативно, поскольку функциональные элементы не более консервативны, чем окружающая нефункциональная последовательность. Если же последовательности очень сильно разошлись, сложно построить удовлетворительное выравнивание. В этом случае совместно с филогенетическим футпринтингом часто используются такие существующие алгоритмы поиска мотивов, как MEME, Consensus или Gibbs sampling.

Клифтен и др. использовали AlignACE для поиска мотивов в сравнительном анализе последовательностей ДНК нескольких видов Saccharomyces, и получили хорошие результаты в тех случаях, когда построить глобальное выравнивание было невозможно. Маккью и др. использовали алгоритм Gibbs sampling совместно с филогенетическим футпринтингом для поиска мотивов в геномах протеобактерий.

Бланшетт и Томпа создали эффективно работающий алгоритм поиска мотивов, использующий филогенетический футпринтинг совместно с динамическим программированием. В своей работе Бланшетт и Томпа отметили, что алгоритмы поиска мотивов часто не учитывают степень эволюционной близости последовательностей и считают их независимыми. Это особенно критично в случае, если производится анализ большого количества геномов, среди которых присутствуют как очень близкие в эволюционном отношении (штаммы), так и далекие (типы и выше). В этом случае построить выравнивание трудно. К тому же, если в наборе исходных данных есть несколько групп разного размера, состоящих из очень близких видов, то с точки зрения алгоритма поиска большую значимость будут иметь мотивы, присутствующие в самой многочисленной группе. Даже использования методов взвешивания последовательностей не достаточно, чтобы решить такую проблему.

Клифтен и др. использовали филогенетический футпринтинг для поиска мотивов в шести геномах Saccharomyces. Авторы применили СLUSTAL W для выравнивания последовательностей. Они обнаружили множество статистически достоверных консервативных участков. Но этот результат был получен только потому, что геномы для исследования были тщательно отобраны так, чтобы расстояние между ними было оптимальным.

1.3.3 Алгоритмы, комбинирующие различные подходы

Комбинированные алгоритмы могут обрабатывать данные, состоящие из последовательностей и перед совместно регулируемыми генами, и перед ортологичными генами. При этом они используют различные методы поиска мотивов, в том числе и филогенетический футпринтинг. Такие алгоритмы учитывают два важнейших аспекта определения значимости мотивов: перепредставленность и межвидовую консервативность.

Алгоритм Келлиса и др. осуществляет поиск в две стадии: сначала в смешанных данных находятся высоко консервативные участки, а уже среди них ищутся перепредставленные. Пракаш и др. разработали алгоритм OrthoMeme, использующий метод максимизации ожидания, в котором два вида поиска в смешанных данных осуществляется одновременно.

Ванг и Стормо разработали алгоритм PhyloCon, основанный на алгоритме Consensus. В нем поиск осуществляется в две стадии. Сначала для областей перед ортологичными генами строится множественные выравнивания и выделяются консервативные области, по которым составляются PWM. Далее с помощью полученных PWM мотивы ищутся перед всеми генами каждого генома. Авторы показали, что PhyloCon имеет очень низкий уровень перепредсказаний. Алгоритм хорошо работает как на коротких, так и на длинных последовательностях.

Особенностью алгоритма, разработанного Синха и др. является то, что при построении множественного выравнивания допускается условие, что мотив может встретиться не во всех последовательностях. Авторы протестировали алгоритм на последовательностях из геномов дрожжей, мухи и даже человека. Сравнение с такими алгоритмами, как MEME, OrthoMEME, PhyloGibbs [85], EMnEm и GIBBS (Wadsworth Gibbs sampler) показало, что алгоритм более эффективен в большинстве случаев.

1.3.4 Сравнительный анализ алгоритмов поиска мотивов

Сейчас доступно огромное количество алгоритмов поиска мотивов. К сожалению, невозможно выделить один универсальный алгоритм: каждый из них имеет свои ограничения. В литературе описано несколько работ по сравнению алгоритмов.

Томпа и др. сравнили 13 алгоритмов поиска мотивов: AlignACE, ANN-Spec [89], Consensus, GLAM [90], Improbizer [91], MEME, MITRA, MotifSampler, Oligo/Dyad-Analysis, QuickScore [92], SeSiMCMC [67], Weeder и YMF. Для этого были созданы наборы данных, состоящие из последовательностей, содержащих известные сайты связывания транскрипционных факторов из базы данных TRANSFAC. Авторам алгоритмов было предложено протестировать свои инструменты на данных наборах с четко установленными параметрами и предоставить для сравнения лучший полученный результат. Сравнение результатов показало общую низкую эффективность работы всех алгоритмов. Однако алгоритм Weeder превзошел другие в большинстве случаев. При этом SeSiMCMC оказался эффективнее при тестировании на последовательностях из геномов мух, а MEME3 (разновидность MEME) и YMF — на последовательностях из геномов мышей. Авторы предположили, что алгоритмы, использующие различные подходы для поиска мотивов, будут более эффективны и универсальны.

Хью и др. также сравнили различные алгоритмы, но при этом использовали последовательности из базы данных регулонов E.coli RegulonDB [93], а также позволили авторам самим подбирать оптимальные параметры поиска и предоставлять не только лучший, но и другие полученные результаты. В эксперименте участвовали пять алгоритмов: AlignACE, MEME, BioProspector, MDScan и MotifSampler. К сожалению, это исследование также показало общую низкую эффективность всех алгоритмов. Основной проблемой практически всех алгоритмов была длина последовательностей — при увеличении значения этого параметра, результаты резко ухудшались. Лидером был признан популярный алгоритм MEME, показавший эффективность в 52% относительно среднего уровня 15%-35%.

1.4 Скрытые марковские модели и вспомогательные алгоритмы

Скрытая марковская модель (СММ, Hidden Markov Model, HMM) — статистическая модель, имитирующая работу процесса, похожего на марковский процесс с неизвестными параметрами. Задача состоит в оценке неизвестных параметров по наблюдаемым данным. Полученные параметры могут быть использованы в дальнейшем анализе. СММ может быть рассмотрена как сеть условных Байесовских вероятностей.

Первые заметки о скрытых марковских моделях опубликовал Баум в 1960;х, и уже в 70-х их впервые применили при распознавании речи. С середины 1980;х СММ применяются при анализе биологических последовательностей, в частности, ДНК.

Элементами скрытой марковской модели являются состояния. Обозначим состояние в момент времени t через x(t) (рис. 9.). Каждое состояние имеет эмиссионные вероятности — распределение среди всех возможных выходных значений. Наблюдаемое значение в момент времени t обозначим через y(t). Между состояниями определены переходные вероятности. Для СММ 1-го порядка состояние x(t) в момент времени t зависит только от состояния x (t? 1) в предыдущий момент времени. Это называется свойством Маркова.

Рис. 9. Общая структура СММ

Вероятность увидеть последовательность значений длины L равна

Здесь сумма пробегает по всем возможным последовательностям скрытых узлов. Метод подсчёта значений P(Y) полным перебором — очень трудоёмкий для многих реальных задач, где количество возможных состояний очень велико. Но существуют эвристические алгоритмы, позволяющие решить эту задачу за приемлемое время. К таким алгоритмам относятся Витерби и вперед-назад (forward-backward).

Алгоритм вперед-назад вычисляет P(Y) для конкретной последовательности наблюдаемых значений и восстанавливает последовательность состояний по последовательности наблюдаемых значений (разметка). СММ можно представить в виде многодольного графа, где узлы соответствуют состояниям, а доли графа — позициям в последовательности наблюдаемых значений. Тогда разметка последовательности наблюдаемых значений в терминах СММ есть путь в данном графе. Основная идея алгоритма заключается в следующем. Вероятность того, что путь пройдет через данный узел графа, складывается из всех путей, которые ведут в данную точку из начальной и из конечной. Алгоритм включает в себя две стадии: проход «вперед» и проход «назад». При проходе «вперед» для каждой точки вычисляется вероятность прийти в эту точку из начальной, а при проходе «назад» — из конечной. Далее для каждой позиции можно сравнить вероятности всех состояний и выбрать наиболее вероятное.

В данной работе для обучения СММ использовался алгоритм Витерби, который находит наиболее вероятную последовательность состояний (путь) для конкретной последовательности наблюдаемых значений. Алгоритм рекурсивен, необходимо наличие начального и конечного состояний.

Пусть дана СММ с набором состояний X, эмиссионными вероятностями e_i, в состоянии x_i и переходными вероятностями a_i,j из i-го состояния в j-ое,. Требуется найти наиболее вероятную последовательность состояний

Пусть вероятности наиболее вероятного пути в состояние k с наблюдаемым значением i известны для всех k. Тогда

Полный алгоритм:

Инициализация (i = 0):, для k > 0.

Рекурсия (i = 1…L):

Терминация:

Обратный проход (i = L…1):

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

2. Материалы и методы

2.1 Базовый алгоритм

2.1.1 Общая схема

В качестве базового алгоритма был использован алгоритм поиска непалиндромных мотивов определенной длины, разработанный в лаборатории биоинформатики ФББ. Алгоритм берет на вход последовательности ДНК, в которых необходимо найти мотив (набор сайтов). Поиск мотива осуществляется в два основных этапа:

1. Поиск перепредставленных слов-кандидатов определенной длины в последовательностях

2. Обучение СММ, созданных на основе слов-кандидатов.

2.1.2 Поиск слов-кандидатов

В качестве исходных параметров на этом этапе используется длина слова-затравки и количество допустимых замен.

Сначала создается словарь из всех возможных слов заданной длины. Соседями в таком словаре считаются слова, количество замен между которыми (расстояние по Хэммингу) не превышает допустимый уровень. Затем для слов из словаря считается количество вхождений в последовательности с заданным количеством замен. Если в последовательности найдены слова из словаря, то счетчик вхождений увеличивается не только для данного слова, но и для всех его соседей. Эта операция проводится для реальных последовательностей и для случайно сгенерированных последовательностей такой же длины с теми же частотами нуклеотидов. Таким образом, мы получаем реальную и фоновую частоты вхождений для каждого слова из словаря с заданным уровнем замен. Для того чтобы получить достоверные фоновые частоты, производится несколько (около 10) генераций случайных последовательностей.

Далее по полученным данным для реального и фонового распределений вычисляется величина hCount (см. ниже), а затем строится частотная гистограмма для этой величины (рис. 10).

где count — количество вхождений слова, min — минимальное значение гистограммы, step — шаг гистограммы, scale — шкала гистограммы, dictsize — размер начального словаря.

Рис. 10. Гистограмма фонового и реального распределений величины hCount. Здесь n(hCount) — количество слов из словаря с данным значением hCount, foreground — реальное распределение, background — фоновое распределение, foreground (trend) и background (trend) — линии тренда для соответствующих распределений