Исследование проводилось в соответствии с Хельсинкской декларацией, а протокол был одобрен Комитетом по этике больницы грудных заболеваний провинции Аньхой (K2025-007) 22 апреля 2025 года. От всех участников исследования было получено информированное согласие.
Извлечение и нормализация данных
Транскриптомные профили и соответствующие клинические наборы данных для LUAD были получены из когорт TCGA и GEO. Набор данных TCGA-LUAD был назначен в качестве обучающей выборки, а GSE72094, GSE31210 и GSE26939 служили когортами для внешней валидации (Таблица 1). Кроме того, 900 MCRGs были собраны из предыдущего исследования12(Дополнительная таблица 1). Данные транскриптома были аннотированы с использованием GENCODE v36 или соответствующих файлов аннотации платформы GPL. Идентификаторы зондов были преобразованы в символы генов, дублирующиеся гены были объединены с помощью функции avereps, и для создания матриц экспрессии на уровне генов были оставлены только белок-кодирующие гены. Для обучающей выборки TCGA-LUAD были отфильтрованы гены с количеством фрагментов на килобазу экзонной модели на миллион картированных фрагментов (FPKM) < 1 в более чем 50% образцов, а оставшиеся значения экспрессии были подвергнуты log2-преобразованию (log2[FPKM+1]). Для валидационных когорт GEO были загружены необработанные данные экспрессии, идентификаторы зондов были сопоставлены с символами генов с использованием соответствующих файлов аннотации платформ, а несколько зондов, соответствующих одному и тому же гену, были объединены путем усреднения их значений экспрессии. Эти наборы данных были подвергнуты log2-преобразованию при необходимости. Коррекция пакетного эффекта между платформами TCGA и GEO не применялась, так как была принята стратегия стандартизации для каждой когорты отдельно для обеспечения относимой сопоставимости. В частности, как для обучающей, так и для валидационных когорт значения экспрессии генов были центрированы и масштабированы (z-преобразование) с использованием среднего значения и стандартного отклонения каждого набора данных по отдельности. Затем те же коэффициенты регрессии Кокса, полученные из обучающей выборки, использовались для расчета показателей риска для всех когорт. Чтобы сохранить клиническую применимость и избежать переобучения под любой валидационный набор, медианный показатель риска обучающей когорты использовался в качестве фиксированного порога для стратификации пациентов на группы высокого и низкого риска во всех внешних валидационных когортах. Была извлечена имеющаяся клиническая информация, включая возраст, пол, патологическую стадию, стадию Tumor-Node-Metastasis (TNM), гистологический тип, время выживаемости, статус выживаемости и тип ткани. Конечной точкой была общая выживаемость (OS). Образцы с неполной информацией о выживаемости или временем выживаемости < 30 дней были исключены. Время выживаемости было переведено в годы, а статус выживаемости был закодирован как 0 для живых и 1 для умерших.
Идентификация и функциональный анализ генов-кандидатов
Пакет Limma использовался для идентификации дифференциально экспрессируемых генов (ДЭГ) между образцами опухоли LUAD и нормальными образцами в обучающей выборке13. ДЭГ определялись по следующим критериям: |log2FC| > 0.5 и скорректированное значение p < 0.05. Впоследствии для разделения ДЭГ на отдельные кластеры экспрессии применялся алгоритм нечеткой кластеризации mfuzz из R-пакета ClusterGVis. Анализ биологических процессов Gene Ontology (GO-BP) проводился для пяти наиболее репрезентативных генов в каждом кластере на основе их показателей принадлежности. Набор общих генов был получен путем нахождения пересечения ДЭГ с MCRG. Для оценки биологической значимости перекрывающихся генов проводился функциональный анализ обогащения с использованием Gene Ontology/Киотской энциклопедии генов и геномов (GO/KEGG). Сети белок-белковых взаимодействий (PPI) были построены на основе базы данных STRING14. Для повышения надежности сети были оставлены только взаимодействия с показателем достоверности > 0.7.
Скрининг прогностических генов
Для проведения одномерного регрессионного анализа Кокса с целью выявления вероятных генов, связанных с общей выживаемостью при LUAD, использовался пакет Survival15. Гены со значением p < 0,05 рассматривались как потенциальные прогностические показатели. Обучающая когорта TCGA-LUAD включала 500 пациентов с полными данными о выживаемости, из которых у 216 (43,2%) в период последующего наблюдения наступила смерть. Соотношение количества генов-кандидатов (n = 108) к количеству событий (n = 216) составило примерно 1:2, что является приемлемым для регрессионного анализа Кокса. Впоследствии для дальнейшего отбора признаков применялся регрессионный анализ LASSO (Least Absolute Shrinkage and Selection Operator) и модель XGBoost (Extreme Gradient Boosting). Модели пропорциональных рисков Кокса строились с параметром family = "cox" с помощью функции cv.glmnet пакета glmnet. Оптимальный параметр регуляризации определялся с использованием 10-кратной перекрестной проверки; в качестве оптимального значения λ было выбрано значение λ.min, соответствующее минимальной ошибке перекрестной проверки. В качестве потенциальных признаков были извлечены гены с ненулевыми коэффициентами регрессии. Для модели XGBoost время выживания и статус выживаемости были объединены в целевую переменную, при этом события смерти обозначались положительными значениями, а цензурированные случаи — отрицательными. Параметры были установлены следующим образом: objective = "survival: cox" и eval_metric = "cox-nloglik", при 100 итерациях и скорости обучения 0,1. После обучения модели показатели значимости генов рассчитывались с использованием значений прироста признака (feature gain). После сортировки показателей значимости в порядке убывания были отобраны 20 наиболее значимых генов для снижения размерности признаков и сложности модели. Гены, совпадающие в результатах LASSO и XGBoost, были определены как гены-кандидаты для прогноза.
Построение и оценка прогностической модели
Прогностическая модель была разработана с использованием многофакторного регрессионного анализа Кокса для выявленных генов-кандидатов. Показатели риска рассчитывались индивидуально следующим образом:
.
где Coefi обозначает коэффициент для гена i, а Expi указывает соответствующее значение экспрессии гена. Впоследствии пациенты были разделены на две группы: высокого и низкого риска, используя медианное значение показателя риска в качестве порогового значения. Затем были построены зависящие от времени ROC-кривые (receiver operating characteristic). Для оценки вероятности переобучения была проведена внутренняя бутстреп-валидация с 1 000 итерациями ресемплирования для расчета скорректированного на смещение C-индекса и зависящих от времени значений AUC с 95% доверительными интервалами. Для оценки соответствия между прогнозируемой и наблюдаемой вероятностью выживаемости через 2, 3 и 5 лет были построены калибровочные кривые. Кроме того, с использованием пакета ggDCA в среде R был проведен анализ кривых принятия решений (DCA) для оценки клинической чистой выгоды модели во временных точках 2, 3 и 5 лет, что позволило количественно определить потенциальную ценность показателя риска при принятии клинических решений при различных пороговых вероятностях. Различия в выживаемости между группами с разным уровнем риска и в других клинических категориях сравнивали с помощью кривых выживаемости Каплана-Мейера (KM) и логрангового критерия. Кроме того, для выяснения вклада отдельных генов в эффективность модели использовался анализ Shapley Additive exPlanations (SHAP) для апостериорной интерпретации результатов.
Разработка и внешняя валидация номограммы
Для оценки клинической применимости модели с помощью тестов Уилкоксона или Краскела-Уоллиса изучались взаимосвязи между рассчитанными показателями риска и различными клиническими характеристиками (включая пол, возраст и стадию по системе TNM). Чтобы определить, является ли показатель риска независимым прогностическим фактором, клинические переменные вместе с показателем риска были включены в многофакторное регрессионное моделирование Кокса. Затем с помощью пакета regplot для языка R была построена прогностическая номограмма, объединяющая независимые клинические факторы риска (например, стадию) и генетический показатель риска, для индивидуализации прогнозов вероятности выживаемости. Калибровочные кривые использовались для оценки соответствия вероятности выживаемости, предсказанной номограммой, фактическим исходам выживаемости. Наконец, итоговая прогностическая способность и обобщаемость интегрированной системы номограмм были строго подтверждены с помощью зависимых от времени ROC-кривых и комплексного анализа клинических подгрупп по методу KM во всех когортах.
Анализ иммунной инфильтрации и иммунных подтипов
Для оценки относительной доли 22 типов иммунных клеток и анализа их инфильтрации у пациентов с LUAD использовали метод CIBERSORT с матрицей сигнатур лейкоцитов (LM22). Взаимосвязь между прогностическими уровнями экспрессии генов и иммунологической инфильтрацией оценивали с помощью корреляционного анализа Спирмена. Показатели иммунного компонента, стромального компонента, чистоты опухоли и общий балл ESTIMATE были получены с помощью алгоритма ESTIMATE, а различия между группами риска оценивали с помощью критерия Вилкоксона. Пациенты с LUAD были распределены по шести иммунным подтипам с использованием пакета ImmuneSubtypeClassifier16. Также применяли критерий Вилкоксона для сравнения распределения иммунных подтипов в различных группах риска.
Анализ контрольных точек иммунитета, иммунофеноскора и цикла противоопухолевого иммунитета
В данном исследовании для оценки экспрессии 21 гена иммунных контрольных точек17 в группах с разным уровнем риска, с целью характеристики иммунного ландшафта LUAD, использовался критерий суммы рангов Уилкоксона. Корреляция по Спирмену позволила связать потенциальные прогностические гены с генами иммунных контрольных точек. Чтобы оценить различия в ответе на ингибиторы иммунных контрольных точек (ICI) у пациентов с LUAD разного уровня риска, данные иммунофеноскора (IPS) для терапии анти-PD-1 и анти-CTLA-4 были получены из Атласа иммуномов рака (The Cancer Immunome Atlas, TCIA)18, а база данных Tracking Tumor Immunophenotype (TIP)19 была использована для оценки активности цикла «рак — иммунитет» путем сравнения соответствующих баллов между группами риска.
Анализ соматических мутаций и чувствительности к препаратам
С помощью инструмента анализа мутаций TCGA были получены профили соматических мутаций для случаев TCGA-LUAD с целью изучения вариаций в паттернах мутаций в разных группах риска. Пакет maftools использовался для обработки и визуализации данных о мутациях. Для каждого образца определялся уровень опухолевой нагрузки мутациями (TMB), который затем сравнивался между двумя категориями риска. Фармакогеномный анализ чувствительности проводился с использованием пакета pRRophetic на основе базы данных Genomics of Drug Sensitivity in Cancer (GDSC)20. Для каждого пациента с LUAD были предсказаны значения полумаксимальной ингибирующей концентрации (IC50) противоопухолевых препаратов, а различия между группами риска оценивались с помощью рангового суммарного критерия Вилкоксона.
Оценка уровней экспрессии прогностических генов
Каждый набор данных использовался для оценки уровней экспрессии выбранных генов-кандидатов, связанных с исходом заболевания. Для установления связи между экспрессией генов и прогнозом пациентов оптимальные пороговые значения были определены с помощью функции surv_cutpoint пакета survminer для языка R. На основании этих порогов случаи LUAD были разделены на подгруппы с высоким и низким уровнем экспрессии для последующего анализа выживаемости.
Кроме того, в больнице грудной клетки провинции Аньхой были получены образцы пяти пар соответствующих опухолей LUAD и прилегающих нормальных тканей, после чего была проведена валидация методом qPCR. Каждый участник предоставил письменное информированное согласие. Для валидации методом qPCR были выбраны шесть потенциальных прогностических генов (PDGFB, LDHA, ZEB2, FKBP4, DMD и S100B). РНК экстрагировали из гомогенизированных образцов тканей с использованием реагента для экстракции РНК с последующей экстракцией хлороформом и осаждением изопропанолом. С помощью спектрофотометра измеряли концентрацию и чистоту РНК. Валидация шести потенциальных прогностических генов методом qPCR проводилась с использованием мастер-микса для ПЦР на основе SYBR Green на системе ПЦР в реальном времени: начальная денатурация при 95 °C в течение 30 s, затем 40 циклов: 95 °C в течение 20 s, 55 °C в течение 20 s и 72 °C в течение 20 s. Относительную экспрессию рассчитывали и стандартизировали по глицеральдегид-3-фосфатдегидрогеназе (GAPDH) с помощью метода 2-ΔΔCt. Подробная информация о всех реагентах и приборах представлена в Таблице материалов.
Статистический анализ
Статистический анализ проводили с использованием программного обеспечения для статистических вычислений и построения графиков. Визуализацию сети белок-белковых взаимодействий осуществляли с помощью программного обеспечения для сетевого анализа. После проверки распределения на нормальность для непрерывных переменных с нормальным распределением использовали t-критерий Стьюдента, а для переменных с ненормальным распределением — U-критерий Манна-Уитни.