CC BY
53
Статистический вывод параметров селекции иммунных рецепторов с помощью алгоритма градиентного спуска
с переменным шагом
Офицеров Е.П., ТулГУ [email protected] Воеводская А.А., РЛЯЕХЕЬ vaa23 9@gmail. сот Назаров В.И., НИУ "Высшая Школа Экономики" [email protected]
ИБХ РАН
Аннотация
Иммунитет - сложнейшая система организма человека. Технологии секвенирова-ния нового поколения позволили получать большие массивы данных Т-клеток и Б-клеток, реализующие самообучающуюся ветвь иммунитета - адаптивный иммунитет. В данной работе мы применяем метод градиентного спуска с переменным шагом для вывода параметров модели селекции - процесса, влияющего на появления и численность клеток в крови организма.
1 Введение
Человеческий организм ежемоментно подвергается атаке самых разнообразных патогенов - вирусов, бактерий, грибов, различных паразитов. Несмотря на то, что только сверхмалая часть патогенов проникает непосредственно в организм, они способны нанести обширный урон. Задачей иммунной системы является распознавание и уничтожение патогенов. Иммуная система условно делится на две ветви, которые выполняет эту задачу разными путями. Первая ветвь - это врожденный иммунитет, который располагает механизмами выявления и уничтожения заранее заданных патогенов. Эти механизмы изначально заложены в геноме человека. Но поскольку геном обладает ограниченным информационным объемом, а патогены обладают способностью мутировать, то в нем невозможно закодировать рецепторы, распознающие абсолютно все возможные патогены. Поэтому в ходе эволюции появилась другая ветвь иммунной системы - адаптивный иммунитет, который позволяет организму обучаться защите от наиболее распространенных патогенов в среде. Мощь адаптивного иммунитета человека реализуется с помощью Т-клеток и Б-клеток, несущие на себе специальные белки-рецепторы, позволяющие рас-
познавать патогены и запускать иммунныи ответ. Такие рецепторы называется, соответственно, Т-клеточные рецепторы (T-cell receptors, TCR) и иммуноглобулины (Immunoglobulins, Ig). Они генерируются случайным образом в стохастическом процессе У(Б)1-рекомбинации, в котором часть генома Т-клеток или Б-клеток перестраивается и транслируется в белок-рецептор. В этом процессе случайным образом выбираются по одному геномному сегменту из множества сегментов (под названиями V(ariable), D(iversity) или J(oining), получившиеся сегменты случайным образом изменяются на концах и конкатенируются в итоговую нуклеотидную последовательность, которая далее транслируется в белок, если это валидная последовательность. В случае TCR получившиеся белки проходят следующий этап - селекцию в специальном органе (тимусе), которая уничтожает аутоиммунные рецепторы (способные атаковать сам организма) или рецепторы, принципиально не способные к распознаванию. Технологии секвенирования нового поколения позволили получать нуклеотидные последовательности рецепторов TCR и Ig из периферической крови человека, что открыло возможность для статистического анализа иммунологических данных и математического моделирования на уровне клеток поведения иммунной системы. В данной работе мы применяем алгоритм градиентного спуска с переменным шагом для вывода параметров математической модели влияния селекции на репертуар рецепторов.
2 Селекция иммунных рецепторов
Потенциальное разнообразие иммунных рецепторов, которые собираются в процессе V(D)J-рекомбинации, огромно, к примеру, для Т-клеточных рецепторов бета-цепей оно оценивается как 1014 [Murugan et al., 2012].
Часть этих рецепторов нефункциональна: не способна среагировать ни с каким возможным патогеном, а часть являются аутоимун-ными (т.е. нападающими на организм хозяина (pre-selection на рис. 1). В силу этого реальное число рецепторов на порядки меньше (naive на рис. 1). Число же рецепторов, которые распознали какой-либо патоген и увеличили свою численность в ответ на него, еще меньше (antigen-experienced на рис. 1). Влияние селекции проявляется в отборе тех, кто потенциально может помочь и не навредить организму (thymic selection), и в отборе тех, кто уже встретился с патогенами и может в будущем создать лучшую защиту против него (clonal selection).
TCR repertoire size
■ Pre-selection Naive
Antigen-experienced
TCR repertoire composition
мим««« ew?« « ее Мише?« «w Ъ в -Сфг в «ос«о. 0 в W9« ««'.'«« э W WHiwwwwii 4* «Ч ШWW WîW 9 W « MW « 9 « в 9 в <3 9(3 9 9SW » UMSO WW«W 9 9 9 999 99 ¡3 çe 9999 9 9 99 99 95 9999 999 999 99 9 95 9999999 > 9 99 99 99 99
Thymic selection Clonal selection
Pre-selection
Antigen-experienced
Рис. 1. Соотношение размеров репертуаров TCR. Наибольшее разнообразие TCR находится в пре-селектированном репертуаре. Селекция в тимусе убирает нефункциональные и аутоиммунные рецепторы, что снижает разнообразие репертуара. После встречи с патогенами, рецепторы увеличивают свою численность, вытесняя наивные рецепторы (которые не встречались с патогенами ранее).
В статье [Elhanati et а!., 2014] была предложена математическая модель селекции, которая позволяет количественно оценивать влияние селекции на клонотипы - определенные последовательности иммунных рецепторов с фиксированными V и I сегментами. Модель описывает селекцию как:
9L4vj
. где Q(T> V'J> - влияние селекции на клоно-тип Тс сегментами VhJ,
'■■■■* - вероятность встретить клонотип Г с сегментами V и J в экспериментальных данных, т.е. его пропорция в данных,
- вероятность генерации клонотипа
с сегментами V и J, X '- нормализующая константа для выполнения условия
£ Q<T,V,J)rjT,VJ)=l
,
- влияние селекции на длину L клонотипа
,
- влияние селекции на сегменты V и J клонотипа 7 ,
: - влияние селекции на i-ую аминокислоту в клонотипе ' с длиной L. Таким образом, для любого клонотипа, зная вероятность его сборки в процессе V(D)J-рекомбинации Р^Л^Л^) и селекционную модель Q, возможно оценить его потенциальную представленность в данных как P^JT, V,J)=P„rjT, VJ)XQ(z, V,J)
3 Методы оптимизации
Задачи оптимизации решаются в основном методами первого порядка (в которых вычисляются первые производные целевых функции) и методами второго порядка (в которых вычисляются и первые, и вторые производные). Одним из самых популярных методов первого порядка является метод градиентного спуска, каждая итерация которого выглядит как
w=w-A ^ V/(4vuJ
,
где - шаг алгоритма.
Этот метод является очень простым в реализации и понимании, но имеет серьезные недостатки: длительность работы, затратность вычислений, необходимость обрабатывать весь массив данных на каждой итерации метода. Существует множество модификаций градиентного спуска, например, метод покоординатного спуска [Wright, 2015]. Его отли-
X
чие заключается в подборе на каждом шаге. В этом алгоритме подбирается так, чтобы в направление антиградиента функция
f достигла минимума, т.е.
Xk=arRmin^JJ(xk-XNf(wtxk))) в общем
случае скорость метода не сильно превышает скорость обычного градиентного спуска.
Ещё одной вариацией градиентного спуска является метод сопряженных градиентов [Nocedal J., Wright S. J.,2006].
Изначально метод применялся для решения квадратичной задачи оптимизации, то влекло за собой решение СЛАУ Ах=Ь Решение задачи высчитывается по формуле 1 ' ■ ' . На каждом шаге вы-
считываем коэффициенты
a*=argm maf (хк_, +al¿>t) _ находим базисный вектор Pk=~f (xí-i)+Pi iPi-j ^ а коэффициент можно высчитывать разными способами, но одним из самых часто употребимых является Метод Флетчера - Ривса по которому
тчкичЬ
.
В общем случае скорость сходимости метода сопряженных градиентов выше скорости покоординатного спуска.
Потомком метода градиентного спуска является метод стохастического градиента [Bot-tou L., 2010]. Он позволяет обрабатывать данные в потоковом режиме, в котором массив данных обрабатывается поэлементно, без необходимости загрузки всего массива в память компьютера. На каждом шаге алгоритма стохастического градиента высчитывается градиент не от всего множества, а от некоторой случайной точки. В этом случае шаг будет выглядеть следующим образом:
Одним из главных минусов метода стохастического градиента является медленная асимптотическая сходимость, связанная с дисперсией. Метод основан на равенстве EfV.f(w~-i>x^',-l] = V.f(w,-l) _ в котором каждое выражение ^ ■■'1 , конечно, от-
личается. На практике это означает, что, если значение V/niW велико, то увеличивает-
ся и дисперсия, а значит, понижается скорость сходимости. Методы, рассмотренные далее, предлагают решение этой проблемы.
Метод усредненного стохастического градиента (SAG) [Schmidt, Roux, Bach, 2013] предназначен для оптимизации суммы конечного числа гладких выпуклых функций. Так же, как и метод стохастического градиента, стоимость итерации не зависит от числа слагаемых. Особенностью данного метода является запоминание предыдущих значений градиента. За счет этого метод сходится гораздо быстрее метода стохастического градиента, а во многих случаях даже быстрее метода градиентного спуска. Шаг алгоритма будет выглядеть следующим образом: d=d-yî
d-ii+yi
.
Переменная d используется для отслеживания количества ^ за счет хранения в памяти последнего вычисленного значения градиента для каждого индекса итерация занимает меньше времени, чем в методе стохастического градиента. Если есть возможность пройти по данным более одного раза, то данный метод быстрее предыдущих методов.
Для уменьшения дисперсии также был предложен метод стохастического градиента с сокращенной дисперсией (SVRG) [Johnson, Zhang, 2013]. Для гладких и строго-выпуклых функций этот метод сходится со скоростью SAG, но интуитивно понятнее и не требует хранения всех градиентов, а следовательно, метод применим и к решению задач с большим число параметров. Также метод применим к задачам невыпуклой оптимизации. Но всё же он требует хранения информации, а значит теряет преимущество потоковой обработки информации, которой обладает метод SGD.
Рис.2. Зависимость логарифма правдоподобия от количества итераций. Алгоритм с переменным шагом
сходится быстрее алгоритмов с постоянным шагом
Алгоритм SVRG метода заключается в следующем. Раз в m шагов мы запоминаем веса и
откладываем их в переменную и и градиент на них — V.1 ■' Тогда шаг итерации будет выглядеть как
.
Метод SVRG-LBFGS [Kolte, Erdogdu, Ozgür, 2015] метод является улучшением метода SVRG и использует BFGS — алгоритм Бройдена - Флетчера - Гольдфарба - Шанно с ограничением по памяти [Liu, Nocedal, 1989]. Алгоритм BFGS заключается в разложении функции в полином второй степени и приблизительную оценку Гессиана. В данном методе применяется LBFGS в начале, используя полный градиент, что приводит к менее частому пересчитыванию Гессиана, увеличению стабильности результатов и скорость сходимости увеличивается.
Из методов второго порядка для решения задач большой размерности применяются вариации метода Ньютона, он основан на ис-
пользовании квадратичной аппроксимации функции в окрестности текущей точки " . Одна итерация в методе Ньютона выглядит
как ^ м . где ■'■ - гессиан, градиент функции в точке Хк .
Метод Ньютона требует вычисления на каждом шагу вторых производных для каждой переменной, что может быть трудоемко, но при этом он сходится гораздо быстрее метода градиентного спуска.
Для решения данной задачи французскими исследователями в [Elhanatl У. et al., 2014] был применен метод градиентного спуска. К примеру, для оптимизации давления селекции на длины используется следующая формула:
Чг^ЛсП^и-Р^ЬО')]
В данной работе предлагается применить метод градиентного спуска с переменным шагом, который в данной задаче требует сопоставимых затрат с алгоритмом градиентного спуска с постоянным шагом, но сходимость в общем случае происходит после меньшего числа итераций.
Рис.3. Совместные распределения параметров моделей селекции, полученные разными алгоритмами. Точки соответствуют параметрам моделей селекции, которые относятся к давлению селекции на аминокислоты в определенных позициях и длинах.
4 Алгоритм наискорейшего спуска для вывода параметров модели селекции
Задача поиска оптимальных параметров заключается в максимизации функции правдоподобия
где
Рthij.^l.;) _
'iV.. .
доля аминокислотных после-
L
в выборке экспе-
рге
s&ata
Это эквивалентно максимизации логарифма правдоподобия:
lnfL)=Y Inf Q(s))+ У Inf Р (s))^ шах
¿W .11 lг_, s^Ldam 1
Так, как от параметров селекции зависит только первая сумма, то в результате требуется решить задачу оптимизации:
Для решения этой задачи применяется метод наискорейшего градиентного спуска. Он заключается в том, чтобы на каждой итерации двигаться вдоль направления градиента. Для рассматриваемой функции частные производные будут иметь следующий вид:
довательностеи длины риментальных данных, ?{ц{ |0)
- вероятность появления последовательности такой длины при условии набора параметров ®. Аналогично вычисляются производные по каждому из параметров селекции.
На практике для оценки значения
Р(с{) |Э;
используется следующий подход. На основе модели генерации инициализируется выборка аминокислотных последовательностей О. После этого, вероятность может быть оценена как,
Совместные распределения параметров моделей селекции, полученные разными алгоритмами (гены)
А1 (211877)
А2 (110505)
С1 (457362)
/о /
4
(0.99841
0.5 1.0 1.5 р (перем. шаг) С2 (396307)
о/
и
/ / [0.9451 ]
у
а
-я
У 0
0.5 1.0 1.5 д (перем. шаг) 01 (350623)
л /
0 /
¡0.9964]
--э
/ [0.995]
0.5 1.0 1.5 Ч (перем. шаг)
02 (397691)
о/' о
•
» в' ® у @
[0.7641]
q (перем. шаг)
0.5 1.0
д (перем. шаг)
q (перем. шаг)
Рис.4. Совместные распределения параметров моделей селекции, полученные разными алгоритмами. Точки соответствуют параметрам моделей селекции, которые относятся к давлению селекции на определенные гены.
Шаг вдоль направления градиента осуществляется пропорционально некоторому коэффициенту "■ . Этот шаг, в свою очередь, при известном градиенте функции определяется так, чтобы максимизировать значение логарифма правдоподобия в направлении градиента. Таким образом, выбор коэффициента представляет собой одномерную задачу оптимизации, которая решается с помощью метода золотого сечения.
5 Результаты
Поскольку в статье [Elhanatl et al, 2014] программный код для вывода параметров был написан на языке Matlab и специлизирован под определенные данные, к которым мы не смогли получить доступ, мы реализовали алгоритмы градиентного спуска с постоянным и переменным шагом на языке С++ и применили их к данным ТСЯ альфа-цепи трех пар однояйцевых близнецов из [Zvyagln et я1, 2014]. Для каждого из людей мы сгенерировали по 2.5 - 3 миллиона рецепторов по вероятностной модели сборки рецепторов [Murugan et я1,
2012]. На (рис. 2) видно, что, как и ожидалось, алгоритм с переменным шагом сходится за гораздо меньшее число итераций, чем с
постоянным. При этом шаг размером в 2 оказывается слишком большим, и алгоритм начинается двигаться вокруг экстремума, не сходясь к нему.
Совместные распределения параметров
(Чч
модели селекции по аминокислотам ( ,
рис. 3), генам (С' '■' , рис. 4) и длинам ( ^ '■ , рис. 5) очень похожи у всех людей, кроме человека С2. Похожесть параметров позволяет сделать вывод, что оптимизация модели с помощью градиентного спуска с переменным шагом является гораздо более эффективной процедурой вывода параметров модели, нежели градиентный спуск с постоянным шагом. В случае человека С2 алгоритм с переменным шагом очень быстро сошелся к оптимуму, при этом алгоритм с постоянным шагом 2 начал расходиться после 40й итерации, что может говорить о сложном ландшафте оптимизируемой функции для человека С2.
Совместные распределения параметров моделей селекции, полученные разными алгоритмами (длина) А1 (211877) А2 (110505) С1 (457362)
1.2
3 1.0
0.6
О'' ✓ О /
аз г О ' У
S/ 9/
У / [ 0.9949 \
1.2
g 1.0
0.6
0.6 0.8 1.0 1.2 L (перем. шаг)
С2 (396307)
,9 у fi
/
Л о ?
[0.9999]
1.2-
3 1.0
0.6
rf О
/ О/ Q /
S ' / в0'
✓ ( 0.9886
0.6 0.8 1.0 1.2 L. (перем. шаг)
D1 (350623)
0.6 0.8 1,0 1.2 L (перем. шаг)
D2 (397691)
3 1.0
✓ 1.2 ¡ 1.0-ь О 1.0.8 0.6
s А 9 У / "b
/ / © О® ш/
/ * & ✓ О у
/ © @ О 0 ( 0.76361 I 0.9911 ]
э 1.0
■h
0,' Л О
© «
f I O.9799]
0.8 1.0 L (перем. шаг)
0.8 1.0 L (перем. шаг)
0.8 10 1,2 L (перем. шаг)
Рис. 5. Совместные распределения параметров моделей селекции, полученные разными алгоритмами. Точки соответствуют параметрам моделей селекции, которые относятся к давлению селекции на определенные длины.
6 Заключение
В данной работе был реализован алгоритм наискорейшего градиентного спуска для статистического вывода параметров модели селекции иммунных рецепторов. Реализованные алгоритмы работают очень быстро (2-5 минут на одного человека) и занимают мало оперативной памяти (500 - 600 мегабайт) для большого объема сгенерированных данных по вероятностной модели сборки рецепторов (2.5 - 3 миллиона рецепторов). Было показано, что для данной задачи метод с переменным шагом сходится гораздо быстрее метода, реализованного в предыдущей работе. В дальнейшем планируются изучить степень влияния размера сгенерированной выборке на точность оценки параметров, и разработать алгоритмы, которые позволят эффективно обрабатывать огромные массивы экспериментальных и сгенерированных данных.
Список литературы
Meriem Attaf, Eric Huseby and Andrew K Sewell.
2015. ab T cell receptors as predictors of health
and disease. Cellular & Molecular Immunology,
12, 391 - 399.
Leon Bottou. 2010. Large-scale machine learning with stochastic gradient descent. Proceedings of COMPSTAT'2010. Physica-Verlag HD, 177-186.
Yuval Elhanati et al. 2014. Quantifying selection in immune receptor repertoires. Proceedings of the National Academy of Sciences,111, 27, 9875-9880.
Rie Johnson, Tong Zhang. 2013. Accelerating stochastic gradient descent using predictive variance reduction. Advances in Neural Information Processing Systems, 315-323.
Ritesh Kolte, Murat Erdogdu, Ayfer Ozgur. 2015. Accelerating SVRG via second-order information. NIPS'15.
Dong C. Liu, Jorge Nocedal. 1989. On the limited memory BFGS method for large scale optimization. Mathematical programming, 45, 1-3, 503-528.
Anand Murugan et al. 2012. Statistical inference of the generation probability of T-cell receptors from sequence repertoires. Proceedings of the National Academy of Sciences, 109, 40, 16161-16166.
Jorge Nocedal, Stephen J. Wright. 2006. Conjugate gradient methods. Numerical Optimization, 101134.
Mark Schmidt, Nicolas Le Roux, Francis Bach F. 2013. Minimizing finite sums with the stochastic average gradient. arXiv preprint arXiv:1309.2388.
Stephen J. Wright. 2015. Coordinate descent algorithms. Mathematical Programming, 151, 1, 3-34.
Ivan Zvyagin et al. 2014. Distinctive properties of identical twins' TCR repertoires revealed by high-throughput sequencing. Proceedings of the National Academy of Sciences, 111, 16, 5980-5985.
