Déclaration du comité d'éthique institutionnel
Cette étude a été menée conformément à la Déclaration d'Helsinki. Le protocole a été approuvé par le Comité d'éthique de l'Hôpital Shenzhen Luohu de médecine traditionnelle chinoise (numéro d'approbation : 2024-LHQZYYYXLL-KY-039), et un consentement éclairé écrit a été obtenu de tous les participants avant leur inclusion. Les détails des outils et matériaux de recherche utilisés dans ce protocole sont fournis dans le Tableau des matériaux.
Source des données et traitement
Les jeux de données d'expression génique liés à la BPCO ont été obtenus à partir du Gene Expression Omnibus (GEO). Le jeu de données GSE54837 a été utilisé comme jeu de données de transcriptome, et le jeu de données GSE112811 a servi de jeu de validation (Tableau 1). Les gènes liés à l'ac4C (ac4C-RGs) ont été recueillis à partir de la littérature18. Les gènes différentiellement exprimés (DEG) entre les groupes BPCO et témoins ont été identifiés à l'aide du package R limma. Les DEG ont été considérés comme statistiquement significatifs si |log2FC| > 0 et p < 0,05. Des graphiques en volcan ont été générés pour visualiser la distribution globale des modifications d'expression génique.
Construction de WGCNA
L'analyse WGCNA a été réalisée sur le jeu de données GSE54837 à l'aide de R afin d'identifier les modules liés à la BPCO. Avant la construction du réseau, des échantillons aberrants ont été identifiés et éliminés par analyse de clustering hiérarchique utilisant la fonction hclust avec la méthode de liaison moyenne et une métrique de distance euclidienne. La puissance optimale de seuil doux (soft-thresholding) (β = 10) a été choisi pour obtenir un indice d'ajustement à une topologie sans échelle de R2 ≥ 0,85, équilibrant la topologie sans échelle et la connectivité moyenne. Une matrice d'adjacence a été construite puis transformée en une matrice de chevauchement topologique (TOM). Les modules de gènes ont été identifiés à l'aide de l'algorithme de découpage dynamique d'arbre hiérarchique (deepSplit = 2, minClusterSize = 50). Les modules présentant des corrélations d'expression des gènes caractéristiques > 0,75 ont ensuite été fusionnés à l'aide de la fonction mergeCloseModules. Les eigengènes des modules ont ensuite été corrélées aux caractéristiques cliniques (statut BPCO, âge, sexe et statut tabagique) à l'aide de coefficients de corrélation de Pearson afin d'identifier les modules associés au BPCO pour les analyses ultérieures.
Étude de criblage, analyse d'enrichissement et analyse du réseau PPI des gènes chevauchants
Un diagramme de Venn a été généré à l'aide du package R ggvenn afin d'identifier les gènes communs parmi les DEG, les gènes du module MEsaumon et les ac4C-RG. Une analyse d'enrichissement fonctionnel des gènes communs a été réalisée à l'aide des bases de données Gene Ontology (GO) et Kyoto Encyclopedia of Genes and Genomes (KEGG) avec le package R clusterProfiler. Les informations sur les interactions entre protéines (PPI) ont été obtenues à partir de la base de données STRING (https://string-db.org/) afin d'analyser les interactions au niveau protéique parmi les gènes communs. Le logiciel Cytoscape a été utilisé pour visualiser le réseau PPI obtenu.
Identification de gènes clés par apprentissage automatique
Trois techniques d'apprentissage automatique ont été appliquées : la régression par l'opérateur de réduction et de sélection des moindres valeurs absolues (LASSO), le boosting par gradient extrême (XGBoost) et la forêt aléatoire (RF). La régression LASSO a été mise en œuvre à l'aide du package glmnet avec une validation croisée à 10 folds pour déterminer le paramètre de pénalité optimal λ. Le paramètre type.measure a été fixé à « deviance », et le paramètre family a été fixé à « binomial ». Le λ optimal a été sélectionné selon le critère λmin, qui minimise la déviance validée par validation croisée, produisant 17 gènes. XGBoost a été réalisé à l'aide du package xgboost avec les hyperparamètres suivants : nrounds = 100, max_depth = 6, eta = 0,3, subsample = 0,8, colsample_bytree = 0,8, et eval_metric = « logloss ». L'importance des variables a été classée selon la métrique gain, et les 30 gènes les plus importants ont été sélectionnés. La forêt aléatoire a été mise en œuvre à l'aide du package randomForest avec ntree = 200. L'importance des variables a été classée selon la diminution moyenne de l'indice de Gini, et les 30 gènes les plus importants ont été sélectionnés. Les gènes sélectionnés par les trois méthodes d'apprentissage automatique ont été croisés afin d'identifier les gènes clés pour les analyses ultérieures.
Construction et évaluation du modèle de régression logistique pour la prédiction du risque
Le jeu de données GSE54837 a été divisé aléatoirement en un ensemble d'apprentissage (70 %) et un ensemble de test (30 %). Un modèle de régression logistique a été construit sur l'ensemble d'apprentissage à l'aide de la fonction glm du package MASS, en utilisant les niveaux d'expression des gènes clés comme variables prédictives. Les performances du modèle ont été évaluées à l'aide de courbes ROC générées avec le package pROC. Les intervalles de confiance à 95 % pour l'AUC ont été calculés par 2 000 répétitions de bootstrap. La calibration du modèle a été évaluée à l'aide de courbes de calibration générées avec 1 000 rééchantillonnages bootstrap (package rms). L'analyse de décision clinique (DCA) a été réalisée à l'aide du package dca pour évaluer le bénéfice clinique net sur une gamme de probabilités seuils. Un nomogramme a été construit à l'aide de la fonction nomogram du package rms afin de faciliter l'estimation personnalisée du risque.
L'équation de régression était la suivante :
logit(P) = 0,5823 + 0,6010 × UPP1 - 0,6563 × PTRF + 0,3853 × B4GALT2 - 0,3972 × FAM168B + 0,1848 × PRKCDBP - 0,4787 × TOR3A. (1)
Ici, P représente la probabilité prédite de BPCO, et chaque coefficient représente la contribution de la valeur d'expression génique correspondante aux cotes logarithmiques de la BPCO.
Analyse d'expression, réseau GeneMANIA et réseau régulatoire moléculaire
Les niveaux d'expression génique entre les groupes BPCO et témoins dans le jeu de données GSE54837 ont été comparés à l'aide du test de Wilcoxon-Mann-Whitney. Des boîtes à moustaches ont été générées à l'aide du package ggplot2 afin de visualiser la distribution des niveaux d'expression, avec la médiane, l'intervalle interquartile (IIQ) et les points de données individuels superposés. GeneMANIA a été utilisé pour construire les réseaux géniques et prédire les interactions fonctionnelles. La recherche a été effectuée avec les paramètres par défaut : espèce = Homo sapiens, nombre maximal de gènes associés = 20. Le réseau obtenu a été téléchargé et visualisé, les couleurs des arêtes indiquant les types d'interaction. Un réseau d'ARN endogènes compétitifs (ceRNA) a été construit afin d'étudier les mécanismes régulateurs post-transcriptionnels. Les miARN ciblant les six gènes clés ont été prédits à l'aide de deux bases de données indépendantes : DIANA-microT (score ≥ 0,8) et miRanda (score ≥ 140, énergie ≤ −20 kcal/mol). L'intersection des miARN identifiés par les deux bases de données a été utilisée pour construire des paires miARN-mARN. Par la suite, les lncARN ciblant ces miARN ont été prédits à l'aide de la base de données StarBase. Un réseau régulateur lncARN-miARN-mARN a été construit et visualisé à l'aide de Cytoscape. Les relations régulatrices transcriptionnelles ont été prédites à l'aide de l'analyse d'enrichissement ChIP-X version 3 (ChEA3). Pour chaque gène clé possédant des facteurs de transcription (FT) prédits, les 10 premiers facteurs de transcription présentant les scores d'enrichissement les plus élevés ont été sélectionnés. Un réseau régulateur FT-cible a été construit dans Cytoscape.
Analyse de l'enrichissement des ensembles de gènes et évaluation de l'infiltration des cellules immunitaires
Une analyse d'enrichissement de gènes (GSEA) a été réalisée à l'aide du package clusterProfiler afin d'étudier les fonctions biologiques de chaque gène clé. Pour chaque gène clé, les échantillons ont été divisés en deux groupes, à expression élevée et faible, selon la valeur médiane. Une analyse d'expression différentielle entre les deux groupes a été effectuée à l'aide de limma, et la liste de gènes obtenue a été classée selon le changement de rapport logarithmique en base 2 signé. L'analyse GSEA a été menée à l'aide de la fonction gseGO pour les termes du processus biologique de l'ontologie générique (GO) et de la fonction gseKEGG pour les voies KEGG, avec les paramètres suivants : minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0,05, et nPerm = 1 000. L'abondance relative de 28 types de cellules immunitaires a été estimée à l'aide de l'analyse d'enrichissement de gènes par échantillon unique (ssGSEA) mise en œuvre dans le package GSVA. Une matrice de signature de gènes soigneusement sélectionnée, comprenant des gènes marqueurs pour 28 types de cellules immunitaires, a été obtenue à partir de travaux antérieurs19. Pour chaque échantillon, la fonction gsva a été appliquée avec method = "ssgsea", ssgsea.norm = TRUE, et kcdf = "Gaussian". Les coefficients de corrélation de Spearman entre les scores d'enrichissement ssGSEA et les niveaux d'expression des six gènes clés ont été calculés à l'aide de la fonction cor.test. Les valeurs de p ont été ajustées pour les tests multiples selon la méthode de Benjamini-Hochberg. La matrice de corrélation a été visualisée sous forme de carte thermique à l'aide du package pheatmap.
Prédiction de médicaments, docking moléculaire et analyse d'association avec les maladies
Des composés thérapeutiques potentiels ciblant des gènes clés ont été identifiés à l'aide de la base de données DrugBank. Un réseau d'interactions « médicament ciblant un gène clé » a été construit dans Cytoscape afin de visualiser les interactions prévues entre médicaments et gènes. Un dockage moléculaire a été réalisé à l'aide de la plateforme CB-Dock2 pour évaluer les affinités de liaison. La structure protéique tridimensionnelle de l'UPP1 humain a été obtenue à partir de la Protein Data Bank (PDB ID : 7B8T). Les structures moléculaires des médicaments (format SMILES) ont été récupérées depuis PubChem. Le dockage a été effectué en utilisant le moteur AutoDock Vina, et les résultats ont été classés selon l'énergie libre de liaison (ΔG, en kcal/mol). Les complexes de dockage ont été visualisés à l'aide de PyMOL. Les associations entre les gènes clés et les maladies humaines liées à des expositions environnementales ont été étudiées à l'aide de la base de données Comparative Toxicogenomics Database (CTD). Chaque gène a fait l'objet d'une requête individuelle, et les dix maladies les plus fortement associées ont été extraites et visualisées à l'aide de diagrammes en radar.
Protocole de RT-qPCR
Des échantillons de sang veineux périphérique ont été prélevés chez huit patients atteints de BPCO et huit témoins sains à l'hôpital Shenzhen Luohu de médecine traditionnelle chinoise. Le diagnostic de BPCO a été établi selon les critères de l'Initiative mondiale pour la maladie pulmonaire obstructive chronique (GOLD), définis par un rapport VEMS/CVF post-bronchodilatateur < 0,70. Le groupe témoin comprenait des volontaires sains appariés par âge et sexe, sans antécédent de maladies respiratoires et présentant des épreuves fonctionnelles respiratoires normales (VEMS % prédit ≥ 80 % et VEMS/CVF ≥ 0,70). Les informations de base concernant les patients sont indiquées dans le Tableau 2. L'ARN total a été extrait des échantillons sanguins de patients BPCO à l'aide d'un kit d'extraction d'ARN sanguin. Pour la synthèse d'ADN complémentaire, 500 ng d'ARN total ont été rétrotranscrits à l'aide d'un kit de synthèse d'ADN complémentaire avec élimination de l'ADN génomique, conformément au protocole fourni. L'ADN complémentaire obtenu a été dilué à 150 ng/μL.
La RT-qPCR a été réalisée à l'aide d'un mélange maître qPCR basé sur le SYBR Green sur un système de PCR en temps réel. Chaque réaction de 10 μL contenait 5 μL de mélange maître 2x SYBR Green, 0,5 μL de chacun des amorces avant et arrière (10 μM), 1 μL d'ADNc dilué (15 ng/μL) et 3 μL d'eau sans nucléase. Les conditions de cyclage comprenaient une dénaturation initiale à 95 °C pendant 5 min, suivie de 40 cycles de 95 °C pendant 10 s et de 60 °C pendant 30 s, avec une analyse finale de courbe de fusion allant de 60 °C à 95 °C afin de vérifier la spécificité de l'amplification. Toutes les réactions ont été effectuées en triplets techniques. L'actine β a été utilisée comme gène de référence interne. L'efficacité des amorces pour chaque gène cible a été validée à l'aide de séries de dilutions par courbe standard et variait entre 90 % et 110 %. Les niveaux d'expression génique ont été normalisés par rapport à l'actine β, et l'expression relative a été calculée selon la méthode 2-ΔΔCt. Les comparaisons statistiques entre les groupes BPCO et témoins ont été effectuées à l'aide du test de Mann-Whitney U.
Analyse statistique
Les visualisations de réseaux ont été créées à l'aide de Cytoscape, et les analyses statistiques ont été effectuées avec le logiciel R. Sauf indication contraire, le test de Mann-Whitney U a été utilisé pour les données non normalement distribuées, et le test t de Student a été utilisé pour les données normalement distribuées afin de comparer deux groupes. Une valeur de p < 0,05 a été considérée comme statistiquement significative.