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

Пусть задана весовая система <M, g, d, b, c>; <S1, P1> и <S2, P2> структуры и A - их выравнивание с набором сопоставленных позиций (склеек) {<pk, qk> | k = 1,…, N}. Пусть m - общее число удаленных фрагментов в A; l - суммарная длина этих фрагментов, k - количество спариваний в P1 and P2; которые не сопоставлены в A другим спариваниям: t = (|P1 | + |P2 |-k)/2 - количество сопоставлений спариваний. Тогда структурным весом выравнивания A называется величина

Score(A) = k=1..NM(S1[pk], S2[qk]) - gm - dl + bt - ck. (3.3.1)

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

Представленный алгоритм сводит задачу построения оптимального выравнивания последовательностей РНК с известными вторичными структурами при заданной весовой системе <M, g, d, b, c> к обобщению задачи выравнивания лесов на случай РНК-лесов. Специфика РНК-лесов заключается (1) в различной интерпретации листьев и внутренних вершин (нуклеотиды и спаривания), в частности, в понятии главных сыновей и (2) в использовании аффинной функции штрафов за удаление символов.

Определение 3.3.4. Выравниванием РНК-лесов F1 и F2 называется такое множество A Vertex(F1) Vertex(F2), что для каждой пары (v1,w1), (v2,w2) A выполнено:

1) v1 = v2 тогда и только тогда, когда w1 = w2 (иными словами, A - взаимно однозначное соответствие между подмножествами Vertex(F1) и Vertex(F2)).

2) вершина v1 предшествует вершине v2 в смысле левого обхода леса F1 тогда и только тогда, когда вершина w1 предшествует вершине w2 в смысле левого обхода леса F2.

3) если (v, w) A, то либо обе вершины v и w - листья, либо обе они - внутренние вершины;

4) Пусть (v, w) Vertex(F1) Vertex(F2); v1 и v2 - главные сыновья вершины v; w1 и w2 - главные сыновья вершины w. В этих условиях (v, w) A тогда и только тогда, когда (v1,w1), (v2,w2) A.

Первые два условия традиционны для выравнивания упорядоченных лесов общего вида. Условия (3) и (4) отражают специфику РНК-лесов. Определение веса выравнивания РНК-лесов очевидным образом следует из (3.3.1).

Алгоритм построения оптимального выравнивания РНК-лесов - это алгоритм динамического программирования, использующий парадигму «левый-правый», предложенную в работе Клейна [Klein P.N.. Computing the edit-distance between unrooted ordered trees // Proceedings of the 6th European Symposium on Algorithms(ESA). 1998. P. 91-102]. В указанной работе вес удаления фрагмента пропорционален его длине, что соответствует посимвольным удалениям. Для того, чтобы использовать аффинные штрафы за удаления фрагментов, для каждого промежуточного РНК-леса вычисляются значения дополнительных (по сравнению с алгоритмом Клейна) характеристик.

Определение 3.3.13. Пусть L, R {1, 2} и A - выравнивание лесов F1 и F2. Пусть kL(A) (соответственно, kR(A) ) - количество таких индексов i L (соответственно, i R) , что лес Fi пуст или при выравнивания A в лесе Fi не удален самый левой (соответственно, правый) лист. Делеционным весом WD(A, L, R) выравнивания A при ограничениях L, R называется величина

WD(A, L, R) = W(A) - (kL+ kR)*g,

где W(A) - вес A относительно выбранной весовой системы <M, g, d, b, c>; g -штраф за удаление фрагмента, аналогичный gap opening penalty алгоритма Смита-Уотермана.

Определение 3.3.14. Пусть <F1 ,F2> - пара лесов, L, R {1, 2}. Через Best(F1 ,F2, L, R) обозначается наибольшее возможное значение делеционного веса WD(A, L, R) для выравниваний A лесов F1 и F2.

Очевидно,

WD(A, , ) = W(A)

и Best(F1 ,F2, , ) - это вес оптимального выравнивания лесов F1 и F2.

Сложность работы алгоритма естественно выражается через количество промежуточных пар лесов K, обработанных в ходе работы алгоритма. Потребная память есть O(K), а время - O(Klog(K)). Аналогично алгоритму Клейна, можно показать, что K = O(m2n lg (n)), где m ?n - длины сравниваемых последовательностей. Таким образом, получаем для памяти оценку O(m2n log (n)), а для времени работы - оценку O(m2n log2 (n)).

В разделе 3.4 представлен алгоритм вычисления энергии внутренних петель в структурах РНК в рамках модели Цукера-Тернера-Мэтьюза (Nearest Neighbour Model, NNM) - общепринятой в настоящее время модели строения РНК. Такой алгоритм является необходимой частью общего алгоритма поиска оптимальной (т.е. имеющей минимальную свободную энергию) структуры для данной последовательности РНК в рамках NNM. Как было показано в предыдущих разделах главы 3, использование сведений о вторичной структуре, в том числе - теоретических предсказаний структуры, важно для сравнения биологических последовательностей.

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

Фиксируем последовательность РНК длины L и множество U, состоящее из M спариваний позиций этой последовательности; cпаривания из U будут называться допустимыми. Множество ROWB = {(x, B) U} называется B-й строкой.

Определение 3.4.2. Пусть P U - структура, (A, B), (p, q) P. Спаривания (A, B), (p, q) образуют внутреннюю петлю (internal loop), если A < p < q < B и ни один из нуклеотидов x, где A < x < p или q < x < B, не участвует в спариваниях структуры P. Спаривание (A, B) называется замыкающим спариванием петли; спаривание (p, q) - внутренним спариванием петли; фрагменты [A+1, p-1] и [q+1, B-1] - крыльями петли.

Определение 3.4.4. Будем говорить, что спаривание (A, B) - замыкающее спаривание (closing pairing) в структуре P, если P не содержит пар (x, y), где x < A или y> B.

Алгоритм, который предсказывает вторичную структуру РНК в рамках модели NNM, просматривает все допустимые спаривания строка за строкой в порядке возрастания номеров строк B [1, L]. При этом, для каждого спаривания (A, B) ROWB вычисляется оптимальная структура каждого типа t {0, 1, 2, 3}, тип структуры определяется типом петли, к которой принадлежит замыкающее спаривание структуры. Далее просмотром этих «частных» оптимальных структур вычисляется структура, имеющая наименьшую энергию ДG(A, B) среди всех структур, не имеющих спаренных нуклеотидов вне сегмента [A, B].

Энергия структуры, в которой замыкающая склейка (A, B) принадлежит внутренней петле, образованной склейками (A, B) и (x, y), задается формулой

GInternal_Loop(A, B; x, y) = = fLen(dA+dB) + fDiff(dA-dB) + ДG(x, y) = = fLen((B-A) - (y-x+2)) + fDiff((B+A) - (y+x)) + ДG(x, y). (3.4.1)

Здесь fLen(t) - фиксированная выпуклая функция, она хорошо аппроксимируется логарифмической функцией. Функция fDiff(t) определяется соотношениями:

fDiff (t) = base_level, если |t| width, fDiff (t) = (base_level / width)*|t|, если |t| < width. (3.4.2)

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

GInternal_Loop(A, B) =

= min{GInternal_Loop(A, B; x, y) |(x, y) U & A<x<y<B} =

= min{ fLen((B-A) - (y-x+2)) + fDiff((B+A) - (y+x)) + ДG(x, y)|

|(x, y) U & A<x<y<B}. (3.4.3 )

Говоря неформально, GInternal_Loop описывает увеличение (т.е. ухудшение - структура с большей энергией менее стабильна) энергии структуры РНК в связи с образованием внутренней петли. Это увеличение зависит от длин dA и dB крыльев петли и представляется в виде суммы двух слагаемых. Одно из слагаемых (fLen(t)) зависит от суммарной длины крыльев, другое (fDiff (t)) - от асимметрии петли, т.е. величины | dA - dB |. Если асимметрия мала, то fDiff (t) растет линейно; для значений асимметрии, превышающих порог width, значение fDiff (t) постоянно. Такой своеобразный вид функции GInternal_Loop(A, B; x, y) и представляет основную трудность при анализе внутренних петель.

В первых алгоритмах построения оптимальной вторичной структуры минимум в (3.4.3) находился прямым перебором, что приводило к общему времени анализа внутренних петель O(L4); позднее был предложен алгоритм [Lyngsш et al. Fast evaluation of internal loops in RNA secondary structure prediction //Bioinformatics. 1999. V. 15. P. 440-445] со сложностью O(L3). Представленный в разделе 3.4 алгоритм AFOLD использует технику динамического программирования для разреженных матриц (sparse dynamic programming, см. [Eppstein et al. Sparse Dynamic Programming II: Convex and Concave Cost Functions J. ACM. 1992. V. 39, P. 546-567]) и имеет временную сложность O(M*log2(L)), где L - длина последовательности, M - количество теоретически допустимых спариваний. Т.к. M L2, эта оценка существенно улучшает оценку O(L3) лучшего из ранее предложенных алгоритмов. Отметим, что использованная алгоритмическая техника близка к технике выравнивания последовательностей и может быть использована при построении выравниваний геномов.

В заключение разбора раздела 3.4 сделаем два замечания. Во-первых, в отличие от предлагавшихся ранее подходов, в оценку времени нашего алгоритма явно входит число M допустимых склеек. Благодаря этому алгоритм хорошо сочетается с различными (как экспериментальными, так и теоретическими) способами предварительного отбора (фильтрации) допустимых склеек. Во-вторых, для некоторых классов относительно небольших молекул РНК известно, что их структуры не содержат ветвящихся петель. Для таких молекул наш подход позволяет построить оптимальную структуру за время O(M*log2(L)).

Рассмотренные в главах 2 и 3 алгоритмы оптимального выравнивания имеют квадратичную временную сложность - их время работы (в худшем случае) пропорционально произведению длин последовательностей. Такое время слишком велико для многих биоинформатических приложений, в частности, - для сравнения синтенных фрагментов геномов (длина ~ 107 нуклеотидов. Глава 4 посвящена выравниванию последовательностей, не связанному с оптимизацией весовой функции, и возникающим в этой связи алгоритмическим задачам. В разделе 4.1 представлен иерархический алгоритм геномного выравнивания OWEN, основанный на построении систем коллинеарных локальных сходств. Разделы 4.2 и 4.3 посвящены поиску локальных сходств с помощью т.н. затравок - наиболее распространенному и быстродействующему из известных методов построения локальных сходств. В разделе 4.2 дано определение затравки в наиболее общем виде и рассмотрен простейший класс «неклассических» затравок - разреженные затравки, предложенные в [Ma, B., Tromp, J, and Li, M. PatternHunter: Faster and More Sensitive Homology Search // Bioinformatics. 2002. V. 18. P. 440-445, 2002]. Для этого класса указан вид оптимальной оптимальной затравки среди затравок с одной несущественной позицией. Раздел 4.3 посвящен введенному нами (совместно с Г.Кучеровым и Л. Ноэ) классу затравок - классификационным затравкам, и некоторым алгоритмическим проблемам, связанным с вычислением чувствительности затравок этого класса. Задача вычисления чувствительности произвольных затравок рассмотрена в разделе 5.2 главы 5.

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

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

Есть три основных источника нарушения коллинеарности локальных сходств или микроколлинеарности (термин «микроколлинеарность» используется, чтобы отличить коллинеарность локальных сходств, как кодирующих, так и некодирующих, от макроколлинеарности - коллинеарности генов в синтенных участков геномов). Этими источниками являются (1) локальная конвергентная эволюция, (2) консервативный повтор, присущий обоим геномам, но по-разному расположенный в них; (3) локальные перестановки в геноме.

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

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

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

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

2. Описанный выше принцип применяется как при сравнении последовательностей в целом, так и при сравнении их ортологичных фрагментов («фрактальность»). Благодаря этому, после того, как выделена цепь локальных сходств, значимых на уровне последовательностей в целом, алгоритм может пополнить эту цепь менее сильными сходствами, расположенными в коллинеарных участках между уже выделенными сходствами. Эти «дополнительные» сходства не являлись значимыми при сравнении исходных последовательностей длиной, скажем, 107, но становятся значимыми при сравнении коллинеарных фрагментов длиной, скажем 103. Иными словами, более «сильные» локальные сходства поддерживают коллинеарные с ними более слабые сходства.

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



Рис. 4.1.1. Матрица сходств (dot-matrix), на которой показаны различные вида конфликтов между локальными сходствами. a) Неустранимый конфликт: сходные фрагменты на разных последовательностях расположены в разном порядке. b) Неустранимый конфликт: проекции сходств на последовательность U существенно перекрываются. c) Есть незначительное перекрытие проекций сходств на последовательность U В этом случае конфликт может быть разрешен путем обрезания концов сходств

Рис. 4.1.2. Сравнение синтенных участков последовательностей мыши (AF139987) и человека (AF045555). Наиболее сильное сходство H - это ортологичное сходство, представляющее экзон 10 локусе Rfc2 (нуклеотиды 104514-104583 по последовательности мыши). Каждое из коллинеарных между собой сходств A1, A2 и A3 не является ортологичным. Однако их общий вес выше, чем вес сходства H

Формальное определение остовной цепи основано на понятии p-достоверности (p [0,1] - порог значимости).

Определение 4.1.1. P-значением сходства H между последовательностями U и V называется вероятность того, что независимые случайные последовательности длин |U|, |V| имеют локальное сходства веса не менее Score(H).

Определение 4.1.3. Пусть F - цепь локальных сходств в последовательностях U, V и H - элемент F. Сходство H p-достоверно в F, если его P-значение относительно фрагментов U, V, ограниченных ближайшими к H элементами цепи F, имеющими больший, чем H, вес, не превышает порога p (см. рис. 4.1.3).

Цепь локальных сходств F в последовательностях U, V называется p-достоверной, если каждое, входящее в F сходство p-достоверно в F.

Остовная цепь для U, V с уровнем значимости p - это максимальная p-достоверная цепь для U, V (уточнение понятия максимальности см. определение 4.1.6 в диссертации).

Построение p-остовной цепи для последовательностей U, V выполняется иерархически.

Рис. 4.1.3. Иллюстрация понятия p-достоверности сходства относительно цепи. Пусть цепь F = {A1, H, A2, A3} и Score(H) < Score(A2) < Score(A3) < Score(A1). Тогда (1) FH = {A1, A2, A3}, FA2 = {A1, A3}, FA3 = {A1}, и FA1 пусто; (2) A2 p-достоверно в F, если в выделенном прямоугольнике U[End(A1,U)+1, Beg(A3,U)-1] V[End(A1,V)+1, Beg(A3,V)-1], ограниченном сходствами A1 and A3; (3) A1 p-достоверно в F, если оно p-достоверно относительно последовательностей U and V в целом.

Сначала строятся все сходства, p-достоверные на уровне последовательностей U, V в целом и для этого множества сходств строится максимальная цепь Затем эта процедура рекурсивно применяется к промежуткам между соседними сходствами уже построенной подцепи искомой остовной цепи. Такой иерархический подход позволяет избежать затратного по времени анализа относительно слабых локальных сходств по последовательности в целом, что и определяет эффективность алгоритма.

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

Время работы описанного алгоритма оценивается величиной O(Z) + O(NlogN) + TimeSetLoc, где Z - количество сходств в остовной цепи, N - общее количество локальных сходств, рассмотренных в ходе работы алгоритма, TimeSetLoc - общее время построения локальных сходств. При обычно используемых малых значения параметра p (p <0.01) значение N = k T, где k мало.

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

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

Практически все используемые в настоящее время алгоритмы поиска локальных сходств основаны на фильтрации пространства поиска c помощью предварительного выделения т.н. затравочных сходств. Приведенные ниже определения обобщают определения из работ [Ma, B., Tromp, J, and Li, M. PatternHunter: Faster and More Sensitive Homology Search. // Bioinformatics. 2002. V. 18. P. 440-445], [Brown D. Optimizing multiple seeds for protein homology search. //IEEE Transactions on Computational Biology and Bioinformatics. 2005. Vol. 2 , P. 29 - 38.] и др.

Определение 4.2.1. Затравкой [seed] будем называть произвольное множество локальных выравниваний (=сходств), эти выравнивания называются затравочными выравниваниями. Затравка Z допускает выравнивание G = <v, w, E>, если существует затравочное выравнивание z0 Z такое, что z0 - подвыравнивание выравнивания G.

Простейший пример затравки - это множество точных совпадений заданной длины m. В этом случае затравка состоит из выравниваний вида <u, u, S0>, где u - слово длины m, S0 - тривиальное выравнивание слова с самим собой.

Пусть р - затравка. Назовем слова w', w'' р-эквивалентными, если

?v((v, w') р (v, w”) р).

Очевидно, множество С(р, v) всех слов w, образующих с данным словом v затравочное сходство в смысле затравки р , есть объединение классов р-эквивалентности. На этом наблюдении основан стандартный способ поиска затравочных сходств для данных последовательностей v, w.. Он состоит из двух этапов.

1. Индексирование. Для каждого класса р-эквивалентными Z строится список всех вхождений элементов этого класса в последовательность v (т.н. индексный список).

2. Поиск. Просматриваем последовательность w слева направо. Чтобы найти все фрагменты последовательности v, образующие затравочное сходство с очередным фрагментом w[x, x+d], достаточно просмотреть индексные списки для всех классов р -эквивалентности Z таких, что Z С(р, w[x, x+d).

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

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

Этап 1. Поиск затравок (см. выше).

Этап 2. Поиск в окрестности затравок локальных сходств, представляющих интерес с точки зрения поставленной биологической задачи.

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

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

Вычисление избирательности затравки определяется вероятностью того, что фрагменты независимых случайных последовательностей, которые начинаются в заданных позициях, образуют затравочное сходство. Например, для рассмотренной выше затравки «точное совпадение длины n» и бернуллиевских случайных последовательностей избирательность равна pn. В то же время вычисление чувствительности затравки представляет определенную сложность. Вычисление чувствительности затравки означает вычисление вероятности (в рамках заданной вероятностной модели) того, что случайное выравнивание из заданного множества целевых выравниваний содержит подвыравнивание, принадлежащее затравке. Алгоритмы вычисления чувствительности затравок рассмотрены в разделе 5.2 главы 5.

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

Первый из этих классов (см. раздел 4.2) - это разреженные затравки (spaced seeds) - первый класс затравок, отличный от «классических» идущих подряд n совпадений; этот класс был предложен в начале 2000-х годов (см. [Burkhardt S. and Karkkainen J. Better Filtering with Gapped q-Grams // Fundamenta Informaticae. 2003. Vol. 56, N. 1-2. P. 51-70], [Ma, B., Tromp, J, and Li, M. PatternHunter: Faster and More Sensitive Homology Search // Bioinformatics. 2002. V. 18. P. 440-445]. Неформально говоря, разреженная затравка, как и классическая затравка, требует наличия нескольких совпадений, но расположенных не подряд, а по описанной в затравке схеме. Более формально разреженные затравки описываются следующими определениями.

Определение 4.2.3. Разреженным термом E называется список положительных целых чисел {p1, …,pd} таких, что p1< …<pd и p1=1. Размером терма E называется величина span(E) = pd; весом терма E называется число weight(E) = d. Отметим, что вес разреженной затравки напрямую связан с ее избирательностью.

Терм E ={p1, …,pd} часто для наглядности представляют словом tE длины pd в двухбуквенном алфавите {#, -} таким, что tE[i] = `#' i {p1, …,pd}. Символ `-` называют джокером.

Определение 4.2.5. Пусть G = безделеционное выравнивание слов v и w; E ={p1, …,pd} - разреженный терм. Выравнивание G соответствует терму E, если

1) |v|=|w\ = pd;

2) ?j {1, …d}( v[k+ pj+1] = w[k+ pj+1]).

Разреженная затравка, задаваемая термом E - это множество всех выравниваний, соответствующих терму E.

Т.к. разреженные затравки трактуют выравнивания только с точки зрения совпадений-несовпадений, то при работе с разреженными затравками выравнивания обычно перекодируются в двухбуквенный алфавит с буквами 1 (совпадение) и 0 (несовпадение).

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

Определение 4.2.9. Выравнивание g{1, 0}m называется (m, k)-выравниванием, если в нем есть ровно k нулей и m-k единиц. Затравка E решает (m, k)-проблему, если она допускает все (m, k)-выравнивания.

Определение 4.2.11. Пусть Q - затравка и k - допустимое количество несовпадений. k-критическая длина для Q - это минимальная длина m, при которой Q решает (m, k)-проблему. Затравка Q называется k-оптимальной в некотором классе затравок, если она имеет минимальную k-критическую длину в этом классе.

Утверждение 4.2.2. Рассмотрим класс затравок веса не менее d с одним джокером, k - натуральное число. Тогда k-оптимальной затравкой будет затравка #d-r-#r, где

r = [ d/2] , если k=1,

r = [ d/3] , если k> 1,

где [x] - ближайшее целое к числу x и [n+1/2] = n.

Этот результат был обобщен на случай затравок с несколькими джокерами в [Farach-Colton, M., Landau, G.M., Sahinalp, S.C., Tsur, D: Optimal spaced seeds for faster approximate string matching // J. Comput. Syst. Sci. 2007. V. 73. P. 1035-1044].

Пример 4.2.1. Затравка #4-#2 - оптимальная в смысле приведенного выше определения среди затравок веса 6 с одним джокером. Это, в частности, означает, что она решает (m, 2) - проблему для всех m ? 16 и это - наилучшая оценка в рассматриваемом классе. Аналогично, эта затравка решает (m, 3) - проблему для всех m ? 20 и это - наилучшая оценка в рассматриваемом классе и т.д.

Естественным обобщением разреженных затравок являются классификационные затравки, представленные в разделе 4.3.

Пусть р - затравка. Слова равной длины w' и w'' называются р-эквивалентными, если для любого слова v той же длины безделеционные выравнивания (v, w') и (v, w”) либо одновременно являются затравочными сходствами для р (см. определение 4.2.1), либо одновременно не являются. Через С(р, v) обозначим множество всех слов w, образующих с данным словом v затравочное сходство в смысле затравки р. Очевидно, множество С(р, v) есть объединение классов р-эквивалентности.

Для произвольной разреженной затравки E каждый класс С(E, v) состоит ровно из одного класса E-эквивалентности (говорят, что разреженные затравки допускают однозначное индексирование), что упрощает поиск соответствующих затравочных сходств. Возможность однозначного индексирования для разреженных затравок связана с тем, что эти затравки различают только два вида отношений между символами - совпадение и несовпадение. В [Brejova B., Brown, D., Vinar T. Vector seeds: an extension to spaced seeds allows substantial improvements in sensitivity and specificity // Proceedings of the 3rd International Workshop in Algorithms in Bioinformatics / Eds. Benson, G., Page, R. Lecture Notes in Computer Science. Vol. 812. Springer. 2003] была предложена более гибкая и общая модель затравок - векторные затравки, в которых критерий соответствия между затравкой и выравниванием базируется на интегральной характеристике - сумме вкладов разных позиций. Платой за эту общность является то, что для векторной затравки R класс С(R, v) может состоять из десятков классов R-эквивалентности, что усложняет как поиск, так и индексирование.

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

Фиксируем алфавит последовательностей П. Пусть M = {(a, a)|a П} - множество всех пар-совпадений.

Определение. 4.2.12. Классификационный алфавит B - это алфавит, каждая буква b B которого соответствует непустому подмножеству Пb П П, причем

1) ?b B (M Пb ) (любая классификационная буква допускает совпадение);

2) В B есть буква #, для которой П# = M (# требует только совпадений, она ниже называется единичной буквой).

Очевидно, алфавит разреженных затравок {#, -} - классификационный, причем джокеру соответствует множество П П.

Определение. 4.3.1. Классификационная затравка - это слово в некотором классификационном алфавите B. Затравка р = b1…bm Bm допускает фрагмент выравнивания a1…am (П П)m, если ?i {1, …, m} ai bi. Размером span(р) затравки р называется ее длина m, #-весом - количество sharp(р) букв `#' среди b1, …, bm

Определение. 4.3.2. Назовем буквы x, y П П эквивалентными относительно алфавита B (B-эквивалентными), если ?bB (x Пb y Пb).

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

Пример 4.3.1. Рассмотрим выравнивания нуклеотидных последовательностей, т.е. выравнивания над алфавитом последовательностей {A, C, G, T}. Транзицией называется одна из четырех пар {(A, T), (T, A), (C, G), (G, C) }. Трансверсией называется любое несовпадение, отличное от транзиции. Известно, что среди замен в нуклеотидных выравниваниях частота транзиций существенно выше, чем частота трансверсий. В качестве алфавита затравок возьмем алфавит B = { #, @, -}; элементы B представляют соответственно только совпадение (#), совпадение или транзицию (@), любое сопоставление (-). Очевидно, порожденный алфавитом B алфавит выравниваний A включает три множества: (1) совпадения (символ: 1); (2) транзиции (символ: h); (3) трансверсии, т.е. несовпадения, отличные от транзиций (символ:0).

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

Пусть B - алфавит классификационных затравок. Будем считать, что каждой букве c B сопоставлена вероятность в(c). Это можно сделать, например, так. Рассмотрим алфавит последовательностей П, связанный с алфавитом B. Пусть на П задано распределение вероятностей в: П -> [0,1], где

Уx П в(x) = 1.

Это распределение индуцирует распределение вероятностей на множестве пар П П по правилу в(x, y) = в(x) в(y). Так как буквы классификационного алфавита соответствуют подмножествам в П П, то, тем самым, зная распределение в, каждому классификационному символу c можно приписать вероятность в(c). Отметим, что, вообще говоря,

Уc B в(c) 1.

Определение 4.3.3. Пусть B - алфавит классификационных затравок и каждой букве c B сопоставлена вероятность в(c) [0, 1]. Пусть # - единичный символ алфавита B, т.е. символ, соответствующий множеству совпадений M = {(x, x)|x П}. Весом произвольного символа c B называется величина

weight(c) = -log(в(c))/log(в(#)).

Весом weight(р) классификационной затравки р называется сумма весов входящих в нее букв.

Пример 4.3.2 (продолжение примера 4.3.1). Пусть на алфавите {A, C, G, T} задано равномерное распределение вероятностей. Тогда weight(#) = 1; weight(@) = Ѕ; weight(-) = 0. Вес weight(р) затравки р = #@-# равен 2.5.

В п.4.3.4 диссертации показано (см. таблицы 4.3.3, 4.3.4), что при фиксированном весе, а следовательно, при фиксированной избирательности, классификационные затравки из рассмотренных выше примеров 4.3.1, 4.3.2 имеют более высокую чувствительность, чем разреженные затравки. В обеих таблицах множество целевых выравниваний - это множество всех безделеционных выравниваний длины 64. В таблице 4.3.1 распределение вероятностей на этом множестве - это бернуллиевское распределение, в таблице 4.3.4 - марковское распределение порядка 3. Параметры этих распределений получены из анализа парных геномных выравниваний для 40 бактериальных геномов. Ниже приведена таблица 4.3.3, данные таблицы 4.3.4 аналогичны.

Таблица 4.3.3.

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

Вес	Разреженные затравки	Чувстви-тельность	Классификационные затравки, два @	Чувстви-тельность
9	##-##-#-#-###	0.4183	###-#-#@#-@##	0.4443
10	##-##-##-#-###	0.2876	###-@#-@#-#-###	0.3077
11	###-#-#-###-###	0.1906	##@#-##-#-#-@###	0.2056
12	###-#-##-#-##-###	0.1375	##@#-#-##-#@-####	0.1481

Оптимальные затравки в классах разреженных затравок и классификационных затравок в алфавите {#, @, -} (см. пример 4.2.2) с не более, чем двумя символами @ для различных весов (от 9 до 12). Целевые выравнивания имеют длину 64. Распределение вероятностей - бернуллиевское, все символы равновероятны.

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

Утверждение 4.3.3. Пусть р = b1…bm Bm - классификационная затравка длины m, w - количество символов # в р и Sр = < A, Q, q0, QF, ш> - классификационный автомат затравки р. Тогда количество состояний автомата Sр не превосходит (w+1)*2m-w.

Эта оценка не улучшаема, что показывает следующее

Утверждение 4.3.4. Пусть A = {#, -} и р = #-r# (между символами # расположено r джокеров). Тогда классификационный автомат Sр неприводим, т.е.

(i) любое состояние Sр достижимо из начального состояния q0;

(ii) никакие два состояния не эквивалентны.

В силу определения 4.2.1, классификационный автомат затравки р эквивалентен автомату Ахо-Корасик для множества слов Mр, соответствующих затравке р. Для двухбуквенного алфавита выравниваний, соответствующего разреженным затравкам, оценка утверждения 4.3.3 совпадает с оценкой количества состояний автомата Ахо-Корасик из работы [Aho, A.V.; Corasick, M.J. Efficient String Matching. // Communications of the ACM. 1975. V. 18. P. 333-340.], которая дается формулой (w+1)*am-w, где a - размер алфавита выравниваний. Таким образом, утверждение 4.3.3 дает верхнюю оценку для числа состояний приведенного автомата, соответствующего автомату Ахо-Корасик для множества слов Mр.

Таблица 4.3.1 представляет данные о среднем количестве состояний классификационного автомата соответственно для разреженных и классификационных затравок из примера 4.3.1. Для каждого веса w было подсчитано среднее количество состояний у автомата Ахо-Корасика (столбцы «Ахо-Корасик») и классификационного автомата, описанного в разделе 4.2.4 (столбцы «Классиф. автомат»). В каждом случае показывалось как среднее количество состояний (столбцы «Разм.»), так и отношение этого среднего к среднему количеству состояний для соответствующего минимального автомата (столбцы «Отн»). Среднее количество состояний дли минимального автомата (это число одинаково для обеих конструкций, т.к. соответствующие автоматы эквивалентны) приведено в столбце «Разм. мин. авт.». Отметим, что классификационный автомат компактнее автомата Ахо-Корасик не только для 3-буквенного алфавита B (см. выше), но и для двухбуквенного алфавита R. При данном весе для регулярных затравок среднее подсчитано по всем затравкам веса w и длины не более w+8; для затравок примера 4.3.1 - по всем затравкам веса w, длины не более w+5 и содержащим не более двух символов @.

Таблица 4.3.1.

Классификационные автоматы и автоматы Ахо-Корасик.

Вес	Разреженные затравки {#, -}	Классификационные затравки {#, @, -}
	Ахо-Корасик	Классиф. автомат	Размер мин. авт.	Ахо-Корасик	Классиф. автомат	Размер мин. авт
	Разм.	Отн.	Разм.	Отн.		Разм.	Отн.	Разм.	Отн.
9	345.9	3.1	146.3	1.3	113.2	1900.7	16.0	167.6	1.4	119.0
10	380.9	3.2	155.1	1.3	120.6	2104.0	16.5	177.9	1.4	127.5
11	415.4	3.3	163.8	1.3	127.6	2306.3	17.0	188.1	1.4	136.0
12	449.5	3.3	172.4	1.3	134.9	2507.9	17.4	198.1	1.4	144.0
13	483.3	3.4	180.9	1.3	141.8	2709.0	17.8	208.1	1.4	152.3

В главах 2-4 представлены алгоритмы построения парного выравнивания биологических последовательностей, ориентированные на широкий спектр постановок задач. Большинство из этих алгоритмов основано на методе динамического программирования. Как известно (см. [Finkelstein, A.V. and Roytberg, M.A. Computation of biopolymers: a general approach to different problems. // BioSystems. 1993. Vol. 30, P.1-19.]), динамическое программирование, наряду с построением оптимальных объектов, позволяет решать, в определенном смысле, двойственную задачу. Это задача вычисления т.н. обобщенных статистических сумм, которые характеризуют множество допустимых объектов в целом. Более формально, речь идет о следующей задаче.

Дан конечный ациклический ориентированный граф G,для простоты будем считать, что у графа G один источник A и один сток Z. Полным путем в графе G называется путь из источника в сток. Ребра графа помечены элементами множества R, на котором определены две ассоциативные операции - «сложение» (обозначается ) и «умножение» (обозначается ), обладающие следующими свойствами:

1) нейтральные элементы 0 и 1:

?x M (0 x) =x; ?x M (1 x) = x;

2) коммутативность сложения:

?x, y M (y x = x y);

3) дистрибутивность:

?x, y, z M ( (x y) z = (xz) (yz) ) ;

Весом пути называется «произведение» (относительно операции ) весов входящих в него ребер, взятых в порядке следования от источника к стоку.

Обобщенной статистической суммой (ОСС) графа G называется «сумма» (относительно операции ) весов всех его полных путей.

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

Значение обобщенной статистической суммы может быть вычислено за время O(Edge(G)), где Edge(G) - количество ребер в графе G.

Глава 5, заключительная глава диссертации, посвящена применению обобщенных статистических сумм в задачах биоинформатики. Наиболее распространенной интегральной характеристикой множества объектов является его вероятность. В разделе 5.1.1 приведены достаточные условия, при которых вычисление вероятности множества символьных последовательностей может быть сведено к вычислению обобщенной статистической суммы. Этот подход применен к вычислению чувствительностей затравок (см. главу 4 и раздел 5.2) и вероятностей обнаружения сигналов в последовательностях (раздел 5.3).

Фиксируем алфавит A и длину последовательностей m. Пусть задано множество последовательностей S Am и распределение вероятностей с на множестве Am. В разделе 5.1 вычисление вероятности Pс(S) сведено к вычислению обобщенной статистической суммы для некоторого графа. Этот граф строится на основе конечного детерминированного автомата (КДА), распознающего язык S и конечно-автоматного генератора вероятностей, задающего распределение с.

Множество S Am, вероятность которого требуется найти - конечное, и, следовательно, допускается конечным детерминированным автоматом K с не более, чем m|S| состояниями. Специальные случаи построения такого автомата K, распознающего множество S рассматриваются ниже. «Конечно-автоматность» распределения с является достаточным условием применимости предлагаемого подхода, она эквивалентна существованию НММ, описывающей распределение с.

Говоря неформально, конечно-автоматный генератор вероятностей - это недетерминированный автомат без заключительных состояний, в котором для каждого перехода (q, a, q'), где q - состояние автомата, a - символ используемого алфавита, q' - новое состояние, задана соответствующая вероятность. Соответствие между конечно-автоматными генераторами вероятностей и НММ аналогично соответствию между автоматами Мили и Мура.

Определение 5.1.1. Пусть A - алфавит. Конечно-автоматный генератор вероятностей над алфавитом A - это четверка G= <Q, q0, A, с >, где Q - множество состояний; q0- начальное состояние; A - алфавит, с: Q A Q -> [0, 1] - функция генерации вероятностей такая, что

ЃНq Ѓё Q У qЎЗЎфQ,aЎфA с(q, a, q') = 1.

Переходом в G называется тройка e =<q,a,q'> такая, что с(q, a, q') >0. Буква a называется меткой перехода и обозначается label(e). Состояния q и q' называются соответственно начальным и конечным состоянием перехода и обозначаются соответственно start(e), end(e). Меткой пути P = (e1, …, en) в G называется слово label(e1)…label(en).

Генератор G называется детерминированным, если для любой пары q Ѓё Q, aЃёA есть не более одного перехода вида <q,a,q'>.

Путь P называется инициальным, если его начальное состояние - это состояние q0 генератора G. Вероятностью с(P) пути P = (e1, …, en) называется произведение с(P) = Р i=1, n с(ei). Вероятностью слова w относительно генератора G называется величина PG(w), равная сумме вероятностей всех инициальных путей с меткой w. Вероятностью конечного набора слов (языка) L A* называется сумма PG(L) вероятностей всех слов w ? L. Очевидно, для любого n выполнено PG(An) = 1.

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

Пусть даны КДА K = <QK, q0K, QFK, A, цK>, - распознающий язык S и конечно-автоматный генератор G= <QG, q0G, A, сG >, задающий распределение вероятностей с на Am. Вычисление вероятности PG(L) сводится к вычислению ОСС для графа т.н. вероятностного распознающего автомата (ВР-автомата), который строится как декартово произведение K и G. Говоря неформально, ВР-автомат - это недетерминированный конечно-автоматный генератор вероятностей, в котором выделено подмножество допускающих состояний.

Определение 5.1.6. Пусть K = <QK, qK0, QFK , A, ц > - детерминированный конечный автомат без выхода; G= <QG, q0G, A, с > - конечно-автоматный генератор вероятностей над алфавитом A. ВР-автомат W = KG - это пятерка W= =<QW, q0W, QFW, A, , сW>, где QW = QK QG ; q0W = q0K q0G ;

QFW = {<k, g>= QK QG | k QFK };

сW(<k, g>, a, <k', g'>) = с( g, a, g'), если ц(k, . a) = k';

сW(<k, g>, a, <k', g'>) = 0 - в противном случае.

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

1. Значение вероятности PG(L) равно значению обобщенной статистической суммы для графа ВР-автомата W.

2. Автомат W имеет |QG|*|QK| состояний и для каждой пары, состоящей из состояния q автомата W и символа a A есть не более |QG| исходящих из q ребер.

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

Утверждение 5.1.4. Рассмотрим алфавит A и конечное множество L Am. Пусть K = <QK, q0K, QFK, A, цK> - КДА, распознающий множество L и G= <QG, q0G, A, сG > - конечно-автоматный генератор, задающий распределение вероятностей с на Am. Тогда вероятность PG(L) множества L относительно распределения сG может быть найдена за время O(|QG|2*|QK|*|A|) с использованием памяти O(|QG|*|QK|). Если генератор G детерминированный, то для времени работы алгоритма верна оценка O(|QG|*|QK|*|A|).

Следствие 1. Пусть распределение вероятностей на множестве Am - бернуллиевское. Тогда для времени и памяти вычисления вероятности PG(L) верны оценки TimeBern? O(|QK|*|A|) и SpaceBern? O(|QK|)

Следствие 2. Пусть распределение вероятностей на множестве Am - это Марковское распределение порядка r. Тогда для времени и памяти вычисления вероятности PG(L) верны оценки TimeMarkov? O(|QK|*|Ar+1|) и SpaceMarkov? O(|QK|*|Ar|)

Утверждение 5.1.6. Рассмотрим алфавит A и конечное множество L Am. Пусть B = <QB, q0B, QFB, A, цB> - автомат, распознающий множество L' такое, что L' Am = L и G= <QG, q0G, A, сG > - конечно-автоматный генератор, задающий распределение вероятностей с на Am. Тогда вероятность Pс(L) множества L относительно распределения с может быть найдена за время O((|QB|*|QG|2*|A|*m) с использованием памяти O(|QB|*|QG|). В случае детерминированного генератора G время работы равно O((|QB|*|QG|*|A|*m).

Следствие 1. Пусть распределение вероятностей на множестве Am - Бернуллиевское. Тогда для времени и памяти вычисления вероятности PG(L) верны оценки TimeBern? O(|QB|*|A|*m) и SpaceBern? O(|QB|)

Следствие 2. Пусть распределение вероятностей на множестве Am - это Марковское распределение порядка r. Тогда для времени и памяти вычисления вероятности PG(L) верны оценки TimeMarkov? O(|QB|*|Ar+1|*m) и SpaceMarkov? O(|QB|*|Ar|)

Утверждение 5.1.8. Рассмотрим алфавит A и конечное множество L Am. Пусть M - конечное множество слов, множество L' - это множество

L' = {w A*| w содержит подслово из M}

и при этом выполнено L = Am L'.

Пусть далее, B = <QB, q0B, QFB, A, цB> - автомат Ахо-Корасик для множества M, распознающий множество L' и с - распределение вероятностей на Am, задаваемое марковской моделью порядка r. Тогда вероятность Pс(L) множества L относительно распределения с может быть найдена за время O((|QB|+|Ar|)*|A|*m) с использованием памяти O((|QB|+|Ar|)).

В разделе 5.2 результаты предыдущего раздела применены для вычисления чувствительности данной затравки, т.е. вероятности (в рамках заданной вероятностной модели) того, что случайное выравнивание из заданного множества целевых выравниваний допускается затравкой. Все выравнивания в этом разделе считаются безделеционными и представляются словами в некотором алфавите A. В качестве множества целевых выравниваний рассматривается множество Am всех выравниваний длины m. На Am задано распределение вероятности с, отражающее вероятности замен в интересующем исследователя классе выравниваний. Как и в разделе 5.1, будем считать, что распределение вероятности с задано конечно-автоматным генератором вероятностей G= <QG, q0G, A, с >. Фиксируем затравку р.

Определение 5.2.1 Через S(р) A* обозначается множество всех выравниваний, допускаемых затравкой р. Автоматом затравки р называется детерминированный конечный автомат без выхода, распознающий множество S(р). Через S(р, m) Am обозначается множество всех выравниваний, допускаемых затравкой р

...