CC BY
9
Горбунов К.Ю.1, Любецкий В.А.2
аИППИ РАН, г. Москва, к.ф.-м.н., в.н.с., gorbunov@iitp . ги 2ИППИ РАН, г. Москва, д.ф.-м.н., зав.лаб., lyubetsk@iitp . ги
РЕКОНСТРУКЦИЯ ПРЕДКОВЫХ ХРОМОСОМНЫХ СТРУКТУР
КЛЮЧЕВЫЕ СЛОВА
Хромосомная структура, эволюционное дерево, предковый вид, реконструкция вдоль дерева, эффективный алгоритм, булево линейное программирование.
АННОТАЦИЯ
Рассматривается задача реконструкции хромосомных структур в нелистовых вершинах данного филогенетического дерева на основе заданных структур в листьях. Структура может иметь паралоги - гены с одинаковыми именами. Рассматриваемая задача сводится к задаче булева линейного программирования, решаемой эффективным алгоритмом.
Определения и постановка задачи. Будем рассматривать реконструкцию хромосомных структур на основе модели хромосомной структуры как произвольного набора цепей и циклов, состоящих из ориентированных рёбер: генов, которым приписаны имена - натуральные числа i. На рис.1 изображены две структуры: а и Ь.
Структура а
Структура Ь
10
-о-
Ф Ф Ф
-о-
7 -о-
-о-
9
-о-
7 -о-
10
-о-
3
-о-
Рис.1. Две структуры а и Ь, каждая с двумя паралогами - генами, помеченными номером 3
Паралоги - гены из одной структуры, имеющие одинаковые номера /: для их идентификации применяется нумерация вида /.у, в которой число i назовём первым номером гена, а у - вторым номером (само выражение /у будем называть полным номером; полные номера внутри одной хромосомной структуры не повторяются). В [1] мы рассматривали задачу реконструкции в предположении весьма существенного ограничения на модель - структуры должны иметь постоянный генный состав и не иметь паралогов. Здесь рассмотрим общий случай. Будем использовать понятие брейкпоинтового расстояния между двумя структурами. Для структур без паралогов (где есть лишь первые номера) оно определяется как число пар соответствующих краёв генов, соседних (или, как говорится в [1], склеенных) в одной структуре и не склеенных или отсутствующих в другой, сложенное с числом генов, присутствующих в одной структуре и отсутствующих в другой (под генами в разных структурах здесь понимаются гены с одинаковыми номерами).
Для структур с паралогами расстояние определяется как минимум расстояний между структурами, которые получаются установлением частичной биекции между паралогами каждого гена, представленного в обеих исходных структурах. Точнее, генов с первым номером / может быть разное число в структурах а и Ь, пусть это - соответственно множества Ра(/) и Рь(/); обозначим /
биекцию между частью множества Ра(/) и частью множества Рь(/); определяются нумерации /./ всех генов в Ра(Г) и Рь(/) вторыми номерами, относительно которых / - тождественная функция; гены, не попавшие в области определения и значений //, имеют все разные вторые номера. Такая нумерация имеет естественную интерпретацию: если у=//(х), то у «наследуется» от х; прочие гены из Ра(/) и из РьМ - самостоятельные, все различные между собой, поэтому они имеют различные значения j. Гены, не принадлежащие области определения или области значений /, можно считать потерянными и соответственно возникшими (на ребре, соединяющем а с Ь в филогенетическом дереве). Теперь расстояние между а и Ь определяется как минимум расстояний между а' и Ь', которые получаются из а и Ь при всевозможных нумерациях вторыми номерами всех генов в каждом множестве генов с одним и тем же первым номером (такую нумерацию будем называть добавлением вторых номеров). Например, для структур на рис.1 минимальное брейкпоинтовое расстояние достигается при нумерации в каждой структуре одного из паралогов гена 3 как 3.1 (неважно какого), а другого как 3.2. Оно равно 18.
Итак, дано корневое (возможно, небинарное) дерево, и в каждом листе своя хромосомная структура. Требуется расставить хромосомные структуры по нелистовым вершинам, чтобы сумма по всем рёбрам расстояний между структурами на концах ребра была минимальна. Покажем что это эквивалентно более «простому» требованию: расставить хромосомные структуры с нумерацией генов полными номерами по нелистовым вершинам и добавить вторые номера генам в листьях так, чтобы сумма по всем рёбрам расстояний между структурами на концах ребра (уже фактически без паралогов) была минимальна. Пусть найдена расстановка, минимизирующая описанную функцию F при первом требовании. Опишем расстановку, дающую не большее её значение для второго требования. Рассмотрим биекции / между генами на концах прикорневых рёбер, соответствующие определению функции К Добавим генам в корне дерева вторые номера произвольным образом. Добавим вторые номера в некорневом конце каждого прикорневого ребра в соответствии с Таким образом, идя от корня, мы добавим вторые номера генам во всех вершинах дерева, включая листья. Поскольку сами биекции / не меняются, получаем расстановку для второго требования с тем же значением К. Обратное соответствие (от второго требования к первому) получается «стиранием» вторых номеров у всех генов. Таким образом, далее рассматриваем задачу в последней формулировке.
Сведение задачи реконструкции к задаче булевого линейного программирования. Легко видеть, что в искомой расстановке имеются лишь гены, полные номера которых присутствуют хотя бы в одном листе. Кроме того, нам удобно считать, что каждая (тождественная) биекция //, определяемая искомой расстановкой, определена на объединении М всех присутствующих в листьях номеров /./, а отсутствующие в структуре гены специальным образом помечены. Очевидно, это определение эквивалентно исходному. «Началом» ребра в дереве считается та вершина, которая ближе к корню; аналогично определяется «конец» ребра.
Добавим произвольным образом генам в данных листовых структурах вторые номера; полученную нумерацию обозначим как I. Для перехода к задаче линейного булевого программирования введём переменные, линейные ограничения и линейную функцию, которую нужно минимизировать.
Введём переменные, указывающие соответствие между нумерацией I и искомой нумерацией генов в листьях. Для каждого листа е дерева и каждой упорядоченной пары к./, к/ элементов из М введём переменную Zk/je (первый индекс будем для краткости опускать); эта переменная равна 1, если ген к/ в е имеет в искомой нумерации номер к./, иначе zк/e=0. Соответствие между генами с первым номером к, определённое переменными z/je, должно быть биективным, что выражается линейными равенствами: для фиксированного / и е сумма по / значений z/je равна 1, и аналогично для фиксированных/ и е.
Введём переменные, определяющие хромосомные структуры в вершинах.
Для каждой вершины V и каждой неупорядоченной пары к, 1 различных краёв генов введём переменную х^; эта переменная равна 1 если эти края склеены в вершине V, иначе она равна 0. Каждый край может быть склеен не более чем с одним краем, что выражается линейными неравенствами: при фиксированных к и V сумма по 1 значений хш не превышает 1, и аналогично для фиксированных 1 и V. В листьях значения переменных хи заданы.
Для каждой вершины V и каждого гена к введём переменную у^; эта переменная равна 1 если ген к отсутствует в V; иначе она равна 0. Края отсутствующих генов ни с чем не склеены, что выражается линейными неравенствами х^<1-ук» где / или / - край гена к. Для каждого гена к и
каждого его края / вместо этих неравенств можно использовать одно неравенство ^ 1
, где/ пробегает все края генов в V. В листьях значения переменныхук заданы.
Введём переменные, вычисляющие брейкпоинтовое расстояние между структурами, приписанными началу и концу ребра дерева. Пусть края к и 1 генов лежат в конце ребра е, а края т, п - в его начале. Назовём четвёрку к, 1, т, п краёв генов брейкпойнтовой, если пара возможных склеек - склейка краёв к и 1 и склейка краёв тип, - учитывается при определении брейкпойнтового расстояния на нелистовом е или может учитываться на листовом е после какого-нибудь (возможно, отличного от I) добавления вторых номеров генам в листьях (эта четвёрка -упорядоченная пара двух неупорядоченных пар: (к,Г) и (т,п)).
Для каждого ребра е и каждой брейкпойнтовой четвёрки к, 1, т, п краёв введём переменную Sklmne; эта переменная равна 1, если эти края принадлежат соответствующим генам (напомним: на листовом е соответствие определяется переменными типа z) и склеены в начале ребра и не склеены (или хотя бы один ген отсутствует) в конце ребра или наоборот. Это условие для Sklmne выражается линейными неравенствами; приведём их ниже для случая, когда е листовое, все четыре края - начала генов (или все - концы) и все четыре гена имеют одинаковые первые номера (другие случаи аналогичны): если е=(ы^), гены /1 и / имеют, соответственно, края к и т, а гены /2 и/2 - края 1 и п, то
Sk1mne — Хk1v-Хmnu-(1-Z/yle)-(1-Z/2j2e), Sk1mne — Хmnu-Хk1v-(1-Z/yle)-(1-Z/2j2e),
Sk1mne — Хk1v—Хтпи-(1-Z/1/'2e) — (1—Z/2/1e), Sk1mne — Хтпц—Хк^—(1 — Z/1/2e) — (1 — Z/2/1e).
Если правая часть неравенства равна 0, переменная Sklmne принимает значение 0 за счёт того, что она входит в качестве слагаемого в минимизируемую функцию.
Для каждого ребра е=(ы^) и каждых двух генов / и /, лежащих, соответственно, в конце и в начале ребра е, введём переменную / эта переменная равна 1, если эти гены соответствуют друг другу и один из них присутствует в структуре, а другой отсутствует. Это условие выражается (для листового е] линейными неравенствами > у,,.,-у/„-(1-2,/„] и% >у,1-у1[-[1-г1,,]. з)
о о о о о о
Ь)
1 1 1
% % % % % %
Рис.2. Пример реконструкции для искусственных данных. а) В каждом листе дерева дана хромосомная структура с тремя паралогами одного и того же гена. Ь) Реконструированные структуры, построенные алгоритмом в нелистовых вершинах, и добавленные вторые номера генов в листьях
Минимизируемая функция - сумма всех переменных Sklmne и / линейные ограничения описаны выше. Решение описанной задачи линейного булева программирования даёт искомую расстановку хромосомных структур по нелистовым вершинам дерева и вторые номера генам в
листьях.
Замечание 1. Описанное сведение задачи реконструкции хромосомных структур к задаче булевого линейного программирования легко обобщается на случай, когда стоимости брейкпоинтовых операций (переход от двух склеенных краёв к расклеенным или наоборот, и переход от присутствия гена к отсутствию или наоборот) различны. Для этого каждую переменную Skimne следует заменить на две переменных slkimne и s2kimne. Для первой из них действуют линейные ограничения вида slkimne ^ Xkiv-Xmnu-—', для второй - ограничения вида s2kimne ^ Хтпи-Хш-.... Аналогично, переменную Sje следует заменить на sljje и s2ye. Теперь можно в минимизируемой функции умножать новые переменные на коэффициенты, отражающие веса событий.
Замечание 2. Задачу вычисления брейкпойнтового расстояния между двумя данными хромосомными структурами с паралогами можно рассматривать как частный случай задачи реконструкции. Здесь "дерево" состоит из одного ребра с данными структурами на концах. Очевидно, остаются лишь переменные типа z и s.
На рис.2Ь приведено решение задачи реконструкции хромосомных структур для искусственных данных, показанных на рис.2 а. Все гены в данных листовых структурах имеют один и тот же первый номер 1 (не показан), в реконструированных структурах показаны лишь вторые номера.
DCJ conversion of genome "Genome A" to "Genome В"'
1 2 -3.x 4 S 6 -7 в 9 lO 3.2
Step (J; ^---* ^ ^ > ^ *-•--7» •-3»
c±>
S 4 XI -J,l i -3-3 -7 9 * Ю
Step*/- » 4 -« * * >
Step 2
X
6 -3.1 14 -7 В 9 1П
)f * x—* » »>, Ж—«..*> » *—>
D; &
6 -3.2 1 Z -3.1 S 4 -7 S 9 ID
Step ^ Ж t ** ^ t ** —*——^ *—
Step i
i -3*2 5 4 -7 i S 1Д -2-1 3,1
4 .-—ж—• > * + » > .-< *-Ж
= CZZZ)
6- -S.2 ! 4
Slep 5:
Step 6:
С
6 -3.1 5 4 -7 10 Sal -2 -1 8 9
-ж- >—-л- ж f 4 * .j» %
i -3.2 i 4 -7 ID 3.1 -2 S 9
Step 7:^ X * » « » > * + ' >
-7 ID 3.1 8 9
-Ы—> >
6 -3.2 E 4
StepS
^ Ж W JT j»^
■7 ID 3.1 S 9 1 2
* * *
Step 9:
S -3-2 5 4 -7 Jit i.l
О о сЪ о о
Рис.3. Кратчайшая последовательность операций, переводящая структуру а в структуру Ь (рис.1) для случая, когда верхний ген 3 в структуре а соответствует нижнему гену 3 в структуре Ь.
Обобщение понятия расстояния. Пусть, для простоты, даны две структуры с одинаковым набором генов. Брейкпоинтовое расстояние между ними можно рассматривать как минимальное число операций расклейки двух склеенных краёв генов или склейки двух свободных краёв, необходимых для перевода одной структуры в другую. Добавим к этим операциям следующие операции над хромосомной структурой: двойная переклейка - расклейка двух склеек краёв генов и переклейка четырёх краёв, приводящая к новой структуре и полуторная переклейка - расклейка двух склеенных краёв и склеивание одного края с каким-то несклеенным краем, второй край остаётся свободным. Получим обобщение брейкпойнтового расстояния, называемое биологическим расстоянием. Для случая двух хромосомных структур с одинаковым набором генов (без паралогов) реализация алгоритма вычисления биологического расстояния описана в [3]. Пример с двумя структурами на рис.1 показывает, что биологическое расстояние "чувствительнее" брейкпоинтового: теперь два возможных соответствия между паралогами гена 3 дают различные значения расстояния. Если верхний ген в структуре а соответствует нижнему гену в структуре Ь, перевести а в Ь можно самое меньшее за девять операций - рис.3; при альтернативном соответствии достаточно семи операций - рис.4 (рисунки получены программой, описанной в [3]). о? депснп* "Сепмпе А'" 1о "О^попие В".
1 2 -3.1 4 I в
■7 е 9 ifl s. г
3bp0:f » * ♦ > » ^ « * » » —* С*}
5 4 12 -9.1 Ь -7 S 9 Ш 3.2
ЗИМ Г > * ♦ ^ * ' > *
6 12 - J.l t -I
S.epZi^Jt * * ■ > ■< » > >
•7 г $ 1Д з. 2
CD
6 -3.1 t * -г I i ло з.г i i
Slep 3:
С
ж—* * + »>
■
CD CDD
■3.1 f
Slip A:
С
-M-
» -
-Г 1Д
3
i.:
> ::*
CDD
6 -3.1 t 4 -7 ID S.J
S|fp 5. ^ Ж * * * * * *
e 9
г ^i
Slep
С
С -3.1 E 4
-Ж—
-7 ID S,i
> > » »-*
6 4
cb О
-3.1 f
Step 7:
С
-Ш-
5
-Г 10 > » »—>
ООО
Рмс.4. Кратчайшая последовательность операций, переводящая структуру а в структуру Ь (рис.1) для случая, когда верхний ген 3 в структуре а соответствует верхнему гену 3 в структуре Ь
Для биологического расстояния задачи его вычисления (в том числе для структур с различным набором генов и с паралогами) и соответствущей реконструкции хромосомных структур рассматриваются в [1, 2], а также в нескольких готовящихся к публикации работах
авторов. Реализации соответствующих алгоритмов доступны на сайте http://lab6.iitp.ru/ru/chromo/.
Работа выполнена за счёт гранта Российского научного фонда (проект № 14-50-00150).
Литература
1. Горбунов К.Ю., Гершгорин Р.А., Любецкий В.А. Перестройка и реконструкция хромосомных структур // Молекулярная биология, 2015, том 49, № 3, стр. 372-383.
2. Gershgorin R.A., Gorbunov K.Yu., Seliverstov A.V., Lyubetsky V.A. Evolution of Chromosome Structures // "Information Technology and Systems 2015" An IITP RAS Interdisciplinary Conference & School (ITaS'15), Sochi, Russia, Sep 7-11 2015, pp. 105-120.
3. Hilker R., Sickinger C., Pedersen C., and Stoye J. UniMoG - a unifying framework for genomic distance calculation and sorting based on DCJ // Bioinformatics, 2012, vol. 28, pp. 2509-2511.
