CC BY
6
Полуян С.В., Ершов Н.М.
DOI: 10.33693/2313-223X-2019-6-4-29-43
С.В. Полуян1, Н.М. Ершов2
1 Государственный Университет «Дубна»,
141982, Дубна, Московская область, Российская Федерация
2 Московский государственный университет имени М.В. Ломоносова, 119991, Москва, Российская Федерация
КВАНТИЛЬНОЕ ПРЕОБРАЗОВАНИЕ В ЗАДАЧАХ СТРУКТУРНОЙ БИОИНФОРМАТИКИ
Аннотация. Настоящая работа посвящена особенностям реализации многомерной эмпирической квантильной функции, для которой элементы выборки распределены на узлах регулярной сетки. В работе приведено краткое «рекурсивное» определение структуры многомерной квантильной функции и предложен алгоритм, позволяющий выполнять как дискретное, так и непрерывное многомерное квантильное преобразование. Приводятся результаты численного исследования представленного алгоритма преобразования и показана вычислительная сложность в зависимости от способа хранения выборки. Продемонстрированы результаты применения эволюционных алгоритмов оптимизации с использованием квантильного преобразования в двух задачах структурной биоинформатики: задаче предсказания пространственной структуры пептида по аминокислотной последовательности и задаче локального пептид-белок докинга с известным линейным интерфейсом связывания и линейной структурой пептида.
Ключевые слова: эмпирическая квантильная функция, конечный автомат, глобальная оптимизация
f Л
ССЫЛКА НА СТАТЬЮ: Полуян С.В., Ершов Н.М. Квантильное преобразование в задачах структурной биоинформатики // Computational nanotechnology. 2019. Т. 6. № 4. С. 29-43. DOI: 10.33693/2313-223X-2019-6-4-29-43
V J
DOI: 10.33693/2313-223X-2019-6-4-29-43
S. Poluyan1, N. Ershov2
1 Dubna State University,
Moscow Region, Dubna, 141982, Russian Federation
2 Lomonosov Moscow State University, Moscow, 119991, Russian Federation
QUANTILE TRANSFORM IN STRUCTURAL BIOINFORMATICS PROBLEMS
Abstract. In this paper we study features of the multivariate empirical quantum function implementation for which sample is distributed at the mesh points of the regular grid. We present an algorithm for continuous and discrete quantile transform based on recursive definition of the multivariate quantile function. We perform numerical study of the presented algorithm and demonstrate it computational complexity according to representation of the sample. We present the results of using evolutionary optimization algorithm with quantile transform for solving the problems in structural bioinformatics: protein structure prediction from amino acid sequence and protein-peptide docking with known binding site and linear peptide structure.
Key words: empirical quantile function, finite-state automata, global optimization
f \
CITATION: Poluyan S., Ershov N. The ray method for solving the dynamics of the elastic viscous-plastic shells. Computational nanotechnology. 2019. Vol. 6. No. 4. P. 29-43. DOI: 10.33693/2313-223X-2019-6-4-29-43
i J
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
1. Введение
Целью настоящего исследования является совершенствование ранее предложенного подхода [1] применения стохастических алгоритмов оптимизации в задачах структурной биоинформатики. В основе предложенного подхода лежит использование многомерной эмпирической квантильной функции [2] (или просто квантильного преобразования). Необходимость совершенствования связана с высокой вычислительной сложностью упомянутого алгоритма квантильного преобразования и невозможностью его применения для задач высокой размерности.
Для того чтобы выполнить квантильное преобразование, необходима выборка и алгоритм преобразования. В настоящей работе предложен подход к вычислению квантильной функции, построенный с применением парадигмы разделяй и властвуй, который позволяет эффективно хранить выборку и выполнять преобразование с учетом особенностей хранения элементов выборки. Предлагаемый подход и его реализация являются предметами настоящего исследования.
Необходимо отметить, что основной особенностью предложенного подхода выполнения квантильного преобразования является способ хранения выборки, на основе которого предложены несколько алгоритмов преобразования. В настоящей работе проведено сравнение с ранее предложенным подходом, приведены теоретические оценки вычислительной сложности алгоритмов и продемонстрированы результаты численных экспериментов. При этом общая идея выполнения преобразования, вне зависимости от алгоритма преобразования, остается неизменной и полностью соответствует определению квантильной функции.
Потребность в реализации квантильной функции возникла при применении непрерывных стохастических алгоритмов оптимизации в двух задачах структурной биоинформатики: задаче предсказания пространственной структуры белка по аминокислотной последовательности и задаче локального пептид-белок докинга с известным линейным интерфейсом связывания и линейной структурой пептида. При этом предполагается, что оптимальное состояние белка или комплекса пептид-белок находится в глобальном минимуме свободной энергии, поиск которого можно осуществить минимизацией энергии посредством глобальной оптимизации.
В силу сложности задач [3] проблематично использовать алгоритмы оптимизации без учета дополнительной информации. Например, в задаче предсказания структуры белка можно значительно упростить пространство поиска, используя низкоэнергетические фрагменты, а также статистически известные распределения для торсионных углов. В задаче пептид-белок докинга возможно исключить из рассмотрения потенциально неподходящие конформации и положения пептида. Чтобы учесть указанные особенности рассматриваемых задач, необходимо наложить взаимозависимые ограничения на большое число контролируемых параметров. Специфика рассматриваемых задач, выражаемая высокой размерностью и большим количеством разнородной дополнительной информации, не позволяет просто использовать ограничения типа неравенств на определенные параметры. Например, в случае пептид-белок докинга положение конечного остатка пептида зависит от вектора смещения, кватерниона вращения и положения углов всех предшествующих остатков, для которых физически доступные значения выражаются известными плотностями распределения. При этом предлагаемые в настоящее время подходы к глобальной
05.13.18
оптимизации, в частности, рассматриваемые эволюционные алгоритмы, требуют непрерывного пространства поиска без ограничений, кроме границ поиска для каждого параметра.
В связи с этим возникает проблема использования дополнительной информации и последующего наложения ограничений, которая разрешима с помощью применения многомерного эмпирического квантильного преобразования. При этом пространство поиска для алгоритмов оптимизации сводится в непрерывный единичный гиперкуб. В предложенном подходе элементы выборки распределены на узлах регулярной сетки, что позволяет с контролируемой для каждого параметра точностью аппроксимировать часть доступного пространства поиска. Следует отметить, что предложенные алгоритмы позволяют производить как непрерывное квантильное преобразование, так и дискретное.
Статья организована следующим образом. В разд. 2 приводится определение квантильной функции и излагается принцип предлагаемого алгоритма преобразования. Демонстрируется псевдокод алгоритма с примерами, которые раскрывают детали предложенного подхода. Разд. 3 посвящен вычислительной сложности рассматриваемых алгоритмов. Использование минимального детерминированного автомата приведено в разд. 4. В разд. 5 и 6 представлены результаты численных экспериментов для рассматриваемых задач биоинформатики. В разд. 7 приводятся результаты выполненной работы и указываются направления дальнейших исследований.
2. Алгоритм квантильного преобразования
Определение многомерной эмпирической квантильной функции (или эмпирического квантильного преобразования) естественно выводится из определения эмпирической функции распределения. Здесь будет приведено наиболее распространенное [4] краткое «рекурсивное» описание структуры квантильной функции.
Пусть дано вероятностное пространство, и на нем определена случайная величина X. Функцией распределения случайной величины X назовем функцию Рх: Я ^ [0, 1], задаваемую формулой Рх (х) = Р (X < х). Для заданной функции квантильное преобразование Р'1 : [0, 1] ^ Я определяется следующей формулой:
FX1 (р) = ^ (х е Я: Р (X < х)> р. (1)
Определим эмпирическую функцию распределения следующим образом:
Рх (а)= * £I(хм< а), (2)
п I=1
где I - индикаторная функция.
Важно отметить, что в предложенном подходе рассматриваются случайные величины, распределенные на узлах регулярной сетки. Причины такого ограничения следующие. Во-первых, в процессе дальнейшего применения построенного преобразования важна только доступность той или иной области поиска. Во-вторых, если значения будут распределены не по сетке, а, например, просто равномерно, то в выражении (2) количество найденных элементов может стать равно нулю, что приводит к неоднозначности распределения.
Необходимо отметить, что каждый узел регулярной сетки определяется несколькими компонентами в зависимости
Полуян С.В., Ершов Н.М.
от размерности. При этом в выборку узел сетки может входить определенное количество раз. Поэтому далее узлы сетки, входящие в выборку, будут называться просто элементами выборки. Каждый узел представляется вектором, а каждый компонент вектора - вещественное значение компоненты узла, либо целочисленный порядковый индекс, который может быть переведен в вещественное значение с использованием информации о регулярной сетке. При этом необходимо хранить всю структуру сетки отдельно, т.е. сеточные узлы в вещественном виде. Для использования рассматриваемых алгоритмов преобразования необходимо указать левую и правую границы и количество сеточных узлов (>1) для каждой размерности.
В настоящей работе для хранения выборки используются различные структуры данных. Для дальнейшего описания алгоритма преобразования необходимо пояснить используемую терминологию. Префиксное дерево будет представлять собой корневое дерево с вершинами и дугами, состоящее из нескольких уровней. Для представления детерминированного конечного автомат (ДКА) будет использоваться диаграмма с переходами и состояниями. Упоминание термина «узел» будет иметь отношение только к сетке.
Первый, простой способ хранения, или просто явный, заключается в хранении элементов выборки в таблице, т.е. представляет собой массив векторов, каждый компонент вектора либо вещественное значение, либо целочисленный индекс узла, который может быть переведен в вещественное значение. Второй способ хранения - префиксное дерево, построенное с использованием целочисленных индексов компонент узлов. Третий способ - модифицированное префиксное дерево. Поскольку каждый узел представим одинаковым количеством компонент, последний уровень префиксного дерева можно объединить. При этом происходит потеря контроля количества элементов на последнем уровне. Такой способ подходит, когда каждый узел входит в выборку только один раз и может использоваться только для «равномерного» распределения. Заключительный способ -представление выборки с помощью минимального ДКА.
Рассмотрим предлагаемый подход к выполнению кван-тильного преобразования на одномерном примере, который представлен на рис. 1. Допустим, имеются два отрезка, представляющие собой непрерывное пространство поиска, к которым нужно организовать доступ посредством кван-тильного преобразования. Если ввести сетку из 12 узлов, можно, с определенной точностью, покрыть данные отрезки. При этом в выборку добавляются только узлы, принадлежащие отрезкам. По выбранным узлам возможно построить эмпирическую функцию распределения дискретной случайной величины. Для выполнения квантильного преобразования необходимо построить функцию, которая производит «обратное преобразование», основываясь на выборке. Если каждый узел будет ассоциирован с двумя границами, которые расположены на расстоянии половины шага сетки от позиции узла, то, используя эти границы, можно произвести линейную интерполяцию и обеспечить непрерывность преобразования.
На данном этапе можно избавиться от явного представления узлов и сформировать выборку, используя порядковые индексы узлов. При этом для каждого индекса необходимо контролировать количество вхождений в выборку. Количество вхождений элементов будет влиять на равномерность распределения. На рис. 1 рассматриваются две выборки, элементы которых, занимая одни и те же узлы, имеют разное количество вхождений. Если проводить индексацию, начиная с нуля, то использованные выборки, содержащие 7 и 16 элементов, представляются следующим образом:
[11, 21, 31, 71, 81, 91, 101], [13, 21, 32, 73, 84, 92, 10,]. (3)
Поскольку предложенный алгоритм квантильного преобразования использует линейную интерполяцию с непрерывным аргументом, возможно оценить плотность распределения. На рис. 1, б, в продемонстрированы ядерные оценки плотности распределения с использованием пакета статистической обработки Я [5], которые демонстрируют зависимость распределения при использовании одной из выборок (3).
1 г
-2 -1
1 х б
-2 -1
yw
Рис 1. Выборка и ядерная оценка плотности распределения:
а - эмпирическая функция распределения; б - равномерное распределение; в - неравномерное распределение
у
У
0
В одномерном случае возможно перейти к представлению с помощью префиксного дерева. При этом каждая вершина дерева будет содержать индекс узла и число вхождений элемента в выборку. Необходимость перехода к дереву будет раскрыта позднее и на данном этапе связана с терминологией описания алгоритма.
Псевдокод процедуры одномерного квантильного преобразования для непрерывного числа из отрезка [0, 1] представлен на рис. 2. Процедура выглядит следующим образом. Пусть дана выборка, состоящая из потомков корневой вершины S, и узлы сетки G, представленные в вещественном виде в массиве. Основной цикл алгоритма представляет
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
собой процедуру двоичного поиска по массиву в. При этом используется вариант двоичного поиска, который выполняет поиск нижней границы, т.е. в результате будет найдена позиция первого узла, который не меньше (равен или больше), чем искомый узел. Необходимо отметить, что искомый узел в явном виде изначально не известен, присутствует только аргумент V е [0, 1], который, по определению эмпирической функции распределения, должен удовлетворять следующему выражению а(т)/п < V < Ь(т)/п, где п - суммарное количество всех элементов; а и Ь - общее число
05.13.18
элементов в узлах меньших т, где т - рассматриваемый с помощью двоичного поиска узел. Таким образом, для каждой операции сравнения требуется выполнить проход по всей выборке и произвести подсчет количества элементов у вершин с меньшим индексом, чем у рассматриваемого. Псевдокод этого действия представлен на рис. 3. После того, как узел найден, необходимо, в случае непрерывного преобразования, выполнить линейную интерполяцию. В случае дискретного преобразования используется индекс найденного узла.
QuantileTransform:
Input: S - текущая вершина, потомки которой определяют выборку (уровень), G - соответствующая уровню сетка, 8 - текущий шаг сетки, v Е [0,1] - аргумент Output: к - состояние для перехода, и - значение после преобразования п S.count > количество элементов (не вершин) на текущем уровне
а 0 Ъ<г- 0 у<- 0
¿«-О f 4- 0 s 0 Count G.size - 1
> узел сетки найден
while Count > 0 do <<-/
s Count/2 i i + s m i
a,b CountLess(5, m x a/n if v > x then y b/n if v < y then
break end
i i + 1 /<-<
CountCount — s — 1 else
Count s end end
if Count = 0 then
y b/n end
26 k 4— 0
27 if a = b then if a = 0 then
k MinimallndexPosition(S') u G[S.children[k].index] + v-8 return u end
if a = n then
k MaximallndexPosition(S') u G[S.children[k].index] + v-8 return u end
k InsideIndexPosition(5) u 4— G[S.children[k].index] + v-8 return £;, u end
for j 1 to S.children.size step 1 do if S.children[j].index = m then k <-j break end end
u G[m] + (v return k, u
> основной цикл поиска
> текущий искомый узел сетки
> число элементов меньше шиш
Компоненты двоичного поиска.
> если все итерации цикла v < х
> индекс узла в выборке для перехода на новый слой
> исключительный случай равенства а и включающий три варианта:
> 1) на левой границе
> поиск минимального узла
> интерполяция на один шаг сетки
> 2) на правой границе
> поиск максимального узла
> интерполяция на один шаг сетки
> 3) на границе между узлами
> интерполяция на один шаг сетки
> поиск узла в выборке для перехода
x)-(G[m + 1] — G[m])/(y — х) > линейная интерполяция
Рис. 2. Псевдокод одномерного квантильного преобразования
Полуян С.В., Ершов Н.М.
Возможна исключительная ситуация равенства а и Ь, которая включает три случая, при которых невозможно использовать линейную интерполяцию. Первые два возникают при указании аргумента на границах (0 или 1). Данные случаи разрешимы нахождением минимального и максимального узла в выборке. Третий случай равенства а и Ь соответствует попаданию на границу между узлами, при этом достаточно выбрать один из узлов. Ввиду очевидности поиска узла в исключительных случаях данная часть в псевдокоде на рис. 2 опущена. Детально с обработкой данных случаев можно ознакомиться в полной реализации [5].
Перейдем к описанию многомерного преобразования. Введем определение многомерной эмпирической функции
^ Х2.....X, ("1/ "1/.../ и ) =
1А (4)
= -X1 ((/ 1] ^ и1/ Х[/ 2] ~ и2/.../ Х[/ ,] ^ )
П I = 1
где d - размерность; п - размер выборки. Теперь можно определить многомерную квантильную функцию [0, 1]d — Я d. Пусть Г - d-мерная функция распределения и X,, ... , Хп - выборка. Используя одномерное квантильное преобразование (1) и выбрав вектор г е [0, можно определить рекурсивно квантильную функцию У = т"1 (г):
У, = Г-1 (г,), (5)
Ук = Fk]l ..., к -1 (.....Ук _ г), 2 < к < d. (6)
Согласно определению (6), одномерное квантильное преобразование применяется d раз для каждой размерности
с использованием только текущей компоненты узла сетки. При этом получаемый после каждого одномерного преобразования компонент сетки используется для отбора элементов из выборки для последующего преобразования.
Ключевая особенность предлагаемого подхода заключается в хранении выборки с помощью префиксного дерева, при котором все элементы выборки разделяются по значению узловых компонент. На этом этапе применяется парадигма разделяй и властвуй, которая естественно сочетается с рекурсивным определением многомерной квантильной функции (6). Важно отметить, что при этом не происходит потери количества элементов на любом уровне дерева, поскольку каждая вершина включает в себя число элементов выборки с тем же порядковым индексом. При этом размер любого уровня дерева и, соответственно, используемая при подсчете (рис. 3) часть выборки, не превышают размера сетки. Вместо поэлементного подсчета, как в случае явного хранения выборки, каждая вершина дерева учитывает количество элементов с аналогичными индексом, что позволяет использовать только находящиеся на уровне вершины. Псевдокод многомерного эмпирического квантильного преобразования представлен на рис. 4.
Следует отметить, что поскольку количество компонент в каждом узле одинаково, для уменьшения количества вершин последний уровень дерева можно объединить. Такое префиксное дерево будем называть модифицированным.
Для оптимизации процесса подсчета возможно хранение уровней дерева в отсортированном виде. Это позволяет использовать двоичный поиск в процессе подсчета. Подробнее этот вариант преобразования представлен в разд. 3.
CountLess:
Input: S - текущая вершина, г - рассматриваемый узел сетки
Output: a, b - общее число элементов с индексом меньше г и г + 1 соответственно
а<- 0 Ь<- 0
1
2 for i i— 0 to S.children.size step 1 do
3 tS.children[i].index
4 if t < r + 1 then
5 b b + S.children[i].count
6 if t < r then
7 | a 4— a + S.children[i\.count end
end
10 end
11 return a, b
> для каждой вершины уровня
> для индексов меньше г + 1
> подсчет количества элементов
> для индексов меньше г
> подсчет количества элементов
Рис. 3. Псевдокод функции подсчета
MultivariateQuantile Transform:
Input: d - размерность, I € [0, l]d - входной вектор, root G - массив сеток для каждой размерности Output: О - входной вектор
1 pf- root
2 for i0 to d step 1 do
k,0[i] QuantileTransform(p.c/n7dren, G[i], /[г]) p <r- p.children[k]
3
4
5 end
6 return О
начальная вершина (состояние),
> начиная с корня дерева
> для каждой размерности
> одномерное преобразование
> переход на следующий уровень
Рис. 4. Псевдокод многомерного квантильного преобразования
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
На рис. 5, б продемонстрировано покрытие непрерывной области поиска выборкой, элементы которой являются узлами регулярной сетки 9 х 10, и показано распределение 2 • 103 точек, полученных с помощью двумерного квантиль-ного преобразования. Основные недостатки предложенного
05.13.18
сеточного подхода, возникающие из-за большого шага сетки, заключаются в том, что, во-первых, часть области поиска оказывается недоступной и во-вторых, присутствует выход за границы непрерывной области поиска. Важно отметить, что в рассматриваемых задачах данные недостатки несущественны.
а б в
Рис. 5. Покрытие непрерывной области поиска и выборка в префиксных деревьях:
а - модифицированное префиксное дерево; б - сетка, выборка и распределенные точки; в - дерево с отсортированными уровнями
На рис. 5 представлены модифицированные префиксные деревья, которые были использованы для распределения точек. В вершинах указаны порядковые индексы узлов сетки с соответствующими цветами. В индексе каждой вершины указывается количество элементов, причем последний уровень содержит каждый элемент в единственном экземпляре, поэтому для него индексация опускается. Важно отметить, что указанное в индексе количество элементов не является числом потомков для каждой вершины и в данном случае совпадает только по причине двухмерности примера. С трехмерным примером, который будет описан в дальнейшем, можно ознакомится на рис. 9. При хранении выборки в модифицированном дереве с отсортированными уровнями в степени указывается количество элементов, которые содержатся в вершинах с меньшим порядковым индексом. Детали применения такого подхода описаны в разд. 3.
3. Вычислительная сложность
Оценим вычислительную сложность предложенного выше алгоритма преобразования, учитывая способ хранения выборки. В первую очередь необходимо отметить, что оценка сложности будет приведена для случая квадратной сетки. При этом в выборку будут добавлены все возможные элементы, т.е. все узлы используемой сетки. Такие допущения не влияют на асимптотику оценок и необходимы только для компактности представления выражений. Следует отметить, что при полном описании сложности необходимо учитывать размеры регулярной сетки (количество узлов на каждом уровне), количество «занимаемых» выборкой узлов на каждом уровне, а также, для случая явного представления, размер выборки.
Рассмотрим квадратную сетку и обозначим за х количество узлов сетки для каждого уровня. Размерность функции обозначим за d. Таким образом, размер выборки составит хd элементов. Вычислительную сложность преобразования определим как функцию двух переменных х и d. Оценка будет произведена по максимальному количеству элементарных операций.
Приведем сложность преобразования при хранении выборки в явном виде [1], т.е. в виде массива векторов или таблицы. Рассмотрим каждый этап работы алгоритма.
1. Для выполнения преобразования необходимо выделить из выборки массив, состоящий из компонент узлов текущей размерности. Изначально массив будет состоять из первых компонент всей выборки. После первого преобразования необходимо перебрать все элементы выборки и выделить с найденным первым компонентом. На следующем шаге поиск происходит среди ранее выбранных для второго компонента и так далее, кроме последнего шага. При этом достаточно каждый раз сохранять выделенные элементы, чтобы избежать полного перебора выборки для каждого преобразования. Таким образом, для выбора необходимых элементов понадобится ^._2X операций, при условии хранения порядковых индексов узлов. Данный этап возможно незначительно ускорить при хранении выборки и преобразовании с использованием графических ускорителей [6].
2. Основной цикл алгоритма (см. рис. 2) выполняется d раз и займет d (х) операций.
3. На каждой операции сравнения основного цикла происходит подсчет элементов с меньшим индексом, который потребует ^2(х)^ = 1 х' операций.
Таким образом, временная сложность преобразования составит:
^x' + d log2 (x) + log2 (xx1 =
i = 2 i = 1 d
= (log2 (x) +1)) x' +(d + x) log2 (x).
(7)
Для преобразования необходим дополнительный вектор длины хd, при этом выборка потребует хранения ^ + 1)хd узловых компонент. Необходимо отметить, что элементы выборки могут быть представлены с использованием вещественного или целочисленного типа данных. Целочисленный способ представления не влияет на сложность преобразования при условии хранения вещественных значений узлов сетки для каждой размерности, т.е. dx значений, и использования индексной адресации. При хранении выборки в вещественном виде выбор элементов потребует в два раза больше операций. Добавление элемента в выборку является
=2
Полуян С.В., Ершов Н.М.
обычным добавлением в массив вектора, и в общем виде потребует сС операций.
Рассмотрим сложность преобразования при хранении выборки в модифицированном префиксном дереве. Необходимо отметить, что сложность преобразования при хранении в префиксном дереве идентична рассмотренной, отличие заключается только в пространственной сложности.
Начиная с корня дерева, производится одномерное преобразование и выполняется переход на следующий уровень. При этом выборка на каждом уровне не превосходит количества узлов сетки для данного уровня. Сложность каждого этапа представлена ниже.
1. Сложность основного цикла алгоритма составит сС 1о§2 (х) операций.
2. Каждая операция сравнения основного цикла потребует С 1о§2 (х) • х операций.
3. Поиск вершины для перехода на следующий уровень потребует ССх операций.
4. Переход между уровнями потребует сС операций.
Таким образом, временная сложность преобразования
составит:
d log2 (x) + d log2 (x) x + dx + d = : (dx + d)(log2 (x) + 1) = d (x + 1) log2 (2x).
(8)
1. На создание частичных сумм потребуется Схопераций.
2. Каждая операция сравнения основного цикла теперь потребует сС 1о§2 (х)(1о§2 (х) + 1) операций.
3. Процедура поиска не изменится и составит сС 1о§2 (х) операций.
4. Поиск вершины для перехода на следующий уровень теперь потребует сС (1о§2 (х) + 1) операций.
5. Переход между уровнями идентичен, сС операций.
Таким образом, временная сложность преобразования
составит:
d (log2 (x) log2 (8x) + x + 2).
(9)
Для преобразования необходимо хранить Сх вещественных значений узлов сетки для каждой размерности. Пространственная сложность преобразования также включает использование дополнительных переменных. В силу их незначительного количества и отсутствия влияния на асимптотику их упоминанием можно пренебречь. Гораздо больший интерес представляет собой пространственная сложность хранения выборки, которая напрямую зависит от реализации модифицированного префиксного дерева и применяемых типов данных. Укажем количество вершин, в каждой из которых необходимо хранить два целочисленных числа и общее число ссылочных типов в дереве.
1. Вершин в дереве будет ^"_0х", а также х вершин на последнем уровне.
2. Общее количество ссылок в дереве составит
У 0х" -1 ■
¿—!" = О
Сложность добавления элемента в дерево составит Сх операций, поскольку на каждом уровне необходимо проверить наличие компоненты узла.
Возможно добиться снижения сложности преобразования. Для этого необходимо хранить вершины в отсортированном виде на каждом уровне префиксного дерева (неважно, модифицированного или нет). Это позволит упростить процедуру выбора элементов с меньшим индексом, и вместо линейного прохода по всем вершинам уровня использовать двоичный поиск. При этом каждая найденная таким образом вершина должна быть ассоциирована с общим количеством элементов с компонентом меньше представленного в вершине индекса. Для этого необходимо перед каждым одномерным преобразованием создавать дополнительный массив длиной х + 1. Данный массив начинается с 0 и представляет собой последовательность частичных сумм ряда, члены которого содержат количество элементов в каждом поддереве. Подробнее отличие такого подхода представлено в примере на рис. 5, в, где элементы данного ряда указаны в степени каждой вершины.
Оценим вычислительную сложность данного подхода.
К пространственной сложности добавится только массив частичных сумм, который можно создавать при преобразовании на каждом уровне.
Для выполнения преобразования необходимо хранить каждый уровень в отсортированном виде, при этом возможны два варианта. При формировании дерева или поддерживать уровни отсортированными или выполнить процедуру сортировки всех уровней после формирования дерева. В первом случае вставка в отсортированную последовательность составит С (1о§2 (х) + 1 + х) операций, при этом, поскольку в текущей реализации использовался векторный список, необходимый для выполнения преобразования, в худшем случае потребуется х довольно затратных операций для смещения элементов на каждом уровне. Во втором случае сложность добавления составит Сх операций, и на заключительном этапе потребуется, например, быстрая сортировка, т.е. х 1о§2 (х) операций для каждого из ^п=0х" уровней (в случае модифицированного дерева). В настоящей реализации предпочтение отдано второму варианту.
На рис. 6 и 7 продемонстрированы теоретические оценки, представленные выражениями (7)-(9), и получаемое усредненное время преобразования 105 равномерно распределенных случайных векторов при изменении размера сетки с фиксированной размерностью 4 и при изменении размерности с фиксированным размером сетки 10. Результаты представлены для следующих реализаций: Е - при хранении выборки в массиве (явный вид), I - при хранении выборки в модифицированном префиксном дереве, / - при хранении в модифицированном префиксном дереве с отсортированными уровнями. Преимущество использования отсортированных уровней достигается только при увеличении размера сетки, причем, как видно на рисунках, в силу особенности реализации и сопутствующих дополнительных затрат у используемого двоичного поиска, на практике необходим больший размер сетки.
Двоичный поиск используется как на этапе поиска узла сетки, так и в процессе подсчета. Вершины на уровнях отсортированы, а узлы сетки распределены равномерно. Возникает вопрос, возможно ли использовать вместо двоичного поиска, например, интерполяционный поиск, и тем самым снизить сложность преобразования? Чтобы дать ответ, рассмотрим каждый этап отдельно.
Во-первых, в рассматриваемых задачах уровни не заполнены целиком. В связи с этим теряется равномерность, и вычислительная сложность интерполяционного поиска значительно деградирует. При этом использовалось только байтовое представление узлов сетки, т.е. уровень мог содержать максимум 256 вершин. При таком размере уровня преимущество использования интерполяционного поиска есть, однако, оно не значительное.
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
107
,g 106
t 105 Q.
<U
55 104
§ 103
O)
s 102
о
* 101 10°
Размер выборки 164
16
Размер сетки а
244
24
324
\ - Е
: - i
32
10-1
Ч 10-2
сс
¡Ig 10-3
CG
§ 10-4 о. ю
ю-5
CL
SI 10-6
s
£ 10-7 10-8
Размер выборки 164
16
Размер сетки б
Рис. 6. Сложность в зависимости от размера сетки (логарифмическая шкала):
а - число элементарных операций; б - усредненное время преобразования
244
24
05.13.18
324
: - Е
: - i
32
4
4
8
8
1
1
Размер выборки
Размер выборки
109 108 107 106 105 104 103 102 101
102
103
104
105
3 4 5
Размерность а
106
107
1
- I
's
—
100
и„ 10-1
СС
1 10-2 а
СО
810-3 а
С
1 10-6 е
at 10-7 10-8
102
103
104
105
345 Размерность б
106
107
Рис. 7. Сложность при изменении размерности (логарифмическая шкала)
а - число элементарных операций; б - усредненное время преобразования
1
2
6
7
1
2
6
7
Во-вторых, основной цикл поиска оперирует с сеткой, но сравниваются не сами значения сетки, а происходит подсчет количества элементов всех вершин с меньшим порядковым индексом, т.е. основной цикл также зависит от распределения вершин на уровне. Необходимо отметить, что в применяемой реализации двоичного поиска сразу происходит поиск нижней границы. Это условие требуется для
обработки исключительных случаев. В данном случае интерполяционный поиск потребует модификации, что частично повысит его вычислительную сложность.
Таким образом, ускорение процедуры преобразования с помощью интерполяционного поиска затруднительно. Необходимо отметить, что преимущество, проявляемое при преобразовании с отсортированными уровнями настолько
Полуян С.В., Ершов Н.М.
незначительно, что в рассматриваемых задачах сортировка не производилась. Поэтому дальнейшая модификация алгоритма представляет только теоретический интерес.
Информацию о количестве элементов в каждой вершине возможно установить как в процессе добавления новых элементов, так и после построения выборки с помощью рекурсивного обхода дерева. При этом, если речь идет о выборке в общем виде с неоднократным вхождением некоторых узлов, возможно только инкрементное пополнение в процессе добавления элементов. Следует отметить, что счетчик в корне дерева содержит количество всех элементов в дереве. Поскольку обход дерева для установления количества элементов в каждом поддереве - довольно простая процедура, к тому же несущественная с точки зрения времени выполнения, в настоящей работе ее описание опущено, с ней можно ознакомиться непосредственно в реализации [5].
На рис. 8, a показано ускорение, получаемое при использовании квантильного преобразования, при котором выборка хранится в модифицированном префиксном дереве относительно явного квантильного преобразования при предсказании структуры белка 3v1a. Приведены усредненные данные по 102 запускам, при фиксированной сетке и размерности 141 изменялось только количество элементов выборки. Следует отметить, что при размере выборки 106 элементов среднее время одного преобразования составило 7,38 с, в то время как при использовании префиксного дерева всего 1,32 • 10-5 с.
Как было указано ранее в данном случае сортировка уровней не дает преимущества в силу слабой заполненности уровней. На рис. 8, б представлена средняя заполненность уровней модифицированного префиксного дерева для выборки размера 106 для пептида 2jsb, длиной 21 остаток. Поскольку выборка построена из структур, которые сформированы из фрагментов длины 7 с учетом только двух углов, кратные уровни демонстрируют большую загруженность. Размер сетки на каждом уровне различается, поэтому приведено нормированное значение, выраженное в процентах. С деталями формирования элементов выборки можно ознакомится в разд. 5. Как видно на рис. 8, б, средняя загруженность уровней минимальна, что приводит к отсутствию преимущества (см. рис. 8, а), что подтверждается теоретической оценкой.
Таким образом, в отличии от явного представления выборки, где сложность преобразования напрямую зависит от количества элементов выборки, предложенный алгоритм преобразования зависит только от размерности функции, сетки и заполненности уровней, встречающихся в процессе преобразования. В рассмотренных задачах применялось квантильное преобразование с использованием модифицированного префиксного дерева и конечного автомата. Детали применения конечного автомата, который позволяет оптимально использовать память для хранения выборки при аналогичной процедуре квантильного преобразования, представлены в следующем разделе.
4. Применение конечного автомата
Недостатком использования модифицированного префиксного дерева является значительный рост количества вершин. Например, в проводимых экспериментах хранение выборки в несколько десятков миллионов занимало более ста гигабайт оперативной памяти, что представляло собой трудности в рамках одного вычислительного узла используемого кластера. В общем случае, при отсутствии полной уникальности компонент у добавляемых узлов, в дереве
встречаются одинаковые и изоморфные [7] поддеревья. При этом возможно перейти к более компактному представлению выборки, руководствуясь теми же принципами, которые лежат в основе перевода префиксного дерева в направленный ациклический граф. Аналогичный подход используется, например, при хранении морфологических словарей.
В настоящей работе вместо направленного ациклического графа будет использован детерминированный ациклический конечный автомат, который, при схожей структуре, позволяет формально поставить задачу минимизации количества состояний. Данный способ компактного представления применим как для хранении выборки, так и для выполнения многомерного квантильного преобразования, которое будет идентично применяемому при хранении выборки в модифицированном префиксном дереве. Однако, с точки зрения реализации, необходимо соблюдать представление, которое «подобно» [8] префиксному дереву. Имеется в виду использование такого представления, где для каждого состояния хранятся исходящие из него ребра. В настоящей работе для построения минимального ДКА использовался инкре-ментный алгоритм [7], производящий минимизацию после добавления каждого элемента.
Определим детерминированный ациклический конечный автомат как набор из пяти элементов (2, О, s, f, б), где 2 - узлы регулярной сетки (алфавит), О - множество состояний, s е О - начальное состояние, f е О - терминальное состояние, б - частичное отображение б : О х 2 ^ О, обозначающее функцию переходов. Переход из одного состояния в другое определяется меткой, содержащей порядковый индекс и число элементов, которые в префиксном дереве хранилась в вершинах.
Рассмотрим пример использования минимального ДКА и трехмерной квантильной функции. На рис. 9, а представлены результаты распределения 104 точек и выборка из 14-ти элементов на регулярной сетке 5 х 5 х 5. Указанный минимальный ДКА (рис. 9, б) соответствует выборке (см. рис. 9, а), где каждый используемый узел представлен с помощью куба для учета границ распределения.
На рис. 9, в для четырехмерного случая с аналогичной квадратной сеткой продемонстрировано изменение количества вершин (состояний) и дуг (переходов) в зависимости от способа хранения выборки. Результаты представлены для модифицированного префиксного дерева (/) и минимального ДКА (А и А'), при случайном (/ и А) и последовательном (А') добавлении всех 54 узлов в выборку. Под последовательным добавлением подразумевается добавление узлов в порядке изменения сначала последней компоненты, а затем наименьшего суффикса. Такой способ добавления имеет частичную аналогию при составлении выборки в задаче пептид-белок докинга. Важно отметить, что, как видно на рис. 9, в, при последовательном добавлении количество состояний и переходов минимального ДКА значительно меньше случайного. Следует отметить, что в данном случае минимальный ДКА, в который добавлены все возможные узлы сетки, содержит только 5 состояний и 20 переходов.
Важно отметить, что уменьшение количества состояний и переходов зависит от конкретного набора добавляемых элементов. Сокращение достижимо только тогда, когда есть схожие и частично совпадающие узлы. При предсказании структуры белка используемые фрагменты схожи по своей структуре. Поскольку для формирования выборки применяются комбинации фрагментов, построение выборки аналогично созданию словаря с использованием одинаковых
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ 05.13.18
словоформ одной и той же лексемы. В задаче пептид-белок докинга, где рассматривается линейное положение пептида, и схожие элементы выборки отличаются по большей части только кватернионом вращения и вектором смещения, преимущество использование ДКА более заметно. На рис. 10 продемонстрировано отношение числа вершин и дуг мо-
дифицированного префиксного дерева к числу состояний и переходов минимального ДКА в рассматриваемых задачах в зависимости от количества элементов, добавленных в дерево и конечный автомат. Как видно на рис. 10, в рассматриваемых задачах использование минимального ДКА снижает расход памяти в разы.
1000 г
5 5 а
161 100 65
10
-\
шРРрт -вершины ^ /ду™ . !состо>ни, - ^Состо»™» : - А' ' переходы переходы
100 200 300 400 500 Количество узлов в
600
Рис. 9. Трехмерное преобразование, используемый ДКА и четырехмерный пример:
а - результат преобразования; б - минимальный ДКА; в - число вершин и состояний
1
1
- Верши - Дуги ны
2-10"
4-10" 6-104 Количество узлов а
8-10"
ф <
^ S-ф ч
^.—.
з Ф
12
10
-Вершины — Дуги
5-104 1-105
Количество узлов б
1,5-105
Рис. 10. Преимущество использования минимального ДКА:
а - при предсказании структуры белка; б - при пептид-белок докинге
6
5
8
4
6
3
4
2
2
1
0
1
1
Основным преимуществом использования минимального ДКА является оптимальное количество состояний и переходов. Причем, по аналогии с предыдущим разделом, если учесть фиксированность сетки, можно указать строго верхнюю границу количества состояний, а именно (согласно [9])
* < * \
£тт(х', 2х -1, (10)
I = о
где х - количество узлов сетки на каждом уровне, С - размерность функции. На рис. 9, в указана граница, равная 65, найденная по формуле (10), и количество вершин, равное 161, в модифицированном префиксном дереве.
Недостатком использования минимального ДКА является время добавления элемента, зависящее от применяемого алгоритма. Поскольку в используемом инкрементном алгоритме построения ДКА минимизация автомата происходит при добавлении каждого элемента [7], сложность этой операции О (С 1о§2 (О)), где О - текущее количество состояний в автомате. Необходимо отметить, что, естественно, минимизацию можно производить поэтапно, например, при достижении определенного количества состояний. Однако, как
было указано ранее, поскольку реализация автомата «подобна» дереву, необходимо использовать соответствующий алгоритм.
Заполнение счетчика количества элементов для каждого перехода выполняется по аналогии с префиксным деревом и реализуется рекурсивной процедурой после построения автомата. Вычислительная сложность преобразования идентична сложности представленной для деревьев и представлена выражениями (8) и (9). При этом присутствует возможность хранения переходов в отсортированном виде.
Квантильная функция может использоваться для аппроксимации распределения на основе известной плотности распределения вероятности. Для этого необходимо добавить в выборку тот или иной узел соответствующее количество раз. При этом в данном случае для хранения выборки подойдет только явное представление или префиксное дерево. В случае префиксного дерева число вхождений элемента в выборку будет контролироваться счетчиком в каждой вершине. Конечный автомат также применим, однако его инкрементное построение с сохранением информации о количестве элементов в каждом поддереве - нетривиальная задача, которая требует разработки специального алгоритма.
Полуян С.В., Ершов Н.М.
Возможно построение минимального ДКА по аналогии с префиксным деревом, где автомат организован таким же образом, как и префиксное дерево, в котором все листья -терминальные состояния, однако, его минимизация с контролем информации о числе элементов также требует разработки специального алгоритма.
На рис. 11 демонстрируется применение квантильно-го преобразования для аппроксимации распределения по представленной плотности распределения с использова-
нием сетки 100 х 100. Благодаря 105 распределенным точкам и ядерной оценки плотности распределения, построенной по точкам с помощью пакета статистической обработки R [5], видно (рис. 11, б) довольно точную аппроксимацию плотности распределения. При этом используются две выборки: первая хранится в префиксном дереве и подходит для аппроксимации, вторая (рис. 11, в) хранится в минимальном ДКА и подходит только для «равномерного» распределения по указанной области.
Необходимо отметить, что форма автомата не влияет на результаты квантильного преобразования. При хранении выборки, элементы которой входят в нее в единственном экземпляре, в явном виде, префиксном дереве, модифицированном префиксном дереве и ДКА результаты преобразования, включая все исключительные ситуации, полностью идентичны друг другу и не зависят от используемого алгоритма. При хранении выборки, элементы которой входят в нее определенное количество раз, в явном виде или префиксном дереве, результаты преобразования также идентичны.
Результаты применения квантильного преобразования представлены в следующих разделах.
5. Предсказание структуры белка
Рассматривается задача предсказания пространственной структуры глобулярного белка по аминокислотной последовательности, которая представляет собой комбинацию из 20-ти стандартных аминокислот, которые при образовании цепи называют остатками. Поскольку в настоящей работе строение пептида и белка идентично в смысле построения главной цепи из остатков, а также из-за отсутствия строгого разграничения белка и пептида, в текущем разделе все структуры будем называть белками.
В соответствии с «термодинамической гипотезой» Ан-финсена предполагается, что нативное уникальное состояние белка находится в глобальном минимуме свободной энергии. Задача оптимизации формулируется как задача поиска этого минимума и ставится в непрерывном пространстве степеней свободы белка. Для полноатомного представления пространственной структуры макромолекулы и вычисления энергии используется фреймворк Rosetta [10]. Детальное описание постановки задачи, степеней свободы белка и силового поля Rosetta представлено в [6; 11].
В настоящей работе обратим внимание только на этап предсказания, непосредственно касающийся использования квантильного преобразования, которое позволяет значительно сократить пространство поиска, исключить невозможные конформации и произвести поиск среди наиболее предпочтительных структур.
Выделяют вторичную структуру белков, такое конформа-ционное расположение главной цепи, которое определяется углами главной цепи белка, т.е. углами каждого остатка. При этом не учитываются боковые цепи остатков и их конформации. Вторичные структуры отличаются регулярной, периодической формой главной цепи. Укладка вторичных структур одной цепи в глобулу определяет пространственное строение белка и называется третичной структурой.
Существуют различные методы предсказания вторичной структуры по аминокислотной последовательности. Различные пакеты программ по предсказанию в большинстве случаев предоставляют унифицированный вывод, представляющий собой предсказание для каждого остатка цепи с указанием, к какой структуре он относится (спираль, лист или ни то, ни другое) и с какой вероятностью (от 0 до 9). В настоящих экспериментах использовался пакет PSIPRED [12], поскольку вывод его формата необходим для генерации фрагментов пакетом Rosetta. Фрагменты - короткие части известных белковых последовательностей, как правило низкоэнергетические, которые выбираются протоколом Rosetta в соответствии с данными PSIPRED. В выборку добавляются значения всех возможных структур, получаемых полной или, в случае большого количества, неполной комбинацией фрагментов.
В настоящей работе поиск глобального минимума производился с помощью адаптивной дифференциальной эволюции JADE [13]. Выбор данного алгоритма обусловлен хорошими результатами [11] в используемом силовом поле.
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
В результате выполненной работы произведена реализация библиотеки для предсказания структур белков. Данная библиотека позволяет применять любой непрерывный стохастический алгоритм оптимизации. Для поиска структуры предоставляется пространство [0,1]п, где п - количество параметров, определяемое предоставленной аминокислотной последовательностью. Данная библиотека реализована на языке С++ и зависит только от пакета Rosetta, который свободно распространяется для академических целей. Для создания пространства поиска также требуется сгенерированный для этой последовательности средствами пакета Rosetta файл с фрагментами, при этом длина фрагментов может быть различной. Регулярная сетка строится в зависимости от вероятности предсказанного значения PSIPRED и может быть настроена вручную. В экспериментах шаг сетки выбирался из отрезка [5°; 10°] и при большей
05.13.00
05.13.18
вероятности предсказания использовался меньший диапазон поиска.
Для боковых цепей выбираются наиболее вероятностные значения углов в зависимости от углов главной цепи. При этом с использованием данных библиотеки Dunbrack [14] также формируется квантильное преобразование для каждого вида остатка.
На рис. 12 представлены основные этапы предсказания структуры. Исходными данными служат последовательность аминокислот в однобуквенном виде (FASTA-формат) и файл с фрагментами. Незначительная часть структур, полученных с помощью комбинации доступных фрагментов, отображена на этапе «пространство поиска». Данные структуры, если быть точнее - значения их углов, порождают выборку. Таким образом, квантильное преобразование позволяет производить поиск в локальной области каждой из этих структур.
Рис. 12. Основные этапы предсказания с использованием квантильного преобразования
На рис. 13 указано среднеквадратичное отклонение атомов главной цепи (Са RMSD) и всех атомов (RMDS) белка относительно его нативной структуры. Представлены предсказанные структуры (светлым) и нативные структуры (темным), взятые из базы данных PDB [15]. Визуализация и пространственное выравнивание выполнено средствами РуМо1 [16]. Для соответствия глобальному минимуму используе-
мого силового поля нативные структуры прошли процедуру локальной оптимизации стандартными средствами пакета Rosetta. Размерность квантильной функции может быть найдена по формуле 3(1 - 1), где I - количество остатков в цепи. Помимо этого, для углов боковых цепей, т.е. для остальных параметров для оптимизации, используются квантильные преобразования с размерностью от 1 до 4.
Рис. 13. Нативные структуры белков (темные) и предсказанные (светлые)
Как видно на представленных рисунках, использование дополнительной информации посредством применения квантильного преобразования позволило успешно предсказывать, с помощью алгоритма оптимизации структуры порядка 20-50 остатков, что в несколько раз больше, чем в проводимых ранее [11] экспериментах с такими же параметрами алгоритма и аналогичном числе эволюционных шагов.
Необходимо отметить, что структуры, образованные с использованием фрагментов, перед добавлением в выборку проходили процедуру фильтрации средствами Rosetta. Физически невозможные структуры с явными «коллизиями» на уровне атомов главной цепи исключались из рассмотрения. Для улучшения предсказаний данный этап может быть значительно модифицирован. Например, можно добавлять только структуры, которые удовлетворяют основным
Полуян С.В., Ершов Н.М.
особенностям строения белка. В общем случае выборка может генерироваться с помощью определенной математической модели, поскольку с точки зрения применения квантильного преобразования не имеет значения, каким образом была сформирована выборка.
6. Пептид-белок докинг
В задаче пептид-белок докинга необходимо найти оптимальное место связывания белка и пептида при взаимодействии друг с другом, а также соответствующую этой связи конформацию комплекса. При этом известен интерфейс связывания, т.е. предполагаемое место связывания пептида и его «линейная» [17] конформация. Задача оптимизации формулируется как задача минимизации энергии связывания, которая вычисляется как разница между энергией комплекса в связанном состоянии и энергией в свободном состоянии, т.е. когда белок и пептид друг с другом не взаимодействуют. К параметрам пептида, указанным в предыдущем разделе, добавляется вектор смещения пептида относительно белка, кватернион вращения и углы боковых цепей белка, находящиеся в интерфейсе связывания в непосредственной близости к пептиду. Для представления структур в полноатомном разрешении и вычисления энергии использовалось силовое поле Rosetta. С детальным описанием задачи можно ознакомиться в [1].
В силу описанных выше особенностей при докинге можно исключить из рассмотрения положения пептида вне места связывания. На данном этапе используется многомерное квантильное преобразование. Вводится предположение, что при связывании концы пептида не должны выходить за определенные границы, например, границ двух сфер, ориентировочно сформированных на границах интерфейса связывания. Рассматриваются оба положения пептида в указанных границах, которые определяются углами главной цепи пептида, параметрами смещения и вращения. На начальном этапе выборка универсальна и зависит только
от длины пептида, диаметра сфер и сетки. Для создания выборки используется многомерный аналог алгоритма заливки [1]. Затем с использованием карт Рамачандрана из выборки исключаются менее предпочтительные структуры.
В вычислительных экспериментах рассматривались комплексы 1jwg и 2fse с длиной пептидов 5 и 15 соответственно. Размерность квантильной функции может быть найдена по формуле 2(1 - 1) + 7, где l - длина пептида. В текущих экспериментах размерность составила 15 и 35 соответственно. В отличии от задачи предсказания структуры бека, где в квантильное преобразование добавлялись все углы главной цепи, в данном случае исключаются углы, которые «стремятся» к планарности (находятся в trans-конформации), их незначительное отклонение существенно не влияет на конформацию пептида. При этом при оптимизации данные параметры учитываются.
Для докинга использовалась адаптивная дифференциальная эволюция (JADE). Сравнение результатов производилось с протоколом Rosetta FlexPepDock (FPD) [18] с двумя стартовыми позициями. Первая позиция соответствует положению пептида, максимально приближенному к нативному, перевернутому вокруг оси межу сферами. Вторая позиция соответствует положению пептида в обратном направлении в интерфейсе связывания. Условием применения данного протокола является присутствие пептида в радиусе пяти ангстрем от места связывания. Работа протокола включает в себя несколько различных этапов, заключительная стадия которого - алгоритм Монте-Карло с локальной оптимизацией.
На рис. 14 и 15 продемонстрированы интерфейсы связывания комплексов и результаты докинга. На графиках представлена энергия комплекса (точнее скоринг-значение [11]) и получаемое среднеквадратичное отклонение атомов главной цепи. Результаты для протокола FPD с двух стартовых позиции отмечены красным и синим цветом соответственно. Количество запусков FPD соответствует времени 10 запусков JADE. Приемлемым результатом докинга является субангстремное значение отклонения.
• FPD (flip + rotate) • FPD (flip) • JADE
0 2 4 6 8 10 12 CA RMSD
а б
Рис. 14. Результаты докинга для комплекса 1jwg:
а - интерфейс связывания и искомое положение пептида; б - среднеквадратичное отклонение атомов и значения скоринг-функции
Результаты численных экспериментов показывают, что с применением квантильного преобразования эволюционный алгоритм оптимизации способен успешно произвести докинг только для случая комплекса 1^. В проводимых ранее экспериментах [6] с тем же комплексом рассматриваемый алгоритм с аналогичными параметрами с поставленной
задачей не справлялся. В случае комплекса 2fse при увеличении оптимизируемых параметров приблизительно в три раза дифференциальная эволюция не способна произвести докинг. Следует отметить, что это указывает на неспособность алгоритма оптимизации преодолеть сложный энергетический ландшафт рассматриваемой задачи.
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ 05.13.18
• FPD (flip + rotate) • FPD (flip) • JADE
Рис. 15. Результаты докинга для комплекса 2fse:
а - интерфейс связывания и искомое положение пептида; б - среднеквадратичное отклонение атомов и значения скоринг-функции
Многоэтапный протокол FPD с поставленными задачами докинга справился. Причем для комплекса 1jwg вне зависимости от стартовой позиции получено субангстремное отклонение. Для комплекса 2fse получены приемлемые результаты только для первой позиции.
В результате выполненной работы произведена реализация библиотеки для выполнения пептид-белок докинга. Данная библиотека зависит только от пакета Rosetta и позволяет применять любой непрерывный стохастический алгоритм оптимизации. Для поиска предоставляется пространство [0, 1]n, где n - количество параметров, определяемое предоставляемым пептид-белок комплексом.
Расчеты проводились на базе гетерогенной вычислительной платформы HybriLIT (ЛИТ, ОИЯИ) [19, 20]. Длительность одного отдельного эксперимента на одном вычислительном узле составило одни сутки. Например, один из стартов FPD или 10 запусков JADE.
7. Заключение
В результате выполненной работы получены следующие результаты. Предложен алгоритм многомерного эмпирического квантильного преобразования. Исследованы вычислительные особенности алгоритма и проанализирована вычислительная сложность преобразования в зависимости от способа представления выборки. Приведенные в работе результаты численных экспериментов подтвердили корректность представленных теоретических оценок.
Предложенный алгоритм квантильного преобразования с использованием конечного автомата позволяет эффективно хранить выборку и выполнять преобразование, что предоставило возможность применять непрерывные стохастические алгоритмы оптимизации для решения задач структурной биоинформатики. На основании проведенного исследования можно заключить, что использование квантильного преобразования позволило учитывать дополнительную информацию для сокращения пространства поиска и успешно применять эволюционные алгоритмы оптимизации для нахождения структур белков. Представленные результаты демонстрируют, что применение квантильного преобразования позволяет рассматривать только потенциальные положения пептида в интерфейсе связывания, что делает возможным производить докинг коротких пептидов.
Разработка алгоритма инкрементного построения минимального ДКА для хранения выборки в общем случае и исследование использования квантильного преобразования
как алгоритма оценки распределения в рамках задач оптимизации могут служить целями дальнейшей работы. Предложенные реализации квантильного преобразования, численные эксперименты, результаты которых были представлены, а также библиотеки для рассмотренных задач биоинформатики с возможностью применения различных алгоритмов оптимизации доступны в открытом доступе в ре-позиториях [5] сервиса GitHub.
Литература / References
1. Полуян С.В., Ершов Н.М. Применение многомерной квантильной функции в задаче пептид-белок докинга // Вестник Южно-Уральского государственного университета. Серия «Вычислительная математика и информатика». 2019. Т. 8, № 2. С. 63-75. [Polu-yan S.V., Ershov N.M. The use of multidimensional quantile function in the peptide-protein docking problem. Bulletin of the South Ural State University. Series "Computational Mathematics and Computer Science". 2019. Vol. 8, No. 2. P. 63-75.]
2. O'Brien G.L. The Comparison Method for Stochastic Processes. The Annals of Probability. 1975. Vol. 3. No. 1. P. 80-88.
3. Rentzsch R., Renard B.Y. Docking small peptides remains a great challenge: an assessment using AutoDock Vina. Briefings in Bioinforma-tics. 2015. Vol. 16. No. 6. P. 1045-1056.
4. Einmahl J.H.J., Mason D.M. Generalized Quantile Processes. The Annals of Statistics. 1992. Vol. 20. No. 2. P. 1062-1078.
5. Репозитории GitHub. URL: https://github.com/poluyan (дата обращения: 10.12.2019). [GitHub repositories. URL: https://github.com/ poluyan (accessed: 12/10/2019).]
6. Poluyan S., Ershov N. Parallel evolutionary optimization algorithms for peptide-protein docking. EPJ Web of Conferences. 2018. Vol. 173. P. 06010-06010.
7. Daciuk J., Mihov S.W., Watson B.W., Watson R.E. Incremental Construction of Minimal Acyclic Finite-State Automata. Computational Linguistics. 2000. Vol. 26. No. 1. P. 3-16.
8. Ferrandez A., Peral J. MergedTrie: Efficient textual indexing. PLoS One. 2019. Vol. 14. No. 4.
9. Kjos-Hanssen B., Liu L. The number of languages with maximum state complexity. Theory and Applications of Models of Computation. TAMC 2019. Lecture Notes in Computer Science. 2019. Vol. 11436. P. 394-409.
10. Alford R.F., Leaver-Fay A., et al. The Rosetta all-atom energy function for macromolecular modeling and design. Journal of Chemical Theory and Computation. 2017. Vol. 13. No. 6. P. 3031-3048.
11. Полуян С.В., Ершов Н.М. Применение параллельных эволюционных алгоритмов оптимизации в задачах структурной биоинформатики // Вестник УГАТУ. 2017. Т. 21, № 4. С. 143-152. [Poluyan S.V., Ershov N.M. Application of parallel evolutionary optimization algorithms in problems of structural bioinformatics. Bulletin of USATU. 2017. Vol. 21. No. 4. P. 143-152.]
Полуян С.В., Ершов Н.М.
12. Buchan D., Jones D. The PSIPRED Protein Analysis Workbench: 20 years on. Nucleic Acids Research. 2019. Vol. 47. No. W1. P. W402-W407.
13. Zhang J., Sanderson A. JADE: Adaptive differential evolution with optional external archive. IEEE Transactions on Evolutionary Computation. 2009. Vol. 13. No. 5. P. 945-958.
14. Shapovalov M., Dunbrack R.L. A Smoothed Backbone-Dependent Ro-tamer Library for Proteins Derived from Adaptive Kernel Density Estimates and Regressions. Structure. 2011. Vol. 19. No. 6. P. 844-858.
15. Berman H.M., Westbrook J. et al. The Protein Data Bank. Nucleic Acids Research. 2000. Vol. 28. No. 1. P. 235-242.
16. Schrödinger LLC. The PyMOL Molecular Graphics System. URL: https://pymol.org (дата обращения: 10.12.2019).
17. Sellers M.S., Hurley M.M. XPairIt Docking Protocol for peptide docking and analysis. Molecular Simulation. 2015. Vol. 42. No. 2. P. 149-161.
18. Raveh B., London N., et al. Rosetta FlexPepDock ab-initio: Simultaneous Folding, Docking and Refinement of Peptides onto Their Receptors. PLoS ONE. 2011. Vol. 6. No. 4.
19. Adam G., Bashashin M. et al. IT-ecosystem of the HybriLIT heterogeneous platform for high-performance computing and training of IT-specialists // Selected Papers of the 8th International Conference "Distributed Computing and Grid-technologies in Science and Education" (GRID 2018). 2018. P. 638-644.
20. Heterogeneous Computing Cluster HybriLIT. URL: http://hybrilit.jinr. ru/en/ (дата обращения: 10.12.2019).
Статья поступила в редакцию 16.12.2019, принята к публикации 26.12.2019 The article was received on 16.12.2019, accepted for publication 26.12.2019
РЕЦЕНЗИЯ
на статью С.В. Полуяна и кандидата физико-математических наук Н.М. Ершова «Квантильное преобразование в задачах структурной биоинформатики»
В представленной работе рассматривается подход к применению непрерывных стохастических алгоритмов глобальной оптимизации для решения двух задачах структурной биоинформатики, основанный на использовании многомерного эмпирического квантильного преобразования. Данное преобразование позволяет принять во внимание известную информацию о рассматриваемой задаче и значительно сократить пространство поиска для применяемого метода оптимизации. Авторы обосновывают необходимость использования данного подхода и предлагают реализацию квантильного преобразования на основе регулярной сетки с применением парадигмы разделяй и властвуй. Предложенный алгоритм преобразования подробно описан, приведены теоретические оценки вычислительной сложности и результаты численных экспериментов. Для рассматриваемых задач демонстрируются результаты использования эволюционного алгоритма оптимизации, подтверждающие эффективность применения квантильного преобразования.
Авторами детально продемонстрирована полная цепочка формирования алгоритма - от определения одномерной эмпирической функции распределения до представления многомерного преобразования с использованием оптимальной структуры хранения данных в форме конечного автомата. В работе рассмотрены две актуальные задачи биоинфор-матки в рамках современного силового поля и приведены иные области применения предложенного подхода.
Реализация алгоритма квантильного преобразования, разработанные для задач программные библиотеки и результаты численных экспериментов представлены авторами в открытом доступе, что подчеркивает верифицируемость проведенного исследования и достоверность полученных результатов. Содержание работы полностью соответствует указанной специальности.
Считаю, что работа может быть рекомендована к публикации в журнале «Computational nanotechnology».
Кандидат физико-математических наук; доцент кафедры СКИ факультета ВМК МГУ им. М.В. Ломоносова
Н.Н. Попова
Сведения об авторах / About the authors
Полуян Сергей Владимирович, ассистент, кафедра распределенных информационно-вычислительных систем, Институт системного анализа и управления, Государственный Университет «Дубна». Дубна, Московская область, Российская Федерация
Sergey Vladimirovich Poluyan, teaching assistant, Dubna State University. Moscow Region, Dubna, Russian Federation РИНЦ AuthorlD: 1024957 E-mail: svpoluyan@gmail.com
Ершов Николай Михайлович, кандидат физико-математических наук, старший научный сотрудник, факультет вычислительной математики и кибернетики, Московский государственный университет им. М.В. Ломоносова. Москва, Российская Федерация
Ershov Nikolay Mikhaylovich, PhD.; senior research associate of Faculty of Computational Mathematics and Cybernetics, Lomonosov Moscow State University. Moscow, Russian Federation E-mail: ershov@cs.msu.ru
