Методы вероятностной кластеризации для интерпретации матричной и графовой информации в разнородных базах данных тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Усольцев Дмитрий Андреевич

• Усольцев Дмитрий Андреевич
• кандидат наук
• 2025

Оглавление диссертации кандидат наук Усольцев Дмитрий Андреевич

Реферат

Synopsis

Глава 1. ОБЗОР ПРЕДМЕТНОЙ ОБЛАСТИ

1.1. Задача кластеризации в общем виде

1.2. Методы вероятностной кластеризации

1.2.1. Метод cмеси Гауссовых распределений

1.2.2. Байесовские смеси (Bayesian Mixture Models)

1.3. Методы графовой кластеризации

1.3.1. Теория графов

1.3.2. Поиск подграфа максимального веса

1.3.3. Вероятностный метод графовой кластеризации

1.4. Методы прочтения генома

1.5. Методы для разработки баз данных с разнородной информацией

1.5.1. Метод уменьшения размерности матриц

1.5.2. Метод выбора подмножества признаков при наличии корреляции между ними

1.5.3. Равновесие Харди-Вайнберга

1.5.4. Метод импутации генетических данных

1.5.5. Метод определения генетического родства

1.5.6. Метод максимального правдоподобия

1.6. Задача линейного объединения статистик

1.7. Основные международные базы данных с разнородной информацией

1.8. Актуальность темы в контексте биомедицины в Российской Федерации

1.9. Анализ выживаемости как дополнительный статистический источник данных

1.10. Задачи, решаемые в диссертационной работе

Выводы по главе

Глава 2. ВЕРОЯТНОСТНЫЙ МЕТОД КЛАСТЕРИЗАЦИИ МАТРИЧНЫХ ДАННЫХ

2.1. Формулировка задачи

2.2. Описание предлагаемого метода

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

2.4. Применение предлагаемого метода к медико-генетической базе данных российской популяции (МГБД-Р)

Выводы по главе

Глава 3. ВЕРОЯТНОСТНЫЙ МЕТОД ГРАФОВОЙ КЛАСТЕРИЗАЦИИ

3.1. Формулировка задачи

3.2. Расчёт многокомпонентных весов вершин

3.3. Описание предлагаемого метода графовой кластеризации

3.4. Эффективность предлагаемого метода

3.4.1. Эффективность предлагаемого метода для активных модулей, генерируемых поиском в ширину

3.4.2. Эффективность предлагаемого метода для активных модулей, генерируемых случайным поиском

3.5. Эффективность предлагаемого метода на реальных данных выживаемости пациентов

Выводы по главе

Глава 4. МЕТОДИКА АНАЛИЗА ДАННЫХ В РАЗНОРОДНЫХ БАЗАХ

4.1. Описание методики

4.2. Разработка первого модуля МГБД-Р (МГБД-Р-1)

4.2.1. Необходимость создания баз данных МГБД-Р

4.2.2. Сбор фенотипических данных для первого модуля МГБД-Р

4.2.3. Генерация и импутация генетических данных для первого модуля МГБД-Р

4.3. Генетическая структура российской популяции

4.3.1. Анализ генетической структуры российской популяции

4.3.2. Уникальные финские и российские ДНК-варианты

4.4. Разработка второго модуля МГБД-Р (МГБД-Р-2)

4.4.1. Проведение полногеномных ассоциативных исследований

4.4.2. Агрегация результатов полногеномных ассоциативных исследований МГБД-Р-2

4.5. Применение разработанного вероятностного метода матричной

кластеризации к МГБД-Р-1

4.6. Применение разработанного вероятностного метода графовой кластеризации ко второму модулю зарубежной МГБД

4.7. Применение разработанного вероятностного метода графовой кластеризации к МГБД-Р-2

Выводы по главе

Заключение

Список литературы

Список рисунков

ПРИЛОЖЕНИЕ А. Результаты интеллектуальной деятельности

ПРИЛОЖЕНИЕ Б. Акты внедрения

ПРИЛОЖЕНИЕ В. Тексты публикаций

Реферат

Общая характеристика диссертации

Актуальность темы исследования определяется экспоненциальным ростом объёмов и разнообразием современных данных, для анализа которых широко применяются методы кластеризации. Эти методы позволяют выявлять скрытые структуры и закономерности, однако классические детерминированные подходы (k-means [1], иерархическая кластеризация [2], DBSCAN [3]) обладают существенными ограничениями. В частности, они не учитывают неопределённость при отнесении объектов к кластерам, что приводит к потенциальным ошибкам и снижению надёжности полученных результатов.

Альтернативой являются вероятностные методы кластеризации, которые основываются на моделировании эмпирических распределений вероятности принадлежности объектов к кластерам [4]. Используя композиции параметризованных распределений, такие подходы явно учитывают априорные знания о признаках и позволяют количественно оценивать степень неопределённости. Благодаря этому становится возможным получение статистических оценок надёжности кластеров, существенно повышающих точность и интерпретируемость результатов анализа.

Для матричных данных такие вероятностные подходы особенно эффективны в условиях высокой размерности, неполноты и зашумлённости. Они обеспечивают устойчивость к выбросам, позволяют корректно формировать эталонные выборки с учетом эмпирических распределений вероятности принадлежности объектов к кластерам и выявлять скрытые группы объектов, основанные на сложных многомерных взаимосвязях между признаками [5].

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

Современные информационные системы часто представляют собой гибридные хранилища, объединяющие табличные, графовые, временные и текстовые форматы данных. Разработка вероятностных методов кластеризации матричных и графовых данных позволяет проводить согласованный анализ подобных разнородных массивов информации. В

результате становится возможным корректное формирование несмещенных наборов данных, выявление скрытых закономерностей и снижение вероятности систематических ошибок.

Примером таких баз данных с разнородной информацией являются медико-генетические базы данных (МГБД) [7], объединяющие геномные матрицы, клинические данные и графы биомолекулярных взаимодействий. Важным аспектом анализа данных МГБД является полногеномный поиск ассоциаций (Genome-wide association study, GWAS) [8]. Вероятностные алгоритмы кластеризации уже были использованы для проведения GWAS без передачи персональных данных [9]. Однако до сих пор важной проблемой остаётся формирование несмещённого набора данных с учётом клинических характеристик для проведения полногеномного ассоциативного анализа. Использование вероятностных подходов матричной кластеризации может позволить эффективно решить эту проблему, моделируя априорные распределения признаков и формируя корректные выборки для анализа.

Вероятностная кластеризация, в свою очередь, актуальна для выявления так называемых функциональных подграфов (активных модулей) [10], отражающих коллективную активность групп объектов (например, генов). Совмещение топологических характеристик графов с данными GWAS, анализа дифференциальной экспрессии или анализа выживаемости пациентов делает эти методы крайне востребованными в интерпретации сложных многокомпонентных данных. Под интерпретацией понимается аннотация данных - присвоение объектам стандартизированных функциональных меток на основе внешних знаний, включая функциональную аннотацию результатов GWAS.

Таким образом, разработка и повышение качества вероятностных методов кластеризации матричных и графовых данных является важным направлением современной науки. Они позволяют учитывать априорные знания, корректно оценивать неопределённость и обеспечивать высокую точность и интерпретируемость результатов анализа, что имеет широкое применение в различных областях науки и техники.

Степень разработанности темы исследования

В области интеграции разнородных данных накоплен значительный объём методов и алгоритмов, направленных на их интерпретацию (аннотацию) и структурный анализ. Активно развиваются различные специализированные базы, включая МГБД, создаваемые в разных странах [11]. Эти ресурсы аккумулируют матричные и графовые данные из множества источников.

В контексте методов кластеризации графовых данных следует отметить подходы на основе метода Монте-Карло по схеме марковских цепей (Markov Chain Monte Carlo - MCMC) [6] и поиска подграфа максимального веса (Maximum Weight Connected Subgraph - MWCS) [12], применявшиеся для анализа и визуализации крупных сетей. При этом существенным направлением остается повышение качества кластеризации при извлечении функциональных подграфов за счет использования априорной информации.

В диссертации [13] был предложен метод графовой кластеризации для совместного анализа результатов генотипирования и экспрессии генов, однако он ограничивался простыми весами вершин (учитывал информацию лишь из одного источника). Кроме того, данный метод базировался на детерминированном поиске подграфа максимального веса. Между тем интеграция данных из разных источников в многокомпонентные веса вершин может повысить точность выявления функциональных подграфов, как и использование вероятностных методов в сочетании с методами машинного обучения для учета априорных знаний.

Целью данной работы является повышение качества кластеризации в задачах матричного и графового анализа разнородных баз данных за счет использования априорных знаний.

Для достижения поставленной цели в диссертации решаются следующие задачи:

1. Обзор современных методов вероятностной кластеризации матричных и графовых данных, а также статистических методов, применяемых для разработки медико-генетических баз данных с разнородной информацией.

2. Разработка вероятностного метода матричной кластеризации на основе априорного распределения признаков для формирования несмещенного набора данных.

3. Разработка вероятностного метода графовой кластеризации на основе методов Монте-Карло по схеме марковских цепей и априорной информации, интегрируемой с помощью машинного обучения, для выявления кластеров, обладающих целевыми признаками.

4. Апробация и внедрение результатов исследования на основе методики, использующей предложенные вероятностные методы и машинное обучение для учёта априорных знаний, с целью

повышения качества кластеризации при анализе баз данных с разнородной информацией.

Теоретическая значимость работы заключается в:

- объединении методов матричного разложения и вероятностного моделирования с учётом априорных распределений, применимой к широкому классу многомерных данных;

- объединении априорных знаний за счет методов машинного обучения и методов Монте-Карло по схеме марковских цепей для задачи кластеризации в графах с взвешенными вершинами. Также показана применимость данного подхода к широкому спектру данных, что значительно повышает потенциал для практического использования в различных отраслях анализа данных.

Практическую значимость работы определяют:

- организация медико-генетической информации в базу данных (МГБД-Р), состоящую из модуля для хранения персональной медико-генетической информации и модуля для хранения и передачи статистической медико-генетической информации в рамках распределённой системы, доступной по адресу https://biobankrus.almazovcentre.ru/ и обеспечивающей надёжный и конфиденциальный обмен информацией;

- разработанная МГБД-Р внедрена в генетические исследования факторов заболеваний, проводимых на базе ФГБУ НМИЦ им. В. А. Алмазова Минздрава России и применяется при метаанализе медико-генетических данных крупнейшими научными консорциумами (Global Biobank Metaanalysis Initiative - GBMI, The Global Lipids Genetics Consortium -GLGC);

- повышение качества кластеризации предложенным методом вероятностной кластеризации графовых данных до 23 % на симулированных данных и на 13 % при анализе реальных данных, в том числе из разработанной МГБД-Р;

- реализация методов вероятностной кластеризации в виде программы на языке R: https://github.eom/11 Dmitriyll/mgc, а также их внедрение в образовательные процессы Университета ИТМО и ФГБУ НМИЦ им. В. А. Алмазова Минздрава России.

Положения, выносимые на защиту, обладающие научной новизной, а также имеющие важное практическое значение:

1. Вероятностный метод кластеризации для формирования набора данных, использующий сингулярное разложение матриц, отличающийся тем, что с целью снижения вероятности ложноположительных результатов при сравнении двух наборов данных, в нём применяются расстояние Махаланобиса и априорное распределение признаков.

2. Вероятностный метод графовой кластеризации, использующий метод Монте-Карло по схеме марковских цепей, отличающийся тем, что с целью повышения точности кластеризации, в нем используется априорная информация, учитываемая с помощью машинного обучения, и несколько источников априорных знаний.

3. Методика анализа данных в разнородных базах, основанная на сравнительном анализе и последующей графовой кластеризации его результатов, отличающаяся тем, что с целью повышения качества кластеризации, она включает предложенные вероятностные методы матричной и графовой кластеризации.

Научной новизной обладают первые два положения, а третье - имеет важное практическое значение.

Методы исследования. В работе используются методы анализа графов, методы теории вероятности и математической статистики, методы линейной алгебры и методы машинного обучения.

Достоверность научных положений, выводов и практических рекомендаций находят отражение в публикациях в рецензируемых журналах и конференциях, в том числе из базы Scopus и Web of Science, включая публикации результатов в высокорейтинговых журналах, например, Nature Communications (IF = 14.7).

Соответствие паспорту специальности. Диссертация соответствует следующему пункту паспорта специальности 2.3.8. «Информатика и информационные процессы (технические науки)»:

Тринадцатый пункт паспорта специальности: «Разработка и применение методов распознавания образов, кластерного анализа, нейро-сетевых и нечетких технологий, решающих правил, мягких вычислений при анализе разнородной информации в базах данных».

Предложен вероятностный метод матричной кластеризации с априорным распределением признаков, сокращающий систематические смещения при формировании эталонных выборок.

Разработан вероятностный метод графовой кластеризации с использованием методов Монте-Карло по схеме марковских цепей и машинного обучения. Разработанный метод был применен автором к различным сравнительным аналитическим исследованиям, в том числе из второго модуля разработанной базы МГБД-Р, и показал более высокое качество кластеризации в сравнении с аналогами.

Апробация результатов работы

Результаты исследований были представлены в качестве выступлений на следующих международных и всероссийских конференциях:

- XIV Конгресс молодых ученых, 2025, Университет ИТМО, Санкт-Петербург, Россия // Биобанк России - онлайн-платформа для хранения и анализа результатов GWAS-исследований полигенных заболеваний, а также разработка фенотипических предсказательных моделей на базе когорты Биобанка России;

- American Society of Human Genetics, 2024, Denver, CO, USA // A meta-predictor for causal gene identification in GWAS overcomes limitations of existing computational approaches;

- Cold Spring Harbor Laboratory Conference, 2024, Laurel Hollow, NY, USA // Understanding complex trait susceptibilities and ethnic diversity in a sample of 4,145 Russians through analysis of clinical and genetic data;

- American Society of Human Genetics, 2023, Washington, D.C., USA // Understanding complex trait susceptibilities and ethnical diversity in a sample of 4,145 Russians through analysis of clinical and genetic data;

- European Society of Human Genetics, 2021, online // ZNF831 R1393Q is a UV-induced mutational hotspot in cutaneous melanoma significantly affecting patient survival.

Награды

Результатам исследований была присуждена следующая награда:

- Reviewers' Choice Abstract (10 % лучших работ на American Society of Human Genetics, 2024, Denver, CO, USA).

Постерный доклад: A meta-predictor for causal gene identification in GWAS overcomes limitations of existing computational approaches

Публикации

Основные результаты, изложенные в настоящей диссертации, опубликованы в шести работах, индексируемых в базе Scopus, три из которых опубликованы в международных изданиях и также содержатся в базе Web of Science, включая публикации результатов в высокорейтинговых журналах, например, Nature Communications (IF = 14.7).

1. Усольцев Д.А. Вероятностный метод матричной кластеризации с априорным распределением признаков для формирования несмещенной контрольной группы // Научно-технический вестник информационных технологий, механики и оптики. 2025. Т. 25, № 5. С

2. Усольцев Д.А., Молотков И.И., Артемов Н.Н., Сергушичев А.А., Шалыто А.А. Применение марковских цепей Монте-Карло и машинного обучения для нахождения активного модуля в биологических графах // Научно-технический вестник информационных технологий, механики и оптики. 2024. Т. 24, № 6. С

3. Усольцев Д.А., Молотков И.И., Артемов Н.Н., Сергушичев А.А., Шалыто А.А. Метод определения активного модуля в биологических графах с многокомпонентными весами вершин // Научно-технический вестник информационных технологий, механики и оптики. 2025. Т. 25. № 3. С

4. Usoltsev D., Kolosov N., Rotar O., Sergushichev A.A., Artomov N.N. et al. Complex trait susceptibilities and population diversity in a sample of 4,145 Russians // Nature Communications. 2024. Т. 15, № 6212. Open access (первый автор).

5. Usoltsev D., Tolkunova K., Moguchaia E. et al. Transgenerational and intergenerational effects of early childhood famine exposure in the cohort of offspring of Leningrad Siege survivors // Scientific Reports. 2023. Т. 13, № 11188. Open access.

6. Zlotina A., Barashkova S., Zhuk S., Skitchenko R., Usoltsev D., Sokolnikova P., Artomov M., Alekseenko S., Simanova T., Goloborodko M., Berleva O., Kostareva A. Characterization of pathogenic genetic variants in Russian patients with primary ciliary dyskinesia using gene panel sequencing and transcript analysis // Orphanet Journal of Rare Diseases. 2024. Т. 19, № 310. Open access.

Личный вклад автора

1. В публикации [1] - разработка вероятностного метода матричной кластеризации (100 %).

2. В публикации [2] - разработка вероятностного метода графовой кластеризации, учитывающего априорные знания (80 %), остальные авторы - консультирование при интерпретации результатов (20 %).

3. В публикации [3] - разработка и апробация вероятностного метода графовой кластеризации, учитывающего многокомпонентные вершины графа и априорные знания (80 %), остальные авторы -консультирование при интерпретации результатов (20 %).

4. В публикации [4] - разработка МГБД-Р, проведение GWAS в российской популяции, выявление генетических особенностей российской популяции, доказательство преимуществ проведения GWAS в российской популяции (60 %); Kolosov N. - импутация генетических данных (5 %); Rotar O. - рекрутинг пациентов и сбор клинических данных (15 %); Artomov N. - дизайн исследования и генерация генетических данных (15%), остальные соавторы -помощь в сборе клинических данных, рекрутинге пациентов и интерпретации результатов (5 %).

5. В публикации [5] - апробация вероятностного метода кластеризации матричных данных для получения несмещенного набора данных для сравнительного анализа, (45 %). Tolkunova K. - рекрутинг пациентов и сбор клинических данных (45 %), остальные соавторы -консультирование при интерпретации результатов (10 %).

6. В публикации [6] - разработка референсной базы данных аллельных частот мутаций в российской популяции (5%); остальные соавторы -разработка дизайна исследования, сбор данных и интерпретация результатов (95 %).

Внедрение результатов работы

1. Работа по созданию базы данных МГБД-Р и анализу генетических особенностей Российской популяции выполнена в рамках гранта Министерства науки и образования Российской Федерации (номер 075-15-2022-301).

2. Результаты диссертации использовались в образовательном процессе в рамках непрерывного медицинского образования на базе НМИЦ им. В. А. Алмазова при чтении лекций и проведении практических занятий, а также в рамках подготовки магистров по программе «Bioinformatics and System Biology» в Университете ИТМО.

3. МГБД-Р-1 используется в эпидемиологических и генетических исследованиях факторов заболеваний в российской популяции, проводимых на базе НМИЦ им. В. А. Алмазова.

4. МГБД-Р-2 используется крупнейшими научными консорциумами (Global Biobank Meta-analysis Initiative - GBMI, The Global Lipids Genetics Consortium - GLGC) при метаанализе генетических факторов заболеваний.

Регистрация программ

Являюсь соавтором двух свидетельств о регистрации программ для ЭВМ на основе разработанной МГБД-Р:

1. Ротарь О. П., Бабенко А. Ю., Колесова Е. П., Бояринова М. А., Могучая Е. В., Ерина А. М., Толкунова К. М., Кван В. В., Усольцев Д. А., Конради А. О., Шляхто Е. В. Электронный калькулятор риска диабета второго типа // Свидетельство о государственной регистрации программы для ЭВМ. № 2024 690421 от

2. Ротарь О. П., Бабенко А. Ю., Колесова Е. П., Бояринова М. А., Могучая Е. В., Ерина А. М., Толкунова К. М., Кван В. В., Усольцев Д. А., Конради А. О., Шляхто Е. В. Электронный калькулятор риска сердечнососудистых инцидентов // Свидетельство о государственной регистрации программы для ЭВМ. № 2024 691771 от

Введение диссертации (часть автореферата) на тему «Методы вероятностной кластеризации для интерпретации матричной и графовой информации в разнородных базах данных»

Объем и структура работы

Диссертация состоит из введения, четырех глав, заключения и трех приложений. Полный объём диссертации составляет 234 страницы, включая 16 рисунков. Список литературы содержит 113 наименований.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

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

В первой главе приведен обзор литературы, охватывающий базовые понятия и методы, используемые в диссертации, включая методы вероятностной кластеризации матричных и графовых данных, а также статистические методы из прикладной области - популяционной генетики, необходимые для разработки медико-генетической базы данных. Даются определения ключевым терминам, включая понятия «ген» и «ДНК-вариант».

В разделе 1.1. описана задача кластеризации в общем виде. Кластеризация представляет собой задачу разбиения множества объектов на непересекающиеся подмножества так, чтобы объекты внутри одного кластера были близки друг к другу по заданной метрике, а объекты из разных кластеров существенно отличались. В классических алгоритмах, таких как k-means отсутствует предположение о вероятностном распределении данных.

В разделе 1.2. описаны вероятностные методы кластеризации матричных данных. Эти методы отличаются от классических тем, что в них явным образом моделируются процесс порождения данных с помощью вероятностный распределений.

Наиболее распространённым примером такого подхода являются модели смесей гауссовых распределений (Gaussian Mixture Model, GMM) [14], в которых данные аппроксимируются суммой K многомерных нормальных распределений, каждое из которых описывает свой кластер: весовые коэффициенты nk задают долю компоненты в смеси, а параметры ^k и Dk определяют положение и форму кластеров. Оценка этих параметров сводится к максимизации логарифма правдоподобия, что эффективно решается алгоритмом Expectation-Maximization (EM-алгоритмом) [15]: на E-шаге вычисляются апостериорные вероятности принадлежности каждого объекта к каждому кластеру, а на M-шаге обновляются оценки параметров смеси.

Байесовская модель гауссовской смеси (Bayesian Gaussian Mixture Model) [16] представляет собой вероятностное обобщение классической GMM, в котором на параметры модели (веса кластеров, средние значения и ковариационные матрицы компонент) накладываются априорные распределения. Однако точное вычисление апостериорного распределения в данной модели оказывается вычислительно неразрешимым из-за экспоненциального числа скрытых конфигураций и невозможности аналитического вычисления многомерных интегралов. В связи с этим на практике применяются методы приближённого байесовского вывода, в частности - метод Монте-Карло по схеме марковских цепей (Markov Chain Monte Carlo, MCMC) [17]. К наиболее распространённым алгоритмам, реализующим методы на основе MCMC, относятся: Gibbs sampling [18], позволяющий поочерёдно сэмплировать параметры из их условных распределений, и более общий алгоритм Metropolis-Hastings [19], в котором предлагаемые значения параметров в зависимости от отношения плотностей принимаются или отклоняются с определённой вероятностью. Применение этих алгоритмов позволяет получать оценки параметров модели и вероятностей принадлежности объектов к кластерам, что особенно актуально в задачах анализа неоднородных данных.

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

Граф определяется как пара G = (V, E), где V - множество вершин, а E -множество рёбер, соединяющих эти вершины. В неориентированных графах рёбра задают симметричные связи, тогда как в ориентированных графах каждое ребро имеет направление. Взвешенные графы дополнительно содержат числовые значения, назначенные рёбрам или вершинам.

Графовая кластеризация сводится к поиску в исходном графе подграфа (активного модуля), обладающего целевым признаком с учётом априорных данных. Вершинам графа присваиваются p-значения, которые моделируются смесью равномерного распределения U(0,1) и бета-распределения B(a,1) с уровнем шума X [20]. По аппроксимации этой смеси каждому p-значению ставится в соответствие вес, отражающий значимость вершины по отношению к искомому признаку.

Существуют два подхода к поиску активного модуля. В детерминированном варианте задача формулируется как поиск связного подграфа максимального суммарного веса [12]. В вероятностном подходе применяют MCMC-семплирование (алгоритм Метрополиса-Гастингса) [21]: начиная со случайного подграфа фиксированного размера, поочерёдно удаляют и добавляют вершины, принимая новые состояния с вероятностью, зависящей от отношения апостериорных шансов и предложенных переходных вероятностей [6]. После серии итераций для каждой вершины оценивают долю попаданий в выборку подграфов, что и даёт её эмпирическую вероятность принадлежности к активному модулю.

В разделе 1.4 описаны основные определения ключевых терминов популяционной генетики. Рассмотрены основные методы прочтения генома. Рассмотрена классификация мутаций.

Геном человека представляет собой упорядоченную последовательность ДНК, где каждая буква четырехбуквенного алфавита - {T, C, G, A}, соответствует одному нуклеотиду. Для представления генома часто используется структура типа словаря, где ключом является пара (chr, pos) -номер хромосомы и позиция нуклеотида (локус), а значением - сам нуклеотид. Поскольку человеческий геном содержит 23 пары хромосом, каждому локусу соответствует пара нуклеотидов (по одному нуклеотиду от каждой из парных хромосом), образующая генотип. Сравнение полученной последовательности с эталонной (референсной) осуществляется посредством выравнивания [22], в ходе которого выявляются ДНК-варианты, наиболее распространёнными из которых являются однонуклеотидные замены (single nucleotide polymorphism, SNP).

Для прочтения генома применяют два основных метода: SNP-генотипирование [23], позволяющее анализировать заранее известные локусы с альтернативным аллелем (ДНК-вариантом), встречающимся более чем у 1 % популяции, и секвенирование [24], охватывающее весь геном или только экзом [25] (~1 % генома, кодирующий РНК и белки). При использовании SNP-генотипирования, благодаря явлению неравновесной сцепленности [26], можно восстанавливать недостающие генотипы методом импутации [27].

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

В разделе 1.5 описаны методы, необходимые для разработки медико-генетической базы данных (МГБД).

Под интерпретацией информации в медико-генетической базе данных понимается аннотирование этой информации.

Полногеномный ассоциативный анализ [8] оценивает связь между числом альтернативных аллелей SNP и фенотипом (признаком, болезнью) с использованием регрессионных моделей (линейной для количественных и логистической для бинарных признаков). Предполагается аддитивный эффект, при котором вклад двух альтернативных аллелей равен удвоенному эффекту одной. Значимость ассоциаций определяется через z-статистику с учетом корректировки на множественные гипотезы (например, методом Бонферрони) [29]. Для избежания ложноположительных результатов учитывается популяционная структура выборки для GWAS с помощью методов уменьшения размерности матриц, таких как анализ главных компонент (principal component analysis, PCA) [30].

Пусть дана матрица X~G Mnxm с числом индивидов n и числом ДНК-вариантов m (закодированных как 0, 1, 2 - число альтернативных аллелей SNP у конкретного индивида). PCA можно вычислить через ковариационную матрицу, тогда собственные векторы этой матрицы задают направления, вдоль которых дисперсия данных максимальна. Собственные значения указывают насколько велика дисперсия вдоль соответствующего собственного вектора. Альтернативно PCA можно вычислить напрямую методом сингулярного разложения (singular value decomposition, SVD) [31]. Так как существует неравновесная сцепленность вариантов внутри локуса, перед проведением PCA необходимо отобрать независимые ДНК-варианты из каждого локуса, что значительно сократит размер матрицы X и ускорит вычисления без потери информации.

Пусть имеется конечное множество S = {б1, S2, ..., Sn} ДНК-вариантов и задана мера связи между любыми двумя вариантами (X, У) через коэффициент детерминации. Необходимо выбрать подмножество в котором любые

два варианта слабо коррелированы - г2 < t для некоторого заданного порога 1 Этот подход, известный как LD-pгuning [32], позволяет исключить сильно скоррелированные варианты и упростить дальнейший анализ. Также могут быть исключены ДНК-варианты, которые вероятно являются техническими ошибками. Для этого применяется информация о закономерностях динамики аллельных частот в популяции, в частности равновесие Харди-Вайнберга [33].

Если в популяции аллели А и а встречаются с частотами р и q (p+q = 1), тогда генотипные частоты АА, Аа, аа составляют соответственно р2, 2pq, q2. Это соотношение известно как равновесие Харди-Вайнберга и сохраняется при соблюдении ряда идеализированных условий. На практике проверка отклонений от этого равновесия позволяет выявлять ДНК-варианты с некорректными генотипными частотами и повышать качество генотипированных данных. После удаления технических ошибок, возможно, импутировать [34] недостающие генетические варианты.

Пусть имеется матрица G (размером М^) с генотипами gij и референсная панель Н с полногеномными данными. Задача импутации заключается в восстановлении недостающих генотипов ^ путем максимизации апостериорной вероятности Р ^ | вхц), Н), где вуц) - все прочие наблюдаемые генотипы, исключая gij. Для этого используется скрытая марковская модель, параметры которой оцениваются с использованием алгоритма «прямого-обратного хода», который учитывает неравновесную сцепленность между локусами [35].

Для проведения полногеномного ассоциативного исследования необходимо, чтобы данные о фенотипах и генотипах были доступны на одном сервере. Однако медицинская информация часто распределена между различными больницами и хранится локально, а генетические данные требуют значительных вычислительных и хранилищных ресурсов из-за своего объёма. Объединение медицинской и генетической информации в единую структуру с использованием унифицированных идентификаторов позволяет создать первый модуль медико-генетической базы данных (МГБД-1).

Этот модуль МГБД служит фундаментом для эпидемиологических и генетических исследований на уровне популяции. Так как медико-генетические базы данных появились сравнительно недавно, то разработка методов для работы с ними является актуальной. МГБД необходимы для проведения полногеномных ассоциативных исследований с использованием как генетических, так и клинических матричных данных. Поэтому особый

интерес представляют методы вероятностной кластеризации для поиска несмещенной контрольной группы для таких исследований.

Основным результатом GWAS является количественная оценка того, насколько конкретный ДНК-вариант связан с определённым фенотипом или признаком, хранящимся в МГБД-1. Для каждого ДНК-варианта строится линейная или логистическая модель и рассчитывают:

- размер эффекта - насколько сильно изменяется признак при наличии альтернативного аллеля по сравнению с референсным;

- статистическую значимость (р-значение), которая определяет, что найденная ассоциация не является случайной;

- z-статистику, которая отражает, насколько сильно и достоверно отличается оцененный эффект от нуля с учётом его дисперсии.

Результаты нескольких GWAS могут быть объединены с целью формирования второго модуля МГБД (МГБД-2). Такой модуль может содержать сотни тысяч связей между ДНК-вариантами и фенотипами.

МГБД-2 хранит ассоциации между ДНК-вариантом и признаком (болезнью), однако, чтобы определить взаимосвязь между признаком и геном необходима дополнительная информации об экспрессиях генов и белок-белковых взаимодействиях.

МГБД-2 всё ещё является относительно новым форматом, для которого активно разрабатываются различные методы интерпретации (аннотации) данных. При этом автором предложен вероятностный метод графовой кластеризации, который, в том числе, может использоваться для анализа данных в МГБД-2. С помощью этого метода возможно выделить активный модуль в биологическом графе по результатам полногеномных ассоциативных исследований, хранящихся в МГБД-2. Выделение активного модуля является разновидностью задачи вероятностной кластеризации.

В разделе 1.6 представлен один из возможных подходов к интерпретации (аннотации) результатов сравнительного анализа, при котором проводится переход от ассоциаций между наблюдаемыми характеристиками и низкоуровневыми переменными к интерпретации (аннотациии) через функционально значимые компоненты системы. Такой переход позволяет учитывать не только прямые зависимости, но и каскадные или опосредованные связи, влияющие на результат.

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

дополнительно учитывают априорные данные о том, как каждый влияет на экспрессию гена (eQTL-информация [36]). Таким образом, для гена g формируется генный z-показатель, который агрегирует z-статистики для всех SNP, участвующих в экспрессии данного гена, с учётом их весов (полученных из eQTL-модели) и нормирования по дисперсиям [37]. Результирующий показатель позволяет судить о совокупном вкладе SNP в фенотип с учетом генной экспрессии. Таким образом, закладывается основа объединения нескольких типов данных для решения задачи по эффективному поиску причинных генов в генетических исследованиях. При этом важным компонентом является учет взаимодействий в биологических системах. Такие взаимодействия могут быть представлены в виде графа.

Таким образом, в главе выполнено описание вероятностных методов для решения задач кластеризации матричных и графовых данных, а также методы, необходимые для разработки базы данных с разнородной информацией.

Во второй главе описан метод вероятностной кластеризации матричных данных с априорным распределением характеристик для получения несмещенного набора данных. Метод применен автором к МГБД-Р-1.

В разделе 2.1 сформулирована задача анализа матричных данных в МГБД.

Предположим, имеется два набора данных - тестовый и контрольный. В тестовом наборе собрана информация о пациентах, у которых наблюдается определённое заболевание или особый фенотип, а в контрольном наборе, который хранится в МГБД, содержатся данные большой популяции без этой болезни. Возникает задача найти в обширном контрольном наборе кластер, который будет максимально похож на тестовый набор по распределению ключевых факторов, чтобы избежать систематические искажения при сравнительном анализе.

Пусть А1 и А2 - матрицы признаков; dim(Al) = т хт, dim(A2) = П2хт, где т - число различных характеристик (например, возраст, пол, различные лабораторные показатели, социальные факторы), т - число участников в тестовой группе, П2 - в контрольной. Требуется подобрать §2 £ S2 -подмножество контрольного набора, чтобы А1 (подмножество строк/столбцов в тестовом наборе) и А2 - (аналогичное подмножество данных МГБД) не отличались значимо по статистическим признакам и корреляциям. Такой подход обеспечивает корректное сопоставление тестового и контрольного набора данных без существенных систематических смещений.

В разделе 2.2 описан предлагаемый метод вероятностной кластеризации матричных данных с априорным распределением характеристик для получения несмещенного набора данных на основе сингулярного разложения (singular value decomposition, SVD) [31] матриц.

Этапы постановки математической задачи

1. Снижение размерности тестового и контрольного наборов данных

На первом этапе задача сводится к приведению двух наборов данных -тестового и контрольного - к единому признаковому пространству. Для этого проводится стандартизация признаков, устраняющая возможные смещения между наборами. Далее применяется такой метод снижения размерности, как сингулярное разложение матриц, позволяющий выделить наиболее информативные направления изменения данных. Оба набора проецируются в это общее подпространство меньшей размерности, где сохраняются главные закономерности, но исключаются второстепенные и шумовые компоненты.

2. Вероятностная кластеризация в пространстве признаков главных компонент

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

В рамках данного метода на основе тестового набора строится априорное распределение признаков. Затем для каждого объекта из контрольного набора вычисляется вероятность его соответствия этой априорной модели. Если объект с высокой вероятностью может быть отнесён к распределению тестового набора, он включается в формируемый кластер. Таким образом, кластер формируется как совокупность объектов, которые статистически близки к тестовому набору по заранее заданным критериям доверия.

Этапы решения математической задачи

1. Снижение размерности тестового и контрольного наборов данных

1.1. В тестовом наборе А1 6 М$1Хт вычисляются средние значения по каждому столбцу. Они затем вычитаются из обеих матриц А1 и А2 6 ^п2Хт, что обеспечивает выравнивание выборок по среднему значению признаков и устраняет глобальные линейные смещения.

1.2. К матрице А1 применяется сингулярное разложение: А1 = и£ Vх, где и 6 М$1ХП1 - матрица левых сингулярных векторов, £ -диагональная матрица сингулярных чисел, V 6 МтХт - матрица правых сингулярных векторов. Из матрицы и выбираются первые к столбцов, соответствующих наибольшим сингулярным числам. Полученная матрица и 6 М$1Хк определяет проекцию в подпространство наибольшей дисперсии.

1.3. Проекции обоих наборов данных вычисляются следующим образом: А1 = ИА1 и А2 = ИА2, соответственно. Каждая строка матрицы А1 и А2 теперь представляет наблюдение в новом пространстве признаков размерности к общем для обеих выборок. Это обеспечивает сопоставимость обоих наборов данных.

2. Вероятностная кластеризация в пространстве признаков главных компонент

2.1. Оценка априорного распределения признаков в тестовой выборке. После проекции на пространство главных компонент А1, компоненты тестовой выборки рассматриваются как реализация случайной величины, подчиняющейся многомерному нормальному распределению. Оцениваются вектор математических ожиданий ^ 6 Мк и ковариационная матрица £ 6 МкХк, задающие априорное распределение признаков: А1 ~

2.2. Формулировка вероятностного критерия принадлежности. Для каждого объекта х 6 А2 вычисляется расстояние Махаланобиса [38] как мера правдоподобия его принадлежности к распределению тестовой выборки: Dм(x) = (х - ^)т£-1(х -Таким образом, каждый объект из контрольной выборки получает

количественную оценку степени принадлежности к априорному распределению.

2.3. Выделение вероятностного кластера. В отличие от детерминированного отнесения объектов к кластерам, используется вероятностный критерий: объект включается в кластер §2, если значение Dм(x) не превышает порог, соответствующий 95 %-квантилю распределения с2 с к степенями свободы. Это соответствует включению в доверительную область уровня 0.95 для многомерного нормального распределения. Таким образом, кластер §2 £ А2 формируется как §2 = {х е А2 | Dм(Uтx) < Хк,1-а}.

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

В разделе 2.3 описана валидация предлагаемого вероятностного метода в двух экспериментах. Была промоделирована тестовая матрица А1 е КтХт, где т = 200 и т = 30, и контрольная матрица А2 е ЖП2Хт, т = 20 000. В первом эксперименте обе выборки генерировались из одного многомерного нормального распределения Ж(0, 2) с ковариацией Тёплица [39]. Во втором эксперименте контрольная группа смещалась на экспериментально подобранную величину 2,7 в пространстве главных компонент. Оба набора проецировались в подпространство размерности девять. Для тестовой группы оценивались вектор средних ц и ковариационная матрица 2, а отбор искомой контрольной подгруппы из матрицы А2 выполнялся по вероятностному критерию Махаланобиса с порогом Л"(х) = (х - т)Т2-1(х - т) < 1 а, где а =

0,05. В качестве сравнения применялся метод ближайшего соседа [38] по схеме «один к одному» по расстоянию Махаланобиса в исходном пространстве признаков. Было показано, что в первом эксперименте предлагаемый метод корректно выделял сопоставимую контрольную группу (N=19 744, рисунок Р.1 (часть А)). Во втором эксперименте предлагаемый метод не отбирал ни одного наблюдения из контрольной группы (ноль ложноположительных наблюдений), тогда как метод ближайшего соседа формировал 200 ошибочных пар ((Рисунок Р.1 (часть Б)). Таким образом, предлагаемый вероятностный метод является более надежным и лучше предотвращает ложноположительные соответствия между тестовой и контрольной группой.

• Тестовая выборка

• Отобранная контрольная выборка Остальная контрольная выборка

Рисунок Р.1. Метод вероятностной кластеризации матричных данных для задачи подбора контрольной группы. А) Первый эксперимент; Б) Второй эксперимент.

В разделе 2.4 описано применение предложенного метода вероятностной кластеризации к реальному когортному исследованию, направленному на выявление фенотипических особенностей потомков жителей блокадного Ленинграда (тестового набора данных). Метод был применен к разработанному МГБД-Р-1, чтобы определить несмещенный контрольный набор данных.

В тестовом наборе были взяты данные о питании и ряде клинических показателей (в том числе работоспособности почек). Изначально был использован простой метод выбора контрольной группы, основанный лишь на согласовании возраста и пола по данным из МГБД-Р-1. Однако этот метод показал, что тестовая группа заметно отличалась от контрольной группы не только по работе почек, но и по ряду пищевых привычек, что делало результаты сравнительного анализа смещёнными (Рисунок Р.2 (часть А)).

Чтобы исключить такое смещение, был применён предложенный метод вероятностной кластеризации (раздел 2.2) для матрицы признаков, содержащей пищевые привычки. Была сформирована матрица А1, содержащая особенности питания для потомков жителей блокадного Ленинграда, и выполнено её сингулярное разложение А1 = и £ Vх. Затем выбраны шесть первых главных направлений (и), которые объясняют основную дисперсию в данных. Матрица И была использована для расчета матрицы А2. Был сформирован кластер §2 из МГБД-Р-1 так, чтобы он максимально совпадал с тестовой группой на основании расстояний Махаланобиса в пространстве главных компонент, тем самым устраняя систематическое различие между группами (Рисунок Р.2 (часть Б)).

В результате такой корректировки с помощью логистической регрессии было показано, что контрольный набор данных, сопоставимый с тестовым набором по факторам питания, демонстрирует статистически значимые

отличия в работе почек - причём эти отличия уже не могут быть объяснены смещением по пищевому рациону. Таким образом, предложенный метод подтвердил независимое отклонение функционального состояния почек у потомков жителей блокадного Ленинграда и позволил исключить гипотезу, что причина выявленного фенотипа обусловлена особенностями питания (Рисунок Р.2 (часть В)).

• Тестовая группа

Рисунок P.2. Применение метода вероятностной кластеризации матричных данных для задачи подбора контрольной группы. A) Первые две главные компоненты, выделенные из матрицы пищевых факторов, показывают расположение тестовой группы (оранжевая маркировка) относительно контрольной группы в пространстве главных компонент; Б) Сравнение плотностей распределения второй главной компоненты для четырёх выборок: (1) тестовой группы, (2) контрольной группы, подобранной только по полу и возрасту, (3) контрольной группы, сформированной предложенным методом SVD, и (4) всей популяционной выборки из МГБД; В) Сравнительный анализ двух фенотипов, связанных с работой почек (скорость клубочковой фильтрации и уровень креатинина), в той же тестовой группе и различных контрольных группах

Список литературы диссертационного исследования кандидат наук Усольцев Дмитрий Андреевич, 2025 год

Литература

1. Huber W„ Carey V.J., Long L., Falcon S.. Gentleman R. Graphs in molecular biology // BMC Bioinformatics. 2007. V. 8. Suppl. 6. P. S8. https://doi.org/l0.1186/l471-2105-8-S6-S8

2. Szczepanski A.P.. Wang L. Emerging multifaeeted roles of BAP1 complexes in biological processes // Cell Death Discovery. 2021. V. 7. N 1. P. 20. https://doi.oig/10.1038/s41420-021-00406-2

3. Carbone M., Yang H.. Pass H.I.. Krausz T.. Testa J.R.. Gaudino G. ВДР1 and cancer//Nature Reviews Cancer. 2013. V. 13. N3. P. 153 159. https://doi.Org/10.1038/nn:3459

4. I.in J.S.. Lai E.M. Protein-protein interactions: Co-Imiminoprecipitation и Methods in Molecular Biology. 2017 V 1615. P. 211-219. https://doi.org/10.1007/978-l-4939-7033-9_i7

5. Tamara S„ den Boer M.A., Heck AJ.R. High-resolution native mass spectrometry//Chemical Reviews. 2022.V. 122. N8. P. 7269- 7326. https://doi.org/10.1021/acs.chemrev. Ie00212

6. Okpara M.O., Hermann C., van (let Walt P.J.. Oarnelt S,, Blackburn J.M., Leaner V.D. A mass spectrometry-based approach for the identification of Kpn(tl binding partners in cancer cells // Scientific Reports. 2022. V. 12.N 1 P. 20171 https://doi.org/10.1038/ s41598-022-24194-6

7. Li T.. Wernersson R-, Hansen R.B., Horn H.. Mercer J.. Slodkowiez G.. Workman C.T., Rigina O.. Rapaeki K., Starfeldt H.H., Brunak S., Jensen T.S.. Lage K. A scored human protein-protein interaction network to catalyze genomic interpretation // Nature Methods. 2017. V. 14. N 1* P. 61-64. https://doi.org/10.1038/ nmeth.4083

8. Zhu Q.M.. Hsu V.H.. Lassen F.H., Mac-Donald B.T.. Stead S„ Malolepsza £., Kim A.. Li T.. Mizoguehi T., Schenone M., Guzman G., Tanenbaum В., Fornelos N.. Сагт S.A., Gupta R.M.. Hllinor P.T., Lage K. Protein interaction networks in die vasculature prioritize genes and pathways underlying coronary arteiy disease // Communications Biology. 2024. V. 7. N 1. P. 87. htlps://doi. org/10.1038/s42003-023-05705-l

9. Nehmc R.. Pietilainen O., Artomov M.. Tegtmeyer M.. Valakh V., Lehumen L.. Bell C., Singh T., Trehan A., Sherwood J., Manning D.. Peirenl E„ Malik R.. Guss E.J., Hawes D.. Beccard A.. Bara A.M.. Hazelbaker D.Z., Zuccaro E., Genovesc G., Loboda A.A., Neumann A,, Lilliehook C., Kuismin O., Hamalainen E.. Kurki M.. Hultman C.M.. Kahler A.K.. Paulo J.A.. Ganna A , Madison J . Cohen В., MePhie D., Adolfsson R., Perlis R., Dolmetsch R., Failli S., McCarrol! S.. Hyman S.. Neale В., BaneU L.E.. Harper W.. Palotie A.. Daly M., Eggan K. The 22ql 1.2 region regulates presynaptic gene-products linked to schizophrenia // Nature Communications. 2022. V. 13, N 1. P. 3690. https://doi.org/10.1038/s41467-022-31436-8

10 Nguyen H., Shrcstha S.. Tran D„ Shafi A.. Draghici S„ Nguyen T. A Comprehensive survey of tools and software for active subnetwork

References

1. Huber W., Carey V.J., Long L., Falcon S.. Gentleman R. Graphs in molecular biology. BMC Bioinformalics. 2007. vol- 8. suppl. 6. pp. S8. https://doi.Org/IO.l 186/1471 -2105-8-S6-SS

2. Szczepanski A.P., Wang L. Emerging multifaeeted roles of BAP1 complexes in biological processes. Cell Death Discovery, 2021, vol. 7, no. 1, pp. 20. https://doi.org/10.l038/s4I420-021-00406-2

3. Carbone M„ Yang H., Pass H.L, Krausz I . Testa J.R., Gaudino G. BAP1 and cancer. Nature Reviews Cancer, 2013, vol. 13, no. 3, pp. 153 159, https://doi.org/10.1038/nre3459

4. Lin J.S.. Lai E.M. Protein-protein interactions: Co-Immunoprecipilation. Methods in Molecular BioioQ', 2017, vol. 1615, pp. 211 -219. hllps://doi.org/10.1007/978-1 -4939-7033-9 J7

5. Tamara S., den Boer M.A., Heck AJ.R. High-resolution native mass spectrometry. Chemical Reviews, 2022, vol. 122, no. 8. pp. 72697326. https://doi.oig/10.1021/acs.chemrev.lc00212

6. Okpara M.O.. Hermann C„ van der Watt P.J., Ciarnetl S.. Blackburn J.M.. Leaner V.D. A mass spectrometry-based approach for the identification of Kpnfl 1 binding partners in cancer cells. Scientific Reports, 2022, vol. 12, no. 1, pp. 20171. htlps://doi. org/10.1038/s41598-022-24194-6

7. Li T-, Wernersson R-. Hansen R.B.. Horn H-, Mercer J.. Slodkowiez G.. Workman C.T., Rigina O, Rapaeki K.. Suerfeldt H.H., Brunak S , Jeusen T.S.. Lage K A scored human protein-protein interaction network to catalyze genomic interpretation. Nature Methods, 2017. vol. 14, no. 1, pp. 61-64. https://doi.org/10.1038/ nmeth.4083

8. Zhu Q.M.. Hsu Y.H., Lassen F.H.. MocDonald B.T., Stead S.. Malolepsza E., Kim A., Li T-, Mizoguchi T-, Schenone M-, Guzman G„ Tanenbaum B.. Fomelos N., Carr S.A., Gupta R.M., Ellinor P.T., Lage K. Protein interaction networks in the vasculature prioritize genes and pathways underlying coronary artery disease. Communications Biology, 2024, vol. 7, no. 1, pp. 87. htlps://doi. org/10.1038/S42003-023-05705-1

9. Nchme R.. Pietilainen O.. Artomov M-, Tegtmeyer M.. Valakh V.. Lehtonen L., Bell C., Singh T., Trehan A., Sherwood J., Manning D.. Peirent E.. Malik R„ Guss E.J., Hawes D.. Beccard A.. Bara A.M.. Hazelbaker D.Z., Zuccaro E., Genovese G.. Loboda A.A.. Neumann A., Lilliehook C., Kuismin O., Hamalainen E., Kurki M., Hultman C.M.. Kahler A.K.. Paulo J.A.. Ganna A.. Madison J.. Cohen B„ MePhie D., Adolfsson R., Perlis R„ Dolmetsch R., Farhi S„ McCarroll S.. Hyman S.. Neale B.. Barrett L.E., Harper W„ Palotie A.. Daly M., Eggan K. The 22ql 1.2 region regulates presynaptic gene-producls linked to schizophrenia. Nature Communications, 2022, vol. 13, no. 1. pp. 3690. https://doi.org/I0.1038/s41467-022-3I436-8

10. Nguyen H„ Shreslha S., Tran D., Shafi A„ Draghici S., Nguyen T. A Comprehensive survey of toots and software for active subnetwork

identification // Frontiers in Genetics. 2019. V. 10. P. 155. https://doi. org/10.3389/fgene.2019.00155

11. Mitra K., Carvunis A.R., Ramesh S.K., Ideker T. Integrative approaches for finding modular structure in biological networks H Nature Reviews Genetics. 2013. V. 14. N 10. P. 719-732. https://doi. org/10.103 8/nrg3 55 2

12. Strauss B.S. Biochemical genetics and molecular biology: The contributions of George Beadle and Edward Tatum // Genetics. 2016. V. 203. N 1. P. 13-20. https://doi.org/10.1534/genetics.116.188995

13. Montecino-Rodriguez E., Casero D., Fice M., Le J., Dorshkind K. Differential expression of PU.l and key T lineage transcription factors distinguishes fetal and adult T cell development // Journal of Immunology. 2018. V. 200. N 6. P. 2046-2056. https://doi. org/10.4049/jimmunol. 1701336

14 Suzuki K., Hatzikotoulas K.. Southam L., Taylor N.J., Yin X., Lorenz K.M. et al. Genetic drivers of heterogeneity in type 2 diabetes pathophysiology // Nature. 2024. V. 627. P. 347-357. https://doi. org/10.1038/s41586-024-07019-6

15. Kim T.K., Park J.H. More about the basic assumptions of t-test: normality and sample size // Korean Journal of Anesthesiology. 2019. V. 72. N4. P. 331 335 https: 'doi.org/10.4097/kja.d.IS.002'»2

16. Barton S.J., Crozier S.R.. Lillycrop K.A, Godfrey K.M.. Inskip H.M. Correction of unexpected distributions of P values from analysis of whole genome arrays by rectifying violation of statistical assumptions// BMC Genomics. 2013. N 14 P. 161. https://doi. org/10.1186/1471-2164-14-161

17 Alexeev N.. lsomurodov J., Sukhov V.. Korotkevich G.. Sergushichev A. Markov chain Monte Carlo for active module identification problem // BMC Bioinformatics. 2020. V. 21. Suppl. 6. P. 261. https://doi.org/10.1186/sl2859-020-03572-9 18. Dittrich M.T., Klau G.W., Rosenwald A., Dandekar T., Müller T. Identifying functional modules in protein-protein interaction networks: an integrated exact approach // Bioinformatics. 2008. V. 24. N 13. P. ¿223-1231. https://doi.org/10.1093/bioinforrnatics/btnl61 19 Zhu Z.. Zhang F., Hu H., Bakshi A.. Robinson M.R., Powell J.E., Montgomery G.W., Goddard M.E., Wray N.R., Visscher P.M., Yang J. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets // Nature Genetics. 2016. V. 48. N 5. P. 481-487. https://doi.org/10.1038/ng.3538

20. Chen f.. Guestrin C. XGBoost: A scalable tree boosting system // Proc. of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. 2016. P. 785 794. https:// doi.org/10.1145/2939*672.2939785

21. Warde-Farley D, Donaldson S.L., Comes O., Zuben K.. Badrawi R.. Chao P., Franz M., Grouios C.. Kazi F., Lopes C.T., Maitland A.. Mostafavi S., Montojo J.. Shao Q., Wright G , Bader G.D., Morris Q. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function // Nucleic Acids Research. 2010. V. 38. Suppl. 2. P. W214-W220. https://doi. org/10.1093/nar/gkq537

identification. Frontiers in Genetics, 2019. vol. 10, pp. 155. https:// doi.org/10.33 89/ fgene.2019.00155

11. Mitra K., Carvunis A.R., Ramesh S.K., Ideker f. Integrative approaches for finding modular structure in biological networks. Nature Reviews Genetics, 2013. vol. 14, no. 10. pp. 719 732. https:// doi.org/10.1038/nrg3552

12. Strauss B.S. Biochemical genetics and molecular biology: The contributions of George Beadle and Fdward Tatum. Genetics, 2016, vol.203, no. 1. pp. 13-20. https://doi.org/10.1534. genetics. 116.188995

13. Montecino-Rodriguez E., Casero D., Fice M., Le J., Dorshkind K. Differential expression of PU. 1 and key T lineage transcription factors distinguishes fetal and adult T cell development. Journal of Immunology, 2018, vol. 200, no. 6, pp. 2046-2056. https://doi. org/10.4049/jimmunol. 1701336

14. Suzuki K., Hatzikotoulas K., Southam L., Taylor H.J , Yin X., Lorenz K.M. et al. Genetic drivers of heterogeneity in type 2 diabetes pathophysiology. Nature, 2024, vol. 627, pp. 347 357. https://doi. org/10.1038/s41586-024-07019-6

15. Kim T.K., Park J.H. More about the basic assumptions oft-test: normality and sample size. Korean Journal of Anesthesiology^ 2019, vol. 72, no. 4, pp. 331-335. https://doi.Org/10.4097/kja.d.18.00292

16. Barton S.J., Crozier S.R., Lillycrop K.A., Godfrey K.M.. Inskip H.M. Correction of unexpected distributions of P values from analysis of whole genome arrays by rectifying violation of statistical assumptions. BMC Genomics, 2G13, no: 14. pp. 161 fops doi.org KM186/H71-2164-14-161

17. Alexeev N., lsomurodov J., Sukhov V., Korotkevich G., Sergushichev A. Markov chain Monte Carlo for active module identification problem. BMC Bioinformatics, 2020. vol. 21. suppl. 6, pp. 261. https://doi.org/10.1186/s 12859-020-03572-9

18. Dittrich M.T., Klau G.W., Rosenwald A., Dandekar T.. Mfiller T. Identifying functional modules in protein-protein interaction networks: an integrated exact approach. Bioinformatics, 2008, vol. 24, no. 13, pp. i223- i231. https://doi.org/10.1093/bioinfonnatics/bml61

19. Zhu Z., Zhang F., IIu H.t Bakshi A., Robinson MR.. Powell J.E., Montgomery G.W.. Goddard M.E., Wray N.R., Visscher P.M., Yang J. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nature Genetics, 2016, vol. 48. no. 5, pp. 481-487. https://doi.org/10.1038/ng.3538

20. Chen T.. Guestrin C. XGBoost: A scalable tree boosting system. Pwc. of the 22'"1 ACM SIGKDD International Conference on Knowledge Discovery and Data Aiming, 2016. pp. 785-794. https://doi. org/10.1145/2939672.2939785

21. Warde-Farley D., Donaldson S.L., Comes O., Zuberi K., Badrawi R., Chao P., Franz M., Grouios C., Kazi F., Lopes C.T., Maitland A., Mostafavi S.. Montojo J.. Shao Q„ Wright G.. Bader G.D.. Morris Q. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Research. 2010, vol. 38, suppl. 2, pp. W214-W220. https://doi. org/10.1093/nar/gkq537

Авторы

Усольцев Дмитрий Андреевич аспирант. Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация: старший научный сотрудник. Институт геномной медицины. Детская больница Нейшеивайд, Колумбус. 43205, США. sc 57279360300. https://orcid. org/0000-0001 -8072-31 OX, [e-mail] Молотков Иван Игоревич старший научный сотрудник, Институт геномной медицины, Детская больница Нейшеивайд, Колумбус. 43205, США; аспирант. Медицинский колледж Университета штата Огайо, Колумбус. 43210, США. sc 58651494600, https://orcid.org/0009-0008-3566-0160, [e-mail] Артемов Никита Николаевич кандидат химических наук, доцент, главный исследователь. Институт геномной медицины. Детская больница Нейшеивайд, Колумбус, 43205, США; профессор. Медицинский колледж Университета штата Огайо, Колумбус. 43210. США, «с 36542095500, https://orcid.oig/0000-0001-5282-8764, mykyta. [e-mail]

Сергушичев Алексей Александрович — кандидат технических наук, доцент. Университет Вашингтона в Сент-Луисе. Сент-Луис. 63110, США. sc 55772694000. https://orcid.org/0000-0003-1159-7220, [e-mail]

Authors

Dmitrii A. Usoltsev PhD Student, ITMO University, Saint Petersburg, 197101. Russian Federation: Senior Researcher, Institute for Genomic Medicine, Nationwide Chldren's Hospital, Columbus, 43205, USA, sc 57279360300, https://orcid.org/0000-0001-8072-310X, dusoltsev.27@ gmail.com

Ivan I. Molotkov Senior Researcher. Institute for Genomic Medicine, Nationwide Chldren's Hospital, Columbus. 43205. USA: PhD Student, The Ohio State University College of Medicine, Columbus, 43210, USA. sc 58651494600, https://orcid.org/0009-0008-3566-0160. ivan. [e-mail]

Mykyta N. Artomov — PhD (Chemistry), Associate Professor. Chief Researcher, Institute for Genomic Medicine. Nationwide Chldren's Hospital, Columbus, 43205, USA; Professor, The Ohio State University College of Medicine. Columbus. 43210, USA. sc 36542095500. https:// orcid.org/0000-0001 -5282-8764, mykyta.artomov@nationwidechiIdrens. org

Alexey A. Sergushichev PhD, Associate Professor, Washington University School of Medicine in St. Louis, St. Louis, 63110, USA. so 55772694000, https://orcid.org/0000-0003-l 159-7220, asergushiehev@ wustl.edu

НАУЧНО-ТЕХНИЧЕСКИЙ ВЕСТНИК ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ. МЕХАНИКИ И ОПТИКИ

1/iTMO

май-июнь 2025

Том 25 № 3

http://ntv.itmo.ru/

SCIENTIFIC AND TECHNICAL JOURNAL OF INFORMATION TECHNOLOGIES. MECHANICS AND OPTICS

ННШЯРМАЦИОИНЫК1ЕХВОЛОГИИ. MtXAHHKM И ОПТИКИ

Mav-June 2025 ISSN 2226-1494 (print)

Vol 25 No 3

http://ntv.Hmo.ru/en/ ISSN 2500-0373 (onhnel

МАТЕМАТИЧЕСКОЕ И КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ

MODELING AND SIMULATION

(1о|: 10.17586/2226-1494-2025-25-3-487-497 УДК 004.961

Метод определения активного модуля в биологических графах с многокомпонентными весами вершин Дмитрии Андреевич Усольцев1 , Иван Игоревич Молотков2, Никита Николаевич Артемов-1, Алексеи Александрович Сергушнчсв4, Анатолий Абрамович Шалыто5

1 - -1 Институт геномной медицины. Детская больница Нейшенвайд, Колумбус, 43205. США 1 * Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация --1 Медицинский колледж Уттверситета штата Огайо, Колумбус, 43210, США

4 Университет Вашингтона в Сент-Луисе, Сент-Луис, 63110, США 1 [e-mail]1, https://orcid.org/0000-0001-8072f310X

: [e-mail], https://orcid.org/0009-0008-3566-0160

5 [e-mail]. https://orcid.org/0000-0001-5282-8764

4 asergushichevfiiwustl.edu, https://orcid.org/0000-0003-1159-7220

5 [e-mail], https://orcid.org/0000-0002-2723-2077

Введение. Активный модуль в биологических графах представляет собой свя занный подграф, вершины которого объединены общей биологической функцией. Для определения активного модуля необходимо сначала построить взвешенный биологический граф. Вес каждой вершины вычисляется на основе биологических экспериментов, исследующих искомую биологическую функцию. Однако результаты одного эксперимента могут не полност ью описывать искомый активный модуль, а лишь его часть, внося, например, неопределенность в веса вершин. В работе показано, что использование метода Фишера для объединения данных нескольких экспериментов, а затем применение метода Монте-Карло по схеме марковских цепей (МКМЦ) и машинного обучения к результатам метода Фишера, позволяет более эффективно выделять акт ивные модули в биологических графах. Метод. В работе используются граф белок-белковых взаимодействий — InWcblM, граф по реконструкции мозга человека из проекта BigBrain и генный граф для вида живых организмов Caenorhabditis elegans. Для объединения результатов нескольких экспериментов в одном графе в единый результат применяется метод Фишера. После этого поиск активных модулей выполняется с использованием метода МКМЦ и машинного обучения. Для валидации предлагаемого метода на реальных данных применяются результаты полногеномного ассоциативного исследования по шизофрении и курению, а также матрица экспрессии генов пациентов с кожной меланомой из проекта The Cancer Genome Allas. Основные результаты. Применение метода Фишера позволяет учитывать результаты песколъких биологических экспериментов одповремешю. Последующее использование метода МКМЦ и машинного обучения повышает точность определения активных модулей по сравнению с ранжированием вершин графа только на основе метода Фишера. Обсуждение. Учет результатов нескольких биологических экспериментов при определении активных модулей играет ключевую роль в повышении точности нахождения веришн активного модуля. Это способствует лучшему пониманию биологических механизмов заболеваний, что может иметь важное значение для разработки новых методов диагностики и терапии. Ключевые слова

!рафы, метод Монте-Карло по схеме марковских цепей, метод Фишера, биологические графы, активный модуль Ссылка для цитировании: Усольцев Д.А., Молотков И.И., Артемов H.H., Сергупитчев A.A., Шалыто A.A. Метод определения активного модуля в биологических графах с многокомпонентными весами вершин // Научно-технический вестник информационных технологий, механики и оптики. 2025. Т. 25, № 3. С. 487^197. doi: 10.17586/2226-1494-2025-25-3-487-497

© Усольцев Д.А., Молотков И.И.. Артемов H.H., Ссргушичев A.A., Шалыто A.A., 2025

Аннотация

Method for identifying the active module in biological graphs with multi-component vertex weights Dmitrii A. Usoltsev1 , Ivan I. Molotkov2, Mvkyta N. Artomov3, Alexey A. Sergushicher4,

Anatoly A. Shalyto5

'■-J Institute lor Genomic Medicine, Nationwide Children's Hospital, Columbus, 43205, USA 15 ITMO University, Saint Petersburg, 197101, Russian Federation --3 The Ohio Slate University College of Medicine, Columbus, 43210, USA 4 Washington University School of Medicine in St- Louis, St. Louis, 63110, USA 1 [e-mail]®, https.7/orcid.org/0000-0001-8072-31 OX : [e-mail], https://orcid.org/0009-0008-3566-0160 [e-mail], https://orcid.org/0000-0001-5282-8764

4 [e-mail], hltps://oreid.org/0000-0003-l 159-7220

5 [e-mail], https://orcid.org/0000-0002-2723-2077

Abstract

An active module in biological graphs is a connected subgraph whose vertices share a common biological liinction. To identity an active module, one must first construct a weighted biological graph. The weight of each vertex is calculated based on biological experiments investigating the target biological function. However, the results of a single experiment may not fully describe the desired active module, covering only part of it and potentially introducing uncertainty into the vertex weights. This work demonstrates that employing Fisher's method to integrate data from multiple experiments followed by applying a Markov chain Monte Carlo (MCMC) and machine learning-based approach to the results of Fisher's method, enables more effective identification of active modules in biological graphs. lite study utilizes the InWebIM protein-protein interaction graph, a human brain reconstruction graph from the BigBrain project, and a gene graph for the organism Caenorhahditis eiegans. To combine the results of several experiments into a single outcome within one graph. Fisher's method is applied. Aflerwards, the search for active modules is conducted using an MCMC and machine learning-based method. To validate the proposed method on real data, results from (ienome-Wide Association Studies on schizophrenia and smoking are used, along with the gene expression matrix of patients with skin melanoma from the TCGA project. Applying Fisher's method makes it possible to consider the results of multiple biological experiments simultaneously. Subsequent use of Ihe MCMC and machine learning-based method improves the accuracy of identifying active modules compared to ranking graph vertices solely by Fisher's method. Considering the results of multiple biological experiments when determining active modules plays a crucial role in increasing the accuracy of identifying the vertices of the active module, 'fhis, in turn, promotes a deeper understanding of the biological mechanisms of diseases, which can be of great significance for the development of new diagnoslic and therapeutic methods. Keywords

graphs, MCMC, Fisher's method, biological gTaphs, active module

For citation: Usoltsev D. A., Molotkov 1.1., Artomov M.N., Sergushichev A. A., Shalyto A. A. Method for identifying the active module inbiological graphs with multi-component vertex weights. Scientific and Technical Journal of information Technologies, Mechanics and Optics, 2025, vol. 25, no. 3, pp. 487^197 (in Russian), doi: 10.17586/2226-1494-2025-25-3-487-497

Введение

Графы часто используются при анализе сложных биологических взаимодействий, таких как белок-белковые взаимодействия, ко-экспрессия генов или регуляция метаболических путей (1-3]. Внутри таких графов существуют связанные подграфы, вершины которых выполняют общую биологическую функцию. Эти подграфы называют активными модулями (4). Модули, например, в графе белок-белковых взаимодействий позволяют понять причины возникновения заболеваний и определить возможные терапевтические мишени [5, 6]. Таким образом, определение активных модулей является одной из важных задач прн анализе биологических графов.

В работах [7. 8] представлены подходы для определения активных модулей на основе метода Монте-Карло по схсмс марковских цепей (МКМЦ) и метода на основе машинного обучения — градиентного бустинга. Однако точность определения активного модуля как с помощью только метода МКМЦ [7]. так и с помощью подхода, предложенного в [8], зависит от точности

определения весов вершин графа. " " \^FDR

Для определения весов в биологическом графе проводятся эксперименты по дифференциальной экспрессии генов пли мутационной нагрузке между больными и здоровыми группами людей. Результаты экспериментов анализируются с помощью статистического теста (например, 1-тест, линейная регрессия) [9]. В итоге формируется список генов с соответствующими /7-значениями, отражающими значимость ассоциации с исследуемым признаком. Распределение /^-значений хорошо аппроксимируется бета-равномерным распределением BUM (а, X) с плотностью, вычисляемой на основе соотношения:

/(.y) = A. + (1 -X.)av«-1,0<x,A., a< 1, (1)

где \ — вес равномерной компоненты; a — параметр формы бета-компоненты 110].

В работах |4, 9] предложено, исходя из р-значс-ний. вычислять априорный вес вершины для заданного уровня ложноположительных результатов (False Discovery Rate, FDR) на основе соотношения:

( р V"1

и (п) = - . О < р. а < 1

где р — /»-значение; /;FDR — вероятность того, что /»-значение принадлежит равномерному распределению. Эта вероятность вычисляется по формуле |4. 9]:

/>FDR =

(Х + (1 - A.)g - FDR * X FDR(1 -X)

0<А, а< 1.

ш.

m

Традиционно веса вершин рассчитываются на основе результатов одного статистического теста по формулам (2) и (3). Однако в исследованиях часто доступны результаты множества статистических тестов, отражающих разные аспекты данных: дифференциальная экспрессия, эпигенетические модификации, ассоциации с фенотипами. Результаты каждого теста распределены на основе формулы (1). Учет нескольких статистических тестов для определения априорного веса вершин графа может позволить повысить точность анализа. Такой вес является многокомпонентным (рис. 1).

Для совмещения результатов нескольких статистических тестов используем метод Фишера, который включает вычисление статистики хи-квадрат:

X2 = -2 £ 1п(р<), 0 </>,< 1, i-i

(4)

где р/ — /^-значение /-го эксперимента; п — число экспериментов (статистических тестов) (рис. 1).

Для вычисления вероятности получения статистики X2 из распределения хи-квадрат (х2) со степенью свободы равной удвоенному числу экспериментов [ 11] используется соотношение:

(5)

где рс„ть — результирующие /»-значение: Р — вероятность; df — степень свободы.

Активные модули часто достоверно не известны. Кроме того, один эксперимент позволяет определить лишь часть активного модуля. В этих условиях моделирование активных модулей с заданными параметрами в реальных биологических графах является оптимальным вариантом.

В настоящей работе выполнено моделирование активных модулей с заданными параметрами на подграфах с числом вершин N = 1000 в белок-белковом графе InWcbIM [12]. Предложен метод определения активного модуля для случая, когда известны результаты п экспериментов, каждый из которых покрывает лишь часть активного модуля. Предлагаемый метод основан на последовательном использовании метода Фишера, метода МКМЦ и машинного обучения (используется градиентный бустинг1). Эффективность предлагаемого метода продемонстрирована сравнением его с простым ранжированием результатов, полученных только методом Фишера для разного числа экспериментов. Показана применимость предлагаемого метода в графе из проекта по реконструкции мозга BigBrain [13, 14] и генном графе (для вида живых организмов Caenorhubditis eleguns (HS-LC) [13, 15]. Оценено качество метода Фишера и предлагаемого метода на реальных данных полногеномного ассоциативного исследования (Genome-Wide Association Study, GWAS) no шизофрении, курению и экспрессии генов пациентов с кожной меланомой.

1 Chen Т., 11с Т., Benesty М., Khotilovich V., Tang Y., Cho 11., Chen К, Mitchell R„ Cano I., Zhou Т., Li M„ Xie J., Lin M., Geng Y., Li Y. xgboost: Extreme Gradient Boosting. R package version 1.5.0.2. [Электронный ресурс]. Режим доступа: https:// CRAN.R-project.org/package=xgboost, свободный. Яз. англ. (дата обращения: 17.08.2024).

w,

I Эксперимент р-значение

Pi

О '

fjï-2 VI /

ч/

Ч'.иЧ/ _

I w.v

Экспертшепт /»-значение

1 Р\

2 Рг

Метод Фишера РсотЬ

Рис. 1. Сравнение однокомпонентных (а) и многокомпонентных (Ь) весов вершин, и/, — вес однокомпонентной /-ой вершины; Wj - вес многокомпонентной /-ой вершины; ij G {\,N)

Fig. J. Comparison of single-component (a) and multi-componenl (b) node weights. wf — weight of the /-th single-component node; Wj — weight of the /-th multi-component node; i,j 6 {1, N]

Предлагаемый метод

Эт апы постановки математической задачи

1. Определение активного модуля. В отсутствии достоверно известных активных модулей для большинства биологических состояний проводится симуляция активных модулей с О" вершинами в белок-белковом графе InWebIM посредством поиска в ширину из одной случайной вершины. Для повышения скорости вычислений этот граф сэмплируется (разбивается) на связанные подграфы — упрощенные белок-белковые графы. Число вершин в каждом подграфе N = 1000.

2. Определение весов активного модули. Для вершин активного модуля генерируются ¿»-значения из бета-распрсдслсния, для остальных вершин — из равномерного распределения. Генерация повторяется ч раз. Для того чтобы промоделировать направленность биологических экспериментов на исследование только части активного модуля, на каждом этапе генерации /^-значений 90 % вершин маскируются случайными значениями из равномерного распределения X - f/[0, 1 ]. Таким образом, каждой вершине приписывается п независимых ^-значений, что соответствует л экспериментам. С помощью метода Фишера для каждой вершины вычисляется одно /»-значение. Затем вычисляется априорное значение принадлежности вершины к активному модулю.

Поиск активного модули и опенка качест ва поиска

Поиск активного модуля выполняется с помощью метода, предложенного в 18], с учетом 20 % известных белков активного модуля. Качество определения активного модуля оценивается на основе метрик «площадь под кривой ошибок» (Receiver Operating Characteristic Area Under the Curve, ROC AUC: эта величина варьируется в диапазоне от нуля до единицы) и «чувствительность» (Recall®) 100 — доля правильно определенных белков активного модуля в первых 100 вершинах, ранжированных по уменьшению предсказанной с помощью модели вероятности вхождения в активный модуль; эта величина варьируется в диапазоне от нуля до единицы).

Этапы решения задачи 1. Симуляции 100 различных графов с известными активными модулями

1.1. Для получения графов с топологией, встречающейся в реальных данных, выбираются случайные подграфы из графа InWebIM с числом вершин равным 1000.

1.2. В каждом подграфе InWebIM из п. 1.1 равновероятно среди всех возможных связных подграфов заданного размера выбирается активный модуль с использованием метода МКМЦ из mcmcRanking (vO.l.O)1 при условии, что все вершины имеют одинаковый вес. Модули выбираются после 1000 итераций МКМЦ. чтобы гарантировать их

1 GitHub — etlab/mcmcRanking: Tool To Solve The Active Module Problem [Электронный ресурс]. Режим доступа: https://github.com/ctlab/mcmcRanking, свободный. Яз. аига. (дата обращения: 15.12.2024).

независимость от подграфов, использованных для инициализации МКМЦ.

2. Определение весов активного модуля

2.1. Для каждой вершины симулируются/)-значения эксперимента, проверяющего принадлежность вершины активному модулю. Для вершин вне активного модуля/»-значения выбираются из равномерного распределения, для вершин из активного модуля — из бета-распределения. Для 90 % вершин активного модуля /»-значения маскируются случайными значениями из равномерного распределения на интервале (0, 1].

2.2. П. 2.1 повторяется п раз. таким образом, чтобы каждой вершине соответствовало п /»-значений.

2.3. Для каждого /»-значения применяется метод Фншера: считается X2 по формуле (4) и результирующее /»-значение по формуле (5).

2.4. Полученное /»-значение аппроксимацией бета-равномерным распределением преобразуется в априорный вес по формулам (2) и (3).

3. Поиск активного модуля

3.1. Вероятность того, что подграф является активным модулем, рассчитывается как произведение вероятностей того, что каждая вершина принадлежит этому модулю. Для получившегося вероятностного пространства на множестве связных подграфов формируется выборка размера 100 с помощью метода МКМЦ из тстсИапктц (\0.1.0), используя 10 000 симуляций МКМЦ. После этого для каждой вершины определяется эмпирическая вероятность вхождения в активный модуль — доля подграфов из выборки, которые включают в себя эту вершину.

3.2. Применяется метод, изложенный в [8]. 20% вершин активного модуля выбираются с помощью генератора случайных чисел и используются как известные вершины. Вычисляются расстояния от каждой вершины до трех ближайших белков активного модуля. Используя полученные расстояния и вероятность, определенную в п. 3.1. обучается модель градиентного бустинга с параметрами, описанными в [8]. Так как общее число активных модулей равняется 100, то 50 используются для тренировки модели, а оставшиеся 50 — для тестирования. Известные белки в активных модулях исключаются. При тренировке модели белкам, включенным в активный модуль, присваивается единица, а остальным белкам — ноль.

Экспериментальная оценка предлагаемого метода

Симуляции были проведены для G = 100 и п G {1,3, 10, 50]. Эффективность метода Фишера и метода, изложенного в [8]. оценивалась с помощью метрик ROC AUC и Recall@100. Метрики определялись для каждого подграфа InWeblM, после этого вычислялось среднее значение каждой метрики из выборки в 50 подграфов. При л = 1 результаты, получаемые предлагаемым методом и методом, изложенным в |8] совпали, при условии маскировки /»-значений для 90 % вершин активного модуля случайными значениями из равномерного распределения на интервале (0, 1 ].

Таблица. Сравнение предлагаемого метода и метода Фишера в зависимости от параметров а и н Table. Comparison of the proposed method and Fisher's method depending vs. the parameters а и n

Параметры Предлагаемый метод Метод Фишера

и и ROC AUC (±std) RecallrälOO (±std) ROC AUC (istd) Recall@100(±std)

0,2 1 0,69 ±0,11 0,38 ±0,19 0,54 ± 0,03 0.15±0.«3

0,2 3 0,68 ±0,09 0,36 ±0,17 0,58 ±0,03 0,21 ±0,04

0,2 10 0,74 ±0,07 0,46 ±0,13 0,70 ±0,03 0,36 ±0,04

0,2 50 0,87 ± 0,05 0,70 ±0,11 «,90 ± 0,02 0,68 ± 0,«6

0,8 1 0,61 ±0,12 0,26 ±0,19 0,51 ±0,04 0,12 ± 0,03

0,8 3 0,65 ±0,11 0,30 ±0,21 0,51 ±0,04 0,12 ±0,05

0,8 10 0,62 ±0,10 0,27 ±0,16 0,52 ±0,03 0,11 ±0,03

0,8 50 0,67 ±0,14 0,37 ± 0,24 0,54 ± «,03 0,14 ± 0.04

Примечание. (±std) — стандартное отклонение. Жирным шрифтом выделены значения для и = 1 и п = 50.

В таблице представлены результаты качества определения вершин активного модуля предлагаемым методом и методом Фишера для а = 0,2 и а = 0,8 и число экспериментов п.

При увеличении п в случае а = 0,2, когда бета-равномерное распределение содержит более информативную бета-компоненту, значительно растут как эффективность предлагаемого метода (при увеличении п от 1 до 50 прирост ROC AUC составил 26 %, a Recaíl(ó> 100 — 84 %), так и эффективность метода Фишера (ROC AUC — 67 %. Recall@100 — 353 %). Такой рост значений метрик качества свидетельствует о том. что дополнительная информация (большее число экспериментов) улучшает определение вершин активного модуля. Предлагаемый метод устанавливает активные модули преимущественно эффективнее, чем метод Фишера по обеим метрикам качества, особенно при и от 1 до 10. При п = 50 результаты были схожи для обоих методов.

При увеличении » в случае а = 0,8, когда бета-равномерное распределение ближе к равномерному, незначительно растет как эффективность предлагаемого метода (при увеличении п от 1 до 50 прирост ROC AUC составил 10 %, Rccall(ffil00 — 42 %), так и эффективность метода Фишера (ROC AUC — 6 %, Rccall@100— 17 %). Предлагаемый метод определяет активные модули эффективнее, чем метод Фишера по обеим метрикам качества для всех п.

Для вариантов параметров а 6 {0,2; 0,8} и и 6 {1, 3, 10) предлагаемый метод демонстрирует более высокие значения ROC AUC и Rccall(í¿)100. чем метод Фишера. Однако при а = 0,2 и п = 50 метод Фишера показал более высокое значение метрики ROC AUC. что указывает на то, что в этих условиях информации, полученной только из экспериментов, уже достаточно, а учет взаимодействий в графе становится избыточным. Дальнейшее повышение числа экспериментов п является нецелесообразным для данного набора параметров G и а в графе InWeblM.

Суммарно, предлагаемый метод демонстрирует более высокую или равную точность определения вершин активного модуля по сравнению с ранжированием вершин на основе метода Фишера, что подтвержда-

ется лучшими значениями ROC AUC (рис. 2, a d) и Recall(fl400 (рис. 2, e-h) в большинстве симуляций с ростом значения а.

Предлагаемый метод может быть применен к различным однородным биологическим графам. В качестве примера были выбраны два независимых графа: граф нейронов в головном мозге человека BigBrain (177 584 вершины и 15 669 037 ребер) и генный граф вида живых организмов Caemtrhabditùs elegans (HS-LC) (4227 вершин и 39 484 ребра).

К обоим графам был применен предлагаемый метод. Приращения метрик ROC AUC и RecalIliflOO для предлагаемого метода были аналогичны приращениям этих метрик на графе InWeblM по сравнению с методом Фишера за исключением большого числа экспериментов (и = 50). При таком числе экспериментов и высокой информативности каждого отдельного эксперимента (о < 0,6) метод Фишера несет достаточную информацию о принадлежности вершины к активному модулю. Информация о топологии графа в таком случае становится избыточной, что приводит к отсутствию прироста ROC AUC предлагаемого метода по сравнению с методом Фишера как в случае графа BigBrain (рис. 3). так и снижению качества ROC AUC в среднем на 6 % для генного графа (HS-LC) при а < 0,6 (рис. 4). Recall(«>100 в среднем ухудшается на 8 % для графа BigBrain и на 10 % для графа (HS-LC) при а < 0,6.

Валилания предлагаемого метода с использованием результатов биологических экспериментов

Генетические исследовании шизофрении

Один из основных методов определения связи между ДНК и болезнью — GWAS. Результатом GWAS является множество локусов (отрезков ДНК. включающих от одного до нескольких генов) и ах статистическую значимость.

1. Возникает задача приоритизации генов внутри локуса по отношению к болезни [16. 17].

2. Систематическое сравнение точности метода Фишера и предлагаемого метода проводилось на примере решения задачи из п. I. Было использовано

$ Метод Фишера I Предлагаемый метол

I

О

о ai

0,6

0,4

0.6

'S 0.4

В ог

0,2

1 3 10 49 Число экспериментов

1 3 10 49 Число экспериментов

Рис. 6. Коробчатые диаграммы дам предлагаемою метода и метода Фишера для данных GWAS по пристрастию к курению ROC AUC {<>). Reeall(a)200 (/>): t-тест использован для сравнения метрик между методами (**♦* —/»-значение < 0,0001;

NS — незначимо)

Fig. 6. Appliealion of the proposed method and Fisher's method to GWAS data on smoking addiction. ROC AUC (a). Recall@200 (/»); tlte t-test was used to compare metrics between methods (***» —/»-value < 0.0001; NS — not signilkant)

лучше метода Фишера в среднем на 12 % для одного эксперимента, на 6 % для трех экспериментов и на 20 % для 49 экспериментов. В случае 10 экспериментов Recall(&)200 предлагаемого метода и метода Фишера не отличался (рис. 5, Ь).

Генетические исследовании в Российской популяции

Состоялся релиз первой базы данных GWAS исследований в России (Биобанк России) [241. Метод Фишера и предлагаемый метод были применены для такого фенотипа, как генетические риски пристрастия к курению из Биобанка России1. 175 известных гена для данного фенотипа были взяты из GWAS-каталога-, В качестве известных генов активного модуля для предлагаемого метода были выбраны 9 случайных генов из 175 известных. Остальные гены использовались для валидации предлагаемого метода. В качестве биологического графа применен граф белок-белковых взаимодействий InWeblM. ROC AUC для предлагаемого метода не отличался от метода Фншсра для одного эксперимента. Rccall@200 для предлагаемого метода был лучше метода Фншера в среднем на 11 % для трех экспериментов, на 9 % для 10 экспериментов и на б % для 49 экспериментов (рис. 6, Ь).

Генетические исследовании кожной меланомы

В [25] описаны результаты исследования, в котором оценивалось влияние экспрессии каждого отдельного гена на выживаемость пациентов с кожной меланомой в двух независимых группах. К этим результатам были применены метод Фишера и предлагаемый метод. В качестве известных генов использованы 11 иммунных генов из [25]. Биологический граф построен из данных

1 Биобанк России | Электронный ресурс]. Режим доступа: https://biobanknis.almazovcentre.ni/pheno/SMNE, свободный. Яз. рус. (дата обращения: 15.12.2024).

2 [Электронный ресурс]. Режим доступа: https://www.ebrac.uk/gwas/, свободный. Яз. англ. (дата обращения: 15.12,2024).

экспрессии генов методом семплирования байесовских сетей [26].

Для того чтобы оценить правильность полученного активного модуля, гены были ранжированы по предсказанным значениям принадлежности к активному модулю предлагаемым методом и методом Фишера независимо. Затем к первой 1000 генов из каждого списка независимо применен метод обогащения биологических путей [27]. Всего было выделено 62 биологических пути. Результаты статистической значимости для метода обогащения биологических путей на основе результатов предлагаемого метода и метода Фишера показаны на рис. 7. Каждая точка соответствует статистической значимости для метода обогащения биоло-

40

g S я

¡г Р 5Г

8 §1 9 « S

ВЗ а А г а «

¡¡I

— v с

20

щ

а ® 5 х =

5 3 ? 5 S Ь Т » 3

Sgl о

2 Р

МI" 0 20 40

—log | о (/»-значение для биологических путей, определенных па основе активного модуля, полученного предлагаемым методом)

Рис. 7. Оценка качества полученного активного модуля кожной меланомы предлагаемым методом и методом Фишера

Fig. 7. The quality of the obtained active module of skin melanoma by the proposed method and the Fisher method

гическнх путей. Смещение точек вправо относительно диагональной кривой (рис. 7, пунктирная кривая), свидетельствует о том. что результаты метода обогащения биологических путей, полученные на основании результатов предлагаемого метода, статистически надежнее, чем на основании результатов, полученных методом Фишера.

Заключение

В работе предложен метод определения активных модулей в биологических графах с использованием многокомпонентных весов вершин. Этот метод основан на последовательном применении метода Фишера для интеграции данных нескольких биологических экспериментов. метода Монте-Карло по схеме марковских цепей для генерации вероятностных выборок связных подграфов, и машинного обучения — градиентного бустинга для учета топологии графа.

Экспериментальная оценка, выполненная на симулированных активных модулях графа белок-белковых взаимодействий 1п\¥сЫМ, а также на независимых биологических графах, таких как граф реконструкции мозга BigBra¡n и генный граф вида живых организ-

мов — СаепоНшМШь е^аня, продемонстрировала значительное повышение точности предлагаемого метода по сравнению с методом Фишера, который не учитывает топологию графа. В частности, было пока!ано, что точность определения активных модулей предлагаемым методом увеличивается с ростом числа учитываемых биологических экспериментов.

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

Важным ограничением предложенного метода является предположение о независимости статистик из разных биологических экспериментов. Если между экспериментами существует сильная корреляция (например. одинаковые образцы или пересекающиеся данные), классический метод Фишера может переоценивать значимость итоговых результатов. В подобных случаях необходимы дополнительные корректировки или расширенные подходы к объединению статистик, учитывающие наличие корреляций между экспериментами.

Литература

1. Wang S.. Wu R.. Lu J.. Jiang Y„ Huang Т.. Cai YD. Protein-protein interaction networks us miners of biological discovery // Proteomics. 2022. V. 22. N 15-16. P. e2100190. https://doi.org/10.1002/ pmic.202100190

2. Rao X., Dixon R.A. Co-expression networks for plant biology: why and how //Acta Biochimica el Biophysica Sinica. 2019. V. 51. N 1IJ. P. 981 988. https://doi.org/10.1093>'abbs/gmz080

3. Rawls K-. Dougherty B.V.. Papin J. Metabolic network reconstructions to predict drug targets and off-lurget effects Methods in Molecular Biology. 2020. V. 2(188. P. 315 330. https:// doi.org/10.1007/978-1 -0716-0159-4J 4

4. Dittrich M.T.. Klau G.W.. Rosenwald A., Dandekar Т.. Mailer T. Identifying functional modules in protein-protein interaction networks: an integrated exact approach // Bioinformatics. 2008. V. 24. N 13. P. ¡223-i231 https://doi.org/Hi.1093/bioinfonnatics/btnl61

5. Zhu Q.M.. Hsu Y.H.H.. Lassen F.H., MacDonald B.T., Slead S.. Malolepsza E., Kim A., l.i T-, Mizoguchi Т., Schenone M.. Guzman G., Tancnbaum В., Fornelos N.. Carr S.A., Gupta R.M.. I:Minor P.T., Lage K. Protein interaction networks in the vasculature prioritize genes and pathways underlying coronary artery disease // Communications Biologv. 2024. V. 7. N 1. P. 87. https://doi. org/10,1038/s42003-023-05705-1

6. Nehme R„ Pietil&nen O., Artomov M.. Tegtmeyer M.. Valakh V., Lehtonen J... Bell C., Singh Т., Trehan A., Sherwood J. et. al. "1Ъе 22ql].2 region regulates presynaptic gene-products linked to schizophrenia // Nature Communications. 2022. V. 13. N 1. P. 3690. https://doi.org/10.1038/S41467-022-31436-8

7. Alcxeev N.. Jsomurodov J.. Sukhov V., Korotkevich G., Sergushiehev A. Markov chain Monte Carlo for active module identification problem // BMC Bioinformatics. 2020. V. 21. Suppl. 6. P. 261. hllps://doi.org/10.1186/s 12859-020-03572-9

8. Усольцев Д.А., Молотков И.И., Артемов H.H., Сергушичев А.А.. Шалыто А.А. Применение марковских цепей Монте-Карло и машинного обучении для поиска активною модуля в биологических графах // Научно-технический вестник информационных технологий, механики и оптики. 2024. Т. 24. № 6. С. 962-971. https://doi.org/10.17586/2226-1494-2024-24-6-962-97I

9. Kim Т.К.. Park .1.11. More about the basic assumptions oft-test: normality and sample size // Korean Journal of Anesthesiology. 2019. V. 72. N 4. P. 331-335. https://doi.oig/10,4097/kj«.d.I8.00292

References

1. Wang S„ Wu R.. Lu J.. Jiang Y. Huang T„ Cai Y.D. Protein-protein interaction networks as miners of biological discovery. Proteomics, 2022, vol. 22, no. 15-16, 1'. e2100190. hltps://doi.org/10.1002/ pmic.202100190

2. Rao X., Dixon R.A. Co-expression networks for plant biology: why and how. Acta Biochimica el Biophysica Sinica, 2019. vol. 51, no. 10. pp. 981 988. htlps://doi.org/l 0.1093/abbs/gmz080

3. Rawls K.. Dougherty B.V.. Papin J. Metabolic network reconstructions to predict drug targets and off-target effects. Methods in Molecular Biolog)'. 2020, vol. 2088, pp. 315 330. hllps://doi. org/10.1007/978-1-0716-0159-4_ 14

4. Dittrich M.T.. Klau G.W., Rosenwald A.. Dandekar T., MUller T Identifying functional modules in protein-protein interaction networks: an integrated exact approach. Bioinformatics. 2008, vol. 24, no. 13. pp. ¡223-Í231. https://doi.org/10.t093/bioinformalics/btnl6l

5. Zhu Q.M.. Hsu Y.H.H., Lassen F.H., MacDonald B.T., Stead S„ Malolepsza E., Kim A., Li T-, Mizoguchi T., Schenone M-, Guzman G., Tanenbaum B.. Fornelos N.. Carr S.A. Gupta R.M . Ellinor P.T., Lage K. Protein interaction networks in the vasculature prioritize genes and pathways underlying coronary artery disease. Communications Biology, 2024, vol. 7, no, 1, pp. 87. hltps://doi. org/10.1038/s42003-023-05 705-1

6. Nehme R., Pietiláinen O.. Artomov M.. Tegtmeyer M.. Valakh V., Lehtonen L., Bell C\, Singh T„ Trehan A., Sherwood J. el. al. The 22qll.2 region regulates presynaptic gene-products linked to schizophrenia. Nature Communications, 2022, vol. 13, no. 1. pp. 3690. https://doi.org/10.I038/s4I467-022-31436-8

7. Alexeev N.. Isomurodov J., Sukhov V.. Korotkevich G-. Sergushiehev A. Markov chain Monte Carlo for active module identification problem BMC Bioinformatics. 2020, vol. 21, Suppl. 6, pp. 261. https://doi.org/10.1186/sl2859-020-03572-9

8. Usoltscv D.A., Molotkov LI., Artomov M.N., Sergushiehev A.A., Shalyto A.A. Application of Mariiov chain Monte Carlo and machine learning for identifying active modules in biological graphs. Scientific and Technical Journal of Information Technologies. Mechanics and Optics, 2024, vol. 24, no. 6. pp. 962-971. (in Russian). hnps://doi. org/10.17586/2226-1494-2024-24-6-962-971

9. Kim T.K., Park J.H. More about the basic assumptions of l-lest: normality and sample size Korean Journal o/'Anesthesiolog)', 2019, vol. 72, no. 4. pp. 331-335. https://doi.Org/10.4097/kja.d, 18.00292

11

12,

10. Pounds S.. Morris S.W. Estimating the occurrence of false positives and false negatives in microarray studies by approximating and partitioning the empirical distribution of p-values // Bioinformatics. 2003. V. 19. N 10. P. 1 236 -1242. https://doi.org/10.1093/ bioinformatics/btgl48

Ham H., Park T. Combining p-values from various statistical methods for microbiome data // Frontiers in Microbiology. 2022. V. 13. P. 990870. https://doi.org/10.3389/fmicb.2022.990870 Li T„ Wernersson R.. Hansen R.B., Horn H., Mercer J., Slodkowicz G., Workman C.T., Rigrna O.. Rapacki K, Sterfeldt H.H., Brunak S., Jensen T.S., Lage K. A scored human protein-protein interaction network to catalyze genomic interpretation // Nature Methods. 2017. V. 14. N 1. P. 61 -64. https://doi.org/10.1038/nmeth.4083

13. Rossi R„ Ahmed N. The network data repository with interactive graph analytics and visualization // Proc. of the 29th AAAI Conference on Artificial Intelligence. 2015. V. 29. N 1. https://doi.org/10.1609/ aaai.v29i 1,9277

14. Amunts K., Lepage C., Borgeat L., Mohlberg H.. Dickscheid T., Rousseau M.E., Bludau S., Bazin P.L., Lewis L.B., Oros-Peusquens A.M.. Shah N.J., Lippert T.. Zilles K, Evans A.C. BigBrain: an ultrahigh-resolution 3D human brain model // Science. 2013. V. 340. N6139. P. 1472 1475. https://doi.org/10.1126/ science. 1235381

15. Cho A., Shin J.. Hwang S.. Kim C„ Shim H., Kim H„ Kim H.. Lee I. WormNet v3: a network-assisted hypothesis-generating server for Caenorhabditis elegans // Nucleic Acids Research. 2014. V. 42. N W1. P. W76-W82. https://doi.org/10.1093/nar/gku367

16. Zhu Z., Zhang F., Hu H., Bakshi A.. Robinson M.R., Powell J.E., Montgomery G.W., Cioddard M.F.., Wray N.R., Visschcr P.M., Yang J Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets // Nature Genetics. 2016. V. 48. N 5. P. 481 487. https://doi.org/10.l038/ng.3538

17. Usoltsev D., Molotkov I., Artomov M. A meta-predictor for causal gene identification in GWAS overcomes limitations of existing computational approaches // American Society of Human Genetics (Complex Traits and Polygenic Disorders Poster Friday Session). 2024.

18. Pardinas AF., Holmans P., Pocklington A.J., Escott-Price V., Ripke S., Carrera N.. I.egge S.E., Bishop S.. Cameron D.. Hamshere M.L., et al. Common schizophrenia alleles are enriched in mutation-intolerant genes and in regions under strong background selection // Nature Genetics. 2018. V. 50. N 3. P. 381 389. https://doi.org/10.1038/ s41588-018-0059-2

19. Barbeira A.N., Dickinson S.P., Bonazzola R., Zheng J., Wheeler H.E.. Torres I.M., Torstenson E.S.. Shah K.P., Garcia T., Edwards T.L., Stahl E.A., Huckins L.M.,Nicolae D.L., Cox N.J., Im II.K Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics // Nature Communications. 2018. V. 9. N 1. P. 1825. https://doi.org/10.1038/ s41467-018-03621-1

20. Urbut S.M., Wang G., Carbonetto P.. Stephens M. Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions // Nature Genetics. 2019. V. 51. N 1. P. 187 195. https://doi.org/10.1038/s41588-018-0268-8

21. Kolosov N., Daly M.J.. Artomov M. Prioritization of disease genes from GWAS using ensemble-based positive-unlabeled learning // European Journal of Human Genetics. 2021. V. 29. N 10. P. 1527

1535. https://d0i.0rg/10.1038/s41431-021-00930-w

22. Lam M.. Chen C-Y., Li Z.. Martin A.R.. Bryois J., Ma X.. Gaspar H., Ikcda M.. Benyamin B., Brown B.C. et al. Comparative genetic architectures of schizophrenia in East Asian and European populations// Nature Genetics. 2019. V. 51. N 12. P. 1670-1678. https://doi.org/10.1038/s41588-019-0512-x

23. Singh T„ Poterba T„ Curtis D„ Akil II., Al Eissa M„ Barchas J.D., Bass N., Bigdeli T.B.. Breen G., Bromet E.J. et al. Rare coding variants in ten genes confer substantial risk for schizophrenia // Nature. 2022. V. 604. N 7906. P. 509-516. https://doi.org/10.1038/ s4|5sf.-(C2-04556-w

24. Usoltsev D.. Kolosov N.. Rotar O., Loboda A.. Boyarinova M.. Moguchaya F... Kolesova E., Erina A., Tolkunova K.. Rezapova V., Molotkov I. et al. Complex trait susceptibilities and population diversity in a sample of 4,145 Russians // Nature Communications. 2024. V. 15. N1. P. 6212. https://doi.org/10.1038/s41467-024-50304-I

25. Usoltsev D., Njauw C.N., Ji Z., Kumar R., Sergushichev A., Zhang S., Shlyakhto E., Daly M. J., Artomov M., Tsao H. Analysis of variants induced by combined ex vivo irradiation and in vivo tumorigenesis

10. Pounds S., Morris S.W. Estimating the occurrence of false positives and false negatives in microarray studies by approximating and partitioning the empirical distribution of p-values. Bioinformatics, 2003, vol. 19, no. 10, pp. 1236-1242 https://doi.org/10.1093/ bioinformatics/btg 148

11. Ham H., Park T. Combining p-values from various statistical methods for microbiome data. Frontiers in Micmbiology, 2022, vol. 13, pp. 990870. https://doi.org/10.3389/fmicb.2022.990870

12. Li T., Wernersson R., Hansen R.B., Horn H., Mercer J.. Slodkowicz G., Workman C.T., Rigina O . Rapacki K., Stasrfeldt H.H.. Brunak S„ Jensen T.S., Lage K. A scored human protein-protein interaction network to catalyze genomic interpretation. Nature Methods, 2017, vol. 14, no. 1, pp. 61 64. https://doi.org/10.1038/nmeth.4083

13. Rossi R., Ahmed N. fhe network data repository with interactive graph analytics and visualization. Proc. of the 29,h AAAI Conference on Artificial Intelligence, 2015, vol. 29, no. 1. https://doi.org/10.1609/ aaai.v29i 1.9277

14. Amunts K., Lepage C., Borgeat L., Mohlberg H„ Dickscheid T., Rousseau M.E., Bludau S., Bazin P.L., Lewis L.B., Oros-Peusquens A.M., Shah N.J.. Lippert T., Zilles K., Evans A.C. BigBrain: an ultrahigh-resolution 3D human brain model. Science, 2013, vol. 340, no. 6139, pp. 1472-1475. https://doi.org/10.1126/ science. 12353 81

15. Cho A., Shin J.. Hwang S„ Kim C., Shim H„ Kim H„ Kim H.. Lee I. WormNet v3: a network-assisted hypothesis-generating server for Caenorhabditis elegans. Nucleic Acids Research, 2014, vol. 42, no. Wl, pp. W76-W82. https://doi.org/10.1093/nar/gku367

16. Zhu Z„ Zhang F., Hu H.. Bakshi A.. Robinson M.R.. Powell J.F.., Montgomery G.W., Cioddard M.E., Wray N.R., Visschcr P.M., Yang J. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets, Natwe Genetics, 2016, vol. 48, no. 5. pp. 481 487. https://doi.oiE/10.1038/ng.3538

17. Usoltsev D., Molotkov I., Artomov M. A meta-predictor for causal gene identification in GWAS overcomes limitations of existing computational approaches. American Society of Human Genetics {Complex Traits and Polygenic Disorders Poster Friday Session), 2024.

18. Pardinas A.F., Holmans P. Pocklington A.J., Escott-Price V., Ripke S., Carrera N.. I.egge S.E., Bishop S.. Cameron D., Hamshere M.I.., et. al. Common schizophrenia alleles are enriched in mutation-intolerant genes and in regions under strong background selection. Nature Genetics, 2018, vol. 50. no. 3, pp. 381 389. https://doi.org/10.1038/ s415 88-018-0059-2

19. Barbeira A.N., Dickinson S.P., Bonazzola R„ Zheng J.. Wheeler H.E., Torres J.M., Torstenson E.S., Shah K.P., Garcia T., Edwards T.L., Stahl E.A., Huckins L.M., Nicolae D.L.. Cox N.J., Im II.K. Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nature Communications, 2018. vol. 9, no. 1, pp. 1825. https://doi. org/10.1038/s41467-018-03621-1

20. Urbut S.M., Wang G., Carbonetto P., Stephens M. Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nature Genetics, 2019, vol. 51, no. 1, pp. 187

195. https://doi.org/10.1038/s41588-018-0268-8

21. Kolosov N.. Daly M.J., Artomov M. Prioritization of disease genes from GWAS using ensemble-based positive-unlabeled learning. European Journal of Human Genetics, 2021, vol. 29, no. 10. pp. 1527-1535. https://doi.orB/10.1038/s41431-021-00930-vv

22. Lam M.. Chen C-Y., Li Z.. Martin A.R., Bryois J.. Ma X., Gaspar H., Ikeda M., Benyamin B.. Brown B.C. et al. Comparative genetic architectures of schizophrenia in East Asian and European populations. Nature Genetics, 2019, vol. 5L no. 12, pp. 1670-1678. https://doi.org/10.1038/s41588-019-0512-x

23. Singh T., Poterba T., Curtis D„ Akil II., Al Eissa M„ Barchas J.D., Bass N., Bigdeli T.B., Breen G., Bromet E.J. et al. Rare coding variants in ten genes confer substantial risk for schizophrenia. Nature, 2022, vol. 604. no. 7906. pp. 509-516. https://doi.org/10.1038/ s41586-02 2-04550-w

24. Usoltsev D.. Kolosov N., Rotar O.. Loboda A.. Boyarinova M., Moguchaya E., Kolesova E., Erina A.. Tolkunova K. Rezapova V., Molotkov I. et al. Complex trait susceptibilities and population diversity in a sample of 4,145 Russians. Nature Communications. 2024, vol. 15, no. 1. pp. 6212. https://doi.org/10.1038/s41467-024-50304-l

25. Usoltsev D., Njauw C.N., Ji Z., Kumar R., Sergushichev A., Zhang S., Shlyakhto E., Daly M.J., Artomov M., Tsao H. Analysis of variants induced by combined ex vivo irradiation and in vivo tumorigenesis

suggests a role for the ZNF831 p.R1393Q variant in cutaneous melanoma development // Journal of Investigative Dermatology. 2024. In Press, Corrected Proof. https://doi.Org/10.1016/j. jid.2024.08.042

26. Лобода А.Л. Метод графовой кластеризации для совместного анализа данных генотипирования и экспрессии генов: диссертация на соискание ученой степени кандидата технических наук. СПб., 2022, 232 с.

27. Subramanian A., Tamayo P., Mootha V.K., Mukheijee S„ Ebert B.L., Gillette М.Л., Paulovich А», Pomeroy S.L., Golub TJC, Lander E.S., Mesirov J.P. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles // Proc. of the National Academy of Sciences of the United States of America. 2005. V. 102. N 43. P. 15545 15550. https://doi.org/10.1073/ pnas.0506580102

suggests a role for the ZNF831 p.R1393Q variant in cutaneous melanoma development. Journal of Investigative Derma tolog\\ 2024. in Press, corrected proof, httpv //doi.Qrg/10.1016/j,jid.2024.0S .042

26. Loboda A.A. A method of graphical clustering for joint analysis of genotyping and expression data. Dissertation for the degree of candidate of technical sciences. St. Petersburg, 2022. 232 p. (in Russian)

27. Subramanian A., Tamayo P., Mootha V.K., Mukherjee S., Fbert B.L., Gillette M.A., Paulovich A., Pomeroy S.L., Golub T.R., Lander E.S., Mesirov J.P. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. of the National Academy of Sciences of the United States of America. 2005. vol. 102, no. 43. pp. 15545 15550. https://doi.oig/10.1073/ pnas.0506580102

Авторы

Усольнсв Дмитрий Андреевич старший научный сотрудник, Институт геномной медицины, Детская больница Пейшенвайд, Колумбус, 43205, США: аспирант, Университет ИТМО. Санкт-Петербург, 197101, Российская Федерация, se 57279360300, https:// orcid.org/0000-0001 -8072-31 OX, dusoltsev.27@gmail .com Молотков Иван Нгирсвич старший научный сотрудник. Институт геномной медицины. Детская больница Нейшенвайд, Колумбус, 43205. CIIIA: аспирант. Медицинский колледж Университета штата Огайо, Колумбус, 43210, США, se 58651494600. https://orcid.org/0009-0008-3566-0160, [e-mail] Артемов Никита Николаевич кандидат химических наук, доцент. главный исследователь. Институт геномной медицины. Детская больница Нейшенвайд, Колумбус, 43205, США: профессор педиатрии, Медицинский колледж Университета штага Огайо, Колумбус, 43210. США, se 36542095500, https://orcid.org/0000-0001-5282-8764, mykyta.artomovfíi'nationwidechildrens.org

C'epiyiUHHeB Алексей Александрович кандидат технических наук, доцент. Университет Вашингтона в Септ-Луисе. Сент-Луис. 63110, CHIA, se 55772694000, https://orcid.org/0000-0003-1159-7220, asergushichev (а>wustl.edu

Шалыто Анатолий Абрамович доктор технических наук, профессор, главный научный сотрудник, профессор. Университет ИТМО, Санкт-Петербург. 197101. Российская Федерация, se 56131789500. https://orcid.org/0000-0002-2723-2077, anatoly.shaIyto(a!gmail.com

Authors

Dmitrii A. Usoltsev Senior Researcher, Institute for Genomic Medicine, Nationwide Children's Hospital, Columbus, 43205, USA; PhD Student, ITMO University, Saint Petersburg, 197101, Russian Federation, se 57279360300, https://orcid.org/0000-(>001 -8072-310X, dusoltsev.27@ gmail.com

Ivan I. Molotkov - Senior Researcher, Institute for Genomic Medicine. Nationwide Children's Hospital, Columbus, 43205, USA; PhD Student. The Ohio State University College of Medicine, Columbus. 43210, USA. sq 58651494600, https://orcid.org/0009-0008-3566-0160, ivan. molotkov(S>nati onwidechildrens.org

Mykyta N. Artomov PhD (Chemistry), Associate Professor, Chief Researcher, Institute for Genomic Medicine, Nationwide Children's Hospital, Columbus, 43205. USA; Professor of Pediatries. The Ohio Slate University College of Medicine, Columbus, 43210, USA, sc 36542095500, https://orcid.org/0000-0001-5282-8764, mykyta. artomov fa1 nat ionwi dech i ldrens. org

Alexcy A. Sergushichev — PhD, Associate Professor, Washington University School of Medicine in St. Louis, St. Louis, 63110, USA. so 55772694000, https://orcid.org/0000-0003-l 159-7220, asergushichevfa wustl.edu

Anatoly A. Shalyto D.Se., Chief Researcher, Full Professor, ITMO University. Saint Petersburg, 197101, Russian Federation, se 56131789500." https://orcid.org/0000-0002-2723-2077, anatoly. shalyto(ogmai 1 .com

Статья поступила в редакцию 09.02.2025 Одобрена после репетирования ¡4.04.2025 Принята к печати 26.05.2025

Received 09.02.2025

Approved after reviewing 14.04.2025

Accepted 26.05.2025

Работа доступна по лицензии Creative Commons »Attribution-NonCommercial"

nature communications

8

Article

https://doi.org/10.1038/s41467-024-50304-1

Complex trait susceptibilities and population diversity in a sample of 4,145 Russians

Received: 27 March 2023

Accepted: 2 July 2024

Published online: 23 July 2024

* Check for updates

Dmitrii Usoltsev ©'-2'3'4'5, Nikita Kolosov'-2-3'45, Oxana Rotar1, Alexander Loboda1'2'3, Maria Boyarinova1, Ekaterina Moguchaya1, Ekaterina Kolesova1, Anastasia Erina1, Kristina Tolkunova1, Valeriia Rezapova1'2'3, Ivan Molotkov4,5, Olesya Melnik1, Olga Freylikhman©1, Nadezhda Paskar1, Asiiat Alieva1, Elena Baranova1, Elena Bazhenova1, Olga Beliaeva1, Elena Vasilyeva1, Sofia Kibkalo1, Rostislav Skitchenko©1, Alina Babenko1, Alexey Sergushichev©2, Alena Dushina6, Ekaterina Lopina6, Irina Basyrova6, Roman Libis6, Dmitrii Duplyakov ©7,a, Natalya Cherepanova8, Kati Donner9, Paivi Laiho©10, Anna Kostareva1'2, Alexandra Konradi©1-2, Evgeny Shlyakhto1, Aarno Palotie©3'9,11, Mark J. Daly©39'11 & Mykyta Artomov•WMM11

The population of Russia consists of more than 150 local ethnicities. The ethnic diversity and geographic origins, which extend from eastern Europe to Asia, make the population uniquely positioned to investigate the shared properties of inherited disease risks between European and Asian ancestries. We present the analysis of genetic and phenotypic data from a cohort of 4,145 individuals collected in three metro areas in western Russia. We show the presence of multiple admixed genetic ancestry clusters spanning from primarily European to Asian and high identity-by-descent sharing with the Finnish population. As a result, there was notable enrichment of Finnish-specific variants in Russia. We illustrate the utility of Russian-descent cohorts for discovery of novel population-specific genetic associations, as well as replication of previously identified associations that were thought to be population-specific in other cohorts. Finally, we provide access to a database of allele frequencies and GWAS results for 464 phenotypes.

Linking inherited DNA variation to the disease risks is one of the main goals in modern predictive medicine. Large-scale projects such as the UK Biobank1, FinnCen'' and Biobank Japan' have made a substantial contribution to the understanding of human biology and the advancement of personalized medicine. The increasing ethnic diversity of genetic studies resulted in the discovery of population-specific

susceptibility loci that could not be identified in other ancestries'. Inclusion of understudied populations into biobank initiatives improves fine-mapping accuracy of previously identified GWAS signals and novel risk gene discovery efforts'.

Genetic studies in the multinational Russian population, which is geographically located at the crossroads of Europe and Asia, could

'Almazov National Medical Research Centre, St Petersburg, Russia. 2ITMO University, St Petersburg, Russia. 3Broad Institute, Cambridge, MA USA. 4The Institute for Genomic Medicine, Nationwide Children's Hospital, Columbus, OH, USA. department of Pediatrics, The Ohio State University College of Medicine, Columbus, OH, USA. 60renburg State Medical University, Orenburg, Russia. 7Samara State Medical University, Samara, Russia. 8Samara Regional Cardiology Dispensary, Samara. Russia, institute for Molecular Medicine Finland (FIMM), Helsinki. Finland. ,0Finnish Institute for Health and Welfare (THL), Helsinki, Finland. "Analytic and Translational Genetics Unit, Massachusetts General Hospital, Boston, MA, USA. e-mail: mykyta [e-mail]

have a potential power to detect historic origins of population-specific variants. Additionally, genetic data from the Russian population could be a powerful source of replication for population-specific associations found in three largest blobanks, UK biobank'. FinnGen\ and Biobank Japan'.

Historically, polygenic trait genetics has been omitted in Russia, resulting in a lack of GWAS studies based on local cohorts. Russian descent individuals were involved primarily either as a small part of consortia datasets or as a basis for population genetic studies lacking phenotypic information"10.

At the same time, there are more than 150 local ethnicities in Russia that would greatly benefit from the local large genetic variation studies. This becomes especially important in light of the lack of transferability of polygenic risk score models between ancestries" r'. The aforementioned variety of ethnicities represents the genetic history of populations between Eastern Europeans, Finns, and Asians. For example, a study of whole genome sequencing data previously showed that populations in the northern regions of Russia, west of the Ural Mountains, not only belong to a Finno-Ugric language group but are also genetically close to Finns, while populations in the central regions of western Russia showed similarities with eastern Europeans". Furthermore, close genetic relation between north-western Russians and Finns was also observed in Balto-Slavic speaking populations comparison1'. In addition, another study showed that Siberian populations separated from other East Asian populations 8800-11,200 years ago and significantly contributed to the formation of Eastern European populations 4700-8000 years ago1". Thus, a gene flow from Asia through the Ural Mountains to eastern Europe was hypothesized. Consistently, a notable genetic relationship between Finns and Mongolian tribes was observed'. Additional evidence of the great diversity of the Russian population linking European, Asian, and Native American populations was presented in a country-wide study -Genome Russia Project, yet no phenotypic information was collected at that time resulting in the lack of biomedical applications for this data"1.

For studies reflecting population history, a small sample set (-100 individuals) is usually sufficient; however, personalized medicine, biobank assembly and GWAS require a larger sample size and extensive data collection, which is impossible without epidemiological studies. In 2012-2013 in 12 regions of Russia a national study "Epidemiology of cardiovascular diseases in different regions of the Russian Federation" (ESSE-RF) was launched. Within the framework of this study, a stratified multistage random sample of approximately 22,000 residents with blood biobanking and detailed phenotyping was collected".

Here, we present the first results of the analysis of clinical and genetic data from three metro areas that participated in the ESSE-RF study: St. Petersburg, Samara, and Orenburg. Our results reflect the genetic structure of the Western Russian population and phenotypic susceptibilities from 4145 participants across 464 phenotypes.

Results

Data collection

A cohort of 4800 residents of three areas in Russia - St. Petersburg UV=1600. N males = 573, N females = 1027), Orenburg (/V = 1600, N males = 656, N females = 944) and Samara (N= 1600, N males = 697, N females = 903) were recruited in 2012-2013 through an ambulatory visit to local hospitals and polyclinics (aged 461 USD years) (Fig. la). Detailed phenotypic data and a blood sample were collected. All participants signed an informed consent. The study protocol was approved by local ethics committees (Almazov National Medical Research Center, St. Petersburg) (Methods, Data collection: Supplementary Note 1).

Each patient was invited for an ambulatory visit for one day to collect phenotypic information (Fig. lb). In 2018-2019,289 out of 1600

original patients from St. Petersburg were invited for an additional ambulatory visit as a part of different local studies (familial hypercholesterolemia, metabolically healthy obesity, premature vascular aging). These patients were subjected to a detailed follow-up data collection protocol (Methods, Data collection: Supplementary Note 2).

For all participants, biannual updates on phenotypic and vital status were recorded through direct contact (phone calls) and indirect contacts (mail/e-mail letters, information from local clinical district databases. Supplementary Data 1). Additional independent cohort of 138 samples were recruited in 2017-2018 for participation as controls in a local study of early childhood starvation effects. They were evaluated using a similar clinical short protocol, and no follow-up data collection was performed (Supplementary Note 3).

Genetic data generation and imputation

A total of 4723 individuals, comprising 4594 population samples (ESSE) and 129 starvation controls were genotyped using FinnGen ThermoFisher Axiom microarray'. Genotype imputation was conducted using Haplotype Reference Consortium (HRC) panel1' '" and Beagle 5.2 software" (Methods. Genetic data generation and imputation). The quality of the imputation was assessed using a masking experiment (Supplementary Note 4.1, Supplementary Fig. 1).

Samples with sex-mismatch (224 ESSE and 28 Starvation Controls) (Supplementary Materials, Supplementary Note 4.2, Supplementary Fig. 2) and duplicated samples (190) were excluded. Additionally. 347 related individuals were marked for further analysis (Supplementary Note 4.3, Supplementary Fig. 3).

The final dataset included 4281 individuals (4183 ESSE and 98 Starvation Controls) and 11,077,763 variants. Overall, 37,439 variants failed Hardy-Weinberg equilibrium (p < 1 x 10 4) and 371 variants were discordant with the HRC imputation panel. Allele frequencies were compared against gnomAD Finnish and non-Finnish European populations (Supplementary Note 4.4, Supplementary Fig. 4).

Population structure analysis

Initially, we performed principal component analysis (PCA) to identify ancestral clusters within the dataset (Methods; Population structure analysis; Fig. 2a; Supplementary Note 5; Supplementary Figs. 5 and 6). Joint PCA with the 1000 Genomes dataset indicated that the western Russian population is represented by individuals of admixed ancestry, spanning from European to East Asian continental ancestry (Fig. 2b. Supplementary Fig. 7a). A separate joint PCA with only European subpopulations from 1000 Genomes demonstrated closer relatedness of Russians to Finnish population. This relationship was particularly evident in the overlap observed at the PC1:PC2 plane (Supplementary Note 6, Fig. 2c. Supplementary Fig. 7b). We were able to distinguish between Russians and Finnish populations at higher principal components (PC3-PC4; Supplementary Fig. 8). Six clusters within the Russian data set were Identified in the PCA space (Supplementary Note 7, Fig. 2d, Supplementary Fig. 9).

We calculated the contribution of haplotypes from geographically proximate ancestries represented in 1000 Genomes to each cluster using the ADMIXTURE" (Supplementary Note 8, Fig. 2e, Supplementary Fig. 10). The first cluster included mostly individuals related to southern Europeans. The second cluster predominantly included a population close to Utah residents (CEU) with Northern and Western European ancestry and the British with a small admixture of Finnish and Asian ancestry. In the next four clusters, the proportion of Finnish and Asian haplotypes increased and the proportion Northern and Western European ancestry haplotypes decreased, which was consistent with the location of these clusters in the principal components space; these clusters co-localized with the Finnish population in PC1-PC2 (Supplementary Note 6, Fig. 2c). Furthermore, in order to establish the independence of our findings from the genotyping platform, we conducted an ADMIXTURE analysis exclusively utilizing HapMap

d 0.12

REGION

• Orenburg, n=1438

• Samara. n=1214

• St. Petersburg. n= 1493

POPULATION =>

• AFR, n=661 AMR. n-34?

Iii

• EAS. n=504 —

• EUR, n=503 O

D.

• RUS, n=4145

• SAS. n=4B9

POPULATION

• CELLGBR, n:

0.025 0 050 0 075 0.100 PC1. (EV=9.6B. 0.23%)

PC1. (EV=110.12, 1.63%)

PC1, (EV=10.47, 0.22%)

0.025 0.050 0.075 0 100 PC1, (EV=9.68, 0.23%)

Fig. 21 Population structure analysis. The followingabbreviations were used: AFR Africans (ACB African Caribbean in Barbados, ASW African Ancestry in SW USA, ESN Esan in Nigeria, CWD Gambian in Western Division - Mandinka, LWK Luhya in Webuye, Kenya: MSL Mende in Sierra Leone: YRI Yoruba in tbadan, Nigeria), AMR native Americans (CLM Colombian in Medellin, Colombia, MXL Mexican Ancestry in Los Angeles CA USA, PEL Peruvian in Lima Peru, PUR Puerto Rican in Puerto Rico), EAS east Asians (CHB Han Chinese, CHS Han Chinese South, JPT Japanese in Tokyo, CDX Chinese Dai in Xishuangbanna, KHV Kinh in Vietnam), RUS Russians, &45south Asians (BEB Bengali in Bangladesh, GIH Gujarati Indians in Houston, Texas, USA: ITU Indian Telugu in the U.K.: PJL Punjabi in Lahore, Pakistan; STU Sri Lankan Tamil

in the UK), CUR Europeans (CEU Utah residents (CEPH) with Northern and Western European ancestry, GBR British, FIN Finnish in Finland, IBS Iberian populations in Spain, TSI Toscani in Italy), a Principal component analysis with labeling indicating sample collection region: bjoint principal component analysis with 1000 Genomes cohort; c Joint principal component analysis with European subsample of 1000 Genomes cohort; d Clustering of the Russian cohort in the PCA space; e Admixture analysis for each of six populational clusters in the Russian cohort, samples are arranged with respect to their PCI coordinate; f Hierarchical clustering of the Russian cohort with 1000 Genomes subpopulations with respect to the sharing of IBD regions.

Genomes populations. We collected all IBD regions with LOD score quality greater than 3 for each pair of individuals. Next, we merged the IBD regions if the gap was not greater than 0.6 cM and if there was not more than 1 discordant homozygous variant. We summarized the length of resulting IBD regions for each pair of individuals and calculated the median IBD length between pairs of populations (Supplementary Fig. 14). The resulting heatmap for all populations is shown in Fig. 2f. According to the clustering analysis, the Finns are closer to the first five Russian clusters than to the European populations from 1000 Genomes (Supplementary Note 10).

Conclusively, the Russian population sampled in the large metro-areas represents a heterogeneous combination of individuals of admixed ancestries between European and Asian populations. In addition, for a subgroup of the sixth cluster with the highest PCI values more than half of the haplotypes are of Asian origin. This finding is expected, as these samples primarily come from Orenburg, a region located near the border with Kazakhstan (Central Asia). It has previously been shown that individuals from Central Asia have a mixture of European and Asian haplotypes-' 2\

To eliminate the potential influence of systematic errors linked to imputation, we conducted a population structure analysis solely utilizing genotyped variants. Notably, the comparison between cluster assessments based on genotyped variants and assessments involving

all variants revealed a robust correlation (R2 = 0.87). This correlation highlighted the strong alignment between these two approaches. As a result, the outcomes of our population structure analysis provided strong validation for the consistency of our observations, underscoring that imputation had no significant effect on the results (Supplementary Note 11, Supplementary Fig. 15).

Enrichment of finnish and russian variants

The analysis of population structure revealed an intricate pattern of relatedness among Russians, Finns, and East Asians. The Finnish population is historically unique; however, the genetic similarity of the Russian population illustrated above suggests that DNA variants enriched in Finns might also be found in Russia.

To assess the population-specific properties of DNA variants, we created the distributions of log2 allele frequency ratios between the target population and non-Finnish Europeans from gnomAD (NFE) for 265,624 Finnish enriched variants (Methods. Enrichment of Finnish and Russian variants; Fig. 3a, Supplementary Fig. 16). The medians of these log-ratio allele frequency distributions (excluding variants that were not observed) increased from cluster 2 (0.56) to cluster 6 (3.05). In the sixth cluster, it reached an even higher value than in the Finnish population (3.032). As the fraction of East Asian haplotypes increases from cluster 2 to cluster 6 of the Russian

which would potentially explain the presence of population-specific variants. We calculated a maximum likelihood tree using the TreeMix model with reference populations: 1000 Genomes: the populations closest to Russia from the Estonian Genome Diversity Panel (EGDP)R; and Russians from the Human Genome Diversity Panel (HGDP)r. Interestingly. Finnish samples were the closest to Russians from HGDP and cluster 4, while cluster 5 was close to Ural populations (Fig. id). The only Russian cluster within the Asian branch is cluster 6, which is further supported by its close proximity in principal component space to Central Asia (Supplementary Fig. 21). Furthermore, we performed Treemix analysis additionally using all HGDP populations (Supplementary Note 14; Supplementary Fig. 22).

Using the ADMIXTURE in unsupervised mode with 8 clusters we built clusterization consistent with our tree, and detected a genetic component (yellow) widely represented in the Siberian population. The percentage of this genetic component decreased with moving from Siberia to East Europe through the Urals. Also, this component is present in the Finnish population, which may be the result of a known gene flow from Asia to Europe through the Ural Mountains (Supplementary Fig. 23).

The Russian population also had enrichment of population-specific variants compared to the Europeans, some of which have been described previously". We identified 44.936 Russian-enriched variants as having a log-ratio between RUS and NFE from gnomAD greater than 2. Among Russian-enriched variants 40,743 (90.61%) had higher AF in East Asians than in Russians and only 2045 (4.55%) were not observed in the East Asian population in gnomAD (Supplementary Fig. 19c, d).

GWAS

We performed GWAS for 464 phenotypes. Results are available at https://biobank.almazovcentre.ru.

Although 4600 samples is a relatively modest cohort size for GWAS, we provide examples of replication of findings from other biobanks as well as newly identified associations specific to the Russian population.

Several associations observed in the UK biobank were replicated, for example, rs7412 and rs4970S34 for LDL, rs4697701 and rs4549940 for uric acid levels (Supplementary Fig. 24a. b).

Interestingly, nominal associations in the UK biobank, such as the association of rsl3266066 with the initiation of smoking (UKBB phe-nocode 20116 0, beta = -0.0043, p = 0.00022)', later confirmed with the MTAG approach using multiple addiction phenotypes (beta = -0.007, p = 1 xio '°, beta was reversed to match models)27, were significantly associated with the "never smoked" phenotype in the Russian cohort (N never smoked = 2391; AF never smoked = 0.414. N controls = 1488; AF controls = 0.475, beta = -0.28, p = 3.74 * 10"', Supplementary Fig. 24c). To reduce the possibility of technical artifacts associated with this observation, we looked at only directly genotyped variants in this locus and confirmed the presence of highly-associated rs11781072 (p = 1.45 * 10"). eQTL properties of rsl3266066 are associated with expression of PTK2 in cerebellum (p = 8.8 * 10'12). We also performed a gene prioritization analysis which indicated the putative causal role of PTK2 (Supplementary Fig. 25).

Several novel genome-wide significant associations were identified: current smoking (rs7972723, with AF of 0.1465 in RUS, 0.1447 in NFE, and 0.1680 in FIN populations; 834 individuals as current smokers, AF of 0.189, and 3045 controls with AF of 0.137; beta = 0.43, p-2.08 x 10"; Supplementary Fig. 24d), abdominal obesity (rs56046524, with AF of 0.3537 in RUS, 0.3938 in NFE, and 0.3632 in FIN populations; 1405 cases with allele frequency of 0.306, and 2462 controls with AF of 0.378; beta = -0.324, p = 3.7 « 10"; Supplementary Fig. 24e),and increased blood pressure in thesecond half of pregnancy (rsll948871, with AF of 0.2024 in RUS, 0.1727 in NFE, and 0.1479 in FIN populations; 366 cases with AF of 0.279, and 1642 controls with AF of 0.185: beta = 0.55, p = 1.4 * 10s; Supplementary Fig. 24f). However,

given the modest size of the discovery cohort, thorough replication is necessary for these findings (Supplementary Note 15).

Due to the presence of six genetic clusters in our population, we analyzed the effect sizes of the discussed variants in each cluster independently. We observed overall consistency in effect sizes. Some clusters were too small to detect a significant deviation of the effect size from 0 (Supplementary Fig. 26).

Population structure in the Russian cohort indicated that it could be feasible to use it for replication of Finnish-specific genetic associations. Despite cohort size limitations, we attempted to illustrate this through a systematic approach. First, we selected only potentially replicable (MAF RUS >0.01) Finnish-enriched variants that were less frequently observed in the East Asian population than in Finns (logratio between EAS and NFE<2). There were 20,050 such variants, which were not clumped together (includes LD-correlates). Among them, we found 142 variants with genome-wide significance associations with 177 traits in FinnGen (total 773 variant-phenotype pairs). Overlaps with the Russian cohort were found for 62 variants in 53 traits (332 variant-phenotype pairs. Supplementary Fig. 27a).

We combined the phenotypes by similarity into 8 phenotypic groups: Diabetes, Sleep apnea, Asthma, Statin, Alzheimer, Hypothyr oidism. Arthritis, Hypertension, and for each group, we independently performed LD-clumping (R2<0.1; Supplementary Data 3). The resulting data set contained 11 independent variants in 8 phenotypic groups. The only variant,rs7480 0719, passed the Bonferroni-adjusted replication threshold (p = 0.05/11 variants/8 traits = 5.7 x 10"*). This variant was associated with an increased risk of Alzheimer disease in FinnGen (p = 1.14 * 10"32, beta = 0.5718). In the Russian cohort, it was associated with an increase in the comorbid phenotype • apolipoprotein B levels (p = 4.4 x 10'4, beta = 0.15) (Supplementary Fig. 27b)".

Furthermore, we conducted an investigation of the genetic correlations between GWAS from the Russian Biobank and their corresponding counterparts in the UK Biobank and FinnGen datasets. Initially, we excluded all GWAS that had heritability estimates outside the range of [0,1] and a heritability standard error more than 50%. This left us with a total of 35 GWAS, 26 of which were matched with corresponding traits from the UK Biobank (Supplementary Data 4). For complex phenotypes, we compared them with all relevant components from the UK Biobank dataset. Thus, we had 34 pairwise comparisons. Out of 34 comparisons, 26 were nominally significant (p< 0.05). and 22 passed the Bonferroni significance threshold (0.05/ 34 = 0.00147). From the FinnGen traits, we specifically chose dyslipi-demia, hypertension, type 2 diabetes (T2D). obesity, myocardial infarction, ischaemic heart disease, anxiety, depression, smoking, and sleep apnoea. We proceeded to examine the genetic correlation between these selected traits and their corresponding counterparts from our pool of 35 traits. This led to the construction of 18 pairs, with 15 of them showing nominal significance and 8 passing the Bonferroni significance threshold (0.05/18 = 0.00278) (Supplementary Data 4).

Discussion

The Russian biobank resource presented here is an essential step towards accessibility of precision medicine for the patients with a wide variety of genetic makeups not currently represented in other major genetic studies.

Our cohort illustrates that the genetic structure of the Russian population, sampled in metropolitan areas in the European part of the country, consists of the number of subpopulations with high related-ness to Finnish and East Asian populations. We also identified a subgroup that has Central Asian origins. This subgroup exhibited the highest proportion of Asian haplotypes and represents a mixed population with significant genetic similarity to Central Asian populations. This finding is supported by the geographic distribution of the samples and previous studies on Central Asian populations^25'2"". The Finno-Ugric subpopulations in Russia are historically found west

of the Uralic mountains, which is in good agreement with previous whole genome studies conducted in this area"'14. The gene flow from Asia through Siberia and the Ural Mountains to eastern Europe, together with the previously found relationship between Finns and Asians, suggests that the unique Finnish variants could also be found on the territory of modern Russia. The bottleneck that Finnish population went through resulted in an increase of the allele frequencies of common Finno-Ugric ancestor population. Moreover, our IBD and ADMIXTURE analyses provide potential explanations to high related-ness between Finnish and Mongolian populations reported previously1.

Replication studies presented here indicate that even with a relatively modest cohort size, previously reported associations from the UK biobank1 and FinnGen; could be directly observed in the Russian cohort. Such a unique genetic structure of the Russian population provides a potential power to discover and replicate associations that often were considered population-specific. Importantly, such replication can be achieved simultaneously with relatively modest cohort size and resources.

This suggests that the susceptibility to polygenic diseases in Russia could potentially be driven by a mixture of variants from multiple ancestral populations. Some of the ancestral populations in this case have not been studied before. Assessing population-dependent contributions of many associated alleles would be critically important for creating informative polygenic risk score models for individual inherited risk evaluation'-.

We continue to monitor participants of this study with the latest data update in the Fall of 2023. Our proof-of-concept study shows that infrastructure, logistics, and research resources are sufficient to create polygenic trait studies in Russia. The major challenge, yet to be resolved, is an outline of how to scale such efforts to the size of other major biobanks.

Finally, we anticipate that the first local resource for polygenic trait genetics studies in Russia, which provides the largest public reference for allele frequencies and genetic associations - Biobank Russia (https://biobank.almazovcentre.ru) - will become a core for further expansion of complex trait genetics research to yet understudied populations.

Methods

Data collection

Overall. 4800 patients were invited for one day ambulatory visits in 2012-2013. Several types of data were collected: responses to the health and medical history questionnaire, dietary behavior, physical activity, social status, depression, anxiety and perceived stress tests, anthropometric parameters, blood pressure and heart rate measurements, blood metabolic panel test, including glucose, creatinine, uric acid, total cholesterol, high-density lipoprotein (HDL), low-density lipoprotein (LDL), triglycerides, insulin and N-terminal prohormone of brain natriuretic peptide (pro-BNP) levels measurements.

For the St. Petersburg cohort, several additional phenotypes were collected: measurements of blood pressure and heart rate in standing position, vessel stiffness measurements and electrocardiogram (ECG), urine albumin and extended blood metabolic panel test with additional test for C-reactive protein (CRP), lipoprotein (a), apolipoproteins A and B, Cortisol, leptin, adiponectin, and vitamin D. (Supplementary Note 1, Fig. lb).

Genetic data generation and imputation

DNA was extracted from blood samples of 4723 individuals (4594 population sample ESSE+ 129 Starvation Controls) using the QIAamp DNA Mini Kit (Qiagen) and genotyped using a custom FinnGen Ther-moFisher Axiom microarray'.

Genotype imputation in the Russian population was carried out with the Haplotype Reference Consortium (HRC)W reference panel in

accordance with optimal solution for the Russian-descent cohorts' (Supplementary Note 4.1). The HRC panel was preprocessed and filtered to meet the data formatting requirements for further imputation procedure. After the filtering there were 37,620,211 variants and 27,165 individuals left in HRC.

Before the imputation, genotypes of the study cohort were filtered, pre-phased, strand-checked, and split into individual chromosomes. In the end, we kept 474,430 variants and 4723 individuals for further SNP imputation.

All imputations were performed by Beagle 5.221 with the default parameters (burnin = 6, iterations = 12, imp-segment - 6, ne = 1000000). The imputation quality for each variant was measured using Dosage-R2 (DR2), as given in Beagle output. All variants with DR2 > = 0.8 were considered well-imputed and kept for further analysis. We confirmed the quality of imputation using a masking experiment (Supplementary Note 4.1 and Supplementary Fig. 1). 224 ESSE and 28 Starvation Controls samples with mismatch between detected and reported sex were excluded (Supplementary Note 4.2, Supplementary Fig, 2).

Relatedness and sample duplications were assessed for individuals using kinship analysis with PLINK2 '. 190 duplicates (187 ESSE and 3 Starvation Controls) were removed, and 347 related individuals were labeled for exclusion in further analysis (Supplementary Note 4.3, Supplementary Fig. 3).

The resulting genetic dataset consisted of 4281 (4183 ESSE and 98 Starvation Controls) individuals and 11,077,763 variants. The genetic data was subjected to quality filtration using python3 "hair (vO.2.85) package". Overall 37,439 variants failing Hardy-Weinberg equilibrium (p < 1 * 10 4) were eliminated from the analyses. Additionally, we removed 371 variants discordant from the HRC imputation panel. The rest of the variants were subjected to comparison of allele frequency against gnomAD Finnish and non-Finnish European populations (Supplementary Note 4.4, Supplementary Fig. 4).

Population structure analysis

PC A was performed with PLINK2 using535,727common autosomal LD-pruned variants (R! < 0.2). Additionally, we excluded 136 individuals from the ESSE cohort as PCA outliers using "adamethods" (vl.2.1) R package with built-in max Euclidean distances approach" " (Supplementary Note 5, Fig. 2a, Supplementary Fig. 5-6). Then, we merged the genotyping dataset of ESSE/Starvation Controls with 1000 Genomes WGS data and performed PCA on common 506,617 LD-pruned variants using the python3 "hail" (vO.2.85) package (Supplementary Note 6, Fig. 2b, Supplementary Fig. 7a). Additionally, we selected ESSE/Star vation Controls and only 1000 Genomes European populations and ran PCA on common 515,649 LD-pruned variants (Supplementary Note 6, Fig. 2c, Supplementary Figs. 7b and 8). Eigenvalues and the corresponding percentage of variance explained for all PC components were calculated using PLINK2 with the number of PCs equal to the number of individuals in a genotyping dataset.

PCA clusters of Russian samples were identified for 4145 individuals using "SVDFunctions" (ulJ) R package'' (Supplementary Note 7, Fig. 2d, Supplementary Fig. 9). Admixture analysis was performed for combined (ESSE/Starvation Controls - 1000 Genomes) LD-pruned genotype matrix excluding 361 relative individuals (347 - ESSE, 14-1000 Genomes). The final dataset for ADMIXTURE (vl.3.0) consisted of 6649 individuals and 506,617 variants" (Supplementary Note 8, Supplementary Figs. 10 and 11, and Fig. 2d).

We calculated Fst statistics using VCFtools (vO.L15)'s (Supplementary Note 9, Supplementary Figs. 12 and 13).