Сбор данных
Матрица экспрессии генов и соответствующие клинические данные для наборов данных GSE142025 и GSE96804, относящихся к диабетической нефропатии (ДН), были получены из базы данных Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/)15. Обучающая выборка (GSE142025) включала образцы тканей почек 27 пациентов с ДН и девять контрольных образцов, секвенирование которых проводилось с использованием платформы GPL20301. Валидационная выборка (GSE96804) состояла из данных секвенирования 41 пациента с ДН и 20 контрольных образцов, обработанных с помощью платформы GPL17586. Оба набора данных сосредоточены на тканях почек; набор данных GSE96804 специально исследует клубочек — основной фильтрационного звено почки (Рисунок 1). Набор данных GSE142025 (обучающая выборка) включает образцы всей ткани почек и обеспечивает обширный обзор транскриптомного ландшафта ДН. С другой стороны, набор данных GSE96804 (валидационная выборка) сфокусирован конкретно на ткани клубочков, которая является основным местом повреждения при гломерулярной фильтрации. Поскольку эти два набора данных не объединялись напрямую из-за различий в платформах и типах тканей, коррекция эффекта партии не применялась. Вместо этого валидация между наборами данных проводилась независимо. Гены с согласованными направленными изменениями и площадью под кривой (AUC) более 0,7 в обоих наборах данных были выбраны в качестве надежных кандидатов, что подтверждает возможность обобщения результатов для различных отделов почек.
В общей сложности 1 136 генов, связанных с митохондриями (MRGs), были извлечены из базы данных MitoCarta3.0 (https://www.broadinstitute.org/mitocarta). Активные ингредиенты LLF были определены с помощью базы данных Traditional Chinese Medicine Systems Pharmacology (TCMSP) (http://sm.nwsuaf.edu.cn/lsp/tcmsp.php) на основе пороговых значений биодоступности при пероральном приеме (OB) ≥30% и лекарственности (DL) ≥0,18. Впоследствии потенциальные гены-мишени для активных компонентов были предсказаны с использованием базы данных Swiss Target Prediction (http://www.swisstargetprediction.ch/).
Анализ дифференциальной экспрессии
Анализ дифференциальной экспрессии GSE142025 (ДН против контроля) проводили с помощью пакета limma (v3.54.1), при этом критерии значимости были установлены как P.adj < 0.05 и |log2FoldChange| > 0.516 . Для визуализации графиков «вулкан» (volcano plots) и тепловых карт использовали пакеты ggplot2 (v 3.3.6) и ComplexHeatmap (v 2.14.0) соответственно17,18. Было выполнено пересечение дифференциально экспрессируемых генов (DEGs), генов, связанных с метаболизмом (MRGs), и потенциальных генов-мишеней активных ингредиентов; перекрывающиеся гены были определены как гены-кандидаты. Сеть, связывающую активные ингредиенты и гены-кандидаты, построили с помощью программного обеспечения Cytoscape (v 3.9.0)19.
Функциональный анализ обогащения и построение сети белок-белковых взаимодействий (PPI)
Анализ обогащения по Gene Ontology (GO) и Kyoto Encyclopedia of Genes and Genomes (KEGG) для исследуемых генов-кандидатов был выполнен с использованием пакета clusterProfiler (v 4.6.2) для изучения их биологических функций и связанных с ними сигнальных путей (P.adjust < 0.05). Затем список генов-кандидатов был введен в базу данных STRING (https://cn.string-db.org/) для получения данных о белок-белковых взаимодействиях (уровень достоверности ≥ 0.4), после чего с помощью Cytoscape (v 3.9.0) была построена сеть PPI20.
Машинное обучение
Четыре алгоритма машинного обучения, включая случайный лес (RF), k-ближайших соседей (KNN), метод частичных наименьших квадратов (PLS) и метод опорных векторов с радиально-базисной функцией ядра (SVM), были реализованы с помощью пакета caret (v6.0-93) на основе набора данных GSE142025. В качестве предикторов использовались гены-кандидаты, выявленные в ходе предыдущего анализа, а в качестве целевой переменной — статус заболевания (DN или контроль). Для модели KNN была применена 10-кратная перекрестная проверка с использованием функции trainControl при tuneLength = 10. Модель RF была построена с использованием 20 деревьев (ntree = 20); модели PLS и SVM были настроены согласно параметрам, реализованным в исходном рабочем процессе caret.
Распределения остатков модели оценивали с помощью пакета DALEX (v2.4.3). Затем была проведена оценка переменных на основе перестановок с использованием функции variable_importance пакета DALEX, где в качестве функции потерь была указана среднеквадратическая ошибка (RMSE). Полученное значение dropout_loss представляет собой RMSE, полученное после перестановки соответствующей переменной21,22. В исходном аналитическом рабочем процессе были сохранены переменные со значением dropout loss < 0.281 для всех моделей, а гены, общие для всех четырех моделей, были определены как консенсусные гены-признаки для последующей валидации.
Четыре модели машинного обучения использовались преимущественно для приоритизации признаков, а не для построения окончательного клинического классификатора. Соответственно, диагностическая разделительная способность впоследствии оценивалась на уровне отдельных генов с помощью ROC-анализа как в открывающем, так и в валидационном наборах данных.
Идентификация потенциальных биомаркеров
В наборах данных GSE142025 и GSE96804 различия в экспрессии исследуемых генов между образцами с диабетической нефропатией (ДН) и контрольными образцами оценивали с помощью критерия Вилкоксона. Для анализа рабочих характеристик приемника (ROC-кривых) отбирали только те гены, которые имели значимую дифференциальную экспрессию (P < 0.05) и согласованные тенденции экспрессии в обоих наборах данных. Для построения ROC-кривых и расчета площади под кривой (AUC) использовали пакет pROC (v1.18.0); гены, демонстрирующие AUC > 0.7 в обоих наборах данных, классифицировали как потенциальные биомаркеры23.
Анализ обогащения наборов генов (GSEA)
Биологические функции и сигнальные пути, связанные с кандидатными биомаркерами, были дополнительно изучены с помощью GSEA на наборе данных GSE142025. Сначала с использованием пакета psych (v2.2.9) был проведен корреляционный анализ Спирмена кандидатных биомаркеров со всеми остальными генами.24Коэффициенты корреляции были рассчитаны и ранжированы (от высокого к низкому). В качестве эталонного набора генов использовался c2.cp.kegg.v2023.1.Hs.symbols.gmt из базы данных молекулярных сигнатур (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/). Затем с помощью пакета clusterProfiler (v4.6.2) был проведен анализ GSEA для оценки обогащения ранжированных генов в фоновом наборе генов. Была применена поправка на множественное тестирование. посредством метод FDR и скорректированные P-значения (обозначенные как P.adjust) считались значимыми, если < 0.05.
Анализ модификации m6A
Для исследования модификаций метилирования РНК потенциальных биомаркеров использовалась база данных SRAMP (http://www.cuilab.cn/sramp/) для предсказания сайтов модификации m6A в последовательностях биомаркеров с акцентом на высокодостоверные позиции в их вторичных структурах. Затем с помощью базы данных ENCORI (https://starbase.sysu.edu.cn/) были идентифицированы белки, модифицированные m6A, которые взаимодействуют с потенциальными биомаркерами; для скрининга ключевых белков использовался параметр |HepG2 (shRNA)| > 1. Впоследствии использовалась база данных RPISeq (http://pridb.gdcb.iastate.edu/RPISeq/) для предсказания вероятности взаимодействий между ключевыми белками и потенциальными биомаркерами. Последовательности РНК обоих компонентов были загружены в формате простого текста для получения оценок предсказания классификаторами RF и SVM. Взаимодействие считалось значимым, если оценка превышала 0,525. Анализ SRAMP проводился с порогом «High» для предсказания сайтов m6A в режиме «Transcript» и с параметрами по умолчанию. Для анализа в ENCORI использовалась функция взаимодействия «miRNA-mRNA» с параметром «HepG2 (shRNA)» > 1. Анализ RPISeq проводился с использованием классификатора RF с параметрами по умолчанию; значения > 0,5 указывали на положительное взаимодействие. Данные результаты являются вычислительными предсказаниями, а не экспериментальными доказательствами модификации m6A или РНК-белковых взаимодействий в тканях почек. Критерий HepG2 shRNA был получен из предварительно рассчитанных наборов данных ENCORI и может не отражать специфическую регуляцию в почках.
Анализ иммунной инфильтрации
Для оценки пропорций 22 типов иммунных клеток в контрольных образцах и образцах с ДН из набора GSE142025 был применен алгоритм CIBERSORT, а визуализация была выполнена с помощью тепловой карты, созданной с использованием пакета ggplot2 (v3.3.6)26. CIBERSORT запускали с использованием сигнатурной матрицы LM22, 1000 перестановок и отключенной квантильной нормализацией (как рекомендуется для данных микрочипов). Для дальнейшего анализа были отобраны образцы с p < 0.05 по результатам CIBERSORT. CIBERSORT оценивает фракции иммунных клеток на основе экспрессии в общем объеме ткани почки, что не позволяет разделить инфильтрацию по компартментам (например, гломерулярную и тубулоинтерстициальную), а также отличить инфильтрирующие лейкоциты от резидентных иммунных клеток. Таким образом, сообщаемые корреляции относятся к уровню ткани и должны быть подтверждены пространственными методами. Корреляционный анализ Спирмена между иммунными клетками с дифференциальным содержанием и потенциальными биомаркерами проводили с использованием пакета psych.
Построение сетей и молекулярный докинг
МикроРНК (miRNAs), взаимодействующие с потенциальными биомаркерами, были предсказаны с использованием базы данных miRNet (https://www.mirnet.ca). Затем с помощью баз данных TarBase (http://www.diana.pcbi.upenn.edu/tarbase) и starbase (http://starbase.sysu.edu.cn/) были предсказаны длинные некодирующие РНК (lncRNAs), мишенями которых являются выявленные микроРНК. Для построения сети были выбраны lncRNAs, общие для обеих баз данных. После этого с помощью программного обеспечения Cytoscape была построена регуляторная сеть lncRNA-miRNA-mRNA. Потенциальные активные ингредиенты, воздействующие на потенциальные биомаркеры, были отобраны для построения сети «активный ингредиент — биомаркер». Кроме того, активные ингредиенты, потенциальные биомаркеры и пути, выявленные с помощью GSEA, были внесены в Cytoscape для создания сети «активный ингредиент — биомаркер — путь».
Для оценки аффинности связывания между потенциальными активными ингредиентами и кандидатными биомаркерами был проведен анализ молекулярного докинга. 3D-структуры белков-биомаркеров были получены из банка данных белков Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB, https://www.rscb.org/pdb) в формате PDB. 2D-структуры потенциальных активных ингредиентов были извлечены в формате SDF из базы данных PubChem (http://pubchem.ncbi.nlm.nih.gov). Молекулярный докинг выполняли с использованием платформы CB-Dock (http://clab.labshare.cn/cb-dock/php/blinddock.php). Значение энергии связывания менее -5 kcal/mol указывало на сильную аффинность связывания27.
Подготовка и аутентификация Ligustri Lucidi Fructus
В данном случае Ligustri Lucidi Fructus (LLF) представляет собой высушенные спелые плоды Ligustrum lucidum W. T. Aiton (Oleaceae). Ботанический материал был аутентифицирован в соответствии с Фармакопеей Китая; гербарные образцы № 20240506, 20240911, 20241103 были депонированы в Университете традиционной китайской медицины провинции Шаньси.
Для приготовления отвара 200 g соответствующих требованиям ломтиков LLF замачивали в 1,000 mL дистиллированной воды в течение 30 мин при комнатной температуре. Смесь доводили до интенсивного кипения, затем томили на слабом огне в течение 60 мин. Собирали фильтрат, а оставшийся растительный материал повторно кипятили с еще 1,000 mL дистиллированной воды в течение 60 мин. Два полученных фильтрата объединяли, фильтровали, центрифугировали и концентрировали при пониженном давлении до достижения конечной концентрации стока 1 g сырого лекарственного вещества/mL (общий объем 100 mL). Приготовленный отвар хранили при 4 °C для краткосрочного использования или при −20 °C для длительного хранения. Качество LLF и полученного отвара строго проверяли и подтверждали в соответствии со стандартами Фармакопеи Китайской Народной Республики для обеспечения надежности и воспроизводимости экспериментов.
Для качественной идентификации была проведена тонкослойная хроматография. Вкратце: соответствующий объем приготовленного отвара центрифугировали, а супернатант экстрагировали метанолом. После фильтрации раствор образца и раствор стандартного образца specnuezhenide наносили на одну и ту же пластину с силикагелем G. После проявления, высушивания и осмотра под ультрафиолетовым светом пятно раствора образца по цвету и положению совпало с пятном стандартного соединения, что подтвердило наличие характерного биоактивного компонента LLF.
Для количественного контроля качества была проведена детекция методом высокоэффективной жидкостной хроматографии. Анализ выполняли с использованием колонки C18, где в качестве подвижной фазы применяли смесь метанола и воды. Длина волны детектирования была установлена на уровне 224 nm. Содержание specnuezhenide в отваре определяли на основании стандартной кривой. Результаты продемонстрировали стабильный и однородный химический состав приготовленного отвара, что обеспечило постоянное качество препарата на протяжении всего эксперимента с применением животных.
Эксперименты на животных
Двенадцать самцов мышей линии db/db (возраст 8–9 недель), соответствующих стандарту SPF, и шесть контрольных самцов мышей линии db/m того же возраста содержались в SPF-виварии Шаньсийского университета традиционной китайской медицины. Перед началом экспериментов животные проходили период акклиматизации в течение 7 дней при 12-часовом цикле свет/темнота с ad libitum доступ к пище и воде. Исследование было одобрено Комитетом по этике Шаньсийского университета традиционной китайской медицины (номер одобрения 2022DW167). У животных, проявлявших >Животных с потерей массы тела на 20%, в состоянии агонии или при невозможности доступа к пище и воде гуманно усыпляли до достижения запланированной конечной точки исследования. По завершении исследования все оставшиеся мыши были усыплены путем внутрибрюшинного введения пентобарбитала натрия с последующим переломом шейных позвонков в соответствии с протоколами учреждения.
После периода акклиматизации создание модели диабетической нефропатии (ДН) у мышей db/db подтверждалось уровнем глюкозы в крови из хвостовой вены ≥ 16 mmol/L и микроальбуминурией, определяемой по положительному результату теста с тест-полоской на микроальбумин в моче. После успешного создания модели ДН мыши db/db были случайным образом разделены на две группы (n = 6 в каждой группе): группу модели ДН (DN) и группу лечения LLF (Treatment). Кроме того, мыши db/m (n = 6) использовались в качестве контрольной группы (Control). Дозировка была выбрана на основе предыдущих фармакодинамических исследований LLF на диабетических крысах и пересчитана на эквивалентную дозу для человека с использованием нормализации по площади поверхности тела28. Группам Control и DN вводили дистиллированную воду, тогда как группа лечения получала LLF в дозе 3,5 g/kg в течение 8 недель. После 8 недель введения всех мышей усыпляли для сбора сыворотки крови, мочи и тканей почек для последующих анализов.
Показатели крови и мочи
Уровни глюкозы в сыворотке крови анализировали с помощью полностью автоматического биохимического анализатора крови. Концентрацию микроальбумина в моче измеряли в соответствии с инструкцией к набору (Дополнительный файл 1).
Патологическое исследование тканей почек мыши
Ткани почек были подготовлены для гистопатологического исследования. После фиксации в 4% параформальдегиде ткани промывали, дегидратировали, заливали в парафин и нарезали на срезы. Затем проводили окрашивание гематоксилином и эозином (HE), после чего окрашенные срезы изучали под оптическим микроскопом для оценки патологических изменений.
Обратная транскрипция и количественная полимеразная цепная реакция (ОТ-кПЦР)
ОТ-кПЦР использовали для определения экспрессии потенциальных биомаркеров в ткани почек мыши. Тотальную РНК экстрагировали согласно инструкциям производителя, после чего оценивали концентрацию и качество РНК (Таблица 1). кДНК синтезировали из экстрагированной РНК с помощью набора для синтеза кДНК. Амплификацию проводили с использованием пар праймеров, перечисленных в Таблице 1>, используя GAPDH в качестве референсного гена. Относительные уровни экспрессии рассчитывали методом 2−ΔΔCt13,26.
Статистический анализ
Все статистические анализы были проведены с использованием программного обеспечения R (версия 4.2.2) совместно с ПО, необходимым для соответствующих экспериментальных процедур. Если не указано иное, все статистические тесты были двусторонними, а различия считались статистически значимыми при P < 0.05. Анализ дифференциальной экспрессии транскриптома выполнялся с использованием пакета limma. Гены с скорректированным значением P < 0.05 и абсолютным значением log2-кратного изменения более 0,5 были определены как дифференциально экспрессируемые.
Для сравнения экспрессии генов-кандидатов между независимыми образцами с DN и контрольными образцами использовались непараметрические критерии суммы рангов Уилкоксона, где это было указано в исходном аналитическом протоколе. Корреляции между биомаркерами-кандидатами и фракциями иммунных клеток оценивали с помощью ранговой корреляции Спирмена.
Экспериментальные данные представлены как среднее ± SD. Сравнения между тремя независимыми группами животных проводились с помощью однофакторного дисперсионного анализа при соблюдении условий параметрического анализа. Для апостериорных сравнений использовался критерий наименьшей значимой разности при условии гомогенности дисперсий, тогда как при неоднородности дисперсий применялся T3-тест Даннетта. Исходные фармакодинамические показатели и показатели на 8-й неделе анализировались и представлялись раздельно; выводы о взаимодействии факторов группы и времени не делались. Скрипты R и исходные данные, использованные для биоинформатического анализа и анализа с применением машинного обучения, представлены в Дополнительном файле 1.