CC BY
16
УДК 519.65
А. С. Подкопаев1, М.Е. Карасиков1'2, Ю.В. Максимов1,3
1 Московский физико-технический институт (государственный университет) 2 Сколковский институт науки и технологий 3Институт проблем передачи информации РАН
Прогнозирование структур белков методами полуопределенного программирования
Задача заключается в предсказании упаковки белковых молекул в мультимерный комплекс в приближении жестких тел. Для решения поставленной задачи предлагается использовать методы выпуклой оптимизации, например, полуопределенные релаксации. Недостатком большинства существующих алгоритмов (жадных алгоритмов и других) является их вычислительная сложность. В данной работе предлагаются алгоритмы меньшей вычислительной сложности. Основным результатом является оценка их качества, сравнение с алгоритмами, использовавшимися ранее.
Ключевые слова: структуры белков, выпуклая релаксация, вычислительная сложность, хордальные графы, максимальные клики, точность прогнозирования.
A. S. Podkopaev1, M. E. Karasikov1'2, Y. V. Maximov1,3
1 Moscow Institute of Physics and Technology (State University) 2Skolkovo Institute of Science and Technology 3Institute of Infomation Transmission Problems RAS
Protein packing by semidefinite programming
The problem is to predict protein packing into the multidimensional complex in the approximation of rigid bodies. This problem can be solved using methods for convex optimization, for example, semidefinite relaxation. The disadvantage of the majority of existing algorithms (greedy algorithms, etc.) is their computational complexity. This study offers the algorithms with lesser computational complexity for specific input data. The main result is the assessment of their properties and performance in comparison with the algorithms used before.
Key words: protein structures, convex relaxation, computational complexity, chordal graphs, maximal cliques, forecast accuracy.
1. Введение
В основе каждого белка лежит полипептидная цепь, которая организована в трехмерную структуру. Различают четыре уровня пространственной организации белка: первичная, вторичная, третичная и четвертичная структуры белковых молекул. Под четвертичной структурой понимают макромолекулярное образование из нескольких полипептидных цепей в составе единого мультимерного белкового комплекса. В работе решается задача предсказания упаковки боковых цепей белковой структуры при помощи выпуклой оптимизации в предположении известного положения основной цепи.
В работах по прогнозированию структур белков [1, 2] основное внимание уделяется взаимодействиям между парами белков. В качестве результата используемые методы выделяют наиболее правдоподобных кандитатов, которые впоследствии при помощи целевой функции ранжироваются. На данный момент среди наиболее популярных подходов к решению данной задачи следует отметить метод прогнозирования структуры в приближении жестких тел с последующим применением алгоритма полярной корреляции Фурье, описанный в статье [1], а также предложенный в работе [2] подход, основанный на методе Монте-Карло.
Особое внимание в литературе уделяется решению задачи предсказания упаковки боковых цепей. В статьях [3,4] описаны одни из лучших по качеству и скорости методов решения задачи упаковки боковых цепей. В ранних работах М.Е. Карасикова и Ю.В. Максимова также рассмотрена задача предсказания упаковки боковых цепей белковой структуры при помощи выпуклой оптимизации в предположении известного положения основной цепи. Важно отметить, что задача прогнозирования структур является ЖР-трудной задачей квадратичной оптимизации. Проблема упаковки белковых молекул в комплекс связана с исследованием большого числа всевозможных вариантов структур белка (большой размер базового пространства поиска).
Перейдем к формальной постановке задачи: предположим, что в цепи имеется п белков, причем каждый г-й белок может располагаться в одной из р позиций: хг = [х\,х2, ...,хгр], где хгк € {0,1}. Для типичных макроструктур предполагается п ~ 10, р ~ 100. Тогда одному комплексу будет соответствовать вектор х = [х1,... , жга]т. Любой паре позиций хгр и Хт можно приписать функцию энергии д, которая обусловлена силовым взаимодействием. Каждой позиции белка также соответствует собственное значение энергии Ьо (эти значения получены из экспериментальных данных). Обозначим через N число всевозможных позиций белков, составляющих один комплекс.
Тогда задача оптимальной упаковки может быть записана в виде задачи минимизации:
хт0ох + Ь^х ^ шт
(1)
при условии, что Ах = 1П, где бинарная матрица А размера Nхп имеет единичное значение
г— 1 г
элемента А^, если и только если ] : ^ ръ < ] ^ ^ Pt. А матрица фо имеет вид
4=1 4=1
Яо =
[0] [Ч12(Р1,Р2)\ [Ч1п(Р1,Рп)\
[Ч21(Р2,Р1)} [0] ••• [д2п(р2,Рп)]
[Яп1(Рп,Р1)] [Яп2(Рп,Р2)} ■■■ [0]
В общем случае, матрица фо не является положительно полуопределенной. Одна из формулировок задачи записывается как невыпуклая оптимизация. Изначально задача задана на булевом кубе: В^ = {х : XI € {0,1}, 1 ^ г ^ п}. Замечая, что хтЯ^х = {Яо,ххт}м,м, переходим к задаче выпуклой оптимизации. С использованием различных методов релаксации можно попытаться получить приближенное решение исходной задачи с высокой точностью за полиномиальное время.
Решение задачи определяется размерностями и свойствами матрицы фо: в случае малых размерностей достаточно воспользоваться жадными алгоритмами, а в случае больших размерностей предлагается воспользоваться методами выпуклой оптимизации: например, полуопределенной релаксацией, релаксацией Лагранжа. Целью данной статьи является построение алгоритма меньшей вычислительной сложности с использованием предложенных методов [5-8].
2. Постановка задачи
В силу введенных в формальной постановке задачи предположений получаем, что вектор х имеет блочную структуру, то есть он может быть поделен на малое количество блоков п ~ 10, а на каждый блок размера р ~ 100 накладывается ограничение вида ||жг||те = 1. Вектор х представляет собой набор вероятностей того, что составные белки занимают определенное положение в пространстве, а матрицы Qi и векторы Ьi предоставляют данные о ранее известных белковых структурах. Пусть задана выборка из N матриц энергий: {(^1,... }. Квадратные матрицы энергий = {д^}ij=1 т имеют размеры порядка
т = пр ~ 1000 и являются сильно разреженными (доля нулевых элементов составляет порядка 99%).
Целью работы является исследование задачи квадратичной оптимизации, возникающей при решении проблемы оптимальной упаковки белков, и повышение эффективности методов выпуклой оптимизации за счет ее структуры.
3. Используемые алгоритмы
Для снижения вычислительной сложности исходной невыпуклой задачи квадратичной булевой оптимизации, являющейся ЖР-трудной, осуществляется переход к графовой структуре данных (построение описано ниже), а затем задача декомпозируется древесным образом. В процессе используются различные эвристические методы - например, работа с хордальным расширением полученного графа. В результате это приводит к решению задачи за меньшее время в сравнении с обычной полуопределенной релаксацией. Перейдем непосредственно к описанию используемых методов.
3.1. Полуопределенная релаксация
При использовании методов выпуклой оптимизации для решения задачи 1 невыпуклые ограничения заменяются на более слабые, но выпуклые ограничения. Тогда мы переходим к выпуклой задаче, методы решения которой хорошо известны. В результате удается получить хорошую оценку в задачах поиска минимума значения целевой функции. Основным методом выпуклой релаксации, рассматриваемым в этой работе, является полуопределенная релаксация.
Полуопределенную релаксацию можно представить в виде
Яц Ху ^ тт ,
г,] + ' - "
где = : X = Xт ^ 0} и, кроме того, X ^ ххт, Хгг = хг, г = 1,..., N и Ах = 1п.
Используя лемму Шура, перепишем задачу в виде эквивалентной:
ЕЯгэХгп ^ Ш1П
г,3
при дополнительных условиях [X х;хт 1] € 5++1, Хгг = хг, г = 1,... , Ах = 1Г Метод реализуется на МаЛаЪ при использовании пакета сух.
3.2. Хордальная структура данных
Для снижения вычислительной сложности алгоритма воспользуемся результатами из теории графов. Заметим, что каждой квадратной матрице ^о порядка п в соответствие можно поставить некоторый граф С с п вершинами, ребра которого строятся по следующим правилам: наличию ненулевого элемента {г^} в матрице ^о соответствует ребро, соединяющее вершины г и Если же элемент {I,]} в матрице ^о равен 0, ребро между элементами г и ] не проводится.
Определение 1. Граф называется хордальным, если каждый из его циклов, имеющий четыре и более дуг, имеет хорду, которая является ребром, соединяющим две вершины, не смежные в цикле.
Поставим задачу поиска максимальных клик в графе С. Поскольку про структуру таблицы смежности исходного графа ничего не известно, то можно считать, что граф С является произвольным. В общем случае задача поиска максимальных клик является ИР-полной. Однако для исходных графов можно построить хордовое расширение и, воспользовавшись свойствами этих графов, получить решение этой задачи за полиномиальное время. В нашей работе для построения расширения использовался элиминационный алгоритм [12], состоящий из трех основных шагов:
1) в графе G = (V., Е) на п вершинах фиксируем некоторый порядок просмотра вершин а. Пусть Gi = G, Е' = 0;
2) добавим ребра к Gi так, чтобы все соседи вершины г стали попарно смежными. Добавленные ребра включим во множество Е' и удалим из графа Gi вершину г получив граф Gi+i;
3) вычислим хордовое расширение G' = (V; Е U Е') графа G.
Для полученного хордового расширения воспользуемся теоремой [9].
Теорема 1. Пусть матрице X соответствует граф G, Ki,...,Kd - максимальные клики графа G. Тогда следующие условия эквивалентны:
1) Можно дополнить X до положительно полуопределенной матрицы;
2) XKitKi ^ 0,г = l,...,d, - подматрица X с номерами строк и стобцов, принадлежащими Ki.
Тогда при решении задачи 1 для тех матриц энергий, для которых размер максимальной клики у соответствующего ей графа существенно меньше размера этого графа, заменим исходную полуопределенную релаксацию на ее модификацию: поиск минимума не по всем X, удовлетворяющим необходимым условиям, а по отдельным подматрицам существенно меньших размеров, соответствующих вершинам, входящим в состав максимальных клик.
Приведем теоретические факты и алгоритмы, используемые для вычисления максимальных клик в хордальных графах за полиномиальное время [11, 12].
Следующий алгоритм (MSC, maximum cardinality search) вычисляет совершенный порядок исключения вершин графа G = (V, Е) на п вершинах:
1) для всех вершин v £ V положим £[v] =0;
2) для всех г, l ^ г ^ п, выберем непронумерованную вершину v с наибольшим значением l[v] и положим &(v) = г. Для всех смежных с v вершин ш в графе G увеличим £[ш] на единицу.
Порядок а на графе G, полученный в результате работы указанного алгоритма, есть совершенный порядок исключения вершин графа G, если граф хордальный.
Теорема 2. Неориентированный граф является хордальным тогда и только тогда, когда алгоритм MSC вычисляет (G,a), где а является совершенным порядком исключения для графа G.
Теорема 3. Каждая максимальная клика хордального графа G = (V, Е) имеет вид {f} U X(v), где X(v) = {w £ adj(v)la(v) < a(w)}, adj(v) - множество вершин, смежных с v в графе G.
Теорема 4. Хордальный граф имеет не более, чем п максимальных клик.
Алгоритм вычисления максимальных клик для хордального графа G = (V, Е) на п вершинах и совершенного порядка исключения а имеет вид:
1) для всех вершин v £ V инициализируем параметр размера s[v] =0;
2) для всех г, l ^ г ^ п, положим v = v-i(i) и
• если степень вершины d(v) = 0, то возвращаем максимальную клику {f},
• вычисляем множество X = {w : (v,w) £ Е и a(v) < a(w)},
• если множество X непусто и s[v] < |Х |, то возвращаем максимальную клику {f} U X. Иначе обновим оценку расстояния s[m[v]] = max{ s[m[v]], |Х | — 1} для вершины m[v] = ст-1(тт{ст(-ш)|-ш £ X});
3) в качестве ответа алгоритм возвращает набор всех максимальных клик.
В результате от ЖР-полной задачи нахождения максимальных клик переходим к задаче для хордальных графов, решаемой за полиномиальное время.
4. Вычислительный эксперимент
Матрицы энергий (379 штук) из PDB (Protein Database Bank) были разбиты на группы по количеству вершин в соответствующих графах: до 500 (малых размеров), 500-1000 (средних размеров), больше 1000 (больших размеров). Для графов из каждой категории (73 - в первой категории, 145 - во второй и 161 - в третьей соответственно) были построены хордальные расширения и вычислены размеры максимальных клик. Приведем результаты в табл. 1.
Т а б л и ц а 1
Соотношение размеров графов и их максимальных клик в исследуемой
выборке
размер клики <12 12-16 16-20 >20
менее 500 вершин 4 29 32 8
500-1000 вершин 1 41 87 16
более 1000 вершин 0 16 120 25
Для исходной матрицы энергий из выборки приведем графики зависимости, на которых по оси абсцисс откладываются фиксированные размеры, а по оси ординат отложим число графов, размер максимальных клик которых не превосходит это фиксированное число.
Рис. 1. Размеры максимальных клик
Учитывая приведенные выше обоснования, делаем вывод, что на исходной выборке следует попробовать применить модификацию исходной полуопределенной релаксации. На исходной выборке из PDB можно сравнить скорость и точность работы полуопределенной релаксации с/без использования хордовой структуры данных, и сопоставить их методу Монте-Карло.
Вычислив максимальные клики для графов, соответствующих матрицам энергий, можно перейти к работе с методами релаксации. Была изучена зависимость времени работы SDP, SDP с учетом хордальной структуры данных и Монте-Карло от размеров исходных матриц.
Сравнивались средние времена работы исходного SDP и модифицированной задачи на всей выборке из 379 матриц. Было получено, что в среднем на любой матрице из выборки модифицированная задача решается в 3.24 раза быстрее.
Рис. 2. Вычислительная сложность для матриц Рис. 3. Вычислительная сложность для матриц всех размерностей малых размерностей
Также можно посмотреть, в скольких случаях модифицированный SDP дает значения целевой функции быстрее метода Монте-Карло. На нашей выборке первый метод давал результаты быстрее на 11% исходной выборки. В среднем метод работал быстрее на квадратных матрицах размера 200-300. Если задать допустимую точность вычисления значения целевой функции 6 = 5%, то расчет показывает, что в сравнении со стохастическим методом Монте-Карло модифицированный SDP дает более точное значение целевой функции на 2.84% матриц, на которых он работает быстрее.
Таблица2
Задачи с наибольшим выигрышем в скорости модифицированного SDP
алгоритма
id SDP Модифицированный SDP Монте-Карло Размер
Значение Время Значение Время Значение Время
1ado 13452 4.91 13624 0.85 12439 1.9 227
1b9w 15995 11.51 16201 2.48 15416 3.36 336
1cei 15830 10.16 15088 2.06 14717 3.66 335
1f94 13937 5.69 14006 1.31 13131 2.51 244
1gvp 16686 9.22 15184 1.91 15696 2.65 322
1oai 12595 6.62 11859 1.28 11599 2.56 258
1uln 14269 7.44 14079 1.61 13699 2.56 281
1ulr 16295 7.42 14294 1.66 15004 2.64 294
1wm3 14278 7.22 13305 1.57 12751 2.5 292
2g30 11522 6.96 11578 1.36 11508 2.33 275
2gqv 11427 4.39 11609 0.68 10371 1.76 190
1yu5 12412 6.19 11704 1.31 12255 2.32 258
1yzm 8757.2 4.54 8197.3 0.73 7995.7 1.7 190
Приведена таблица с указанием времени работы трех методов на части выборки, для которой соответствующие квадратные матрицы энергий имеют размеры до 500 строк (столбцов), и полученным значением минимизируемого функционала. Под id понимается метка данной белковой структуры, а под размером - размер соответствующей матрицы энергий . В табл. 2 приводятся результаты, в которых хордальная модификация показала наиболее высокие результаты по сравнению с другими методами.
Вместе с тем следует отметить, что в ряде случаев качество методов полуопределенной релаксации существенно выше, чем качество методов Монте-Карло. Пять задач с наи-
большим выигрышем в качестве решения методами полуопределенного программирования приведены в табл. 3.
Таблица3
Задачи с наибольшим выигрышем в качестве полуопределенной релаксации
id SDP Модифицированный SDP Монте-Карло Размер
Значение Время Значение Время Значение Время
2hpl 18169 13.52 17641 3.62 19199 3.95 398
2o0q 20229 14.06 17866 3.35 18621 3.41 390
2i49 16099 10.83 14929 2.73 15683 3.23 350
2fl4 18492 11.91 17055 2.97 18310 3.41 368
1ulr 16295 7.42 14294 1.66 15004 2.64 294
5. Заключение
В данной работе рассмотрены алгоритмы прогнозирования четвертичных структур белков, построенных на основе предположения о минимизации суммарной энергии. Для невыпуклой задачи оптимизации были рассмотрены методы выпуклой релаксации как без учета хордальной структуры данных, так и с учетом этой структуры. Был проведен вычислительный эксперимент на основе данных из PDB. Проведено сравнение результатов работы этих двух алгоритмов и метода Монте-Карло. Показано, что в задачах малой размерности использованние выпуклой релаксации с хордальной структурой позволяет повысить скорость решения. В целом точность решения методами полуопределенной релаксации и методом Монте-Карло примерно равна, что по всей видимости связано с особенностями порождения данных задачи.
Использованная в работе идея из [10] заключалась в разбиении исходной задачи полуопределенной релаксации на подзадачи определенным способом и последующим решением подзадач также методами полуопределенной релаксации. Дальшейшее повышение точности результатов возможно за счет явного представления оптимума в небольших кликах в виде последовательности максимумов, решаемой явно и встраивании данной конструкции в двойственную задачу.
В работе было показано, что в ряде задач структурной биологии данные обладают необходимыми свойствами для проведения эффективной декомпозиции исходной задачи на задачи меньшей вычислительной сложности (например, использованная в статье хор-дальная структура данных). В ряде случаев удается получить существенный выигрыш в скорости и качестве решения поставленной задачи.
Авторы выражают благодарность С. В. Грудинину за постановку задачи и предоставление данных. Работа выполнена при финансовой поддержке грантов РФФИ 14-07-31241 мол_а (Ю. В. Максимов) и 15-07-09121_а (М. Е. Карасиков и А. С. Подкопаев).
Литература
1. Moughon G., Wang, Schueler-Furman, Kuhlman B. Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations // Journal of Molecular Biology. 2003. P. 281-299.
2. Topf M., Lasker K., Webb B., Wolfson H., Chiu W., Sali A. Protein structure fitting and refinement guided by cryo-EM density // Structure. 2008. V. 16, N. 2. P. 295-307.
3. Krivov G., Shapovalov M., Dunbrack R. Improved prediction of protein side-chain conformations with SCWRL4 // Proteins. 2009. P. 778-795.
4. Zhichao Miao, Yang Cao, Taijiao Jiang RASP: rapid modeling of protein side chain conformations // Bioinformatics. 2011. P. 3117-3122.
5. Freund R. Introduction to Semidefinite Programming. 2009. P. 1-49.
6. Boyd S., Vandenberghe L. Convex Optimization // Cambridge University Press. 2004. P. 1716.
7. Nesterov Y. Introductory lectures on convex optimization // Springer Science & Business Media. 2004. V. 87.
8. Nesterov Y. Quality of semidefinite relaxation for nonconvex quadratic optimization // Universite catholique de Louvain, Center for Operations Research and Econometrics (CORE). 1997. N 1997019.
9. Grone R., Charles R., Eduardo M.,Wolkowicz H. Positive definite completions of partial Hermitian matrices // Linear algebra and its applications. 1984. V. 58. P. 109-124.
10. De Klerk E. Exploiting special structurein semidefinite programming: A survey of theory and applications // European Journal of Operational Research. 2010. V. 201, N 1. P. 1-10.
11. Rose D., Lueker G., Tarjan R. Algorithmic aspects of vertex elimination on graphs // SIAM Journal on Computing. 1976. V. 5. P. 266-283.
12. Tarjan R. and Yannakakis M. Simple Linear-time Algorithms to Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs // SIAM Journal on Computing. 1984. — V. 13, N 3. P. 566-579.
References
1. Moughon G, Wang, Schueler-Furman, Kuhlman B. Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations. Journal of Molecular Biology. 2003. P. 281-299.
2. Topf M., Lasker K., Webb B., Wolfson H., Chiu W., Sali A. Protein structure fitting and refinement guided by cryo-EM density. Structure. 2008. V. 16, N 2. P. 295-307.
3. Krivov G., Shapovalov M., Dunbrack R. Improved prediction of protein side-chain conformations with SCWRL4. Proteins. 2009. P. 778-795.
4. Zhichao Miao, Yang Cao, Taijiao Jiang RASP: rapid modeling of protein side chain conformations. Bioinformatics. 2011. P. 3117-3122.
5. Freund R. Introduction to Semidefinite Programming. 2009. P. 1-49.
6. Boyd S., Vandenberghe L. Convex Optimization. Cambridge University Press. 2004. P. 1716.
7. Nesterov Y. Introductory lectures on convex optimization. Springer Science & Business Media. 2004. V. 87.
8. Nesterov Y. Quality of semidefinite relaxation for nonconvex quadratic optimization. Universite catholique de Louvain, Center for Operations Research and Econometrics (CORE). 1997. N 1997019.
9. Grone R., Charles R., Eduardo M.,Wolkowicz H. Positive definite completions of partial Hermitian matrices. Linear algebra and its applications. 1984. V. 58. P. 109-124.
10. De Klerk E. Exploiting special structure in semidefinite programming: A survey of theory and applications. European Journal of Operational Research. 2010. V. 201, N 1. P. 1-10.
11. Rose D., Lueker G., Tarjan R. Algorithmic aspects of vertex elimination on graphs. SIAM Journal on Computing. 1976. V. 5. P. 266-283.
12. Tarjan R. and Yannakakis M. Simple Linear-time Algorithms to Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs. SIAM Journal on Computing. 1984. V. 13, N 3. P. 566-579.
Поступила в редакцию 12.11.2015.
