Die Studie wurde gemäß der Deklaration von Helsinki durchgeführt, und das Protokoll wurde am 22. April 2025 vom Ethikkomitee des Anhui Chest Hospital (K2025-007) genehmigt. Von allen an der Studie beteiligten Probanden wurde informierte Einwilligung eingeholt.
Datenextraktion und Normalisierung
Transkriptomprofile und entsprechende klinische Datensätze für LUAD wurden aus den TCGA- und GEO-Kohorten bezogen. Der TCGA-LUAD-Datensatz wurde als Trainingsset festgelegt, wobei GSE72094, GSE31210 und GSE26939 als Kohorten zur externen Validierung dienten (Tabelle 1). Zusätzlich wurden 900 MCRGs aus einer früheren Studie gesammelt12 (Ergänzende Tabelle 1). Transkriptomdaten wurden mithilfe von GENCODE v36 oder den entsprechenden GPL-Plattform-Annotationsdateien annotiert. Die Probe-IDs wurden in Gen-Symbole umgewandelt, doppelte Gene mittels der avereps-Funktion zusammengeführt, und nur protein-kodierende Gene wurden beibehalten, um Gen-Ebene-Expressionsmatrizen zu erstellen. Für das TCGA-LUAD-Trainingsset wurden Gene mit Fragments per kilobase of exon model per million mapped fragments (FPKM) < 1 in mehr als 50 % der Proben herausgefiltert, und die verbleibenden Expressionswerte wurden log2-transformiert (log2[FPKM+1]). Für die GEO-Validierungskohorten wurden rohe Expressionsdaten heruntergeladen, die Probe-IDs mithilfe der jeweiligen Plattform-Annotationsdateien den Gen-Symbolen zugeordnet, und mehrere Proben, die demselben Gen entsprachen, wurden durch Mittelung ihrer Expressionswerte zusammengefasst. Diese Datensätze wurden gegebenenfalls log2-transformiert. Zwischen TCGA und GEO wurde keine korrekte Bereinigung von plattformübergreifenden Batch-Effekten vorgenommen, da wir eine kohortenweise Standardisierungsstrategie anwandten, um eine relative Vergleichbarkeit sicherzustellen. Konkret wurden für Trainings- und Validierungskohorten die Genexpressionswerte mittels des jeweiligen Mittelwerts und der Standardabweichung jedes Datensatzes zentriert und skaliert (z-Wert-Transformation). Die gleichen aus dem Trainingsset abgeleiteten Cox-Regressionskoeffizienten wurden anschließend verwendet, um Risikoscores für alle Kohorten zu berechnen. Um die klinische Anwendbarkeit zu bewahren und eine Überanpassung an irgendeinen Validierungssatz zu vermeiden, wurde der mediane Risikoscore der Trainingskohorte als fester Schwellenwert verwendet, um Patienten in Hoch- und Niedrigrisikogruppen über alle externen Validierungskohorten hinweg einzuteilen. Klinische Informationen, einschließlich Alter, Geschlecht, pathologisches Stadium, Tumor-Node-Metastase-(TNM-)Stadium, histologischer Typ, Überlebenszeit, Überlebensstatus und Gewebetyp, wurden bei Verfügbarkeit extrahiert. Der Endpunkt war das Gesamtüberleben (OS). Proben mit unvollständigen Überlebensinformationen oder einer Überlebenszeit < 30 Tagen wurden ausgeschlossen. Die Überlebenszeit wurde in Jahre umgerechnet, und der Überlebensstatus wurde als 0 für lebendig und 1 für verstorben kodiert.
Identifizierung und funktionelle Analyse von Kandidatengenen
Das Limma-Paket identifizierte differentiell exprimierte Gene (DEGs) zwischen LUAD-Tumor- und Normalproben im Trainingsset13. Die folgenden Kriterien definierten die DEGs: |log2FC| > 0,5 und adjustierter p-Wert < 0,05. Anschließend wurde der mfuzz-Fuzzy-Clustering-Algorithmus im R-Paket ClusterGVis verwendet, um die DEGs in unterschiedliche Expressionscluster zu unterteilen. Die Genontologie–biologischer Prozess (GO-BP)-Analyse wurde an den fünf repräsentativsten Genen jedes Clusters basierend auf ihren Zugehörigkeitsscores durchgeführt. Eine Menge gemeinsam exprimierter Gene ergab sich aus dem Durchschnitt der DEGs mit den MCRGs. Die funktionelle Anreicherungsanalyse mittels Genontologie/Kyoto Encyclopedia of Genes and Genomes (GO/KEGG) bewertete die biologische Relevanz der überlappenden Gene. Protein-Protein-Interaktionsnetzwerke (PPI) stammten aus der STRING-Datenbank14. Nur Interaktionen mit Konfidenzwerten > 0,7 blieben erhalten, um die Netzwerkzuverlässigkeit zu verbessern.
Screening prognostischer Gene
Mit dem Survival-Paket wurde eine univariate Cox-Regressionanalyse durchgeführt, um mögliche Gene zu identifizieren, die mit dem Gesamtüberleben bei LUAD15 assoziiert sind. Gene mit einem p-Wert < 0,05 wurden als potenzielle prognostische Indikatoren betrachtet. Die TCGA-LUAD-Trainingskohorte umfasste 500 Patienten mit vollständigen Überlebensdaten, von denen 216 (43,2 %) während der Nachbeobachtungszeit ein Todesereignis erlitten. Das Verhältnis von Kandidatengenen (n = 108) zu Ereignissen (n = 216) betrug etwa 1:2, was für eine Cox-Regressionanalyse akzeptabel ist. Anschließend erfolgte eine weitere Merkmalsauswahl mittels Least-Absolute-Shrinkage-and-Selection-Operator-(LASSO-)Regressionsanalyse und eines Extreme-Gradient-Boosting-(XGBoost-)Modells. Die Cox-Proportional-Hazards-Modelle wurden mit family = „cox“ mithilfe der cv.glmnet-Funktion des glmnet-Pakets erstellt. Der optimale Regularisierungsparameter wurde mittels 10-facher Kreuzvalidierung bestimmt, wobei der λ.min-Wert den minimalen Kreuzvalidierungsfehler repräsentierte und als optimaler λ-Wert ausgewählt wurde. Gene mit nicht-null Regressionskoeffizienten wurden als Kandidatenmerkmale extrahiert. Für das XGBoost-Modell wurden Überlebenszeit und Überlebensstatus zur Zielvariablen kombiniert, wobei positive Werte den Todesereignissen und negative Werte den zensierten Fällen zugeordnet wurden. Die Parameter wurden wie folgt festgelegt: objective = „survival:cox“ und eval_metric = „cox-nloglik“, mit 100 Iterationen und einer Lernrate von 0,1. Nach dem Modelltraining wurden die Wichtigkeitsscores der Gene anhand der Feature-Gain-Werte berechnet. Die 20 am höchsten bewerteten Gene wurden nach absteigender Sortierung der Wichtigkeitsscores beibehalten, um die Merkmalsdimensionalität und die Modellkomplexität zu reduzieren. Gemeinsame Gene aus den Ergebnissen von LASSO und XGBoost wurden als prognostische Kandidatengene identifiziert.
Konstruktion und Bewertung eines prognostischen Modells
Ein prognostisches Modell wurde mithilfe einer multivariaten Cox-Regressionanalyse der identifizierten Kandidatengene entwickelt. Die Risikoscores wurden individuell wie folgt berechnet:
.
Dabei bezeichnet Coefi den Koeffizienten für das Gen i und Expi den entsprechenden Genexpressionswert. Die Individuen wurden anschließend anhand des medianen Risikoscores als Trennwert in zwei Gruppen eingeteilt: hohe Risikogruppe und niedrige Risikogruppe. Danach wurden zeitabhängige Receiver-Operating-Characteristic-(ROC)-Kurven erstellt. Um das Potenzial für eine Überanpassung zu bewerten, wurde eine interne Bootstrap-Validierung mit 1.000 Resampling-Durchgängen durchgeführt, um den biaskorrigierten C-Index und die zeitabhängigen AUCs mit 95 %-Konfidenzintervallen zu berechnen. Kalibrierungskurven wurden generiert, um die Übereinstimmung zwischen den vorhergesagten und beobachteten Überlebenswahrscheinlichkeiten nach 2, 3 und 5 Jahren zu bewerten. Außerdem wurde eine Entscheidungskurvenanalyse (DCA) mit dem ggDCA-Paket in R durchgeführt, um den klinischen Nettovorteil des Modells zu den Zeitpunkten 2, 3 und 5 Jahre zu bewerten und den potenziellen Nutzen des Risikoscores in klinischen Entscheidungsprozessen über verschiedene Schwellenwahrscheinlichkeiten hinweg zu quantifizieren. Überlebensunterschiede zwischen den risikostratifizierten Gruppen sowie zwischen anderen klinischen Kategorien wurden mittels Kaplan-Meier-(KM)-Überlebenskurven mit Log-Rank-Test verglichen. Um die Beiträge einzelner Gene zur Modellleistung aufzuklären, wurde eine Shapley-Additive-exPlanations-(SHAP)-Analyse zur nachträglichen erklärbaren Interpretation herangezogen.
Entwicklung und externe Validierung eines Nomogramms
Die Zusammenhänge zwischen den berechneten Risikoscores und verschiedenen klinischen Merkmalen (einschließlich Geschlecht, Alter und TNM-Stadium) wurden mithilfe des Wilcoxon-Rangsummentests oder des Kruskal-Wallis-Tests untersucht, um die klinische Anwendbarkeit des Modells zu bewerten. Um zu prüfen, ob der Risikoscore als unabhängiger prognostischer Faktor fungierte, wurden klinische Variablen zusammen mit dem Risikoscore in eine multivariate Cox-Regression eingebunden. Anschließend wurde ein prognostischer Nomogramm, das die unabhängigen klinischen Risikofaktoren (z. B. Stadium) und den genetischen Risikoscore vereint, mithilfe des regplot-R-Pakets erstellt, um die Vorhersage der Überlebenswahrscheinlichkeit zu individualisieren. Kalibrierungskurven dienten der Beurteilung der Übereinstimmung zwischen den vom Nomogramm vorhergesagten und den tatsächlich beobachteten Überlebensraten. Schließlich wurde die endgültige Vorhersagekraft und Generalisierbarkeit des integrierten Nomogramm-Systems durch zeitabhängige ROC-Kurven und umfassende KM-Analysen klinischer Subgruppen in den Kohorten streng validiert.
Analysen der Immuninfiltration und der Immunsubtypen
CIBERSORT wurde unter Verwendung der Leukozyten-Gen-Signaturmatrix (LM22) verwendet, um die relativen Anteile von 22 Immunzelltypen zur Beurteilung der Immunzellinfiltration bei LUAD-Patienten zu schätzen. Die Beziehungen zwischen prognostischen Genexpressionsniveaus und immunologischer Infiltration wurden mittels Spearman-Korrelationsanalyse bewertet. Immune-, Stromal-, Tumorreinheits- und ESTIMATE-Scores wurden über den ESTIMATE-Algorithmus abgeleitet, und Unterschiede zwischen Risikogruppen wurden mittels Wilcoxon-Test analysiert. LUAD-Patienten wurden mithilfe des ImmuneSubtypeClassifier-Pakets16 sechs Immunsubtypen zugeordnet. Der Wilcoxon-Test wurde zusätzlich verwendet, um die Verteilung der Immunsubtypen zwischen Risikogruppen zu vergleichen.
Analysen der Immun-Checkpoint, des Immunphänotypscores und des Krebsimmunzyklus
In dieser Studie diente der Wilcoxon-Rangsummentest zur Bewertung von 21 Immun-Checkpoint-Genen17 in risikostratifizierten Gruppen, um die immunologische Landschaft von LUAD zu charakterisieren. Die Spearman-Korrelation verknüpfte vielversprechende prognostische Gene mit Immun-Checkpoint-Genen. Um Unterschiede in der Reaktion auf Immun-Checkpoint-Inhibitoren (ICIs) bei LUAD-Patienten mit unterschiedlichem Risiko zu bewerten, wurden Immunphänotypisierungsscores (IPS) für Anti-PD-1- und Anti-CTLA-4-Behandlungen aus dem Cancer Immunome Atlas (TCIA)18 bezogen, und die Datenbank Tracking Tumor Immunophenotype (TIP)19 wurde verwendet, um die Aktivität des Krebs-Immunzyklus zu evaluieren und die entsprechenden Scores zwischen den Risikogruppen zu vergleichen.
Analyse von somatischen Mutationen und Arzneimittelsensitivität
Das TCGA-Mutations-Tool erhielt somatische Mutationsprofile für TCGA-LUAD-Fälle, um die Variation der Mutationsmuster zwischen Risikogruppen zu untersuchen. Das maftools-Paket verarbeitete und visualisierte die Mutationsdaten. Die Tumor-Mutationslast (TMB) wurde für jede Probe bestimmt und zwischen den beiden Risikokategorien verglichen. Eine pharmakogenomische Sensitivitätsanalyse wurde mithilfe des pRRophetic-Pakets gemäß der Genomics of Drug Sensitivity in Cancer (GDSC)-Datenbank20 durchgeführt. Die halbmaximale inhibitorische Konzentration (IC50) von Antikrebsmedikamenten wurde für jeden LUAD-Patienten vorhergesagt, und die Unterschiede zwischen den Risikogruppen wurden mittels des Wilcoxon-Rangsummentests quantifiziert.
Beurteilung der Expressionsniveaus prognostischer Gene
Jeder Datensatz diente der Beurteilung der Transkriptmengen ausgewählter potenzieller, ergebnisbezogener Gene. Um die Genexpression mit der Patientenprognose zu verknüpfen, wurden optimale Cutoff-Werte mithilfe der Funktion surv_cutpoint im R-Paket survminer ermittelt. Basierend auf diesen Schwellenwerten wurden die Fälle von LUAD in Subgruppen mit hoher und niedriger Expression eingeteilt, die anschließend einer Überlebenszeitanalyse unterzogen wurden.
Außerdem wurden fünf gepaarte Proben von Lungenadenokarzinom (LUAD)-Tumoren und angrenzendem Normalgewebe vom Anhui Chest Hospital bezogen, und anschließend erfolgte die Validierung mittels qPCR. Jeder Teilnehmer gab eine schriftliche, informierte Einwilligung ab. Sechs prognostische Kandidatengene (PDGFB, LDHA, ZEB2, FKBP4, DMD und S100B) wurden für die qPCR-Validierung ausgewählt. RNA wurde aus homogenisierten Gewebeproben mittels RNA-Extraktionsreagenz, gefolgt von Chloroform-Extraktion und Isopropanol-Fällung, isoliert. Mithilfe eines Spektrophotometers wurde die RNA-Konzentration und -reinheit gemessen. Die qPCR-Validierung der sechs prognostischen Kandidatengene erfolgte mit einem SYBR-Green-basierten PCR-Mastermix auf einem Echtzeit-PCR-System: zunächst eine Denaturierung bei 95 °C für 30 s, gefolgt von 40 Zyklen mit jeweils 20 s bei 95 °C, 20 s bei 55 °C und 20 s bei 72 °C. Die relative Expression wurde berechnet und mit der 2-ΔΔCt-Methode auf Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) standardisiert. Alle Einzelheiten zu Reagenzien und Instrumenten sind in der Tabelle der Materialien angegeben.
Statistische Analyse
Die statistischen Analysen wurden mit Software für statistisches Rechnen und Visualisierung durchgeführt. Das Protein-Protein-Interaktionsnetzwerk wurde mithilfe von Software zur Netzwerkanalyse dargestellt. Nach der Beurteilung der Normalverteilung wurde der Student-t-Test für normalverteilte kontinuierliche Variablen und der Mann-Whitney-U-Test für nicht normalverteilte Variablen verwendet.