Выравнивание последовательностей

В мышцах человека есть миоглобин, в клубеньках на корнях сои — леггемоглобин. Оба держат гем, группу с атомом железа, и связывают кислород, а их пространственные структуры похожи: цепь свёрнута в спирали вокруг гема. Но если сопоставить последовательности буква к букве, в лучшем куске совпадёт 28 букв из 140, и лучшее сопоставление набирает не больше очков, чем у пары, где одну цепочку случайно перемешали. Выравнивание последовательностей — способ найти лучшее сопоставление двух цепочек букв и понять, говорит ли оно о родстве. Ниже — α- и β-цепи гемоглобина человека. В α-цепи нет одной из спиралей β-цепи, и выравнивание находит это место, ничего не зная о спиралях.

В списке есть и другие пары: миоглобины кашалота и человека, леггемоглобин, лизоцим из яичного белка с α-лактальбумином из молока и пара, у которой родства нет. Сравните, как растёт доля плюсов, когда одинаковых букв становится меньше: у дальней родни буквы чаще заменены похожими, чем случайными.

Идея

Выравнивание ставит две строки друг под другом и вставляет в них пропуски — дефисы, — чтобы похожее встало напротив похожего. Каждый столбец получает очки: пара букв — оценку из матрицы замен, буква против дефиса — штраф. Лучшее выравнивание — с наибольшей суммой.

Это та же таблица, что у редакционного расстояния, только вместо цены правок, которую уменьшают, здесь сходство, которое увеличивают. Клетка (i, j) хранит лучший счёт выравнивания первых i букв строки a с первыми j буквами строки b, и последний столбец этого выравнивания бывает трёх видов:

F(i, j) = max { F(i − 1, j − 1) + s(a_i, b_j);  F(i − 1, j) − d;  F(i, j − 1) − d },   F(i, 0) = −i·d,  F(0, j) = −j·d.

Здесь s — оценка пары букв, d — штраф за букву против пропуска. Ответ лежит в правом нижнем углу, а само выравнивание восстанавливают обратным ходом: из каждой клетки идут туда, откуда пришёл максимум. Так выравнивают две последовательности целиком, от первой буквы до последней, — это глобальное выравнивание, алгоритм Нидлмана — Вунша.

Почему сумма? Очки в матрицах замен — логарифмы. Оценка пары букв говорит, во сколько раз чаще эта пара стоит в одном столбце у родственных белков, чем у случайно взятых, а сумма логарифмов — логарифм произведения. Тогда счёт выравнивания показывает, насколько вероятнее, что такие столбцы оставила общая история, а не случай, и лучшее выравнивание — самая правдоподобная история замен. Для пропусков такого точного смысла нет: их штрафы подбирают так, чтобы находились настоящие родственники.

Откуда берутся очки

В 1992 году Стивен и Джорджа Хеникофф собрали около 2000 блоков — кусков выравнивания без пропусков из более чем 500 групп родственных белков — и посчитали, как часто каждая пара аминокислот стоит в одном столбце. Чтобы почти одинаковые последовательности не перевесили остальные, внутри блока те, что совпадают не меньше чем на 62 %, сливали в одну группу; отсюда число в названии BLOSUM62. Оценка пары x, y — это 2 · log₂(p_xy / q_xq_y), округлённое до целого: частота пары в блоках против ожидаемой, если бы буквы стояли рядом случайно, в половинках бита.

Из этого видно, почему числа такие. Совпадение двух триптофанов W стоит +11: это одна из самых редких аминокислот, и случайно две W друг напротив друга встают редко. Два лейцина L — только +4: лейцин частый. Замена изолейцина I на валин V стоит +3, аспарагиновой кислоты D на глутаминовую E +2, лизина K на аргинин R +2: это пары похожих по свойствам аминокислот, и белок такую замену обычно переносит. Лейцин против аспарагиновой кислоты — −4. Средняя оценка двух случайных букв отрицательна, −0,52: в файле матрицы это строка «Expected = −0.5209». Это важно для локального выравнивания.

Раньше пользовались матрицами Маргарет Дейхофф (PAM, 1978): их считали по близким родственникам и пересчитывали на большие расстояния моделью эволюции. Хеникоффы считали замены прямо по далёкой родне, и, как они пишут, поиск и выравнивание заметно улучшились. Для ДНК матрица обычно проще: одна оценка за совпадение и одна за замену.

Пропуски: дорого открыть, дёшево продлить

Одна мутация может вставить или удалить сразу несколько букв. Три пропуска по одной букве — следы трёх событий, а пропуск из трёх букв подряд — одного. Штраф d за каждую букву их не различает. Поэтому пропуск длины k обычно оценивают как o + e·k: дорого открыть, дёшево продлить. Такой штраф называют аффинным.

Для произвольной зависимости штрафа от длины таблицу можно считать, перебирая в каждой клетке длину последнего пропуска, — это уже O(m·n·(m + n)). В 1982 году Осаму Гото показал, что для аффинного штрафа хватает трёх таблиц и прежнего времени O(m·n). Таблица M — лучшие выравнивания, которые кончаются парой букв, X — буквой строки a против пропуска, Y — пропуском против буквы строки b:

M(i, j) = s(a_i, b_j) + max { M, X, Y в клетке (i − 1, j − 1) }
X(i, j) = max { M(i − 1, j) − o − e;  X(i − 1, j) − e;  Y(i − 1, j) − o − e }
Y(i, j) = max { M(i, j − 1) − o − e;  Y(i, j − 1) − e;  X(i, j − 1) − o − e }

Таблица X помнит, что пропуск уже открыт: продолжить его из X стоит только e, а войти в X из другой таблицы — o + e. Одной таблицы для этого мало: в ней на клетку приходится одно число, и уже не узнать, каким столбцом кончилось лучшее выравнивание.

В первой сцене это видно на гемоглобине. Когда каждая буква пропуска стоит 8, пропуски в α-цепи рассыпаются на куски 1, 1, 2 и 3, счёт 259. При штрафе 10 + 1 за букву их становится три — 1, 1 и 5, счёт 281, — и пять подряд встают напротив остатков 50–54 β-цепи. Это спираль D, по структуре — остатки 50–56, которой нет ни у одной известной α-цепи позвоночных. Алгоритм ничего не знает о спиралях: он нашёл только, что один длинный пропуск дешевле нескольких коротких. Цветные полосы в сцене показывают, что и остальные спирали, от A до H, легли друг напротив друга.

Цены пропусков — часть модели, а не закон природы. В BLAST для BLOSUM62 по умолчанию открытие стоит 11, а не 10. При 11 + 1 за букву оба варианта набирают ровно 277: пропуск 1 + 5 и пропуск из шести букв подряд, при котором буквы D-L-S α-цепи, совпадавшие с D-L-S β-цепи, сдвигаются на пять позиций. Наш листинг при равенстве выбирает шесть подряд. Оптимум модели — не обязательно истина; модель говорит лишь, что эти два варианта для неё одинаково хороши.

Здесь три таблицы для α- и β-цепей сложены стопкой. Клетка светится, если через неё проходит выравнивание почти такое же хорошее, как лучшее: так видно не только оптимальный путь, но и равноценные и близкие к нему. Чтобы это посчитать, таблицы считают дважды — с начала строк и с конца — и складывают: лучший счёт через клетку равен сумме двух половин пути. Путь идёт по листу M, а там, где у α-цепи пропуск, спускается на лист Y и возвращается. Двигайте ползунок вдоль диагонали и переключите цену пропусков на 8 за букву: путь спустится на лист Y уже трижды, а на листах X и Y засветится почти вся полоса вдоль диагонали. Когда несколько коротких пропусков стоят столько же, сколько один длинный, почти таких же хороших вариантов с пропусками в разных местах становится много.

Локальное выравнивание

Глобальное выравнивание тянет обе последовательности от начала до конца. Если у двух белков общий только один домен, а у двух геномов — только часть генов, ему приходится выравнивать и всё чужое. В 1981 году Темпл Смит и Майкл Уотерман изменили одну строку: к максимуму в клетке добавляется ноль.

M(i, j) = s(a_i, b_j) + max { 0;  M, X, Y в клетке (i − 1, j − 1) }

Ноль значит «начать заново здесь»: если всё, что было левее и выше, в сумме хуже нуля, его отбрасывают. Ответ — наибольшее значение во всей таблице, а обратный ход идёт от этой клетки до первой, где выравнивание началось с нуля. Так находится самый похожий кусок одной последовательности и самый похожий кусок другой.

Здесь и нужна отрицательная средняя оценка случайной пары. Если бы случайные буквы в среднем приносили очки, выгодно было бы тянуть любое выравнивание до самых концов, и локальное ничем не отличалось бы от глобального. У BLOSUM62 среднее −0,52, поэтому за пределами настоящего сходства счёт тает и кусок обрывается.

Фаги λ, HK97 и HK022 — вирусы кишечной палочки из одной группы, лямбдоидных фагов. Сцена выравнивает каждый белок одного фага с каждым белком другого: для λ и HK97 это 4380 пар и 177 миллионов клеток таблиц. Родственными считаются пары, у которых E-value меньше 10⁻⁴, — таких 28, и у одиннадцати совпадает 90 % букв и больше.

Картина вышла не такой, какую ждёшь от родни. Из десяти белков головки λ — упаковка ДНК, портал, капсид — ни у одного нет родственника в HK97. Из пятнадцати белков хвоста родню нашли десять, но далёкую: совпадает 22–39 % букв. А рядом стоят почти одинаковые белки: интеграза, которой фаг встраивается в хромосому бактерии, совпадает в 354 буквах из 356, репрессор cI — в 235 из 237, белки выхода из клетки S и R — в 102 из 107 и 156 из 158. Переключите на HK97 и HK022: у них, наоборот, общая головка — все шесть белков, у основного белка капсида совпадают все 385 букв, — а репрессоры не родственны вовсе.

Отрезок показывает, какой кусок белка нашёл пару. У большого белка хвоста J из 1132 букв общий с HK97 кусок кончается за 146 букв до конца. Как показали в 2000 году, последние 249 букв белка J — то место, которым λ цепляется за рецептор LamB на поверхности клетки. Коснитесь этого отрезка: у HK97 для этого конца пары нет.

Так выглядит мозаика. В 2000 году группа Роджера Хендрикса прочитала геномы HK97 и HK022 и, сравнив их с λ и с фагом P22, пришла к выводу: лямбдоидные фаги собраны из кусков с резкими границами между ними, и границы чаще всего проходят по краям генов, иногда — по границам доменов внутри гена. Фаги обмениваются генами, и каждый — сборка из деталей разного происхождения. Глобальное выравнивание целых геномов этого не увидело бы: оно растянуло бы одно выравнивание на сорок тысяч букв и смешало родственные куски с чужими.

Родство или случайность

У любых двух последовательностей есть лучшее выравнивание, и у него есть счёт. Вопрос в том, больше ли этот счёт, чем даёт случай. Нидлман и Вунш в 1970 году проверяли так: перемешивали одну последовательность пары, сохраняя состав букв, выравнивали заново и смотрели, на сколько стандартных отклонений настоящий счёт выше среднего по перемешанным. Перемешиваний было десять — больше не позволяло машинное время, — а значимым считали отклонение больше трёх.

Сцена делает то же, только перемешиваний сотни. У α- и β-цепей настоящий счёт 286, а у перемешанных — 29,5 ± 4,9 на 500 попыток: настоящий выше на 52 стандартных отклонения. У миоглобина кашалота и β-цепи человека — пары, которую сравнивали Нидлман и Вунш, — 102 против 29,8 ± 5,0. Лизоцим из яичного белка и α-лактальбумин из молока делают разную работу, но счёт 226 против 27,0 показывает общее происхождение. Рибонуклеаза и лизоцим — пара, которую Нидлман и Вунш взяли как заведомо неродственную, — набирают 31, а 128 перемешиваний из 500 набрали не меньше.

Леггемоглобин сои против миоглобина человека — 26, меньше среднего у перемешанных. По буквам эти белки неотличимы от чужих, хотя их структуры похожи. Тот же леггемоглобин против β-цепи гемоглобина набирает 63, и это уже за пределами случая: родство в семействе видно через одних его членов и не видно через других. Буркхард Рост в 1999 году сравнил больше миллиона пар белков с известной структурой и описал сумеречную зону — 20–35 % одинаковых букв, где сигнал размывается: выше 30 % девять пар из десяти оказывались родственными по структуре, ниже 25 % — меньше одной из десяти. У леггемоглобина и миоглобина в лучшем куске совпадает 20 %.

Перемешивать каждую пару сотни раз дорого, и для локальных выравниваний есть формула. Самуэль Карлин и Стивен Альтшуль в 1990 году показали, что у случайных последовательностей длин m и n, выровненных локально без пропусков, лучший счёт подчиняется распределению Гумбеля — распределению максимумов, — а ожидаемое число кусков со счётом не меньше S равно

E = K · m · n · e^−λS.

Это и есть E-value: сколько раз такой счёт ожидается у пары случайных последовательностей тех же длин. Параметры λ и K зависят от матрицы и частот букв. Для выравниваний с пропусками точной формулы нет, и параметры оценивают по случайным выравниваниям; в исходниках BLAST для BLOSUM62 со штрафом 10 + 1 за букву записано λ = 0,243 и K = 0,024. Кривая в сцене построена по этим числам, а не подогнана под столбики, и средний счёт перемешанных пар она предсказывает с ошибкой меньше двух очков. У миоглобина кашалота с β-цепью E ≈ 9 · 10⁻⁹, у леггемоглобина с миоглобином — 0,95: такое совпадение у случайных белков этих длин бывает почти всегда.

Код

Листинг выравнивает α- и β-цепи с двумя ценами пропусков, проверяет пять пар перемешиванием и сравнивает репрессор cI фага λ с репрессорами HK97 и HK022. Таблицы — три, как у Гото; функция score держит только две строки каждой таблицы и нужна для перемешиваний.

Тридцати перемешиваний хватает, чтобы отличить α- и β-цепи от случайной пары, но оценки среднего и разброса у них неточные: сравните с пятью сотнями из сцены. Последний пункт — пример того, что ищет локальное выравнивание и как выглядит его отсутствие: репрессор HK97 совпадает с репрессором λ почти целиком, а у HK022 лучший кусок — 55 столбцов с 12 одинаковыми буквами и счётом 28, обычным для случайной пары.

Сколько это стоит

Одна пара — m·n клеток в каждой из трёх таблиц, время O(m·n). Если нужен только счёт, хватает двух строк каждой таблицы, и память — O(min(m, n)). Для обратного хода листинг хранит таблицы целиком: три таблицы по четыре байта на клетку — это 12 МБ для двух белков по тысяче букв и около 19 ГБ для двух геномов фагов по сорок тысяч.

Само выравнивание тоже можно получить в линейной памяти. Дэн Хиршберг в 1975 году сделал это для наибольшей общей подпоследовательности, Юджин Майерс и Уэбб Миллер в 1988 году — для аффинных пропусков. Таблицу считают с начала до средней строки и с конца до неё же, находят столбец, через который проходит лучший путь, — это та же сумма двух половин, которой светятся клетки в сцене с листами, — и решают две половины задачи отдельно. Времени уходит примерно вдвое больше.

Для поиска по базе это всё равно медленно: белок против базы из сотен миллионов букв — сотни миллиардов клеток. BLAST, который в 1990 году описали Стивен Альтшуль, Уоррен Гиш, Уэбб Миллер, Юджин Майерс — тот же, что у diff, — и Дэвид Липман, сначала ищет короткие слова с высоким счётом, общие для запроса и последовательности из базы, расширяет их и только вокруг удачных затравок считает таблицу; значимость находок он оценивает по формуле Карлина — Альтшуля. Точный Смит — Уотерман тоже ускоряют: Майкл Фаррар в 2007 году разложил столбцы таблицы по векторным регистрам процессора полосами и получил в 2–8 раз больше, чем у прежних векторных реализаций, — больше трёх миллиардов клеток в секунду на двухъядерном процессоре 2,0 ГГц. Сцена с мозаикой считает таблицы на обычном JavaScript; сколько это заняло у вас, видно под ней.

Варианты

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

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

Множественное выравнивание. Для k последовательностей длины n таблица k-мерная, в ней n^k клеток, и точно её считают только для нескольких коротких строк. На практике выравнивают постепенно: сначала самые близкие пары, потом к ним добавляют следующие по дереву сходства.

Выравнивание структур. Когда буквы уже ничего не говорят, как у леггемоглобина с миоглобином, сравнивают пространственные структуры — координаты атомов цепи.

История

Таблицу для сравнения строк находили несколько раз независимо, в распознавании речи, биологии и исправлении опечаток, — об этом рассказывает глава курса. Статья Сола Нидлмана и Кристиана Вунша поступила в журнал Journal of Molecular Biology 21 июля 1969 года и вышла в 1970-м. Программы на Фортране для машины CDC 3400 сравнивали миоглобин кита с β-цепью гемоглобина человека и, как пару, от которой родства не ждали, рибонуклеазу быка с лизоцимом куриного яйца. Очки за пару аминокислот авторы брали из генетического кода: сколько букв могут совпадать у их кодонов, а за пропуск вычитали штраф. Все семь вариантов очков показали значимое сходство миоглобина и β-цепи и ни один — у рибонуклеазы и лизоцима. Любопытно замечание в тексте: миоглобин человека против β-цепи человека с самыми простыми очками значимого сходства не дал. С BLOSUM62 и аффинными пропусками эта пара набирает 111 против 29,5 ± 5,2 у перемешанных.

Темпл Смит и Майкл Уотерман описали локальное выравнивание в 1981 году в заметке на три страницы «Identification of common molecular subsequences», Осаму Гото годом позже — три таблицы для аффинных пропусков. Матрицы PAM Маргарет Дейхофф вышли в 1978 году, BLOSUM Хеникоффов — в 1992-м. Статистику локальных счетов Карлин и Альтшуль опубликовали в 1990 году, и в том же году вышла статья о BLAST.

Задачи