Компьютерный анализ сплайсинга Неверов Алексей Дмитриевич

Содержание к диссертации

Введение

1. Введение и обзор литературы 6

1.1 Введение 6

1.2 Алгоритмы предсказания генов, основанные на статистических свойствах геномной последовательности 8

1.3 Альтернативный сплайсинг (АС) 34

2. Предсказание генов в геномах низших эукариот 61

3. Элементарные альтернативы в генах эукариот 80

3.1 Выявление элементарных альтернатив 82

3.2 Координированный альтернативный сплайсинг 92

4. Альтернативный сплайсинг и функции белков 99

4.1 Белковые изоформы альтернативно сплайсируемых генов 101

4.2 Частота ошибок сплайсинга 119

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

5.1 Динамическое программирование (ДП) 123

5.2 Построение базы данных альтернативного сплайсинга (EDAS) 130

5.3 Функциональные группы элементарных альтернатив 133

6. Основные результаты и выводы 134

Приложения 135

Благодарности 138

Список литературы.

• Алгоритмы предсказания генов, основанные на статистических свойствах геномной последовательности
• Предсказание генов в геномах низших эукариот
• Координированный альтернативный сплайсинг
• Частота ошибок сплайсинга

Введение к работе

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

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

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

Молекулярные механизмы сплайсинга.

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

Рисунок 1.1.1 Нуклеотидные последовательности участвующие в процессе

сплайсинга.

донорный сайт акцепторный сайт

сплайсинга сайт сплайсинга

Ul snfl

AGGU
экзон 1 интрон экзон 2

В процессе сплайсинга происходит распознавание трех участков пре-мРНК (Рисунок 1.1.1): 5' сайта сплайсинга (AGGURAGU донорный сайт), сайта ветвления (CTRAYY), и 3' сайта сплайсинга (YYYYYYYNCAGG акцепторный сайт).

Рисунок 1.1.2. Сплайсинг пре-мРНК.

В процессе распознавания сайтов

_*4л

U4/6 Ш

\

ftylAGGU

\

_4&

+

AS|GU Соединенные экзоны

сплайсинга происходит узнавание донорного сайта комплексом малых ядерных РНК сначала U1, затем U6/U4, а также узнавание сайта ветвления малой ядерной РНК U2. На следующем этапе происходит сближение 5' и 3' сайтов сплайсинга и образование комплекса U2, U4/U6. Далее происходит разрезание мРНК по донорному сайту сплайсинга и замыкание 5' конца интрона на 2' положение рибозы нуклеотида точки ветвления (Рисунок 1.1.2А). Малая ядерная РНК U5 сводит вместе донорный и акцепторный сайт сплайсинга. В результате последующей реакции происходит сшивка 5' и 3' сайта и полное вырезание интрона (Рисунок 1.1.2Б).

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

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

Алгоритмы поиска генов можно разделить на два вида - 1) статистические алгоритмы, источником информации для которых служит сама последовательность ДНК и 2) алгоритмы, использующие сходство с белками, EST, кДНК или последовательностями ДНК, содержащими гомологичные гены До недавнего времени такое разделение было абсолютным, однако в настоящее время граница стала сильно размытой. Как показано в [61] программы статистического распознавания генов имеют хорошую чувствительность, однако специфичность таких программ сильно зависит от длины межгенных интервалов -спейсеров. Ложные экзоны, предсказанные внутри длинных спейсеров, свойственных человеческому геному, заметно снижают специфичность. В значительной степени это связано с тем, что в настоящее время не существует хороших статистических моделей, позволяющих определить границы генов. Программы, основанные на сходстве, решают обратную задачу восстановления структуры гена по продукту сплайсинга того же самого гена (в простейшем случае) либо родственного гена, известного в другом организме. Эти алгоритмы имеют высокую специфичность, которая зависит только от степени близости продукта гена, структуру которого надо установить, и известного из базы данных гомолога. Серьезным недостатком такого подхода является то, что он не позволяет картировать ранее неизвестные гены. Специальным случаем является предсказание генов, основанное на том факте, что кодирующие белок экзоны значительно более консервативны по сравнению с некодирующей ДНК. Поэтому можно найти экзоны сравнением данного фрагмента с геномной ДНК, содержащей гомологичные гены. Предсказательная сила такого подхода для аннотации генов с известными ортологами в других организмах очевидно очень высокая. Экзон-интронная структура ортологичных генов млекопитающих часто в высокой степени консервативна. Например, среди ортологов человека и мыши 86% имеют одинаковое количество кодирующих экзонов, 46% имеют одинаковую длину белок-кодирующей области, причем только 1,5% имеют одинаковую длину белок-кодирующей области и разное количество экзонов [58]. Однако,

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

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

Подробный обзор алгоритмов, использующих сходство для предсказания генов, приведен в работе [49]. Такие алгоритмы позволяют идентифицировать 50-70% генов в новом геноме [47,48] и эта доля будет возрастать по мере секвенирования новых геномов. Тем не менее, применение статистических алгоритмов распознавания генов в процессе аннотации является обязательным этапом, особенно для организмов, имеющих большое значение для медицины и биотехнологии. Современные статистические программы распознавания генов эукариот и многие программы, предсказывающие гены прокариот, используют Скрытую Марковскую Модель (НММ) генома. Алгоритмы, основанные на НММ, превосходят в чувствительности и специфичности алгоритмы, использующие альтернативные подходы нейронные сети, линейный дискриминантный анализ, лингвистический анализ [47,50]. Для НММ существуют эффективные алгоритмы

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

Скрытая Марковская Модель эукариотического гена.

Успешное применение НММ к задаче распознавания генов человека в работе Burge и Karlin [51], предложивших программу GenScan, дало начало целой серии программ, использующих аналогичную модель геномной ДНК: Genie [52], FGENESH [53], GeneMark.hmm [54] и другие [47]. Перечисленные программы различались в деталях -статистических моделях сайтов, кодирующей и некодирующей ДНК и организмах, для которых было произведено обучение. В настоящее время все перечисленные программы обучены для предсказания генов широкого спектра организмов. Программа GeneMark.hmm, например, первоначально была разработана для бактериальных геномов, затем модель была расширена на геномы эукариот. При сравнении качества предсказаний различные программы демонстрируют близкие результаты [55].

Прежде чем обсуждать скрытые Марковские модели, дадим определение простой цепи Маркова [56]. Последовательность случайных величин х/, .. х,... называется цепью Маркова с состояниями S={1,2,. N}, Vi, х,є S, если выполняются условия: 1) для любого момента времени /, ^P(x_t = s) = 1 и 2) для любой последовательности отметок времени

l^<h^< .<;«<*'s,t є S и любых 5/,5^, В„ подмножествах S

выполняется P(Xj=t\ х,ієВі,х,2єВ2, ,x_l„eB_nx,=s)=P(Xj=t\x,=s) Свойство (2) называется свойством марковости и означает, что при любом фиксированном положении системы в данный момент времени, будущее поведение системы (/>/) не зависит от поведения системы в прошлом. Марковская цепь может быть представлена в виде набора (S,A,n) состояний S, матрицы вероятностей переходов между состояниями A=fa_st=P(x,=t\x_t i=s)} и вектором начальных состояний п=-( P(xi=s)), s=l N. Если вероятности переходов a_st не зависят от времени і, то цепь называется однородной. Вероятность любой последовательности Х= х;, х,.. , хт определяется выражением

Р(Л = Лх„...х₁...,х_т) = ^Па_ЯііЛ

1=2

Будем называть Марковской цепью порядка к>1 последовательность случайных величин хо, xi, х, ...,хт (Т>к), обладающих свойством: для любого момента времени i>k вероятность состояния зависит от последовательности к предшествующих состояний Р(х,\х,і, ,х,ь хо)= Р(х,\х,і,...,х,.і(). Вектор начальных состояний определяется для последовательностей хо, xi, Xki Легко заметить, что Марковская цепь порядка к сводится к Марковской цепи первого порядка путем определения новых состояний как Sr .

Дадим теперь формальное определение Скрытой марковской модели [57]. Теперь мы будем различать последовательность символов Х= Xi, . х, , хт, х,єГ (алфавит) и последовательность состояний Z=zi,...,zt, z,e S. Пусть задана обычная марковская цепь состояний (S,A,n), для каждого состояния зададим распределение эмиссионных вероятностей символов e_s(b)=P(x,=b\ z, =s), х,єГ, г,є S. Глядя на последовательность символов X, мы более не можем однозначно сказать какой последовательностью состояний Z, она была сгенерирована - состояния "скрыты" от наблюдателя. Мы можем записать совместную вероятность последовательностей символов и состояний:

1=2

Дальнейшим обобщением модели будет НММ с длительностями состояний (GHMM) [51,54]. Мы предполагали ранее, что в каждом состоянии модели генерируется единственный символ, теперь мы будем считать, что в каждом состоянии s может быть сгенерировано п символов с вероятностью d_s(n), В общем случае длина наблюдаемой последовательности символов Т не равна длине последовательности состояний L.

Обозначим D=ni,ri2, пц Vw =Г - последовательность длин в каждом состоянии.

/-і

Обозначим /, = ^ п_к длину последовательности символов, сгенерированных

к-\

последовательностью первых / состояний. Так же как для обычной НММ мы можем
записать совместную вероятность наблюдаемой последовательности X,
последовательности состояний Z с длинами D: (1)

^^2^) = ^^^^1,^)^^^)173^,^^^/^+1,/,-1]^^).

1=2

Здесь мы обозначили фрагмент последовательности символов, сгенерированный в состоянии і - X[l,.i+1,1,-1]. Мы будем называть простыми состояниями состояния Марковской цепи, в которых генерируется единственный символ и обобщенными состояниями - те состояния, в которых генерируется несколько символов. Скрытая Марковская модель может содержать как простые, так и обобщенные состояния.

Состояния цепи в задаче распознавания генов соответствуют функциональным элементам гена или геномной ДНК: сигналам, регулирующим экспрессию гена, трансляцию и сплайсинг; кодирующим и некодирующим фрагментам генома. В алгоритме Genscan [51] структура генома моделируется с помощью 27 состояний НММ. Гены на прямой и обратной цепи ДНК моделируются с помощью 13 состояний для каждой цепи соответственно, и одно состояние моделирует межгенный интервал - спейсер. Диаграмма состояний модели приведена на рисунке 1.2.1.

Рисунок 1.2.1. Структура НММ эукариотического генома, используемая в программе Genscan [51]. На рисунке используются обозначения: кодирующие состояния (экзоны) Emu - начальный, Е* - последний экзоны, Е„ 1=0,1,2 внутренние экзоны в соответствующей рамке считывания; E_s,_ngie - одноэкзонный ген; некодирующие состояния

I„ i=0,1,2 интроны, N -межгенная область (спейсер), Р- промотор, А - сайт полиаденилирования, F и Т -5' и З'-нетранслируемая область

Большинство переходов
между состояниями имеют
тривиальные вероятности: 0
(состояния не связаны на
диаграмме - переходы
запрещены) или 1 (состояние
имеет единственный вариант
перехода в новое состояние).
В большинстве моделей
генома игнорируется

существование нетранслируемых экзонов, за исключением модели гена в алгоритме Doublescan [58], где такие экзоны могут быть предсказаны за счет сильного сигнала в сайтах сплайсинга. Здесь необходимо отметить, что в отличие от

рассматриваемых алгоритмов, модель Doublescan описывает выравнивание двух

последовательностей ДНК, содержащих гомологичные гены. Одноэкзонные гены, а так же первый и последний экзоны гена, содержащего интроны, моделируются отдельными состояниями НММ. Первый экзон гена начинается всегда со СТАРТ кодона (ATG), последний экзон всегда заканчивается одним из СТОП кодонов (ТАА, TAG, TGA). Интроны рассматриваются как вставки в кодирующую последовательность мРНК. Они могут приходиться между соседними кодонами, а так же после первого либо второго нуклеотида внутри кодона. В зависимости от числа нуклеотидов разорванного кодона, которые переносятся на следующий экзон на каждой цепи ДНК, существует по три интронных состояния модели. Іо, І/, І2. Внутренний экзон может транслироваться в одной из трех возможных рамок 1-0,1,2 в зависимости от того, сколько нуклеотидов от начала экзона относятся к кодону мРНК, разорванному предшествующим интроном. Таким образом, для каждого интронного состояния возможен единственный переход в соответствующее состояние внутреннего экзона - /, -> Е, , i=0,l,2. Обратная цепь ДНК моделируется состояниями соответствующими один к одному состояниям прямой цепи, а переходы осуществляются в обратном 3' 5' направлении.

Эмиссионные вероятности и распределения длин.

Мы определили состояния модели и матрицу переходов, описывающих геномную ДНК, теперь кратко рассмотрим статистические модели, лежащие в основе расчета эмиссионных вероятностей. По оценкам [59], для задач различения кодирующей и некодирующей ДНК наилучшими моделями являются однородная (не зависящая от позиции) Марковская цепь пятого порядка для некодирующей ДНК и трех-периодическая Марковская цепь пятого порядка для кодирующей. Трех-периодическая Марковская цепь является неоднородной Марковской цепью, где вероятности зависят от позиции генерируемого нуклеотида в кодоне (0,1,2). Для обучения этих моделей требуется оценить достаточно много параметров: 4⁶ для некодирующей модели и 3*4⁶ для кодирующей. Для многих геномов можно использовать модели с меньшим числом

параметров, при сохранении качества предсказания. Для кодирующей модели часто используются статистика кодонов, Марковские цепи порядка к=2 и 3, для некодирующих состояний эффективны Марковские модели первого порядка и частоты нуклеотидов [64]. Если геном имеет сильно неравномерное распределение GC состава, как, например, геном человека, где существуют протяженные области - изохоры, то часто используются неоднородные Марковские модели, зависящие от GC-состава В программе Genscan выделяются четыре группы с диапазоном GC: <43, 43-51, 51-57, >57% . Наблюдаются существенные различия статистических свойств организации генома между группами [51]: средних длин спейсеров, интронов, транскрипта (пре-мРНК), кодирующей части мРНК, доли одноэкзонных генов. В программе Genscan от группы GC состава, к которой относится аннотируемая последовательность, зависят начальные вероятности состояний модели, вероятности переходов и эмиссионные вероятности кодирующих и некодирующих состояний.

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

Сайт ветвления человека обладает очень слабым консенсусом - только 30% интронов имеют сайт ветвления, удовлетворяющий максимально вырожденному консенсусу YYRAY в области [-40, -21] от акцепторного сайта. В программе Genscan для

моделирования области интрона, содержащей сайт ветвления, используется неоднородная Марковская цепь второго порядка. Важной особенностью модели является способ обучения - условные частоты нуклеотидов в позиции і вычисляются как средние условных частот в пяти соседних позициях: i-2, i-l, і, i+l, i+2. Авторы отмечают, что такая модель имеет слабую, но выраженную тенденцию генерировать в избытке YYY триплеты, а так же триплеты, характерные для сайтов ветвления: TGA, ТАА, GAC, и ААС.

Эмиссионные вероятности нуклеотидов донорного сайта могут быть получены из позиционной матрицы частот. Однако донорный сайт млекопитающих более эффективно распознается моделями, учитывающими зависимости между нуклеотидами в различных позициях. Статистически значимые зависимости в донорном сайте характерны как для соседних, так и удаленных друг от друга позиций сайта [51], которые учитываются в модели максимальных декомпозиций MDD (Maximum Dependence Decomposition). Отличительным свойством используемой модели является то, что безусловные позиционные вероятности заменяются на условные вероятности, зависящие от нуклеотида в позиции, обладающей максимальным совокупным влиянием на остальные позиции сайта.

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

требует специальных алгоритмов кластеризации, которые пока не интегрированы в программы распознавания генов

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

Алгоритмы декодирования для НММ.

Алгоритмы для НММ делятся на две группы: алгоритмы декодирования и алгоритмы автоматического обучения параметров модели [57,51,54]. Мы будем описывать алгоритмы для обобщенной НММ с длинами состояний, ровно те же самые алгоритмы применимы для обычной НММ. Алгоритмы декодирования, так же называемые алгоритмами разметки, ставят в соответствие наблюдаемой последовательности символов последовательность скрытых состояний. Существует три класса алгоритмов декодирования: алгоритм Витерби, алгоритм апостериорного декодирования (forward-backward algorithm) и случайный выбор последовательности состояний (НММ sampling).

Алгоритм Витерби.

Алгоритм Витерби является алгоритмом поиска оптимальной последовательности скрытых состояний и их длин, оптимизирующий совместную вероятность (1 см. Скрытая Марковская модель эукариотического гена). Пусть v(z_k,n_k,l) - максимальная

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

последовательности наблюдаемых Х[1,1] : / = ^и, . Пусть так же вероятностьu(z_k,n_k,l)

известна для всех состояний z*eS и для всех щ

ф*₊1>л*₊1>^/+и*₊1)^:

^u(h,un_k+i,l + n_k+]) = P(X[l,l + n_k+l])d (n_k+l) max (v(z_k,n_k,l)a ).

Запишем теперь рекурсивный алгоритм, вычисляющий оптимальное значение совместной вероятности наблюдаемых X и разметки Z*,D*:

1. инициализация (1=0): u(z_],0,0) = л_г, z/eS;

2. рекурсия 1=Т 1,п_к<1:

v(z_k,n_k,l) = P(X[l-n_k +l,l])d_it(n_k)max_{h i65}^ ^Mz^n^J-n^a^) ptr{z*_k,n*_k,l) = argmax_2j ^ ₍(u(z,_,,n_k__x,l-n_k)a₂^_Zt)

3. завершение: P(X, Z*, D*) = max_{2j ч} (v(z_k, n_k, L)); ptr (z *_k, n *_k, L) = arg max _{ч ч} (v(z_k, n_k, L));

4. обратный ход (1=T l):(z*_ki,n_k.i,l-n^)=ptr(z*_k,n_k,l)

позволяет восстановить оптимальные последовательности состояний Z* и длин D*.

Алгоритм апостериорного декодирования (forward-backward algorithm).

Алгоритм апостериорного декодирования для любого фрагмента последовательности наблюдаемых символов Х[1,1+п] позволяет вычислить вероятность того, что он был сгенерирован в некотором состоянии zeS Применительно к задаче распознавания генов, можно вычислить вероятности того, что некоторый фрагмент ДНК, является экзоном/интроном/.../. Для этого нам надо уметь вычислять величину Р(Х)= ^P(X,Z,D) - вероятность последовательности в рамках принятой НММ. Для

(ZD)

того чтобы вычислить величину Р(Х), введем сначала «прямые» переменные

a(l,z) = P(X[\,l],z), которые являются вероятностью того, что последовательность

наблюдаемых до координаты /, сгенерирована так, что последний символ сгенерирован в

состоянии zeS. Запишем рекурсивный алгоритм, подобный алгоритму Витерби, для

вычисления «прямых» переменных:

1. инициализация (/=7): а(\,д) = л_?, ceS,

2. рекурсия(1=Т. 1). a(/,z) = ^a_g2a(l-n,g)P(X[l-n,!-\]\g)d_g(n-\)

я=1 geS g*z

3) завершение: Р(Х) = ^а(Т,д)

Запишем теперь выражение для вероятности сгенерировать последовательность X, так чтобы символы в интервале а<Ь были сгенерированы в состоянии г.

P(X,X[a,b],z) = P(X[l,a-l],z)P(X[a,b] | z)d₂{b - a + \)P(X[b + \,Т] | X[\.b],z) = Р{Х[\,а -1],z)d_z(b-a + \)Р(Х[а,Ь] | z)P(X[b +1],z) = а(а,z)d_z(b-a + \)P(X[a,b] \ z)p(b,z) Вторая строка выражения является следствием условия марковости - все, что

сгенерировано после момента Ъ зависит только от состояния цепи в этот момент Так же

было введено обозначение для «обратных» переменных /3(l,z). Алгоритм вычисления

«обратных» переменных полностью аналогичен алгоритму вычисления «прямых»

переменных:

3. инициализация(1=Т): /?(r,^) = l,?eS;

4. рекурсия (№1,1): /3(/,2) = а_г/(Х[/Н/ + л]кК(")/?('+ ",<Г);

я-1 geS $*i

3) завершение: Р(Х) = ^^(1,д)тг_?.

Если вычислены «прямые» и «обратные» переменные, то мы легко можем рассчитать апостериорную вероятность состояния для фрагмента наблюдаемых Х[а,Ь]:

a(a,z)d_z(b-a + l)P(X[a,b]\z)№z)
r(z,a,b\X) = .

Р(Х) Апостериорное декодирование имеет приблизительно такую же предсказательную силу, как и алгоритм Витерби, однако существует очень полезное применение этого алгоритма, когда нужно приписать меру надежности предсказания экзона к фрагменту ДНК [51].

Случайный выбор последовательности состояний (НММ sampling).

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

Лемма 1.

Пусть задан ориентированный ациклический граф с источником s и стоком t и веса на ребрах. Будем обозначать вес ребра w(e), или w(u,Uj), вес вершины - w(u_t). Предположим, что выполняется условие нормировки на веса путей

Ир)
Z ^11,)1^1/,.,,11,)4-(1/,) = 1

Ур-(щи₂, н_Нр)) 1-1

Где сумма берется по всем путям из источника в сток. Можно случайным образом за время 0(к(р)) выбрать путь/? = (щ.Щ, . ,щ(р)) с вероятностью

Нр) Рг(^) = ф)\[ Ф,._х, и, )ф,)

1=2

Действительно, пусть на каждой вершине заданы «обратные» переменные {t(u,) -являющиеся суммой весов всех путей из вершины и, в сток. Пусть утверждение верно для

максимального пути из источника в сток длиной п ребер. Если n = 1 утверждение

очевидно. Докажем для случая, когда длина максимального пути равна п+1 ребер. Рассмотрим вершины и/, .щ связанные с источником 5 соответствующими ребрами в], ,ек Выберем ребро і с вероятностью w(s)w(e;)/?,. Обозначим а вес пути из вершины и, в сток t. Расстояние до стока из выбранной вершины не больше п ребер. Выберем путь в сток из вершины и, с вероятностью а//?, Вес пути из источника в сток через выбранную вершину равен w(s)aw(e,), и этот путь был выбран с вероятностью w(s)w(e,)P, cr/fl, = w(s)aw(ej.

Лемма определяет алгоритм случайного выбора пути на графе:

Вычисляем «обратные» переменные р(и^ для каждой вершины.

Затем, начиная от источника, выбираем ребро в следующую вершину и, с вероятностью - w(u,.i)w(e^)P(u), как в доказательстве.

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

Оценка качества предсказания.

Качество предсказания оценивается на основании множества последовательностей ДНК с известной аннотацией - тестового множества. Здесь мы не будем обсуждать подробно проблемы построения тестового множества, отметим, что чаще всего в него включаются фрагменты ДНК, содержащие гены, подтвержденные экспериментальной процедурой или маркером экспрессии (мРНК, EST, белок). Существует несколько мер того, насколько две разметки (предсказание и реальная структура) тестовой последовательности близки друг к другу. Для предсказания генов существует несколько уровней оценивания качества предсказания: нуклеотидный, экзонный и уровень гена. Качество предсказания на нуклеотидном уровне, является мерой того, насколько разметка нуклеотидов на кодирующие РР (predicted positive) и некодирующие PN (predicted

negative) близка к аннотации: АР (actual positive) и AN (actual negative). Нуклеотид может быть предсказан в соответствии с аннотацией, тогда он учитывается как верно предсказанный кодирующий TP=PPnAP (true positive) или некодирующий TN =PNnAN (true negative). Возможны два типа ошибок - нуклеотид неверно предсказан в качестве кодирующего FP (false positive, перепредсказание) или некодирующего FN (false negative). Предсказание кодирующей части гена характеризуется специфичностью Sp⁺=(TP)/(PP) и чувствительностью Sn⁺=(TP)/(AP), обе меры оценивают относительный вклад соответствующих типов ошибок (FP и FN). Величина в скобках обозначает число нуклеотидов в соответствующем множестве. Существуют несколько широко используемых интегральных мер, учитывающих относительные вклады обоих типов ошибок:

(TP)(TN) - (FP)(FN)

- коэффициент корреляции СС = і = (correlation coefficient),

J(PP)(PN)(AP)(AN)

приближенный коэффициент корреляции А С = 1 / 2[Sn⁺ + Sp* +Sn~ + Sp' ] -1 (approximate correlation),

- среднее между чувствительностью и специфичностью (Sn⁺ + Sp⁺)/2 и пересечение

_QQ=(PPnAP) (PPkjAP)

Чувствительность и специфичность по отношению к некодирующим нуклеотидам обозначены как Sn" и Sp' соответственно. Обычно для каждой последовательности из тестового множества мера качества предсказания вычисляется отдельно, затем общая характеристика вычисляется как среднее по всем последовательностям. Если тестовое множество включает в себя последовательности ДНК, подавляющая часть нуклеотидов в которых являются кодирующими, то коэффициенты корреляции СС и АС для таких последовательностей получаются очень низкими либо их невозможно вычислить. Это бывает в том случае, когда интроны в генах в этих последовательностях много короче экзонов, либо гены не содержат интронов, и межгенные участки не включены в тестовое

множество. Для таких последовательностей обычно пренебрегают предсказанием некодирующих нуклеотидов и вместо коэффициентов корреляции СС и АС используют среднее (Sn⁺ + Sp⁺)/2 и пересечение QQ. Очевидно, что альтернативным решением проблемы плохой аннотации некодирующих нуклеотидов является вычисление мер качества предсказания для всей совокупности последовательностей, вместо вычисления для каждой в отдельности с последующим усреднением. Проблема тестирования алгоритмов предсказания генов на длинных последовательностях ДНК, содержащих несколько генов, рассмотрена в работе Guigo [61] где было предложено формировать тестовое множество из псевдо геномных последовательностей, содержащих несколько генов на обеих цепях ДНК, разделенных и искусственно сгенерированными межгенными интервалами в соответствии с распределением длин, характерным для генома. Современные статистические алгоритмы предсказания генов обладают высокой нуклеотидной чувствительностью и специфичностью - более 90% [51,55,62].

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

специфичностью. Эти две меры вычисляются аналогичным образом - также как на уровне нуклеотидов, с той разницей, что учитываются предсказания экзонов или генов. На экзонном уровне используются полезные меры перепредсказания - доля предсказанных экзонов, которые не пересекаются в соответствующей рамке ни с одним из экзонов аннотации WE (wrong exons), а так же недопредсказания - доля экзонов в аннотации, не пересекающихся ни с одним из предсказанных экзонов ME (missed exons). Для алгоритмов предсказания генов чувствительность и специфичность на экзонном уровне составляет 80-

90% для разных геномов, доля пропущенных экзонов ME составляет около 10%, доля ложных экзонов WE около 5% [51,55,61,62] Здесь следует заметить, что игнорируется существование альтернативного сплайсинга Чувствительность на геномном уровне составляет около 40-50% [51,55]

Обучение НММ.

Скрытая Марковская модель содержит большое число параметров - вероятности начальных состояний и переходов, параметры моделей эмиссионных вероятностей. Качество предсказания генов критично зависит от того, насколько хорошо оценены (обучены) эти вероятности. Параметры эмиссионных моделей в значительно большей степени влияют на результат, чем начальные вероятности и вероятности переходов [63]. Если существует обучающее множество последовательностей с известной разметкой (аннотацией), то достаточно просто оценить параметры модели. Вероятности эмиссии могут быть оценены на основании статистики встречаемости олигонуклеотидов внутри соответствующих состояний модели. Оценивание вероятностей переходов осуществляется на основании статистики встречаемости самих состояний модели в обучающем множестве. Например, вероятности переходов в Genscan оцениваются на основании статистических свойств генов, включенных в обучающее множество: среднего числа интронов в генах, содержащих более одного экзона, и доли одноэкзонных генов Распределения длин так же извлекаются из обучающего множества, для состояний с геометрическим распределением достаточно оценить среднюю длину состояния, для экзонных состояний строятся гистограммы, которые затем сглаживаются и усредняются с периодом 3 нуклеотида, для того чтобы убрать трех периодический компонент в распределении. Обучающее множество обычно содержит слишком мало одноэкзонных генов из-за их относительно небольшой частоты в геноме. В качестве приближенной оценки средней длины одноэкзонного гена используется средняя длина кодирующей последовательности (кодирующей части мРНК, CDS).

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

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

Всегда ли сложные модели лучше, чем простые? Сколько параметров должно быть в модели для достижения хорошего качества различения кодирующей и некодирующей ДНК? Для бактериальных геномов можно построить эвристическую модель, число параметров в которой равно трем, позволяющую предсказывать более 93% генов. Для сравнения, трех периодическая Марковская модель пятого порядка (3*4⁶=12288 параметров) позволяет увеличить долю предсказанных генов менее 1% процента. В статье [64] авторы, проанализировав 17 геномов бактерий различного GC состава от 28,6%

(Borrelia burgdorferi) до 65,6% (Mycobacterium tuberculosis), показали, что существует выраженная зависимость позиционной частоты нуклеотида в кодоне от GC состава генома. Частоты 10 из 20 аминокислот так же существенно зависят от GC состава. Четыре из этих 10 аминокислот кодируются SSN (S = {C,G}; N = {A,C,G,T}) типом кодонов: аланин (А), глицин (G), пролин (Р) и аргинин (R). Для аргинина помимо четырех кодонов SSN типа существует еще два кодона: AGA и AGG. Частоты встречаемости этих аминокислот в белках увеличиваются с увеличением GC состава генома Пять других аминокислот имеют WWN (W={A,T}) тип кодонов: фенилаланин (F), изолейцин (I), лизин (К), аспарагин (N), тирозин (Y). Частоты WWN аминокислот уменьшаются с увеличением GC. Метионин - редкая аминокислота в бактериальных геномах (1,8% состава белка), так же является аминокислотой WWN типа, но ее частота не зависит от GC состава. Единственная аминокислота валин (V) из всех аминокислот, кодоны которых содержат в первых двух позициях смесь S и W нуклеотидов, обладает существенной зависимостью частоты от GC состава генома. Зависимости частоты каждого из четырех типов нуклеотидов от GC% в позициях 0, 1 и 2 кодона были приближены линейными функциями. Частоты каждой из 10 аминокислот, зависимых от GC состава, так же были приближены линейной регрессией. Частоты аминокислот, для которых не наблюдалось значимых зависимостей от GC состава, были зафиксированы на уровне, характерном для генома Е coh. Частоты кодонов могут быть восстановлены по известному GC% составу в соответствии с выражением:

/,(w,) = /,(gc%) vlT'*²⁾ л'

/ ,JІ\^ХП^Х\^Х7)

^х<}*\^х1 ~*A

где щп]П2 обозначает один из кодонов аминокислоты А, /а - частота аминокислоты А (зависит от GC%), f\ - частоты кодонов, полученные на основании позиционных частот соответствующих нуклеотидов (зависит от GC%),fR - частота кодона, скорректированная на частоту аминокислоты А и частоты синонимичных кодонов //. Существенными

достоинствами эвристической модели являются: малое число параметров (частоты трех типов нуклеотидов), а, следовательно, малый размер последовательности ДНК, необходимой для обучения (>400 нук.), применимость к геномам с неоднородным распределением GC состава. Малый объем обучающего множества позволил с успехом применять эту модель для распознавания генов в геномах вирусов ВИЧ-1 и Т-клеточного лимфотропного вируса тип I.

Применима ли эвристическая модель для геномов эукариот? Геномы бактерий устроены много проще, чем эукариот. Гены бактерий не имеют интронов, и плотность кодирующей ДНК много больше, чем в геномах эукариот. За счет большой доли некодирующей ДНК влияние GC состава на статистику кодонов может быть выражено слабее, чем в геномах бактерий. Тем не менее, эвристическая модель, использованная для вычисления эмиссионных вероятностей кодирующих состояний, позволила предсказать 88,5% экзонов и 65% генов из тестирующего множества С reinhardtn. Остальные состояния Скрытой марковской модели эукариотического генома были обучены на основании множества аннотированных генов С reinhardtn.

Автоматическое обучение НММ.

Рассмотрим формально задачу выбора значений параметров для скрытой марковской модели. Обозначим набор параметров НММ - в Параметрами марковской цепи являются вероятности переходов (а/сі), вероятности эмиссии ek(b)=P(x=b\z=k) и распределение длин обобщенных состояний. Будем предполагать, что имеется множество из п последовательностей ДНК у=(у, „у"), сгенерированных независимо одной марковской цепью, причем соответствующие последовательности состояний неизвестны (скрыты) Будем максимизировать log правдоподобия по всем возможным значениям параметров в

/-і

Таким образом, задача обучения НММ формулируется как задача максимального правдоподобия (МП) модели. Существует три основных алгоритма решения МП задачи для НММ* алгоритм Баума-Велтча [57], обучение на основе алгоритма Витерби (Viterbi training) [57,62] и выборочный метод Гиббса (Gibbs sampling) [65]. Два первых алгоритма относятся к алгоритмам максимизации математического ожидаемого (expectation maximization), гарантируют нахождение локального максимума целевой функции (функционала). Если в пространстве параметров существует несколько локальных максимумов, то будет существовать зависимость от начальных значений параметров ё⁰⁾. Поэтому для сложной модели следует прилагать специальные усилия, чтобы привести алгоритм в окрестность глобального максимума. Метод Гиббса позволяет находить глобальный максимум за конечное число итераций. Все перечисленные алгоритмы имеют интересные приложения к задаче распознавания генов.

Для объяснения алгоритмов введем сначала дополнительные обозначения. Пусть вектор z обозначает множество оптимальных последовательностей скрытых состояний для совокупности последовательностей наблюдаемых у, h(y^) - число вхождений свободной переменой і єв в подмножество данных у„ определяемое вектором Z. Свободной переменной будем называть параметр модели, независящий от остальных параметров. Подмножество данных у, (фрагментов ДНК) - это подмножество исходных данных, на котором можно подсчитать число вхождений свободной переменной, например, набор экзонов для подсчета вхождений олигонуклеотидов. Рассмотрим два частных случая свободных переменных: вероятности перехода между состояниями аы и эмиссионные вероятности олигонуклеотидов еф). Соответствующие значения h(y,) обозначим Аи и Еф), эти обозначения будут использованы в описании алгоритма Баума-Велтча.

Если известно количество вхождений для соответствующих параметров Аи и Ек(Ь), то можно записать оценки максимального правдоподобия этих параметров. Оценки максимального правдоподобия (МП) вероятностей переходов и эмиссии определяются

выражениями: а_и =^- и е_к = ^к .

На каждой итерации алгоритм Баума-Велтча вычисляет математическое ожидание числа вхождений свободной переменной. На следующей итерации новое значение свободной переменной приравнивается оценке максимального правдоподобия. Алгоритм Баума-Велтча имеет вероятностное обоснование. Запишем вероятность того, что в последовательности наблюдаемых X в позиции b происходит переход из состояния z, в

состояние z_J+i P(X[a,b],Zj =k,z_J+l =l\X,6) = * —— = ,

Г (Л)

где «прямые» а и «обратные» (5 переменные определены в разделе "Алгоритм
апостериорного декодирования". Рассматривая совокупность последовательностей у,
получим ожидаемое число событий перехода к->1:

J Г\У )y>[aJ>W

где верхний индекс «прямых» и «обратных» переменных обозначает номер
последовательности в совокупности наблюдаемых у, выражение у'^Ъ] обозначает
фрагмент последовательности j. Аналогично вычисляется ожидаемое число эмиссий
символа b в состоянии к:

j *\У )yJ[a,b]sy^Jie[a,b}y^J[i]=b

где внутренняя сумма берется только по тем позициям /, в которых генерируется символ - Ъ. Ожидаемое распределение длин в состоянии к так же является параметром модели и может быть вычислено, как сумма вероятностей сгенерировать п=Ь-а+1 символов для п=1,2 Так как в данном итерационном процессе мы движемся к максимуму в

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

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

Алгоритмы Баума-Велтча и обучение Витерби использовались для оценивания параметров в программах GeneMarkS [66] и EasyGene [67] для прокариотической модели гена. Из-за сложной структуры моделей эукариотической ДНК долгое время считалось затруднительным использование автоматического обучения. Как уже было отмечено, алгоритмы Баума-Велтча и алгоритмы обучения по Витерби не гарантируют нахождение оптимальных значений параметров модели. В работе [62] авторы применили обучение Витерби в программе GeneMark.hmm ES-3. По результатам тестирования автоматическое обучение превосходит традиционное обучение. Для того чтобы алгоритм обладал высокой чувствительностью и специфичностью, авторы применили ряд эвристик. На начальных итерациях обучаются только эмиссионные вероятности экзонов и интронов, после того, как алгоритм начнет достаточно хорошо различать кодирующую и некодирующую ДНК, обучаются вероятности эмиссии сайтов, вероятности переходов и распределения длин. Точка в пространстве параметров, к которой сходится итерационный процесс, не зависит от модели начальных эмиссионных вероятностей. Наиболее быстрая сходимость достигается при использовании эвристической модели для кодирующих нуклеотидов [64], которая была описана в предыдущем разделе. Благодаря тому, что при таком способе обучения специфичность предсказания увеличивается быстрее, чем

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

Алгоритмы предсказания генов долгое время четко подразделялись на статистические алгоритмы и алгоритмы, использующие информацию о сходстве с другими последовательностями. В настоящее время разработаны новые алгоритмы, где на основании формализма скрытых Марковских моделей объединяются различные источники информации. Интересный алгоритм предсказания ортологичных генов в совокупности геномных последовательностей различных организмов предложен в работе Chatterji и Pachter [65]. В основе алгоритма лежит следующая идея - гомологичные гены в каждой последовательности сгенерированы одной скрытой Марковской моделью, параметры которой оцениваются методом Гиббса.

В случае скрытой Марковской модели отсутствует информация о последовательностях состояний, поэтому необходимо выбирать параметры 9 вместе с последовательностями состояний из апостериорного распределения p(z, 9\у). Однако было бы удобно разделить выбор z и в так, чтобы итеративно выбирать разметку каждой последовательности ДНК из условного распределения p(z/ \t}'^jI, в,у) (1 ^ < п) (1), а параметры из апостериорного распределения р(9\г,у) (2), где t!'^jI - обозначает вектор с исключенной компонентой ]. Мы уже рассматривали эффективный алгоритм случайного выбора последовательности скрытых состояний (разметки) НММ [44] из распределения:

(3)p(zl,..^\yl,...y_T,^\y^[-^>])

где Т и L обозначают длины соответствующих наблюдаемых последовательностей и последовательностей состояний. Формулу (3) можно получить, проинтегрировав формулу (1) по априорному распределению параметров р(в), тем самым выбор последовательностей скрытых состояний не будет зависеть от первоначального выбора параметров НММ. Можно показать, что распределение (3) пропорционально правдоподобию разметки і для последовательностиУ"

p(z ^J = к I y^J,z^[-^jl , y^["^j]) oc Yl h{y\-^j])^Hyl)

Другими словами, алгоритм состоит из следующих шагов* (А) последовательность j изымается из совокупности у, (Б) пересчитываются вхождения свободных переменных V і ев - h(y^Jy, (В) вычисляются МП оценки параметров модели V / ев на основании h(/^j},) (Г) для исключенной последовательности случайно выбирается разметка z¹, затем последовательность у возвращается в совокупность у Шаги А-Г повторяются рекурсивно, пока не будет выполнен критерий сходимости Напомним, что /_; обозначает совокупность фрагментов последовательности j, зависящих от разметки z", на которой вычисляется число вхождений /z(yV свободного параметра і ев Свободными параметрами (переменными) являются независимые параметры модели. Предсказание генов.

Будем предполагать, что последовательности ДНК, для которых мы хотим получить экзон-интронную структуру, связаны друг с другом через эволюционное дерево в виде звезды, где все ветви, содержащие листья, встречаются в одной точке. Это серьезное ограничение, однако, справедливо для равно удаленных друг от друга последовательностей млекопитающих. По методу Гиббса обучались частоты олигонуклеотидов в экзоне и распределение длин экзонов. Программа предсказания генов SLAM [101] модифицировалась так, чтобы предсказывать гены только в одной последовательности ДНК, наподобие GENSCAN [51]. Отличие состоит в том, что для экзонов используется неоднородная трех периодическая марковская цепь пятого порядка. Каждое экзонное состояние состоит из к -состояний, где к- число моделей гена. Р(е) -вероятность экзона е в такой модели описывается выражением

P(e) = f_ja'_lP(e\M_l)

1=0

где а', - вероятности выбрать модель гена М, в позиции ДНК - /, P(e\MJ - вероятность экзона в соответствующей модели гена. Важно заметить, что вероятность выбрать модель а', - зависит от позиции экзона /. Выбирать число моделей к нужно совместно с другими параметрами марковской цепи. Для выбора параметров модели в используется следующая процедура белки, предсказанные на каждом шаге обучения, сравниваются с помощью программы быстрого локального выравнивания (BLAT [99]). Для этого конструируется неориентированный граф G=(V,E), где вершинами являются предсказанные белки. Если существует существенное сходство между двумя белками, то соответствующие вершины соединяются ребром. Связанная компонента графа G определяет модель гена М„ число связанных компонент определяет число моделей к. Параметры для экзонов определяются для каждой компоненты отдельно. Для того чтобы избежать нулевых частот, для некоторых олигонуклеотидов в экзоне используются псевдоотсчеты, вычисленные на обучающем множестве аннотированных генов. Вероятности моделей генов а', выбирались так, чтобы для каждого экзона использовалась ровно одна модель, так чтобы Р(е)=тах, Р(е\М). Последовательность состояний скрытой марковской цепи для каждой последовательности ДНК выбирается случайно в соответствие с предсказательным распределением (3), как в алгоритме из предыдущего раздела [44]. Метод, предложенный в работе [65] для предсказания генов, обладает высокой чувствительностью (0.9 нуклеотидная) и специфичностью (0.89 нуклеотидная) для последовательностей генома человека, которые сравнивались с гомологичными последовательностями из геномов млекопитающих, секвенируемых в рамках проекта NISC [102] Качество предсказания алгоритма не зависит от геномных перестановок.

1.3 Альтернативный сплайсинг (АС).

Появление новых мощных методов исследования, таких как секвенирование маркеров экспрессии (EST) и микронаборы гибридизационных зондов (microarrays, MA), изменило относительно простую картину интерпретации генетического кода: от ДНК через мРНК к

первичной структуре белка («один ген - один белок»). Оказалось, что для многих генов существует несколько альтернативных продуктов сплайсинга первичного транскрипта пре-мРНК. В результате альтернативного сплайсинга образуются изоформы - мРНК содержащие разные наборы экзонов. Изоформы кодируют функциональные варианты белка или являются «ошибками» сплайсинга, в этом случае они быстро деградируют с помощью специального механизма. Для некоторых генов показано, что "ошибочные" изоформы могут выполнять регуляторные функции. Секвенирование полных геномов совместно с данными по экспрессии генов позволили оценить количество генов в геномах различных организмов. Современные оценки числа генов в геномах человека и мыши составляют около 30 - 40 тысяч [1], согласно более строгим оценкам - 20 - 25 тысяч [3,4] и не отличаются существенно между видами. Из этих оценок следует, что геномы даже сложных организмов, таких как высшие млекопитающие, содержат относительно небольшое число генов. Альтернативный сплайсинг является фундаментальным механизмом, так как он обнаружен в различных таксономических группах эукариот: у растений, дрожжей, беспозвоночных, млекопитающих. При относительно небольшом числе генов альтернативный сплайсинг увеличивает разнообразие протеома -совокупность белков кодируемых в геноме.

Типы элементарных альтернатив.

Рисунок 1.3.1. Основные типы элементарных альтернатив. А -альтернативный донорный сайт вариантами, или альтернативная терминация транскрипции; Б -альтернативный акцепторный сайт или альтернативная инициация транскрипции; В - удержанный интрон; Г - кассетный экзон; Д - чередующиеся экзоны

Принято выделять следующие типы элементарных альтернатив (рис 1.3.1):

кассетный экзон, альтернативные 5'- и 3'- сайты, чередующиеся экзоны, удержанный интрон, альтернативная инициация транскрипции и альтернативное полиаденилирование. У млекопитающих наиболее распространенной альтернативой является кассетный экзон [24,69]

Частота альтернативного сплайсинга.

Достаточно долго альтернативный сплайсинг считался исключением из общего правила По мере накопления EST оценка доли альтернативно сплайсируемых генов в геноме человека составила: от 35% [2] до 60% [103] Однако, есть основания считать, что доля генов, для которых характерно более одной изоформы, существенно недооценена. На основании экспериментов с микронаборами гибридизационных зондов альтернативный сплайсинг показан для 74% генов [см. 5], 79 - 88% генов на хромосомах 21 и 22 [6]. В настоящее время применение методов гибридизационного анализа к исследованиям альтернативного сплайсинга является перспективным направлением. Данный метод лишен нескольких существенных недостатков, свойственных EST, таких как неравномерность покрытия гена и зависимость покрытия от уровня экспрессии гена. Большая часть EST имеет тенденцию лучше покрывать концы мРНК: 5'- и в большей степени 3'- конец, в то время как середина гена может быть недостаточно покрыта для обнаружения альтернативного сплайсинга. Первоначально предполагалось, что доля альтернативного сплайсируемых генов увеличивается с усложнением организма, однако это верно лишь отчасти. Более 95% генов Saccharomyces cerevmae не содержат интронов, и всего для нескольких генов показан альтернативный сплайсинг. В геноме дрожжей обнаружено очень мало генов, продукты которых регулируют сплайсинг, по сравнению с геномами Drosophila melanogaster, С elegans и человека [7]. Доля альтернативно сплайсируемых генов в геноме Arabidopsis thaliana составляет от 1,5 до 6,5% [9]. На основе анализа EST Brett с соавторами [8] показали, что частота альтернативного сплайсинга в геноме Arabidopsis существенно меньше, чем в геномах других 6 видов

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

Влияние на функцию белка.

Один из наиболее интригующих вопросов, как АС влияет на структуру и функцию белка. К сожалению, в настоящее время практически нет данных о пространственной структуре альтернативных изоформ В базе данных PDB, хранящей информацию о белках, для которых известна пространственная структура, представлено менее 10 сплайсированных изоформ. Для некоторых белков экспериментально показано различие 3D структуры их изоформ. Однако для большинства белков пока еще не существует прямого способа предсказать влияние того или иного события альтернативного сплайсинга на их структуру и функцию. Исследование функции альтернативного сплайсинга как явления в целом на уровне протеома требует моделирования 3D структуры изоформ или использования непрямых методов, таких как анализ корреляций альтернативного сплайсинга и известных характеристик белка. Среди непрямых методов наиболее широко используются два. Первый из них заключается в определении областей в доминирующей (базовой) изоформе, в которые альтернативный сплайсинг вносит изменения. Положение таких областей относительно элементов вторичной структуры или доменов в базовой изоформе позволяет делать предположения относительно их влияния на белок.

Второй метод состоит в определении степени выраженности альтернативного сплайсинга в генах, кодирующих белки с общими свойствами. Например, по данным Johnson с соавторами [5] среди генов с наибольшим числом тканеспецифичных альтернативных экзонов преобладают гены из следующих функциональных групп: передача сигналов посредством рецепторов, связанных с тирозин-киназами, и регуляция посредством малых АТФ-аз. Среди генов, для которых характерен тканеспецифичный

альтернативный сплайсинг, локализованный вблизи доменов в белке, преобладают гены, участвующие в формировании цитоскелета [16].

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

В работе Wang с соавторами [10] рассматривали положение альтернативных областей относительно вторичной структуры белка базовой изоформы. Выборка содержала 1209 гена из SWISS-PROT, для которых были указаны альтернативные изоформы. Для 351 гена была известна вторичная структура (PDB), для остальных вторичная структура базовой изоформы была предсказана с очень высокой достоверностью (Р « 10⁸). Обнаружено, что границы делеций и вставок имеют статистически значимую тенденцию избегать спирали и /?-складки и приходятся в основном в область петель. Для замен в белке такой тенденции не выявлено. Это может быть связано с отбором на замены, которые не нарушают вторичную структуру. Около 40% границ делеций, приходящихся на элемент с известной вторичной структурой, совпадают с границами этого элемента Это согласуется с тем, что альтернативный сплайсинг имеет тенденцию к удалению доменов целиком [11]. Для аминокислот была определена степень доступности для растворителя, которая зависит от их положения внутри или на поверхности вторичной структуры. Большинство событий альтернативного сплайсинга (75%) приходиться на поверхность белка. Авторы [10] предполагают, что альтернативный сплайсинг подвержен давлению отбора против сильных изменений вторичной структуры белка. Альтернативный сплайсинг лишь «адаптирует» структуру базовой изоформы для выполнения новой функции. Для того, чтобы проверить эту гипотезу, было проведено моделирование вторичной структуры

альтернативных изоформ и сравнение со структурой базовой изоформы. Так как большинство альтернативных областей (60%) имели небольшой размер (<50 аминокислот), изменения структуры белка были в основном локальными.

Однако даже незначительное изменение аминокислотной последовательности может приводить к существенным изменениям вторичной структуры. Garsia с соавторами [12] исследовал альтернативный сплайсинг в белке Piccolo (Aczonin) Альтернативный сплайсинг встраивает в в домен СгА с известной вторичной структурой девять аминокислот. Моделирование длинной изоформы показало, что вставка приходится в область петли на поверхности белка вдали от Са связующего домена. С помощью ЯМР было установлено, что вставка приходится внутрь Р-цепи в ядре базовой изоформы. В результате фрагмент аминокислотной последовательности со вставкой сильно изменяет Са²⁺ связующий домен так, что меняется белок-белковое взаимодействие. Длинная изоформа, в отличие от базовой, образует Са- зависимый димер. Это очень показательный пример относительно того, как малые изменения последовательности белка могут значительно изменять его труктуру и функцию, и насколько осторожно надо относиться к предсказаниям пространственной структуры альтернативных изоформ. В работе Offinan с соавторами [13] показано на выборке из 21 гена отсутствие корреляции между альтернативно сплайсируемыми областями белка и положением сайтов белок-белкового взаимодействия. Однако это не означает, что альтернативный сплайсинг не влияет на эти сайты непрямым способом, наподобие того, как это было показано в [12] для белка Piccolo.

Тем не менее, существует множество примеров, когда альтернативный сплайсинг удаляет домены в белке. Некоторые типы доменов альтернативный сплайсинг удаляет из базовой формы более часто, чем это могло бы быть по случайным причинам [14], например KRAB домены и анкириновые повторы. Наиболее частый случай, когда удаляется домен, участвующий в белок-белковом взаимодействии, так что изменяется

характер взаимодействия с другими белками. Так, секреторная форма белка может образовываться в результате удаления трансмембранного домена в базовой изоформе [15].

Одно из направлений исследования альтернативного сплайсинга, привлекающее пристальное внимание биологов, является тканеспецифичный альтернативный сплайсинг. Основываясь на EST данных, было выявлено множество примеров тканеспецифичного сплайсинга [104] в различных генах. Однако результаты сильно зависят от объема доступных для анализа кДНК/EST транскриптов из разных тканей. Существуют примеры, когда экспрессия экзонов, считавшихся тканеспецифичными, позднее наблюдалась и в других тканях по мере появления новых EST последовательностей (BRCA1)[104].

Эксперименты с использованием микронаборов гибридизационных зондов могут давать более надежную информацию относительно специфичности включения альтернативных экзонов и позволяют взглянуть на явление в целом. В настоящее время предложены схемы экспериментов, позволяющие делать количественную оценку уровня экспрессии гена и отдельных событий альтернативного сплайсинга, обладающие хорошей чувствительностью и специфичностью. Так, в работе Pan с соавторами [16] были выбраны гибридизационные зонды на 3126 альтернативных кассетных экзонов, известных по кДНК/EST данным, которые были идентифицированы в 2647 различных генах. Был проведен анализ на 10 зрелых тканях мыши. Авторы определяли такие параметры, как доля включения в мРНК каждого кассетного экзона и уровень экспрессии каждого гена в различных тканях. Было показано, что для большинства кассетных экзонов доля их включения в мРНК не значительно изменяется от ткани к ткани. Близкие по происхождению ткани обладают близкими профилями экспрессии генов и профилями включения альтернативных экзонов. Тканеспецифичная экспрессия и альтернативный сплайснг независимо действуют на различные гены и на функциональные группы генов. Нет значимого перекрытия, по сравнению со случайным, между множеством генов для которых характерна тканеспецифичная экспрессия и множеством генов, содержащих

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

Деградация мРНК, вызванная преждевременной терминацией трансляции (NMD).

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

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

Детали процесса пока не известны, мРНК подвергается экзонуклеазному расщеплению, если СТОП-кодон находится на расстоянии больше 50-55 нуклеотидов от положения 5'-сайта последнего экзона [26]. У Drosophila экзонуклеазное расщепление инициируется эндонуклеазой, разрезающей мРНК в окрестности СТОП-кодона [105]. На расстоянии приблизительно в 20-24 нуклеотида от каждой границы между экзонами находится стабильный комплекс белков (EJC от exon junction complex, включающий белки Upf2, Upf3, Upf3x), который участвует в распознавании преждевременного СТОП-кодона Комплекс белков EJC остается связанным с мРНК. Процесс распознавания СТОП кодона происходит во время, так называемой предварительной трансляции. Рибосома, двигаясь вдоль мРНК, счищает белковые комплексы, маркирующие сайты сплайсинга, РНК геликаза Upfl, непосредственно взаимодействует с белком Upf2 в составе комплекса EJC, если он не был удален рибосомой. Как предполагают, Upfl так же учавствует в терминации трансляции. Процесс деградации мРНК начинается с удаления 5'-кэпа и поли-А хвоста, после чего мРНК разрушается экзонуклеазами в обоих направлениях 5'-3' иЗ'-5\

Альтернативный сплайсинг совместно с NMD может участвовать в регуляции экспрессии генов. Для ряда генов такая регуляция была показана экспериментально: рибосомальные гены RP3, RPlOa, и RP12 С Elegans [27] и человека [30], глютаминаза, рецептор-2 фактора роста фибробластов, фактор сплайсинга SC35 и другие. По оценки Lewis с соавторами [28] около 35% изоформ мРНК различных генов, 144 альтернативных мРНК изоформы в базе данных SWISS-PROT [29] и 4,3% экспериментально проверенных мРНК в базе данных RefiSeq являются изоформами, подверженными NMD-деградации. Авторы предположили, что регуляция экспрессии генов, в основе которой лежит альтернативный сплайсинг сопряженный с NMD, является широко используемым,

эволюционно консервативным механизмом. Однако данные, полученные в эксперименте Pan с соавторами [31] с использованием микронаборов гибридизационных зондов противоречат этой гипотезе. В эксперименте было показано, что доля вариантов сплайсинга, приводящих к появлению преждевременных СТОП- кодонов (ПСК), очень мала. С помощью малых интерферирующих РНК (siRNA) экспрессия Upfl была снижена до 5% от нормальной. Однако менее 10% вариантов сплайсига, приводящих к ПСК, увеличили свою долю среди мРНК транскриптов более чем на 15%. Снижение экспрессии Upfl не влияло на экспрессию генов Таким образом, варианты сплайсинга, содержащие ПСК, являются редкими и возможно ошибочными вариантами мРНК.

Эволюционное значение альтернативного сплайсинга.

Современное представление об альтернативном сплайсинге как механизме, ускоряющем эволюцию, сложилось на основании значительного количества наблюдений. Эволюционный процесс может создавать новые функции двумя способами, дуплекацией генов или созданием новых экзонов для существующих генов. Мобильные элементы являются хорошо известным источником новых экзонов в геноме человека. Например, Л/и-повтор входит в состав около 4% кодирующих экзонов человека. Процесс превращения Л/и-повтора в экзон достаточно позднее явление, характерное для приматов [7] На основании данных кДНК/EST Sorek с коллегами [17] показали, что все экзоны, содержащие Alu, являются альтернативно сплайсируемыми В большинстве случаев мобильный элемент встраивается в интрон и затем постепенно в результате случайных мутаций превращается в экзон. В своей работе Singer с соавторами [18] продемонстрировали эволюцию альтернативного первого экзона 1а гена p75TNFR. Выравнивание этого экзона и 5'-фланкирующей последовательности соответствующего интрона с консенсусом семейства повторов Alu Jo свидетельствует, что нуклеотидная замена А --> G привела к возникновению старта трансляции (ATG), замена С Т привела к формированию сайта сплайсинга (GT), а удаление 7 нуклеотидов к появлению открытой

рамки считывания. Филогенетический анализ последовательности первого экзона 1а гена p75TNFR в приматах демонстрирует, что вставка Alu произошла 40-58 миллионов лет назад до расхождения ветвей обезьян старого и нового света, затем произошла мутация А --> G (возникновения ATG). Появление сайта сплайсинга и открытой рамки произошло 25-40 миллионов лет назад уже после расхождения ветвей обезьян старого и нового света Основной функцией молодого мобильного элемента до появления открытой рамки считывания могла быть регуляция экспрессии гена, сопряженная с механизмом деградации мРНК, вызванной преждевременным СТОП-кодоном (NMD).

Другим источником альтернативных экзонов является дупликация существующего экзона Было показано большое число генов человека и животных, содержащих дуплекацию экзонов. Примерно 10% альтернативных чередующихся экзонов произошли в результате дуплекации [19]. Одновременное включение обоих чередующихся экзонов в мРНК обычно нарушает рамку считывания.

Связь альтернативного сплайсинга с возникновением / утратой экзонов в процессе эволюции, по-видимому, не ограничивается примерами дуплекации и превращением мобильных элементов в экзоны. Например, более 98% константных экзонов, присутствующих во всех мРНК гена, консервативны между ортологичными генами человека и мыши. Кассетные экзоны, которые включаются в более 50% кДНК/EST, консервативны практически в той же степени, как и константные экзоны. Только 30% кассетных экзонов, которые включаются в транскрипт реже чем в половине случаев, являются консервативными [20]. Эсперименты с использованием микронаборов гибридизационных зондов подтверждают эти результаты [16]. Более того, эти эксперименты показали, что частота включения экзона в мРНК положительно связана с локализацией соответствующей альтернативной области вблизи или непосредственно в области консервативных доменов. Редкие экзоны (с частотой включения меньше 33%) вызывают изменения в областях консервативных доменов в 6% случаев, в то время как

экзоны с высокой частотой включения ассоциированы с доменами в 15% случаев. Таким образом, экзоны с высокой частотой включения в большей степени связаны с регуляцией функции белков, по сравнению с редкими экзонами [16].

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

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

Таким образом, альтернативный сплайсинг позволяет открыть нейтральный или почти нейтральный путь для эволюции посредством появления новых экзонов [7]. Нейтрализующий эффект альтернативного сплайсинга по степени влияния в некоторой степени близок к влиянию диплоидии. Замечено, что гены на аутосомных хромосомах млекопитающих имеют большую долю мРНК изоформ (9%), содержащих преждевременный СТОП-кодон, по сравнению с генами на X хромосоме (3%) [23]. Так же наблюдалось [23], что редкие изоформы значительно более часто подвергаются NMD деградации по сравнению с доминирующими изоформами.

Алгоритмы анализа EST данных.

До недавнего времени, изучение явления альтернативного сплайсинга практически полностью основывалось на данных кДНК/EST. Экспериментальные методики, использующие микронаборы гибридазационных зондов, позволяют получить данные относительно того, является ли некоторый экзон альтернативным, и в какой доле мРНК он присутствует. К сожалению, с помощью микронаборов гибридизационных зондов нельзя определить для двух экзонов присутствуют ли они в составе одной мРНК или разных. Кроме того, при разработке зондов часто используются данные об экзон-интронной структуре гена, полученные на основе анализа EST [5,16].

EST могут рассматриваться как случайные фрагменты мРНК [25] средней длинной ~300 - 500 нуклеотидов [32], имеющие более высокую плотность покрытия на концах гена [16].

Последовательности EST, полученные современными методами, содержат 1-2% ошибок, основными источниками которых являются: (1) загрязнение последовательностью вектора, в котором производилось клонирование, (2) плохо прочитанные нуклеотиды, (3) химерные кДНК клоны. Проблема загрязнений вектором решается относительно просто. Так как последовательность вектора известна, то ее легко идентифицировать и отрезать, тем более что вставки вектора находятся на концах кДНК. Проблема идентификации химерных EST, содержащих сплайсированные вместе части мРНК разных генов, решается выравниванием с геномной последовательностью. Если геномная последовательность недоступна, то химеры являются мостами между непересекающимися кластерами EST, однако мостом может быть EST, объединяющая два кластера, принадлежащих разным концам одной мРНК [32].

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

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

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

Граф сплайсинга где вершины соответствуют экзонам, а ребра графа - интронам, использовался в работах Lee [34,35]. Для построения графа сплайсинга необходимо множественное выравнивание EST, относящихся к одному кластеру, или выравнивание с геномной последовательностью. Экзоны в разных последовательностях EST, совпадающие по начальной и конечной координате, представляются одной вершиной на графе. Ребра соединяют вершины в том случае, если, хотя бы в одной последовательности, соответствующие экзоны являются соседними. Задача выделения элементарных альтернатив (рис. 1.3.1) - взаимно исключающих вариантов сплайсинга некоторой области пре-мРНК, формулируется как задача поиска подграфов определенного вида.

В работах Heber [36] и Maddle [33] с соавторами используют графы слов - другой способ построения графов сплайсинга. Вершинами графа слов являются все слова длины к-1, а ребрами все слова длины к (к-слова), встречающиеся в кластере EST. Если некоторое k-слово встречается в последовательности более одного раза, то в графе образуется цикл. Специальные усилия следует предпринять для коррекции плохо прочитанных или неправильно прочитанных оснований в последовательностях EST. Heber [36] с соавторами решают эту задачу с помощью построения множественного

выранивания. Ошибки секвенирования выявляются как единичные замены в достаточно длинных (>30 нуклеотидов) выровненных областях EST. Выравнивание внутри таких областей должно содержать более 95% константных позиций. Небольшие делеции/вставки внутри выравнивания пары последовательностей EST могут так же быть признаны ошибками секвенирования, а не вариантами альтернативного сплайсинга. Алгоритм сборки EST заключается в перечислении всех возможных путей на графе.

Алгоритм Maddle [33] решает проблему повторяющихся k-слов и генерирует набор консенсусов - путей, наиболее подтвержденных набором EST. Для этого на ребрах графа слов задаются веса равные числу последовательностей EST, в которых встречается соответствующее k-слово. Также на ребре хранится ассоциативный массив, отображающий идентификаторы последовательностей во множество координат к - слова в соответствующей последовательности. Путь на графе строится жадным способом. Для этого сначала алгоритм выбирает ребро с максимальным весом. Стартуя с этого ребра в обоих направлениях, алгоритм присоединяет ребра, которые соответствуют словам, содержащимся в максимальном числе последовательностей из множества последовательностей на предыдущем ребре. Алгоритм избегает зацикливания следующим способом: при присоединении нового ребра выбираются только те последовательности, для которых это ребро увеличивает текущую координату. Таким образом, последовательность консенсуса может содержать повторяющиеся слова. Перед выбором первого ребра алгоритм уменьшает веса на ребрах, представляющих слова, повторяющиеся в последовательностях по нескольку раз, так как для таких слов нельзя однозначно определить текущую координату. Для построения следующего консенсуса выбирается ребро графа с максимальным весом, которое не включено ни в один из предыдущих.

Третий способ построения графа сплайсинга ClusterMerge использовался в работе Eyras с коллегами [37]. Последовательности EST, выровненные с геномом, были

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

1. Последовательности EST упорядочиваются относительно генома по 5' концу по возрастанию. Если 5' концы совпадают, то по убыванию координаты 3' конца

2. Пусть у расширяет х Ребро х->у присутствует в графе, если не существует EST w, которая расширяет х, при этом сама расширяется последовательностью у (х-> w-> у), и если не существует EST w, которая включает в себя х, при этом сама расширяется последовательностью у (w=>x & w->y)

3. Пусть л: включает в себя у. Ребро х=> у присутствует в графе, если не существует EST w, котороая включает в себя у, при этом сама включена в х (х=> w =>у), и если не существует EST w, которая расширяется л;, при этом включает в себя у (w->x & w=>y).

Определенные подграфы в принципе невозможны, например

А

х с=> у

^w С=> Z

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

> z

_Г-п

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

не расширяется другими вершинами.

_г

s _

На рисунке вершины х, у, z являются корнями, а / и s являются краевыми вершинами. Сборка EST осуществляется следующим образом - алгоритм конструирует минимальный набор множеств EST, которые могут быть частями одной мРНК, при этом каждое такое множество должно включать в себя максимальное число EST. Формально задача решается перечислением всех путей из корневых вершин в краевые по ребрам расширения, при

этом последнее ребро на пути не должно быть частью треугольника. Для каждой вершины на пути необходимо добавить все достижимые вершины по ребрам включения. Для графа Г нерасширяемыми множествами EST являются {x,t}, {y,tj, {z,s,tj. Множество {z,tj может быть расширено за счет s (ребро z-> t является частью треугольника z, s, t).

Сборка EST - одно из наиболее очевидных применений графа сплайсинга Либеральный подход Heber [36, 38], в результате которого генерируются все изоформы, гарантирует, что ни один из возможных вариантов не будет пропущен. При таком подходе для некоторых генов генерируется очень большое число изоформ (более 5000). В основе либерального подхода лежит предположение о том, что интроны вырезаются, преимущественно, независимо друг от друга Другой подход, когда генерируется минимальное число изоформ, как в алгоритме ClusterMerge, оптимизирует специфичность.

Эвристический алгоритм HeaviestBundle, предложенный Lee [35], в отличие от ClusterMerge [37] не гарантирует нахождение точного решения задачи построения минимального числа изоформ, объясняющих входные данные. Строится граф сплайсинга ребрами которого являются интроны, а вершинами экзоны. Веса на ребрах равны числу последовательностей EST. Среди EST выбирается максимально длинная последовательность - шаблон Веса ребер графа, которые соответствуют интронам шаблона, увеличиваются в f (f=1000) раз. Для всех EST, хорошо выравнивающихся с шаблоном (более 93% совпадающих символов и делеции короче 6 нуклеотидов), веса ребер так же увеличиваются в f раз. Алгоритм динамического программирования HeaviestBundle ищет путь на графе с максимальным весом, на каждом шаге, совершая локально оптимальный выбор. Для каждой вершины выбирается входящее ребро с максимальным весом, если существует несколько таких ребер, то из них выбирается ребро, исходящее из вершины с наибольшим весом. Текущей вершине приписывается сумма весов выбранного входящего ребра и веса в предыдущей вершине. Для наглядности

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

Динамическое программирование подразумевает, что ребра графа выбираются независимо (независимый сплайсинг). Алгоритм HeaviestBundle использует увеличение весов ребер в соответствие с последовательностью шаблона, для того, чтобы заставить оптимальный путь проходить через шаблон, или максимально "близко" к шаблону. Таким образом, шаблон расширяется в обоих направлениях в соответствии с принципом максимального правдоподобия. В отличие от алгоритмов Eyras[37] и Maddle [33], HeaviestBundle [35] не гарантирует, что последовательность консенсуса не будет содержать экзонов, сочетание которых никогда не наблюдалось во входных данных. По времени работы алгоритм HeaviestBundle более производителен, по сравнению с похожим алгоритмом Maddle, так как не требуется вычислять пересечение множеств идентификаторов последовательностей EST на ребрах графа.

Несмотря на то, что HeaviestBundle успешно применялся для исследования
альтернативного сплайсинга [39, 23, 14], Lee с соавторами [40] предложили новый
алгоритм, который корректно обрабатывает скоординированные события

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

наблюдаемым фрагментам. В рамках, предложенного подхода может быть объединена информация различных источников данных: EST и микронаборов гибридизационных зондов. Рассмотрим более подробно данный алгоритм. Пусть на графе сплайсинга перечислены все возможные пути - изоформы: //, І2, , h- Обозначим множество наблюдений (последовательностей мРНК) Oi, О2, , Ojv , тогда зададим матрицу сопряженности изоформ и наблюдений Z=(z,k\ z, к~0,1), i=l N (число наблюдений), к=1 К (число изоформ). Предположим, что мы наблюдаем полноразмерные последовательности. Элемент матрицы z, * равен 1, только в том случае, если наблюдаемая последовательность О, совпадает с изоформой 7*. В каждой строке матрицы Z возможна только одна единица, остальные значения равны нулю. Вероятность наблюдать ту или

^{_^}_Л^{,.,„„,„.„„ ,.,±fe.,..}

действительности матрица Z не является наблюдаемой, мы наблюдаем другую матрицу сопряженности Y=(y,k\ у,к~0,1), i=l N, к=1 К между EST и изоформами. Если какой-либо элемент у, к равен нулю, то соответствующий элемент z,k тоже равен нулю, однако если у, к равен 1, то z,k может быть равен 1, а может быть нет. Задача состоит в определении частот изоформ по матрице F. Обозначим множество искомых частот в=(рі,

N К

Р2, > Рк) ^и запишем правдоподобие наблюдаемых данных \(Y\0) = ^\og(^y_lkp_k).

1-І k=\

Задача состоит в определении методом максимального правдоподобия элементов множества частот изоформ в, решение которой не выражается в аналитическом виде. Авторы [40] применили итерационный алгоритм максимизации математического ожидания.

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

последовательностей (около 10) наблюдается значительная вариация расстояния между изначальными частотами и восстановленными. При числе последовательностей около 250 алгоритм хорошо восстанавливает изначальное распределение/?* (вариация менее 0,1)

Имея набор предполагаемых изоформ, для каждого гена можно построить множество кодируемых белков Для этого в каждой изоформе мРНК определяется положение кодирующей последовательности, чаще всего как наиболее длинной открытой рамки считывания. Например, такой подход использовался при построении базы ASP альтернативных белковых изоформ [39], однако он может давать ошибочный результат для некоторых генов. Авторы приводят пример мРНК, кодирующей короткий белок, для которой правильная рамка считывания не является самой длинной рамкой. Для генов с экспериментально известными мРНК можно определять рамку считывания в изоформах собранных из EST по пересечению с кодирующей последовательностью мРНК, как было сделано в работе Lewis с соавторами [28]. Сгенерированная изоформа может быть кандидатом для деградации с помощью NMD, в этом случае она, скорее всего, не кодирует биологически активную форму белка. При построении базы ASP, изоформы, рамка считывания в которых была короче, чем половина от кодирующей последовательности известной мРНК, или расстояние от 5' границы последнего экзона до СТОП-кодона более 50 - 55 нуклеотидов считались изоформами, подверженными NMD.

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

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

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

В принципе, для аннотации генома, с помощью сплайсового выравнивания EST, можно использовать последовательности, полученные из других достаточно близких организмов. Так же структура гена может быть восстановлена с помощью EST генов, относящихся к тому же семейству параллогов. Основная проблема заключается в том, как определить границы экзонов, если в выравнивании существуют делеции/вставки и замены нуклеотидов. Usuke и Brendel [41,42] предложили алгоритм GeneSeqer, использующий для аннотации геномов растений EST ортологичных генов или параллогов. В основе алгоритма лежит скрытая Марковская цепь парных состояний, которая генерирует сплайсированное выравнивание с геномной последовательностью. Для оценивания сайтов сплайсинга используются статистические модели "правильных" и "ложных" сайтов. Ложные сайты встречаются в окрестностях правильных сайтов в кодирующей и некодирующей последовательности ДНК. Кандидаты в сайты классифицируются с помощью функции правдоподобия в одну из категорий: "правильный" сайт в рамке 0,1,2; "ложный" сайт в рамках 0, 1, 2 и ложный сайт в некодирующей последовательности [42]. Состояния Марковской цепи соответствуют классификации нуклеотидов геномной последовательности ДНК на принадлежащих экзонам и интронам. Переходы между состояниями цепи возможны только по кандидатам в сайты сплайсинга, соответствующие вероятности переходов зависят от статических свойств кандидатов. Для выровненных последовательностей EST используется алгоритм сборки полноразмерных транскриптов.

В работе Foissac и Schiex [43] предложен алгоритм EuGene, в основе которого лежит модель гена в виде ациклического ориентированного графа аналогичная моделям, принятым в статистических алгоритмах распознавания генов, основанных на Скрытой Марковской Модели (НММ). Все возможные разметки последовательности ДНК на экзон-интроные структуры генов представлены на графе. Источником графа называют вершину в которую не входит ни одного ребра, стоком - вершину из которой не выходит не одного ребра. Оптимальный путь на таком графе может быть найден с помощью алгоритма Витерби Алгоритм EuGene может предсказывать гены, основываясь только на геномной последовательности, а так же учитывать информацию о экзон-интронной структуре, полученную из сплайсового выравнивания EST. Для каждой последовательности EST, выровненной с геномом, строится дополнительная компонента графа, представляющая собой все возможные разметки геномной последовательности согласованные с выравниванием. Источники и стоки дополнительной компоненты, соответствуют концам выравнивания, соединяются с соответствующими вершинами основного графа, ребрами, ориентированными в направлении дополнительной компоненты. Алгоритм Витерби применяется дважды: по (5' 3') и против (3' 5') ориентации ребер. В результате первого прохода алгоритм находит оптимальные пути из источника в 3' стоки дополнительного графа и оптимальный путь через основной граф. Во время обратного прохода вычисляются оптимальные пути из стока графа в 3'-стоки дополнительной компоненты Таким образом, для каждой дополнительной компоненты графа определяется оптимальная экзон-интронная структура, согласованная с выравниванием EST. Применение EuGene для аннотации генома позволяет обойти ограничения, связанные с недостаточным покрытием EST. Для генов, не покрытых EST,, будет найдена оптимальная структура, кодирующая единственный белок. Для генов с альтернативными экзон-интронными структурами, полученными выравниванием EST, EuGene определяет соответствующий набор белков.

Среди алгоритмов, предсказывающих альтернативные экзон-интронные структуры без использования последовательностей кДНК и EST, существует два основных подхода: поиск субоптимальных экзон-интронных структур с помощью НММ и использование парной НММ для поиска альтернативного сплайсинга по отношению к ранее известной структуре гена. В работе [44] Cawly и Pachter предложили метод, позволяющий эффективно выбирать к- случайных экзон-интронных разметок в соответствии с апостериорным распределением вероятностей, заданном на Скрытой Марковской Модели геномной последовательности. Формальная постановка задачи - пусть задана парная Марковская Модель, содержащая N скрытых состояний (рНММ), пусть X? и Yi^U две геномные последовательности, кодирующие ортологичные белки, Т и U длины соответствующих последовательностей. Алгоритм позволяет случайно выбрать к последовательностей скрытых состояний Zi^L(l), с длительностями di^L(,) , i-l к из распределения вероятностей P(Z]^L,di^L\ Х? ,Yi^U) за время 0(NTU+kN(T+U)) и требует 0(N(T+U)) затрат памяти.

Эта работа имеет скорее теоретический интерес, так как для того, чтобы среди выбранного набора были различающиеся разметки, размер выборки должен быть достаточно большим к~1(Ґ Более серьезным возражением является то, насколько справедливо утверждение, что реальные альтернативные экзон-интронные структуры находятся в окрестности глобально оптимальной разметки? Однако авторы приводят пример применения алгоритма, когда был предсказан консервативный альтернативный сплайсинг (человеко/ мышь) гена тропомиазина ТРМ2.

Алгоритмы, в основе которых лежит парная Скрытая Марковская Модель, используются для предсказания экзон-интронной структуры гомологичных генов одновременно в паре геномных последовательностей [44,45,58]. Существует ряд характеристик, которые позволяют различать эволюционно консервативные константные

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

альтернативного сплайсинга [7]

Burge с соавторами применили рНММ, реализованную в алгоритме UNCOVER [46],

для построения глобального выравнивания пары ортологичньк интронов человека и мыши. Для обнаружения альтернативного сплайсинга, не встречавшегося в кДНК/EST, алгоритм применяется к парам интронов, которые находятся между константными экзонами. Последовательность скрытых состояний позволяет моделировать все возможные элементарные альтернативы: 5' и 3' альтернативные сайты (удлинение фланкирующих экзонов), кассетный экзон и удержанный интрон (рис.1.3.2).

Рисунок 1.3.2 Упрощенная схема НММ программы Uncover [46], предназначенной для предсказания консервативного альтернативного сплайсинга.

удержанный интрон

Фланкирующие сайты кассетного экзона, альтернативные сайты, удлиняющие константные экзоны, и константные сайты являются различными состояниями модели. Схожие и видоспецифичные кодирующие фрагменты, консервативные и неконсервативные части интронов так же представлены различными состояниями модели. Эмиссионные вероятности состояний были обучены на соответствующих выборках сайтов, экзонов и интронов ортологичных генов. Алгоритм UNCOVER обладает высокой чувствительностью (87%) и специфичностью (68%) к консервативным кассетным экзонам. Тестирование проводилось на проверенной вручную, выборке пар ортологичных интронов. Наличие консервативных кассетных экзонов, проверялось выравниванием с EST/кДНК человека и мыши. В экспериментах с использованием микронаборов гибридизационных зондов, показано, что из-за недостаточно покрытия на основании EST

невозможно удовлетворительно оценивать долю транскриптов, содержащих тот или иной вариант альтернативного сплайсинга [16]. Доля альтернативно сплайсируемых генов в геноме так же недооценена на основе только лишь EST данных [5]. По оценкам авторов, около 1900 интронов в геноме могут содержать неизвестные консервативные кассетные экзоны с низкой долей включения в мРНК Одно из возможных применений алгоритма UNCOVER - проверка паттерна альтернативного сплайсинга на консервативность при наличии EST только в одном из двух организмов.

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

Алгоритмы поиска генов можно разделить на два вида - 1) статистические алгоритмы, источником информации для которых служит сама последовательность ДНК и 2) алгоритмы, использующие сходство с белками, EST, кДНК или последовательностями ДНК, содержащими гомологичные гены До недавнего времени такое разделение было абсолютным, однако в настоящее время граница стала сильно размытой. Как показано в [61] программы статистического распознавания генов имеют хорошую чувствительность, однако специфичность таких программ сильно зависит от длины межгенных интервалов -спейсеров. Ложные экзоны, предсказанные внутри длинных спейсеров, свойственных человеческому геному, заметно снижают специфичность. В значительной степени это связано с тем, что в настоящее время не существует хороших статистических моделей, позволяющих определить границы генов. Программы, основанные на сходстве, решают обратную задачу восстановления структуры гена по продукту сплайсинга того же самого гена (в простейшем случае) либо родственного гена, известного в другом организме. Эти алгоритмы имеют высокую специфичность, которая зависит только от степени близости продукта гена, структуру которого надо установить, и известного из базы данных гомолога. Серьезным недостатком такого подхода является то, что он не позволяет картировать ранее неизвестные гены. Специальным случаем является предсказание генов, основанное на том факте, что кодирующие белок экзоны значительно более консервативны по сравнению с некодирующей ДНК. Поэтому можно найти экзоны сравнением данного фрагмента с геномной ДНК, содержащей гомологичные гены. Предсказательная сила такого подхода для аннотации генов с известными ортологами в других организмах очевидно очень высокая. Экзон-интронная структура ортологичных генов млекопитающих часто в высокой степени консервативна. Например, среди ортологов человека и мыши 86% имеют одинаковое количество кодирующих экзонов, 46% имеют одинаковую длину белок-кодирующей области, причем только 1,5% имеют одинаковую длину белок-кодирующей области и разное количество экзонов [58]. Однако, при аннотации генов в составе быстро эволюционирующих параллогичных семейств, в том числе видоспецифичных генов приходится использовать статистические алгоритмы

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

Подробный обзор алгоритмов, использующих сходство для предсказания генов, приведен в работе [49]. Такие алгоритмы позволяют идентифицировать 50-70% генов в новом геноме [47,48] и эта доля будет возрастать по мере секвенирования новых геномов. Тем не менее, применение статистических алгоритмов распознавания генов в процессе аннотации является обязательным этапом, особенно для организмов, имеющих большое значение для медицины и биотехнологии. Современные статистические программы распознавания генов эукариот и многие программы, предсказывающие гены прокариот, используют Скрытую Марковскую Модель (НММ) генома. Алгоритмы, основанные на НММ, превосходят в чувствительности и специфичности алгоритмы, использующие альтернативные подходы нейронные сети, линейный дискриминантный анализ, лингвистический анализ [47,50]. Для НММ существуют эффективные алгоритмы решения задач оптимизации. Одним из достоинств НММ является вероятностная интерпретация - алгоритм ищет экзон-интронную структуру, которая с наибольшей вероятностью соответствует последовательности ДНК При построении модели может быть использована любая вероятностная весовая функция сайтов, нуклеотидного состава экзонов и некодирующей ДНК. Далее мы подробно рассмотрим Скрытую Марковскую Модель эукариотического гена, которая лежит в основе статистических алгоритмов. Так же мы коснемся построения парной Скрытой Марковской Модели (рНММ) для выравнивания гомологичных последовательностей и задачи оценивания параметров модели - обучения.

Скрытая Марковская Модель эукариотического гена. Успешное применение НММ к задаче распознавания генов человека в работе Burge и Karlin [51], предложивших программу GenScan, дало начало целой серии программ, использующих аналогичную модель геномной ДНК: Genie [52], FGENESH [53], GeneMark.hmm [54] и другие [47]. Перечисленные программы различались в деталях -статистических моделях сайтов, кодирующей и некодирующей ДНК и организмах, для которых было произведено обучение. В настоящее время все перечисленные программы обучены для предсказания генов широкого спектра организмов. Программа GeneMark.hmm, например, первоначально была разработана для бактериальных геномов, затем модель была расширена на геномы эукариот. При сравнении качества предсказаний различные программы демонстрируют близкие результаты [55].

Предсказание генов в геномах низших эукариот

Низшие эукариоты имеют большое значение для науки, как модельные микроорганизмы, для медицины и сельского хозяйства, как важные патогены человека, животных и растений, продуценты антибиотиков и биологически-активных веществ. Так же они широко используются в пищевой промышленности и биотехнологии. К настоящему времени были секвенированы геномы некоторых низших грибов, в частности, Aspergillus mdulans, Neurospora crassa, Magnaporte grizea, всего около 20 наименований [Fungal Genetics Stock Center]. В будущем интерес к секвенированию грибных геномов вряд ли ослабеет. Анализ нового генома всегда начинается с автоматической аннотации, целью которой служит предсказание генов для последующего детального анализа. Программы обнаруживающие гены на основе информации о гомологии геномной последовательности и последовательностей белков, мРНК, EST позволяют идентифицировать 50-70% генов [47,48]. Статистические программы предсказания генов, в основе которых лежит скрытая марковская модель (НММ), как правило, нуждаются в переобучении для каждого нового генома. Применение таких программ является обязательным этапом аннотации, так как позволяет идентифицировать гены специфичные для организма. Обучающее множество, необходимое для оценивания параметров НММ, может быть построено на основании генов, предсказанных по гомологии. К моменту начала работы над проектом существовала потребность в программе предсказания генов, которая могла бы использовать широкий набор статистических моделей, подбираемых под каждый аннотируемый геном, и могла бы легко переобучаться для аннотации новых геномов. Нами была разработана программа GipsyGene [80,82,86,94], в основе которой лежит Скрытая Марковская Модель, которая обладает требуемыми свойствами к гибкости обучения и интегрируется в пакет программ автоматической аннотации геномов Применение модели сайта ветвления в программе GipsyGene позволило улучшить распознавание коротких интронов Aspergillus spp и Neurospora crassa. Так же в программе GipsyGene существует возможность использовать прокариотическую модель гена, которая предусматривает отсутствие сплайсинга, альтернативные сайты инициации трансляции, оперонную модель организации генома [86]. Для аннотации некоторых геномов низших грибов целесообразно использовать модель гена без сплайсинга, так геномы некоторых дрожжей содержат на 6000 генов немногим более 200 интронов [66,67]. Программа тестировалась на генах Aspergillus spp и Neurospora crassa, качество предсказания сравнимо с результатами тестирования GENSCAN [51] на одногенных фрагментах ДНК. Дополнительно, мы предприняли тестирование на фрагментах ДНК, содержащих несколько генов. GipsyGene статистическая программа распознавания генов. Программа GipsyGene может использовать как модель эукариотического генома, так и прокариотического. Скрытая Марковская модель генома эукариот подобна моделям Genscan [51], Fgenesh [53], GeneMark.hmm [54]. Модель генома прокариот может рассматриваться как ограничение эукариотической модели. В модель заложены следующие кодирующие и некодирующие состояния: одноэкзонный ген; начальный, внутренний и терминальный экзоны; интрон; межгенный участок (спейсер).

Вероятности донорного и акцепторного сайтов вычисляются как в GENSCAN [51]. В зависимости от того, имеет ли сайт значимые корреляции между позициями, выбирается статистическая модель Если анализ сайтов в обучающей выборке не обнаруживает корреляций между позициями, или эти корреляции статистически не значимы, например, из-за небольшого объема обучающего множества, то вероятность вычисляется по профилю сайта. Если большинство позиций сайта значимо коррелируют в основном только с соседними позициями, то применяется WAM модель - вероятность основания в данной позиции зависит от основания в предыдущей позиции Эта модель применяется только для акцепторного сайта Если донорный сайт имеет значимые корреляции как с соседними, так и с удаленными позициями, его вероятность вычисляется по профилю в зависимости от оснований, стоящих в позициях, которые наиболее сильно взаимодействуют с другими позициями сайта - MDD модель.

Так как для многих грибных геномов сложно сформировать обучающую выборку достаточного объема, для вычисления вероятности кодирующей (некодирующей) части мы реализовали несколько различных моделей. Для вычисления автоматически выбирается наиболее «сильная» модель из доступных. Кодирующая ДНК может быть смоделирована с помощью 1) статистики кодонов; 2) марковской цепи первого порядка для аминокислот и статистики синонимичных кодонов, соответствующих аминокислотам; 3) трех-периодичных марковских цепей третьего-пятого порядков. Модели для некодирующей ДНК - марковские цепи первого, третьего и пятого порядков.

Для улучшения предсказаний 3 -концов интронов была разработана модель сайта ветвления, учитывающая ряд перечисленных свойств. Потенциальный сайт ветвления оценивается по профилю, распределению расстояний до акцепторного сайта, распределению количества динуклеотидов AG между сайтом ветвления и акцептором. Профиль и распределения конструируются на множестве слов, наиболее похожих на сайт ветвления Saccharomyces. Процедура построения профиля сайта ветвления предложена в [84], мы увеличили число позиций, включенных в профиль с 5 до 7, и добавили к профилю псевдоотсчеты. Псевдоотсчеты корректируют профиль, в предположении, что при его построении некоторые слова мало представлены или совсем не представлены. Для того чтобы построить профиль сайта ветвления, в окне [-/, -г] ([-40,-2] для Aspergillus spp и Neurospora crassa) левее акцепторного сайта осуществлялся поиск наилучшего слова, удовлетворяющего консенсусу [CT]T[AG]a[CT], идеально - СТАаС. Для кандидатов допускались отклонения от идеального слова только в двух из трех вырожденных позиций На множестве кандидатов строилось распределение расстояний до акцепторного сайта. Применение этой же процедуры к случайным последовательностям задает уровень шума - среднюю частоту встречаемости слова в случайной последовательности. Далее, окно поиска уточнялось так, чтобы частота в окне превышала уровень шума На множестве кандидатов, найденных в уточненном окне поиска, строилась частотная матрица размером 7 позиций (добавлены две позиции 5 -левее консенсуса), распределение расстояний до акцептора и распределение динуклеотидов AG. К каждому элементу полученной частотной матрицы добавлялось число, пропорциональное квадратному корню из числа последовательностей - псевдоотсчеты. В состоянии НММ "5 сайт экзона" генерируются нуклеотиды: собственно акцепторного сайта сплайсинга и соответствующего сайта ветвления. Сайт ветвления определяется как мотив с наибольшим весом в окне [-1.3/,-/ ] от каждого динуклеотида AG, где / и г - левая и правая координаты окна в котором сайт ветвления встречается чаще, чем в случайной последовательности. Вес мотива определяется как логарифм отношения правдоподобия модели сайта ветвления к интроннои модели: WBPS = log( -). Были Р{х Introri) введены обозначения: х - мотив, P(x\BPS) и P(x\Intron) вероятности мотива по профилю сайта ветвления и модели внутренней части интрона соответственно, d и NAG - количество нуклеотидов и число динуклеотидов AG между сайтом ветвления и кандидатом в акцепторный сайт соответственно.

Координированный альтернативный сплайсинг

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

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

Мы применили алгоритм выделения соседних альтернатив для анализа базы данных EDAS-I [92,93, 95] Для 630 генов было получено 916 матриц сопряженности вариантов сплайсинга в соседних альтернативах. Для того чтобы обнаружить координированный альтернативный сплайсинг для каждой таблицы был применен точный тест Фишера. Для 60 генов был обнаружен координированный альтернативный сплайсинг с уровнем значимости pval 0,1. Пороговое значение уровня значимости 0,1 было выбрано как ожидаемое, исходя из случайного распределения. Нами было обнаружено несколько соседних альтернатив, сплайсинг в которых был значимо координирован даже после учета поправки Бонферони на множественное применение теста Фишера. Поправка Бонферони заключается в том, что выполняется следующее условие: minfpval,}N, 0,l. Где минимум берется по всем значениям уровня значимости (pval,), полученным при однократном применении теста к каждой таблице i Nt. На рисунке 3.2.1 приведено пять примеров генов, обладающих координированным альтернативным сплайсингом достоверным на уровне pval 0,007. Для трех генов (CNOT8, РСВР2, BRRN1) различие длин вариантов сплайсинга кратно трем. Экзон в 5 -альтернативе гена CNOT8, длинной 189 нуклеотидов находится в 5 -нетранслируемой области. Экзоны в 5 -альтернативах генов РСВР2, BRRN1 находятся в кодирующей области, но не приводят к появлению преждевременных СТОП-кодонов (ПСК). Включение экзона в 5 -альтернативе гена MGAT4B приводит к появлению ПСК внутри предпоследнего экзона на расстоянии меньше 50 нуклеотидов от донорного сайта, что исключает активацию механизма NMD. Таким образом, механизм NMD не может быть объяснением недостатка числа изоформ в ячейках матриц сопряженности вариантов сплайсинга в соседних альтернативах этих четырех генов. В принципе, действие тканеспецифических факторов на области соседних 5 - и 3 -альтернатив может объяснить наблюдаемое распределение вариантов сплайсинга. Мы исключили эту возможность, проанализировав информацию о тканях из соответствующих EST клонотек (таб. 3.S1). Для всех генов кроме BRRN1 не наблюдается зависимости между вариантами сплайсинга и типами тканей, из которых были получены последовательности EST. Рисунок 3.2.1 Примеры генов, для которых был обнаружен координированный альтернативный сплайсинг с высоким уровнем достоверности. В каждой альтернативе содержится по два пути a, b в 5 и с, d в З . В таблицах показано число EST в EDAS-I, подтверждающее каждую пару путей в соседних альтернативах. ( ) - Экзоны в 5 нетранслируемой области гена, ( )- экзоны в кодирующей области гена не приводящие к появлению СТОП-кодона, ( ) - экзоны, приводящие к появлению СТОП-кодона, который не может активировать механизм NMD.

Применение алгоритма идентификации соседних альтернатив позволило обнаружить гены, для которых наблюдается координированный сплайсинг пре-мРНК. Из 630 генов, содержащих соседние альтернативы с достаточным для анализа покрытием EST, только для 60 генов был показан координированный сплайсинг. Мы оценили, что около 25% генов человека содержат более одной альтернативно сплайсируемой области, следовательно, сплайсинг этих областей может быть координирован Оценка была получена с помощью программы IsoformCounter [87] (см. следующую главу) В настоящей работе использовались EST, которые содержат в среднем четыре экзона, кроме того, 3 и 5 концы генов, кодирующих длинные транскрипты мРНК, обычно имеют значительно более высокое покрытие по сравнению с серединой. Мы идентифицировали лишь незначительную часть существующих соседних альтернатив для которых EST покрытие было достаточно для анализа. В работе [89] мы рассмотрели пять генов (рис. 3.2.1), содержащих области координированного альтернативного сплайсинга, в которых, совместное распределение вариантов в соседних альтернативах с высокой значимостью отличалось от ожидаемого (pval 0,007). Мы показали, что альтернативный сплайсинг четырех генов (CNOT8, РСВР2, BRRN1, MGAT4B) не может быть объяснен активацией механизма деградации изоформ, содержащих преждевременный СТОП-кодон. Действие тканевых факторов так же не может объяснить паттерны сплайсинга для всех генов кроме BRRN1.

Частота ошибок сплайсинга

Как было показано в главе 3.1, около 59% всех альтернатив, имеющих EST покрытие, приходится на редкие варианты сплайсинга (см. таб. 3.1.2). Для выявления редких альтернатив в главе 3.1 был использован биноминальный тест, предложенный Кап с соавторами. [25] с параметром 0,01. Мы предполагаем, что с некоторой вероятностью сплайсосома способна выбрать в качестве сайта произвольный ДНК мотив, обладающий минимально возможным сходством с реальным сайтом (например, произвольный GT или AG динуклеотиды). Параметр, используемый в биноминальном фильтре, отражает нашу оценку вероятности ошибки сплайсосомы [87] (см. главу 4.1 нормализация на EST покрытие). Здесь будет приведена более точная оценка, по сравнению с опубликованной ранее [87].

Программа IsoformCounter для каждого интрона позволяет определить существование хотя бы одной изоформы, не попадающей под действие механизма NMD и кодирующей, достаточно длинную аминокислотную последовательность (глава 4.1: фильтры (1)-(5)) К сожалению, из факта существования такой изоформы никак не следует, что кодируемый в ней белок способен выполнять какую-либо функцию. Мы полагаем, что если для интрона не существует изоформы, удовлетворяющей фильтрам IsoformCounter, то с высокой долей вероятности, она не может кодировать функциональный белок. Рассмотрим множество (error jntrons) состоящее из интронов, подтвержденных только последовательностями EST через которые IsoformCounter не может провести ни одной изоформы. Кроме того, каждый из интронов, входящих в состав множества error jntrons пересекается с каким-либо интроном, достоверным на уровне белка. Рассмотрим множество последовательностей EST (Good) из EDAS-II, которые сплайсированы только в тех интронах, которые имеют достоверность на уровне белка, обозначим Ng00d - число таких EST. Обозначим Error множество EST, сплайсированных хотя бы в одном из ошибочных интронов из множества error jntrons , при этом эти EST должны быть обязательно сплайсированы хотя бы в одном белковом интроне, пусть Nerror - число таких EST. EST из множества Error соответствуют транскриптам, в которых произошла ошибка при сплайсинге интронов в достоверно кодирующих изоформах. Обозначим а - вероятность ошибки сплайсинга одного интрона. Пусть EST в среднем содержит «/ интронов, тогда (1- а) т - вероятность, что EST не содержит ни одной ошибки сплайсинга. Вероятность ошибки сплайсосомы может быть получена, как решение уравнения: (1-(1- аут )(Ngood + Nermr) = Nenor (см. таблицу 4 2 1)

Оценка вероятности ошибки сплайсинга одного интрона - а. Первая оценка получена без ограничения на число клонотек, подтверждающих EST интроны, входящие в состав множества error jntrons. Вторая оценка получена только с использованием тех EST интронов, которые наблюдались только в одной клонотеке. любое число клонотек 1 EST клонотека Naood (ЧИСЛО EST) 1848958 1848958 Nerror (число EST) 38845 17926 a 0,0069 0,0032

Для того, чтобы понять насколько наша оценка близка к истине, будем рассуждать следующим образом. Рассмотрим ген содержащий пі интронов Если а - вероятность ошибки при сплайсинге одного интрона, то (1- а)" /?, где /? - доля нормальных мРНК необходимых для эффективной экспрессии гена. Мы предполагаем, что ошибка сплайсинга с большой вероятностью вызовет преждевременный СТОП. Оценкой «здравого смысла» можно считать следующую величину a 1- /? " . В таблице 4.2.2 приведено число интронов, которое может содержать ген, так чтобы при вероятности ошибки сплайсинга а обеспечивалась эффективность экспрессии /?.

Несмотря на то, что среднее число интронов в генах человека относительно невелико (5,5) [51], некоторые гены имеют очень большое число интронов. Экстремальный пример - ген титин содержит 363 экзона, длина ORF, проходящей через все экзоны, составляет 114414 пар оснований (4200 кД). Экспериментально была показана коэкспрессия изоформ титана -700 кД и 3000 кД [79]. Из таблицы 4.2 2 следует, что величина ошибки сплайсосомы должна находится в интервале 0,001 - 0,003, чтобы обеспечить достаточную эффективность экспрессии длинной изоформы титана. В главах 3.1 и 4.1 была использована завышенная оценка ошибки сплайсомы 0,01. Однако, ошибки сплайсинга, с высокой вероятностью сохраняют рамку считывания, особенно в случае делеции в базовой изоформе (не происходит вставки СТОП кодонов из интронной последовательности см. рисунок 3.S1). Так как мы рассматривали только те EST интроны, через которые с точки зрения IsoformCounter, не проходит функциональная открытая рамка считывания, то наша оценка ошибки сплайсинга (0,003-0,007) может быть заниженной. Это не противоречит рассуждениям об эффективности сплайсинга - гены с большим числом интронов составляют незначительную долю всех генов, и для них за счет регуляции может быть существенно снижена ошибка сплайсинга. Вероятность ошибки сплайсинга может зависеть от контекста сайтов сплайсинга базовой изоформы и варьировать в достаточно широком диапазоне: 10 3 - 10 2. Представляет интерес сравнение вероятности ошибки сплайсинга с ошибкой РНК-полимеразы и ошибкой трансляции. В отличие от ДНК полимеразы, которая имеет очень низкий уровень ошибок - 1 на 109 нуклеотидов, РНК полимераза, не имеющая корректирующей активности, ошибается с частотой 1 на 104. Ошибки трансляции происходят с частотой 1 на 104 аминокислотных остатков [96]. Длина кодирующей части мРНК среднего гена 900 нуклеотидов, тогда вероятность того, что при транскрипции в ней содержится хотя бы одна ошибка 0,09! Таким образом, около 10% мРНК содержат ошибки в кодирующей части гена, внесенные РНК полимеразой, примерно такой же уровень ошибок вносится при сплайсинге пре-мРНК генов с 10 (верхняя оценка а) - 100 (нижняя оценка) нитронами.

Похожие диссертации на Компьютерный анализ сплайсинга

Анализ соответствия молекулярных и морфологических данных при анализе филогении на примере семейств бобовые (Leguminosae) и зонтичные (Umbelliferae)Дегтярева Галина Викторовна

Антирестрикционный белок ArdA, кодируемый трансмиссивной плазмидой IncI группы: мутационный анализ и изучение механизма действияНекрасов Сергей Вячеславович

Структурно-функциональный анализ генов, кодирующих специфические субъединицы РНК-полимераз II и III делящихся дрожжейПрошкина Галина Михайловна

Анализ регуляции транскрипции генов дыхания в гамма-протеобактериях методами сравнительной геномикиГерасимова Анна Викторовна

Сравнительный анализ и консервация геномов штаммов вируса натуральной оспы из российской коллекцииБабкина Ирина Николаевна

Экспрессия генов белков NS5B и NS5A вируса гепатита C, ингибиторный анализ и белок-белковое взаимодействиеИванов Александр Владимирович

Структура и анализ генома вируса КаршиЛяпина Ольга Викторовна

Разработка молекулярных маркеров для эпидемиологического анализа клинических штаммов Mycobacterium TuberculosisЛихошвай Екатерина Юрьевна

Молекулярно-генетический анализ эндогенных ретровирусов у некоторых пород свиней и диких кабановАйтназаров Руслан Бейшеналиевич

Получение и анализ трансгенных растений Nicotiana tabacum со смысловой и антисмысловой формами агробактериальных генов ipt и iaaMАлексеева Валерия Витальевна