Алгоритмы биоинформатики ФББ 2013 г., весенний семестр, 3-й курс. Миронов Андрей Александрович Спирин Сергей Александрович. - презентация

1 Алгоритмы биоинформатики ФББ 2013 г., весенний семестр, 3-й курс. Миронов Андрей Александрович Спирин Сергей Александрович

2 Информатика и Биоинформатика Биологическая задача Формализация Алгоритм Тестирование Параметры Определение области применимости

3 Пример: сравнение последовательностей Тестирование: алгоритм должен распознавать последовательности, для которых известно, что они биологически (структурно и/или функционально) сходныТестирование: алгоритм должен распознавать последовательности, для которых известно, что они биологически (структурно и/или функционально) сходны

4 Сравнение последовательностей Формализация1: глобальное выравнивание Алгоритм1: Граф выравнивания, динамическое программирование Алгоритм1а: Граф выравнивания, динамическое программирование, линейная память Параметры: Матрица сходства, штраф за делецию

5 Сравнение последовательностей Формализация2: локальное выравнивание Алгоритм2: Граф локального выравнивания, динамическое программирование Параметры: Матрица сходства, штраф за делецию

6 Сравнение последовательностей Формализация3: локальное выравнивание с аффинными штрафами Алгоритм3: Расширенный граф локального выравнивания, динамическое программирование Параметры: Матрица сходства, штраф за открытие делеции, штраф за расширение делеции

7 Сравнение последовательностей Алгоритм4: BLAST. Формальная задача плохо определена. Параметры: Размер якоря, матрица сходства, штраф за делецию

8 Выравнивания

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

10 Сколько существует выравниваний? Дано: две последовательности S 1 и S 2 длиной m и n. Сколько есть способов написать одну последовательность под другой (со вставками)? Построим выборочную последовательность S длиной m+n следующим образом: возьмем несколько символов из последовательности S 1, потом несколько символов из последовательности S 2 потом опять несколько символов из S 1, потом опять несколько из S 2. (Доказать!) –Каждой выборочной последовательности S соответствует выравнивание и по каждому выравниванию можно построить выборочную последовательность. (Доказать!) (Доказать!) –Количество выборочных последовательностей равно N sel = C n+m m = (m + n)! / (m!n!) (Доказать!)

11 Динамическое программирование для редакционного расстояния Граф редакционного расстояния для последователь- ностей S 1,S 2 : вершина v i,j соответствует префиксам последовательностей {S 1 1..i }, {S 2 1..j }. На вершине записано редакционное расстояние между префиксами. (красные стрелки соответствуют вставкам и удалениям) d i,j d i+1,j d i,j+1 d i+1,j+1 d i+1,j+1 = min{d i+1,j +1, d i,j+1 +1, d i,j +e i,+1,j+1 } e i,j ={0, S 1 i = S 2 j ; 1, S 1 i S 2 j }

12 Подмена задачи и обобщение Заменим расстояния d i,j на – d i,j. Тогда операцию min надо заменить на max. Прибавим к – d i,j ½ (w i,j = ½ – d i,j ), тогда получим функцию сходства: совпадение = ½, замена = –½, делеция = –1. Функцию сходства W легко обобщить, варьируя штрафы за замену и делеции. Новая задача: написать одну последовательность под другой так, чтобы максимизировать сходство Алгоритм Нидлмана – Вунша решает эту задачу, используя динамическое программирование.

13 Граничные условия w i,j w i+1,j w i,j+1 w i+1,j+1 w1,1w1,1 начало w1,2w1,2 d2,1d2,1 w n,m-1 w n,m w2,1w2,1 w n-1,m конец При таких граничных условиях начальные и концевые делеции штрафуются

14 Как не штрафовать за концевые делеции w i,j w1,1w1,1 начало w1,2w1,2 w2,1w2,1 w n,m-1 w n,m w3,1w3,1 w n-1,m конец w n,m-2 w n-2,m w1,3w1,3 0 0 В граф добавляются ребра веса 0, ведущие из начала во все граничные вершины (i=1 | j=1) и из граничных вершин (i=n | j=m) в конец

15 Алгортим Нидлмана – Вунша: оценка времени работы и необходимой памяти Алгоритм просматривает все вершины графа В каждой вершине делается 3 сравнения Количество необходимых операций (время работы алгоритма): T=O(nm). Говорят, что алгоритм выравнивания квадратичен по времени работы. Для запоминания весов и восстановления оптимального выравнивания надо в каждой вершине запомнить ее вес и направление перехода. Таким образом, алгоритм квадратичен по памяти.

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

17 Линейный по памяти алгоритм Миллера – Маерса Разбиваем одну из последовательностей на две равные части Для каждой точки x линии раздела находим веса оптимальных выравниваний из начала в x и из конца в x: W + (x), W – (x) Вес оптимального выравнивания, проходящего через точку x равен W(x)=W + (x) + W – (x) Вес оптимального выравнивания равен W = max x (W(x)) Таким образом, за время T=Cn 2 найдена одна точка, через которую проходит оптимальное выравнивание S1 S2 x

18 Алгоритм Миллера – Маерса Найденная точка x разбивает матрицу выравнивания на четыре квадранта, два из которых заведомо не содержат оптимального выравнивания Для двух квадрантов, содержащих оптимальный путь можно применить тот же прием, и запомнить точки x' и x". Просмотр оставшихся квадрантов требует времени T=Cn 2 /2 (почему?) Продолжая процедуру деления пополам найдем все точки, через которые проходит оптимальный путь. Время работы алгоритма T=Cn 2 +Cn 2 /2+Cn 2 /4+…= = Cn 2 (1+1/2+1/4+1/8+…); T=2Cn 2 S1 S2 x x' x" Важно, что при просмотре мы не запоминали обратных переходов! x x' x" x x' x" x x' x"

19 Еще один способ сэкономить время и память Ясно, что выравнивания D1 и D2 не представляют интереса, поскольку содержат в основном делеции Разумные выравнивания (A) лежат в полосе Алгоритм: задаемся шириной полосы w и просматриваем только те вершины графа, что лежат в указанной полосе. D2 D1 A w

20 Локальное выравнивание Локальным оптимальным выравниванием называется такое оптимальное выравнивание фрагментов последовательностей, при котором любое удлинение или укорочение фрагментов приводит только к уменьшению веса Локальному оптимальному выравниванию отвечает путь с наибольшим весом, независимо от того, где он начинается и где кончается Локальное оптимальное выравнивание может иметь бóльший биологический смысл, чем глобальное, но только если математическое ожидание веса сравнения букв, случайно взятых из последовательностей, отрицательно (почему?) Например, для алфавита из 4 букв, встречающихся с одинаковой частотой, годятся параметры 1 за совпадение, –1 за замену или 5 за совпадение, –2 за замену, но не имеет смысла использовать 5 за совпадение и –1 за замену.

21 Алгоритм Смита – Ватермана w i,j w1,1w1,1 начало w1,2w1,2 w2,1w2,1 w n,m-1 w n,m w3,1w3,1 w n-1,m конец w n,m-2 w n-2,m w1,3w1,3 0 0 В граф добавляются ребра веса 0, ведущие из начала во все вершины и из всех вершин в конец

22 Алгоритм Смита – Ватермана Пусть есть какой-то путь с неотрицательными весами Построим график веса вдоль пути Абсолютный максимум на этом графике определит точку окончания пути w max

23 Алгоритм Смита – Ватермана Точка конца пути (от нее начинаем обратный просмотр и восстановление пути) определяется так: (i max, j max ) = argmax (w i,j ) w i,j = max { w i-i,j-1 + e i,j, i > 1, j > 1 w i-1,j – d, i > 1 w i,j-1 – d, j > } Пусть (при одинаковых параметрах) мы получили вес глобального выравнивания S glob и вес локального выравнивания S loc. Какая величина больше?

24 Более общая зависимость штрафа за делецию от величины делеции Простейшая модель делеции: элементарное событие – удаление одного символа. Протяженная делеция – несколько независимых событий удаления одного символа. Работает плохо. Более реалистичная модель: делеция нескольких символов происходит за одно элементарное событие, а размер делеции является некоторой случайной величиной. Поэтому в качестве штрафа хорошо бы взять что-нибудь вроде Δ ( l ) = a log( l + 1 ), где l – длина делеции В любом случае функция Δ ( l ) должна быть выпуклой – должно выполняться неравенство треугольника: Δ ( l 1 + l 2 ) Δ ( l 1 ) + Δ ( l 2 )

25 Более общая зависимость штрафа за делецию от величины делеции. Алгоритм. wwwwwwwwwwwwwwww Теперь надо просматривать все возможные варианты делеций. Поэтому в каждую вершину входит не 3 ребра, а примерно (n+m)/2 ребер, где n, m – длины последовательностей Поэтому время работы алгоритма становится кубичным: T = O ( nm (n+m) )

26 Аффинные штрафы за делецию Δ ( l ) = d open + l d extВместо логарифмической зависимости используют зависимость вида: Δ ( l ) = d open + l d ext d opend open – штраф за открытие делеции d extd ext – штраф за длину делеции Δ l d open d ext a log( l + 1 )

27 Алгоритм для аффинных штрафов w v w v w v w v Веса на ребрах e i,j сопоставление d open d open открытие делеции d ext d ext продолжение делеции e i,j закрытие делеции Модификация стандартного графа: 1. 1.В каждой ячейке вводится дополнительная вершина (v), отвечающая делеционному пути 2. 2.Вводятся делеционные ребра для открытия и закрытия делеции (из вершин типа w в вершины типа v и обратно) 3. 3.Ребра, отвечающие продолжению делеции переносятся на новые вершины Число вершин графа равно 2mn число ребер равно 5mn Трудоемкость алгоритма равна: T = O (mn)

28 Рекурсия для аффинных штрафов w i, j = max ( w i-1, j-1 +e i j, v i-1, j-1 +e i j, 0 ); v i, j = max ( w i, j – d open, v i-1, j – d ext, v i,j-1 – d ext ); (i max, j max ) = argmax (w i,j )

29 Статистика выравниваний

30 Параметры выравнивания В простейшем случае есть три параметра: –премия за совпадение (match) –штраф за несовпадение (mism) –штраф за делецию (indel) match=1.Если все параметры умножить на одну и ту же положительную величину, то само оптимальное выравнивание не изменится, а вес выравнивания умножится на ту же величину. Поэтому можно положить match=1. Если mism > 2 indel, то оптимальное выравнивание не будет иметь замен (почему?)

31 Статистика выравниваний Допустим мы выровняли две последовательности длиной 100 и получили вес 20. Что это значит? Может быть при выравнивании двух случайных последовательностей будет тот же вес? А что такое случайные последовательности?

32 Модели случайных последовательностей Базовая (вообще говоря неправильная) модель – бернуллиевские последовательности: символы генерируются независимо друг от друга с заданной вероятностью. Для этой модели математика проще и проще получить оценки. Уточненная модель (лучше, но тоже неправильная) – марковская цепь (вероятность появления следующего символа зависит от нескольких предыдущих символов). Математика значительно сложнее. Почти ничего не известно.

33 Частные случаи локального выравнивания mism = 0, indel = 0 – максимальная общая подпоследовательность mism =, indel = – максимальное общее подслово

34 Наибольшая общая подпоследовательность Длина общей подпоследовательности есть случайная величина r(n), зависящая от длины последовательностей. Пусть две последовательности длиной n разбиты каждая на два фрагмента длиной n 1 и n 2 (n 1 + n 2 = n) Ясно, что общая подпоследовательность будет не короче, чем объединение общих подпоследовательностей для фрагментов: r(n) r 1 (n 1 )+r 2 (n 2 ) (попробуйте понять смысл неравенства) Отсюда следует, что математическое ожидание E(r(n)) E(r(n 1 )) + E(r(n 2 )), или E(r(n)) cn Можно показать, что E(r(n)) – ( E(r(n 1 )) + E(r(n 2 )) ) 0 :Поэтому: r(n) r 1 (n 1 ) r 2 (n 2 ) E(r(n)) cn (n E(r(n)) cn (n )

35 Наибольшее общее слово Наложим одну последовательность на другую. Будем идти вдоль пары последовательностей и, если буквы совпали, то будем считать успехом, иначе – неудача. Имеем классическую схему испытаний Бернулли. Наибольшему общему слову при таких испытаниях будет соответствовать максимальная серия успехов. Известно, что средняя величина максимальной серии успехов равна: E(l) = log 1/p (n) Возможных наложений много (порядка длины последовательности). Максимальное общее слово есть максимум от максимальных серий успехов при всех возможных наложениях. Показано (Waterman), что: E(l) log 1/p (nm) + log 1/p (1-p) + γlog 1/p (e) – ½ log 1/p (Knm), (m, n 0,577 E(l) log 1/p (nm) + log 1/p (1-p) + γlog 1/p (e) – ½ = log 1/p (Knm), (m, n, γ 0,577) (l) [ π log 1/p (e) ] 2 / 6 + ½, (l) [ π log 1/p (e) ] 2 / 6 + ½, (не зависит от n !)

36 Зависимость от параметров ожидаемогоПоказано, что зависимость ожидаемого веса выравнивания от длины последовательности может быть либо логарифмической, либо линейной в зависимости от параметров. Все пространство параметров разбивается некой поверхностью на две области поведения. mism indel W = O(log(n)) W = O(n) При безделеционном выравнивании поведение логарифмическое, если мат.ожидание вéса сравнения двух случайных сегментов отрицательно.

37 Матрицы замен

38 Откуда берутся параметры для выравнивания? Пусть у нас есть выравнивание. Если последовательности случайные и независимые (модель R), то вероятность увидеть букву α против β p(α, β | R) = p(α) p(β) а вероятность выравнивания (x,y) будет равна p(x,y | R) = Π p(x i ) Π p(y i ) Если выравнивание не случайно (модель M), то p(x,y | M) = Π p(x i, y i ) Отношение правдоподобия: p(x,y | M) Π p(x i, y i ) = p(x,y | R) Π p(x i ) Π p(y i ) Логарифмируя, получаем log( p(x,y|M)/p(x,y|R) ) = s(x i,y i ); Матрица замен: s(α, β) = log(p α β /p α p β )

39 Серия матриц BLOSUM База данных BLOCKS ( Henikoff & Henikoff ) – безделеционные фрагменты множественных выравниваний (выравнивания получены специальной программой). В каждом блоке отбираем подмножество последовательностей такое, что для каждой пары в нём процент идентичных аминокислот не больше заданного значения ID. В урезанном блоке в каждой колонке подсчитываем число пар аминокислот n bl col (α, β) Усредняем по всем колонкам и по всем блокам: f (α, β) = n bl col (α, β) / N col Элемент матрицы BLOSUM ID : BLOSUM ID (α, β) = log( f (α, β) / f (α) f ( β) )

40 Серия матриц PAM Point Accepted Mutation – эволюционное расстояние, при котором произошла одна замена на 100 остатков. Эволюционный процесс можно представить как марковский процесс. Если в начальный момент времени t = 0 в некоторой позиции был остаток α, то через время Δt в этой позиции с некоторой вероятностью будет остаток β: p(β| α, Δt) = M Δt (β, α) M Δt – эволюционная матрица Через время 2Δt p(β| α, 2Δt) = γ M Δt (β, γ) M Δt (γ, α) = M Δt 2 (β, α) Через время NΔt p(β| α, NΔt)= M Δt N (β, α)

41 Серия матриц PAM Находим выравнивания, отвечающие расстоянию PAM1 Находим частоты пар и вычисляем частоты пар. Поскольку расстояние мало, то один из символов соответствует предку: p(αβ) = p(α β) p(α)+ p(β α) p(β) полагая эволюцию равновесной: p(α β) p(α) = p(β α) p(β) получаем p(α β) = 2p(αβ) / p(α) p(α α) = 1 – βα p(α β) Марковский процесс: p N (αβ)= p N (α β) p(α) PAM N (αβ) = log (p N (α β) / p β )

42 Распределение экстремальных значений Пусть вес выравнивания x (случайная величина) имеет распределение G(S) = P(x < S) Тогда при N независимых испытаниях распределение максимального значения будет G N (x) = G N (x); NМожно показать, что для нормально распределенного G(x) при больших N G N (x) exp(–KN e –λx )

43 E-value и P-value Для бернуллиевских последовательностей длин m и n математическое ожидание количества независимых локальных выравниваний с весом >S описывается формулой ( Karlin &Altschul ) : E(S) = Kmn e –λS где λ – положительный корень уравнения αβ p α p β e λs(αβ) = 1 s(αβ) – матрица замен (с отрицательным матожиданием сравнения случайных букв: αβ p α p β s(αβ) < 0 ) K – константа, зависящая от p α и s(αβ) Для поиска по банку n – суммарная длина банка. E-value: E(S)E-value: E(S) ожидаемое количество выравниваний с таким или большим весом P-value: p(x > S) = 1– e –E(S)P-value: p(x > S) = 1– e –E(S) вероятность встретить (хоть одно) выравнивание с таким или большим весом

44 Проблема малой сложности В формулах для E-value и P-value присутствуют вероятности букв p α Эти формулы работают удовлетворительно, если частоты букв в последовательностях примерно равны вероятностям букв в модели Но в биологических последовательностях встречаются участки, в которых частоты букв сильно сдвинуты (например, богатые пролином петли в белках). Выравнивание таких участков при использовании стандартной матрицы s(αβ) получит высокий вес даже при отсутствии реальной гомологии. Близкая проблема: участки вида … atgcatgcatgcatgcatgcatgcatgc …

45 Проблема малой сложности: подходы к решению «Маскировка» участков малой сложности: –применялась в BLAST до 2005 года как единственный вариант; сейчас применяется в основном для нуклеотидных последовательностей; –для белковых последовательностей применяется программа seg, маскирующая участки с частотами букв, сильно отличающимися от «базовых»; –для нуклеотидных последовательностей применяется программа dust, которая маскирует участки с сильно «сдвинутым» составом триплетов (подслов длины 3); –основной недостаток «всё или ничего»: можно замаскировать биологически осмысленное выравнивание, с другой стороны участок чуть выше «порога сложности» может дать много бессмысленных выравниваний с последовательностями банка. Корректировка матрицы замен s(αβ) в соответствии с частотами букв в последовательностях: –предложена в статьях Yu (Юй) и Альтшуля 2003–2005; –с 2005 включена по умолчанию в BLAST для белков.

46 Алгоритм dust На входе: нуклеотидная последовательность (в алфавите A, T, G, C). На выходе: маскированная последовательности, отличающаяся от исходной тем, что участки малой сложности заменены буквами N. Для каждой подпоследовательности a и для каждого триплета (слова длины 3) t определяется число c t (a), равное числу подслов в a, совпадающих с t. Вес S(a) определяется как (n–3) –1 t c t (a)(c t (a) – 1)/2, где n – длина a. При более или менее равных частотах триплетов вес мал, а если какие-то триплеты встречаются сильно чаще других, велик (почему?). Алгоритм: – рассматриваются все подпоследовательности a длины W=64 – рассматриваются все подпоследовательности a длины W=64 – для каждого a находится префикс a' максимального веса S(a') – для каждого a находится префикс a' максимального веса S(a') – если S(a') > T = 2, то в a' находится суффикс a'' максимального веса S(a'') – если S(a') > T = 2, то в a' находится суффикс a'' максимального веса S(a'') – все найденные так подпоследовательности a'' маскируются. – все найденные так подпоследовательности a'' маскируются.

47 Алгоритм seg Вход: белковая или нуклеотидная последовательность. Выход: маскированная последовательность (нуклеотиды в участках малой сложности заменяются на N, аминокислоты на X). В алгоритме SEG используется две меры сложности последовательности: K 1 = 1/L log N Ω здесь N – размер алфавита (4 или 20), L – длина последовательности, а Ω равно числу различных последовательностей, имеющих в точности те же частоты букв, что и данная. Чем больше частоты букв различаются между собой, тем меньше K 1 (почему?) K 1 = 1/L log N Ω здесь N – размер алфавита (4 или 20), L – длина последовательности, а Ω равно числу различных последовательностей, имеющих в точности те же частоты букв, что и данная. Чем больше частоты букв различаются между собой, тем меньше K 1 (почему?) K 2 = α f (α) log 2 f (α), где f – частоты букв. K 2 представляет собой приближение K 1 : K 2 K 1 при L K 2 = α f (α) log 2 f (α), где f – частоты букв. K 2 представляет собой приближение K 1 : K 2 K 1 при L

48 Алгоритм seg (окончание) Проходим по последовательности окном длины W (по умолчанию W = 12) и определяем сложность (K 2 ) для каждого окна. Проходим по последовательности окном длины W (по умолчанию W = 12) и определяем сложность (K 2 ) для каждого окна. Помечаем все буквы, попавшие в окна, чья сложность ниже первого порога, по умолчанию равного 2,2. Помечаем все буквы, попавшие в окна, чья сложность ниже первого порога, по умолчанию равного 2,2. Помечаем все буквы, попавшие в окна, пересекающиеся с уже помеченными участками и чья сложность ниже второго порога, по умолчанию равного 2,5. Помечаем все буквы, попавшие в окна, пересекающиеся с уже помеченными участками и чья сложность ниже второго порога, по умолчанию равного 2,5. Для каждого непрерывного отрезка из помеченных букв уточняем границы участка малой сложности, используя меру K 1. Для каждого непрерывного отрезка из помеченных букв уточняем границы участка малой сложности, используя меру K 1.

49 Корректировка матрицы замен Матрица замен (BLOSUM или PAM) может быть представлена в виде: s(αβ) = λ –1 ·log 2 ( f (α, β) / f (α) f ( β) ), где f (α, β) – частоты замен в эталонных выравниваниях, а f (α) – частоты букв f (α) = p α, вероятностям бернуллиевской модели в формуле для E-value Имеется очевидное равенство: α f (α, β) = f (β) для любой буквы β Пусть теперь имеется пара последовательностей, в которых частоты букв f ' (α) f (α). Задача: подобрать f '(α, β) так, чтобы: а) α f ' (α, β) = f ' (β), т.е. частоты замен соответствовали частотам букв; б) матрица f ' (α, β) была (в каком-нибудь смысле) «самой близкой» к матрице f (α, β) изо всех матриц, удовлетворяющих условию (а). б) матрица f ' (α, β) была (в каком-нибудь смысле) «самой близкой» к матрице f (α, β) изо всех матриц, удовлетворяющих условию (а). Когда f ' (α, β) подобраны, для выравнивания используется матрица замен s' (αβ) = λ –1 ·log 2 ( f '(α, β) / f '(α) f '( β) )

50 Какая матрица самая близкая? Используется один из двух вариантов. Вариант 1. Самая близкая матрица f '(α, β) та, для которой минимальна величина D = αβ f ' (α, β) log( f ' (α, β) / f (α, β) ) Вариант 2. Самая близкая матрица та, для которой, во первых, выполнено равенство: αβ f ' (α, β) log( f ' (α, β) / f '(α) f '(β) ) = αβ f (α, β) log( f (α, β) / f (α) f (β) ) а во-вторых, величина D минимальна среди всех матриц, для которых выполнено это равенство. Именно второй вариант реализован в программе BLASTP. Смысл равенства в том, что величина в правой и левой его части («энтропия матрицы») характеризует отличие частот пар в выравниваниях от произведений частот букв. Для матриц PAM энтропия уменьшается с ростом номера (для PAM10 она больше, чем для PAM100), для матриц BLOSUM наоборот (почему?)

51 Поиск по банку

52 Поиск по банку. Постановка задачи На входе: последовательность-запрос (query) и банк из большого количества последовательностей. Нужно найти в банке кандидатов в гомологи запроса. Кандидатами в гомологи считаются последовательности с достаточно высоким весом локального выравнивания с запросом. Достаточно ли высок вес, оценивается по P-value или E- value. Главная проблема – время: алгоритм Смита – Ватермана слишком медленный. Нужны приёмы быстрого выравнивания. Нахождение всех достаточно хороших выравниваний не гарантируется (эвристические алгоритмы).

53 Поиск по банку. Хеширование. l-tuple l-граммаПодготовка банка – построение хэш-таблицы. Хэш- функция – номер слова заданного размера (l-tuple, l-грамма). В хэш-таблице хранятся списки ссылок на последовательности и на позиции в последовательностях, где встречается соответствующая l-грамма. При поиске запроса (query) в последовательности запроса последовательно находятся l-граммы, далее, по хэш-таблице для них находятся соответствующие документы и позиции. затравкой, якорем, seed.Пара совпадающих l-грамм в запросе и в банке называется затравкой, якорем, seed.

54 Поиск по банку. BLAST1. Ищем якоря с помощью хэш-таблицы. Длина якоря: l = 3 или 4 для белков, l 7 для н.к. Каждый якорь расширяем с тем, чтобы получить сегмент совпадения наибольшего веса (HSP – high scoring pair). Оцениваем его статистическую значимость (E-value), и, если она больше порога, то выдаём (Altschul, Gish, Miller, Myers, Lipman, 1990)

55 Поиск по банку. BLAST2. Расширяются не одиночные якоря, а пары якорей, которые находятся недалеко и на близких диагоналях При расширении пары якорей допускаются делеции T-соседней l-граммой L T для l-граммы L называется такая l- грамма, что вес ее сравнения с L не меньше заданного T: s(L i, L i T ) T Для аминокислотных последовательностей при просмотре запроса формируем не только те l-граммы, которые встретились в нем, но также все T-соседние l-граммы. Используются l = 2 или 3 для белков, l 4 для н.к. (Altschul et al., 1997)

56 Поиск по банку. FASTA. Используются якоря длины l = 1 или 2 (для белков); l = 3, …, 6 (для нк). Два якоря (i 1,j 1 ), (i 2,j 2 ) принадлежат одной диагонали, если i 1 – j 1 = i 2 – j 2 Мощностью диагонали называется количество якорей, принадлежащих диагонали. Иногда в мощность диагонали включают мощности соседних диагоналей (чтобы учесть возможность небольших делеций) Отбираем n (например n = 10) самых мощных диагоналей и для каждой строим локальное выравнивание в полосе заданной ширины вокруг диагонали (Wilbur, Lipman, Pearson)

57 Ещё один алгоритм быстрого выравнивания Ищем якоря Якорь (i 1,j 1 ) предшествует якорю (i 2,j 2 ), если i 1 < i 2 & j 1 < j 2 & i 2 – i 1 < d & j 2 – j 1 < d Получаем ориентированный граф с небольшим количеством вершин и ребер Можно найти оптимальную цепочку якорей методом динамического программирования

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

59 δ-функция Определение: δ (x) = 0, x 0; + δ (x) dx = 1 δ (x) dx = 1- Свойство:+ δ (x – a) f(x) dx = f(a) δ (x – a) f(x) dx = f(a)- Символ Кронекера. Определение: 0, m n 0, m n δ(m, n) = δ m, n = δ(m, n) = δ m, n = 1, m = n 1, m = n

60 Γ-функция Определение: + Γ(x) = e –t t x–1 dt Γ(x) = e –t t x–1 dt 0 0 Свойства: Γ(0) =1; Γ(1) =1; Γ(x+1) =x Γ(x) Следствие: Γ(n+1) = n! Формула Стирлинга Γ(x+1) 2π x x x e – x, x +

61 Модели последовательностей Подсчитаем частоту встречаемости букв в русском языке. Будем генерировать символы с подсчитанными частотами Получим ли мы что-либо похожее на слова русского языка? Способ генерации случайных слов (последовательностей) называется (статистической) моделью языка.

62 Марковские цепи Очевидно, что ничего хорошего не получится. Есть наблюдение, что вероятность появления буквы зависит от предыдущей буквы. Например, 'ь' никогда не появляется после гласной буквы; две гласные подряд встречаются редко, и т.п. НЕ НЕЗАВИСИМЫСимволы в последовательности НЕ НЕЗАВИСИМЫ. Марковская модель первого порядка: p(S i ) = β p(β S i ) δ (β, S i-1 ) p(β α) – матрица переходных вероятностей: p(β α) = P(α | β) вероятность появления символа α при условии, что предыдущий символ – β.

63 Матрица переходных вероятностей Размер матрицы – Σ × Σ, где Σ – число исходов (например, размер алфавита). Вообще говоря, количество исходов может быть бесконечным. Сумма значений в строке : β p(β α) = 1 β p(β α) = 1 Пусть первый символ был a. Тогда распределение вероятностей для второго символа будет p 2 (α) = p (a α) = β p(β α) δ(β,a); p 2 = Pδ(β,a) Здесь P означает матрицу переходных вероятностей, δ(β,a) – одномерный вектор (столбец) из нулей и единицы. Распределение вероятностей для n-го символа в генерированной последовательности: p n (α) = β p(β α) p n–1 (β); p n = P p n–1 = P 2 p n–2 =…= P n–2 p 2 = P n–1 δ(β,a)

64 Марковские цепи и эволюция Пусть происходит эволюция некоторого белка. Ясно, что некоторые замены (например, IV) фиксируются часто, а некоторые – редко. Если в некоторой позиции в предке была аминокислота G, то можно построить распределение частот аминокислоты через некоторое время. Марковский процесс: изменение аминокислотного остатка в данной позиции.

65 Марковские цепи и эволюция Матрица P переходных вероятностей имеет вид: AGVI… A 1-S A dt P GA dt P VA dt P IA dt … G P AG dt 1-S G dt P VG dt P IG dt … V P AV dt P GV dt 1-S V dt P IV dt … I P AI dt P GV dt P VI dt 1-S I dt … ……………… S α = βα P(β α) Распределение вероятностей через время t: p(α, t) = exp (t Q) p(β,t) Загадочный объект – экспонента от матрицы. Это матрица, которая определяется через ряд Тейлора, только вместо степени числа пишется соответствующая степень матрицы. P = = I + Q dt

66 Марковские цепи высших порядков Вероятность появления очередного символа зависит не от одного, а от нескольких предыдущих символов: p(S i ) = α1 α2 … αm δ(S i-k,α 1 ) δ(S i-k+1,α 2 )… δ(S i-1,α m ) P(S i | α 1 α 2 … α m )

67 Оценка порядка марковской цепи в модели последовательностей Оценка переходной вероятности: p*(α 1,…, α k ; α k+1 ) = n(α 1,…, α k, α k+1 ) / n(α 1,…, α k,*); n(α 1,…, α k,*) = β n(α 1,…, α k, β) Информационный критерий Байеса: лог- правдоподобие: L k = все слова n(α 1,…, α k, α k+1 ) ln p*(α 1,…, α k ; α k+1 ) BIC( k ) = –2L k + a ln (N) a = A k (A–1) – число независимых параметров цепи (A – размер алфавита), N – число последовательностей в обучении. истинный порядок цепи k * = argmin k (BIC( k ) )

68 Задача Испытания лекарства: М: 50/90=5/9 > 4/12=3/9 Ж: 10/12=5/6 > 80/120=4/6 Вместе: Л: 60/102=20/34 К: 84/132=7/11=21/33 20/34 < 21/33 !!!! МужЖен лекКнтрлеккнтр К чему бы это?

69 Введение в байесову статистику Задача. Мы 3 раза бросили монету и 3 раза выпал орел. Какова вероятность выпадения орла у этой монеты? –Если мы уверены, что монета не кривая, то p = ½ –Допустим, что мы взяли монету из мешка, а в мешке монеты разной кривизны. Но при этом мы знаем как распределена кривизна монет P a (p) (априорное распределение). распределений вероятностей данной –Мы хотим на основе наблюдения 3о и априорного распределения распределений вероятностей оценить вероятность выпадения орла у данной монеты.

70 Введение в байесову статистику P(3o | p) = p 3 ; P(3o, p) = P(3o | p) P a (p) = P(p | 3o) P(3o); P(p | 3o)= {P(3o | p) P a (p)} / P(3o); Загадочный объект P(3o) – безусловная вероятность трех орлов. Определяется из условия нормировки: P(p | 3o) = 1; Окончательно, распределение вероятностей вероятности орла будет: P(p | 3o)= p 3 P a (p) / p 3 P a (p) ;

71 Введение в байесову статистику P(p | 3o)= p 3 P a (p) / p 3 P a (p) dp; В качестве оценки для искомой вероятности удобно иметь число, а не распределение: –Максимальное значение p ML =argmax p ( P(3o | p)) – максимальное правдоподобие (max likelihood, ML) –Среднее значение p E =E( P( p | 3o))= p P( p | 3o) dp;

72 Введение в байесову статистику ML оценка (максимальное правдоподобие): p ML = argmax (p 3 ) = 1; E оценка (матожидание апостериорной вероятности) p E = p 4 P a (p) dp / p 3 P a (p) dp; –Если мы уверены, что монета правильная, то P a (p)=δ(p – ½); p E = ½ ; –Если мы ничего не знаем о распределении P a (p), то положим P a (p) = const. Тогда p E = p 4 P a (p) dp / p 3 P a (p) dp = (1/5) / (1/4) = 4/5 ; p E (no) = (n+1)/(n+2); В более общем случае p E (no) = (n+1)/(n+2); MAP оценка (максимум апостериорной вероятности) p MAP = argmax { P(p | 3o)};

73 Определения Пусть у нас есть несколько источников Y событий X (например, несколько монет). Тогда : условная вероятность P(X | Y) – условная вероятность совместная вероятность P(X,Y) = P(X | Y) P(Y) – совместная вероятность полная вероятность P(X) = Y P(X,Y) = Y P(X |Y) P(Y) – полная вероятность – апостериорная вероятность выбора источника (правдоподобие гипотезы) P(Y | X) – апостериорная вероятность выбора источника (правдоподобие гипотезы) априорная вероятность выбора источника P(Y) – априорная вероятность выбора источника Теорема Байеса: P(X | Y)= P(Y | X) P(X) / P(Y)

74 Пример 1 Пусть есть две кости – правильная и кривая (с вероятностью выпадения шестёрки, равной ½ вместо 1/6). Пусть нам подсовывают кривую кость с вероятностью 1%. Мы бросили кость 3 раза и 3 раза получили 6. Какова вероятность того, что нам дали кривую кость? PP( кривая кость | 3 шестерки ) = PP P( 3 шестерки | кривая кость ) P( кривая кость ) P P( 3 шестерки ) PPP P( 3 шестерки )=P( 3 шестерки | кривая кость ) P( кривая кость ) + PP P( 3 шестерки | правильная кость ) P( правильная кость ) = (1/6) = = P P( кривая кость | 3 шестерки ) = / =0.21 Вывод – кость скорее правильная!Вывод – кость скорее правильная! Сколько шестерок подряд надо, чтобы мы поняли, что нас обманывают?Сколько шестерок подряд надо, чтобы мы поняли, что нас обманывают?

75 Пример 2 Есть редкая болезнь, P(б.)=10 -6 Имеется тест со свойствами: если больны, то вероятность ошибки теста P(–|б.) = 0, если здоровы, P(+|з.)=10 -4 Стоит ли проходить тест? P(б.|+)=P(+|б.)·P(б.) / P(+); P(+)= P(+|б.)·P(б.) + P(+|з.)P(з.)10 -4 P(б.|+)=10 -6 /10 -4 =10 -2 P(з.|+)=1 – P(б.|+)=0,99

76 Пример 3 В последовательности A нашли взаимно-комплементарную структуру. Последовательность B имеет степень сходства id. В выроненных участках последовательности B нашли аналогичную шпильку (буквы в ней не обязательно такие же, важно, что шпилька) Какова значимость этого наблюдения? α1α1α1α1 α2α2α2α2 β1β1β1β1 β2β2β2β2 Символы α 1 и взаимно-комплементарны; Символы β 1 и β 2 также взаимно-комплементарны; Найдем вероятность такого события.

77 Пример 3 (продолжение)

78 Пример 4 Пусть ORF начинается всегда с ATG и кончается стоп-кодоном. Найти распределение длин ORF. P(ORF длины L | c поз. i)=P(start)·P L-2 (codon)P(stop) P(ORF длины L)= P(ORF длины L | c поз. i)/P(ORF c поз. i)

79 Оценка параметров по результатам Пусть у нас есть наблюдение D и некоторый набор параметров распределения θ, которые мы хотим оценить (см. пример про 3 орла). Кроме того, у нас есть представление о том, как эти параметры распределены (prior) Апостериорное распределение вероятностей параметров получаем из теоремы Байеса: PP P P(θ) P(D |θ) P(θ | D) = PP θ P(θ') P(D |θ')

80 Распределение Дирихле Определение: D(θ|α)=Z -1 θ i αi δ( θ i – 1); –Z – нормировочный множитель –α i – параметры распределения – θ i 0 – область определения распределения –δ – дельта-функция (δ(x)=0, x0; δ(x)dx=1;) θ1θ1 θ2θ2 θ3θ3 Симплекс Задача: найти объем симплекса в n-мерном пространстве

81 Оценка по максимуму апостериорной вероятности (MAP) Пусть есть модель с L исходами. Пусть есть наблюдения n 1,n 2,…,n L. Пусть априорное распределение – распределение Дирихле с параметрами α 1, α 2,…, α L : Найдем максимальную апостериорную вероятность Условие максимума при ограничении θ=1

82 MAP-оценка

83 prior = распределение Дирихле псевдо-отсчетами (pseudo counts).Часто в качестве prior используют распределение Дирихле. Параметры этого распределения α i называют псевдо-отсчетами (pseudo counts). Они определяют степень нашего доверия к результатам На графиках показаны распределения для случая 4-х орлов при 4-х бросаниях монеты. θ – вероятность орла P(D | θ) –Синяя линия – P(D | θ) P(θ) –Красная линия – распределение Дирихле P(θ) P(θ | D) –Желтая линия – апостериорная вероятность выпадения орла P(θ | D) α 1 =1, α 2 =1 α 1 =3, α 2 =3

84 Скрытые Марковские модели (HMM)

85 Пример Пусть некто имеет две монеты – правильную и кривую. Он бросает монету и сообщает нам серию результатов. С некоторой вероятностью он может подменить монету. Моменты подмены монеты нам неизвестны, но известно: –результаты бросков –вероятность с которой он заменяет монету –степень кривизны каждой монеты Задача: определить моменты смены монеты

86 Биологические примеры Дана аминокислотная последовательность трансмембранного белка. Известно, что частоты встречаемости аминокислот в трансмембранных и в растворимых частях белка различаются (аналог разных монет). Определить по последовательности где находятся трансмембранные участки. Дана геномная последовательность. Статистические свойства кодирующих областей отличаются от свойств некодирующих областей. Найти кодирующие области.

87 Описание HMM конечного автоматаПример с монетой можно представить в виде схемы конечного автомата: состояния –Прямоугольники означают состояния эмиссии –Кружки означают результат бросания (эмиссии) переходы –Стрелки – возможные переходы между состояниями –Числа около кружков – вероятности эмиссии e i –числа около стрелок – вероятности переходов между состояниями a ik –Есть начальное и конечное состояния Сумма весов исходящих стрелок равна 1 Сумма весов эмиссии в каждом состоянии рана 1 (Конечный автомат Мура – в алгоритмах был автомат Мили) B E 0.1

88 Решение задачи о монете Пусть нам известна серия бросков: Этой серии можно поставить в соответствие граф переходов: –Красные вершины соответствуют эмиссии соответствующих значений правильной монетой –Синие вершины – эмиссия значений кривой монетой –на ребрах – вероятности переходов –на вершинах – вероятности эмиссии Каждому пути по графу соответствует одна из гипотез о порядке смены монеты BE a -- e 0- a +- a -+ e 1- e 1+ e 0+

89 Решение задачи о монете Для любого пути можно подсчитать вероятность того, что наблюденная серия соответствует этому пути (порядку смены монет) P = a 0,1 a i,i+1 e i+1 Найдем путь, отвечающий максимуму P. log является монотонной функцией, поэтому можно прологарифмировать формулу для вероятности. (почему?) π*= argmin {– log a 01 – π (log(a i,i+1 ) + log(e i+1 )} π*= argmin {– log a 01 – π (log(a i,i+1 ) + log(e i+1 )} Это задача поиска оптимального пути на графе. Решается динамическим программированием ViterbiАлгоритм динамического программирования для поиска наиболее вероятного пути называется Viterbi

90 Viterbi рекурсия ОбозначенияОбозначения –v k (i) –v k (i) – наилучшая вероятность пути, проходящего через позицию i в состоянии k. π k (i) – π k (i) – наилучший переход из позиции i в состоянии k в предыдущую позицию (предыдущее состояние) π*(i) – π*(i) – наилучшее состояние в позиции i ИнициацияИнициация v k (0) = δ(0,k); v k (0) = δ(0,k); k – номер состояния РекурсияРекурсия v k (i) = e k (x i ) max m ( v m ( i – 1 ) a mk ); π(i,k) = argmax m ( v m ( i – 1 ) a mk ); π(i,k) = argmax m ( v m ( i – 1 ) a mk ); обратный переход ЗавершениеЗавершение P(x,π*)= max m ( v m ( L ) a m0 ); π*(L) = argmax m ( v m ( L ) a m0 ); Оптимальный путьОптимальный путь π*( i – 1 ) = π ( i, π* ( i ) ) ;

91 Другая постановка задачи Для каждого наблюденного значения определить вероятность того, что в этот момент монета была правильной. i +Для этого надо просуммировать по всем путям, проходящим через точку i + вероятности этих путей. Для решения этой задачи достаточно вспомнить динамическое программирование над полукольцом с использованием операции сложения и умножения. Нас интересует вероятность P(π i =k |x) = P(x, π i =k) / P(x) P(x, π i =k) = P(x 1 …x i, π i =k) P(x i+1 …x L | π i =k)Оцениваем значение P(x, π i =k) = P(x 1 …x i, π i =k) P(x i+1 …x L | π i =k) f k (i) –Первый сомножитель f k (i) = P(x 1 …x i, π i =k) определяем просмотром вперед b k (i+1) –Второй сомножитель b k (i+1) = P(x i+1 …x L | π i =k) определяем просмотром назад

92 Алгоритм Forward / backward Forward: по определению f k (i) = P(x 1 …x i, π i =k) f 0 (0)=1, f k (0)=0, k>0 f 0 (0)=1, f k (0)=0, k>0 f l (i) = e i (x i ) k f k (i-1) a kl P(X)= k f k (L)a k0 Backward: b k (i) = P(x i+1 …x L | π i =k) b k (L) = a k0 b k (i) = l a kl e l (x i+1 ) b l (i+1) P(X)= l a 0l e l (x 1 ) b l (1)

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

94 Оценка параметров HMM при наличии обучающей выборки Здесь используется техника оценки параметров методом наибольшего правдоподобия. Пусть –x n – набор независимых наблюдений –θ – набор параметров, которые надо оценить θ * =argmax θ l(x 1 … x n | θ) = argmax θ { j log P(x j | θ)}Тогда надо максимизировать θ * =argmax θ l(x 1 … x n | θ) = argmax θ { j log P(x j | θ)}

95 Оценка параметров HMM при наличии обучающей выборки При условиях Метод неопределенных множителей Лагранжа aab aa

96 Оценка параметров HMM при наличии обучающей выборки Можно показать, что при большом количестве наблюдений справедливы оценки a kl = A kl / l' A kl' ; e k (b) = E k (b) / b' E k (b'); –A kl –A kl – наблюденное количество переходов между моделями –E k (b) –E k (b) – количество порожденных символов в соответствующих моделях При малых размерах выборки используют технику псовдоотсчетов, добавляя к наблюденным значениям некоторое количество шума.

97 Если нет обучающей выборки Итеративный алгоритм обучения Витерби.Итеративный алгоритм обучения Витерби. 1.Выберем некоторые наборы параметров HMM (обычно они генерируются случайно). 2.Найдем для них оптимальные пути во всех представленных примерах 3.По найденным оптимальным путям определим новые параметры, подсчитывая частоты эмиссии и переходов. 4.Перейдем к шагу 2. Итеративный алгоритм Баума-ВелчаИтеративный алгоритм Баума-Велча – то же самое, но параметры оцениваются с помощью Forward-Backward. –Показано, что алгоритм сходится (отношение правдоподобия растет на каждой итерации) Есть опасность нахождения локального, а не глобального экстремума.

98 Оценки параметров по Бауму – Велчу Имея заданные параметры модели можно определить вероятность перехода между состояниями: где f k (i) = P(x 1 …x i, π i =k), b l (i+1) P(x i+1 …x L | π i+1 =l) – значения, полученные при прямом и обратном проходе. Тогда для переходных и эмиссионных вероятностей получим оценки для количества переходов и порожденных символов: где x j – j-последовательность в выборке, f j k, b j l – результаты прямого и обратного прохода по последовательности x j

99 Предсказание кодирующих областей в прокариотах Реальная схема HMM для поиска кодирующих областей сложнее: –Включает в себя SD сайт –Учитывает неравномерность следования кодонов A C G T eAeA eCeC eGeG eTeT AT 11 AAC AAG AAT AAA 1 1 Кодоны p codon G T GA 1 1 T AA T AG Стоп p start p stop 1 1-p start Старт некодирующая последовательность

100 Оценка качества обучения Выборку разбивают на два подмножества – обучающую и тестирующую На первой выборке подбирают параметры На второй – тестируют и определяют качество обучения: –TP –TP – количество правильно определенных позитивных позиций (например, кодирующих) –TN –TN – количество правильно определенных негативных позиций (например, некодирующих) –FP –FP – количество неправильно определенных позитивных позиций (некодирующих, предсказанных как кодирующие) –FN –FN – количество неправильно определенных негативных позиций (кодирующих некодирующих, предсказанных как некодирующие)

101 Оценка качества обучения Специфичность: Sp = TP / (TP + FP) Чувствительность: Sen =TP / (TP + FN) Качество (пересечение/объединение) QQ =TP/(TP+FP+FN) Коэффициент корреляции CC=(TP*TN–FP*FN) / ((TP+FP)*(TN+FN)*(TP+FN)*(TN+FP)), реальность предсказание TP TN FNFP TN

102 Казалось бы … Построим модель с миллионом параметров, включая учет притяжения Луны. Можно ожидать, что в этом случае мы получим очень точную модель, которая будет правильно все предсказывать. НО… точность модели точность оценки параметровНО… При этом для оценки каждого параметра будет использовано примерно одно наблюдение. Поэтому хоть точность модели и велика, точность оценки параметров очень мала, и ее предсказательная сила будет также очень мала.

103 HMM и парное выравнивание

104 Конечный автомат для парного выравнивания M IXIXIXIX IYIYIYIY Порождение пары символов Символ в Y и делеция в X

105 HMM для выравнивания Парная HMMПарная HMM Состояния:Состояния: –Начало –Сопоставление (генерация пары сопоставленных символов) e ij =s(x i,y j ) –Генерация символа в X и делеция в Y e i =q(x i ) –Генерация символа в Y и делеция в X e i =q(y i ) –Конец M IXIXIXIX IYIYIYIY Begin End τ τ τ εε δ δ 1-ε-τ 1-2δ-τ δ δ τ 1-ε-τ 1-2δ-τ

106 Viterbi для выравнивания M IXIXIXIX IYIYIYIY Begin End τ τ τ εε δ δ

107 Случайная модель: независимое порождение последовательностей X Y Begin End 1-η Отношение правдоподобия: Вероятность для случайного независимого порождения последовательностей η η η η 1-η

108 Viterbi для log отношения правдоподобия Завершение:

109 Если есть несколько слабых выравниваний Можно оценить полную вероятность Для этого можно использовать Forward - алгоритм вычисления полной вероятности Доменная перестройка

110 ForwardM IXIXIXIX IYIYIYIY Begin End τ τ τ εε δ δ Инициация: Рекурсия: Завершение (полная вероятность):

111 Вероятностная генерация выравниваний На обратном пути мы выбираем переходы не по максимуму, а с вероятностями:

112 Вероятность того, что x i и y j выравнены y j и x i - выравнены

113 BackwardИнициализация Рекурсия Искомая вероятность

114 Информация и энтропия

115 Микро- и макросостояния (кое-что из статистической физики) Пусть у нас есть p состояний. Числом заполнения n i состояния i называется число частиц, находящихся в состоянии i. Микросостоянием системы из N частиц называется распределение (размещение) частиц по состояниям. Макросостоянием называется набор чисел заполнения. Одному Макросостоянию отвечает набор микросостояний Энтропией Макросостояния называется логарифм количества микросостояний, отвечающих данному макросостянию.

116 Энтропия По определению: S(N) = log( N! / (n 1 ! n 2 ! … n p !)); используем приближение n! = n n e – n. S(N) = N log N – n 1 log n 1 – n 2 log n 2 - … n p log n p + (- N + n 1 + n 2 + … +n p ) = (n 1 + n 2 + … +n p ) log N – n 1 log n 1 – n 2 log n 2 - … n p log n p ; окончательно получаем: S(N) = – N i f i log f i ;

117 Энтропия и информация Для источника символов энтропия равна: H источника = - α P(α) log 2 P(α) если P(x) = 0, то вклад этого члена равен 0. P(α) – вероятность генерации символа Энтропия – степень неопределенности при генерации символов Энтропия аддтивна: энтропия неопределенной последовательности X равна сумме энтропий позиций: H(X) = i H i = N H i H(X) = i H i = N H i Энтропия максимальна, если все символы равновероятны При генерации последовательности неопределенность становится определенностью. Полное Информационное содержание – потеря энтропии: I(X) = H before – H after = – i α P(α) log 2 P(α)

118 Информация Информация при генерации очередного символа: I = α P(α) log 2 P(α)= α P(α) I(α) I(α) –I(α) – частная информация Частная информация (информационное содержание) последовательности: I(X) = i I(x i ) = – i log 2 P(x i )

119 Информация выравнивания (bit-score) S 1 AFGILVQRSTASGNMFLC A|G| Q||TA|GN F|C S 2 AYGVLVQKTTATGNWYIC Информационное содержание выравнивания bit-score = – log 2 p(s 1 i, s 2 i );

120 Взаимная энтропия Вероятность макросостояния: Взаимная энтропия:

121 Взаимная информация Для двух распределений взаимная информация (расстояние Кульбака): Свойство: если f i p i, то I(f | p) > 0. Простое доказательство:

122 Профили

123 Способы описания множественного выравнивания Дано:Дано: множественное выравнивание. Задача:Задача: определить принадлежит ли некая последовательность данному семейству. Простейший способ описания множественного выравнивания – консенсус – все просто и ясно – пишется наиболее часто встречающаяся буква L[ST]XX…Регулярное выражение (используется в Pro-Site): L[ST]XX… Матрица частот встречаемости аминокислот в колонке LSPADKTNVKAAWGKV LTPEEKSAVTALWGKV LSEGEWQLVLHVWAKV LSADQISTVQASFDKV LSAAEKTKIRSAWAPV LTESQAALVKSSWEEF LSAAQRQVIAATWKDI Ls......v.a.W.kv L S T P E A G D Q

124 Энтропия колонки Пусть колонка содержит n α букв типа α. Тогда вероятность появления такой колонки при случайных независимых последовательностях будет определяться мультиномиальным распределением: N! P column = α p α nα ; p α – вероятность появления α α n α ! Логарифм этой величины равен: log ( P column ) = log N! + α (n α log p α – log n α !) Заменим n на N f α (f α – частота) и применим оценку для факториала n! (n/e) n. Получим полную энтропию колонки доказать! H column = log( P column ) = N α f α (log p α – log f α ); доказать! Величина I = – α f α (log p α – log f α ) информационным содержанием называется информационным содержанием колонки

125 HMM профиль Модель: каждая последовательность множественного выравнивания является серией скрытой Марковской модели. Профиль – описание Марковской модели. Каждой позиции соответствует свое состояние. Вероятности переходов между соседними состояниями равны 1. Вероятность того, что некоторая последовательность x соответствует профилю M: P( x | M)= e i (x i ); Значимость определяется отношением правдоподобия: сравнением с P( x | R) – вероятностью, что последовательность сгенерирована случайной моделью R: S = log (P( x | M) / P( x | R)) = log {e i (x i ) / q (x i )}; Величины w i (α)= log {e i (α) / q (α)} называют позиционной весовой матрицей (PSSM, PWM) B eAecef…eAecef… eAecef…eAecef… eAecef…eAecef… eAecef…eAecef… E M1M1 M2M2 M3M3 MnMn

126 HMM с учетом возможности вставок Делеция в профиле и в последовательности могут идти подряд (в отличие от парного выравнивания) Делеционные состояния – молчащие (не имеют эмиссии) Вероятность перехода в делеционное состояние зависит от позиции B E M1M1 M2M2 M3M3 MnMn D1D1 D2D2 D3D3 DnDn I0I0 I1I1 I2I2 I3I3 InIn Делеция в профиле Делеция в последовательности

127 Определение параметров модели Для начала надо определиться с длиной модели. В случае, если обучающее множественное выравнивание не имеет вставок/делеций это тривиально. Наличие же вставок/делеций требует различать вставки и делеции. Простейшее правило если колонка содержит больше половины вставок, то она не включатся в модель, а события вставок трактуются как вставки в последовательность с соответствующими эмиссионными вероятностями. Если выравнивание толстое, то для параметров можно использовать обычные оценки: a kl = A kl / l' A kl' ; e k (a) = E k / a' E k (a');

128 Для тонких выравниваний Простейшие варианты псевдоотсчетов: –Правило Лапласа: к каждому счетчику прибавить 1: e k (a) = (E k (a) +1) / ( a' E k (a')+ N α ); где N α – размер алфавита (20) – Добавлять псевдоотсчеты, пропорционально фоновым частотам: e k (a) = (E k (a) +Aq a ) / ( a' E k (a')+ A); A N α ; Такие псевдоотсчеты соответствуют Байесовой оценке P(θ | D) = P(D | θ) P(θ) / P(D) ; при априорном распределении P(θ) – распределение Дирихле с параметром α a = Aq a.

129 Смеси Дирихле Представим себе, что на распределение вероятностей влияют несколько источников – частота встречаемости символа в белках вообще, частота встречаемости символа в петлях, частота встречаемости символа в трансмембранных сегментах и т.п. Каждое такое распределение дает свои псевдоотсчеты α k. Тогда для вероятности эмиссии можно написать: e k (a) = d P(d| E k ) (E k (a) + α d a ) / ( a' E k (a')+ α d a' ); где P(d| E k ) – вероятность выбора распределения d при условии наблюдаемых частот: P(d| E k ) = P(E k | d) P(d) / d' P(E k | d') P(d') ; Для оценки P(E k | d) используют простую формулу: ( a E k (a) )! Γ( a (E k (a) + α d a ) ) Γ( α d a ) P(E k | d)= a E k (a)! a Γ (E k (a) + α d a ) a Γ (α d a )

130 Использование матрицы замен Еще один способ введения псевдоотсчетов. У нас есть матрица замен аминокислотных остатков (например, PAM120). Матрица замен может трактоваться как то, что каждая аминокислота является немножко другой аминокислотой. Поэтому в качестве псевдоотсчетов используют величину α ia = A b f ib P(a | b), где f ib – частота встречаемости в колонке буквы b, P(a | b) – вероятности замены буквы b на a

131 Использование предка Все последовательности x k в выравнивании произошли от общего предка y. P(y j =a | alignment)=q a k P(x k j |a) / a' q a k P(x k j |a) Тогда для оценки эмиссионной вероятности e j (a) = a' P j (a| a') P(y j =a' | alignment) где P j (a| a') – матрица замен. Матрица замен зависит от скорости эволюции соответствующей колонки. Для выбора матрицы можно использовать принцип максимального правдоподобия: P(x j 1, x j 2,…, x j N ) = a' q a k P(x k j | a, t) max ; Для матрицы замен можно использовать выражение: P(a|b, t) = exp( t P(a|b, 1) )

132 А чему же равно A? Для компенсации малости выборок используют псевдоотсчеты. Разные подходы дают разные распределения псвдоотсчетов α i, но не определяют величину коэффициента A при α i. Часто предполагают, что псевдоотсчеты должны быть сопоставимыми с точностью определения частот Δ, которая пропорциональна Δ N, где N – количество испытаний (толщина выравнивания) поэтому полагают: A=κ N, κ 1 (0.5…1);

133 Множественное выравнивание

134 Способ написать несколько последовательностей друг под другом (может быть с пропусками) так, чтобы в одной колонке стояли гомологичные позиции. "Золотой стандарт" – совмещенные пространственные структуры гомологичных белков. Соответствующие позиции в разных последовательностях отвечают гомологичным позициям Задача. Найти способ (алгоритм и параметры), выравнивающий последовательности "золотого стандарта" правильно. Есть надежда, что в случаях, когда пространственные структуры неизвестны, этот алгоритм правильно выровняет последовательности.

135 Оценка качества множественного выравнивания Энтропийная оценка Обычно считают, что колонки в выравнивании независимы. Поэтому качество выравнивания можно оценить как сумму качеств колонок: S = G + columns S(m k ) G – веса делеций, S(m k ) – вес колонки Пусть с ia – количество появлений аминокислоты a в колонке i. Вероятность колонки можно описать как P(m i ) = a p ia cia Вероятность выравнивания = i P(m i ); В качестве веса можно использовать логарифм вероятности: S = columns S(m k ); S(m k ) = – a c ia log p ia = H(m i ) H(m i ) – энтропия колонки; для вероятностей остатков принимают: p ia = c ~ ia / a' c ~ ia' где c ~ ia – количество остатков в колонке с поправкой на псевдоотсчеты

136 Оценка качества множественного выравнивания Сумма пар Другой традиционный способ оценки – сумма весов матрицы соответствия аминокислотных остатков SP: S(m i ) = k

137 Если есть функционал, то его надо оптимизировать Элементарные переходы: –Сопоставление трех –Сопоставление двух и одна делеция –Делеция в двух последовательностях Seq1 Seq2 Seq3

138 Динамическое программирование для множественного выравнивания Количество вершин равно посл. L i = O(L N ) (почему ?)Количество ребер из каждой вершины = 2 N -1 (почему ?) Количество операций равно T = O(L N ) Надо запоминать обратные переходы в L N вершинах. Если количество последовательностей > 4, то задача практически не разрешима.

139 Прогрессивное выравнивание Строится бинарное дерево (guide tree, путеводное дерево) – листья = последовательности Дерево обходится начиная с листьев. При объединении двух узлов строится парное выравнивание супер- последовательностей (профилей) и получается новая суперпоследовательность x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 x1x1 x2x2 x3x3 x4x4 x5x5 x6x6 x7x7 Final alignment Путеводное дерево строится приближенно – главное быстро. Обычно это кластерное дерево

140 Выравнивание профилей Выравнивание одной стопки последовательности относительно другой – обычное динамическое программирование. Оптимизируется сумма парных весов: i S(m i ) max, S(m i ) = k< lN s(x k i, x l i ) Если мы выравниваем две стопки – 0 < i n и n < i N, то сумму разбиваем на три части: S(m i ) = k< ln s(x k i, x l i ) + n< k< lN s(x k i, x l i ) + kn, n< lN s(x k i, x l i ) Две первые суммы являются внутренним делом стопок, последняя сумма отвечает за сравнение стопок (профилей) При сравнении используем расширенную матрицу сходства, добавив в нее сравнение делционного символа '-' : s(-,-)=0, s( a,-) = -d ; При множественном выравнивании обычно используют линейные штрафы за делеции

141 Взвешивание последовательностей

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

143 Взвешивание последовательностей Способ учета неравномерной представленности последовательностей в выборке называется взвешиванием последовательностей. Каждой последовательности в выравнивании присваивается свой вес β k. Тогда частота каждого символа a в колонке k подсчитывается по формуле: E a k = i β i δ(S i k, a) / β i где S i k – буква в последовательности i в колонке k, β i – вес последовательности i.

144 Взвешивание последовательностей Метод Герштейна – Сонхаммера – Чотьи Пусть нам известно филогенетическое дерево с расстояниями на ветвях. На листьях – последовательности. В начале все веса последовательностей приравниваются длинам веток Далее веса определяем итеративно, внося поправки в веса по ходу движения вверх по дереву: Δw i =t n w i / k-листья ниже узла n w i Смысл заключается в том, что длина ветки распределяется по дочерним узлам w 1 =2 w 2 =2 w 3 =5 w 4 = w 1 =2+3/2=3.5 w 2 =2+3/2=3.5 w 3 =5 w 4 = w 1 = /12=4.4 w 2 = /12=4.4 w 3 =5+35/12=6.25 w 4 = w 1 =4.4 w 2 =4.4 w 3 =6.25 w 4 =8

145 Взвешивание последовательностей Многогранники Вороного Поместим объекты в некоторое метрическое пространство. Каждый объект хочет иметь "поместье" – некоторую область пространства. Отнесём точку пространства x к "поместью" объекта A, если A – самый близкий к x объект. Тогда границы между "поместьями" будут отрезками прямых, проходящих посредине между объектами. В результате все "поместья" будут иметь форму многогранника. Эта конструкция называется многогранниками Воронова. Можно определить вес последовательности как объем поместья. Вопрос только в том, как и в какое метрическое пространство помещать последовательности.

146 Взвешивание последовательностей Многогранники Вороного Один из вариантов метрического пространства – большое количество случайных последовательностей Обычно при генерации случайных последовательностей для взвешивания по методу Вороного i-ая буква каждой последовательности выбирается равновероятно из букв, представленных в i-ой колонке входного выравнивания Метод часто используется, если время работы не слишком важно.

147 Взвешивание последовательностей Максимизация энтропии – метод Хеникофф Пусть k(i,a) – количество остатков типа a в колонке i, m i – количество типов остатков в колонке i. Выберем вес для последовательности k равным w k (i)=1/(m i k(i,a)). Такой вес обеспечивает наиболее равномерное распределение частот остатков в колонке. Чтобы задать вес для последовательности в целом, просуммируем соответствующие веса: w k = i w k (i) = i 1/(m i k(i,a)). Такой вес работает достаточно хорошо и считается быстро. Используется, например, в PSI-BLAST.

148 Обобщенный подход: i H i (w) max, k w k =1; где H i (w) = a p ia log p ia ; p ia – вероятности встречаемости аминокислоты a в колонке i, подсчитанные с учетом весов последовательностей: p ia = k w k δ ( x k i, a); Задача максимизации приводит к системе уравнений: k w k =1; i H i (w)/ w k – λ = 0; Здесь неизвестные w k и неопределенный множитель Лагранжа λ Взвешивание последовательностей Максимизация энтропии

149 ClustalW 1.Строится матрица расстояний с использованием попарных выравниваний. 2.По матрице расстояний строится дерево. 3.Строится прогрессивное выравнивание. Используются дополнительные эвристики: –Взвешивание последовательностей (с учетом только топологии дерева) –На разных уровнях дерева используются разные матрицы сходства – Используется контекстно-зависимые штрафы за открытие делеции –Если при построении выравнивания появляются очень низкие веса, то дерево корректируется Сравните время работы первого и третьего этапов

150 Улучшение выравнивания Недостаток прогрессивных методов: если для некоторой группы последовательностей выравнивание построено, то оно уже не перестраивается. Алгоритм итеративного улучшения 1.Вынимаем из выравнивания одну последовательность 2.По оставшимся последовательностям строим профиль 3.Выравниваем вынутую последовательность с профилем. Фиксируем, иначе ли подровнялась эта последовательность. 4.Переходим к этапу 1. 5.Останавливаемся, если после перебора всех последовательнсотей ничего не изменилось.

151 Улучшение выравнивания Более мощный алгоритм итеративного улучшения 1.Построим по выравниванию дерево 2.Выберем ветвь дерева. Выбор ветви делит выравнивание на две части (последовательности по каждую сторону от ветви). 3.Строим два профиля и выравниваем их друг с другом. Фиксируем, если выравнивание изменилось. 4.Переходим к этапу 2. 5.Заканчиваем, если при переборе всех ветвей ничего не изменилось. Этот алгоритм применён в программе Muscle, за счёт чего достигается преимущество в качестве над ClustalW. Преимущество в скорости достигается за счёт построения матрицы расстояний (см. первый этап ClustalW) не из парных выравниваний, а из сравнений частот слов в последовательностях.

152 Поиск сигналов

153 Постановка задачи Дано несколько (например, 20) последовательностей. Длина каждой последовательности равна 200 В каждой последовательности найти короткий (длиной 20) фрагмент (сайт), такой, что все сайты между собой похожи. Например, даны регуляторные области совместно регулируемых генов. Найти сайты связывания белков-регуляторов.

154 Источник данных ChIP-Chip или ChIP-seq эксперименты SELEX Регуляторные области ортологичных генов Регуляторные области генов, принадлежащих общему метаболическому пути или регуляторной системе.

155 Графовая постановка задачи. Дан многодольный граф: –Каждой доле соответствует последовательность –Вершины – сайты –Ребра проводятся между всеми сайтами, или если эти сайты между собой похожи. На каждой клике графа определено число. Например, информационное содержание безделеционного множественного выравнивания сайтов Задача:Задача: Найти клику наибольшего веса attcgctgac catcgctaac ctttgcaatg

156 HMM-постановка задачи Найти HMM, описывающую наилучший сайт. Для описания сайта используют следующую модель: Start Не сайт x 1 a e a 1 c e c 1 … x 2 a e a 2 c e c 2 … x L a e a L c e c L … End Сайт – p end 1 – p site – p end 1 – p site p site p end p site p end

157 Алгоритм максимизации ожидания (MEME) Допустим, нам приблизительно известна структура сайта. Применяем алгоритм Баума – Велча. Получаем структуру сайта. MEMEАлгоритм MEME: –В качестве исходной модели выбираем модель, индуцированную первым словом в первой последовательности (с учетом псевдоотсетов). –Находим HMM –Берем в качестве исходной следующее слово из первой последовательности. –Так перебираем все слова во всех последовательностях –Отбираем наилучшие HMM

158 Гиббс сэмплер Задача: найти набор позиций сайтов в последовательностях Инициация: В качестве решения выбираем произвольный набор позиций. Итерации: –Удаляем из выборки одну последовательность. –По позициям, определенным для остальных последовательностей строим профиль (HMM). –Для каждой позиции в удаленной последовательности рассчитываем вероятность того, что сайт находится там. –Разыгрываем позицию сайта в удаленной последовательности в соответствии с рассчитанными вероятностями. –Повторяем процедуру много раз для всех последовательностей

159 Вероятности для Гиббс сэмплера Вероятности для Гиббс сэмплера. Позиция разыгрывается с вероятностью, пропорциональной отношению: s(k) s(k) – символ в позиции k f i сигнал (α) f i сигнал (α) – частота появления символа α в позиции i сигнала. Часто используют поправки псевдоотсчетов и взвешивания последовательностей. f фон (α) f фон (α) – фоновая частота появления символа α

160 Дополнительные замечания Сигнал часто имеет структуру – палиндром, повтор. Обычно длина сигнала должна быть заранее известна. Стартуя со случайных сайтов мы можем получить: –Неправильное решение –Решение (сайты), которые по случайным причинам сдвинуты относительно настоящих сайтов

161 RNA

162 Вторичная структура РНК Вторичной структурой называется совокупность спаренных оснований Биологическая роль вторичной структуры: –Структурная РНК – рибосомная, тРНК –Регуляция – Рибопереключатели аттенюация микроРНК –Рибозимы –Стабильность РНК

163 Элементы вторичной структуры Шпилька Спираль Внутренняя петля Множственная петля петля Выпячивание Псевдоузел 5'5'5'5' 3'

164 Способы представления вторичных структур Топологическая схема Круговая диаграмма Массив спаренных оснований Список спиралей from 1 to 1 from 2 to 2 A B C

165 Задача ДанаДана последовательность. НайтиНайти правильную вторичную структуру. Золотой стандарт:Золотой стандарт: тРНК, рРНК. Количество возможных вторичных структур очень велико. Дополнительные ограничения: –Нет псевдоузлов. (На самом деле они очень редки и энергетически невыгодны) Количество возможных структур все равно очень велико оптимальнуюНадо найти оптимальную структуру. А что оптимизировать? Как оптимизировать?

166 Комбинаторный подход Построим граф: –вершины – потенциальные нуклеотидные пары (или потенциальные спирали) –Ребро проводится, если пары совместимы (не образуют псевдоузлов и не имеют общих оснований) Допустимая вторичная структура – клика в этом графе a b c d e fg h a b c d ef g h a b c d ef g h a b c d ef g h a b c d ef g h

167 Структуры без псевдоузлов Структура без псевдоузлов = правильное скобочное выражение Может быть представлено в виде дерева Оценка количества возможных структур: T(L) 1.8 L (очень много) ( ( ( ( ( ( ( ) ) ) ) ( ( ) ) ) ) ) 28,3 22,26 8,12 23,27 6,14 7,13 29,2 30,1

168 Оптимизация количества спаренных оснований Обозначим |s| - мощность структуры (количество спаренных оснований) Пусть s 1 и s 2 две непересекающиеся структуры (структуры без общих оснований) Тогда |s 1 +s 2 | = |s 1 | + |s 2 | s2s2 s1s1 s 1 +s 2

169 Оптимизация количества спаренных оснований Пусть нам известны оптимальные структуры S rt для всех фрагментов i r t j Тогда можно найти оптимальную структуру для сегмента [i, j+1] Для этого нам надо понять, спаривать ли основание j+1, и, если спаривать, то с кем i+1 k S k+1,i S 1,k-1

170 Динамическое программирование для количества спаренных оснований (Нуссинофф) Количество спаренных оснований в оптимальной структуре S * i,j+1 определяется как максимум: S * i,j+1 = max { S * i,j ; (нет спаривания) max k (S * i,k-1 + S * k, j )+1; (k спаривается с j+1) }; Время работы алгоритма: TO(L 3 ) j+1 k S k+1,j S i,k-1 i

171 Динамическое программирование для количества спаренных оснований При поиске оптимального количества спаренных оснований заполняется треугольная матрица весов S i,j, i < j. Обозначим π ij – номер основания, с которым надо спарить основание j при анализе сегмента [i, j], или 0, если не надо спаривать. При оптимизации запоминаем треугольную матрицу спаривания (аналог матрицы обратных переходов)

172 AUCG AU CG GC Энергия вторичной структуры Энергия спиралей Энергия петель (энтропия) A – U C – G A – U G – C C – G ΔG = Энергия спирали рассчитывается как сумма энергий стэкингов =

173 Энергия петель Энергия свободной цепи ΔG = B + 3/2 kT ln L Для шпилек при L=3..5 кроме энтропии есть некоторое напряжение структуры. Для внутренних петель и для мультипетель L – суммарная длина петель + количество ветвей. Параметр B зависит от типа петли Для выпячивания сохраняется стэкинг. Обычно используют не формулу, а таблицы.

174 Минимизация энергии Обычное динамическое программирование не проходит – нет аддитивности. Определения –нуклеотид h называется доступным для пары ij, если НЕ существует спаривания kl, такого, что i < k < h < l < j –Множество доступных нуклеотидов для пары ij называется петлей L ij, а пара ij называется замыкающей парой. Частный случай петли – стэкинг. –Энергия структуры рассчитывается как сумма энергий петель (в том числе и стекингов): ΔG = e(L ij )

175 Алгоритм Зукера Введем две переменные: –W(i,j) – минимальная энергия для структуры на фрагменте последовательности [i, j]; –V(i,j) – минимальная энергия для структуры на фрагменте последовательности [i, j] при условии, что i и j спарены; Рекурсия:

176 Алгоритм Зукера Рекурсия для W требует времени TO(L 3 ) Рекурсия для V требует гораздо большего времени TO(2 L ) Причина – мультипетли. Можно: –Ограничить размер или индекс мультипетель –Применить упрощенную формулу для их энергии –Просматривать мультипетли только если i+1, j-1 не спарены. –Применить приближенную эвристику i,j i l,j l

177 Проблемы минимизации энергии 1.Только около 60% тРНК сворачиваются в правильную структуру 2.Энергетические параметры определены не очень точно. Более того, в клетке бывают разные условия, и, соответственно, реализуются разные параметры. 3.Находится единственная структура с минимальной энергией, в то время как обычно существует несколько структур с энергией, близкой к оптимальной.

178 Решение проблем Искать субоптимальные структуры Искать эволюционно консервативные структуры. –структуры тРНК и рРНК определены именно так