Поиск биологического обоснования наследуемости фенотипов является одной из основных задач современной медицинской генетики. Исторически все известные подходы к поиску патогенных геномных вариантов делят на две группы: биологические и математические. К биологическим методам относятся подходы объяснения фенотипа на основе изученных биохимических процессов, а к математическим — различные статистические подходы. В случае, когда невозможно напрямую проследить биохимическую цепь формирования фенотипа, но развитие заболевания явно отягощено семейным анамнезом, применяют различные статистические подходы: полногеномные ассоциативные исследования (GenomeWide Association Studies, GWAS) [1], методы оценки полигенных рисков (Polygenic Risk Score, PRS) [2], подходы по выделению гаплотипов [3] и другие методы. Однако у каждого из методов есть свои особенности и ограничения. При вычислении по алгоритмам GWAS полиморфизм оценивают независимо, а затем производят фильтрацию по критерию значимости. Таким образом, результат любого заболевания, являющегося совокупностью редких наследуемых мутаций, с помощью данного подхода будет некорректным, так как максимально значимыми окажутся не относящиеся к фенотипу полиморфизмы или вообще ни один полиморфизм не преодолеет критерии фильтрации. Подход PRS представляет собой расширение методологии GWAS, однако он тоже основан на предположении, что все используемые в модели однонуклеотидные варианты (single nucleotide polymorphisms, SNP) независимы. К тому же PRS-модели не воспроизводятся между популяциями при воспроизводимости локусов точек, на которых они построены, т. е. использования одиночных вариантов явно недостаточно для построения моделей и классификации пациентов. Кроме того, существует значительное количество заболеваний, генетическую структуру которых не удается определить ни биологическими, ни популярными статистическими подходами. Чаще всего в основе алгоритма наследования таких фенотипов находятся гаплотипные структуры. Под гаплотипом подразумевается линейно сцепленная совокупность из нескольких (в том числе до нескольких сотен) вариабельных вариантов, которые совместно образуют небольшое (менее 100, в среднем 10–20) количество вариантов аллелей. Подход с выявлением конкретных вариантов гаплотипов активно используют в фармакогенетике при анализе цитохромов P450. Для CYP2D6, например, существует более 120 гаплотипных вариантов, которые обусловлены более чем 500 полиморфизмами [4]. Однако подобный подход очень редко обсуждают для большинства локусов полигенных заболеваний.

Особой популярностью GWAS-исследования пользуются при изучении заболеваний с вовлечением нервной системы, в том числе нарушений развития и нейродегенеративных патологий с полигенной природой, таких как боковой амиотрофический склероз, болезнь Паркинсона, шизофрения, патологии аутического спектра. Они позволяют выявлять участки генома, изменения в которых перепредставлены у больных по сравнению с общей популяцией. Интересным свойством GWAS является возможность выявления локусов, патогенным вариантом в которых является не однонуклеотидный вариант из структуры чипа, а структурная вариация, напрямую чипом не детектируемая, но выявляемая по сцепленным с ней SNP. В частности, в GWAS бокового амиотрофического склероза детектируется локус в гене C9orf72, содержащий увеличенный шестинуклеотидный повтор G4C2 (GGGGCC) [5], при этом варианты повтора непосредственно чипом не детектируются. В GWAS болезни Хантингтона выявляется локус гена HTT в 15-й хромосоме, содержащий тринуклеотидные повторы [6], при этом в чипе нет проб, соответствующих непосредственно участку повтора.

Эндотелиальная дистрофия роговицы Фукса (ДФ) — наследственное заболевание глаз, характеризующееся уменьшением количества эндотелиальных клеток роговицы, поддерживающих водный баланс стромы роговицы. ДФ является полигенным заболеванием и представляет значительный интерес для исследования генетики [7]. Существуют две формы ДФ: ранняя и поздняя. Они различаются генетическими причинами. Ранняя форма диагностируется в возрасте до 50 лет, встречается очень редко и связана с патогенными вариантами гена COL8A2 [8]. Поздний дебют ДФ проявляется в возрасте старше 50 лет и является наиболее распространенной формой ДФ. Было показано, что поздний дебют ДФ ассоциирован с экспансией тринуклеотидных повторов CTG18.1 в интроне гена TCF4 [9]. По нашим данным и данным зарубежных авторов, экспансия тринуклеотидных повторов CTG18.1 в гене TCF4 является наиболее частым вариантом, ассоциированным с ДФ, среди популяций европейского происхождения. Экспансия по крайней мере одного аллеля тринуклеотидного повтора CTG18.1 была обнаружена примерно у двух третей пациентов с ДФ в когортах европейского происхождения. Позднее Афшари и др. [10] предприняли попытку с помощью GWAS в большой когорте найти другие варианты, связанные с ДФ. Они подтвердили ассоциацию установленного локуса TCF4, а также идентифицировали три новых локуса в генах KANK4, LAMC1 и вблизи гена ATP1B1, однако их независимость от экспансии тринуклеотидного повтора не была проверена [10]. Обсуждается также роль мутаций в генах ZEB1 [11], SLC4A11 [12], AGBL1 [13] и LOXHD1 [14] как каузальных и независимых от экспансии повторов для реализации ДФ. Вопрос о том, является ДФ набором фенокопий или полигенным заболеванием, до сих пор остается открытым. Но известное бессимптомное носительство экспансии повторов [9, 15] и спорность четкой моногенной связи ДФ с некоторыми другими генами позволяют предположить, что ДФ с поздним началом — это набор полигенных фенокопий, что делает ее схожей с другими болезнями экспансий с поздней манифестацией.

В связи с этим представляется интересным вопрос — можно ли и с какой точностью разделить пациентов на группы внутри локусов по результатам GWAS. И если речь идет о ДФ позднего начала, можно ли по данным микрочипового анализа определить пациентов с фенокопией с экспансией и без нее? Сравнимы ли для этих целей использование гаплотипного подхода/комбинаций вариантов и носительство минорных аллелей отдельных генетических вариантов? Целью данной работы было разработать и проверить подход разделения пациентов на группы на основе результатов чип-генотипирования и полногеномного ассоциативного исследования (GWAS).

ПАЦИЕНТЫ И МЕТОДЫ

Изложенный в статье анализ был выполнен на образцах следующих двух выборок базы данных dbGaP: мультицентрового исследования по изучению дистрофии Фукса (FECD Genetics Multi-center Study) [16] и исследования возрастных заболеваний глаз (Age-Related Eye Disease Study AREDS, Refractive Error Substudy) [17–18] с использованием чипов Illumina HumanOmni2.5-4v1.

Клинические проявления заболевания были классифицированы с использованием модифицированной оценочной шкалы Krachmer J., основанной на результатах биомикроскопии щелевой лампой [19].

Контроль качества данных осуществляли на уровне образцов и на уровне вариантов. Для подготовки данных генотипирования использовали программное обеспечение PLINK 1.9 [20], GRAF 2.4 [21–22], а также собственный программный код, написанный на языке R версии 4.1.0.

Для контроля качества вариантов на первом шаге были удалены генотипы с низкими показателями качества прочтения GenCall (GC) путем удаления вариантов с GC < 0,3. В анализ были  включены маркеры, удовлетворяющие следующим критериям: доля пропущенных генотипов < 10%, частота минорного аллеля > 1%, число ошибок Менделевского наследования < 2, p-value точного теста на равновесие Харди–Вайнберга > 1 × 10^–10 для контрольных образцов и p-value > 1 × 10^–15 — для пациентов с ДФ. Маркеры, являющиеся дупликатными, т. е. маркеры с различными ID, но идентичными геномными позициями и кодировкой аллелей, были детектированы и проанализированы отдельно. Пары дупликатных маркеров с дискордантными генотипами более чем в 10 образцах были полностью удалены из рассмотрения. Из каждой пары дубликатов с 10 отличиями и менее было исключено по одному маркеру с наименьшим показателем доли пропущенных генотипов (genotype missing rate). После применения всех фильтров в анализе осталось 1580746 SNP-маркеров.

На основе метаданных были сформулированы следующие критерии включения в группу пациентов с ДФ: возраст 47 лет и старше; кератопластика не менее, чем на одном глазу, или наличие 2-й или выше степени развития заболевания по модифицированной шкале Krachmer J., манифестировавшего не менее чем на одном глазу.

Критерии включения в группу контроля: возраст 60 лет и старше; нормальное развитие роговицы глаза (без зафиксированных в метаданных аномалий в состояниях эпителия, эндотелия и стромы), за исключением травм роговицы.

Критерии исключения: образцы с ошибками менделевского наследования; образцы с несовпадениями между заявленным полом и полом, установленным на основе генетических данных (по гетерозиготности Х-хромосомы с одновременным  обнаружением вариантов Y-хромосомы по данным PLINK); образцы с отсутствием более 5% генотипов; родственные образцы до второй линии родства включительно (по данным GRAF-rel).

С целью получения генетически однородной выборки в популяционном приближении с использованием GRAF-pop была проведена оценка популяции образцов. Были отфильтрованы образцы, являющиеся выбросами в пространстве генетических координат. Участники исследований были рекрутированы на базе нескольких научных центров, использующих разные критерии включения участников в выборку пациентов с ДФ и в контрольную выборку. Чтобы унифицировать данные критерии, мы дополнительно провели фильтрацию образцов по фенотипическим признакам, руководствуясь единой для всех образцов логикой формирования групп сравнения.

Далее разбиение пациентов на подгруппы с предполагаемым наличием или отсутствием экспансии повторов осуществляли в несколько этапов:

1-й этап: ранжирование и отсечение наиболее значимых вариантов локуса;

2-й этап: кластеризация участников исследования по гаплотипам/комбинациям отобранных вариантов;

3-й этап: проверка соответствия полученного результата доле вероятного носительства экспансии повторов по данным, известным из публикаций. Аллель повторов считали экспансированной, если число повторов составляло более 40 копий, и не экспансированной, если число повторов было менее 40.

Для реализации первого этапа был проведен полногеномный поиск ассоциаций с использованием модели логистической регрессии c добавлением пола и первых шести главных компонент в качестве ковариат. Значения p-value были скорректированы на множественное тестирование методом Бенджамини–Хохберга по FDR. Набор вариантов из 18 хромосомы (где расположен локус с повторами) сначала был отфильтрован по p-value < 1 × 10^–15. Для сравнения с гаплотипным подходом из полученного набора вариантов были рассмотрены три SNP с самыми низкими значениями p-value в качестве возможных маркеров увеличенного числа повторов. Затем все варианты локуса были проведены через процедуру прунинга — сокращения числа анализируемых вариантов путем удаления высокоскоррелированных вариантов с попарной корреляцией: r² > 0,6. Матрица генотипов перед дальнейшим анализом была закодирована в соответствии с доминантной моделью наследования.

На втором этапе использовали предположение, что пациенты с экспансией повторов в гене TCF4 будут обладать определенным набором SNP. Мы ожидали, что образцы ДФ из анализируемого сета сгруппируются в локусе гена TCF4 на основании гаплотипов и комбинаций одиночных вариантов. При этом люди с фенокопией, обусловленной экспансией, на основании сходства комбинаций минорных вариантов попадут в общие кластеры. В эти же кластеры попадут люди из контрольной выборки, которые являются бессимптомными носителями повторов (2–10%), а также небольшая часть контрольной выборки, несущая минорные гаплотипы, не содержащие повторов. При этом комбинации мажорных вариантов и гаплотипы с преимущественно мажорными вариантами сформируют кластеры преимущественно из представителей контрольной выборки. Но сюда попадет и часть пациентов с ДФ с фенокопией без экспансии повторов, и часть пациентов с экспансией, более не сцепленной с минорными гаплотипами (у 7% пациентов с ДФ частично гаплотипы расцепляются с повторами, что ранее было показано для варианта rs613872 [23]). Поэтому доля пациентов с ДФ и участников без ДФ может быть использована как суррогатный маркер носительства повторов в конкретных кластерах.

Для разделения пациентов на группы к выборке был применен метод иерархической агломерационной кластеризации, реализованный в функции hclust пакета stats. Данный алгоритм выстраивает данные в древовидную структуру, последовательно объединяя пары кластеров с наименьшим расстоянием до тех пор, пока все элементы не будут объединены в один кластер.  Алгоритм принимает на вход матрицу попарных расстояний между точками (образцами), на начальном этапе каждая точка представляет отдельный кластер. Из-за неидентичности гаплотипов мы ожидаем больше двух кластеров, а оптимальное их количество определяется метрикой  Silhouette.

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

Для подобранных трех SNP носительство минорного аллеля по группе контролей считали маркером носительства увеличенного числа повторов.

Для подтверждения правильности полученного разбиения было вычислено отношение шансов обнаружить экспансию в группе пациентов с ДФ по сравнению с группой контроля. Для отношения шансов был дополнительно рассчитан 95%-й доверительный интервал.

Для сравнения с экспериментальными данными нами был проведен подбор исследований по следующим критериям.

В исследовании методом фрагментного анализа или ПЦР с праймированием триплетов определяли число повторов в интроне гена TCF4.

Исследовали выборку людей европейского происхождения.

Число исследованных людей в группе больных ДФ или в группе контроля — не менее 50 человек.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

После реализации процедуры контроля качества была сформирована комбинированная рабочая выборка из 3660 образцов европейского происхождения (табл. 1).

После применения всех фильтров в анализе осталось 1 580 746 SNP-маркеров.

Так как GWAS был проведен на тех же наборах данных, что были использованы в исследовании Afshari et al, 2017 [10], его результаты (рис. 1) сопоставимы с результатами, описанными в статье. Коэффициент генетической инфляции составил 1.05, что свидетельствует о незначительном расслоении населения.

Для дальнейшего анализа рассматривали только локус 18-й хромосомы. В результате фильтрации вариантов по значениям p-value осталось 134 варианта, из которых три имеют наиболее низкие значения p-value: rs784257, rs72932578 и rs618869 (и частоты минорных аллелей С, T и С в европейской популяции по данным проекта gnomAD v3.1.2, равны 0,17932, 0,05649 и 0,13451 соответственно) и были дополнительно проанализированы на предмет разделения пациентов на группы.

Размер гаплотипного блока составил 50 вариантов, оставшихся после процедуры прунинга. В результате кластеризации образцы были распределены по 10 подгруппам (рис. 2). Из них в трех подгруппах преобладали образцы пациентов с ДФ, в семи подгруппах преобладали образцы контрольной группы.

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

Наш анализ показал, что доля людей из кластеров с предположительным носительством экспансированных повторов в группе образцов с ДФ значительно выше, чем в контрольной группе (табл. 2). При этом рассчитанная доля людей с предположительным носительством экспансии повторов сильно варьирует в зависимости от выбора метода — предсказание экспансии по гаплотипам/ комбинациям вариантов  или по генотипам отдельных вариантов с низкими показателями p-value).

Для проверки полученных результатов нами был проведен подбор исследований с результатами экспериментального определения экспансии повторов. Традиционно число повторов определяют каноническим фрагментным анализом или ПЦР с праймированием триплетов с последующим фрагментным анализом. Было подобрано пять публикаций с подходящими выборками (табл. 3).

Для сравнения предсказанной частоты носительства экспансии с частотой носительства экспансии в исследованиях было проведено объединение выборок из публикаций. Сравнительный анализ  показал, что маркеры в разной степени воспроизводят частоту носительства экспансии повторов в группах сравнения (рис. 3).

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

Интересно отметить, что рассмотренные нами одиночные варианты дают сильно дискордантные между собой результаты (рис. 4), т. е. во многом разные люди являются носителями минорных аллелей в данных вариантах, что делает неустойчивой используемую метрику. rs784257 сильнее всего отличается по носительству аллелей от гаплотипного подхода и одновременно является наиболее значимым по результатам GWAS. При этом он показывает максимальное несоответствие по долям потенциального носительства экспансий в группе контроля и не лучшее соответствие с группой ДФ. Скорее всего данный аллель менее сцеплен с носительством повторов, чем два других.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Молекулярно-генетическая классификация пациентов с полигенными заболеваниями является полезным инструментом для изучения генетической природы заболеваний. При этом в рамках одного фенотипа могут существовать пациенты с группами каузальных вариантов, сцепленных с разными гаплотипами. И хотя ген явно связан с заболеванием, p-value вариантов будет выше из-за обилия групп сцепленных вариантов, т. е. изза особенностей генетической структуры заболевания варианты, значимые для каждой из групп в отдельности, не преодолевают общепринятый критерий значимости (p-value < 5 × 10^–8 ). К тому же для многих локусов, приближающихся по значимости к общепринятому порогу, характерна заметная разреженность по вариантам в области низких значений p-value, при которой только единичные варианты показывают сильные ассоциации с заболеванием. Таким образом, без дополнительного анализа не представляется возможным выбрать между отнесением геномных вариантов к популяционным выбросам (значимость которых обусловлена случайным популяционным смещением частот) или к потенциально значимой группе вариантов. Поэтому методы структуризации генетических данных, включающие анализ взаимодействий как вариантов в блоках гаплотипов, так и гаплотипных блоков между собой, являются перспективным инструментом для уточнения генетики заболевания.

GWAS позволяет получить больше информации о генетической структуре заболевания, чем исключение вариантов по p-value и формирование локусов с привязкой к ближайшему гену и по сути переходом с уровня «вариант» на уровень «ген». При этом остаются вопросы о неравноценности вклада разных локусов в генетическую природу заболевания у конкретных групп людей с одним фенотипом. Это связано в том числе с отсутствием развитых подходов к формированию комбинаций вариантов, т. е. к работе на промежуточном уровне между уровнями «ген» и «вариант». Из-за неполной сцепленности вариантов между собой целесообразно рассматривать совокупности гаплотипов/комбинаций вариантов, определяющие дифференцированный риск заболевания, а не конкретные рисковые или протективные гаплотипы в отдельности, т. е. вариант с самым высоким популяционным атрибутивным риском (сочетание частоты аллелей и относительного риска), вероятно, и будет самым значимым в локусе.

До сих пор анализ групп гаплотипов, сцепленных с каузальными вариантами, является сложной задачей, однако все чаще превосходит методологию GWAS, даже несмотря на отсутствие высокопроизводительного единого подхода. В 2005 г. в результате GWAS было показано, что ген CFH ассоциирован с возрастной макулодистрофией сетчатки [27]. А позже появилось сообщение, что эта ассоциация не ограничивается одиночными вариантами, и видна в том числе в группах пациентов со структурными перестройками, в частности с частичными делециями генов CFHR1-5 [27]. Более того, было выявлено, что большая часть вариабельности, которую приписывали одиночным вариантам, на самом деле являлась маркером гаплотипов с крупными структурными перестройками в этом регионе. И именно гаплотипные варианты структуры локуса, в том числе с разной популяционной представленностью, уже имеют намного более сильную связь с риском дегенерации сетчатки, чем большинство одиночных вариантов в этом локусе [28].

В данной работе мы реализовали кластеризацию образцов по вариантам региона с экспансией, опираясь на данные об ассоциации с повторами отдельных вариантов [14, 23, 29], в частности, аллеля G варианта rs613872, и гаплотипных блоков [29]. В результате кластеризации образцы из группы пациентов с ДФ и контрольной группы оказались неравномерно распределены по кластерам, что является косвенным свидетельством кластеризации по гаплотипам, сцепленным с экспансией. Все кластеры за исключением одного (кластера № 3), обладали однозначным статусом. Неопределенность при выявлении статуса была обусловлена равным соотношением в кластере долей пациентов и лиц контроля. В дальнейшем предстоит принять решение, что делать с такими кластерами — перекластеризовать людей из них отдельно или оставлять с неопределенным статусом. Необходимо также подобрать другую метрику качества кластеризации, а для этого нужны дополнительные данные, содержащие одновременно и результаты генотипирования образцов, и информацию о длине повторов. Несмотря на данные ограничения, результаты, полученные с использованием гаплотипного подхода, превзошли результаты, показанные одиночными вариантами, однако оказались недостаточно точными, чтобы считать использованную нами реализацию оптимальной.

Данная работа позволила реализовать две цели.

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

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

ВЫВОДЫ

В данном исследовании была продемонстрирована возможность применения подхода, основанного на гаплотипах, к молекулярно-генетической классификации пациентов по причинам генетического заболевания, в данном случае, по наличию экспансии повторов. Результаты работы позволили нам сформулировать следующие выводы: 1) гаплотипный подход лучше, чем отдельные варианты, подходит для выявления связи локусов с отдельными группами пациентов; 2) для получения более точной картины следует пересмотреть подход к определению гаплотипного состава и моделированию матрицы данных для кластеризации. В частности планируется проанализировать ряд методов по вычислению генетического сходства (генетических расстояний) между образцами, а также применить более специфичные методы для начального отбора вариантов; 3) полученные результаты указывают на то, что кластеризация разделяет пациентов с ДФ и группу контроля с учетом групп гаплотипов/комбинаций вариантов, ассоциированных с экспансией повторов. Для дальнейших этапов тестирования подхода необходима дополнительная доказательная база, требующая подключения большего объема валидационных данных.