Содержание к диссертации
Введение
Глава 1. Первичный анализ регуляции генов 14
1.1. Основные понятия и определения . 14
1.2. Создание базы знаний о факторах и предсказание сайтов связывания 19
1.3. Обработка экспериментальных данных и предсказание регулирующих факторов 30 CLASS Глава 2. Модель регуляторного комплекса 34 CLASS
2.1. Постановка задачи 34
2.2. Построение модели 37
2.2.1. Конструирование и анализ алгоритма оконного класса моделей 37
2.2.2. Конструирование булева класса моделей 44
2.3. Целевая функция 49
2.4. Поиск оптимального комплекса 54
2.4.1. Выбор эффективного алгоритма поиска 54
2.4.2. Операторы создания и мутации 58
2.4.3. Оператор кроссовера , 59
2.5. Обобщённый класс моделей: унификация знаний о сайтах и факторах 62
2.5.1. Общая структура 62
2.5.2. Подкомплекс: семантика и структура 64
2.5.3. Вес обобщённого комплекса , 73
2.5.4. Эффективная программная реализация математической модели обобщенного комплекса 76
2.5.5. Операторы создания, мутации и кроссовера в обобщённом классе 81
2.6. Обобщённая целевая функция 86
2.6.1. Компоненты целевой функции ..86
2.6.2. Ограничения обобщённой целевой функции 97
Глава 3. Взаимодействие с пользователем и оценка качества 100
3.1. Вычисление комплекса и вывод результата 100
3.2. Методы оценки качества результата 103
3.2.1. Мультизапуск: проверка устойчивости поведения генетического алгоритма ,.103
3.2.2. Запуск с кластеризацией 114
3.2.3. Запуск с расщеплением выборки: достоверность результата и переобучение .117
3.3. Значимость компонентов подкомплекса 117
Глава 4. Реализация и тестирование 121
4.1. Реализация СМА 121
4.2. Тестирование СМА 128
4.2.1. Тестирование на искусственных данных ...128
4.2.2. Тестирование на экспериментальных данных І36
4.3. Система ExPlain 142
Заключение 148
Список литературы
- Создание базы знаний о факторах и предсказание сайтов связывания
- Конструирование и анализ алгоритма оконного класса моделей
- Мультизапуск: проверка устойчивости поведения генетического алгоритма
- Тестирование на искусственных данных
Введение к работе
Согласно [1], современный путь создания лекарства в общем случае включает следующие этапы: выбор молекулярной мишени для действия лекарства, нахождение базовой структуры нового лекарства, оптимизация базовой структуры, доклинические испытания, клинические испытания, производство препарата. При этом в течение длительного времени на стадиях экспериментального тестирования подавляющее большинство прототипов лекарств отбрасывались как бесперспективные по ряду причин: низкая активность, высокая канцерогенность, сложность производства и т. д. Приблизительно одно соединение из 100 000 выходило на фармацевтический рынок. В [2] отмечено, что производство нового лекарства требует 12-15 лет и свыше 800 млн. долларов.
В последнее десятилетие ситуация изменилась. Благодаря внедрению компьютерных технологий и достижениям биоинформатики первые три этапа были рационализированы: поиск новых прототипов происходит более направленно, и значительно больше заведомо бесперспективных соединений отбрасывается на ранних этапах. Хотя, безусловно, компьютерное моделирование не заменит собой клинические испытания, согласно [1] уже сегодня оно позволяет снизить на два порядка количество прототипов, которые необходимо синтезировать и проверить.
Данная работа связана с выбором молекулярной мишени. Лекарство призвано активировать или дезактивировать какой-либо белок в организме человека либо в микроорганизмах. Во многих случаях конечной целью является воздействие на работу генов. К примеру, в результате болезни или плохой наследственности с некоторого гена не производится белок или, напротив, производится чрезмерно много белка. В этих случаях цель лекарства — восстановить нужный уровень активности. Поэтому зачастую мишенями становятся транскрипционные факторы (белки, регулирующие работу генов в ор-
ганизме) либо белки, передающие сигналы факторам. Однако для того, чтобы выбрать оптимальные мишени, необходимо тщательно изучить регуляторные процессы, действующие в клетках, характер их отличий в здоровом и больном организме (либо в здоровых и больных клетках), динамику изменений регуляции во времени и т. д.
Таким образом, сегодня перед исследователем стоит актуальная задача: определить набор активных элементов — факторов, которые регулируют тот или иной процесс в организме, а также характер взаимодействия между ними, отличить нормальную регуляцию от нарушенной и по возможности найти пути воздействия на этот процесс, способы восстановить регуляцию в больном организме. На данный момент биоинформатика ещё далека от того, чтобы дать точные и ясные ответы на такие вопросы в большинстве случаев. Не открыто ещё общих законов, которые достаточно хорошо описывали бы регуляцию генов и при этом работали бы в большинстве ситуаций.
Тем не менее, обнаружено достаточно много частных закономерностей и накоплен огромный экспериментальный материал. Это позволило нам создать комплексный инструмент, который может значительно помочь исследователю в поиске ответов на его вопросы. Хотя разработанные нами программные средства не дают абсолютно точного ответа, умелое комбинирование созданного нами инструментария с экспериментальной работой и творческий подход к решению задачи позволяют решить её в самых различных вариациях. Иногда можно найти регулирующие элементы исключительно экспериментально, не прибегая к компьютерному моделированию, но с помощью последнего это можно сделать быстрее и дешевле. Так, вместо проверки всевозможных факторов разработанное нами программное обеспечение позволяет смоделировать процесс регуляции и получить подмножество факторов, которые действуют в данном эксперименте с наибольшей вероятностью, после чего проверять уже только их.
Хотелось бы также отметить потенциальную важность исследования генной регуляции не только для генетики, но и для изучения новых вычислительных технологий и принципов построения алгоритмов и программ. Рассматривая представление ДНК как блока данных, а транскрипционных факторов — как программы, обрабатывающей эти данные, и изучая механизмы их взаимодействия, можно получить новые эффективные методы построения программ или решения некоторого класса задач. Напомним, что, исследуя устройство и функционирование генетического кода, учёные уже описали семейство эволюционных алгоритмов (включая генетический алгоритм), базирующихся на тех же принципах. Также можно упомянуть, что в 1994 году благодаря работе Эдлмана [3] появился метод решения NP-полных задач с помощью фрагментов ДНК и так называемые биокомпьютеры [4]. Изучение законов генной регуляции, возможно, позволит в будущем развить эту область и увеличить возможности биокомпьютеров1.
Цель работы заключается в разработке и совершенствовании математических и программных методов для анализа работы регуляторной системы при взаимодействии с промоторами генов в общем случае и в отдельных частных экспериментах. Основные задачи в рамках достижения этой цели включают:
Изучение и формализацию знаний, накопленных в области генной регуляции, описание поведения регуляторной системы в виде алгоритмов.
Построение параметризованной модели регуляторного комплекса, включающей в себя набор транскрипционных факторов и характер их взаимодействия.
Мнение о потенциальной пользе нашей работы для биокомпьютерных исследований высказал Danny van Noort (PhD, ранее старший научный сотрудник кафедры биомолекулярной обработки информации Института Фраунгофера, Германия, ныне профессор Национального Университета Сеула, читает курс лекций по биокомпьютерным вычислениям) на конференции ICNC05 в Китае. Уже сейчас проводятся эксперименты по использованию генной регуляции в биокомпьютерах.
Реализацию способа оптимизации (подбора параметров) этой модели для достижения наибольшего соответствия модели экспериментальным данным.
Разработку программной системы, упрощающей и автоматизирующей обработку экспериментальных данных, связанных с генной регуляцией.
Создание программного инструментария, который выполняет процесс оптимизации и представляет подробные результаты, а также позволяет гибко управлять видом модели и процессом её оптимизации и оценить качество результата.
Проектирование и реализация методов тестирования программных систем на искусственных и экспериментальных данных для оценки качества полученных результатов.
В результате проведённой работы разработано два программных средства — СМА и ExPlain. СМА [66,68,69,70,75] представляет собой неинтерактивное приложение, основанное на библиотеке GRESA [5] и управляемое множеством параметров командной строки. Большая часть проектирования, разработки и отладки СМА выполнена автором данной работы, но в создании этого продукта принимали участие и другие разработчики. ExPlain [77,80] — это мощная интерактивная система анализа генной регуляции в целом, которая, в частности, является оболочкой для СМА, хотя удобный запуск СМА и визуальное представление результатов этого запуска — лишь малая доля функциональности этой системы. Она сопряжена с различными базами данных, выпускаемыми BIOBASE GmbH, и призвана сделать обработку данных по генной регуляции более эффективной. ExPlain разрабатывается автором этой работы в составе интернациональной группы, в которой около 10 человек.
Существуют некоторые конкурирующие с СМА разработки, такие как TOUCAN [6,7] и TeLIS [10]. В работе [73] мы проводили сравнение их функциональности с СМА, но с тех пор проект СМА ушёл далеко вперёд, хотя некоторые полезные особенности конкурентов до сих пор не имеют аналогов в СМА (например, разбиение весовых матриц на классы с возможностью исключения из регуляторного комплекса матриц из одного класса). С другой стороны, СМА содержит много уникальных особенностей, таких как предсказание композиционных элементов, комплекса, состоящего из нескольких подкомплексов, учёт нескольких сайтов одной матрицы и многое другое.
Система ExPlain на мировом рынке аналогов не имеет. Хотя для большинства отдельных компонентов системы существуют подобные конкурирующие разработки, до сих пор никто не пытался интегрировать их в цельную систему.
Некоторые идеи, методы и алгоритмы, описанные в работе, нашли также применение в программном комплексе CisSearch [71,72,76,78,79,81], представляющем собой интерактивное графическое приложение под Windows, которое облегчает обработку разнообразных экспериментальных данных, связанных с генной регуляцией. Впрочем, мы практически не будем говорить об этом продукте.
Практическая ценность
Разработанные программные средства успешно использовались различными специалистами для обработки данных по генной регуляции включая представителей международных биотехнологических и фармацевтических компаний (например, Serono Group); научных институтов и центров (например, Германского центра изучения рака, DKFZ); компании, занимающейся выпуском биологических баз данных BIOBASE GmbH.
Апробация работы
Результаты работы докладывались на международных научных конференциях, включая Ist Intl. Conf. on Natural Computations (ICNC05) в г. Чанша, Китай; German Conf. on Bioinformatics (GCB'05) в г. Гамбург, Германия; 3rd Annual RECOMB Satellite Workshop on Regulatory Genomics в г. Сингапур, Сингапур; 3rd Intl. Conf. "Genomics, Proteomics, Bioinformatics and Nanotech-nologies for Medicine" в г. Новосибирск и др. Работа была представлена на рабочем семинаре «Наукоёмкое программное обеспечение» конференции памяти академика А. П. Ершова «Перспективы систем информатики», на различных встречах, семинарах. Система ЕхРїаІп демонстрировалась на пленарных докладах международных конференций, на встречах с представителями свыше десятка биотехнологических и фармацевтических компаний.
Структура и объем работы
Диссертационная работа состоит из введения, четырёх.глав, заключения и списка литературы. Объем диссертации — 163 стр. Список литературы содержит 81 наименование. Работа включает 26 рисунков и графиков, полученных в результате расчётов на ЭВМ, в том числе с использованием разработанного программного обеспечения.
В первой главе вводятся некоторые необходимые понятия и описан процесс первичной обработки данных — поиск сайтов связывания транскрипционных факторов на промоторах. Раздел 1.1. посвящен описанию некоторых необходимых биологических понятий и их математической формализации.
В разделе 1.2. описана проблема предсказания сайтов связывания транскрипционных факторов и приведён один из способов её решения, основанный на весовых матрицах. По известным экспериментально сайтам составляется модель, которая представляет собой матрицу размера 4 х L, где L — длина сайта (количество нуклеотидов). Каждой позиции внутри сайта для
нуклеотидов А, С, G, Т соответствует вероятность появления данного нуклео-тида в данной позиции сайта с некоторой перенормировкой. Имея весовую матрицу, для любой последовательности длины L можно вычислить степень соответствия (вес) между последовательностью и моделью; в случае, если вес превышает некоторый порог, последовательность можно считать предсказанным сайтом, соответствующим данной матрице. Далее раздел включает описание работы программы Match, созданной ранее с использованием библиотеки GRESA. Так как почти всегда в качестве входных данных СМА принимает результат работы Match, эта часть является значимой для понимания дальнейшего материала, поэтому она описана достаточно подробно.
В конце первой главы приведён простейший способ предсказания регулирующих факторов по значениям экспрессии, полученным экспериментально. Промоторы с низкой экспрессией помещаются в категорию 'No', а промоторы с высокой — в категорию 'Yes', после чего находится отношение плотности предсказанных сайтов в этих категориях для каждой матрицы. В разделе приведён пример экспериментальных данных, где для ряда матриц некоторых факторов такое отношение весьма заметно (то есть предсказанные сайты этих матриц встречаются в 'Yes' значительно чаще, чем в 'No'), это даёт возможность предположить, что в данном эксперименте эти факторы действительно регулировали активность.
Вторая глава содержит основной объем теоретического материала. В ней описана общая постановка задачи, решаемой СМА, вводится и подробно описывается регуляторная модель, описывается способ её оптимизации, а также целевая функция. Часть изложенной теории оформлена в виде алгоритмов, затронуты вопросы быстродействия и оптимальной реализации.
В разделе 2.1. описана общая постановка задачи, решаемой программой СМА. Вводится понятие модели, как параметризованной функции J?5:57 —> [0,1], где SI — множество предсказанных сайтов для некоторого
промотора. Модель задаёт набор правил и определяет, насколько данный промотор (а точнее — список сайтов, предсказанных на нём) соответствует этим правилам. Приводится пример простейшей модели с одним параметром:
RSv(Sl) =
1, если 3 сайт матрицы X в SI О, иначе
Здесь параметр X— имя матрицы. Такой модели соответствуют промоторы, на которых присутствует как минимум один сайт матрицы X. При оптимизации параметров модели получается регуляторный комплекс, который наилучшим образом соответствует экспериментальным данным. В качестве алгоритма оптимизации был выбран генетический алгоритм, а в качестве целевой функции — линейная комбинация различных компонент, включающих следующие: соответствие профиля вычисленных значений комплексов RS экспериментальным значениям экспрессии (2Л); степень различия значений RS для промоторов категорий 'Yes' и 'No' по критерию Стьюдента (2Г); количество ошибок перепредсказания и недопредсказания {ZF)\ соответствие распределения значений RS для категории 'Yes' нормальному (Zv); степень сложности комплекса (Zp). В разделах 2.2. и 2.5. детально описываются различные классы (типы) моделей. В разделах 2.3. и 2.6. вводится целевая функция. В разделе 2.4. обосновывается выбор алгоритма оптимизации и приводятся операторы мутации и кроссовера.
Третья глава посвящена особым вспомогательным возможностям, заложенным в СМА, которые позволяют не просто получить оптимальный комплекс, но и провести дополнительные исследования: оценить надёжность результата и склонность к переобучению для заданных параметров и набора входных данных; оценить важность отдельных компонентов полученного комплекса; получить детальную информацию как по наилучшему комплексу, так и по другим хорошим комплексам из популяции или по комплексу, сконструированному вручную. Один из методов оценки результата — это мульти-
запуск, при котором оптимальный комплекс находится генетическим алгоритмом несколько раз, и затем результаты сравниваются и определяется надёжность по следующей формуле:
Здесь Syc^cA — функция схожести двух комплексов, полученных в
результате соответствующих запусков генетического алгоритма. Она принимает значения от 0 до 1, причём большие значения соответствуют большей степени схожести. В главе конструируется эта функция для всех классов моделей, введённых во второй главе, после чего приводится пример использования надёжности.
Четвёртая глава посвящена реализации и тестированию. В ней приведена внутренняя структура СМА с точки зрения программирования; структура основных пакетов, краткое описание классов. В раздел 4.2. включены некоторые результаты тестирования на искусственных и реальных данных. Результаты тестирования на искусственно сгенерированных промоторах, в которые были внедрены сайты определённых факторов и композиционных элементов, показали, что СМА успешно находит внедрённые факторы. Результаты тестирования на двух экспериментальных наборах данных, для которых правильный результат частично известен, показали, что комплекс, находимый СМА, хорошо согласуется с известными данными. Кроме того, описана система ExPlain, которая также создавалась при участии автора. Она интегрирует пользовательские экспериментальные результаты с базами данных биологической информации и позволяет проводить различные типы анализа, преобразуя необходимым образом входные и выходные данные.
В заключении описаны основные результаты работы, вклад автора в проделанную работу и идеи дальнейшего развития проекта.
Для описания алгоритмов в тексте работы используется псевдокод. Мы не используем какой-либо реальный язык программирования, чтобы не отвлекаться в описании алгоритмов на формальности (к примеру, объявление переменных), а сосредоточиться на главном, Большинство конструкций псевдокода должны быть интуитивно понятны любому программисту. Алгоритм начинается с заголовка, состоящего из слова alg, названия алгоритма и набора параметров. Сам алгоритм и блоки операторов циклов заключены в операторные скобки begin...end. Конструкция if, циклы for и while, оператор преждевременного выхода из цикла break, оператор присваивания := подобны существующим в языке Pascal. Возвращаемое алгоритмом значение задаётся с помощью оператора return, как в языке Си. Помимо скалярных типов данных существуют списки, для которых определены операции вставки элемента и списка элементов того же типа (insert...into), удаления элемента (delete...from), получения первого элемента (pop), цикла с перебором элементов (foreach...in) и проверки на отсутствие элементов (empty). Набор переменных, объединённый в круглые скобки, соответствует записи (record) в Pascal или структуре (struct) в Си. Для упрощения сама переменная-структура обычно не именуется, а обращение к полям производится, как к независимым переменным. Кроме того, существуют массивы. Под массивом понимается список элементов, доступ к которым осуществляется по ключу. Массив обозначается символом с индексом в фигурных скобках: \Р;\, а его элемент, соответствующий ключу / — тем же символом, но без фигурных скобок: рг. Добавление элемента в массив выполняется через присваивание вида pt :=1; если элемента рг ещё не было, то он добавляется в массив со значением 1,
иначе значение обновляется. В программном коде такой массив может быть реализован как хэш, AVL-дерево или массив в общепринятом смысле (если ключи можно представить как последовательные целые числа). В начале выполнения алгоритма скалярные численные переменные, не являющиеся вход-
ными параметрами, предполагаются равными нулю, а списки и массивы — не содержащими элементов. С символа # начинается комментарий, который обычно объясняет параметры алгоритма.
Создание базы знаний о факторах и предсказание сайтов связывания
Знание о том, какие сайты связывания присутствуют на том или ином промоторе, является весьма полезным для понимания процессов регуляции того или иного гена. Однако экспериментальный поиск сайтов in vitro — это достаточно трудоёмкий и дорогостоящий процесс. Так, в базе данных TRANSFAC 9.4 (декабрь 2005 года) зафиксировано 2409 экспериментально подтверждённых сайтов для Homo sapiens, что сравнительно немного. В ДНК человека 15000-20000 генов, то есть в среднем известен один сайт на 6-8 генов. Как правило, один открытый сайт — это одна статья в журнале уровня Nucleic Acids Research или Journal of Biological Chemistry либо дипломная работа одного студента. В действительности количество сайтов связывания выше на порядок, в связи с чем становится актуальной задача предсказания сайтов in silico. Такой метод быстрее и дешевле, хотя существующие на сего дняшнии день алгоритмы не могут предсказывать сайты надёжно: количество ошибок первого и второго рода достаточно велико.
Мы рассмотрим хорошо известный метод предсказания сайтов, основанный на весовых матрицах [18]. Хотя этот метод был реализован без участия автора, он будет описан достаточно подробно, так как он является фундаментом для дальнейшего анализа. Метод также был описан нами в [67]. Реализация этого метода выполнена в программе Match [21], которая используется на первом этапе дальнейшего анализа регуляторных областей. Существуют и другие методы предсказания сайтов, к примеру, основанные на построении меры информации [19], нейронных сетях [20] и т. д. Мы не будем останавливаться на этих методах.
Метод весовых матриц позволяет построить модель сайтов определённого фактора и затем для произвольной подцепочки определить, является ли она сайтом данного фактора. Пусть предикат BSn (х),х є S определяет, является ли подцепочка х сайтом связывания фактора TF. Метод основан на следующих предположениях:
Значение BSTF(x) зависит только от последовательности подцепочки и не зависит от положения и стрэнда х, а также какой-либо другой информации (например, нуклеотидов, окружающих данную подцепочку).
Чем больше похожа последовательность подцепочки х на некоторый образец, тем больше степень уверенности, что х — сайт.
Вообще говоря, эти предположения могут и не быть верными, однако они позволяют построить довольно простую и вместе с тем биологически обоснованную модель.
Для примера рассмотрим несколько сайтов фактора AhR (агу! hydrocarbon (dioxin) receptor), известных из экспериментальных данных. Так как мы решили, что положения и стрэнды не представляют интереса, рассмотрим только последовательности нуклеотидов этих сайтов (слева приведены идентификаторы из базы TRANSFAC):
Заметим, что в результате экспериментов установлено лишь, что внутри каждого из этих фрагментов есть сайт связывания, то есть, вероятно, не все нуклеотиды из приведённых последовательностей участвуют в связывании регуляторной области с фактором, а только некоторая подпоследовательность. Процедура выравнивания, в детали которой мы углубляться не будем, смещает последовательности друг относительно друга таким образом, чтобы пересекающиеся части были наиболее похожи. К некоторым последовательностям также применена операция комплементарности (они выделены курсивом). В результате получается следующее:
Конструирование и анализ алгоритма оконного класса моделей
Из предметной области известно, что факторы нередко действуют в группе и увеличивают экспрессию только в том случае, если все необходимые факторы оказались связаны с промотором, причём сравнительно недалеко друг от друга. Мы воспользуемся этой информацией для конструирования оконного класса моделей Q.
Назовём окном ширины / диапазон в цепочке нуклеотидов, в котором должны находиться все сайты, учитываемые при подсчёте значения модели RS. Будем рассматривать не более одного сайта для каждой матрицы (выберем сайт с максимальным значением rw, если их несколько). Общее количество матриц, значимых для модели, обозначим N. Будем считать, что матрицы не повторяются. Зададим функцию Ex(Sl,id,c), похожую на (2.2): Ex(Sl,id,c) = \,eaiu3p,rw,cw, 4morw cu{p,rw,cw,id} ESl (2-3) О, иначе Модель без учёта окна для заданного N может выглядеть следующим образом: SV vC S/) = l ( ,,C,) (2.4)
Здесь следует объяснить необходимость вынесения в параметры модели значений порогов С,. Дело в том, что хотя пороги для матриц, использованные в профайле при выполнении Match, имеют под собой некоторое основание, по большому счёту исходные соображения могут быть и ошибочными. Оптимизировав же пороги в результате поиска комплекса, мы можем не только найти комплекс, который лучше соответствует данному эксперименту, но и использовать полученные пороги в дальнейшем в других профайлах, обосновывая это тем, что именно эти значения улучшают разделение в некотором эксперименте. Кроме того, из-за плохо выбранных порогов для Match может оказаться, что некоторой матрице соответствует чрезмерно большое количество сайтов; соответственно, она присутствует на каждом промоторе, и её наличие в комплексе никак не способствует разделению промоторов.
Подбор порога для веса ядра не проводится, так как он обычно в профайле достаточно велик, а во многих случаях просто равен единице, что означает требование полного соответствия ядра во всех сайтах.
Введём предикат InSIM[p),который принимает подцепочку р, позиции начала и конца st, en и возвращает истину в случае, если подцепочка лежит в данном диапазоне: , , истина, если Pos[p) st и Pos{p) + Len(p)-[ en І ложь, иначе Добавим в (2.3) проверку на попадание сайта в заданный диапазон с помощью функции In: Ex (Sl,id,c,st,en) = 1, если 3p,nv,cw:nv c, Inslen(p) и (p,rw,cw,id)eSl (2.6) О, иначе
Пусть Lmn и Z, — диапазон координат, в котором находятся все сайты данного промотора, то есть V(p,nv,cw,/rf)eS/ In, , (р). В качестве этих значений можно взять границы промотора (если исходный промотор известен как подцепочка р, то Lmm =Pos(p);Lmux = Pos(p) + Len(p)-\). Но, вообще говоря, модели необязательно знать о промоторе, ей необходим только набор сайтов. Поэтому эти значения можно получить как максимальную и минимальную позиции сайтов из SI: puSl (2.7) = ? Р0$ (р}+ Len р} [ Теперь мы можем полностью записать оконную модель:
Параметры N и / будем считать параметрами класса Q. Заметим, что при / = Ьтгх - 1тш +1 (2.8) превращается в (2.4), то есть модель без учёта окна является частным случаем модели с окном. Запишем алгоритм вычисления оконного комплекса по (2.8):
Проанализируем этот алгоритм. В первом цикле вычисляются значения Lmn и Lmax. Затем для каждого положения окна вычисляется Ех для всех матриц модели и выбирается положение окна, при котором сумма Ех максимальна.
Как видно, в таком виде сложность алгоритма определяется количеством итераций наиболее вложенного цикла и составляет ОШ\Sl\-{Lnm-Lmm l) \. Такая сложность неудовлетворительна, учитывая, что даже для одного комплекса RS вычисляется для всех промоторов. Попытаемся оптимизировать этот алгоритм.
Чтобы избавиться от необходимости многократно пробегать по списку 5/, преобразуем его в массив списков {Slld}, где ключом является имя матрицы, к которой относится тот или иной сайт. Кроме того, вычислим заранее значениями
Мультизапуск: проверка устойчивости поведения генетического алгоритма
Этот вариант возник, когда мы задались следующими вопросами. Каким образом определить надёжность найденного комплекса? Достаточно ли поколений отработал генетический алгоритм, хватает ли размера популяции? Насколько найденное решение далеко от глобального максимума?
Проще всего ответить на эти вопросы, если заранее известен оптимальный ответ. К сожалению, его у нас нет, а если бы он был, то необходимость в использовании генетического алгоритма отпала бы вместе с вопросами.
Тем не менее, мы можем сделать некоторые выводы, если выполним запуск генетического алгоритма несколько раз и сравним найденные результаты. Так, сильное различие в полученных оптимальных комплексах скажет о том, что результат неустойчив, а значит, необходимо либо увеличивать размер популяции и количество поколений, либо упрощать модель, так как она перегружена параметрами и не может дать устойчивого решения за разумное
время. Напротив, похожие или одинаковые результаты покажут, что параметры запуска алгоритма подобраны достаточно хорошо, чтобы найти оптимальный комплекс. Хотя, разумеется, остаётся вероятность, что генетический алгоритм постоянно оказывается заперт в одном и том же локальном минимуме.
Для сравнения найденных комплексов можно сравнивать их веса. Дополнительно мы решили, что представляет интерес введение численной меры для сравнения двух комплексов одинаковой структуры. Пусть дана функция сходства S[c]tc2), которая принимает два комплекса (точнее, набор параметров модели для этих комплексов) и возвращает значение от 0 до 1, причём 5(с,,с:) = 1 означает, что комплексы в точности совпадают, а S(c,,c2) = 0 означает, что между ними нет ничего общего. Потребуем также симметричности: 5(с,,с2) = 5(с2,с,). Предположим, мы прогнали генетический алгоритм п раз и получили в результате п комплексов cvc2,...,cn —оптимальные комплексы каждого прогона. Определим сходство /-го результата с прочими посредством усреднения сходства с, со всеми прочими комплексами: [ _ , , 5,=-2: s(Cl,c,) (за) П lJ = lr J-ll + l, ,Я
Вычислим теперь величину R, как усреднённое значение 5( (воспользовавшись симметричностью Syc Cjj): й45 =та с ) (3-2) Я 1=1 п\п l) i 2 j=l
Назовём R мерой надёжности полученного результата. Она находится в [0,1], причём большие значения соответствуют большей надёжности.
Способ вычисления S{cvc2) зависит от класса моделей. Рассмотрим их по очереди.
Пусть даны наборы параметров с, -ІХ1п..,,Х\,СІ,...,С1Л и с2 = lxf,...,X\,,Cf t...,Cl- J, которые описывают два комплекса класса Q. Введём функцию (с,,с2), которая определяет вклад в сходство матрицы X): фі Сг) = - , если 3/, Ч/ИО X = X, 2(і-яшг,) (3.3) О, иначе
Здесь mss , — порог матрицы Jf, в профайле (1-mjw , — максималь но возможное различие порогов для этой матрицы). Если такая матрица есть в с2, то st принимает значение от 0,5 до 1 в зависимости от различия порогов, в противном случае st=0. S(cvc2) будет усреднённым значением st{cvc2): 1 ,v v 1=1 S{cvc2) симметрична в силу того, что внутри комплекса оконного класса матрицы не повторяются. Значения этой функции лежат в [0, l], что и требовалось.
Так как мы сконцентрировались на работе над обобщённым классом, то в булевом классе S{cvc2) была проработана плохо. Здесь 5(с],с2) = 1-А,/тах(«1,й2), где щ и и2 — количество непустых групп в сх и с2% а к — количество совпадающих групп в комплексах, причём группа G0 в сх сравнивается с группой G0 в с2, а для всех прочих групп из с, перебираются группы с2 в поисках такой же. Сравнивается только набор матриц в группе, значения порогов игнорируются полностью.
Тестирование на искусственных данных
Для проверки работоспособности СМА на начальной стадии мы создали генератор искусственных подцепочек с внедрёнными на них сайтами. Генератор действует следующим образом: изначально создаётся набор случайных последовательностей, после чего на части из них (категории Yes ) вставляются сайты некоторых заданных матриц, а точнее подцепочки, которые могут быть в дальнейшем предсказаны с помощью Match как сайты этих матриц. Для увеличения сходства сгенерированных промоторов с настоящими генератор может принять на вход набор настоящих промоторов и на их основании сформулировать марковскую модель порядка 0-3, после чего генерировать промоторы в соответствии с созданной моделью. Марковские модели активно применяются для анализа закономерностей в геноме [56]. Сайты внедряются с заданной частотой и с заданным минимальным весом rw в окне заданной ширины. Также есть возможность внедрения пар с заданной ориентацией и расстоянием между сайтами. По сути дела, предоставлена возможность внедрения произвольных комплексов обобщённого класса, состоящих из одного подкомплекса. Количество сайтов одной матрицы, внедряемых в окно, не фиксировано, а варьируется в соответствии с нормальным распределением.
Для сгенерированных последовательностей создаётся также файл со значениями экспрессии, причём средняя экспрессия на промоторах категории Yes превышает экспрессию на промоторах категории No (хотя диапазоны могут перекрываться; степень перекрытия задаётся пользователем). После этого СМА предлагается найти внедрённый комплекс по сгенерированным последовательностям и значениям экспрессии.
Для тестирования мы взяли достаточно большой профайл VertGood-MinSum, входящий в поставку базы данных TRANSFAC. В него входит 358 наиболее достоверных матриц, сконструированных по факторам позвоночных.
В качестве порогов mss и ess в нём используются значения, минимизирующие сумму ошибок недопредсказания и перепредсказания Match. Было сгенерировано 440 искусственных регуляторных подцепочек длиной 1100 нуклео-тидов каждая, половина из которых была отнесена к категории Yes и содержала внедрённые сайты. Для генерации использовалась модель Маркова третьего порядка, построенная по набору промоторов человека. Для начала мы внедрили следующий комплекс из трёх матриц: Kl=[V?aRPl_01@0.88/4V?MYOGNFl_0ie0.77/21V$ZID_0ie0.91/2] RS=K1
Матрицы и пороги были выбраны методом проб и ошибок с таким расчётом, чтобы их сайты могли случайно возникать на последовательностях, где не были внедрены специально. В противном случае в процессе генетического алгоритма предпочтение было бы вполне логично отдано комплексу, состоящему только из одной матрицы внедрённого комплекса, как более простому, но хорошо разделяющему выборки. Мы же хотим сгенерировать последовательности, где важны все составляющие внедрённого комплекса.
В первом эксперименте значения экспрессии в категории Yes были строго выше значений экспрессии в категории No .
При поиске оптимального комплекса мы задали в параметрах класса наличие строго одного подкомплекса и запретили композиционные элементы. Однако количество матриц было установлено достаточно свободно: Nmmn = 1; ty«wg =3; Wmmax -5. Значения v мы ограничили диапазоном от 1 до 4. После прогона Match выяснилось, что у каждой матрицы в среднем 80 уровней порога, и общее число различных комплексов заданной модели превышает 1025.
Веса компонентов целевой функции мы выбрали так: а = 0; 6 = 5/11; с = 5/11; d - 0; е = 1/11. Мы не учитываем компоненту нормальности Z v, так как сгенерированные комплексы заведомо не удовлетворяют нормальному распределению. Соответствие профиля экспрессии профилю RS также трудно ожидать, потому что внутри категорий экспрессия генерировалась как равномерно распределённая случайная величина. По этой причине ZR также не учитывалось. Влияние штрафной компоненты выставлено достаточно слабым, так как можно предположить, что любая матрица из комплекса по отдельности даст неплохое значение Zr и ZE, и большой вклад Zp не позволит найти комплекс из всех трёх матриц.
Сначала мы высчитали значение внедрённого комплекса с помощью CMAcalc. Получилось Z = 0,9624 , причём значения компонент составили ZT = 0,9357; ZE = 1,0; Z = 0,7445. Хотя сайты отдельных матриц из комплекса встречаются на 16% промоторов категории No , они не образуют комплексов, поэтому мы имеем оптимальное разделение и полное отсутствие ошибок перепредсказания и недопредсказания. Затем мы запустили генетический алгоритм на 200 поколений с популяцией из 200 хромосом.
Похожие диссертации на Алгоритмы и программный инструментарий для исследования процессов генной регуляции
