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

Компьютерный анализ сплайсинга

ВАК РФ 03.00.03, Молекулярная биология

Автореферат диссертации по теме "Компьютерный анализ сплайсинга"

На правах рукописи

Неверов Алексей Дмитриевич

компьютерный анализ сплаисинга

03.00.03 - Молекулярная биология

АВТОРЕФЕРАТ Диссертации на соискание ученой степени кандидата биологических наук

- Москва 2007 -

003054430

Работа выполнена в лаборатории биоинформатики Государственного научно-исследовательского института генетики и селекции промышленных микроорганизмов ФГУП "ГосНИИ генетика".

Научный руководитель:

доктор биологических наук, кандидат физико-математических наук Миронов Андрей Александрович Официальные оппоненты:

Доктор биологических наук КОрягина A.C. Кандидат биологических наук Боринская С.А

Ведущая организация: Институт Математических Проблем Биологии РАН (ИМПБ РАН)

Защита диссертации состоится «06» марта 2007 г. в 14-00 на заседании Диссертационного совета Д 217.013.01 при Государственном научно-исследовательском институте генетики и селекции промышленных микроорганизмов по адресу: 117545, г. Москва, 1-й Дорожный проезд д., 1

С диссертацией можно ознакомиться в библиотеке "ГосНИИ генетика".

Автореферат разослан "6" февраля 2007 г.

Ученый секретарь Диссертационного совета,

Кандидат биологических наук

Заиграева Г.Г.

Общая характеристика работы

Актуальность темы

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

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

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

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

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

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

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

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

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

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

• Проверка гипотезы о независимости сплайсинга нитронов в пре-мРНК.

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

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

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

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

Научная новизна и практическое значение

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

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

• С помощью оригинальной модели сайта ветвления в программе GipsyGene было улучшено распознавание коротких интронов в геномах грибов: Aspergillus spp. и Neurospora crassa.

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

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

• Была оценена доля генов, содержащих соседние альтернативы, в которых сплайсинг в 3'-альтернативе зависит от варианта сплайсинга в 5'-альтернативе.

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

• С помощью IsoformCounter, было показано, что альтернативно сплайсируемые гены из категорий "рибосома" и "передача сигналов посредством малых ГТФаз" имеют меньше изоформ, а гены из категории "репликация ДНК и хромосомный цикл" больше изоформ, чем в среднем по всем генам. Было показано, что среди генов, кодирующих белки, участвующие в образовании взаимодействий с другими белками, больше альтернативно сплайсируемых генов, чем в белках, не участвующих в таких взаимодействиях.

• Была оценена вероятность ошибки при сплайсинге одного интрона. Было проведено сравнение вероятности ошибки сплайсинга с вероятностями ошибок других биологических процессов - синтеза РНК и трансляции белка.

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

Апробация работы

Результаты работы представлены на международных конференциях: "Third International conference on bioinformatics of genome regulation and structure (BGRS'2002, Новосибирск, 2002);

"Moscow Conference on Computational Molecular Biology" (МССМВ'ОЗ, Москва, 2003); "Moscow Conference on Computational Molecular Biology" (MCCMB'05, Москва, 2005); "5lh Int. Conf. Bioinformatics of Genome Regulation and Structure" (BGRS'2006, Новосибирск, 2006). По материалам диссертации опубликовано 10 печатных работ.

Объем и структура диссертации

Диссертация изложена на 138 страницах и состоит из 6 глав. В главах 2-4 представлены оригинальные результаты. Список литературы, приведенный в конце диссертации, содержит 106 наименований. Работа содержит 19 рисунков и 7 таблиц. Глава 1 содержит введение и обзор литературы.

Глава 2 посвящена разработке статистической программы распознавания генов в геномах низших эукариот. Была предложена статистическая модель сайта ветвления значительно улучшающая качество предсказания коротких нитронов в геномах низших грибов Aspergillus spp. и Neurospora crassa.

Глава 3 посвящена разработке методов анализа альтернативного сплайсинга, показанного выравниванием EST с геномной последовательностью. В этой главе был предложен алгоритм выявления и определения типов областей альтернативного сплайсинга пре-мРНК. Было показано, что пре-мРНК некоторых генов содержат несколько областей альтернативного сплайсинга. При этом сплайсинг в альтернативной 3' области, синтезируемой позже, зависит от варианта сплайсинга ранее синтезированной части пре-мРНК.

Глава 4 содержит описание разработанного алгоритма сборки EST в полноразмерные

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

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

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

Глава 5 описывает материалы и методы.

Глава б содержит сводку основных результаты и выводы.

Содержание работы

Глава 2. Предсказание генов в геномах низших эукариот Низшие эукариоты имеют большое значение для науки, как модельные микроорганизмы, для медицины и сельского хозяйства, как важные патогены человека, животных и растений, а также как биотихнологические продуценты антибиотиков и биологически-активных веществ. К настоящему времени были секвенированы геномы некоторых низших грибов, в частности, Aspergillus niclulans, Neurospora crassa, hiagnaporte grizea, всего около 20 наименований [Fungal Genetics Stock Center]. В будущем интерес к секвенированию геномов грибов вряд ли ослабеет. Анализ нового генома всегда начинается с автоматической аннотации, целью которой служит предсказание генов для последующего детального анализа. Программы обнаруживающие гены на основе информации о сходстве геномной последовательности и последовательностей белков, мРНК, EST позволяют идентифицировать 50-70% генов. Статистические программы предсказания генов, в основе которых лежит скрытая марковская модель (НММ), как правило, нуждаются в новом обучении для каждого нового генома. Применение таких программ является обязательным этапом аннотации, так как позволяет идентифицировать гены, специфичные для организма. Обучающее множество, необходимое для оценивания параметров НММ, может быть построено на основании генов, предсказанных по сходству с белками. К моменту начала работы над проектом существовала потребность в программе предсказания генов, которая могла бы использовать широкий набор статистических моделей, подбираемых под каждый аннотируемый геном, и могла бы легко переобучаться для аннотации новых геномов.

GipsyGcne статистическая программа распознавания генов

Программа GipsyGene может использовать модель как эукариотического, так и прокариотического генома. Модель генома прокариот может рассматриваться как ограничение эукариотической модели. Для решения задачи распознавания генов строится граф, вершинами которого являются кандидаты в сайты сплайсинга, СТАРТ и СТОП кодоны. Ребра этого графа

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

В модель заложены следующие кодирующие и некодирующие состояния: одноэкзонный ген; начальный, внутренний и терминальный экзоны; нитрон; межгенный участок (спейсер). В ^рвуСепе реализованы наиболее широко используемые модели донорпого и акцепторного сайтов. Статистическая модель выбирается в зависимости от того, имеет ли сайт значимые корреляции между позициями. Если анализ сайтов в обучающей выборке не обнаруживает корреляций между позициями, или эти корреляции статистически не значимы, например, из-за небольшого объема обучающего множества, то вероятность вычисляется по профилю сайта. Другие две модели учитывают корреляции между позициями сайта. (1) \УАМ модель - неоднородная Марковская модель первого порядка, учитывающая корреляции между соседними нуклеотидами. (2) МОО модель - применяется для сайтов, которые имеет значимые корреляции, как с соседними, так и с удаленными позициями. Вероятность вычисляется по профилю в зависимости от оснований, стоящих в позициях, которые наиболее сильно взаимодействуют с другими позициями сайта.

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

Для улучшения предсказаний З'-концов нитронов была разработана модель сайта ветвления. Потенциальный сайт ветвления оценивается по профилю, распределению расстояний до акцепторного сайта, распределению количества динуклеотидов Ав между сайтом ветвления и

акцептором. Профиль и распределения конструируются на множестве слов, наиболее похожих на сайт ветвления Saccharomyces.

Для того чтобы построить профиль сайта ветвления, в окне [-/, -/*] ([-2,-40] для Aspergillus spp. и Neurospora crassa) левее акцепторного сайта осуществлялся поиск наилучшего слова, удовлетворяющего консенсусу [CT]T[AG]a[CT], идеально - СТАаС. Для кандидатов допускались отклонения от идеального слова только в двух из трех вырожденных позиций. На множестве кандидатов строилось распределение расстояний до акцепторного сайта. Применение этой же процедуры к случайным последовательностям задает уровень шума - среднюю частоту встречаемости слова в случайной последовательности. Далее, окно поиска уточнялось так, чтобы частота в окне превышала уровень шума. На множестве кандидатов, найденных в уточненном окне поиска, строилась частотная матрица размером 7 позиций, распределение расстояний до акцептора и распределение динуклеотидов AG. К каждому элементу полученной частотной матрицы добавлялось число, пропорциональное квадратному корню из числа последовательностей - псевдоотсчеты. В состоянии НММ "5' сайт экзона" генерируются нуклеотиды: собственно акцепторного сайта сплайсинга и соответствующего сайта ветвления. Сайт ветвления определяется как мотив с наибольшим весом в окне [-1.3/,-г] от каждого динуклеотида AG, где / и г - левая и правая координаты окна в котором сайт ветвления встречается чаще, чем в случайной последовательности. Вес мотива определяется как логарифм отношения правдоподобия модели

, ,P{x\BPS)P{d)P(N

сайта ветвления к интроннои модели: WBrs = log(-1-. Были использованы

Р{х | Iniron)

следующие обозначения: * - мотив, P(x\BPS) и P(x\lntron) вероятности мотива по профилю сайта ветвления и модели внутренней части интрона соответственно, d и Nac, - количество нуклеотидов и число динуклеотидов AG между сайтом ветвления и кандидатом в акцепторный сайт соответственно. Результаты тестировании

Программа была обучена и протестирована на генах Aspergillus spp. и Neurospora crassa на одногенных и мультигенных фрагментах. Качество предсказания GipsyGene к кодирующим

нуклеотидам в одногенных фрагментах составило 96% чувствительность и 95% специфичность для Aspergillus spp, 93% и 95% соответственно для Neurospora crassa. Качество предсказания экзонов в одногенных фрагментах Aspergillus spp. (и Neurospora crassa) составило: чувствительность - 80 (75)%, специфичность - 73 (69)%, потерянных экзонов - 4 (8)%, полностью ложных экзонов - 12 (12)%, что соответствует характеристикам аналогичных программ. Вклад модели сайта ветвления в качество предсказания незначителен для статистик кодирующих нуклеотидов (-1%), однако, увеличение чувствительности и специфичности распознавания экзонов может достигать 10%, что объясняется улучшением распознавания коротких нитронов.

Как и у всех программ этого класса, применение GipsyGene на многогенных фрагментах приводит к потере в специфичности, в зависимости от длины межгенных спейсеров, что объясняется тем, что в настоящее время не существует надежных моделей, позволяющих предсказывать границы генов. При тестировании программы на искусственных многогенных фрагментах ДНК, со средней длиной межгенного интервала (2000 оснований), свойственного Aspergillus, снижение специфичности составило ~5%. Чувствительность программы не зависит от длины анализируемой последовательности и длинны спейсеров. Комплекс программ автоматической аннотации геномов

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

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

2) Предсказание генов по сходству. На этом этапе применяются программы сплайсового выравнивания последовательностей белков и геномной ДНК, например, PROFRAME. Для эффективной работы программ этого класса используется предыдущий этап.

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

4) Применение обученной программы.

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

Sciccharopolyspora erythraea (предварительная сборка, draft) длиной около 2-106 пар оснований. Качество предсказания генов программным комплексом оценивалось по отношению к тестирующему множеству 43 аннотированных последовательностей в GenBank, для которого нуклеотидная чувствительность составила 97%, специфичность - 86%. Относительно низкая специфичнос+ь объясняется неполной аннотацией последовательностей, так как для всех предсказанных генов длинной более 300 нуклеотидов с помощью BlastX было обнаружено существенное сходство с каким-либо белком.

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

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

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

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

Д

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

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

Для фильтрации редких альтернатив, был применен биноминальный тест с параметром 0,01. Значение параметра теста подбиралось исходя из оценки частоты ошибок сплайсосомы (см. далее). Вариант сплайсинга, через который проходит наибольшее количество EST, называется базовым. Пусть через источник и сток альтернативы проходит N EST, а через минорный вариант сплайсинга К EST. Сумма соответствующих биномиальных частот:

является вероятностью того, что событие с количеством EST в минорном пути меньше наблюдаемого (К) встречаются с частотой 0,01. Если вероятность Р>0,95 то отвергается нулевая гипотеза, мы считаем, что частота минорного варианта больше пороговой, иначе альтернатива считается редкой.

Свойства элементарных альтернатив

Наиболее распространенным типом элементарных альтернатив является кассетный экзон 57% - 64% (последнее значение получено без учета редких альтернатив). Около 51% кассетных экзонов являются пропусками экзона в базовой изоформе. Доля удержанных интронов составляет 8% (как с учетом, так и без учета редких делеций и вставок). В противоположность кассетным экзонам, 58% удержанных интронов являются вставками в базовую изоформу. Около 59% всех альтернатив, имеющих EST покрытие, приходится на редкие варианты сплайсинга. Несмотря на значительное число EST человека ~ 4-10й последовательностей, более половины элементарных альтернатив, по-видимому, являются ошибками сплайсинга в генах с высоким уровнем

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

Один из двух вариантов сплайсинга в элементарных альтернативах соответствует более длинной мРНК, чем другой вариант. Были построены распределения частоты включения в мРНК длинного варианта сплайсинга для каждого типа альтернатив. Диапазон частот включения от 0 до 0,4 соответствует вставкам нуклеотидов в базовую изоформу, диапазон от 0,4 до 0,6 соответствует равнозначным альтернативам и диапазон от 0,6 до 1 - делениям в базовой изоформе. Распределения частот включения для альтернативных сайтов (рис. 1 А и Б) и кассетных экзонов (рис. 1 Г) обладают большим сходством что определяется общим механизмом, лежащим в основе альтернативного сплайсинга, - конкуренцией сайтов базовой и альтернативной изоформ. В случае кассетного экзона, по-видимому, решающим является распознавание какого-нибудь одного из альтенативных сайтов.

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

Алгоритм выделения альтернатив позволяет обнаружить координированный альтернативный сплайсинг. Две альтернативы (не обязательно элементарные) называются соседними, если между стоком 5'- альтернативы и источником З'-альтернативы па графе сплайсинга существует единственный путь. Для соседних альтернатив проверялась нулевая гипотеза Н0 о независимости сплайсинга, в этом случае распределение наблюдаемых путей в 3'-альтернативе не зависит от пути следования через 5'- альтернативу. Для проверки гипотезы нужно рассмотреть совместное распределение числа EST па всех возможных парах путей в соседних альтернативах. Это распределение может быть представлено в виде матрицы смежности (рис. 2). Гипотеза о независимости проверяется с помощью точного теста Фишера.

Применение алгоритма идентификации соседних альтернатив позволило обнаружить гены, для которых сплайсинг в 3'- альтернативе зависит от варианта сплайсинга в 5'-альтерпативе (координированный сплайсинг). Из 630 генов, содержащих соседние альтернативы с достаточным для анализа покрытием EST, только для 60 генов был показан координированный сплайсинг. Мы

оценили, что около 25% генов человека содержат более одной альтернативно сплайсируемой области, следовательно, сплайсинг этих областей может быть координирован. Оценка была получена с помощью программы IsoformCounter (см. следующую главу). EST содержат в среднем четыре экзона, кроме того, 3' и 5' концы генов, кодирующих длинные транскрипты мРНК, обычно имеют значительно более высокое покрытие по сравнению с серединой. Мы идентифицировали лишь незначительную часть существующих соседних альтернатив, для которых EST покрытие было достаточно для анализа. Для пяти генов, содержащих области координированного альтернативного сплайсинга, совместное распределение вариантов в соседних альтернативах с высокой значимостью отличалось от ожидаемого исходя из модели Но (pval<0,007). Было показано, что распределение наблюдаемых частот вариантов сплайсинга в этих генах не может быть объяснено деградацией изоформ, содержащих преждевременный СТОП-кодон (NMD). Действие тканевых факторов для четырех из пяти генов так же не может объяснить наблюдаемое распределение.

Наши коллеги на модели минигенов, содержащих два идентичных кассетных экзона EDI гена фибронектина (FN) человека, разделенные фрагментом, из трех последовательных константных экзонов и соответствующими нитронами, показали, что наблюдаемый координированный сплайсинг в соседних альтернативах зависит от скорости транскрипции. Высокая скорость транскрипции уничтожает координационный эффект. Координация сплайсинга в соседних альтернативах зависит от промотора, под управлением которого инициируется транскрипция, т.е. координация исчезает под управлением промотора а-глобинового гена и восстанавливается при ингибировании полимеразной активности pol-11. Альтернативный сплайсинг одного из пяти генов (РСВР2), для которых координированный сплайсинг был показан в EST с высоким уровнем значимости, был исследован экспериментально. Было показано, что существует высокая корреляция между распределением частот изоформ, оцененным по EST-данным и по результатам ОТ-ПЦР.

Рисунок 2. Координированный альтернативный сплайсинг гена РСВР2. В каждой альтернативе содержится по два пути a, b в 5' и с, d в 3'. В таблицах показано число EST в EDAS, подтверждающее каждую пару путей в соседних альтернативах.

с d

а 0 12

b 44 48 F

Распределение вариантов сплайсинга гена РСВР2 приведено на рисунке 2. Если 5' кассетный экзон пропускается (вариант а) то, 3' кассетный экзоп включается в состав 100% мРНК. Если 5' кассетный экзон включается (вариант Ь), то доля экзона d п мРНК снижается до 52%. Таким образом, введение задержки транскрипции РНК должно приводить к увеличению доли включения 5' экзона Ь и как следствие снижать долю мРНК, содержащих экзон d. Этот эффект был обнаружен после обработки трех клеточных линий ингибитором полимеразиой активности ОЯВ, что привело к снижению доли включения экзона d примерно на 50%.

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

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

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

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

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

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

(1) Идентификация изоформ.

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

(2) Путциация трансляции.

В генах эукариот трансляция инициируется по механизму линейного сканирования. 40S субъединица рибосомы связывается с кэп-структурой на 5'-конце мРНК, затем начинает скольжение в 3' направлении до ближайшего ATG кодона, где инициирует трансляцию, если этот

кодон находится в подходящем контексте. Стабильные шпильки или ATG кодоны 5' левее сайта инициации трансляции снижают эффективность линейного сканирования. Небольшая доля мРНК (2-8%) содержит внутренний сайт инициации трансляции. IsoformCounter допускает существование не более 2-х ATG кодонов 5' левее ORF-СТЛРТ-вершины. Данный фильтр не применяется, если СТАРТ-вершина соответствует началу выравнивания ДНК - белок.

(3) Короткие изоформы.

Алгоритм допускает изоформы, если длина ORF превышает 50% от средней длины белка в EDAS. Под средней длинной белка понимается средняя длина RefSeq изоформ данного гена. Если для гена не существует RefSeq изоформ, то его средняя длина равна средней длине всех известных белков данного гена. Изоформы длины которых меньше 33 аминокислот не рассматриваются.

(4) Согласованность с белками

Необходимо чтобы изоформа имела как минимум одну общую аминокислоту с известной последовательностью белка, кодируемой данным геном. Основное назначение данного фильтра -удалять длинные ORF в 3' нетранслируемой мРНК.

(5) Терминация трансляции.

Трапскрипты, содержащие преждевременный СТОП-кодон возникающие, например, в результате ошибочного сплайсинга, разрушаются с помощью механизма NMD. Так как не существует правила, по которому можно различать функциональные мРНК изоформы и ошибки сплайсинга, мы применили специальный фильтр, имитирующий механизм NMD. IsoformCounter допускает сплайсированные изоформы в которых СТОП-кодон находится на расстоянии не более 55 нуклеотидов от последней границы между экзонами. Фильтр не применяется к ORF-СТОП-вершинам, соответствующим З'-концам выравниваний с последовательностями белков. Процедура нормализации достоверности экзонов и нитронов

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

небольшим числом EST клонотек. Для того чтобы сравнивать достоверность альтернативного сплайсинга различных генов необходимо произвести нормировку на общее число EST. IsoformCounter использует вероятностную модель, которая позволяет оценить минимальную достоверность экзонов и нитронов, необходимую для исключения ошибочного сплайсинга на данном уровне экспрессии гена. Пусть для некоторого гена известен уровень экспрессии \ (среднее число транскриптов на ген). Мы предполагаем, что ^=f(N), где N - наблюдаемое число EST. Обозначим P(N) вероятность, что клетка содержит хотя бы одну неверно сплайсированную EST, Р=1-(1-аД где а - вероятность ошибки при сплайсинге одного интрона. Вероятность, что ошибочный интервал будет иметь подтверждение к библиотеками клонов (клонотеками) равна Рк<Р, где Р - уровень значимости - вероятность, с которой мы допускаем ошибку сплайсинга. Решая неравенство относительно к, получаем выражение: к(фМпф)/1п(1-(1-а)(). Уровень экспрессии был оценен как £=N/5. Оценка а составила 0,01 (обоснование будет дано далее). В результате было использовано следующее выражения для определения порогового числа EST-клонотек для того чтобы принять или отклонить экзон или интрон: k(N)=[-l/ln(l-0,99N3)].

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

Программа IsoformCounter была применена для анализа альтернативного сплайсинга (АС) генов базы данных EDAS. Среди всех генов в EDAS, которые содержат интроны, 77% генов имеют более одной функциональной изоформы. Около 20% генов в геноме человека являются одноэкзонными, следовательно, оценка доли альтернативно сплайсируемых генов составляет примерно 60% (75% из 80%). Большинство генов (91%) имеют относительно небольшое число изоформ (от 1 до 15). Для каждого гена было подсчитано число альтернативных и константных сегментов в изоформе, кодирующей самый длинный белок. Было получено следующее

распределение: 52-72% генов содержат единственный альтернативный сегмент, нижняя оценка получена на основании EST данных, верхняя - только белков, в промежутке (62%) находится оценка, полученная с использованием мРНК и белков вместе; 20-29% содержат два сегмента; 612% три; 1-6% более трех. Таким образом, около 25% генов содержат более одного альтернативного сегмента. Сплайсинг в альтернативных сегментах этих генов может быть координирован.

Были рассмотрены следующие функциональные категории онтологии генов GO [http://wwvv.geneontology.org/GO.doc.html]: "передача сигналов посредством малых ГТФаз" (145 генов), "катаболизм" (512 генов), "репликация ДНК и хромосомный цикл" (99 генов), "рибосома" (123 гена). Число функциональных изоформ было оценено с применением нормализации на общее число EST. Было получено значимое отличие распределений числа изоформ генов из категорий "рибосома" и "передача сигналов посредством малых ГТФаз" от распределения по всем генам (р = 0,003 U-тест Манн-Уитни). В обеих категориях наблюдалось меньшее, чем в среднем, число изоформ. Категория "рибосома" содержит 46% константных генов, несмотря на очень большое EST покрытие, что является результатом нормализации. Гены, принадлежащие категории "репликация ДНК и хромосомный цикл", имеют больше изоформ, чем в среднем (р = 0,07 U-тест Манн-Уитни), среди них более высокий процент генов являются альтернативно сплайсируемыми (имеют две и более изоформы). Для изучения связи между АС и белок-белковыми взаимодействиями из базы данных МРР1 было выбрано 198 пар взаимодействующих белков, обе составляющих которых кодируются генами, представленными в EDAS (всего 262 гена). Доля альтернативного сплайсинга была выше среди генов, участвующих хотя бы в одном белок-белковом взаимодействии. Дефицит константных генов составил 17-30%, избыток альтернативно сплайсируемых - 10-25%. Несмотря на относительно небольшие различия наблюдаемых и ожидаемых значений, эти различия статистически значимы (р < 0,1-1%) па всех рассмотренных уровнях достоверности: только белки, белки и мРНК, EST с различным числом клонотек, нормализованные EST.

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

Программа IsoformCounter для каждого нитрона позволяет определить существование хотя бы одной изоформы, не попадающей под действие механизма NMD и кодирующей достаточно длинную аминокислотную последовательность (фильтры (1)-(5)). Мы полагаем, что если для интрона не существует ни одной изоформы, то с высокой долей вероятности, она не может кодировать функциональный белок. Рассмотрим множество (errorJntrons), состоящее из нитронов, подтвержденных только последовательностями EST через которые IsoformCounter не может провести ни одной изоформы. Кроме того, каждый из интронов, входящих в состав множества error_introns пересекается с каким-либо нитроном, достоверным на уровне белка. Рассмотрим множество последовательностей EST (Good) из EDAS, которые сплайсированы только в тех нитронах, которые имеют достоверность на уровне белка, обозначим Ngooj - число таких EST. Обозначим Error множество EST, сплайсированных хотя бы в одном из ошибочных интронов из множества error Jntrons , при этом эти EST должны быть так же сплайсированы хотя бы в одном белковом нитроне, пусть Nermr - число таких EST. EST из множества Error соответствуют транскриптам, в которых произошла ошибка при сплайсинге интронов в достоверно кодирующих изоформах. Обозначим а - вероятность ошибки сплайсинга одного интрона. Пусть EST в среднем содержит <ni> интронов, тогда (1- а) - вероятность, что EST не содержит ни одной ошибки сплайсинга. Вероятность ошибки сплайсосомы может быть получена, как решение уравнения: (1 -(1 -cc)<nl>)(Neoml + Nrrmr) = Nerrlir, результаты приведены в таблице 1.

Первая оценка получена без ограничения на число клонотек, подтверждающих EST интропы, входящие в состав множества error Jntrons. Вторая оценка получена с ограничением одна EST-клонотека.

Таблица 1 Оценка вероятности ошибки сплайсинга одного интрона - а.

любое число клонотек 1 EST клонотека

Nqood (ЧИСЛО EST) 1848958 1848958

Nerror (ЧИСЛО EST) 38845 17926

а 0,0069 0,0032

Для того, чтобы понять насколько приведенаая оценка близка к истине, будем рассуждать следующим образом. Рассмотрим ген содержащий ni нитронов. Если а - вероятность ошибки при сплайсинге одного интрона, то (I- а)"' > Р, где /? - доля нормальных мРНК необходимых для эффективной экспрессии гена. Несмотря на то, что среднее число нитронов в генах человека относительно невелико (5,5), некоторые гены имеют очень большое число нитронов. Экстремальный пример - ген титин содержит 363 экзона. Величина ошибки сплайсосомы должна находится в интервале 0,001 - 0,003, чтобы обеспечить достаточную эффективность экспрессии длинной изоформы титина (/i=fl,5 - 0,8). В главах 3 и 4 использована верхняя оценка ошибки сплайсомы -0,01, так как ошибка сплайсинга может зависеть от контекста сайтов базовой изоформы. Представляет интерес сравнение вероятности ошибки сплайсинга с ошибкой РНК-полимеразы и ошибкой трансляции. Частота ошибок РНК полимеразы, не имеющей корректирующей активности, составляет 1 на 104 нуклеотидов. Ошибки трансляции происходят с частотой 1 на 104 аминокислотных остатков. Длина кодирующей части мРНК среднего гена -900 нуклеотидов, тогда вероятность того, что при транскрипции в ней содержится хотя бы одна ошибка - 0,09! Таким образом, около 10% мРНК содержат ошибки в кодирующей части гена, внесенные РНК полимеразой, примерно такой же уровень ошибок вносится при сплайсинге пре-мРНК генов с 10 (верхняя оценка а) - 100 (нижняя оценка) нитронами

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

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

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

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

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

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

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

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

Список работ опубликованных по теме диссертации

1. A. D. Neverov, М. S. Gelfand, A. A. Mironov 2003. GipsyGene: A Statistics-Based Gene Recognizer for Fungal Genomes. Biophysics (Moscow), Vol. 48, Suppl. I, 2003, pp. S71-S75.

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

3. J. Fededa, E. Petrillo, M. Gelfand, A. Neverov, S. Kadener, G. Nogues, F. Pelisch, F. Baralle,

A. Muro, A. Kornblihtt. 2005. A polar mechanism coordinates different regions of alternative splicing within a single gene. Mol Cell. 2005, Vol. 19, N. 3, pp.393-404.

4. P.H. Нуртдинов, А.Д. Неверов, Д.Б. Малько, И.А. Космодемьянский, Е.О. Ермакова,

B.Е. Раменский, А.А. Миронов и М.С. Гельфанд. 2006. EDAS, база данных альтернативно епланенруемых генов человека. Биофизика, 2006, т. 51, номер 4, стр. 589-592.

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

6. A.D. Neverov, M.S. Gelfand, A.A. Mironov. 2002. Gene prediction in genomics DNA of Aspergillus. Proceedings of the Third International Conference on Bioinformatics of Genome Regulation and Structure. (BGRS'2002), Novosibirsk, Russia. Vol. 1, p. 116.

7. A.D. Neverov. 2003. GipsyGene: A HMM-based gene recognitional tool for lower fungi. . Proceedings of the First International Moscow Conference on Computational Molecular Biology (MCCMB'03), Moscow, Russia, p. 161.

8. A.D. Neverov, L. Milanesi. 2005. A pipeline for computational gene recognition in the Sacharopolyspora erythraea genome. Proceedings of the Second International Moscow Conference on Computational Molecular Biology (MCCMB'05), Moscow, Russia, p. 246.

9. R.N. Nurtdinov, A.D. Neverov, D.B. Malko, I.A. Kosmodemyansky, A.A. Mironov, M.S. Gelfand. 2005. EDAS: EST-derived alternative splicing database, Proceedings of the Second International Moscow Conference on Computational Molecular Biology (MCCMB'05), Moscow, Russia, p. 259.

10. V. Ramensky, R. Nurtdinov, A. Neverov, A. Mironov, M. Gelfand. 2006. Proceedings of the 5th International. Conf. Bioinformatics of Genome Regulation and Structure (BGRS'2006), Novosibirsk, Russia, p. 211.

Заказ № 6/02/07 Подписано в печать 06.02.2007 Тираж 80 экз. Усл. пл. 1,5

ООО "Цифровичок", тел. (495) 797-75-76; (495) 778-22-20 www.cfr.ru; е-таИ:info@cfr.ni

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

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

Апробация

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 случаев ошибочного сплайсинга на интрон.

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

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.