| Результаты этапа: ===== Задача 1. Разработка методик огрубленного моделирования нуклеосомных фибрилл с учетом процессов активной перестройки хроматина. =====
В ходе выполнения проекта производились работы над библиотекой для моделирования двухцепочечной ДНК PyNAMod. PyNAMod - прикладная библиотека для языка программирования Python. Библиотека представляет собой реализацию дискретной модели гибкой ДНК в комплексе с белками (белки представлены в огрубленной форме). В общих чертах, ДНК в данной модели представлена в виде гибкой цепи из анизотропных бусин, где каждая бусина соответствует паре нуклеотидов. Геометрию такой модели можно описывать как в лабораторной системе координат, так и во внутренних координатах цепи ДНК. В первом случае определяют координаты центров бусин. Во втором случае положение и ориентацию анизотропных бусин задают относительно друг друга: положение следующей пары нуклеотидов задают через 6 параметров Tilt, Roll, Twist, Shift, Rise, Slide [1] (Рисунок 1 А). Такой способ описания цепи ДНК позволяет рассчитывать энергию изгиба ДНК из статистического потенциала от ряда кристаллических структур [2] или продолжительных траекторий молекулярной динамики [3]. Данное координатное пространство нельзя использовать для расчета перекрывания частиц в модели. Положение белков в модели задается относительно “ключевых“ пар нуклеотидов (визуализация модели показана на Рисунке 1 Б). Ядром библиотеки PyNAMod является набор быстрых алгоритмов пересчета геометрии ДНК из лабораторной системы координат во внутреннюю систему координат.
Для расчета конформационных ансамблей в PyNAMod используется метод Монте-Карло. Реализация метода устроена следующим образом: для нуклеотидного шага между нуклеосомами (линкерные участки, участки ДНК, не принадлежащие нуклеосомам) выбирается случайное отклонение параметров шага во внутренних координатах ДНК. После изменения шага рассчитываются новые позиции всех частиц в лабораторных координатах, вычисляется изгибная энергия и энергия парных потенциалов, проводится тест Метрополиса [4]. Один проход по всей нуклеосомной фибрилле называется макроциклом (схема работы алгоритма показана на рисунке 1 В). В ходе одного макроцикла (в данной модели) не происходит активных процессов. Перемещения нуклеосом, отворачивание ДНК и другие процессы проходят между макроциклами, в ходе которых происходит релаксация напряженных конформаций нуклеосомных цепочек.
На втором этапе проекта, в PyNAMod были добавлены новые функции, а также исправлены алгоритмические ошибки. В частности была произведена замена модуля расчета геометрии ДНК по атомистическим структурам, для удаления зависимости ПО от программы 3DNA. Изменены и добавлены новые алгоритмы моделирования перемещения нуклеосом по ДНК и частичного раскручивания нуклеосом. Добавлен алгоритм расчета траекторий протягивания хроматиновых петель. Проведена оптимизация производительности ряда функций. Описание работ согласно пунктам плана приведено ниже.
====Разработка алгоритмов для описания процессов перестроения хроматина: частичного разворачивания нуклеосом и протягивания петель из нуклеосомных фбирилл====
----Моделирование процесса частичного разворачивания нуклеосом----
Несмотря на то, что большинство структур нуклеосом в PDB не содержат открученной ДНК, экспериментально показано, что частичное откручивание возможно в физиологических условиях [5,6]. Руководителем проекта ранее было показано, что такое откручивание значительно изменяет характеристики регулярных нуклеосомных фибрилл [7]. В PyNAMod был добавлен алгоритм для реализации частичного откручивания ДНК от нуклеосом в процессе конформационного поиска. Во время конформационного поиска каждые несколько макроциклов (частота задается исследователем) производится попытка открутить ДНК от нуклеосомы с последующей релаксацией. В случае, если конформация с открученной ДНК принимается за следующий макроцикл, расчет продолжается, в обратном случае откручивание отвергается. Изменения происходящие с фибриллой при изменении степени откручивания ДНК от нуклеосомы показаны на рисунке 2А. Так как нуклеосомы наиболее стабильны в компактном состоянии, для нуклеосом в открученном состоянии характерна большая потенциальная энергия. Для учета этого вклада, в полную энергию системы был дополнительно введен компонент, зависящий от степени отделения ДНК от нуклеосомы. В PyNAMod доабвлено два набора параметров для описания энергии отгиба ДНК в единицах кТ (в тех же единицах, что и энергия изгиба ДНК). Первый набор параметров адаптирован из [8], второй представляет откручивание в виде линейной функции с наклоном 0,28 кТ/п.н. (адаптировано из [9]). Конформационные ансамбли, рассчитанные с данными наборами моделей не показали значимых различий, по этой причине, для дальнейшей работы была выбрана линейная модель раскручивания нуклеосом. При моделировании изменения степени откручивания, важно корректно подбирать число шагов релаксации (число макроциклов моделирования) между попытками открутить ДНК от нуклеосомы (скорости продвижения волны раскручивания нуклеосом). Таким образом, для моделирования продвижения РНК полимеразы (продвижение фронта частично раскрученных нуклеосом), необходимо корректно выбрать число шагов релаксации после откручивания ДНК. Для фибрилл длиной 50 нуклеосом энергия моделей выходит на плато после 20-30 макроциклов моделирования после изменения степени раскручивания (Рисунок 2 Б), что определяет максимальную частоту (скорость) изменения степени откручивания в данной модели.
----Моделирование процесса протягивания петель----
Одним из интересных проявлений активного перестроения хроматина является выпетливание участков генома при помощи комплексов когезина и транскрипционного фактора CTCF. Формирование таких петель может как активировать, так и репрессировать транскрипцию, при этом сам процесс транскрипции может быть ассоциирован с протягиванием петель [10].
Был разработан протокол моделирования протягивания петель из нуклеосомных цепочек при помощи PyNAMod. Для реализации данного протокола в библиотеку была добавлена возможность введения виртуальных частиц, невзаимодействующих с ДНК и нуклеосомами. Данные частицы были выбраны так, чтобы находиться вблизи выхода из каждой нуклеосомы. В систему была добавлена система потенциалов, протягивающих одну часть фибриллы относительно другой (Рисунок 3А, потенциалы прикладывались между виртуальными частичками, потенциал позволяет нечетным связям увеличиваться, а четным уменьшаться, вытягивая петлю). Данный потенциал позволяет петле из нуклеосом на нити ДНК либо затягиваться, либо расслабляться (Рисунок 3В). Была оценена скорость релаксации фибриллы после акта протягивания цепочки. Энергия выходит на плато в среднем после 30-40 макроциклов для цепочек из 50 нуклеосом (Рисунок 3Б). По этой причине для дальнейших расчетов использовалось 50 макроциклов на одно изменение уровня откручивания единичной нуклеосомы в модели или одно изменение состояния протягивания цепи нуклеосом.
===Разработка оптимизированных алгоритмов расчета геометрии фибриллы и создания геометрических выборок===
Одним из важнейших элементов при построении фибрилл является выбор шаблонов для создания нуклеосомных цепочек. Для оценки геометрического разнообразия структур, а также их качества был проведен анализ всех доступных на данный момент структур нуклеосом в PDB. По результатам этого анализа в PyNAMod был загружен набор всех известных на данный момент мононуклеосом с разрешением выше 3 Å и числом н.п. более 140 (всего 161 структура). В ходе работы над курируемым набором структур была опубликована статья о нуклеосомном разнообразии и перспективах развития структурной биологии хроматина (Armeev GA, Gribkova AK and Shaytan AK (2022) Nucleosomes and their complexes in the cryoEM era: Trends and limitations. Front. Mol. Biosci. 9:1070489. doi: 10.3389/fmolb.2022.1070489)
--- Оптимизация работы алгоритмов в PyNAMod ---
В PyNAMod были введены оптимизированные алгоритмы, позволяющие проводить параллельные вычисления конформационного поиска и экспорта конформационных ансамблей. Параллелизация была реализована на двух уровнях: при моделировании активного перестроения фибрилл и для эффективного создания конформационных выборок. При моделировании активного передвижения нуклеосом, после окончания макроцикла конформационного поиска, запускается N (по числу доступных процессорных ядер) процессов релаксации для случайного перемещения N нуклеосом (по одной нуклеосоме на процесс). Для результата каждого из процессов релаксации производится проверка по критерию Метрополиса относительно последнего шага до перестроения. Отбирается один из результатов с наименьшей энергией. Для моделирования эффективного создания конформационных выборок в PyNAMod был реализован аналог “параллельного темперированния” [11]. При старте конформационного поиска запускается N систем при разных температурах (в методе Монте-Карло температура управляет жесткостью критерия Метрополиса). Через заданное число шагов моделирования системы случайным образом попарно обмениваются температурой (вероятность обмена также определяется по критерию Метрополиса). Реализация данного подхода не увеличивает скорость конформационного поиска, однако значительно увеличивает вероятность обнаружения допустимого перемещения нуклеосомы. Такой подход позволяет расширить конформационный поиск для систем в которых не происходит активного изменения топологии, и он используется между актами перемещения, раскручивания нуклеосом, по завершении акта протягивания фибрилл относительно друг друга. Также в PyNAMod была добавлена возможность простого параллельного запуска для удобства набора повторностей эксперимента. Разработанные методы параллелизации показали линейное масштабирование на число процессоров, так как не зависят друг от друга. Все приведенные алгоритмы реализованы при помощи встроенной в Python библиотеки multiprocessing.
Помимо параллелизации алгоритма конформационного поиска был дополнительно оптимизирован алгоритм пересчета пространства координат ДНК в лабораторные координаты, так как это требуется на каждом шаге моделирования. В основном алгоритме для расчета координат каждой пары нужно использовать положение и ориентацию предыдущей пары. Из-за этого алгоритм может работать только последовательно, что значительно увеличивает время его работы. Однако, так как за один шаг моделирования изменения вносят только в одном положении ДНК, изменения для всех последующих шагов будут аналогичными относительно точки сдвига и поворота. Такой подход не идеален: при больших размерах фибрилл, ошибка численной математики приводит к расхождению между значениями стандартного алгоритма и ускоренного. Было показано, что ошибка возрастает с ростом длины фибриллы, однако для фибрилл длиной как минимум 10^6 н.п. остается пренебрежительно малой (Рисунок 4Б). Ошибка также нарастает при внесении последовательных изменений в геометрические параметры, поэтому необходимо периодически рассчитывать значения для реального пространства с помощью стандартного последовательного алгоритма (раз в 100 шагов для фибрилл длиной 1000000 н.п, так как ошибка начинает достигать сотых долей Å) (Рисунок 4В). При условии проверки ошибки каждые 100 шагов для фибрилл длиной 10000 н.п. он позволил увеличить скорость создания выборок конформаций приблизительно в 5 раз (Рисунок 4А).
Дополнительно был изменен алгоритм поиска конфигураций после перемещения нуклеосом. Введен управляющий параметр для ограничения области ДНК подверженной релаксации. Данный подход позволил увеличить долю успешных попыток отбора перемещений нуклеосом приблизительно в 15 раз; хотя такой подход не снижает время общего расчета, он позволяет моделировать больше актов перемещения нуклеосом по цепи ДНК.
Возможности библиотеки, а также подробности реализации были доложены на тематических конференциях, а также во время приглашенного доклада руководителя проекта на летней школе по эпигенетике на биологическом факультете МГУ.
===== Задача 2. Исследование структуры и динамики модельных хроматиновых фибрилл в процессе активной перестройки хроматина. =====
Нуклеосомы в геноме распределены неслучайно. Существуют области как с высокой упорядоченностью нуклеосом, так и участки генома лишенные четко определенной нуклеосомной организации. В работе [12] впервые было показано, что в областях вблизи промоторов S. сerevisiae нуклеосомы точно позиционированы. Такое позиционирование связано с высокой транскрипционной активностью и работой транскрипционных факторов [13]. На данном этапе проекта были исследованы модели нуклеосом на нити ДНК в процессе протягивания петель ДНК с возможностью раскручивания и транслокации нуклеосом.
--- Разработка моделей фибрилл нуклеосом с высоким и низким уровнем транскрипционной активности, сопряженной с протягиванием петель ДНК ---
Для исследования свойств модельных фибрилл были созданы 2 группы фибрилл по 50 нуклеосом, отличающиеся отступами между нуклеосомами.
- Первая группа содержала регулярные или случайно расположенные нуклеосомы таким образом, чтобы средний отступ между ними составлял 160 н.п. (расстояние поддерживаемое ремоделлером CHD1)
- Вторая группа содержала регулярные или случайно расположенные нуклеосомы таким образом, чтобы средний отступ между ними составлял 170 н.п (расстояние поддерживаемое ремоделлером ISW1)
Все фибриллы строились на случайной последовательности ДНК на основании структуры с кодом PDB 3LZ0 (145 н.п. на нуклеосому), а на старте все нуклеосомы были соединены прямыми участками ДНК. Для всех структур были рассчитаны конформационные выборки в режиме протягивания петель (применение потенциала, перемещающего нити фибриллы относительно друг друга показано на рисунке 3 В, отличия петель с регулярным и нерегулярным расположением нуклеосом показаны на рисунке 5АБ). Данные выборки были посчитаны в трех режимах: режиме фиксированных нуклеосом, режиме транслокации нуклеосом и режиме частичного раскручивания нуклеосом перед продвигающейся петлей. Каждая модель была запущена в трех повторностях, и для каждой модели было рассчитано не менее 100000 конформаций. Было показано, что при протягивании петель хроматина средняя энергия деформации ДНК отличается между разными режимами моделирования. Возможность транслокации нуклеосом снижает деформационную энергию фибриллы после протягивания петли в случае регулярной стартовой конфигурации, но в случае нерегулярной конфигурации, возможность транслокации увеличивает энергию модели рисунке 5Г. Такой эффект связан с тем, что петли, сформированные из нерегулярных нитей нуклеосом деформированы, что затрудняет перемещение внутри них. Частичное разворачивание позволяет как регулярным так и нерегулярным фибриллам быстрее релаксировать, так как в этом случае высвобождаются дополнительные гибкие сегменты ДНК (Рисунок5 ВГ). Таким образом было показано, что возможность активного изменения структуры нуклеосом при моделировании процесса протягивания петель позволяет модели обнаруживать более энергетически выгодные состояния.
--- сравнение данных моделирования цепочек нуклеосом с участками генома дрожжей ---
Полученные модели фибрилл с разрешенной транслокацией нуклеосом были сравнены с участками генов дрожжей с высокой и низкой транскрипционной активностью (аннотация уровнем экспрессии была взята из эксперимента [14], позиционирование нуклеосом из работы [13] идентификатор эксперимента в SRA - SRR1802189). Размер гена в геноме дрожжей в среднем составляет 2000 н.п и содержит порядка 10 нуклеосом. По этой причине было проведено сравнение для центральных нуклеосом в петеле из фибрилл длиной 50 нуклеосом. В геноме дрожжей присутствует достаточно мало участков, для которых характерна строгая периодичность позиционирования нуклеосом. Даже при моделировании транслокации нуклеосом, профили заселенности для моделей плохо коррелируют с экспериментальными профилями для разных участков генома. Однако, так как результатом экспериментов по определению позиций нуклеосом является усреднение по популяции клеток, в результирующем экспериментальном профиле содержат фибриллы с разной степенью отступа. Например, для гена TAH11 (выбран как один из наиболее покрытых в эксперименте), была показана хорошая корреляция (0,63) между экспериментальным профилем заселённости нуклеосомами и комбинацией активных моделей нуклеосомных фибрилл в соотношении 30% фибрилл с отступом 160 и 70% фибрилл с отступом 170 (Рисунок 6А). Такие отступы поддерживаются ремоделлерами CHD1 и ISW1 соответственно. Интересно, что последний является репрессором транскрипции и ген TAH11 аннотирован как мало экспрессирующийся.
Также были построены все возможные комбинаторные модели позиционирования нуклеосом на 6 разных генах дрожжей (см следующий пункт отчета, Рисунок 6 Б), отличающихся по уровню экспрессии, однако зависимости между параметрами моделей нуклеосомных цепочек (число доступных конформаций, радиус инерции) и уровнем экспрессии гена обнаружено не было (Рисунок 6 ВГ). Данный результат может являться следствием того, что эксперимент по определению уровня экспрессии и эксперимент по картированию нуклеосом были проведены в разных условиях и на разных клеточных культурах.
--- Исследование конформационного пространства всех возможных пентануклеосом для уточнения позиционирования нуклеосом по экспериментальным данным ---
Эксперименты по определению позиций нуклеосом по большей части получены для культур клеток и содержат сигналы позиционирования для большой совокупности геномов. На итоговых профилях заселенности мы видим результат суперпозиции множества альтернативных вариантов позиционирования, который дополнительно смазан неоднородностью длин сигналов. Определение возможных расположений нуклеосом по таким данным - сложная вычислительная задача. В ходе выполнения проекта был разработан метод отбора возможных позиций нуклеосом по данным MNase секвенирования. Метод основан на том, что не все взаимные расположения нуклеосом допустимы, так как могут приводить к стерическим перекрываниям.
Для поиска доступных комбинаций позиций нуклеосом были созданы модели всех возможных комбинаций цепочек из 5 нуклеосом и длиной соединяющей ДНК от 0 до 100 н.п. (отступы от 145 до 245). Всего было создано 10^8 комбинаций нуклеосом. Расчеты проводились в параллельном режиме с использованием 10 узлов вычислительного кластера, всего 80 процессорных ядер. Ряд комбинаций линкерных областей приводит к появлению заслоненных структур, а на тепловой карте (для первых трех нуклеосом) видна периодичность рисунке 7 А. Для всех комбинаций фибрилл из 5 нуклеосом чаще всего встречаются линкеры длиной 10N+5 (рисунок 7 В), что хорошо согласуется с литературными данными для регулярных цепочек [15]. Было показано, что с ростом длины нуклеосомных цепочек, доля допустимых комбинаций длин линкерных областей значительно падает (рисунок 7 Б). Таким образом, для плотных цепочек нуклеосом доступно значительно меньшее конфигурационное пространство чем для цепочек с большим межнуклеосомным расстоянием. Опираясь на данное наблюдение, был создан алгоритм отбора и фильтрации нуклеосомных цепочек произвольной длины по экспериментальным данным. Алгоритм позволяет уточнять позиции нуклеосом на генах, а также отбрасывать сигналы нуклеосом, несовместимые с остальными сигналами. Алгоритм основан на представлении кандидатов позиций нуклеосом в виде графа с последующим поиском путей в графе от начала гена до конца с отсевом недопустимых взаимных ориентаций из 5 нуклеосом. Алгоритм позволяет уточнять позиции нуклеосом по экспериментальным картам позиционирования низкого разрешения, а также по сигналам от одиночных нуклеосом (Рисунок 6 Б). Показано, что неравномерности малого масштаба в профилях заселенности нуклеосомами воспроизводятся при расчете профилей от отобранных моделей. Данный подход показал эффективность на наборе генов для дрожжей. По результатам моделирования всех цепочек из 5 нуклеосом был создан профиль вероятности обнаружения соседней нуклеосомы, относительно заданной (рисунок 7 Г). Данный профиль можно применять в качестве оконного фильтра для данных о позиционировании нуклеосом низкого разрешения. Такой подход позволяет предположить наиболее вероятные положения нуклеосом на участках генома с учетом контекста соседей, особенно в случае появления широких пиков заселенности нуклеосомами (рисунок 7 Д). Результаты данной работы приняты к публикации в журнале Биофизика (В.А. Васильев, Д.М. Рябов, А.К. Шайтан, Г.А. Армеев, Уточнение позиций нуклеосом внутри отдельных генов с использованием методов молекулярного моделирования и данных MNase секвенирования. Биофизика. 2023). |
