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

Нетрудно заметить, что в правой части гистограммы (в области около 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 задаются в соответствии с буквой в последовательности слова-кандидата и «размываются» псевдокаунтами. Такая СММ позволяет найти в последовательности произвольное число сайтов.

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

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

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

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

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

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

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

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

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

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

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



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

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

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

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

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

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

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

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

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

2.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.2.6 Уточнение мотива

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

2.3 Веб-ресурс

Для обеспечения открытого доступа к алгоритму был разработан веб-интерфейс http://bioinf.fbb.msu.ru/HMMSiteFinder/. Приложение реализовано с помощью технологии Apache Struts (http://struts.apache.org/) на языке Java в среде разработки Eclipse. Эта технология является расширением Java Servlet API и в архитектурном плане реализует идеологию MVC (Model-View-Controller).

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



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

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

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

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

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

Транскрипционный фактор	Информационное содержание	Logo
MetJ_Shewanella (Repeat 8-0-8)	14.519
BirA_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. Разработанный в лаборатории биоинформатики ФББ алгоритм поиска мотивов в нуклеотидных последовательностях дополнен возможностью идентификации мотивов сложной структуры с вариабельным спейсером. Для этого модифицированы все стадии базового алгоритма, в том числе разработаны новые схемы СММ. Также добавлены: начальная обработка данных, кластеризация слов-кандидатов и возможность уточнения мотивов.

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

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

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

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

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

Размещено на Allbest.ru