По данным организации мониторинга онкологических заболеваний GLOBOCAN, в 2020 г. зарегистрировано около 2,3 млн новых случаев и 684 996 смертей от рака молочной железы (РМЖ), который занимает первое место по частоте среди онкологических заболеваний во всем мире [1] и представляет собой высокогетерогенное заболевание с различными молекулярными и клиническими характеристиками [2].

В настоящее время в широкой клинической практике подтип РМЖ определяют в опухолевой ткани иммуногистохимическим (ИГХ) методом [3] в том числе по экспрессии в опухоли белков-рецепторов эстрогена, прогестерона, HER2, а также скорости деления опухолевых клеток. Разработка методов анализа экспрессии генов с помощью ДНК-микрочипов сыграла важную роль в определении молекулярных подтипов РМЖ. С помощью использования классификатора на основе экспрессии 50 генов PAM50 удается четко выделить люминальный A (LumA), люминальный B (LumB), HER2-обогащенный (HER2+) молекулярные подтипы и базальноподобный или тройной негативный РМЖ (ТН РМЖ) [4]. ТН РМЖ включает в себя 15–20% всех случаев РМЖ и характеризуется агрессивным течением, высоким уровнем метастазирования, частым возникновением рецидивов и низкой выживаемостью по сравнению с другими подтипами РМЖ [5].  Мультигенные микрочиповые диагностические системы позволяют получить важную прогностическую информацию для онкологических больных, особенно в случае двусмысленного прогноза по результатам клинических характеристик и иммуногистохимических маркеров. К таким системам относят Mammaprint/Blueprint и Prosigna/ PAM50, которые в дополнении к их прогностической и предикторной ценности предоставляют возможность разбиения на молекулярные подтипы [6]. Данные системы можно использовать для определения оценки высокого и низкого рисков рецидива для пациенток с РМЖ, однако для ТН РМЖ и HER2+-молекулярного подтипа такая возможность пока отсутствует ввиду недостатка клинических исследований.

Эпигенетические изменения модулируют использование генома с помощью модификаций гистонов, вариантного состава гистонов, ремоделинга хроматина, метилирования ДНК, расположения нуклеосом и некодирующих РНК (экспрессия специфичных микроРНК). Для проявления своего эффекта вышеперечисленные эпигенетические изменения действуют совместно. Метилирование ДНК является одним из наиболее известных факторов регуляции экспрессии генов. Оно происходит за счет ковалентной модификации цитозина путем присоединения метильной группы в контексте CpG-динуклеотида к 5'-углероду пиримидинового кольца [7]. CpG-динуклеотиды обычно концентрируются в CG-богатых участках ДНК под названием CpG-островки, значительная часть которых сосредоточена в промоторных участках гена и в областях с длинными повторами, например в ретротранспозонных элементах или центромерных повторах. Реакцию метилирования цитозинов опосредует класс ферментов под названием метилтрансферазы (DNMT) [7]. Всего у млекопитающих идентифицировано пять членов семейства DNMT: DNMT1, DNMT2, DNMT3a, DNMT3b и DNMT3L. DNMT3a и DNMT3b являются метилтрансферазами de novo, которые взаимодействуют с неметилированными CpGдинуклеотидами. DNMT1 выполняет функцию поддержания метилирования при репликации в S-фазе. Было показано, что DNMT3L стимулирует метилирование de novo, которое происходит с помощью DNMT3a и опосредует репрессию транскрипции с помощью гистондеацетилазы 1 (HDAC1) [7]. Аберрантное метилирование ДНК связано с широким спектром заболеваний и наиболее выражено в злокачественных опухолях [8]. Исследования последних лет продемонстрировали, что каждая эпителиальная опухоль содержит приблизительно 10–15 генов, инактивированных структурными изменениями генома, и несколько сотен, инактивированных гиперметилированием ДНК, что демонстрирует значение этой модификации в развитии опухоли. Другой особенностью опухолевых геномов является их тотальное гипометилирование. Это полногеномное гипометилирование, происходящее главным образом из-за потери метилирования повторяющихся элементов, ведет к геномной нестабильности и хромосомным перестройкам [8]. Важную роль в патогенезе РМЖ играет повышенное промоторное метилирование в генах-супрессорах опухолевого роста, которые подавляют различные механизмы опухолевого прогрессирования, приводящее к их эпигенетическому «умолканию» и обратимой инактивации [8]. Идентификация опухолеспецифичных паттернов аберрантного метилирования ДНК может быть полезным для ранней диагностики рака, дифференциальной диагностики злокачественных новообразований, в качестве прогностических и предикторных маркеров [9]. Изучение специфичных паттернов метилирования ДНК, выявляемых с использованием широкогеномных методов анализа, вносит важный вклад в понимание молекулярного патогенеза рака молочной железы [10]. Как отмечалось выше, каждый тип рака подразделяется на подтипы, и существуют геномные паттерны, в том числе эпигенетические, характерные для этих подтипов. Таким образом, необходимо проведение специфического широкогеномного профилирования метилирования ДНК при онкологических заболеваниях, наряду с традиционными исследованиями точечных событий гиперметилирования промоторов отдельных генов [11].

Прогнозирование подразумевает под собой предсказание вероятного течения и исхода онкологического заболевания. Анализ выживаемости основан на математическом подходе к прогнозированию онкологического заболевания и позволяет предсказать вероятность жизни после определенного времени. Маркеры метилирования ДНК благодаря их биологической важности и стабильности являются эффективным прогностическим фактором [12]. В одной из работ на основе данных результатов широкогеномного анализа метилирования ДНК в образцах РМЖ в составе базы The Cancer Genome Atlas Breast Cancer (TCGA-BRCA) была построена модель из семи CpG-динуклеотидов, в которой хорошо различают опухоли молочной железы всех подтипов и нормальные ткани, а также идентифицировано шесть сайтов метилирования, которые высоко коррелировали с общей выживаемостью (ОВ) [13]. В ходе анализа данных метилирования из открытых источников с помощью LASSO-регрессии и бустинга были обнаружены соответственно 29 и 11 CpG-динуклеотидов, связанных с ОВ [14]. Исследование данных из открытого источника TCGA-BRCA тоже позволило выявить три гена (TDRD10, PRAC2 и TMEM132C), состояние метилирования которых имеет прогностическую ценность, но преимущественно для эстроген-положительных опухолей молочной железы [15]. Для ТН РМЖ была разработана прогностическая модель на основании данных из источника TCGA-BRCA, состоящая из пяти генов (TGFBR2, EIF4EBP1, FOSB, BCL2A1, ADRB2), которая одинаково хорошо прогнозирует ОВ и безрецидивную выживаемость (БРВ) [16].

Необходимость в проведении исследования обусловлена отсутствием таких сигнатур для HER2-обогащенного подтипа и довольно ограниченным количеством для остальных молекулярных подтипов РМЖ. Диагностический потенциал существующих моделей для прогнозирования выживаемости тоже остается неясным, поэтому мы использовали модифицированный алгоритм для поиска CpG-динуклеотидов, ассоциированных со всеми доступными клиническими конечными точками из базы данных TCGA-BRCA.

Целью данного исследования было получить различные сигнатуры на основании открытых данных ДНК-метилирования в РМЖ The Cancer Genome Atlas Breast Cancer для прогнозирования различных типов выживаемости пациенток (общей выживаемости, безрецидивной выживаемости и выживаемости без прогрессирования) в молекулярных подтипах РМЖ и проверить зависимость клинико-патологических характеристик от полученных сигнатур.

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

Публично доступные клинические параметры, а также данные широкогеномного профилирования метилирования ДНК, полученные с помощью гибридизационных чипов HumanMethylation450 (HM450) (Illumina Inc.; США) проекта The Cancer Genome Atlas Breast Cancer (TCGA-BRCA) (https:// portal.gdc.cancer.gov/projects/TCGA-BRCA), получали и обрабатывали с помощью пакета программ TCGAbiolinks [17]. Критерии включения пациентов для дальнейшего отбора кандидатных CpG-пар были следующими: наличие соответствующего молекулярного подтипа РМЖ, наличие доступной клинико-патологической информации, наличие данных о профилировании метилирования ДНК с чипов Illumina HumanMethylation450. Критерии исключения: отсутствие данных о показателях времени для клинических конечных точек, возрасте пациента, классификации TNM и стадии. В дальнейшем из матрицы данных профилирования были исключены результаты, полученные с FFPE (formalin fixed paraffin-embedded) блоков пациенток, а также кроссгибридизационные зонды.

Отбор CpG-пар, ассоциированных с ОВ, БРВ или выживаемостью без признаков прогрессирования заболевания (ВБП), производили с помощью одномерной регрессии Кокса (Univariate Cox regression) [18]. Из отобранных CpG-пар в дальнейший анализ были включены прошедшие поправку на множественное тестирование (скорректированный показатель p < 0,05, использовали тест Уальда) методом средней доли ложных отклонений гипотез (False Discovery Rate, FDR). В целях отбора наиболее стабильных CpG-пар использовали метод LASSO-регрессии Кокса (Cox LASSO-regression) [19] пакета программ SurvHiDim [20]. Многомерную регрессию Кокса [21] применяли для расчета сигнатур из CpG-пар и проверки независимости клинических параметров пациенток на эти сигнатуры. С помощью логистической регрессии определяли способность к классификации различных исходов. Для стратификации на высокий и низкий уровень риска использовали медиану. Методом cvROC (cross-validated receiver operative curve) [22] определяли качество построенных моделей и строили ROC-кривые. Наилучшие показатели чувствительности и специфичности определяли с помощью индекса Юдена. Кривые Каплан–Майера были построены с помощью пакета survminer [23]. Для сравнения двух кривых выживаемости использовали критерий Мантеля–Кокса. Для всех этапов выбора маркеров и вычисления сигнатур использовали 10-кратную рандомизированную кросс-валидацию. Все вышеуказанные вычисления выполняли с помощью языка статистического программирования R [24].

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

Исследуемый массив данных TCGA-BRCA состоял из профиля ДНК-метилирования, полученного с помощью чипов HM450 и клинико-патологических характеристик 735 образцов первичных опухолей РМЖ. После исключения образцов из парафиновых блоков суммарно осталось 555 образцов LumA+B подтипа (LumAB), 134 образца ТН-подтипа и 46 образцов HER2-обогащенного подтипа (табл. 1). Перед отбором признаков для дальнейшего анализа из матрицы метилирования исключили кроссгибридизационные зонды, вследствие чего число зондов снизилось с 485 577 до 456 344 соответственно.

Дальнейшим этапом анализа было применение одномерной регрессии Кокса для поиска сайтов метилирования, которые коррелировали с продолжительностью ОВ, БРВ и ВБП в разных молекулярных подтипах РМЖ. Всего, после первоначального отбора, с учетом p-value с поправкой на множественное тестирование, были выделены:

• в подтипах LumAB — 10 433 зондов, ТН — 3214 зондов, HER2-обогащенном — 6471 зонд, ассоциированный с ОВ;
• в подтипах LumAB — 4419 зондов, ТН — 168 зондов, HER2-обогащенном — 483 зонда, ассоциированных с БРВ;
• в подтипах LumAB — 2345 зондов, ТН — 43 зонда, HER2-обогащенном — 3216 зондов, ассоциированных с ВБП.

Для каждого из этих наборов была применена LASSOрегрессия Кокса, что позволило отобрать наиболее важные для анализа CpG-динуклеотиды. На каждом этапе кросс-валидации выявляли разное число таких CpG-пар, и были отобраны такие CpG-пары, которые встречались больше чем в 50% кросс-валидационных разбиений (табл. 2).

Для выбора сочетаний CpG-динуклеотидов, обеспечивающих значимую связь для различных типов выживаемости, были оценены все возможные комбинации (сигнатуры) таких CpG-динуклеотидов в различных молекулярных подтипах РМЖ. Для каждой клинической конечной точки и молекулярного подтипа РМЖ были получены результаты cvAUC (cross-validated area under curve, усредененная площадь под кривой на каждом этапе кросс-валидации), чувствительности, специфичности и точности для различных комбинаций. Первые 10 комбинаций с высокими показателями cvAUC были проверены на независимость от клинико-патологических характеристик, и диагностические характеристики этих комбинаций вместе с числом зондов и указанием принадлежности генов к зондам представлены в табл. 3.

Наибольшей установленной нами комбинацией является комбинация из двенадцати CpG-динуклеотидов для прогнозирования ОВ в LumAB подтипе, а меньшей — сочетание из двух CpG-динуклеотидов для прогнозирования БРВ в HER2-экспрессирующем подтипе. Для каждой сигнатуры были построены кривые cvROC (cross-validated receiver operative characteristics, на каждом этапе кросс-валидации строится ROC-кривая, в дальнейшем строится усредненная кривая) и кривые Каплан– Майера для отображения диагностического потенциала и оценки функции выживаемости. Комбинации для LumAB имели меньший показатель cvAUC (от 0,76 до 0,83), тогда как комбинации для ТН и HER2-экспрессирующего подтипа показывали высокий уровень cvAUC при меньших количествах сочетаний (от 0,83 до 1) (рис. 1).

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

Анализ кривых Каплана–Майера показал статистически значимое (p < 0,05) снижение ОВ, БРВ и ВБП в группе пациенток с высоким риском летального исхода, рецидива и прогрессирования болезни по сравнению с группой пациенток с низким риском при всех молекулярных подтипах РМЖ для каждой из выбранных комбинаций (рис. 2).

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

В настоящем исследовании с помощью методов анализа выживаемости, а также с использованием данных ДНК-метилирования была рассмотрена возможность идентификации сайтов дифференциального метилирования CpG-динуклеотидов для прогнозирования видов выживаемости в различных молекулярных подтипах РМЖ. Подход для расчета дифференциального метилирования с помощью одномерной регрессии Кокса широко используют в различных работах. Так, этот метод применяли для идентификации 249 810 и 249 811 зондов в данных ДНК-метилирования при раке яичников и при РМЖ соответственно [12] и для идентификации зондов в данных ДНК-метилирования при меланоме кожи [25]. 

Нами показано, что при использовании различных комбинаций (от 2 до 12 CpG-динуклеотидов) можно добиться приемлемого (cvAUC между 0,7 и 0,8), хорошего (от 0,8 до 0,9) и очень хорошего (от 0,9 до 1) качества классификации высокого и низкого риска летального исхода, рецидива и прогрессирования. В рамках данной работы были идентифицированы 47 зондов/генов (SLC30A7, EXTL2, C15orf41, MIA3, NIPAL3, HEY2, HK1, DIRC3, TMEM41A, SH3BP5L, RFX2, SLC25A39, BAT2, ZNF417, PSMA6, RG9MTD3, ZNF827, ABCC5, HLA-DRB5, HIST3H2A, RERE, SPAG5, cg13447284, cg03512997, LIN54, RASGRP2, LDLRAD3, ZNF643, PKNOX1, KCNMB2, ZFAND1, HDAC9, cg13745678, DPPA5, cg02927111, PKNOX1, SSU72, CADPS2, PEX5L, GSTM4, cg26290926, BIRC5, cg10660854, SLC43A1, BOD1, cg00297843, KCNN1), метилирование которых ассоциировано с ОВ, БРВ и ВБП, при этом для семи из них (cg13447284, cg03512997, cg13745678, cg02927111, cg26290926, cg10660854, cg00297843) не удалось установить принадлежность к генам, а пять из них (BIRC5, PKNOX1, SPAG5, HDAC9, PSMA6) ранее были описаны в научной литературе как молекулярные маркеры выживаемости больных с РМЖ на основе данных об количественной экспрессии генов [26–31].

Примечательно, что в рамках данного исследования мы обнаружили такое же число CpG-динуклеотидов (зондов HM450) для прогнозирования ОВ и БРВ (по пяти зондам), что и в другой работе [16], но зонды различались по принадлежности к генам. Согласно нашим данным, сигнатура для прогнозирования ВБП состоит из шести зондов; в работе [16] сигнатуру ВБП не вычисляли. Другие исследователи предлагают использовать для прогноза ОВ и БРВ больных с ЭР + РМЖ индивидуальные маркеры промоторного гиперметилирования семи генов (RASSF1, BRCA1, PITX2, RARB, PGR, CDH1 и PCDH10), а также рассматривают использование панели из трех генов (GSTP1, RASSF1 и RARB) для предсказания ОВ, основываясь на анализе литературных данных (по результатам систематического обзора публикаций) [26], в то время как в нашем исследовании применена стратегия формирования панелей из шести, девяти и 12 маркеров метилирования с учетом диагностического потенциала маркеров, установленного статистическим анализом массива экспериментальных данных.

Среди генов, вошедших в полученные нами комбинации, внимание привлекает ген BIRC5 (кодирует белок — бакуловирусный ингибитор мотива апоптозных повторов 5), который гиперэкспрессируется в большинстве опухолей, в том числе при РМЖ, связан с худшим прогнозом общей, безрецидивной и метастатической выживаемости. Показано, что использование химиотерапевтических препаратов из группы таксанов может увеличивать экспрессию данного гена [27]. Ген PKNOX1 (кодирует одноименный белок, расположен на коротком плече 21-й хромосомы, играет важную роль в эмбриональном развитии) является онкосупрессором, а его повышенная экспрессия ассоциирована с худшей выживаемостью [28]. Прогностическим фактором является также повышенная экспрессия гена SPAG5 (кодирует белок, связанный с аппаратом митотического веретена), ассоциированная с худшим прогнозом ОВ, БРВ и выживаемости без метастазирования только для эстроген-положительных (ЭР+) опухолей молочной железы [29], что также подтверждается нашим исследованием.

Результаты исследования, проведенного на образцах ЭР+ РМЖ, показывают, что повышенная экспрессия гена эпигенетического фермента HDAC9 (кодирует белок фермента гистондеацетилазы 9) в опухолях связана с худшим прогнозом БРВ [30]. В нашей работе показана ассоциация аномального метилирования этого гена с выживаемостью для больных с эстроген-отрицательными (ЭР–) опухолями, точнее с ТН РМЖ. Сниженная БРВ больных с ЭР+ РМЖ была показана при повышенной экспрессии гена PSMA6 (кодирует белок протеасомной субъединицы альфа-типа-6) [31], что также подтверждают результаты нашего исследования.  

ВЫВОДЫ

С помощью методов анализа выживаемости, сочетания различных сайтов метилирования и оценки диагностических показателей были обнаружены молекулярные эпигенетические сигнатуры для различных подтипов РМЖ. Данную методику можно рекомендовать для поиска сигнатур, характерных для РМЖ, а также для других опухолевых заболеваний. В дальнейшем обнаруженные нами эпигенетические сигнатуры можно использовать для разработки диагностических тестсистем на основе метилчувствительной количественной ПЦР. Такие тест-системы при прохождении клинических испытаний могут стать более дешевой и практичной альтернативой экспрессионным микрочипам без ущерба для диагностической эффективности.