Компьютерный анализ сплайсинга тема диссертации и автореферата по ВАК РФ 03.00.03, кандидат биологических наук Неверов, Алексей Дмитриевич

• Неверов, Алексей Дмитриевич
• кандидат биологических наук
• 2007

Оглавление диссертации кандидат биологических наук Неверов, Алексей Дмитриевич

Список использованных сокращений и обозначений.

Апробация

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

1.1 Введение.

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

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

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

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

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

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

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

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

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

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

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

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

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

Введение диссертации (часть автореферата) на тему «Компьютерный анализ сплайсинга»

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

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

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

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

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

VIstM Mi

AGGU (STnA&GU экзон 1 интрон экзон2

I wiiiiarivn ния

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

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

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

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

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

AG4 (VJoAGGU uWiE;

GU

PJ3 В

GO m

AS|SU JF 0/У

Соединенные экзоны

Алгоритмы поиска генов можно разделить на два вида - 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(xt = s) = 1 и 2) для любой последовательности отметок времени seZ l<i}<i2< .</n</</, любых состояниях s,t е S и любых В/,В2, В„ подмножествах S выполняется P(Xj=t\ х^еВ^хаeB?, ,xmeB„ x,=s)=P(xj=t\x,=s) Свойство (2) называется свойством марковости и означает, что при любом фиксированном положении системы в данный момент времени, будущее поведение системы (j>i) не зависит от поведения системы в прошлом. Марковская цепь может быть представлена в виде набора (S,A,n) состояний S, матрицы вероятностей переходов между состояниями A=fasl=P(x,=t\x, i=s)} и вектором начальных состояний тт=( P(xi=s)), 5=7 N. Если вероятности переходов as, не зависят от времени i, то цепь называется однородной. Вероятность любой последовательности Х= Х\, х,. , xj определяется выражением

1=2

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

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

PiXtZ^x^wfl^e^x,). i=2

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

Обозначим D=ni,ri2, riL, ^nt=T ~ последовательность длин в каждом состоянии.

-1 I

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

P(X,Z,D) = xzezt (ЛЦ К, (И,)ГК №. + Ц - 1К, (И,).

1=2

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

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

Рисунок 1.2.1. Структура НММ эукариотического генома, используемая в программе Genscan [51]. На рисунке используются обозначения: кодирующие состояния (экзоны) Einit - начальный, Eterm - последний экзоны, Е„ 1-0,1,2 внутренние экзоны в соответствующей рамке считывания; Esingie - одноэкзонный ген; некодирующие состояния рассматриваемых алгоритмов, модель Doublescan описывает выравнивание двух последовательностей ДНК, содержащих гомологичные гены. Одноэкзонные гены, а так же первый и последний экзоны гена, содержащего интроны, моделируются отдельными состояниями НММ. Первый экзон гена начинается всегда со СТАРТ кодона (ATG), последний экзон всегда заканчивается одним из СТОП кодонов (ТАА, TAG, TGA). Интроны рассматриваются как вставки в кодирующую последовательность мРНК. Они могут приходиться между соседними кодонами, а так же после первого либо второго нуклеотида внутри кодона. В зависимости от числа нуклеотидов разорванного кодона, которые переносятся на следующий экзон на каждой цепи ДНК, существует по три интронных состояния модели, /о, //, h• Внутренний экзон может транслироваться в одной из трех возможных рамок 1-0,1,2 в зависимости от того, сколько нуклеотидов от начала экзона относятся к кодону мРНК, разорванному предшествующим интроном. Таким образом, для каждого интронного состояния возможен единственный переход в соответствующее состояние внутреннего экзона - /, Е, , i=0,l,2. Обратная цепь ДНК моделируется состояниями соответствующими один к одному состояниям прямой цепи, а переходы осуществляются в обратном 3' 5' направлении.

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

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

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

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

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

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

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

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

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

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

Алгоритм Витерби является алгоритмом поиска оптимальной последовательности скрытых состояний и их длин, оптимизирующий совместную вероятность (1 см. Скрытая Марковская модель эукариотического гена). Пусть o(zk,nk,l) - максимальная вероятность пути из к состояний, который заканчивается в состоянии z* , при этом щ наблюдаемых сгенерированы в этом состоянии. Полная длина сгенерированной к последовательности наблюдаемых X[l,l] : l = . Пусть так же вероятностьv{zk,nk,t) 1 известна для всех состояний zи для всех щ </, тогда легко вычислить вероятность u{zk+l,nk+i,l + nk+]) = P(X[l,l + nk+l])d Jnk+l) max (v(zk,nk,l)a). zkeb l+nktl<.l

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

1) инициализация (/=0): u(z,,0,0) = я- , z/eS;

2) рекурсия l=T 1, Пк<1\ u{zk,nk,l) = P(X[l-nk +\,l])d2k(nk)maxZk ^ ^Пк(р{гкх,пкА,1-пк)а1к Л) ■ ptr{z *k,n*k,l) = argmax2t ^ ({v{zkx,nkx,l-nk)a2^2i)•

3) завершение: P(X, Z*,D*) = maxZj щ (u(zk, nk, L)); ptr(z *k,n*k>L) = arg max Zt „t (u(zk, nk, L));

4) обратный ход (l=T 1): (z\ !,nk.,,l-nk)=ptr (z\,n\l) позволяет восстановить оптимальные последовательности состояний Z* и длин D*.

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

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

ZD) того чтобы вычислить величину Р(Х), введем сначала «прямые» переменные a(l,z) = P(X[\,l],z), которые являются вероятностью того, что последовательность наблюдаемых до координаты /, сгенерирована так, что последний символ сгенерирован в состоянии zeS. Запишем рекурсивный алгоритм, подобный алгоритму Витерби, для вычисления «прямых» переменных:

1) инициализация (1=1): а(\,д) = жд, ceS,

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

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

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

P(X,X[aMz) = P(X[ha-\],z)P(X[aMz)d2(b-a + \)P(X[b + \,T]\X[\.b],z) = Р(Х[ 1,а-1],z)dz(b-a +1 )Р(Х[а,Ь] | z)P(X[b +1],z) = а(а,z)dz(b-a +1 )P(X[a,b] \ z)p{b,z) Вторая строка выражения является следствием условия марковости - все, что сгенерировано после момента Ъ зависит только от состояния цепи в этот момент Так же было введено обозначение для «обратных» переменных p(l,z). Алгоритм вычисления обратных» переменных полностью аналогичен алгоритму вычисления «прямых» переменных:

1) инициализация (1=Т): Р{Т,д) = 1, ceS;

2) рекурсия(№1,1): /?(/,*) = £ ^а1дР{Х[1 + \,1 + п}\дКдШ{1 + п,дУ,

Л-1 f€S f *Z

3) завершение: Р(Х) = . eS

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

14= a b\X)- a{a'2)d'{Ь~а + 1)Р{Х[а'Ь] 1

Р{Х)

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

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

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

Лемма 1.

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

Ир)

Z VV<M,)fIw(MM,M,)w(M() = l

Vp-(u,u2, иНр)) 1-2

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

Hp)

Рг (р) = н<5)[] w(uiA, и, Щи,) 1=2

Действительно, пусть на каждой вершине заданы «обратные» переменные являющиеся суммой весов всех путей из вершины и, в сток. Пусть утверждение верно для максимального пути из источника в сток длиной п ребер. Если n = 1 утверждение очевидно. Докажем для случая, когда длина максимального пути равна п+1 ребер. Рассмотрим вершины и/, .щ связанные с источником s соответствующими ребрами е/, ,вк Выберем ребро i с вероятностью w(s)w(e,)(il. Обозначим о вес пути из вершины и, в сток t. Расстояние до стока из выбранной вершины не больше п ребер. Выберем путь в сток из вершины и, с вероятностью а/Д Вес пути из источника в сток через выбранную вершину равен w(s)aw(et), и этот путь был выбран с вероятностью w(s)w(el)fil cr/f], = w(s)ow(eJ.

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

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

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

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

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

Качество предсказания оценивается на основании множества последовательностей ДНК с известной аннотацией - тестового множества. Здесь мы не будем обсуждать подробно проблемы построения тестового множества, отметим, что чаще всего в него включаются фрагменты ДНК, содержащие гены, подтвержденные экспериментальной процедурой или маркером экспрессии (мРНК, 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[5и+ + Sp+ + Sn~ + Sp~] -1 (approximate correlation),

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

РР и АР)

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

Меры, оценивающие качество предсказания на экзонном и геномном уровнях более строгие, чем нуклеотидные. В качестве правильного предсказания на экзонном уровне (TP) рассматриваются экзоны, сайты которых точно совпадают с экспериментальной аннотацией. На уровне генома правильными предсказаниями считаются предсказания генов, все сайты которых совпадают с сайтами кодирующих экзонов аннотации, включая СТАРТ кодон. Качество предсказания оценивается чувствительностью и специфичностью. Эти две меры вычисляются аналогичным образом - также как на уровне нуклеотидов, с той разницей, что учитываются предсказания экзонов или генов. На экзонном уровне используются полезные меры перепредсказания - доля предсказанных экзонов, которые не пересекаются в соответствующей рамке ни с одним из экзонов аннотации 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*46=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 состава, были зафиксированы на уровне, характерном для генома Е coll. Частоты кодонов могут быть восстановлены по известному GC% составу в соответствии с выражением: где ЩП]П2 обозначает один из кодонов аминокислоты А, /а - частота аминокислоты А (зависит от GC%), fi - частоты кодонов, полученные на основании позиционных частот соответствующих нуклеотидов (зависит от GC%),/r - частота кодона, скорректированная на частоту аминокислоты А и частоты синонимичных кодонов fi. Существенными и0и,и2) достоинствами эвристической модели являются: малое число параметров (частоты трех типов нуклеотидов), а, следовательно, малый размер последовательности ДНК, необходимой для обучения (>400 нук.), применимость к геномам с неоднородным распределением GC состава. Малый объем обучающего множества позволил с успехом применять эту модель для распознавания генов в геномах вирусов ВИЧ-1 и Т-клеточного лимфотропного вируса тип I.

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

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

Рассмотрим формально задачу выбора значений параметров для скрытой марковской модели. Обозначим набор параметров НММ - в Параметрами марковской цепи являются вероятности переходов (aid), вероятности эмиссии ek(b)=P(x=b\z=k) и распределение длин обобщенных состояний. Будем предполагать, что имеется множество из п последовательностей ДНК у=(у', сгенерированных независимо одной марковской цепью, причем соответствующие последовательности состояний неизвестны (скрыты) Будем максимизировать log правдоподобия по всем возможным значениям параметров в maхД^ДУ m = max*\-b°zp(yJ m j-1

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

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

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

Аи Ек(Ь) выражениями: аы = и ек = * . lalAkl lab

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

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

DfYr ,п , ,\v т a(a,k)dk(b-a +1 )Р(Х[а,Ь]\к)ашр(Ъ +1,/) состояние zJ+i Р(Х[а,Ь],г; = k,zJ+l =1\Х,в) =-*-——-=-,

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

Лш =1-БГТТ Та'(а>к)анПу'[а,Ь]\к)с1к(Ь-а + 1)^(Ь + 1,1),

J Р\У )y>[aJ>W где верхний индекс «прямых» и «обратных» переменных обозначает номер последовательности в совокупности наблюдаемых у, выражение У/"а,Ь] обозначает фрагмент последовательности j. Аналогично вычисляется ожидаемое число эмиссий символа b в состоянии к:

Еь (*)=Е^гтт I 1У w & w* (b~a+vp(yJ № I*), j )yJ[a,b]eyJie[a,b]yJ[i]=b где внутренняя сумма берется только по тем позициям /, в которых генерируется символ - Ь. Ожидаемое распределение длин в состоянии к так же является параметром модели и может быть вычислено, как сумма вероятностей сгенерировать п=Ь-а+1 символов для «=7,2 Так как в данном итерационном процессе мы движемся к максимуму в непрерывном пространстве, мы никогда не достигнем максимума Поэтому, необходимо установить критерий сходимости, который останавливает процесс, если изменения логарифма правдоподобия совокупности данных достаточно малы.

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

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

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

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

Ъ) p{z\,. где Т и L обозначают длины соответствующих наблюдаемых последовательностей и последовательностей состояний. Формулу (3) можно получить, проинтегрировав формулу (1) по априорному распределению параметров р(9), тем самым выбор последовательностей скрытых состояний не будет зависеть от первоначального выбора параметров НММ. Можно показать, что распределение (3) пропорционально правдоподобию разметки i для последовательности У' p(z J = к I yj,z[-jl , y[-j]) oc Yl h(y[-j])Hyn i

Другими словами, алгоритм состоит из следующих шагов- (А) последовательность j изымается из совокупности у, (Б) пересчитываются вхождения свободных переменных V ied- h()/j]J, (В) вычисляются МП оценки параметров модели V \евш основании (Г) для исключенной последовательности случайно выбирается разметка z?, затем последовательность j возвращается в совокупность у Шаги А-Г повторяются рекурсивно, пока не будет выполнен критерий сходимости Напомним, что /, обозначает совокупность фрагментов последовательности у, зависящих от разметки г*, на которой вычисляется число вхождений свободного параметра / ев Свободными параметрами переменными) являются независимые параметры модели. Предсказание генов.

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

P(e) = fja'!P(e\Ml)

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

Заключение диссертации по теме «Молекулярная биология», Неверов, Алексей Дмитриевич

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

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

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

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

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

4. С помощью разработанного метода выявления альтернатив было показано, что для ряда генов выбор варианта сплайсинга 3'-альтернативы зависит от способа сплайсинга 5'-альтернативы. Была дана оценка доли таких генов - не более 25%.

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

6. Дана оценка вероятности ошибки сплайсинга - 10"2 -10'3 случаев ошибочного сплайсинга на интрон.

Список литературы диссертационного исследования кандидат биологических наук Неверов, Алексей Дмитриевич, 2007 год

1. Eyras Е, Caccamo М, Curwen V, Clamp M: ESTGenes: Alternative splicing from ESTs in Ensembl. Genome Res 2004,14(5).976-987.

2. Mironov AA, Fickett JW, Gelfand MS: Frequent alternative splicing of human genes Genome Res 1999,9:1288-1293.

3. International Human Genome Sequencing Consortium: Finishing the euchromatic sequence of the human genome. Nature 2004,431:931-45.

4. Southan C: Has the yo-yo stopped? An assessment of human protein-coding gene number. Proteomics 2004,4:1712-1726.

5. Johnson JM, Castle J, Garrett-Engele P, Kan Z, Loerch PM, Armour CD, Santos R, Schadt EE, Stoughton R, Shoemaker DD: Genomewide survey of human alternative pre-mRNA splicing with exon junction microarrays. Science 2003,302:2141-2144.

6. Kampa D, Cheng J, Kapranov P, et al: Novel RNAs identified from an in-depth analysis of the transcriptome of human chromosomes 21 and 22. Genome Res 2004,14(3):331-342.

7. Xing Y, Lee C: Alternative splicing and RNA selection pressure evolutionary consequences for eukaryotic genomes. Nat Rev Genet. 2006 Jul;7(7):499-509

8. Brett D, Pospisil H, Valcarcel J, Reich J, Bork P: Alternative splicing and genome complexity. Nature Genet 2002,30:29-30.

9. Foissac S, Schiex T: Integrating alternative splicing detection into gene prediction. BMC Bionformatics 2005,10:6-25.

10. Wang P, Yan B, Guo J, Hicks C, Xu Y: Structural genomics analysis of alternative splicing and application to isoform structure madelling. PNAS 2005.

11. Lee C, Wang Q: Bioinformatics analysis of alternative splicing. Brief Bioinform 2005, 6(1): 23-33

12. Garsia J, Gerber S, Sugita S et al: A conformational switch in the Piccolo C2A domain regulated by alternative splicing. Nat. Struct. Mol. Biol. 2004,11(1): 45-53.

13. Offman M, Nurtdinov R, Gelfand M, Frishman D: No statistical support for correlation between the positions of protein interaction sites and alternative spliced regions. BMC Bioinformatics 2004,5(1): 41.

14. Resch A, Xing Y, Modrek В et al: Assessing the impact of alternative splicing on domain interaction in human proteome J. Proteome Res. 2004,3(1): 76-83.

15. Xing Y, Xu Q, Lee C: Widespread production of novel soluble protein isoforms by alternative splicing removal of transmembrane anchoring domains. FEBS Lett. 2003, 555(3): 572-578.

16. Pan Q, Shai O, Misquitta C, et al: Revealing global regulatory features of mammalian alternative splicing using a quantitative microarray platform. Mol. Cell 2004,16(6): 929-941.

17. Sorek R, Ast G, Graur D: Alu-containing exons are alternativly spliced. Genome Res. 2002, 12: 1060-1067.

18. Singer S, Mannel D et al: From "junk" to gene : curriculum vitae of a primate receptor isoform gene. J. Mol. Biol. 2004,341: 883-886.

19. Kondrashov F, Koonin E: Origin of alternative splicing by tandem exon duplication. Hum. Mol. Genet. 2001,10: 2661-2669.

20. Modrek B, Lee C: Alternative splicing in the human, mouse and rat genomes is associated with an increased rate of exon creation/loss. Nature. Genet. 2003,34: 177-180.

21. Cusak B, Wolfe K: Changes in alternative splicing of human and mouse genes are accompanied by faster evolution of constitutive exons. Mol. Biol. Evol. 2005,22: 2198-2208.

22. Wang W et al: Origin and evolution of new exons in rodents. Genome Res 2005, 15: 12581264.

23. Xing Y, Lee C: Negative selection pressure against premature protein truncation is reduced by alternative splicing and diploidy. Trends Genet. 2004,20:472-475.

24. Sugnet C, Kent W, Ares M, Haussler D: Transcriptome and genome conservation of alternative splicing events in human and mice. Рас. Symp Biocomput. 2004,66-77.

25. Kan Z, States D, Gish W: Selecting for Functional Alternative Splices in ESTs Genome Res 2002 12: 1837-1845

26. Lewis B, Green R, Brenner S: Evidence for widespread coupling of alternative splicing and nonsense-mediated mRNA decay in humans. Proc. Natl. Acad. Sci. USA 2003,100(1): 189-192.

27. Hillman R, Green R, Brenner S: An appreciated role for RNA surveillance. Genome Biol 2004,5(2): R8.

28. Cuccurese M, Russo G, Russo A, Pietropaolo C: Alternative splicing and nonsense-mediated mRNA decay regulate mammalian ribosomal gene expression. Nucleic Acids Research 2005, 33(18): 5965-5977.

29. Pan Q, Saltzman A, Kim Y et al: Quantitative microarray profiling provides evidence against widespread coupling of alternative splicing with nonsense-mediated mRNA decay to control gene expression. Genes & Dev 2006 20: 153-158

30. Kalyanaraman A, Aluru S: Expressed sequence tags: Clustering and applications. Handbook on Computational Molecular Biology.2006 S.Alluru, ed. CRC Press:

31. Malde K, Coward E, Jonassen I. A graph based algorithm for generating EST consensus sequences. Bioinformatics 2005,21(8):1371-5.

32. Lee C, Grasso C, Sharlow M: Multiple sequence alignment using partial order graphs. Bioinformatics 2002,18(3): 452-464.

33. Lee C: Generating consensus sequences from partial order multiple sequence alignment graphs. Bioinformatics 2003,19(8): 999-1008.

34. Heber S, Alexeyev M, Sze S et al: Splicing graphs and EST assembly problem. Bioinformatics 2002,18 Suppl. 1: S181-188.

35. Eyras E, Caccamo M, Curwen V, Clamp M: ESTGenes: alternative splicing from ESTs in Ensembl. Genome Res 2004,14(5). 976-87.

36. Leipzig J, Pevzner P, Heber S: The alternative splicing gallery (ASG): bridging the gap between genome and transcriptome. Nucleic Acids Res. 2004,32(13): 3977-3983.

37. Xing Y, Resch A, Lee C: The multiassembly Problem: Reconstructing multiple transcript isoforms from EST fragment mixtures. Genome Res. 2004,14. 426-441.

38. Xing Y, Yu T, Wu YN, Roy M, Kim J, Lee C: An expectation-maximization algorithm for probabilistic reconstructions of full-length isoforms from splice graphs. Nucleic Acids Res. 2006,34(10):3150-60.

39. Usuka J, Zhu W, Brendel V: Optimal spliced alignment of homologues cDNA to a genomic DNA template. Bioinformatics 2000,16:200-211.

40. Brendel V, Xing L, Zhu W: Gene structure prediction from consensus spliced alignment of multiple ESTs matching the same genomic locus. Bioinformatics 2004,20:1157-1164.

41. Foissac S, Schiex T: Integrating alternative splicing detection into gene prediction. BMC Bioinformatics 2005,6:25.

42. Ohler U, Shomron N, Burge C: Recognition of unknown conserved alternatively spliced exons. PLoS Comp Biol 2005,1(2): e5.

43. Mather C, Sagot M, Schiex T,Rouze P: Current methods of gene prediction, their strength and weaknesses. Nucl. Acids Res. 2002,30:4103-4117.

44. Pertea M, Salzberg S: Computational gene finders in plants. Plant Mol. Biol. 2002,48:39-48.

45. Neverov A, Mironov A, Gelfand M: Splice alignment and similarity-based gene recognition. Handbook of computational molecular biology, Aluru S, Chapman & Hall 2006,2:1-18.

46. Do J, Choi D: Computational approaches to gene prediction. The J. of Microbiol. 2006, 4:137-144.

47. Burge C, Karlin S: Prediction of Complete Gene Structures in Human Genomic DNA. J. Mol. Biol. 1997,268:78-94.

48. Kulp D, Haussler D, Reese M, Eeckman F: A generalized Hidden Markov Model for the recognition of human genes in DNA. Intell. Sys. for Mol. Biol., 4:134-142.

49. Salamov A, Solovyev V: Ab initio gene finding in Drosophila genomic DNA. Genome Res. 2000,10:516-522.

50. Lukashin A, Borodovsky M: GeneMark.hmm: new solution for gene finding. Nucl. Acids Res 1998,26(4): 1107-1115.

51. Yao H, Guo L, et al: Evaluation of five ab initio gene prediction programs for the discovery of maize genes. Plant Mol. Biol. 2005, 57:445-460.

52. Чистяков В. П: Курс теории вероятностей. "Агар", Москва 2000,161-174.

53. Durbin R, Eddy SR, Krogh A, G. Mitchison: Biological sequence analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press 2002.

54. Mayer I, Durbin R: Comparative ab initio prediction of gene structures using pair HMMs. Bioinformatics 2002,18:1309-1318.

55. Fickett J, Tung C: Assessment of protein coding measures. Nucl. Acids Res. 1992, 20:64416450.

56. Kriventseva E, Gelfand M: Statistical analysis of the exon-intron structure of higher and lower eukaryote genes. J.of Biomol. Struct, and Dynamics 1999,17:281-288.

57. Guigo, R., Agarwal, P., et al: An assessment of gene prediction accuracy in large DNA sequences. Genome Res. 2000,10: 1631-1642.

58. Lomsadze A, Ter-Hovhannisyan V, Chernoff Y, Borodovsky M: Gene identification in novel eukaryotic genomes by self-training algorithm. Nucl. Acids Res. 2005,33(20): 6494-6506.

59. Mitrofanov A, Lomsadze A, Borodovsky M: Sensitivity of hidden Markov models. J. Appl. Probab. 2005,42:632-642.

60. Besemer J, Borodovsky M: Heuristic approach to deriving models for gene finding. Nucl. Acids Res. 1999,27(19):3911-3920.

61. Chatteiji S, Pachter L: Large multiple organism gene finding by collapsed Gibbs sampling", J. Сотр. Biol 2005,12(6):599-608.

62. Besemer J, Lomsadze A, Borodovsky M: GeneMarkS: a self-training method for prediction of gene starts in microbial genomes, implications for finding sequence motifs in regulatory regions. Nucleic Acids Res. 2001,29: 2607-2618.

63. LarsenT, Krogh A: EasyGene a prokaryotic gene finder that ranks ORFs by statistical significance. BMC Bioinformatics 2003,4: 21.

64. Altschul, Stephen F, Madden T, Schaffer A, Zhang J, Zhang Z, Miller W, Lipman J: Gapped BLAST and PSI-BLAST: a new generation of protein database search programs, Nucleic Acids Res. 1997,25:3389-3402.

65. Boue S, Vingron M, Kriventseva E, Koch I: Theoretical analysis of alternative splice forms using computational methods. Bioinformatics 2002,18(Suppl. 2): s65-s73.

66. Webster N, Evans L, Caples M, Erker L, Chew S: Assembly of splicing complexes on exon 11 of the human insulin receptor gene does not correlate with splicing efficiency in-vitro. BMC Mol Biol. 2004,5:7.

67. Cote J, Dupuis S, Wu J: Polypyrimidine track-binding protein binding downstream of caspase-2 alternative exon 9 represses its inclusion. J Biol Chem. 2001,276(11):8535-43.

68. Kozak M: Pushing the limits of the scanning mechanism for initiation of translation. Gene 2002,299:1-34.

69. Pestova TV, Kolupaeva VG, Lomakin IB, Pilipenko EV, Shatsky IN, Agol VI, Hellen CU: Molecular mechanisms of translation initiation in eukaryotes. Proc Natl Acad Sci USA 2001, 98:7029-7036.

70. Kan Z, Rouchka EC, Gish WR, States DJ: Gene structure prediction and alternative splicing analysis using genomically aligned ESTs. Genome Res 2001,11:889-900.

71. Zavolan M, van Nimwegen E, Gaasterland T: Splice variation in mouse full-length cDNAs identified by mapping to the mouse genome. Genome Res 2003,12:1377-1385.

72. Resch A, Xing Y, Modrek B, Gorlick M, Riley R, Lee C: Assessing the impact of alternative splicing on domain interactions in the human proteome. J Proteome Res 2004, 3:76-83.

73. A. D. Neverov, M. S. Gelfand, A. A. Mironov 2003. GipsyGene: A Statistics-Based Gene

74. Recognizer for Fungal Genomes. Biophysics (Moscow), Vol 48, Suppl 1, 2003, pp S71-S75

75. Mironov A, Novichkov P, Gelfand M:. Pro-Frame: similarity-based gene recognition in eukaryotic DNA sequences with errors. Bioinformatics. 2001,17(l):13-5.

76. Neverov AD, Gelfand MS, Mironov AA: Gene prediction in genomics DNA of Aspergillus, Third International conference on bioinformatics of genome regulation and structure. (BGRS'2002), 1:116, Новосибирск, Россия.

77. Kupfer D, Drabenstot S, et al: Introns and splicing elements of five diverse fungi, Euk. Cell 2004,3(5):1088-1100.

78. Clark F, Thanaraj T: Categorization and characterization of transcript-confirmed constitutively and alternatively spliced introns and exons from human. Hum Mol Genet. 2002, ll(4):451-64.

79. Neverov AD, Milanesi L: A pipeline for computational gene recognition in the Sacharopolyspora erythraea genome. Moscow Conference on Computational Molecular Biology (MCCMB'05), Москва, Россия, 1:246.

80. Neverov A, Artamonova I, Nurtdinov R, Frishman D, Gelfand M, Mironov A: Alternative splicing and protein function. BMC Bioinformatics 2005,6:266

81. Menta C, Patel, N: Algorithm 643: FEXACT: A FORTRAN Subroutine for Fisher's Exact Test on Unordered Contingency Tables. ACM Trans Math Softw 1986, 12(2):154-161.

82. Fededa J, Petrillo E, Gelfand M, Neverov A, Kadener S, Nogues G, Pelisch F, Baralle F, Muro A, Kornblihtt A: A polar mechanism coordinates different regions of alternative splicing within a single gene. Mol Cell. 2005,19(3):393-404.

83. Price D. P-TEFb, a cyclin-dependent kinase controlling elongation by RNA polymerase II. Mol. Cell. Biol. 2000,20:2629-2634.

84. Ramensky V, Nurtdinov R, Neverov A, Mironov A, Gelfand M: 5th Int. Conf. "Bioinformatics of Genome Regulation and Structure" BGRS'2006 2006, 1:211.

85. Нуртдинов PH, Неверов АД, Малько ДБ, Космодемьянский ИА, Ермакова ЕО, Раменский BE, Миронов АА и Гельфанд МС. (2006): EDAS, база данных альтернативно сплайсируемых генов человека. Биофизика, 2006,51(4): 589-592.

86. Ramil N Nurtdinov, Ilya Kosmodemyansky: The EDAS (EST-derived alternative splicing) database. Moscow Conference on Computational Molecular Biology (MCCMB'03), 1:171, Москва, Россия, 2003.

87. Neverov AD: GipsyGene: A HMM-based gene recognitional tool for lower fungi. Moscow Conference on Computational Molecular Biology (MCCMB'03), 161, Москва, Россия, 2003.

88. Nurtdinov RN, NeveroM AD, Malko DB, Kosmodemyansky IA, Mironov AA, Gelfand MS: EDAS: EST-derived alternative splicing database, Moscow Conference on Computational Molecular Biology (MCCMB'05), Москва, Россия, 1:259.

89. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molecular biology of the cell 4th edition, Garland Publishing, Inc. 2002, New York & London.

90. Д. Гасфилд: Строки, деревья и последовательности в алгоритмах: Информатика и вычислительная биология, СПб., Невский диалект, Петербург 2003.

91. Florea L, Hartzell G, Zhang Z, Rubin GM, Miller W: A computer program for aligning a cDNA sequence with a genomic DNA sequence. Genome Res 1998,8(9): 967-974.

92. Kent WJ: BLAT-the BLAST-like alignment tool. Genome Res 2002,12(4): 656-64.

93. Ogasawara J, Morishita S: Fast and sensitive algorithm for aligning ESTs to humangenome. Proc IEEE Comput Soc Bioinform Conf 2002,1: 43-53.

94. Alexandersson M, Cawley S, Pachter L: SLAM: cross-species gene finding and alignment with a generalized pair hidden Markov model. Genome Res. 2003,13(3):496-502.

95. Thomas JW, Touchman JW, Blakesley RW, Bouffard GG, Beckstrom-Sternberg SM, Margulies EH, Blanchette M, Siepel AC, Thomas PJ, McDowell JC: Comparative analyses of multi-species sequences from targeted genomic regions. Nature 2003,424: 788-793.

96. Imanishi T, Itoh T, Suzuki Y, O'Donovan C, Fukuchi S, Koyanagi КО, Barrero RA, Tamura T, Yamaguchi-Kabata Y, Tanino M, et al: Integrativeannotation of 21,037 human genes validated by fulllength cDNA clones. PloS Biology 2004,2:1-20.

97. Gelfand MS: Computational analysis of alternative splicing. Handbook of Computational Molecular Biology (Chapman & Hall/CRC). Ed. Alluru S. 2006,16-1-16-33.

98. Gatfield D, Izaurralde E: Nonsense-mediated messenger RNA decay is initiated by endonucleolytic cleavage in Drosophila. Nature, 2004,429(6991): 575-578.

99. De Hoon MJ, Imoto S, Kobayashi K, Ogasawara N, Miyano S: Predicting the operon structure of Bacillus subtilis using operon length intergene distance, and gene expression information Рас. Symp. Biocomput. 2004; :276-278.