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

 

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

Введение

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

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

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

Цель работы

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

Задачи

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

a.Палиндром

b.Повтор

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

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

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

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

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

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

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

.1.1 Общие сведения

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

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

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

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

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

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

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

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

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

Палиндром

- GACTGCGCAGTC-3

3 - GACTGCGCAGTC-5Прямой повтор

- GACTGCagtGACTGC-3

3 - GCAGTCactGCAGTC-5Инвертированный повтор

- GACTGCagtGCAGTC-3

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

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

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

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

ДНК-связывающий домен стероидных рецепторов содержит аминокислоты, специфично связывающиеся с гормон-чувствительнным элементом (сайтом связывания) на последовательности ДНК. Этот участок состоит из 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)

Размер спейсера между полусайтами гормон-чувствительных элементов определяет взаимодействие с ДНК димерных ядерных рецепторов. Чем больше длина спейсера, тем более специфичен гормон-чувствительный элемент [8] (рис. 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 [11] с использованием соответствующей PWM и отличающихся наиболее высоким весом, было экспериментально показано связывание с HNF-1 in vitro [12].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Переборные алгоритмы обеспечивают нахождение глобального оптимального решения, но при этом на больших выборках работают довольно долгое время. К переборным алгоритмам относятся: 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] и дополнения к нем.

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

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

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

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

Тогда

,

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

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

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

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

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

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

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

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

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

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

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

Метод Gibbs sampling

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Первые заметки о скрытых марковских моделях опубликовал Баум в 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, эмиссионными вероятностями ei, в состоянии xi и переходными вероятностями ai,j из i-го состояния в j-ое, . Требуется найти наиболее вероятную последовательность состояний

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

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

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

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

Терминация:

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

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

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

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

.1.1 Общая схема

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

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

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

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

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

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

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

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

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

Нетрудно заметить, что в правой части гистограммы (в области около 60-70) значения n(hCount) реального распределения в какой-то момент начинают превышать соответствующие значения фонового. В этой области находятся часто встречающиеся слова (то есть слова с большой величиной count и, соответственно, hCount), которых в реальных последовательностях больше, чем в сгенерированных, то есть эти слова перепредставлены. Пики гистограммы в области глобального максимума соответствуют словам с самопересекающейся структурой (например, ATGCGATG)

Далее строится отношение гистограмм реального и фонового распределений (рис. 11). Пороговое значение (xCount) для отбора перепредставленных слов-кандидатов выставляется в наименьшей точке, где отношение реального распределения к фоновому превышает 1.

В данном подходе при построении фонового распределения учитываются только частоты нуклеотидов в исходных данных, но не другие особенности, которые должны моделироваться СММ более высокого порядка. Предположим, существует слово, вероятность встретить которое в сгенерированных данных очень мала. В результате, для того, чтобы это слово оказалось перепредставленным, ему и его соседям достаточно встретиться в исходных данных всего несколько раз. Такие редко встречающиеся слова (то есть слова с маленькой величиной count и, соответственно, hCount) выделяются в левой части реального распределения (область около 45-50), где отношение распределений также превышает 1. Но они не являются искомыми мотивами. Поэтому было введено условие, что пороговое значение xCount должно всегда быть правее максимума фонового распределения.

Рис. 11. Гистограмма отношения реального распределения частот вхождений слов к фоновому. Пороговое значение выбрано равным 61

Далее составляется список реальных слов-кандидатов (т.е. участков начальных последовательностей), соответствующих словам из словаря, величина count которых больше порогового значения xCount. Если два перепредставленных слова перекрываются на реальной последовательности, то в список слов-кандидатов добавляется и их объединение (большей длины). Эта операция позволяет строить словарь из коротких слов (например, 8 нуклеотидов), что существенно экономит память и ускоряет работу алгоритма, а затем объединять их до нужной длины. Кроме того, такой подход в какой-то мере позволяет решить проблему поиска мотива неизвестной длины.

Каждое слово-кандидат используется для определения эмиссионных и переходных вероятностей СММ.

2.1.3 Структура СММ

Структура СММ базового алгоритма содержит следующие состояния:

·Begin - «молчащее» (т.е. не имеющее эмиссионных вероятностей) состояние, определяющее начало работы.

·Background (bckg) - состояние фона. Начальные эмиссионные вероятности для всех нуклеотидов равны 0.25.

·Site - «молчащее» состояние, определяющее начало сайта на последовательности.

·Pos 1, pos 2 … pos n - состояния, определяющие позиции сайта.

·End - «молчащее» состояние, определяющее конец работы.

Рис. 12. Схема СММ базового алгоритма

На рисунке 12 изображена схема СММ для базового алгоритма, стрелками показаны возможные переходы между состояниями СММ. Переходные вероятности между состояниями pos 1, pos 2…pos n (позиции в сайте) равны 1 и не меняются в процессе обучения. Начальные эмиссионные вероятности состояний pos 1 … pos n задаются в соответствии с буквой в последовательности слова-кандидата и «размываются» псевдокаунтами. Такая СММ позволяет найти в последовательности произвольное число сайтов.

.1.4 Обучение СММ и применение для поиска сайтов

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

.2 Модификации алгоритма

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

.2.1 Начальная обработка последовательностей

Известно, что в последовательностях ДНК существуют так называемые области низкой сложности [94] - чаще всего это высоко периодичные AT-богатые участки. Наличие таких участков в начальных последовательностях приводит к появлению множества бессмысленных перепредставленных слов-кандидатов и завышению порогового значения. В нашей программе мы применили процедуру маскировки областей низкой сложности до начала основной работы алгоритма. После составления словаря и получения слов-кандидатов обучение СММ проводится на полных исходных последовательностях.

.2.2 Модификация процедуры поиска слов-кандидатов

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

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

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

В процессе тестирования было замечено, что для маленьких выборок начальных данных (около 5 последовательностей длины 200), пороговое значение иногда выбирается неверно по случайным причинам из-за недостаточного количества информации. Для решения этой проблемы пороговое значение уменьшается на величину , где N - количество последовательностей в начальной выборке, - коэффициент поправки, вычисленный экспериментально - методом Монте-Карло с 10000 итераций. В случае больших выборок сдвиг на такую величину не оказывает большого влияния на результат, но помогает восполнить недостаток информации в случае маленьких выборок (рис. 13).

Рис. 13. Гистограммы реального (foreground) и фонового (background) распределений со сдвинутым пороговым значением (xCount) для выборки из 25 последовательностей длины 200 нуклеотидов

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

В результате слово-кандидат выбирается по следующим критериям:

1.В исходных последовательностях присутствует повтор, соответствующий данному слову

.Слово перепредставлено в исходных последовательностях (частота вхождений соответствующего ему слова из словаря выше порогового значения)

.В составе слова должна быть хотя бы одна буква С или G.

Перекрывающиеся слова-кандидаты не рассматриваются.

.2.3 Кластеризация слов-кандидатов

Так как слова-кандидаты - это реальные участки начальных последовательностей, среди них встречаются идентичные или очень похожие. Чтобы не создавать несколько одинаковых СММ, слова-кандидаты кластеризуются с помощью алгоритма, использующего бинарное дерево. В одном кластере могут находиться слова-кандидаты, отличающиеся в одной или двух конкретных позициях. Кластеры представлены в виде PWM, составленной по входящим в него последовательностям. Рассматриваются только кластеры, содержащие более двух слов.

.2.4 Модификация СММ

Помимо уже известных состояний, модифицированная СММ содержит также состояния:

·Spacer - «молчащее» состояние, определяющее начало промежутка между плечами структуры

·Spacer 1, spacer 2 … spacer n - состояния, определяющие позиции спейсера

·Pos 1, pos 2 … pos n - состояния, определяющие позиции второго плеча мотива

Рис. 14. Схема модифицированной СММ

На рисунке 14 показаны возможные переходы между состояниями модифицированной СММ. Такая схема СММ подходит для любой структуры мотива. Переходные вероятности между состояниями site, pos 1 … pos n и pos 1 … pos n равны 1 и не меняются в процессе обучения. Эмиссионные вероятности состояний pos 1 … pos n задаются в соответствии с PWM, составленной по кластеру слов-кандидатов. В случае повторов эмиссионные вероятности состояний pos 1 … pos n равны вероятностям состояний pos 1 … pos n соответственно, в случае палиндромов и инвертированных повторов - вероятностям симметричных (комплементарных) состояний (pos n … pos 1). Из молчащего состояния spacer возможен переход сразу в следующее плечо сайта, что позволяет найти мотивы без спейсера. Из конечного состояния второго плеча сайта возможен обратный переход в спейсер, что позволяет найти мотив, состоящий больше, чем из 2-х плеч. Такое иногда встречается в случае повторов.

2.2.5 Модификация процесса обучения

В процессе обучения СММ, эмиссионные вероятности в состояниях pos 1 … pos n и pos 1 … pos n обучаются одновременно: после каждой итерации обучения эмиссии соответствующих позиций двух плеч усредняются.

.2.6 Уточнение мотива

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

.3 Веб-ресурс

Для обеспечения открытого доступа к алгоритму был разработан веб-интерфейс #"justify">Страница формы (рис. 15) содержит следующие поля: файл с начальными данными, структура искомого сайта (прямой повтор, инвертированный повтор или палиндром), минимальная и максимальная длина спейсера, длина искомого сайта, количество допустимых замен между плечами сайта, порог на вес итоговых сайтов, количество генераций случайных последовательностей для расчета фонового распределения, требование о том, что каждый сайт должен иметь хотя бы одну букву C или G, и возможность маскировать области низкой сложности.

Рис. 15. Страница формы

На странице результатов (рис. 16) отображена необходимая информация о мотиве: информационное содержание соответствующей PWM, logo, полученное на ее основе, и список найденных сайтов.

Рис. 16. Страница результатов

3. Результаты

.1 Тестирование на сгенерированных данных

Первый тест состоял из выборок, содержащих 5, 15, 25 или 50 последовательностей, в каждой из которых содержался ровно один искомый сайт. Выборки варьировались по количеству допустимых замен между плечами сайта (от 0 до 3). Для каждого случая было сгенерировано 10 тестовых файлов. Тестирование было проведено для сайтов со структурами типа прямой и инвертированный повтор с вариабельным спейсером.

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

Рис. 17. Результаты тестирования алгоритма на выборках разной. Здесь по оси абсцисс отложено количество последовательностей, по оси ординат - средняя по 10 экспериментам доля верно найденных сайтов. Mism - количество допустимых замен между плечами искомых сайтов в начальных данных

Также в процессе данного тестирования было оценено время работы алгоритма (Таблица 2). Технические характеристики компьютера, на котором производилось тестирование: процессор Intel® Core™ Duo 2.0 ГГц, 3Гб оперативной памяти.

Таблица 2. Время работы алгоритма на выборках разного объема

Величина выборкиВремя работы алгоритма5 последовательностей10-15 сек.15 последовательностей<1 мин.25 последовательностей1-2 мин.50 последовательностей5-6 мин.100 последовательностейОколо 40 мин.200 последовательностейБолее 2 ч.

Во втором тесте оценивалась чувствительность алгоритма относительно количества сайтов в каждой из исходных последовательностей. Начальные выборки состояли из 15 последовательностей длиной 200 нуклеотидов и варьировались по количеству сайтов в каждой последовательности и уровню допустимых замен между плечами сайта. Первая серия выборок содержала по два искомых сайта в каждой последовательности, вторая - по одному, а каждая следующая - на два сайта меньше, чем предыдущая. Каждая серия состояла из 10 сгенерированных выборок

Результаты тестирования (рис. 18 и рис. 19) показали прямую зависимость средней доли найденных сайтов от количества сайтов в последовательности: чем больше сайтов в среднем приходится на последовательность, тем выше доля найденных сайтов.

Рис. 18. Чувствительность алгоритма в зависимости от процентного содержания сайтов в последовательности при поиске прямых повторов. По оси абсцисс - доля сайтов в последовательности в процентах. По оси ординат - средняя доля верно найденных сайтов. Mism - количество допустимых замен между плечами искомого сайта

Рис. 19. Чувствительность алгоритма в зависимости от процентного содержания сайтов в последовательности при поиске инвертированных повторов. По оси абсцисс - доля сайтов в последовательности в процентах. По оси ординат - средняя доля верно найденных сайтов. Mism - количество допустимых замен между плечами искомого сайта

Тестирование показало высокую эффективность алгоритма. В случае прямых повторов доля найденных сайтов снижается постепенно при уменьшении их количества в последовательности. При этом доля верно найденных сайтов держится выше 0.8 вплоть до 40% сайтов в последовательности. В случае же инвертированных повторов доля найденных сайтов практически равна 1 вплоть до 60% сайтов в последовательности, а затем эффективность резко падает. Такая разница для прямых и инвертированных повторов обусловлена особенностями СММ. Предположим, что в последовательности до или после искомого сайта есть подпоследовательность, похожая на одно из плеч сайта. В таком случае при поиске сайтов типа прямой повтор существует возможность по случайным причинам выбрать сайт, имеющий в качестве одного плеча плечо искомого сайта, а в качестве второго - случайную последовательность. При поиске сайтов типа инвертированный повтор такое невозможно, так как плечи сайта имеют не идентичную, а обратно комплементарную последовательность, и для СММ это будет совершенно другой сайт.

3.2. Тестирование на реальных данных

Алгоритм был применен для поиска сайтов связывания двух транскрипционных факторов в геномах рода Shewanella - MetJ и BirA. Метиониновый репрессор MetJ уменьшает экспрессию метионинового регулона и белков, вовлеченных в синтез S-аденозилметионина. Сайт связывания представляет собой прямой повтор с плечом дины 8 без спейсера. Выборка содержит 62 последовательности длиной около 250 нуклеотидов. Биотиновая голоэнзимная синтетаза, кодируемая геном BirA, отвечает за репрессию биотинового оперона, адсорбцию биотина и удержание его внутри клетки. Сайт связывания представляет собой инвертированный повтор с плечом дины 8 и спейсером длины 15. Выборка содержит 15 последовательности длиной около 200 нуклеотидов. В обоих случаях в качестве лучшего результата был найден верный мотив (Таблица 3).

Таблица 3. Результаты тестирования нашего алгоритма на последовательностях из геномов Shewanella

Транскрипционный факторИнформационное содержаниеLogoMetJ_Shewanella (Repeat 8-0-8)14.519BirA_Shewanella (Inverted repeat 8-15-8)10.455

3.3 Сравнение с другими алгоритмам

Мы сравнили наш алгоритм с другими алгоритмами поиска мотивов в последовательностях ДНК: MEME [55] и SeSiMCMC [67]. Оба алгоритма показали высокую эффективность в проведенных ранее сравнительных экспериментах (см. обзор литературы, сравнительный анализ алгоритмов поиска мотивов), кроме того они общедоступны, удобны для пользования и работают достаточно быстро.

Алгоритм MEME не приспособлен для поиска двойных мотивов, разделенных спейсером, поэтому тестирование проводилось для палиндромных сайтов. Сравнение алгоритмов осуществлялось на двух тестовых выборках: «пуриновой» и «аргининовой».

«Пуриновая» выборка изначально состояла из 19 последовательностей, содержащих 20 экспериментально определенных сайтов связывания фактора транскрипции PurR. Большинство сайтов в выборке одинарные, то есть по одному сайту в каждой последовательности. Сайты при этом довольно консервативные, с четко выраженной палиндромностью.

«Аргининовая» выборка изначально состояла из 8 последовательностей, содержащих 19 сайтов связывания факторов транскрипции ArgR. Сайты в выборке двойные и слабые (два сайта в одном фрагменте ДНК, каждый слабо похож на исходный мотив). Палиндромность мотива также очень слабо выражена.

Далее, один за другим мы вырезали сайты из последовательностей выборки, то есть каждый следующий тест содержал на один сайт меньше, чем предыдущий. Таким образом, все меньшее число последовательностей в тесте содержало искомый мотив. Каждый раз мы убирали самый сильный сайт, то есть сайт, имеющий наибольший вес относительно PWM.

При данном тестировании оценивалась чувствительность и специфичность рассматриваемых алгоритмов. Чувствительность (sensitivity) - это доля верно найденных сайтов из всех реально существующих сайтов в начальной выборке. Специфичность (specificity) - это доля верно найденных сайтов из всех найденных.

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

3.3.1 «Пуриновый» тест

При поиске сайтов в начальном файле с полным набором сайтов все три алгоритма нашли искомый мотив в качестве лучшего.

Чувствительность нашего алгоритма сопоставима с чувствительностью MEME. Алгоритм SeSiMCMC находит верный мотив в выборках с низким содержанием сайтов, но в целом чувствительность этого алгоритма ниже, чем у других.

Специфичность алгоритма MEME выше, чем у других алгоритмов. На выборках с большим количеством исключенных сайтов наш алгоритм показал более высокую специфичность, чем SeSiMCMC

3.3.2 «Аргининовый» тест

При поиске сайтов в начальном файле с полным набором сайтов все три алгоритма нашли искомый мотив в качестве лучшего. MEME нашел сайт со сдвигом на один нуклеотид.

Чувствительность и специфичность нахождения мотива на аргининовом тесте у всех алгоритмов сильно ниже, чем на пуриновом. Это объясняется тем, что мотив в «аргининовой» выборке значительно слабее, чем в «пуриновой». На выборках с большим количеством сайтов наш алгоритм показал чувствительность, сравнимую с MEME, а на выборках с меньшим количеством сайтов превосходит его. SeSiMCMC показал относительно низкую общую чувствительность на данном тесте, однако он смог найти верный мотив в выборках с наименьшим количеством сайтов в исходных последовательностях. По специфичности (рис. 27) наш алгоритм превосходит MEME и SeSiMCMC на выборках с низким содержанием сайтов.

Заключение

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

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

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

Выводы

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

.Алгоритм реализован на языке Java в среде разработки Eclipse.

.Тестирование на сгенерированных и реальных данных показало высокую эффективность алгоритма.

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

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

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

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

Больше работ по теме: