Заявление институционального обзора
Данное исследование было проведено в соответствии с Хельсинкской декларацией. Протокол был одобрен Комитетом по этике больницы традиционной китайской медицины района Луху города Шэньчжэнь (номер одобрения 2024-LHQZYYYXLL-KY-039), и от всех участников было получено письменное информированное согласие перед включением в исследование. Подробная информация об инструментах и материалах, использованных в данном протоколе, представлена в Таблице материалов.
Источник и обработка данных
Наборы данных об экспрессии генов, связанных с ХОБЛ, были получены из базы данных Gene Expression Omnibus (GEO). Набор данных GSE54837 использовался в качестве набора данных транскриптома, а набор данных GSE112811 служил валидационным набором (Таблица 1). Список ac4C-RGs был собран на основе литературных данных18. Дифференциально экспрессируемые гены (ДЭГ) между группами с ХОБЛ и контрольной группой были определены с помощью пакета лимма (limma) для языка R. ДЭГ считались статистически значимыми при |log2FC| > 0 и p < 0,05. Для визуализации общего распределения изменений экспрессии генов были построены графики «вулкан» (volcano plots).
Построение WGCNA
Для выявления модулей, связанных с ХОБЛ, был проведен анализ WGCNA набора данных GSE54837 с использованием языка R. Перед построением сети были выявлены и удалены образцы-выбросы с помощью иерархического кластерного анализа с использованием функции hclust с методом среднего сцепления (average linkage) и метрикой евклидова расстояния. Для достижения индекса соответствия безмасштабной топологии R2 ≥ 0.85 при балансе между безмасштабной топологией и средней связностью была выбрана оптимальная мощность мягкого порога (β = 10). Была построена матрица смежности, которая затем была преобразована в матрицу топологического перекрытия (TOM). Генетические модули были определены с помощью алгоритма динамического разрезания дерева (deepSplit = 2, minClusterSize = 50). Модули с корреляцией собственных генов > 0.75 впоследствии были объединены с помощью функции mergeCloseModules. Затем собственные гены модулей были сопоставлены с клиническими признаками (наличие ХОБЛ, возраст, пол и статус курения) с использованием коэффициентов корреляции Пирсона для выявления модулей, ассоциированных с ХОБЛ, для последующего анализа.
Скрининг, анализ обогащения и анализ сети PPI для перекрывающихся генов
Для выявления общих генов среди дифференциально экспрессируемых генов (DEGs), генов модуля MEsalmon и ac4C-RGs с помощью пакета ggvenn для языка R была построена диаграмма Венна. Функциональный анализ обогащения общих генов проводили с использованием баз данных Gene Ontology (GO) и Kyoto Encyclopedia of Genes and Genomes (KEGG) с помощью пакета clusterProfiler для R. Информацию о белок-белковых взаимодействиях (PPI) получали из базы данных STRING (https://string-db.org/) для анализа взаимодействий на уровне белков между общими генами. Для визуализации полученной сети PPI использовали программное обеспечение Cytoscape.
Идентификация ключевых генов с помощью машинного обучения
Были применены три метода машинного обучения: регрессия методом наименьших абсолютных отклонений (LASSO), экстремальный градиентный бустинг (XGBoost) и случайный лес (RF). Регрессия LASSO была реализована с использованием пакета glmnet с 10-кратной кросс-валидацией для определения оптимального параметра штрафа λ. Параметр type.measure был установлен на "deviance", а параметр family — на "binomial". Оптимальное значение λ было выбрано по критерию λmin, который минимизирует отклонение при кросс-валидации, что позволило выделить 17 генов. XGBoost выполнялся с использованием пакета xgboost со следующими гиперпараметрами: nrounds = 100, max_depth = 6, eta = 0.3, subsample = 0.8, colsample_bytree = 0.8 и eval_metric = "logloss". Важность признаков ранжировалась по метрике gain, и были отобраны 30 наиболее значимых генов. Случайный лес был реализован с помощью пакета randomForest с ntree = 200. Важность признаков ранжировалась по среднему снижению индекса Джини, и были выбраны 30 наиболее значимых генов. Гены, отобранные тремя методами машинного обучения, были подвергнуты пересечению для выявления ключевых генов для последующего анализа.
Построение и оценка модели логистической регрессии для прогнозирования риска
Набор данных GSE54837 был случайным образом разделен на обучающую выборку (70%) и тестовую выборку (30%). На основе обучающей выборки была построена модель логистической регрессии с использованием функции glm из пакета MASS, где в качестве входных признаков использовались уровни экспрессии ключевых генов. Эффективность модели оценивали с помощью ROC-кривых, построенных с помощью пакета pROC. 95% доверительные интервалы для AUC были рассчитаны с помощью 2 000 бутстрэп-репликаций. Калибровку модели оценивали с помощью калибровочных кривых, построенных на основе 1 000 бутстрэп-выборок (пакет rms). Для оценки чистой клинической выгоды при различных пороговых значениях вероятности был проведен анализ кривой принятия решений (DCA) с использованием пакета dca. Для облегчения индивидуальной оценки риска с помощью функции nomogram из пакета rms была построена номограмма.
Уравнение регрессии имело следующий вид:
logit(P) = 0.5823 + 0.6010 × UPP1 - 0.6563 × PTRF + 0.3853 × B4GALT2 - 0.3972 × FAM168B + 0.1848 × PRKCDBP - 0.4787 × TOR3A. (1)
Здесь P представляет собой прогнозируемую вероятность ХОБЛ, а каждый коэффициент отражает вклад соответствующего значения экспрессии гена в логарифмический шанс развития ХОБЛ.
Анализ экспрессии, сеть GeneMANIA и молекулярная регуляторная сеть
Уровни экспрессии генов в группах с ХОБЛ и контрольной группе в наборе данных GSE54837 сравнивали с помощью критерия суммы рангов Вилкоксона. Для визуализации распределения уровней экспрессии с помощью пакета ggplot2 строили диаграммы размаха (box plots), на которых отображали медиану, межквартильный размах (IQR) и отдельные точки данных. Для построения генных сетей и прогнозирования функциональных взаимодействий использовали GeneMANIA. Поиск проводили с параметрами по умолчанию: вид = Homo sapiens, максимальное число связанных генов = 20. Полученную сеть загружали и визуализировали, при этом цвета ребер указывали на типы взаимодействий. Для изучения механизмов посттранскрипционной регуляции была построена сеть конкурентных эндогенных РНК (ceRNA). МикроРНК (miRNA), нацеленные на шесть ключевых генов, предсказывали с использованием двух независимых баз данных: DIANA-microT (score ≥ 0.8) и miRanda (score ≥ 140, энергия ≤ −20 kcal/mol). Пересечение микроРНК, выявленных обеими базами данных, использовали для построения пар miRNA-mRNA. Впоследствии с помощью базы данных StarBase предсказывали длинные некодирующие РНК (lncRNA), нацеленные на эти микроРНК. Регуляторная сеть lncRNA-miRNA-mRNA была построена и визуализирована с помощью Cytoscape. Транскрипционные регуляторные связи предсказывали с помощью ChIP-X Enrichment Analysis Version 3 (ChEA3). Для каждого ключевого гена с предсказанными TF выбирали 10 транскрипционных факторов с наивысшими показателями обогащения. В Cytoscape была построена регуляторная сеть TF-мишень.
Анализ обогащения наборов генов и оценка инфильтрации иммунных клеток
Анализ обогащения наборов генов (GSEA) был выполнен с помощью пакета clusterProfiler для изучения биологических функций каждого ключевого гена. Для каждого ключевого гена образцы были разделены на группы с высокой и низкой экспрессией на основе медианного значения. Анализ дифференциальной экспрессии между двумя группами проводили с помощью limma, и полученный список генов был ранжирован по знаковому log₂ fold-change. GSEA проводили с использованием функции gseGO для терминов биологических процессов GO и функции gseKEGG для путей KEGG со следующими параметрами: minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05 и nPerm = 1,000. Относительная численность 28 типов иммунных клеток оценивали с помощью анализа обогащения наборов генов для одного образца (ssGSEA), реализованного в пакете GSVA. Курируемая матрица сигнатур наборов генов, включающая маркерные гены для 28 типов иммунных клеток, была получена из предыдущих публикаций19. Для каждого образца применяли функцию gsva с параметрами method = "ssgsea", ssgsea.norm = TRUE и kcdf = "Gaussian". Коэффициенты корреляции Спирмена между показателями обогащения ssGSEA и уровнями экспрессии шести ключевых генов рассчитывали с помощью функции cor.test. Значения p были скорректированы с учетом множественного тестирования по методу Бенджамини-Хохберга. Матрица корреляции была визуализирована в виде тепловой карты с помощью пакета pheatmap.
Прогнозирование лекарственных препаратов, молекулярный докинг и анализ ассоциации с заболеваниями
С помощью базы данных DrugBank были определены потенциальные терапевтические соединения, воздействующие на ключевые гены. В программе Cytoscape была построена сеть взаимодействий «ключевой ген — таргетный препарат» для визуализации прогнозируемых взаимодействий между лекарственными средствами и генами. Для оценки аффинности связывания был проведен молекулярный докинг с использованием платформы CB-Dock2. Трехмерная структура белка UPP1 человека была получена из Protein Data Bank (PDB ID: 7B8T). Молекулярные структуры препаратов (в формате SMILES) были взяты из PubChem. Докинг выполнялся с помощью движка AutoDock Vina, а результаты ранжировались по свободной энергии связывания (ΔG, в kcal/mol). Комплексы докинга были визуализированы с помощью PyMOL. Связи между ключевыми генами и заболеваниями человека, связанными с воздействием факторов окружающей среды, исследовались с использованием базы данных Comparative Toxicogenomics Database (CTD). По каждому гену был проведен отдельный поиск, после чего десять наиболее сильно связанных с ними заболеваний были извлечены и визуализированы с помощью лепестковых диаграмм.
Протокол ОТ-кПЦР
Образцы периферической венозной крови были собраны у восьми пациентов с ХОБЛ и восьми здоровых добровольцев (контрольная группа) в Шэньчжэньской больнице традиционной китайской медицины района Луоху. Диагноз ХОБЛ ставили в соответствии с критериями Глобальной инициативы по хронической обструктивной болезни легких (GOLD), определяемыми как ОФВ1/ФЖЕЛ после применения бронхолитиков < 0,70. Контрольная группа состояла из здоровых добровольцев соответствующего пола и возраста без заболеваний дыхательных путей в анамнезе и с нормальными результатами тестов функции легких (ОФВ1% от должного ≥ 80% и ОФВ1/ФЖЕЛ ≥ 0,70). Исходные данные пациентов представлены в Таблице 2. Общую РНК экстрагировали из образцов крови пациентов с ХОБЛ с помощью набора для выделения РНК из крови. Для синтеза кДНК 500 ng общей РНК подвергли обратной транскрипции с использованием набора для синтеза кДНК с удалением геномной ДНК согласно прилагаемому протоколу. Полученную кДНК развели до концентрации 150 ng/μL.
RT-qPCR проводили с использованием мастер-микса для qPCR на основе SYBR Green на системе ПЦР в реальном времени. Каждая реакция объемом 10 μL содержала 5 μL 2x мастер-микса SYBR Green, по 0,5 μL прямого и обратного праймеров (10 μM), 1 μL разведенной кДНК (15 ng/μL) и 3 μL воды, свободной от нуклеаз. Условия циклирования включали начальную денатурацию при 95 °C в течение 5 min, затем 40 циклов: 95 °C в течение 10 s и 60 °C в течение 30 s, с последующим анализом кривых плавления от 60 °C до 95 °C для проверки специфичности амплификации. Все реакции проводили в трех технических повторностях. В качестве внутреннего референсного гена использовали β-actin. Эффективность праймеров для каждого целевого гена была подтверждена с использованием серий разведений по стандартной кривой и находилась в диапазоне от 90% до 110%. Уровни экспрессии генов нормировали по β-actin, а относительную экспрессию рассчитывали методом 2-ΔΔCt. Статистическое сравнение групп с ХОБЛ и контрольной группы проводили с помощью U-критерия Манна-Уитни.
Статистический анализ
Визуализация сетей была создана с помощью Cytoscape, а статистический анализ выполнялся с использованием программного обеспечения R. Если не указано иное, для данных с ненормальным распределением использовался U-критерий Манна-Уитни, а для данных с нормальным распределением — t-критерий Стьюдента для сравнения двух групп. Значение p < 0,05 считалось статистически значимым.