'

Алгоритмы биоинформатики

Понравилась презентация – покажи это...





Слайд 0

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


Слайд 1

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


Слайд 2

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


Слайд 3

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


Слайд 4

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


Слайд 5

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


Слайд 6

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


Слайд 7

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


Слайд 8

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


Слайд 9

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


Слайд 10

Сколько существует выравниваний? Дано: две последовательности S1 и S2 длиной m и n. Сколько есть способов написать одну последовательность под другой (со вставками)? Построим выборочную последовательность S длиной m+n следующим образом: возьмем несколько символов из последовательности S1, потом несколько символов из последовательности S2 потом опять несколько символов из S1, потом опять несколько из S2. Каждой выборочной последовательности S соответствует выравнивание и по каждому выравниванию можно построить выборочную последовательность. (Доказать!) Количество выборочных последовательностей равно Nsel = Cn+mm =(m+n)! / (m!*n!) (Доказать!) (2n)! 22n Nalgn = C2nn = ? (n!)2 v?n


Слайд 11

Динамическое программирование для редакционного расстояния Граф редакционного расстояния для последователь-ностей S1,S2: вершина vi,j соответствует префиксам последовательностей {S11..i}, {S21..j}. На вершине записано редакционное расстояние между префиксами. (красные стрелки соответствуют вставкам и удалениям) di,j di+1,j di,j+1 di+1,j+1 di+1,j+1= min{ di+1,j+1, di,j+1+1, di,j+ei,+1,j+1} ei,j={ 0, S1i = S2j ; 1, S1i ? S2j }


Слайд 12

Подмена задачи и обобщение Заменим расстояния di,j на -di,j. Тогда операцию min надо заменить на max. Прибавим к -di,j ? (wi,j= ? - di,j ), тогда получим функцию сходства: совпадение = ?, замена = -?, делеция = -1. Функцию сходства W легко обобщить, варьируя штрафы за замену и делеции. Новая задача: написать одну последовательность под другой так, чтобы максимизировать сходство


Слайд 13

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


Слайд 14

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


Слайд 15

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


Слайд 16

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


Слайд 17

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


Слайд 18

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


Слайд 19

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


Слайд 20

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


Слайд 21

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


Слайд 22

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


Слайд 23

Алгоритм Смита-Ватермана Точка конца пути (от нее начинаем обратный просмотр и восстановление пути) определяется так: (imax, jmax) = argmax (wi,j) wi,j = max { wi-i,j-1 + ei,j , i > 1, j > 1 wi-1,j – d , i > 1 wi,j-1 – d , j > 1 0 } Пусть (при одинаковых параметрах) мы получили вес глобального выравнивания wglob и вес локального выравнивания wloc. Какая величина больше?


Слайд 24

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


Слайд 25

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


Слайд 26

Аффинные штрафы за делецию Вместо логарифмической зависимости используют зависимость вида: ? ( l ) = dopen+ l dext dopen – штраф за открытие делеции dext – штраф за удлинение делеции ? l dopen dext a log( l + 1 )


Слайд 27

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


Слайд 28

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


Слайд 29

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


Слайд 30

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


Слайд 31

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


Слайд 32

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


Слайд 33

Серия матриц PAM Находим выравнивания, отвечающие расстоянию PAM1 Находим частоты пар и вычисляем частоты пар: p(??) = p(? > ?) p(?)+ p(? > ?) p(?) полагая p(? > ?) = p(? > ?) получаем p(? > ?) = p(??) / (p(?)+ p(?)) p(? > ?) = 1 – ? ??? p(? > ?) PAMN(??) = log (p (? > ?)N / p?p? )


Слайд 34

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


Слайд 35

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


Слайд 36

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


Слайд 37

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


Слайд 38

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


Слайд 39

Наибольшая общая подпоследовательность Длина оптимальной подпоследовательности есть случайная величина r(n), зависящая от длины последовательностей. Пусть две последовательности длиной n разбиты на два фрагмента длиной n1 и n2 (n1+n2=n) Ясно, что оптимальная подпоследовательность будет не хуже, чем объединение оптимальных подпоследовательностей для фрагментов: r(n) ? r1(n1)+r2(n2) (попробуйте понять смысл неравенства) Отсюда следует, что математическое ожидание M(r(n)) ? M(r(n1)) + M(r(n2)), или M(r(n)) ? c •n Можно показать, что M(r(n)) – M(r(n1)) + M(r(n2) > 0 Поэтому: r(n) r1(n1) r2(n2) M(r(n)) ? c • n, (n > ?)


Слайд 40

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


Слайд 41

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


Слайд 42

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


Слайд 43

e-value & p-value Количество независимых локальных выравниваний с весом >S описывается распределением Пуассона (Karlin &Altschul) : E(S) = Kmn e –?S где ? – положительный корень уравнения ? p?p? e ?s(??) = 1, s(??) – матрица замен K – константа, зависящая от p? и s(??). e-value: E(S) – ожидаемое количество выравниваний с заданным весом p-value: p(x>S) = 1- e –E(S) – Вероятность встретить выравнивание с таким или большим весом


Слайд 44

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


Слайд 45

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


Слайд 46

Поиск по банку. FASTA. Используется техника поиска якорей с помощью хэш-таблицы. Два якоря (i1,j1), (i2,j2) принадлежат одной диагонали, если i1 – j1 = i2 – j2 Мощностью диагонали называется количество якорей, принадлежащих диагонали. Иногда в мощность диагонали включают мощности соседних диагоналей (чтобы учесть возможность делеций) Отбираем n* (n*=10) самых мощных диагоналей и для них пытаемся построить цепочки якорей, или строим S-W выравнивание в полосе (Wilbur-Lipman-Pearson) Для оценки стат.значимости используют z-score


Слайд 47

Поиск по банку. BLAST1. Ищем якоря с помощью хэш-таблицы Каждый якорь расширяем с тем, чтобы получить сегмент совпадения наибольшего веса (HSP – high scoring pair). Оцениваем его статистическую значимость, и, если она больше порога, то репортируем Для оценки значимости используется формула Альтшуля (Altschul, Lipman, Pearson)


Слайд 48

Поиск по банку. BLAST2. T-соседней l-граммой LT для l-граммы L называется такая l-грамма, что вес ее сравнения с L не меньше заданного T: ?s(Li, LiT) ? T Для аминокислотных последовательностей при просмотре запроса формируем не только те l-граммы, которые встретились в нем, но также все T-соседние l-граммы. Характерное количество l-грамм для белка длиной 300 остатков – 15000. Расширяются только те якоря, которые принадлежат мощной диагонали (как в FASTA), причем мощность диагонали должна быть ? 2T При расширении диагонали допускается небольшое количество делеций


Слайд 49

Быстрое выравнивание Ищем якоря с помощью хэш-таблицы Якорь (i1,j1) предшествует якорю (i2,j2), если i1 < i2 & j1 < j2 & i2 – i1 < d & j2 – j1 < d Получаем ориентированный граф с небольшим количеством вершин и ребер Можно найти оптимальную цепочку якорей методом динамического программирования


Слайд 50

Введение в Байесову статистику


Слайд 51

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


Слайд 52

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


Слайд 53

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


Слайд 54

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


Слайд 55

Определения Пусть у нас есть несколько источников 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)


Слайд 56

Пример Пусть есть две кости – правильная и кривая (с вероятностью выпасть 6 – 0.5). И пусть нам подсовывают кривую кость с вероятностью 1%. Мы бросили кость 3 раза и 3 раза получили 6. Какова вероятность того, что нам дали кривую кость? P(кривая кость | 3 шестерки) = P(3 шестерки | кривая кость) • P(кривая кость) P(3 шестерки) P(3 шестерки)=P(3 шестерки | кривая кость) • P(кривая кость) + P(3 шестерки | правильная кость) • P(правильная кость) = 0.53 • 0.01 + (1/6)3 •0.99 = 0.00125+0.0046 = 0.00585 P(кривая кость | 3 шестерки) = 0.00125/0.00585=0.21 Вывод – кость скорее правильная!


Слайд 57

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


Слайд 58

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


Слайд 59

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


Слайд 60

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


Слайд 61

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


Слайд 62

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


Слайд 63

Описание HMM Пример с монетой можно представить в виде схемы конечного автомата: Прямоугольники означают состояния Кружки означают результат бросания (эмиссии) Стрелки – возможные переходы между состояниями Числа около кружков – вероятности эмиссии ei числа около стрелок – вероятности переходов между состояниями aik 0+ 1+ 0.5 0.5 0- 1- 0.9 0.1 0.9 0.1 0.8 0.2 Сумма весов исходящих стрелок равна 1 Сумма весов эмиссии в каждом состоянии рана 1


Слайд 64

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


Слайд 65

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


Слайд 66

Viterbi рекурсия Обозначения vk(i) – наилучшая вероятность пути, проходящего через позицию i в состоянии k. ?k(i) – наилучший переход из позиции i в состоянии k в предыдущую позицию (предыдущее состояние) ?*(i) – наилучшее состояние в позиции i Инициация vk(0) = ?(0,k); k – номер состояния Рекурсия vk(i) = ek(xi) maxm( vm( i-1 ) amk); ?(i,k) = argmaxm( vm( i-1 ) amk); обратный переход Завершение P(x,?*)= maxm( vm( L ) am0); ?*(L) = argmaxm( vm ( L ) am0); Оптимальный путь ?*( i-1 ) = ? ( i, ?* ( i ) ) ;


Слайд 67

Другая постановка задачи Для каждого наблюденного значения определить вероятность того, что в этот момент монета была правильной. Для этого надо просуммировать по всем путям, проходящим через точку i+ вероятности этих путей. Для решения этой задачи достаточно вспомнить динамическое программирование над полукольцом с использованием операции сложения и умножения. Оцениваем значение P(x, ?i=k) = P(x1…xi, ?i=k) •P(xi+1…xL | ?i=k)/P(x); Первый сомножитель fk(i) = P(x1…xi, ?i=k) определяем просмотром вперед Второй сомножитель bk (i+1) = P(xi+1…xL | ?i=k) определяем просмотром назад


Слайд 68

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


Слайд 69

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


Слайд 70

Оценка параметров HMM при наличии обучающей выборки Можно показать, что при большом количестве наблюдений справедливы оценки akl = Akl / ?l'Akl' ; ek(b) = Ek(b) / ?b'Ek(b); Akl – наблюденное количество переходов между моделями Ek(b) – количество порожденных символов в соответствующих моделях При малых размерах выборки используют технику псовдоотсчетов, добавляя к наблюденным значениям некоторое количество шума.


Слайд 71

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


Слайд 72

Оценки параметров по Бауму-Велчу Имея заданные параметры модели можно определить вероятность перехода между состояниями: P(?i=k, ?i+1=l | x,?) = fk(i)aklei(xi+1)bl(i+1)/P(x), где fk(i) = P(x1…xi, ?i=k), bl(i+1) •P(xi+1…xL | ?i+1=l) – значения, полученные при прямом и обратном проходе. Тогда для переходных и эмиссионных вероятностей получим оценки для количества переходов и порожденных символов: Akl= ?j1/P(xj) ?i f jk(i)aklei(xi+1)b jl(i+1); Ek(b) = ?j1/P(xj) ?i f jk(i) b jk(i) ?(x ji,b); где x j – j-последовательность в выборке, f jk , b jl – результаты прямого и обратного прохода по последовательности x j


Слайд 73

Предсказание кодирующих областей в прокариотах Реальная схема HMM для поиска кодирующих областей сложнее: Включает в себя SD сайт Учитывает неравномерность следования кодонов A C G T eA eC eG eT A T 1 1 A A C A A G A A T A A A 1 1 Кодоны pcodon G 1 1 1 T G A 1 1 T A A T A G Стоп pstart pstop 1 1-pstart Старт некодирующая последовательность


Слайд 74

Оценка качества обучения Выборку разбивают на два подмножества – обучающую и тестирующую На первой выборке подбирают параметры На второй – тестируют и определяют качество обучения: TP – количество правильно определенных позитивных позиций (например, кодирующих) TN – количество правильно определенных негативных позиций (например, некодирующих) FP – количество неправильно определенных позитивных позиций (некодирующих, предсказанных как кодирующие) FN – количество неправильно определенных негативных позиций (кодирующих некодирующих, предсказанных как некодирующие) Специфичность: Sp = TP / (TP + FP) Чувствительность: Sen = Качество QQ = Коэффициент корреляции CC =


Слайд 75

Профили


Слайд 76

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


Слайд 77

Энтропия колонки Пусть колонка содержит 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 ? ) называется информационным содержанием колонки


Слайд 78

HMM профиль Модель: каждая последовательность множественного выравнивания является серией скрытой Марковской модели. Профиль – описание Марковской модели. Каждой позиции соответствует свое состояние. Вероятности переходов между соседними состояниями равны 1. Вероятность того, что некоторая последовательность x соответствует профилю M: P( x | M)= ? ei (xi); Значимость определяется отношением правдоподобия: сравнением с P( x | R) – вероятностью, что последовательность сгенерирована случайной моделью R: S = log (P( x | M) / P( x | R)) = ? log {ei (xi) / q (xi)}; Величины wi(?)= log {ei (?) / q (?)} называют позиционной весовой матрицей (PSSM) B eA ec ef … eA ec ef … eA ec ef … eA ec ef … E M1 M2 M3 Mn


Слайд 79

HMM с учетом возможности вставок Делеция в профиле и в последовательности могут идти подряд (в отличие от парного выравнивания) Делеционные состояния – молчащие (не имеют эмиссии) Вероятность перехода в делеционное состояние зависит от позиции B E M1 M2 M3 Mn D1 D2 D3 Dn Делеция в профиле Делеция в последовательности


Слайд 80

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


Слайд 81

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


Слайд 82

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


Слайд 83

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


Слайд 84

Использование предка Все последовательности xk в выравнивании произошли от общего предка y. P(yj=a | alignment)=qa?kP(xkj|a)/?a' qa?kP(xkj|a) Тогда для оценки эмиссионной вероятности ej (a) = ?a' Pj(a| a') P(yj=a' | alignment) где Pj (a| a') – матрица замен. Матрица замен зависит от скорости эволюции соответствующей колонки. Для выбора матрицы можно использовать принцип максимального правдоподобия: P(xj1, xj2,…, xjN) = ?a' qa?kP(xkj| a, t) > max ; Для матрицы замен можно использовать выражение: P(a|b, t) = exp( t P(a|b, 1) )


Слайд 85

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


Слайд 86

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


Слайд 87

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


Слайд 88

Взвешивание последовательностей Метод Герштейна-Сонхаммера-Чотьи Пусть нам известно филогенетическое дерево с расстояниями на ветвях. На листьях – последовательности. В начале все веса последовательностей приравниваются длинам веток Далее веса определяем итеративно, внося поправки в веса по ходу движения вверх по дереву: ?wi=tn wi/ ?k-листья ниже узла n wi Смысл заключается в том, что длиня ветки распределяется по дочерним узлам


Слайд 89

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


Слайд 90

Вероятность модели при условии, что последовательность x принадлежит модели M: P( x | M) P(M) P(M | x) = P( x | M)P(M) + P( x | R)(1-P(M) Вероятность модели при условии нескольких последовательностей (дискриминатор): D=?k P( M | xk) Веса последовательностей подбираем так, чтобы максимизировать D. Но: чтобы вычислить D нам необходимы параметры модели, которые зависят от того, как мы взвешивали последовательности. Применяют итеративную процедуру. Взвешивание последовательностей Максимально дискриминирующие веса


Слайд 91

Взвешивание последовательностей Максимизация энтропии Пусть k(i,a) – количество остатков типа a в колонке i, mi – количество типов остатков в колонке i. Выберем вес для последовательности k равным wk(i)=1/(mi k(i,a)). Такой вес обеспечивает наиболее равномерное распределение частот остатков в колонке. Чтобы задать вес для последовательности в целом, просуммируем соответствующие веса: wk = ?i wk(i) = ?i 1/(mi k(i,a)).


Слайд 92

Обобщенный подход: ?i Hi(w) > max, ?kwk=1; где Hi(w) = ?a pia log pia; pia – вероятности встречаемости аминокислоты a в колонке i, подсчитанные с учетом весов последовательностей: pia= ?k wk ? ( xki, a); Задача максимизации приводит к системе уравнений: ?kwk=1; ?i ? Hi(w)/ ?wk – ? = 0; Здесь неизвестные wk и неопределенный множитель Лагранжа ? Взвешивание последовательностей Максимизация энтропии


Слайд 93

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


Слайд 94

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


Слайд 95

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


Слайд 96

Оценка качества множественного выравнивания Сумма пар Другой традиционный способ оценки – сумма весов матрицы соответствия аминокислотных остатков SP: S(mi) = ?k<l s(xki , xli); Способ не совсем правильный. Более правильная оценка для трех последовательностей S(mi) = log (pabc / qaqbqc), а не log (pab/qaqb) + log (pbc/qbqc) + log (pac/qaqc); (вспомним определение матрицы замен)


Слайд 97

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


Слайд 98

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


Слайд 99

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


Слайд 100

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


Слайд 101

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


Слайд 102

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


Слайд 103

Множественное выравнивание с помощью HMM Каждому множественное выравнивание соответствует скрытая Марковская модель. Можно применить алгоритм максимизации ожидания Баума-Велча: Порождаем случайные параметры HMM. Выравниваем все последовательности с этой моделью Переоцениваем параметры. Проблема: легко попасть в локальный максимум Обход проблемы: время от времени параметры HMM возмущаются. Другой вариант – использование искусственного отжига. Достоинство подхода: одновременно анализируются все последовательности. Нет проблемы необратимости, характерной для прогрессивного выравнивания.


Слайд 104

Блочное выравнивание


Слайд 105

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


Слайд 106

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


Слайд 107

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


Слайд 108

HMM-постановка задачи Найти HMM, описывающую наилучший сайт. Для описания сайта используют следующую модель: Start Не сайт x1 a ea1 c ec1 … x2 a ea2 c ec2 … xL a eaL c ecL … End Сайт 1 1 1 1- pend 1 – psite- pend 1 - psite psite pend psite pend


Слайд 109

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


Слайд 110

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


Слайд 111

Вероятности для Гиббс сэмплера Вероятности для Гиббс сэмплера


Слайд 112

Комбинаторные методы


Слайд 113

RNA


Слайд 114

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


Слайд 115

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


Слайд 116

Способы представления вторичных структур Топологическая схема Круговая диаграмма Массив спаренных оснований Список спиралей


Слайд 117

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


Слайд 118

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


Слайд 119

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


Слайд 120

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


Слайд 121

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


Слайд 122

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


Слайд 123

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


Слайд 124

Восстановление структуры по матрице спаривания - - 17 - 2 2 16 - 3 15 - 11 11 11 - - - - - - - 14 - 9 9 9 9 9 9 13 - 10 10 10 10 10 10 - 12 - 11 - 10 - 6 9 - 4 8 - 5 7 - 6 - 5 - 4 - 3 - 2 - 1 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1


Слайд 125

Восстановление структуры по матрице спаривания SearchStruct (int i, int j) { int i0=i, j0=j; do{ if(i >= j) return; if(?ij == 0) j--; if(?ij != i) i++; if(?ij == i) { StorePair(i,j); SearchStruct (i0, i-1); SearchStruct (i+1, j0); return; } }while(true) }


Слайд 126

Энергия вторичной структуры Энергия спиралей Энергия петель (энтропия) A – U C – G A – U G – C C – G ?G = -3.2 -3.2 -3.7 -4.5 Энергия спирали рассчитывается как сумма энергий стэкингов = - 14.6


Слайд 127

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


Слайд 128

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


Слайд 129

Алгоритм Зукера Введем две переменные: W(i,j) – минимальная энергия для структуры на фрагменте последовательности [i, j]; V(i,j) – минимальная энергия для структуры на фрагменте последовательности [i, j] при условии, что i и j спарены; Рекурсия: V(i,j) = mini?i1<j1<i2<j2<…<ik<jk?j ?lk V(il ,jl) W(i,j)=min{ W(i+1,j), i не спарено W(i,j-1), j не спарено V(i,j), i и j спарены mini<k<j(W(i,k) + W(k+1,j)) } i j спарены с кем-то.


Слайд 130

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


Слайд 131

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


Слайд 132

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


Слайд 133

Поиск субоптимальных структур и структурных элементов Статистическая сумма Z = ? exp(-?Gi /kT) Если мы просуммируем по всем структурам, содержащим данную пару, то мы можем оценить ее значимость (чем Z больше, тем более значимым является спаривание) Для подсчета Z можно использовать тот же алгоритм динамического программирования, заменив min на суммирование, а сложение на умножение. Больцмановская вероятность того, что нуклеотиды i,j спарены равна: P(i,j) = exp(-?Gij /kT) /Zij ; Разыгравыем пары оснований в соответствии с этой вероятностью и восстанавливаем соответствующие субоптимальные вторичные структуры.


Слайд 134

Консенсусные вторичные структуры РНК


Слайд 135

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


Слайд 136

Метод ковариаций Пусть дано множественное выравнивание последовательностей Взаимная информация двух колонок: I(A,B) = ? ?? fAB(??) log2{fAB(??) / (fB(?))} fAB(??) – частоты одновременной встречи буквы ? в колонке A и буквы ? в колонке B. fA(?) – частота встречаемости буквы ? в колонке A. fB(?) – частота встречаемости буквы ? в колонке B. Пары колонок с высоким значением взаимной информации с большой степенью вероятности образуют комплементарную пару (если высоки совместные частоты для пар букв AT, CG) Для восстановления вторичной структуры можно использовать алгоритм Нуссинофф, приписывая в качестве весов пар значение взаимной информации.


Слайд 137

Грамматики Определения Терминальным символом называется символ, который может получаться в строке (обозначается малыми буквами) Нетерминальный символ – символ для обозначения промежуточной подстроки Грамматика – набор правил генерации слов Пример: W2 > aW1 , W1 > bW2 , W1 > ? ; Порождает слова вида "ababababab"


Слайд 138

Стохастические контекстно-свободные грамматики Контекстно-свободные грамматики имеют правила вида W > ? ? – терминальные и/или нетерминальные исключая нулевую строку Правила преобразования могут быть снабжены вероятностями Обобщает скрытые Марковские модели. Позволяет описывать вторичные структуры РНК. Пример. Грамматика S > aW1t | cW1g | gW1c | tW1a; W1 > aW2t | cW2g | gW2c | tW2a; W2 > aW3t | cW3g | gW3c | tW3a; W3 > gaaa | gcaa Порождает шпильки с длиной спирали 3 и с последовательностью в петле gaaa или gcaa


Слайд 139

Задача выравнивания СКСГ с последовательностью СКСГ для описания вторичной структуры. Есть шесть типов преобразований:


Слайд 140

Общая модель для выравнивания вторичной структуры с последовательностью


Слайд 141

Общая идея алгоритма разбора последовательности Заполняется трехмерная матрица A(i,j,v): i,j – рассмотренный сегмент, v – тип вершины (MP …) Просмотр начинается "изнутри" с коротких фрагментов. для каждого сегмента вероятность того, что этот сегмент может быть порожден соответствующей грамматикой. (Вариант динамического программирования) Затем просмотром "внутрь" находится способ вложения последовательности в грамматику MP MP ML MP MP MP S


Слайд 142

Поиск генов


Слайд 143


×

HTML:





Ссылка: