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

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

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

03.00.28 Биоинформатика

АВТОРЕФЕРАТ

диссертации на соискание ученой степени доктора физико-математических наук

Ройтберг Михаил Абрамович

Москва - 2009

Работа выполнена в Учреждении Российской Академии Наук

Институт математических проблем биологии РАН

Официальные оппоненты: доктор биологических наук,

профессор В.В. Поройков

доктор физико-математических наук

О.В. Галзитская

доктор физико-математических наук

А.М. Райгородский

Ведущая организация: Учреждение Российской академии наук Институт теоретической и экспериментальной биофизики РАН

Защита состоится «______» ______________2009 г. в ________ч. на заседании диссертационного совета Д 002.077.02 при Учреждении Российской академии наук Институт проблем передачи информации РАН им. А.А. Харкевича по адресу 127994, Москва, ГСП-4, Большой Каретный переулок, д.19, стр.1.

С диссертацией можно ознакомиться в библиотеке Учреждения Российской академии наук Институт проблем передачи информации РАН им. А.А. Харкевича

Автореферат разослан «______»________________2009г.

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

доктор биологических наук, профессор Рожкова Г.И.

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

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

Практически одновременно с накоплением данных о биологических последовательностях (в 60-х - 70-х годах XX века) происходило развитие прикладной теории алгоритмов - разработка базовых алгоритмов анализа символьных последовательностей и связанных с ними алгоритмов сортировки и алгоритмов на графах. Начиная с 70-х годов ХХ века аппарат прикладной теории алгоритмов начал применяться для анализа (прежде всего - сравнения) биологических последовательностей. При этом достаточно быстро стало понятно, что постановки задач должны максимально учитывать специфику предметной области, а собственно алгоритмы поиска (построения) искомых объектов должны дополняться исследованием статистической значимости полученного результата и/или его соответствия эталонным (экспериментально подтвержденным) результатам - для тех случаев, когда эти результаты известны.

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

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

Основные понятия выравнивания биологических последовательностей были сформированы в конце 70-х - начале 80-х годов XX века в работах Нидлмана, Вунша, Смита, Уотермана, Селлерса, Санкоффа, Эриксона, Туманяна и др. В частности, были предложены основанные на методе динамического программирования алгоритмы, которые строят оптимальное выравнивание последовательностей для различных классов весовых функций. Алгоритм для линейных штрафов имеет время работы O(m*n); алгоритмы построения оптимального выравнивания для выпуклых весов делеций и для произвольных весов делеций имеют временную сложность соответственно O(m*n*(m+n)) и O(m2*n2) (m и n - длины последовательностей). Однако оставался открытым вопрос - существует ли класс весовых функций, более широкий, чем линейные функции, допускающий построение оптимального выравнивания за время O(m*n)? Кроме того, указанные алгоритмы имеют два существенных недостатка с точки зрения биологических приложений. Во-первых, оптимальное выравнивание аминокислотных последовательностей белков в ряде важных случаев по внутренним причинам не может воспроизвести выравнивание этих белков, согласованное с их пространственной структурой. Во-вторых, не разработана теория выбора значений параметров штрафов за делеции. Как правило, выбор параметров осуществлялся эмпирически. Это определяет важность поиска путей инкорпорирования содержательной биологической информации, как при собственно выравнивании, так и при подборе параметров штрафов.

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

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

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

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

Цель исследования

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

Задачи исследования

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

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

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

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

Методы исследования

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

Научная новизна

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

Теоретическая и практическая ценность

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

Основными из этих алгоритмов являются:

1) алгоритм построения оптимального выравнивания нуклеотидных последовательностей при штрафах за делеции, задаваемых кусочно-линейными функциями;

2) алгоритм построения множества Парето-оптимальных выравниваний относительно векторной весовой функции;

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

4) алгоритм предсказания внутренних циклов во вторичной структуре РНК;

5) алгоритм выравнивания вторичных структур РНК при заданной вторичной структуре;

6) алгоритм выравнивания геномов;

7) алгоритм вычисления чувствительности затравок для поиска локальных сходств;

8) алгоритм оценки значимости кластеров регуляторных сайтов

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

Апробация результатов

Материалы диссертации докладывались на международных и всероссийских конференциях и семинарах, в том числе: Московском семинаре по компьютерной генетике; отчетных конференциях программы «Геном человека»; II Съезде биофизиков России (Москва, 1999), симпозиуме «The Informatics of Protein Classification» (Университет Ратгерс, США, 2000), III, IV,V,VI международных конференциях по биоинформатике и регуляции структуры генома (Новосибирск, 2002, 2004, 2006, 2008), I, II и III Московских международных конференциях по вычислительной биологии (2003, 2005, 2007), международной конференции Combinatorial Pattern Matching-2004 (Стамбул, Турция), международной конференции Workshop on Algorithms in Bioinformatica (Пальма де Майорка, Испания, 2005), международной конференции «Implementation and Application of Automata» (Прага, Чехия, 2007), международном симпозиуме NETTAB (Варенна, Италия, 2008).

С использованием материалов диссертации автором сделаны доклады в NCBI USA (2000), Georgia Tech (2005), INRIA (2003, 2007), на семинарах в Институте белка РАН, Московском физико-техническом институте, факультете биоинженерии и биоинформатики МГУ и ряде других учреждений.

Публикации

Основные материалы диссертации изложены в 37 статьях в реферируемых научных изданиях (из них 33 в соавторстве).

Структура и объем работы

Диссертация состоит из введения, пяти глав, заключения и списка литературы (345 наименований). Полный объем диссертации составляет 223 страницы, количество рисунков - 29, количество таблиц - 17.

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

Во введении дана общая характеристика работы и приведены основные определения и обозначения. Глава 1 посвящена обзору литературы, отмечена связь проанализированных работ с предметом исследования диссертации. Результаты работы представлены в главах 2 - 5.

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

Напомним основные понятия. Выравнивание символьных последовательностей u и v - это тройка <u, v, S>, где S ={<i1, j1 >, …<in, jn >} - набор пар позиций в словах u и v соответственно, таких, что 1 ? i1< …< in ? |u|; 1 ? j1< …< jn ? |v|. Имеется в виду, что ik-я позиция слова u сопоставлена с jk-й позицией слова v (k = 1, …, n), а непустые фрагменты вида u[ik+1, ik+1-1] и v[jk+1, jk+1-1] удалены (здесь k = 0, …, n; i0 = j0 = 0; in+1=|u|+1; jn+1=|v|+1). Пары позиций <is, js > называются склейками. Вес P(A) выравнивания A определяется как функция (обычно - разность) весов сопоставлений букв и штрафов за удаление фрагментов. Предполагается, что нулевой (начальной) буквой любой символьной последовательности является начальный маркер, один и тот же для всех последовательностей.

В разделе 2.1 введено понятие кусочно-линейных функций штрафов (весовых функций) за удаление фрагмента и представлен соответствующий алгоритм построения оптимального выравнивания. Предложенный класс функций является наиболее широким из известных классов весовых функций, для которых построен квадратичный алгоритм построения оптимального выравнивания (под квадратичным алгоритмом понимается алгоритм со временем работы, пропорциональным длинам сравниваемых последовательностей). Класс весовых функций, рассмотренный в разделе 2.2, наоборот, - простейший класс, при котором задача построения оптимального выравнивания двух последовательностей не сводится к хорошо изученной задаче построения максимальной общей подпоследовательности (Longest Common Subsequence, LCS). Для указанного класса функций предложен адаптивный алгоритм построения оптимального выравнивания. Хотя оценка времени алгоритма в худшем случае является квадратичной, время работы алгоритма для близких последовательностей почти линейно. В частности, если для последовательностей v1 и v2 длины m существует оптимальное выравнивание, содержащее t несовпадений и s удаленных символов, то время работы предлагаемого алгоритма составляет O((t+s)m).

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

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

1) штрафы за удаления фрагментов на концах последовательностей могут быть произвольными;

2) штраф за удаление фрагмента может зависеть от граничных позиций удаляемого фрагмента;

3) зависимость штрафа от длины фрагмента может задаваться произвольной кусочно-линейной функцией;

4) штрафы за удаления фрагментов в каждой из сравниваемых последовательностей могут задаваться по-своему;

5) вес сопоставления символов u[i] и v[j] задается произвольной функцией з(i, j, u, v).

Вес выравнивания последовательностей определяется как разность S - D, где S - сумма весов сопоставлений букв, D - сумма штрафов за удаление фрагментов.

Идея алгоритма построения оптимального выравнивания символьных последовательностей (слов) v1 и v2 при заданной кусочно-линейной системе весовых функций состоит в следующем. Последовательно (в цикле по i, а внутри него - в цикле по j) вычисляются т.н. оптимальные (i, j)-выравнивания A[i, j] и их веса P[i, j]; здесь (i, j)-выравнивание слов u и v - это такое выравниваний префиксов u[1, i], v[1, j], в котором сопоставлены позиции u[i] и v[j]. Одновременно вычисляется выравнивание M[i, j] - оптимальное выравнивание слов u и v, в котором последнее сопоставление позиций (x, y) удовлетворяет условиям x?i; y ? j. Выравнивание M[|u|, |v|] будет искомым. Для вычисления веса P[i, j] и соответствующего ему выравнивания A[i, j] алгоритм поддерживает k1 * k2 рекурсивно вычисляемых величин m[s, t], где ki. - количество интервалов линейности у функции штрафов за удаление фрагментов внутри слова vi. При вычислении веса P[i, j] величина m[s, t] хранит (в приведенной форме, см. подробнее п.2.1.5 диссертации) максимальный вес среди весов таких (x, y)-выравниваний, что длина i-x-1 попадает в s-й интервал линейности слова u, а длина j-y-1 - в t-й интервал линейности слова v. При этом каждая из величин m[s, t] для очередной пары (i, j) перевычисляется за конечное время, что и обеспечивает указанное время работы алгоритма O(k1*k2*|v1|*|v2|). Память, нужная алгоритму, оценивается формулой (см. раздел 2.1.10 диссертации)

,

где c1, …,c4 - константы. Однако большая часть этой памяти нужна лишь при восстановлении оптимального выравнивания как пути в специальном дереве (дереве опорных склеек). Поэтому для потребной оперативной памяти получаем оценку O(*m2).

В разделе 2.2 вес выравнивания G последовательностей v1 и v2 длины соответственно m1 и m2 задается формулой

W(G) = p - f*r - d*( m1-p-r) - d*( m2-p-r)

где p - количество совпадений в выравнивании; r - количество несовпадений; m1-p-r - количество символов, удаленных в первой последовательности; m2-p-r - количество символов, удаленных во второй последовательности; f, d - весовые коэффициенты. Таким образом, система весовых функций описывается двумя коэффициентами f и d.

Базовыми понятиями алгоритма, описанного в разделе 2.2, являются понятие опорного множества и специального веса начального выравнивания. Начальным выравниванием слов v1 и v2 называется выравнивание префиксов v1[0, x1] и v2[0, x2]. Начальное выравнивание G = < v1[0, x1], v2[0, x2], Q> слов v1 и v2 называется граничным, если x1 = |v1| или x2 = |v2|.

Определение 2.2.2. Множество S начальных выравниваний слов v1 и v2 называется опорным множеством, если в S найдется выравнивание, которое является префиксом некоторого оптимального выравнивания слов v1 и v2.

Например, множество S0, состоящее из единственного начального выравнивания A0 = < v1[0], v1[0], {(0,0)}>, является опорным.

Определение 2.2.1. Пусть G = < v1 [0, x1], v2 [0, x2], Q> - начальное выравнивание слов v1 и v2 длин m1 и m2. Специальным весом выравнивания G называется величина

SP(G) = P(G)+ c* min( m1-x1, m2-x2 ) - d*|( m1-x1)-( m2-x2 ) |.

Начальное выравнивание G называется максимальным в множестве S, если в S нет выравнивания, которое имеет специальный вес больше, чем SP(G). Опорное множество S называется терминальным, если некоторое граничное выравнивание является максимальным в нем.

Специальный вес выравнивания обладает следующими свойствами.

1. Пусть G - выравнивание слов v1 и v2 и начальное выравнивание B является префиксом выравнивания G. Тогда

P(G)? SP(B).

2. Пусть S - опорное множество для слов v1 и v2; B S - граничное начальное выравнивание, которое является максимальным в S. Тогда выравнивание B' слов v1 и v2, имеющее тот же набор сопоставлений, что и выравнивание B - оптимальное выравнивание слов v1 и v2.

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

Построение терминального опорного множества и соответствующего граничного максимального выравнивания выполняется с помощью построения т.н. канонической последовательности опорных множеств {Si}. Множество S0 состоит из одного выравнивания <v1[0,0], v2[0,0], {(0,0)} >. Пусть опорное множество Si построено. Если Si - терминальное множество, то процесс построения канонической последовательности закончен. В противном случае выбирается некоторое начальное выравнивание Li Si (лидер множества Si) и это выравнивание заменяется таким набором продолжений D(Li), что любое оптимальное выравнивание, продолжающее выравнивание Li, является продолжением одного из выравниваний множества D(Li). Это гарантирует то, что полученное множество Si+1 будет опорным. Далее из построенного множества, возможно, удаляются некоторые избыточные выравнивания, т.е. такие выравнивания, удаление которых оставляет множество Si+1 опорным. Таким образом,

Si+1 (Si - { Li }) D(Li).

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

Время работы алгоритма Time и потребная память Space описываются соотношениями Space = O(K+ |v1| + |v2|), где K - длина канонической последовательности; Time = O(Klog(|v1| + |v2|). При этом для величины K имеет место следующая оценка. Пусть для слов v1 и v2 длин соответственно m1 и m2 существует оптимальное выравнивание, содержащее t несовпадений и s удаленных символов. Тогда

K ? ( 2t +s + 2.5*|m1 - m2| ) *min(m1 , m2).

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

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

Определение 2.3.1. Пусть k ?2 - целое число. Векторная весовая функция - это функция, сопоставляющая каждому выравниванию A k-мерный вектор V(A), называемый (векторным) весом выравнивания A.

Примером векторного веса является вес

V(A) = (NumMatch(A), -NumDel(A)),

где k = 2; NumMatch(A) и NumDel(A) - это соответственно число совпадений и число удаленных символов в выравнивании A. Компоненты векторного веса (в данном случае NumMatch(A) и -NumDel(A)) называются элементарными весовыми функциями. Другими примерами элементарных весовых функций являются:

-NumGap(A) - количество удаленных сегментов последовательностей («делеций», “gaps”);

WeightMatch(A) - суммарный вес сопоставлений символов относительно выбранной матрицы замен;

-MisMatch(A) - количество несовпадений (соответствует использованию единичной матрицы замен);

Определение 2.3.4. Пусть S1 и S2 - последовательности; V - векторная весовая функция. Выравнивание A последовательностей S1 и S2 называется Парето-оптимальным относительно весовой функции V, если V(A) является Парето-оптимальным вектором в множестве векторных весов всех возможных выравниваний последовательностей S1 и S2.

Важность множества Парето-оптимальных выравниваний определяется следующим наблюдением. Пусть V(A) = < V1(A), …, Vk(A) > - векторная весовая функция; g(x1, …,xk) - функция k переменных, монотонно неубывающая по каждому из аргументов, ц(A) = g(V1(A), …, Vk(A)) - скалярная весовая функция и B - оптимальное выравнивание некоторых последовательностей относительно функции ц. Тогда B - оптимальное выравнивание этих последовательностей относительно вектор-функции V. В частности, если B - оптимальное выравнивание некоторых последовательностей относительно линейной комбинации функций V1(A), …, Vk(A) с положительными коэффициентами, то B - оптимальное выравнивание этих последовательностей относительно вектор-функции V(A).

В разделе 2.3 представлен алгоритм построения множества всех Парето-оптимальных выравниваний данных последовательностей S1 и S2 относительно данной весовой функции V(A). Этот алгоритм является алгоритмом динамического программирования и основан на соотношении дистрибутивности между операцией сложения векторов и операцией взятия Парето-подмножества объединения двух Парето-оптимальных множеств. Более точно, пусть c - число, множества T, T1, T2 Rk, причем T, T1, T2 - Парето-оптимальные множества; через Pareto(T) обозначается Парето-подмножество множества T;. Рассмотрим следующие операции:

c T = Pareto({< x1+c, …, xk +c>| < x1, …, xk > T}),

T1 T2 = Pareto(T1 T2).

Тогда

c ( T1 T2) = (c T1 ) (c T2).

Доказательство непосредственно следует из определения операций и и определения Парето-множества.

Для последовательностей длины m1 и m2 время работы такого алгоритма оценивается как O(с(m1, m2)*m1*m2), где коэффициент с(m1, m2) зависит от выбранной весовой функции и определяется временем выполнения операций и для этой функции. В некоторых случаях приведенная оценка может быть улучшена. В частности, оценку времени для случая упомянутой выше весовой функции VD(A) = (NumMatch(A), -NumDel(A)) дается следующей леммой:

Утверждение 2.3.4. Пусть S1 и S2 - последовательности; их длины равны соответственно n и m. Пусть p - длина наибольшей общей подпоследовательности S1 и S2; d = m + n - 2p и r = min(p, d).

Алгоритм Pareto_Align_Del_Dmax строит множество Парето-оптимальных весов и Парето-оптимальных выравниваний для последовательностей S1 и S2 относительно весовой функции VD(A) = (Match(A), -NumDel(A)) за время O(min(n, 2d)*m*n*log(r)).

Такая же оценка верна и для весовой функции VG(A) = (Match(A), -NumGap(A)).

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

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

Предлагаемый подход основан на следующем наблюдении. Рассмотрим две последовательности u и v и фиксируем некоторую весовую матрицу замен. Для произвольного целого g, через S(g) обозначим наибольший суммарный вес сопоставлений среди выравниваний, содержащих не более g удаленных фрагментов. Через D(g) обозначим «производную» S(g):

D(g) = S(g+1) - S(g).

Значения S(g) могут быть найдены путем построения Парето-оптимальных выравниваний относительно векторной весовой функции

WG(A) = (WeightMatch(A), -NumGap(A))

Допустим, что удалением некоторого количества фрагментов из последовательностей u и v можно получить достаточно похожие последовательности равной длины; эти фрагменты назовем «правильными». Удаление «правильного» фрагмента восстанавливает соответствие между (обычно, протяженными) гомологичными фрагментами исходных последовательностей и, тем самым, ведет к существенному увеличению веса S(g) (значение D(g) велико). Когда все «правильные» фрагменты удалены, новые удаления уже не мотивированы биологически и, следовательно, не могут привести к существенному увеличению S(g) (значение D(g) мало). Наибольшее значение параметра g, при котором значение D(g) велико, назовем критическим значением, а соответствующее ему выравнивание - критическим выравниванием. Мы полагаем, что критическое выравнивание наиболее (среди других Парето-оптимальных выравниваний) соответствует биологически корректному выравниванию данных последовательностей. Это приводит к следующему методу выравнивания последовательностей.

1. Построить множество Парето-оптимальных выравниваний относительно векторной весовой функции WG(A) = (WeightMatch(A), -NumGap(A)).

2. Определить критическое значение количества удаленных фрагментов g и взять критическое выравнивание в качестве результата.

Успешность применения этого метода зависит от алгоритма, который определяет порог, отделяющий «малые» и «большие» значения D(g). Компьютерные эксперименты показали, что для нуклеотидных последовательностей со степенью сходства не менее 30% и для аминокислотных последовательностей со степенью сходства не менее 20% выбор такого порога возможен. Получена верхняя оценка для величины Dcr «малых» значений D(g):

Dcr ? log(L)/log(1/(1 - p)) -1, (2.4.1)

где L - средняя длина последовательностей, p - доля совпадающих букв. Эта оценка хорошо согласуется с данными компьютерных экспериментов при p ? 0.4 (см. таблицу 2.4.8). Величина ошибки увеличивается с убыванием p и возрастанием L.

В свою очередь, ожидаемое значение M величины D(g) при удалении «правильного» фрагмента (что соответствует «докритической» области значений D(g)) зависит не только от L и p, но и от других характеристик сравниваемых последовательностей. Такими характеристиками являются: длина d удаляемого фрагмента, длина D каждого из гомологичных фрагментов, которые оказались сопоставленными правильно в результате удаления; доля совпадений p0 при сопоставлении негомологичных (случайных) фрагментов последовательностей. Для среднего значения M имеем формулу:

M = D(p - p0) - pd.

Таблица 2.4.8.

Оценки значений величины Dcr

L	P
	0.3	0.4	0.5	0.6	0.7	0.8	0.9
200	FORM	13.8	9.4	6.6	4.8	3.4	2.3	1.3
	EXP	10	8	6	4	3	2	1
300	FORM	15	10.2	7.2	5.2	3.7	2.5	1.5
	EXP	10	8	6	5	3	2	1
700	FORM	17.3	11.8	8.5	6.1	4.4	3.07	1.8
	EXP	11	9	7	5	3	3	1
1000	FORM	18.3	12.5	9	6.6	4.7	3.3	2
	EXP	11	9	7	6	4	3	2

Верхняя граница для величины Dcr, вычисленная по формуле (2.4.1), строки FORM, и максимальные значения величины Dcr, полученные в ходе компьютерных экспериментов (см. п.2.4.2 диссертации, строки EXP. Данные приведены для различных длин L и уровней сходства p.

Если указанное значение будет менее Dcr, то удаление таких «правильных» фрагментов не может быть диагностировано на фоне статистического шума. Это показывает статистические пределы применимости алгоритмических методов восстановления биологически корректных выравниваний. Другие ограничения возможностей алгоритмических методов рассмотрены в главе 3.

Таким образом, в главе 2 рассмотрена общая задача построения оптимального выравнивания символьных последовательностей, в частности, - проблема выбора штрафов за удаление фрагментов. В главе 3 рассматривается более специальная задача - выравнивание биологических последовательностей. Центральная тема этой главы - учет пространственной структуры сравниваемых молекул при сравнении их первичных структур (последовательностей). Имея это ввиду, мы ограничиваемся рассмотрением наиболее распространенного в современной биоинформатике класса весовых функций удаления фрагментов - аффинными функциями. Это позволяет существенно упростить изложение, хотя приведенные в разделах 3.2 и 3.3 алгоритмы могут быть обобщены на случай рассмотренных в главе 1 кусочно-линейных функций. В разделе 3.4. представлен оригинальный алгоритм предсказания вторичной структуры РНК.

Раздел 3.1 посвящен изучению связи между биологически корректными выравниваниями аминокислотных последовательностей и выравниваниями, полученными с помощью алгоритма Смита-Уотермана (SW) - наиболее распространенного в настоящее время алгоритма построения оптимального выравнивания последовательностей. В качестве «биологически корректных» (эталонных) выравниваний использованы выравнивания, основанные на наложении пространственных структур белков. Адекватность такого подхода обоснована существенно большей консервативностью пространственной структуры белков по сравнению с их первичной структурой. В качестве источника эталонных выравниваний использовалась база данных BaliBase [http://www-igbmc.u-strasbg.fr/BioInfo/BAliBASE/]. Большая часть сравниваемых пар последовательностей (всего 583 пары) имела процент совпадений %ID от 15% до 40%.

Исследована зависимость степени сходства между алгоритмическими и структурными выравниваниями от степени сходства между сравниваемыми последовательностями и выявлены причины расхождений между этими выравниваниями. В качестве количественной оценки качества алгоритмически полученных выравниваний использовались две взаимодополняющих меры: (1) Точность выравнивания (обозначение: Ali_Acc) равна отношению количества одинаковых сопоставлений (I) в обоих выравниваниях к общему количеству сопоставлений в эталонном выравнивании (G): Ali_Acc = I/G*100 %. (2) Достоверность выравнивания (обозначение: Ali_Conf), равна отношению количества одинаковых сопоставлений в обоих выравниваниях (I) к общему количеству сопоставлений в алгоритмически построенном выравнивании (A): Ali_Conf = I/A*100 %. Неформально говоря, точность Ali_Acc показывает, какую долю эталонного выравнивания удалось восстановить, а достоверность Ali_Conf - насколько можно доверять построенному выравниванию.

Как видно из рис. 3.1.2, алгоритм SW может строить выравнивания, хорошо совпадающие с эталонными, только при уровне сходства сравниваемых белков более 30-40%. Этот диапазон уровня гомологии (ID > 30%) примерно совпадает с известным порогом %ID, выше. которого можно достоверно восстановить структурное выравнивание, зная только последовательности [%56, %58-60]

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

Рис. 3.1.2. Зависимость точности восстановления эталонных выравниваний методом SW от уровня сходства (%ID) сравниваемых белков. Каждая точка соответствует паре эталонных белков. Х-координата точки равна %ID пары, а Y-равный координата - значению точности. Зависимость для достоверности практически такая же

При уровне гомологии меньше 10% метод SW не может восстановить правильное выравнивание даже частично. Для диапазона уровня гомологии от 10% до 30% выравнивания Смита-Уотермана показывают очень широкий разброс точности и достоверности. Для разных пар последовательностей с одинаковым уровнем сходства построенные SW-выравнивания могут иметь очень различные значения точности и достоверности. Это означает, что в этом диапазоне %ID качество алгоритмических выравниваний определяется не только уровнем сходства сравниваемых белков, но и «внутренними свойствами» их эталонных выравниваний. Эти «внутренние» свойства удобно формулировать в терминах т.н. «островов».

Определение 3.1.1. Островом в выравнивании A = <u, v, Q> называется непродолжаемая последовательность сопоставлений, не разделенных удалениями фрагментов. Весом острова называется сумма весов входящих в остров сопоставлений. Выравнивание можно представить как цепочку островов, разделенных делециями.

Рисунок 3.1.5. (а) Гистограмма распределения количества островов в эталонных выравниваниях (белый) и выравниваниях SW (черный) по весу острова. (б) Суммарная длина островов, имеющих вес в пределах диапазонов. Эталонные острова - белый, SW - черный

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

Рисунок 3.1.6. Более детальное представление данных рис. 3.1.5а. Гистограммы суммарных длин островов с весом в пределах указанных по оси Х диапазонов отдельно для каждой их 3_х областей %ID.Данные по эталонным островам показаны белым, SW - черным

Эталонные выравнивания и SW-выравнивания имеют существенно различную структуру островов с точки зрения веса островов (см. рис. 3.1.5). Неожиданно много эталонных островов имеют очень низкий или даже отрицательный вес, в то время как алгоритмические выравнивания совсем не содержат островов малого веса. Стоит отметить, что суммарная длина таких «слабых» островов в эталонных выравниваниях достаточно велика (см. рис. 3.1.5 б)

Эталонные острова веса меньше 5 составляют 32% от всех островов и покрывают 20% всей длины эталонных выравниваний. Только 5% островов такого малого веса были восстановлены алгоритмом. Для выравниваний из серой зоны (10 < %ID 30) картина еще более критическая - восстановлено всего 2.5% островов веса меньше 5. Эти «слабые» острова обычно не имеют шансов быть восстановленными любым алгоритмом, использующим данную матрицу замен.

С увеличением степени сходства сравниваемых белков (см. рис. 3.1.6) различие в весе эталонных и построенных островов уменьшается. Однако даже для высокого уровня сходства белков (ID > 30%) встречаются эталонные острова отрицательного веса.

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

Представленный в разделе 3.2 алгоритм STRUSWER является модификацией алгоритма Смита-Уотермана. Отличие состоит в том, что при сопоставлении i-го аминокислотного остатка одной последовательности и j-го остатка другой к весу сопоставления добавляется бонус. Этот равен произведению коэффициента SBON, определяющего вклад вторичной структуры в вес сопоставления, на величину сходства элементов вторичной структуры.

Как указывалось выше, предложенный метод может использоваться как с экспериментально полученными, так и с теоретически предсказанными вторичными структурами. Для предсказания вторичной структуры использовалась программа PSIPRED [Bryson K, McGuffin LJ, Marsden RL, Ward JJ, Sodhi JS. & Jones DT. Protein structure prediction servers at University College London. Nucl. Acids Res. 2005. Vol. 33 (Web Server issue): W36-38] в двух режимах: совместного предсказания структуры для группы гомологичных белков (“full version”) и предсказание структуры только по аминокислотной последовательности (“single version”). При этом точность предсказания для использованного набора белков составила соответственно 82% и 65%, что согласуется с результатами, приведенными на сервере EVA (http://cubic.bioc.columbia.edu/eva/).

Для каждой из этих версий использовалось два способа представления предсказанной вторичной структуры. В первом случае («тип_структуры»), каждому остатку аминокислотной последовательности приписывается определенный символ вторичной структуры (H - спираль; Е - бета-структура; L - петля). Во втором («вероятность_структуры»), каждому остатку приписываются вероятности принадлежности остатка к каждому из трех типов вторичной структуры, которые также рассчитываются программой PSIPRED.

При тестировании качество выравниваний, полученных методом STRUSWER с различными способами разметки вторичной структуры (см. выше) сравнивалось с качеством выравниваний, полученных методом SW, а также методом WFMFL, представленный в работе [Wallqvist A, Fukunishi Y, Murphy L.R., Fadel A, Levy R.M. Iterative sequence/secondary structure search for protein homologs: comparison with amino acid sequence alignments and application to fold recognition in genome databases. Bioinformatics. 2000. V. 16. P. 988-1002]

Таблица 3.2.2.

Точность (Асс) и достоверность (Conf) различных методов выравнивания при тестировании на тестовой выборке

Метод	SBON	GOP	GEP	Acc	Conf	ID < 30%
						Acc	Conf
SW	не исп.	7	1	0.525	0.585	0.353	0.429
a) предсказание вторичной структуры по последовательности
STRUSWER_SIN_S	2	10	1	0.578	0.622	0.428	0.482
STRUSWER_SIN_%	7	8	2	0.602	0.618	0.461	0.477
WFMFL_SIN	не исп.	13	1	0.399	0.488	0.263	0.346
б) предсказание вторичной структуры с привлечением данных о гомологичных белках
STRUSWER_PSI_S	8	9	1	0.659	0.683	0.546	0.573
STRUSWER_PSI_%	17	6	2	0.683	0.695	0.579	0.589
WFMFL_PSI	не исп.	16	1	0.631	0.672	0.503	0.56
в) экспериментально известная структура
STRUSWER_EXP	8	10	1	0.677	0.7	0.577	0.601
WFMFL_EXP	не исп.	15	1	0.638	0.698	0.527	0.602

Параметры (SBON, GOP, GEP) были подобраны на обучающей выборке для получения максимальной точности (Асс) каждым из методов. Представлены данные как для всей тестовой выборки (288 пар белков) так и для «серой зоны» (белки с гомологичностью ниже 30%, 182 пары). Для методов STRUSWER и WFMFL указаны способы разметки вторичной структуры: Exp - экспериментально полученная структура; PSI - предсказание структуры по гомологии (“full version” PSIPRED); SIN. Суффиксы _S и _% указывают способ представления предсказаний: «тип_структуры» (_S) или «вероятность_структуры» (_%). Метод WFMFL не приспособлен для использования предсказаний, представленных вероятностями структур

В качестве эталонных выравниваний, как и в разделе 3.1, использовались выравнивания из базы BaliBase (см. выше); корпус эталонных парных выравниваний был разделен на обучающий и тестовый наборы. Источником данных об экспериментально определенных вторичных структурах белков служила база данных DSSP [Kabsch W, Sander C. Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers. 1983. V. 22, P. 2577-2637].

Результаты тестирования представлены в таблице 3.2.2. Параметры каждого из методов, были предварительно подобраны на обучающей выборке выравниваний. Отметим, что параметр SBON имеет смысл только для метода STRUSWER, а параметры аффинных штрафов за удаление (GOP, gap opening penalty, и GEP, gap elongation penalty) используются всеми методами. Значения параметров в таблице 3.2.2 были получены путем максимизации средней точности выравниваний по обучающей выборке. Результаты тестирования с параметрами, полученными при оптимизации достоверности, носят сходный характер.

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

Таким образом, в качестве основного метода предсказания вторичной структуры следует рассматривать метод, делающий предсказание вторичной структуры по отдельной аминокислотной последовательности, в данной работе - программой PSIPREDsingle. Обладая меньшей, чем полный вариант PSIPREDfull, точностью предсказаний, он, тем не менее, имеет ряд преимуществ. Первое преимущество состоит в том, что программа PSIPREDsingle не основана на поиске гомологов, и поэтому STRUSWER в соответствующих режимах (STRUSWER_SIN_S и STRUSWER_SIN_%) можно использовать в случаях, когда произвести поиск гомологов по тем или иным причинам невозможно. Второе преимущество вытекает из первого, и заключается в существенно меньших затратах времени, необходимого для предсказания вторичной структуры. Подобный факт может оказаться решающим в больших вычислительных проектах. Выравнивания, сделанные алгоритмом STRUSWER_SIN, с использованием вторичной структуры, предсказанной по аминокислотной последовательности, превосходят по качеству аналогичные выравнивания, полученные алгоритмом SW и алгоритмом WFMFL_SIN, как по точности, так и по достоверности. Все указанные соотношения остаются верными, если мы ограничимся только слабогомологичными парами белков. При этом относительный выигрыш от использования вторичной структуры существенно возрастает.

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

Определение 3.3.3. РНК-дерево - это такое корневое упорядоченное дерево, что

(1) листья помечены буквами РНК-алфавита {A, C, G, U};

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

РНК-лесом F = <T1, …, Tn > называется упорядоченное множество РНК-деревьев.

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

Представленный в разделе 3.3 алгоритм строит оптимальное выравнивание РНК-лесов при аффинных штрафах за удаление фрагментов РНК. Уточним постановку задачи.

Пусть <S1, P1>, <S2, P2> - РНК-структуры без псевдоузлов, здесь Si - последовательность, Pi - множество спариваний (i=1, 2); G - выравнивание S1 и S2. Будем говорить что в выравнивании G спаривание (x1, y1) P1 сопоставлено спариванию (x2, y2) P2, если x1 сопоставлено x2 и y1 сопоставлено y2.

Весовая система для выравнивания структур РНК - это пятерка <M, g, d, b, c>, где M - весовая матрица замен; g и d - это коэффициенты аффинной весовой функции удалений фрагментов; b - бонус за сопоставление спариваний, c - штраф за каждое спаривание, не участвующее в сопоставлениях.

...