Реконструкция генных сетей на основе данных микрочиповых экспериментов: отбор генов мишеней и сравнение различных моделей регуляции

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

1. Генная экспрессия, основные определения

2. Микрочиповые эксперименты

3. Выявление генов изменивших уровень экспрессии

4. Гипергеометрический метод

5. Анализ данных по раку молочной железы

6. Моделирование генной экспрессии

7. Линейная модель

8. Методы сглаживания

9. Стохастическая модель

10. Сравнение моделей

11. Программная реализация в системе biouml

ЗАКЛЮЧЕНИЕ

СПИСОК ЛИТЕРАТУРЫ

ПРИЛОЖЕНИЕ 1

ПРИЛОЖЕНИЕ 2



ВВЕДЕНИЕ

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

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

Одним из наиболее развитых подходов для решения этой проблемы в современной биохимии является математическое моделирование.

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

Целью данной работы являются:

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

Задачи, решаемые в работе:

· Рассмотрение и сравнение различных математических методов анализа данных микрочиповых экспериментов, предложение своих методов.

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

· Написание программного продукта, реализующего рассмотренные методы.



1. Генная экспрессия, основные определения

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

Ген -- это участок ДНК, кодирующий одну полипептидную цепь (белок) или одну молекулу тРНК, мРНК или рРНК. Это - смысловая единица генетического кода.

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

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

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

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

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

Рисунок 1. Процессы генной экспрессии и регуляции

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

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

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

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

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

2. Микрочиповый эксперимент, исходные данные

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

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

С помощью микрочипов можно изучать:

- экспрессию генов в различных тканях.

- изменение экспрессии генов в ходе клеточного цикла и в ходе развития.

- - экспрессию генов в норме и при патологии (в нормальных и раковых клетках).

- изменение экспрессии генов с течением времени как результат внешнего воздействия (взаимодействие клетки с патогеном, лекарством).

- наличие чужеродной ДНК/РНК вирусов для диагностики различных заболеваний (например, ЗППП)

- анализ мутаций и полиморфизмов.

Основываясь на характере числовых данных в чипе, можно построить две независимые классификации микрочипов:

- статические (данные замерены единовременно).

- динамические (данные распределены во времени).

А также

- с конкурентной гибридизацией (тестовый образец против контроля).

- без конкурентной гибридизации (только тестовый образец).

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

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

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

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



3. Выявление генов изменивших уровень экспрессии

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

Напомним, что исходными данными являются матрицы значений экспрессии генов:

Строка соответствует отдельной пробе (гену, либо его участку).

- экспрессия i-той пробы в условиях эксперимента j. Основными типами данных являются повременные замеры экспрессии ( - экспрессия в момент времени ) и замеры у разных пациентов ( - экспрессия у j-того пациента).

Задача ставится следующим образом: необходимо выделить из общей группы те пробы, которые значимо изменили свою экспрессию либо с течением клеточного цикла, либо в результате воздействия на пациентов. Часто данные представляют собой отношение экспрессии объекта у пациента к экспрессии данного объекта в норме.



В этом случае задача сводится к выявлению того, насколько значимо отклоняется от единицы в ту или иную сторону.

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

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

Замер экспрессии у j-того пациента, соответственно - реализация случайной величины.

В случае повременных замеров имеем случайный процесс:

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



(

либо сравнить экспрессию в норме и патологии, при этом нужно иметь несколько реализаций объекта в норме для успешного статистического анализа. Т.е. здесь сравниваются две выборки, и оценивается различие между ними. Рассмотрим некоторые стандартные методы, которые могут применяться к этой задаче.

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

1. Статистика Стьюдента

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

- первая выборка значений экспрессии для гена,

- вторая выборка значений экспрессии для гена,

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

где

,

-- соответствующие выборочные средние,

,

-- соответствующие выборочные дисперсии,

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

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

Таким образом, используя двухсторонний критерий, P-значение вычисляем по формуле:

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

Окончательно, получаем:



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

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

2. Критерий ранговых сумм Вилкоксона

Свободный от распределения критерий суммы рангов для гипотезы об отсутствии эффекта обработки.

Имеем наблюдений и .

Положим:

где и -- наблюдаемые, а -- ненаблюдаемые случайные величины; интересующий нас параметр -- это неизвестный сдвиг в положении, обусловленный “обработкой”.

Все случайных величин взаимно независимы.

Все извлечены из одной генеральной совокупности.

Проверяем гипотезу , для этого:

Объединим выборочные значения обоих выборок и расположим их в порядке возрастания, т.е. образуем общий вариационный ряд, и каждой величине из этого ряда сопоставим ее ранг , равный порядковому номеру величины в общем вариационном ряду. Заметим, что если нулевая гипотеза имеет место, то любое распределение рангов по этим двум выборкам равновероятно.

В качестве статистики критерия берется сумма рангов одной (например, первой выборки):

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

Отсюда промежуточное P-значение вычисляем по формуле:

Далее, если , то в качестве итогового P-значения мы берем:

Если же наше , тогда подсчитываем все различные способы группирования рангов, при которых статистика принимает значения, равные или больше наблюдаемого , аналогично с предыдущим подсчитываем итоговое P значение:



3. Тест Колмогорова

Свободный от распределения критерий для проверки однородности двух выборок .

Имеем наблюдений и .

Предположим, что:

Все наблюдений и взаимно независимы.

Все эти извлечены из одной непрерывной совокупности .

Все извлечены из одной непрерывной совокупности .

Рассматриваем гипотезу о том, что совокупности и идентичны, т.е. о том, что обе выборки извлечены из одной и той же совокупности. Её можно переписать так:

Будем тестировать эту гипотезу против альтернативы

.

Для проверки делаем:

Будем считать, что наши выборки упорядочены по возрастанию:

(если это не так просто упорядочиваем их).

Статистика критерия Колмогорова-Смирнова измеряет различия между эмпирическими функциями распределения, построенными по выборке:

Однако на практике значение данной статистики удобнее вычислять при помощи формул:

Если гипотеза справедлива, то при неограниченном увеличении объемов выборок:

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

В этой связи можно предложит простую модификацию статистики :

у которой практически отсутствует последний недостаток.

Таким образом, вычислив статистику , P-значение вычисляем по формуле:

.

4. Критерий однородности Лемана-Розенблатта

Свободный от распределения критерий однородности Лемана-Розенблатта представляет собой критерий типа .

Имеем наблюдений и .

Предположим (точно так же как и в критерии Колмогорова-Смирнова), что:

Все наблюдений и взаимно независимы.

Все эти извлечены из одной непрерывной совокупности , распределенной согласно закону .

Все извлечены из одной непрерывной совокупности , распределенной согласно закону .

Рассматриваем гипотезу о том, что совокупности и идентичны, т.е. о том, что обе выборки извлечены из одной и той же совокупности. Её можно переписать так:



Будем тестировать эту гипотезу против альтернативы

.

Для проверки делаем:

Будем считать, что наши выборки упорядочены по возрастанию:

(если это не так просто упорядочиваем их).

Статистика критерия имеет вид:

Где

-- эмпирическая функция распределения, построенная по объединенному вариационному ряду двух выборок. Статистика используется в форме:

где -- порядковый номер (ранг) , -- порядковый номер (ранг) в объединенном вариационном ряду.

Статистика в пределе распределена как :

где и -- модифицированные функции Бесселя:

В отличие от статистики Колмогорова-Смирнова распределение статистики быстро сходится к предельному.

Таким образом, вычислив статистику , P значение вычисляем по формуле:

.

5. Статистическая регрессия

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

- выборка значений экспрессии для гена.

- временные точки, в которых замерена экспрессия

Т.е. значение замерено в момент времени

Строим полиномиальную регрессию:

Где

- сумма квадратов отклонений

Из требования , получаем систему для нахождения МНК оценок для b:

Также нам потребуется оценка для дисперсии коэффициентов:

Оценка для :

Проверяя гипотезу будем использовать статистику



в случае справедливости нулевой гипотезы, она имеет распределение Стьюдента с степенями свободы. Отсюда, формула для P-значения i-того коэффициента:

Где находится из (2), а

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

4. Гипергеометрический тест

Отдельно рассмотрим метод, специально разработанный для решения исходной задачи - обнаружения объектов (клонов и генов) с повышенным и пониженным уровнем экспрессии. Здесь в отличие от предыдущих методов анализ отдельных объектов проводится при учете всей совокупности данных в таблице.

Пусть:

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

При этом мы предполагаем, что для подавляющего количества генов в таблице нулевая гипотеза верна.

Параметр метода BV (Boundary Value) -- предназначен для разграничения между повышенной и пониженной экспрессией.

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

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

Для удобства, построим вариационный ряд:

Ясно, что в качестве не уменьшая общности можно рассматривать только значения из ряда .

Допустим, мы фиксировали , введем обозначение:

В случае верности нулевой гипотезы, отклонения от произошли по случайным причинам и в этом случае имеет гипергеометрическое распределение с параметрами , где:

- количество значений в таблице больше либо равных .

- общее количество данных в таблице.

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

Таким образом, вероятность ошибиться при отвержении гипотезы c помощью критерия может быть вычислена по формуле:

Таким образом, мы получили серию критериев с критическими статистиками где k пробегает все значения, для которых Логично из этих критериев выбрать критерий с наименьшей ошибкой первого рода:

Аналогичным образом вводим альтернативную гипотезу и обозначения:

- набор критических статистик.

- количество значений в таблице меньше либо равных .



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

И, окончательно

Причем, если то, отвергнув нулевую гипотезу, целесообразно объявить исследуемый объект повысившим свою регуляцию. Если же , то ген объявляется понизившим регуляцию.

Введем обозначения:

Тогда алгоритм можно описать следующим образом:

Упорядочиваем значения

Идем по массиву значений и находим величины



Итоговый ответ для каждого объекта - это скор, знак которого указывает на то, вероятность чего больше: того, что ген повысил или понизил экспрессию, а абсолютное значение, на то какова вероятность ошибиться при отвержении нулевой гипотезы. Скор близкий к нулю означает, что нулевая гипотеза, скорее всего, верна. Десятичный логарифм здесь взят для удобства. Т.е. например скор, равный указывает нам, что если мы отвергнем гипотезу о том, что объект не изменил свою экспрессию и примем альтернативу то мы ошибемся с вероятностью 0,01. Вероятность ошибки при принятии альтернативы

Важным моментом является то, что мы предполагаем (а обычно так и есть), что большая часть генов действительно не изменила своей экспрессии и отклоняется от некоторого среднего уровня только случайно или из-за погрешностей в измерении и поэтому правомерно использовать значения как некоторую “кучу” из которой можно брать случайные значения, характерные для экспрессии генов в данном эксперименте. Хотя количество данных в строке может быть каким угодно, но чем больше значений, тем достовернее результаты теста. При анализе данных, в которых для каждого объекта есть всего одно измерении экспрессии, тест вырождается в простое ранжирование объектов и выборе в качестве изменивших свою экспрессию объектов с большими значениями экспрессии.

В качестве BV можно выбирать, например 1 (когда данные уже представляют собой отношение паталогии к норме), уровень экспрессии в норме (для каждого объекта выбирается своё BV), усредненная экспрессия контрольной группы или средняя экспрессия самого гена (тогда исследуется отклонения гена от его среднего сзначения).

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

Мета-анализ

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

В качестве входных данных мы имеем несколько микрочиповых таблиц . Ключами в которых являются объекты, а значениями - замеры экспрессии для них.

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

На выходе для каждой имеем таблицу .

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

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

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

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

5. Анализ данных по раку молочной железы

Был проведен анализ данных микрочипового эксперимента, призванного выявить реакцию организма больного раком пациента на потенциальные лекарства. Данные представляют собой замеры уровня экспрессии генов пациентов больных раком молочной железы, после введения препаратов. В качестве потенциальных лекарств тестировались вещества RITA в объеме 0,1мM и 1мM, Nutlin в объеме 10мM. Замеры проводились в следующих временных точках (единицы измерения в часах):

Количество исследованных объектов - 22276.

Одной из задач анализа была оценка того, как изменится экспрессия генов в первые часы после приема препаратов. В частности очень важно проанализировать поведение генов, которые влияют на процесс апоптоза - программируемой клеточной смерти, которая наступает при нарушении правильного функционирования клетки (например, в результате мутации ДНК). В клетках раковых опухолей процесс апоптоза подавлен, что позволяет раковым клеткам продолжать делиться. Одним из самых важных противораковых генов является ген tp53 (tumor protein 53), который активируется при повреждениях генетического аппарата, а также при стимулах, которые могут привести к подобным повреждениям, или являются сигналом о неблагоприятном состоянии клетки (стрессовом состоянии). Его функция состоит в удалении из пула реплицирующихся клеток тех клеток, которые являются потенциально онкогенными (отсюда образное название белка р53 - хранитель генома). Сам процесс апоптоза, конечно, регулируется не одним геном, а целой генной сетью, однако решающим фактором является именно tp53, тогда как остальные гены заменяемы (существует несколько сигнальных путей, ведущих к tp53). Числовые значения в таблице варьировались от 0.1 до 70390. Перед началом анализа, данные подвергнулись следующей обработке: все значения экспрессии меньшие заданного порога, равного 200 были приравнены к этому порогу. Это было сделано ввиду стохастической природы данных: все значения меньше 200 представляют собой в основном случайный шум.

На первом этапе анализа данные были разбиты на три таблицы в соответствии с исследовавшимся препаратом: RITA_0.01, RITA_1 и Nutlin, далее каждая из этих таблиц анализировалась отдельно. В каждой временной точке кроме контроля (0ч.) было проведено только одно измерение экспрессии для каждого препарата, поэтому временные точки были сгруппированы следующим образом: и каждая из этих групп сравнивалась с контрольной группой в качестве которой выступали 3 замера уровня экспрессии в 0ч.



Рисунок 2. Эмпирические плотности распределения экспрессий генов в разные моменты времени. Вертикальной линией обозначен порог по которому отсекались несущественные значения

Статистический анализ проводился с помощью статистик Стьюдента, Вилкоксона и гипергеометрического теста. Граничным было выбрано значение 0,01. Т.е. отбирались только те гены, скор которых превышал 2, что является довольно сильным ограничением. Стоит отметить, что корреляции между скорами, полученными по одним и тем же данным различными методами во всех 27 случаях был выше 0,95, что позволяет говорить о надежности результата. Количество генов, отобранных различными методами для препаратов RITA и Nutlin приведены в таблицах 1 и 2 соответственно.

Таблица 1

Количество объектов изменивших свою экспрессию под воздействием препарата RITA в объеме 1uM

	2-4-6 часов	4-6-8 часов	6-8-10 часов
Статистика Стьюдента	2684	4297	4625
Статистика Вилкоксона	5340	6137	6293
Гипергеометрический тест	6933	7240	7220
Пересечение	2684	4296	4625



Таблица 2

Количество объектов изменивших свою экспрессию под воздействием препарата Nutlin в объеме 10uM

	2-4-6 часов	4-6-8 часов	6-8-10 часов
Статистика Стьюдента	981	2180	3040
Статистика Вилкоксона	3109	4450	5356
Гипергеометрический тест	5581	6355	7030
Пересечение	981	2180	3040

По этим таблицам видно, что хотя корреляции между результатами большие, но сами методы дают различные результаты. Это связано с тем, что методы Вилкоксона и гипергеометрический имеют дискретный характер, в отличие от метода Стьюдента. В то же время видно, что все гены отобранные методом Стьюдента также отобраны и другими методами. Однако, помимо пересечения с результатами анализа методом Стьюдента необходимо исследовать также результаты других анализов. Так, например, играющий важную роль в процессе апоптоза, онкоген JUN при гипергеометрическом анализе показал довольно большое значение скора: 5,033741591, в то время как при анализе методом Стьюдента скор не превысил порога в 2.

Итого было получено 27 таблиц скоров (для каждого из трех препаратов было проведено по три сравнения каждое тремя методами). Окончательно отбирались только те гены, скор которых превысил пороговое значение для всех методов анализа.

На втором этапе, после импорта полученных списков в базу данных Explain (http://www.biobase-international.com/), где объекты были соотнесены с генами, эти гены, с помощью имеющегося в Explain метода анализа, были привязаны к соответствующим сигнальным путям. Далее исследовались гены, имеющие отношение к генной сети, регулирующей апоптоз. В результате исследования было выявлено, что большая часть генов, составляющих сигнальный путь p53, значительно изменила свою экспрессию во всех девяти подгруппах. Количество генов из p53-сигнального пути, сильно изменившие свою экспрессию после приема RITA и Nutlin приведено в таблицах 1 и 2 соответственно.

Таблица 3

Количество генов p53-pathway, значительно изменивших свою экспрессию под воздействием препарата RITA

	Уровень экспрессии	2-4-6 часов	4-6-8 часов	6-8-10 часов
Активаторы апоптоза (антионкогены)	Повышен	3	10	11
	Понижен	10	12	14
Репрессоры апоптоза (онкогены)	Повышен	4	6	7
	Понижен	8	12	10

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

Таблица 4

Количество генов p53-pathway, значительно изменивших свою экспрессию под воздействием препарата Nutlin

	Уровень экспрессии	2-4-6 часов	4-6-8 часов	6-8-10 часов
Активаторы апоптоза (антионкогены)	Повышен	1	7	4
	Понижен	2	5	8
Репрессоры апоптоза (онкогены)	Повышен	2	2	2
	Понижен	4	3	7

Выводы

Стоит отметить, что сам ген tp53 даже при малых отклонениях от нормы сильно влияет на процесс апоптоза, поэтому при исследовании действия Nutlin он не проявил себя среди наиболее изменивших уровень экспрессии генов, в то же время RITA является значительно более мощным препаратом, уже на втором этапе tp53 появляется в результатах. Кроме того на третьем этапе стабильно сильнее всего активизируется ген MDM2, который связан обратной негативной связью с tp53: наработанный белок р53 стимулирует экспрессию гена MDM2, который в свою очередь подавляет активность tp53. Кроме того, активность белка р53 в отношении mdm2 и остальных генов-мишеней р53 задаёт тот временной промежуток, в течение которого р53 может выполнять свои функции. Таким образом, в результате действия как Nutlin так и RITA ген tp53 активизируется настолько, что включается природный механизм его подавления.

Из таблиц 3,4 видны некоторые особенности действия препаратов: RITA действует практически сразу и очень сильно, с большим количеством побочных эффектов. В то время как Nutlin действует постепенно и значительно слабее. Относительно результатов RITA в концентрации 0.01µM можно сказать, что действие сходно с 1µM, однако значительно слабее (действует даже слабее чем Nutlin). О более детальном действии препаратов на раковые опухоли говорить пока рано, требуется более детальное исследование полученных результатов совместно с исследованием реакции тканей на уровне белков. При анализе результатов были использованы следующие материалы [16,17], а также базы данных:

· http://www.genecards.org/.

· http://www.ihop-net.org/.

· http://www.ncbi.nlm.nih.gov/ (Entrez, Unigene).

· https://www.affymetrix.com/.

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

6. Моделирование генной экспрессии

Фиксируем некоторый интересующий нас объект . Для простоты будем называть его ген-мишень. Обозначим результат измерения уровня экспрессии гена-мишени в момент времени t через .

Обозначим набор генов, прямо или косвенно влияющих на ген-мишень: . Это - гены-регуляторы. Обозначим через - результат измерения уровня экспрессии i-того регулятора в момент времени t.

Как именно связан уровень экспрессии гена-мишени с уровнями экспрессии генов-регуляторов? В самом общем виде можно представить зависимость следующим образом:

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

здесь:

- фоновый уровень экспрессии, независящий от регуляторов.

- коэффициент деградации.

- уровень влияния i-того регулятора на ген-мишень.

- временная задержка нужная для завершения процесса трансляции транскрипционных факторов.

- стандартный винеровский процесс, представляющий собой ошибку измерения.

- квадратный корень из дисперсии ошибки.

Моделирует процесс “насыщения” гена мишени. Начиная с некоторого предела, наработка новых транскрипционных факторов уже практически не влияет на процесс транскрипции гена-мишени.

Задача построения модели регуляции гена-мишени состоит из нескольких этапов:

1. Найти набор генов регуляторов .

2. Выбрать вид модели, описывающей регуляцию гена-мишени.

3. Вычислить коэффициенты модели.

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

7. Линейная модель

Рассмотрим линейную дифференциальную модель регуляции гена-мишени одним геном-регулятором со сдвигом по времени.



Известными данными являются замеренные экспрессии гена-мишени и гена-регулятора в фиксированные моменты времени.

Необходимо найти параметры и оценить насколько выбранный регулятор объясняет поведение гена-мишени. Рассмотрим далее только первую модель, для второй все действия аналогичны. Исходя из биологических соображений, параметр должен принадлежать определенному заранее интервалу . В частности может быть , где T - период клеточного цикла. Выберем достаточно большое и разделим этот интервал на равные части. Обозначим:

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

1. Интерполирующий сплайн.

2. Сглаживающий “коридорный” сплайн.

3. Сглаживающий сплайн с постоянным параметром.

4. Ядерное сглаживание с помощью ядра Епанечникова.

5. Сглаживание методом наименьших квадратов.

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

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

Где - коэффициент детерминации. Чем больше тем лучше модель аппроксимирует , значения близкие к 1 говорят об отсутствии связи между X и Y. Таким образом, выбирается модель с максимальной величиной .

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

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

Алгоритм

1. Фиксируется ген-мишень.

2. Выбирается потенциальный ген-регулятор .

3. Полученные профили экспрессии сглаживаются одним из методов.

4. Выбирается величина сдвига Строится линейная регрессия:

Находятся параметры , коэффициент детерминации.

5. Проделываем шаги 4-5 для . Находим оптимальную задержку .

6. Повторяем шаги 2-5 для каждого потенциального гена-регулятора.

7. Выбираем ген-регулятор с наилучшим коэффициентом детерминации.

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

Возможны модификации метода

1. С биологической точки зрения целесообразно добавить коэффициент деградации

В этом случае считается множественная регрессия.

2. Можно рассматривать модель с несколькими регуляторами:

8. Методы сглаживания

Часто в работе с микрочиповыми таблицами имеет смысл преобразовывать дискретные данные в непрерывные. Это может быть связано с необходимостью вычислить значение в промежуточной точке, вычислить численную производную профиля экспрессии или в связи с тем, что некоторые значения уровня экспрессии могут быть пропущены из-за погрешности эксперимента. К сожалению данные полученные с помощью микрочипового эксперимента неизбежно несут некоторую иногда довольно значительную погрешность, так что точная интерполяция смысла практически не имеет. Рассмотрим несколько методов сглаживания данных, предложенных в литературе [7-10]. Всюду ниже, аналогично предыдущему, будем обозначать:

сетка.

дискретно заданная функция, которую требуется сгладить.



1. Сглаживающий “коридорный” сплайн

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

Сглаживающие сплайны возникают при решении задачи о минимизации функционала:

Где и -- заданные величины. Видно, что чем меньше коэффициенты , тем ближе проходит функция минимизирующая функционал к заданным значениям . Доказано, что среди всех функций , кубический сплайн , удовлетворяющий условиям , минимизирует функционал . [7] Исходя из этого требования сам сплайн строится итерационным способом.

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

2. Сглаживающий сплайн с постоянным весовым параметром

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

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

Основным вопросом при построении такого сплайна является выбор параметра . Можно попытаться выбирать значение , минимизируя некоторую оценку величины . В [9] показано, что в случае равномерной сетки с шагом значение дает минимум оценке величины . Однако, обратим внимание, что при таком значении достигается лишь минимум некоторой оценки этой величины, поэтому в конкретных случаях минимум погрешности приближения может достигаться при других значениях . Значение , может служить лишь некоторым ориентиром при выборе параметра сглаживания.

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

3. Метод наименьших квадратов

Введем непрерывную функцию для аппроксимации дискретной зависимости . Обозначим отклонения в узлах и . Метод построения аппроксимирующей функции из условия минимума величины Q называется методом наименьших квадратов (МНК).

Широко распространен выбор базисных функций в виде степеней :

.

Однако такой выбор базисных функций не является оптимальным с точки зрения решения системы нормальных уравнений с наименьшими погрешностями. При выборе такого базиса построение аппроксимирующего полинома выше 7 степени приводит к значительным погрешностям при решении системы нормальных уравнений.

Намного лучшие результаты может дать применение ортогональных полиномов (Чебышева, Лежандра, Лаггера, Якоби). Выбор таких полиномов в качестве базисных функций обеспечивает лучшую обусловленность системы нормальных уравнений и позволяет строить качественные аппроксимации больших степеней.

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

4. Ядерное сглаживание

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

Соответственно для производной (в случае дифференцируемости весовых функций):

Ядерный подход состоит в том, что весовая функция задается посредством функции плотности , которая регулирует размер и форму весов вблизи точки . Эту функцию принято называть ядром. Ядро -- это непрерывная ограниченная симметричная функция с единичным интегралом:

Вид весовых функций, использующий ядро был предложен в работах (Nadaraya, 1964) и (Watson, 1964):

где

- Ядерная оценка Розенблатта-Парзена для маргинальной плотности переменной x.

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

Согласно эти формулам окончательная оценка значения сглаживаемой функции и ее первой производной:



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

Чаще всего из-за некоторых оптимальных свойств используется ядро Епанечникова:

где

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

Для корректного вычисления значений первой производной сглаженной функции во всех точках интервала [ должно выполняться:

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

Выводы

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