Смиряев а. В., Панкина л. К. Основы биоинформатики

Самый простой вид штрафа – так называемый линейный штраф, пропорциональный длине разрыва:

где g – длина разрыва; d – штраф за одиночный разрыв.

Другой вид, так называемый аффинный штраф за разрыв определяют по формуле:

R(g)= – d – (g – 1)e,

где g – длина разрыва; d – штраф за открытие разрыва; е – штраф за его продолжение.

Обычно штраф за продолжение разрыва (е) меньше штрафа за открытие (d); тогда длинные вставки и делеции аффинной функцией штрафа наказываются меньше, чем линейной. Это желательно, когда известно заранее, что ожидаемая частота разрывов в один и несколько остатков примерно одинакова. Типичные значения штрафов за разрывы, используемые на практике, равны d = 8 для линейного штрафа, или d = 12, e = 2 для аффинного случая.
3.4. Алгоритмы выравнивания.

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

глобальных выравниваний двух последовательностей длиной n. Например, при n=9 N≈10⁹. Перебрать все варианты выравниваний для выявления оптимального физически невозможно, даже при небольших значениях п. Поэтому для нахождения оптимального выравнивания были разработаны специальные алгоритмы.

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

Для того чтобы проиллюстрировать работу разных методов выравнивания, возьмем две короткие аминокислотные последовательности, x) HEAGAWGHEE и y) PAWHEAE. Для вычисления веса выравнивания будем использовать матрицу BLOSUM50 (табл. 3.2). Воспользуемся формулой (3.1) – линейной функцией штрафа за разрывы. Штраф за удаление одного остатка примем d = 8. В таблице 3.3 две тестовые последовательности x и y изображены так, чтобы показать соответствующие веса матрицы BLOSUM50 для всех возможных пар аминокислот этих последовательностей. Веса выравнивания идентичных (консервативных) остатков выделены жирным шрифтом.

y	H	E	A	G	A	W	G	H	E	E
P	–2	–1	–1	–2	–1	–4	–2	–2	–1	–1
A	–2	–1	5	0	5	–3	0	–2	–1	–1
W	–3	–3	–3	–3	–3	15	–3	–3	–3	–3
H	10	0	–2	–2	–2	–3	–2	10	0	0
E	0	6	–1	–3	–1	–3	–3	0	6	6
A	–2	–1	5	0	5	–3	0	–2	–1	–1
E	0	6	–1	–3	–1	–3	–3	0	6	6

Рассмотрим построение оптимального глобального выравнивания двух последовательностей x – длиной n и y – длиной m. Проблема пошаговой процедуры выравнивания в том, что на каждом ее «шаге» накопленная сумма весов зависит от выбранных вариантов (замены, вставки – разрывы) на предыдущих «шагах». Поэтому даже при небольших значениях n и m прямой перебор и сравнение эффективности всех сочетаний вариантов на всех «шагах» общего выравнивания затруднен – слишком велико число сочетаний. Алгоритм решения этой задачи основан на методе динамического программирования. Его суть в следующем. Для каждого сочетания i, j длин коротких фрагментов последовательностей x и y, начиная с i = 1, j = 1, подбирается условно оптимальное выравнивание – с максимальным весом. Условно оптимальное потому, что рассматриваются все пары длин i, j, среди которых наверняка будут лишние – они не войдут в полное оптимальное выравнивание последовательностей x и y. Всего таких условно оптимальных выравниваний n x m. То есть алгоритм динамического программирования состоит в пошаговом построении оптимального выравнивания большей длины, используя полученные на предыдущем «шаге» условно оптимальные выравнивания фрагментов меньшей длины.

Для этого строится матрица F размера (n+1) x (m+1). Элемент F(i, j) этой матрицы содержит вес наилучшего условного выравнивания между фрагментом x_1...i длины i последовательности х и фрагментом y_1...j длины j последовательности у.

Например, примем F(0, 0) = 0 для исходного (нулевого) шага алгоритма.

Заполняется матрица с верхнего левого угла к нижнему правому. Для того чтобы определить значение F(i, j) для очередного шага необходимо рассмотреть квадрат из четырех ячеек (рис. 3.1).

Рис 3.1. Квадрат из трех «предыдущих» и одной «последующей» ячейки для определения значения F(i, j).

Рассматриваем каждый квадрат из 4–х ячеек и определяем «последующее» значение F(i, j), применяя уравнение (3.2) пока вся матрица F(i, j) не будет заполнена. После каждого вычисления максимального значения F(i , j) ставим только одну стрелку из трех, обозначенных на рис 3.1, направленную в ту «предыдущую» ячейку, которую использовали для получения этого значения F(i , j). Среди прочих будет правая нижняя ячейка матрицы для i=n, j=m. Наилучшим весом всего оптимального выравнивания x и y является значение F(n, m) в этой ячейке.

Для того чтобы прочитать само оптимальное выравнивание, необходимо восстановить последовательность выбора вариантов (замен, вставок – разрывов), которая привела к максимальному весу F(n, m). Процедура восстановления выборов называется процедурой обратного прохода. Она осуществляет восстановление выравнивания, от правой нижней ячейки матрицы (n, m), двигаясь по шагам строго по указателям – стрелкам, которые были оставлены при построении матрицы F(i, j). На каждом шаге обратного прохода добавляют пару символов аминокислот или «–» слева к текущему выравниванию. А именно добавляем:

аминокислоты x_i и y_j (из i–ой и j–ой позиций x и y), если стрелка указывает на ячейку (i – 1, j – 1);

«–» и y_j, если – на (i, j – 1).

В конце концов, достигнем левого верхнего угла матрицы, где i=j=0. Оптимальное глобальное выравнивание построено.

Рассмотрим пример. Построим глобальное выравнивание для последовательностей x) HEAGAWGHEE и y) PAWHEAE. Линейный штраф за удаление одного остатка (разрыв) d = 8.

Построим матрицу F(i, j) динамического программирования (рис 3.2).


x

y

1

42

Рис. 3.2. Матрица F(i,j), динамического программирования глобального выравнивания двух модельных последовательностей x и y. Под матрицей – полученное оптимальное выравнивание с весом 1 (по Р Дурбин и др.,2006).

		i=0	i=1
		x	H
j=0	y	0	–8
j=1	P	–8	F(1, 1)

Определим значение F(i, j)= F(1, 1).

Для нашего случая F(i – 1, j – 1)= F(0,0)=0;

S(1, 1)= S(Н, Р)= –2 (из матрицы табл. 3.3);

F(i – 1, j)= F(0, 1)= –8;

F(i, j – 1)= F(1, 0)= –8.

Тогда, используя формулу 3.2, получаем:

Рассмотрим следующий квадрат из четырех ячеек.

		i=0	i=1	i=2
		x	H	E
j=0	y	0	–8	–16
j=1	P	–8	–2	F(2, 1)

Определим значение F(2, 1).

Значение F(i – 1, j – 1)= F(1, 0)= –8;

S(2, 1)= S(Е, Р)= –1 (из матрицы табл. 3.3);

F(i – 1, j)= F(1, 1)= –2;

F(i, j – 1)= F(2, 0)= –16.

Применяя формулу 3.2, получаем:

Максимальное значение F(2, 1) равно –9.

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

На рис.3.2 показана матрица динамического программирования глобального выравнивания двух последовательностей с указателями (стрелки выделены жирным шрифтом) для процедуры обратного прохода. Значения F(i, j), соответствующие оптимальному выравниванию (на каждом оптимальном шаге) выделены жирным шрифтом. Оптимальное глобальное выравнивание имеет вес 1 (на рис. 3.2 в кружке). Оно выстраивается справа налево, читая по шагам указания жирных стрелок. Результат указан под матрицей.

Процедура обратного прохода, описанная здесь, находит только одно выравнивание с наилучшим весом; но если в некий момент (шаг) максимум величины F(i, j) достигается двумя или даже тремя способами (это в принципе возможно), то делается случайный выбор оптимального шага. В такой ситуации оптимальные выравнивания могут отличаться, но их вес при этом не изменится.
Локальное выравнивание: алгоритм Смита–Уотермана.

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

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

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

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


Выбор нуля соответствует началу нового локального выравнивания, начиная с пары остатков с номерами i, j двух последовательностей x и y. В алгоритме заложено правило перехода к новому локальному выравниванию. Например, если наилучшее выравнивание на очередном шаге приобретает отрицательный вес, лучше начать новое выравнивание, чем продолжать старое. Тогда, в частности, из-за гарантированно отрицательных значений F(–id, 0) и F(–jd, 0) матрицы глобального выравнивания (включение разрывов), элементы верхней строчки и левого столбца таблицы F(i, j) локального выравнивания определяют равными нулю (рис.3.3.):

F(0, 0) = 0; F(i, 0)= 0, a F(0, j)= 0

Второе отличие локального по сравнению с глобальным выравниванием состоит в том, что в первом случае выравнивание может заканчиваться в любом месте таблицы F(i, j). Наилучшим весом локального выравнивания считается наибольшее значение весов всей матрицы F(i, j). Для использованных выше двух последовательностей x и y максимальным значением веса во всей матрице является 28 (помечено кружком на рис.3.3.). Стандартную процедуру обратного прохода нужно начинать с F(9, 5), т.е. с аминокислот Е, Е последовательностей x и y. Обратный проход заканчивается, когда встречаем нулевой элемент таблицы (помечен квадратом), что соответствует началу лучшего локального выравнивания.

Рис. 3.3 Матрица динамического программирования для поиска оптимального локального выравнивания последовательностей x и y (по Р Дурбин и др.,2006).

Отдельное от первого «менее оптимальное» локальное выравнивание длиной 3 можно выделить двигаясь от F(3, 6) до F(1, 4). Его вес равен 21.