In silico предсказание взаимодействия транскрипционного фактора и энхансера как первого этапа регуляции аксонального роста

Методика машинного обучения начала широко применяться в области биомедицины для поиска потенциальных лигандов, в дальнейшем используемых как фармакологические мишени [1]. Классически можно определить энхансеры как цис-действующие последовательности дезоксирибонуклеиновой кислоты (ДНК), функцией которых является увеличение транскрипции (экспрессии) гена. Как правило, они функционируют независимо от ориентации и на различных расстояниях от промотора-мишени (или промоторов) [2]. Одним из необходимых условий для функционирования энхансеров является доступность их сайта связывания транскрипционного фактора (transcription factor-binding site, TFBS), то есть сам энхансер и промотор, как и непосредственно последовательность ДНК какого-либо гена, должны находиться в деконденсированном состоянии [3]. Длина обычной энхансерной последовательности (ЭП) составляет примерно 100–1000 bp (base pairs, пары оснований).

Функционально оценить работу энхансера можно путем использования нескольких метрик, например типа предполагаемого для связывания транскрипционного фактора (ТФ), ориентацию (местоположение) энхансера, его аффинность; порядок, число и расстояние между несколькими TFBS на протяжении всего энхансера и, в конечном счете, лежащую в основе топологию ДНК, в совокупности именуемой «энхансерной архитектурой» [4]. Возможность потенциального управления энхансерами как одними из элементов регуляции экспрессии белка является важной задачей в современной молекулярной биологии.

Центральной проблемой в понимании регуляции генов остается объяснение того, как отбираются специфические наборы генов для экспрессии во время роста клеток, дифференцировки или в ответ на сигналы окружающей среды. Основная цель: определить, как геном фиксированного размера устанавливает огромный набор различных программ развития. Вся информация, необходимая для регуляторных функций, опосредованных энхансерами и промоторами, закодирована в последовательности ДНК с помощью их уникальной комбинации модулей [5]. То есть каждый отдельный модуль связывает один или более ТФ, выполняя таким образом одну из функций всего регуляторного элемента. Кроме того, определенные модули могут служить центральными блоками переключения, соответственно реагируя на входные данные от других модулей того же элемента [5].

Структурные и функциональные исследования ТФ показали, что они представляют собой модульные белки, несущие отдельные области, предназначенные для различных функций: ДНК-связывающий домен, который направляет белок к определенному сайту ДНК, домен мультимеризации, позволяющий собирать гомо- или гетеро-мультимеры, и эффекторный домен, который может модулировать скорость транскрипции (активацию или репрессию) [5]. Модульная природа ТФ вместе с модульной архитектурой энхансеров и промоторов обеспечивает основу для комбинаторного режима экспрессии генов, а неограниченные возможности смешивания и сопоставления энхансерных и белковых модулей позволяют предположить, что может существовать бесконечное количество уникальных программ экспрессии генов, встроенных в геномы относительно ограниченного размера. Подразумевается, что энхансеры и промоторы лежат в конечных точках путей передачи сигнала, которые модифицируют ТФ. Таким образом, определенный ген экспрессируется только в том случае, если клетка воспринимает и соответствующим образом интерпретирует специфические сигналы [5].

Аксональное наведение (axonal guidance) наряду с клеточной миграцией и синаптогенезом – один из ключевых процессов, необходимый для правильного формирования нейронной цепи. Его регулирует широкий ряд сигнальных каскадов, происходящих как в самих нейронах, так и в других клетках, включая нейроглию [6].

С техническим улучшением визуализации процессов аксонального наведения как in vivo, так и in vitro показано, что с развитием нейродегенеративных заболеваний ассоциированы дефекты аксонального транспорта. Эти дефекты обусловлены генетическими мутациями, приводящими к отсутствию связывания моторных белков (кинезин, динеин) либо нарушению функции или развитию нестабильности микротрубочек [6].

Среди регуляторов аксонального наведения можно выделить рецептор DCC (deleted in colorectal cancer / colorectal cancer suppressor, супрессор колоректального рака), который конститутивно экспрессируется на поверхности аксонов. Связывание нетрина-1 с DCC вызывает хемоаттракцию [7]. Сам нетрин-1 секретируется клетками пластинки дна спинного мозга, диффундирует во внеклеточный матрикс и создает градиент, привлекающий растущие комиссуральные аксоны к вентральной срединной линии спинного мозга. Различные направляющие сигнал рецепторы на поверхности конуса роста нейрона постоянно контактируют с окружающей средой, взаимодействуя с соответствующими сигналами, высвобождаемыми клетками-мишенями, тем самым позволяя аксону правильно перемещаться по точной траектории среди многих возможных маршрутов [8]. По некоторым данным, нарушение функции работы DCC-рецептора вместе с нетрином-1 способствует гипоплазии зрительного нерва [9].

В качестве внутриклеточного регулятора аксонального направления и роста можно выделить белок SHH (sonic hedgehog protein, «сверхзвуковой ёжик»). Он играет значимую роль в опосредованном направлении аксонов комиссуральных нейронов в развивающемся спинном мозге через неканонический транскрипционно-независимый путь, взаимодействуя с рецептором Boc и представителем семейства рецепторов, сопряженных с G-белком, – рецептором SMO (Smoothened receptor). Активация киназы семейства Src через SMO под воздействием SHH направляет аксоны комиссуральных нейронов, играя ключевую роль во всем процессе роста и направления [10][11].

Цель работы – предсказание взаимодействия определенных энхансерных областей генов DCC и SHH с аннотированными для них факторами транскрипции.

МАТЕРИАЛЫ И МЕТОДЫ

Подготовка данных

В качестве белков-мишеней выбраны нетриновый рецептор DCC и белок SHH, принимающий опосредованное участие в процессе направленного роста аксонов.

ЭП генов белков DCC и SHH взята из открытой базы данных NCBI (National Center for Biotechnology, Национальный центр биотехнологической информации США) в FASTA-формате^1,2. Для картирования энхансеров использовалась база Ensembl, для отбора потенциальных энхансеров и ТФ к ним – GeneCards.

Структуры факторов транскрипции, а также их ДНК-связывающие домены были взяты из базы данных UniProtKB/Swiss-prot.

Поиск энхансерных последовательностей и предсказание их силы

Для выполнения эксперимента отбирались последовательности, активность которых, по данным Ensembl, наблюдалась в клетках нейрального происхождения; для генов DCC и SHH общей клеточной популяцией стали звездчатые нейроглиальные клетки астроциты. Другими важными критериями отбора стало наличие ЭП и в базе Ensembl, и в базе GeneHencer (представлена на сайте GeneCards), близкое ее расположение относительно самого гена (высокое значение Gene Association Score по GeneHencer) и наличие желаемого гена в качестве аннотированной мишени (например, DCC).

Алгоритмы для оценки силы энхансерной последовательности

Для оценки силы ЭП были выбраны два алгоритма с относительно высокой точностью предсказания, удобные в использовании.

Первым алгоритмом выбран iEnhancer-2L³ [12]. Его входными данными является предполагаемая ЭП в FASTA-формате, а выходными – список подстрок, полученных путем фрагментации исходной последовательности рамкой считывания длиной в 200 bp (модельная длина энхансера) и шагом в 1 bp, а также сила полученной подстроки. Посредством iEnhancer-2L проанализирована сила всей последовательности и в случае, если выходные данные предсказания свидетельствовали о том, что все 200 bp-подстроки являлись сильными, то и саму последовательность считали сильной. При получении разнородных данных: перемежающихся сильных, слабых и не функциональных участков энхансера делалось предположение, что структура последовательности не полностью вовлекается в связывание с факторами транскрипции, а лишь с потенциально сильными. Вычислительные возможности iEnhancer-2L ограничены при анализе очень длинных последовательностей, поэтому ЭП для гена DCC (3399 bp) была разбита на несколько подстрок длиной в 500 bp, эти подстроки проанализированы через iEnhancer-2L, и в качестве выходных данных отобраны те из них, которые являются сильными. Это помогает оптимизировать время, отведенное для докинга данной ЭП с фактором транскрипции.

Вторым алгоритмом для анализа силы ЭП служил ES-ARCNN [13]; его входными данными также является ЭП в FASTA-формате, а выходными – сила всей ЭП. Согласно исследованию, проведенному Т. Zhang и соавт. [13], относительные точности предсказаний iEnhancer-2L и ES-ARCNN составляют 60,5 и 65,5% соответственно.

Молекулярный докинг

Для анализа взаимодействия ТФ с ЭП был выбран метод молекулярного докинга. Это один из современных методов в вычислительной биологии, сутью которого является предсказание оптимального взаиморасположения двух биомолекул, обеспечивающее их стабильное связывание. Сам процесс представляет собой генерацию потенциально возможных конформаций/ориентаций лиганда в сайте связывания белка [14].

В качестве докинг-алгоритма для нашей работы выбран HDOCK⁴. Он отличается от своих типичных аналогов, так как в качестве лиганда принимает последовательность ДНК в FASTA-формате, из которой сам выстраивает структурную модель; в качестве второго субстрата связывания выступает сам ТФ, который загружается в pdb-формате. Из дополнительных настроек были указаны аннотированные на UniProtKB ДНК-связывающие сайты ТФ, если таковые имелись [15]. Интерпретация метрик связывания осуществлялась в соответствии с правилами, указанными на сервере HDOCK. Score (скор, оценка) представляет собой результат расчета оценочных функций ITScorePP или ITScorePR и является собственной разработкой авторов, а confidence score (оценка достоверности) служит эмпирически выведенной оценочной функцией докинга и рассчитывается по формуле:

Важно уточнить, что мы не приравнивали показатель оценки к реальной энергии связывания, поскольку первая не оптимизирована с учетом экспериментальных данных.

Анализ энхансерной последовательности на наличие сайта связывания транскрипционного фактора

Выбранные участки ЭП были проанализированы на наличие в них соответствующих TFBS при помощи модулей Find TFBS with SITECON и Find TFBS with matrices в свободном биоинформатическом программном обеспечении Unipro UGENE⁵ [16]. Данные модули сопоставляют последовательности ДНК с уже имеющейся библиотекой SITECON, а также с позиционными матрицами весов/частот. Отобранные ЭП были загружены в FASTA-формате в Unipro UGENE, далее производился поиск TFBS: для CEBPA (CCAAT/enhancer-binding protein alpha, CCAAT/энхансер-связывающий белок альфа) были использованы модули SITECON⁶ и матрицы весов JASPAR ID: MA0102.2 (технически применялась позиционная матрица частоты – файл в pfm-формате), для NANOG (homeobox protein NANOG, гомеобоксный белок NANOG) использовались только позиционные матрицы весов (файл в pwm-формате), взятые из базы данных HOCOMOCO⁷ (взяты данные с quality A) и адаптированные под структуру, используемую Unipro UGENE [17]. Используемые численные показатели рассчитаны на основе эмпирических данных SELEX (SITECON для CEBPA), а также ChIP-Seq (matrix для CEBPA и NANOG). В качестве выходных данных были получены картированные участки потенциальных TFBS, среди которых учитывались имеющие наивысшую метрику score.

Визуальная интерпретация результатов осуществлялась при помощи PyMol ver. 2.5.4 (Schrödinger, LLC, США).

РЕЗУЛЬТАТЫ

Полученные результаты TFBS представлены в таблице 1.

Использованные модели ТФ, ЭП, ее локализация, длина и предсказанная сила суммированы в таблице 2.

Исследование ЭП с помощью разных алгоритмов показало разнородные результаты: с одной стороны, и iEnhancer-2L, и ES-ARCNN предсказали ЭП гена SHH как сильную, с другой, ЭП гена DCC была оценена iEnhancer-2L как сильная, а ES-ARCNN – как слабая.

Результаты связывания ЭП генов с указанными в таблице 2 факторами транскрипции приведены в таблице 3.

В соответствии с полученными в таблице 3 данными можно утверждать, что связывание ЭП гена DCC с ТФ NANOG на промежутках 1–206 bp и 686–885 bp является наиболее вероятным, связывание ЭП гена SHH с ТФ CEBPA на промежутке 1–500 bp (ограничение HDOCK в 500 bp) является потенциально возможным.

На рисунке представлена визуальная интерпретация результатов.

ОБСУЖДЕНИЕ

По нашим результатам, полученным с помощью комбинации нескольких методов in silico, можно утверждать, что взаимодействие NANOG с ЭП гена DCC и взаимодействие CEBPA с ЭП гена SHH является потенциально возможным и может служить релевантным предметом для будущих исследований. Поскольку в нашей работе предложен оригинальный протокол исследования, а метрики score и confidence score привязаны к конкретному докинг-алгоритму (HDOCK), сравнение полученных данных с результатами других исследователей ограничено. Тем не менее метод молекулярного докинга уже успел зарекомендовать себя в качестве независимого. Так, в исследовании P. Giri и соавт. [18], используя метод стыковки in silico, идентифицированы паттерны (AtMAPK3P) у Arabidopsis thaliana. Из 131 изученного транскрипционного фактора лишь MYB 41 демонстрировал взаимодействие с AtMAPK3P. С использованием поиска минимальной последовательности мотивов, а также с помощью метода докинга также сообщалось о нескольких новых белках, взаимодействующих с MYB, которые необходимо подтвердить киназным анализом in vitro [18]. Данная работа, как и наша, демонстрирует успешное применение биоинформатического подхода на начальных этапах эксперимента в качестве хорошего прогностического моделирования.

В работе отечественных ученых А.М. Андрианова и соавт. [19] методами виртуального скрининга и молекулярного моделирования идентифицированы шесть потенциальных пептидомиметиков кросс-реактивного нейтрализующего анти-ВИЧ-1 (вирус иммунодефицита человека) антитела N6, способных имитировать фармакофорные свойства этого иммуноглобулина путем специфических и эффективных взаимодействий с CD4-связывающим сайтом белка gp120 оболочки вируса ВИЧ [19]. Показано, что ключевую роль в связывании этих соединений с белком gp120 играют ван-дер-ваальсовы взаимодействия с консервативными остатками Phe43-полости гликопротеина, критическими для присоединения ВИЧ-1 к клеточному рецептору CD4, а также водородная связь с остатком Asp-368gp120, образование которой увеличивает химическое сродство без активации нежелательного аллостерического эффекта [19]. На основе полученных результатов авторами был сделан вывод о том, что идентифицированные соединения могут рассматриваться в качестве перспективных кандидатов для проведения детальных экспериментальных исследований с целью их дальнейшего использования в работах по созданию новых противовирусных препаратов. Таким образом, можно утверждать, что in silico исследование А.М. Андрианова и соавт. [19], как и наше исследование, также прошло успешно с точки зрения выявления основных доменов в качестве кандидатов на роль фармакологических молекул.

Также стоит упомянуть, что разработанный нами протокол не является совершенным по нескольким причинам: метод не учитывает конкретные сайты связывания ТФ, т.е. докинг проводился «вслепую»; сам метод межмолекулярного докинга является «статичным» по своей природе, поскольку показывает одну из возможных конформаций ДНК-белкового комплекса (которая вполне может быть неустойчива); значения score и confidence score не являются универсальными.

Указанные ранее недостатки можно исключить включением в пайплайн ab initio метода молекулярной динамики. Молекулярно-динамическое моделирование использует силовые поля, представляющие собой параметрические уравнения с описанием компонентов для различных сил (растяжения связей, Ван-дер-Ваальса и др.), действующих между атомами внутри молекул и между ними [20] и представляет собой вычислительный метод, применяемый для стабильности конформаций и динамики различных биомолекул, включая белки и нуклеиновые кислоты [21].

Преждевременное использование такого сильного метода без предварительной апробации было бы нерациональным, поскольку подобные симуляции вычислительно (задействуются GPU, или графические процессоры) и финансово затратны (аренда сервера для вычислений). На сегодня нашей группой уже ведется изучение нового метода и разработка обновленного пайплайна, который дополнит и, возможно, скорректирует результаты полученные путем докинга.

ЗАКЛЮЧЕНИЕ

In silico предсказание взаимодействия ТФ с ЭП продемонстрировало удовлетворительные результаты и достигло поставленной цели. Несмотря на тот факт, что примененные методики имеют большую неточность и не позволяют делать каких-либо точных выводов, исследование показало, что применение таких методов может служить полезным инструментом для скрининга потенциального взаимодействия ТФ с релевантной последовательностью ДНК (с энхансером, промотором и др.).

ВКЛАД АВТОРОВ

Д.Д. Котельников внес основной вклад в разработку концепции и дизайна исследования, а также подготовку иллюстраций. И.А. Синякин внес вклад в поиск литературных источников, написание и редактирование текста. Е.А. Бородин и Т.А. Баталова руководили процессом статистической обработки, написания и редактирования статьи. Все авторы одобрили окончательный вариант статьи и готовы взять на себя ответственность за все аспекты представленной публикации.

AUTHOR CONTRIBUTIONS

Danil D. Kotelnikov made a major contribution to the development of the concept and design of the study, as well as the preparation of illustrations. Ivan A. Sinyakin contributed to the search for literary sources, writing and editing the text. Evgeny A. Borodin and Tatyana A. Batalova supervised the process of statistical processing, writing and editing the text. All authors approved the final version of the article and are ready to take responsibility for all aspects of the submitted publication.