Dichiarazione del comitato etico
Lo studio è stato condotto in conformità con la Dichiarazione di Helsinki. Il protocollo è stato approvato dal Comitato Etico dell'Ospedale Shenzhen Luohu di Medicina Tradizionale Cinese (numero di approvazione 2024-LHQZYYYXLL-KY-039), e il consenso informato scritto è stato ottenuto da tutti i partecipanti prima dell'arruolamento. I dettagli degli strumenti e dei materiali di ricerca utilizzati in questo protocollo sono riportati nella Tabella dei Materiali.
Fonte dei dati ed elaborazione
I set di dati sull'espressione genica correlati alla BPCO sono stati ottenuti dal Gene Expression Omnibus (GEO). Il dataset GSE54837 è stato utilizzato come dataset del trascrittoma, mentre il dataset GSE112811 ha funto da set di validazione (Tabella 1). I geni ac4C-RG sono stati raccolti dalla letteratura18. I geni differenzialmente espressi (DEG) tra i gruppi BPCO e controllo sono stati identificati utilizzando il pacchetto R limma. I DEG sono stati considerati statisticamente significativi se |log2FC| > 0 e p < 0,05. Sono stati generati grafici a vulcano per visualizzare la distribuzione complessiva delle variazioni dell'espressione genica.
Costruzione di WGCNA
È stata eseguita un'analisi WGCNA sul set di dati GSE54837 utilizzando R per identificare moduli associati alla BPCO. Prima della costruzione della rete, i campioni outlier sono stati identificati ed eliminati mediante analisi di clustering gerarchico utilizzando la funzione hclust con metodo di collegamento medio e una metrica di distanza euclidea. È stata selezionata la potenza ottimale di soglia morbida (β = 10) per ottenere un indice di adattamento alla topologia libera da scala R2 ≥ 0,85, bilanciando topologia libera da scala e connettività media. È stata costruita una matrice di adiacenza e trasformata in una matrice di sovrapposizione topologica (TOM). I moduli genici sono stati identificati utilizzando l'algoritmo dinamico di taglio dell'albero (deepSplit = 2, minClusterSize = 50). I moduli con correlazioni di eigengene > 0,75 sono stati successivamente uniti utilizzando la funzione mergeCloseModules. Gli eigengene dei moduli sono stati quindi correlati con le caratteristiche cliniche (stato di BPCO, età, sesso e abitudine al fumo) mediante coefficienti di correlazione di Pearson per identificare i moduli associati alla BPCO da sottoporre ad analisi successive.
Selezione, analisi dell'arricchimento e analisi della rete PPI dei geni sovrapposti
È stato generato un diagramma di Venn utilizzando il pacchetto R ggvenn per identificare i geni sovrapposti tra i DEG, i geni del modulo MEsalmon e gli ac4C-RG. L'analisi di arricchimento funzionale dei geni sovrapposti è stata eseguita utilizzando le basi di dati Gene Ontology (GO) e Kyoto Encyclopedia of Genes and Genomes (KEGG) con il pacchetto R clusterProfiler. Le informazioni sulle interazioni proteina-proteina (PPI) sono state ottenute dal database STRING (https://string-db.org/) per analizzare le interazioni a livello proteico tra i geni sovrapposti. Il software Cytoscape è stato utilizzato per visualizzare la rete PPI risultante.
Identificazione di geni chiave mediante apprendimento automatico
Sono state applicate tre tecniche di machine learning: regressione con operatore di riduzione e selezione della norma L1 (LASSO), boosting estremo del gradiente (XGBoost) e foresta casuale (RF). La regressione LASSO è stata implementata utilizzando il pacchetto glmnet con validazione incrociata a 10 ripetizioni per determinare il parametro di penalizzazione ottimale λ. Il parametro type.measure è stato impostato su «deviance», e il parametro family su «binomial». Il valore ottimale di λ è stato selezionato mediante il criterio λmin, che minimizza la devianza ottenuta dalla validazione incrociata, producendo 17 geni. XGBoost è stato eseguito utilizzando il pacchetto xgboost con i seguenti iperparametri: nrounds = 100, max_depth = 6, eta = 0.3, subsample = 0.8, colsample_bytree = 0.8 e eval_metric = «logloss». L'importanza delle caratteristiche è stata ordinata in base alla metrica gain, selezionando i 30 geni più significativi. La foresta casuale è stata implementata utilizzando il pacchetto randomForest con ntree = 200. L'importanza delle caratteristiche è stata valutata in base alla riduzione media dell'indice di Gini, selezionando i 30 geni più importanti. I geni selezionati dai tre metodi di machine learning sono stati incrociati per identificare i geni chiave da utilizzare nelle analisi successive.
Costruzione e valutazione del modello di regressione logistica per la previsione del rischio
Il dataset GSE54837 è stato suddiviso casualmente in un insieme di addestramento (70%) e un insieme di test (30%). Un modello di regressione logistica è stato costruito sull'insieme di addestramento utilizzando la funzione glm del pacchetto MASS, con i livelli di espressione dei geni chiave come caratteristiche di input. Le prestazioni del modello sono state valutate mediante curve ROC generate con il pacchetto pROC. Gli intervalli di confidenza al 95% per l'AUC sono stati calcolati tramite 2.000 replicati bootstrap. La calibrazione del modello è stata valutata utilizzando curve di calibrazione generate con 1.000 ricampionamenti bootstrap (pacchetto rms). L'analisi decisionale clinica (DCA) è stata eseguita utilizzando il pacchetto dca per valutare il beneficio clinico netto in un intervallo di probabilità soglia. Un nomogramma è stato costruito utilizzando la funzione nomogram del pacchetto rms per facilitare la stima del rischio individualizzata.
L'equazione di regressione era:
logit(P) = 0,5823 + 0,6010 × UPP1 - 0,6563 × PTRF + 0,3853 × B4GALT2 - 0,3972 × FAM168B + 0,1848 × PRKCDBP - 0,4787 × TOR3A. (1)
In questo caso, P rappresenta la probabilità predetta di BPCO, e ciascun coefficiente rappresenta il contributo del corrispondente valore di espressione genica alle probabilità logaritmiche di BPCO.
Analisi dell'espressione, rete GeneMANIA e rete regolatoria molecolare
I livelli di espressione genica tra i gruppi COPD e di controllo nel dataset GSE54837 sono stati confrontati utilizzando il test della somma dei ranghi di Wilcoxon. I grafici a scatola (box plot) sono stati generati con il pacchetto ggplot2 per visualizzare la distribuzione dei livelli di espressione, con mediana, intervallo interquartile (IQR) e punti dati individuali sovrapposti. GeneMANIA è stato utilizzato per costruire reti geniche e prevedere interazioni funzionali. La ricerca è stata eseguita con parametri predefiniti: specie = Homo sapiens, geni correlati massimi = 20. La rete risultante è stata scaricata e visualizzata, con i colori dei collegamenti (edge) che indicano i tipi di interazione. È stata costruita una rete di RNA endogeno competitivo (ceRNA) per indagare i meccanismi regolatori post-trascrizionali. Le miRNA che bersagliano i sei geni chiave sono state previste utilizzando due database indipendenti: DIANA-microT (punteggio ≥ 0,8) e miRanda (punteggio ≥ 140, energia ≤ −20 kcal/mol). L'intersezione delle miRNA identificate da entrambi i database è stata utilizzata per costruire coppie miRNA-mRNA. Successivamente, le lncRNA che bersagliano queste miRNA sono state previste utilizzando il database StarBase. È stata costruita e visualizzata una rete regolatoria lncRNA-miRNA-mRNA mediante Cytoscape. Le relazioni regolatorie trascrizionali sono state previste utilizzando l'analisi di arricchimento ChIP-X Versione 3 (ChEA3). Per ciascun gene chiave con fattori trascrizionali (TF) previsti, sono stati selezionati i primi 10 fattori trascrizionali con i punteggi di arricchimento più elevati. Una rete regolatoria TF-target è stata costruita in Cytoscape.
Analisi dell'arricchimento di geni e valutazione dell'infiltrazione delle cellule immunitarie
È stata eseguita un'analisi di arricchimento dei set di geni (GSEA) utilizzando il pacchetto clusterProfiler per indagare le funzioni biologiche di ciascun gene chiave. Per ogni gene chiave, i campioni sono stati suddivisi in gruppi ad alta e bassa espressione in base al valore mediano. L'analisi dell'espressione differenziale tra i due gruppi è stata effettuata utilizzando limma, e l'elenco di geni risultante è stato ordinato in base al log₂ fold-change con segno. La GSEA è stata condotta utilizzando la funzione gseGO per i termini del processo biologico GO e la funzione gseKEGG per i percorsi KEGG, con i seguenti parametri: minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0,05 e nPerm = 1.000. L'abbondanza relativa di 28 tipi di cellule immunitarie è stata stimata mediante analisi di arricchimento dei set di geni su singolo campione (ssGSEA) implementata nel pacchetto GSVA. Una matrice di firma di set di geni curata, comprendente geni marcatore per 28 tipi di cellule immunitarie, è stata ottenuta dalla letteratura precedente19. Per ciascun campione, è stata applicata la funzione gsva con method = "ssgsea", ssgsea.norm = TRUE e kcdf = "Gaussian". I coefficienti di correlazione di Spearman tra i punteggi di arricchimento ssGSEA e i livelli di espressione dei sei geni chiave sono stati calcolati utilizzando la funzione cor.test. I valori p sono stati corretti per i test multipli mediante il metodo di Benjamini-Hochberg. La matrice di correlazione è stata visualizzata come mappa termica (heatmap) utilizzando il pacchetto pheatmap.
Previsione di farmaci, docking molecolare e analisi dell'associazione con malattie
I composti terapeutici potenziali che bersagliano geni chiave sono stati identificati utilizzando il database DrugBank. È stata costruita in Cytoscape una rete di interazioni tra "farmaci che bersagliano geni chiave" per visualizzare le interazioni previste tra farmaci e geni. È stato eseguito il docking molecolare utilizzando la piattaforma CB-Dock2 per valutare le affinità di legame. La struttura proteica tridimensionale dell'UPP1 umano è stata ottenuta dal Protein Data Bank (PDB ID: 7B8T). Le strutture molecolari dei farmaci (formato SMILES) sono state recuperate da PubChem. Il docking è stato effettuato utilizzando il motore AutoDock Vina e i risultati sono stati ordinati in base all'energia libera di legame (ΔG, in kcal/mol). I complessi di docking sono stati visualizzati utilizzando PyMOL. Le associazioni tra geni chiave e malattie umane correlate a esposizioni ambientali sono state analizzate utilizzando il Comparative Toxicogenomics Database (CTD). Ogni gene è stato interrogato singolarmente ed è stato estratto il gruppo delle dieci malattie con associazione più forte, i cui dati sono stati rappresentati mediante grafici radar.
Protocollo di RT-qPCR
Sono stati raccolti campioni di sangue venoso periferico da otto pazienti affetti da BPCO e da otto controlli sani presso l'Ospedale Shenzhen Luohu di Medicina Tradizionale Cinese. La diagnosi di BPCO è stata stabilita in base ai criteri della Global Initiative for Chronic Obstructive Lung Disease (GOLD), definita come un rapporto VEMS/FVC post-broncodilatatore < 0,70. Il gruppo di controllo era composto da volontari sani abbinati per età e sesso, privi di anamnesi di malattie respiratorie e con test di funzionalità polmonare normali (VEMS% previsto ≥ 80% e VEMS/FVC ≥ 0,70). Le informazioni di base dei pazienti sono riportate nella Tabella 2. L'RNA totale è stato estratto dai campioni ematici dei pazienti con BPCO utilizzando un kit per l'estrazione di RNA da sangue. Per la sintesi di cDNA, 500 ng di RNA totale sono stati retrotrascritti utilizzando un kit per la sintesi di cDNA con rimozione del DNA genomico, seguendo il protocollo fornito. Il cDNA ottenuto è stato diluito a 150 ng/μL.
La RT-qPCR è stata eseguita utilizzando una miscela maestra per qPCR a base di SYBR Green su un sistema di PCR in tempo reale. Ogni reazione da 10 μL conteneva 5 μL di miscela maestra 2x SYBR Green, 0,5 μL ciascuno di primer diretto e inverso (10 μM), 1 μL di cDNA diluito (15 ng/μL) e 3 μL di acqua priva di nucleasi. Le condizioni di termociclaggio prevedevano una denaturazione iniziale a 95 °C per 5 min, seguita da 40 cicli di 95 °C per 10 s e 60 °C per 30 s, con un'analisi finale della curva di denaturazione da 60 °C a 95 °C per verificare la specificità dell'amplificazione. Tutte le reazioni sono state eseguite in triplicato tecnico. Il gene β-actina è stato utilizzato come gene di riferimento interno. L'efficienza dei primer per ciascun gene bersaglio è stata validata mediante serie di diluizioni per la curva standard ed è risultata compresa tra il 90% e il 110%. I livelli di espressione genica sono stati normalizzati rispetto a β-actina e l'espressione relativa è stata calcolata utilizzando il metodo 2-ΔΔCt. I confronti statistici tra i gruppi COPD e controllo sono stati effettuati utilizzando il test di Mann-Whitney U.
Analisi statistica
Le visualizzazioni della rete sono state create utilizzando Cytoscape e le analisi statistiche sono state eseguite con il software R. Se non indicato diversamente, il test U di Mann-Whitney è stato utilizzato per dati non distribuiti normalmente, mentre il test t di Student è stato utilizzato per dati distribuiti normalmente al fine di confrontare due gruppi. Un valore di p < 0,05 è stato considerato statisticamente significativo.