Identification of prognostic biomarkers and immunotherapy response predictors in lung adenocarcinoma: integrative Mendelian randomization and machine learning analysis

Introduction

In 2022, lung cancer was the most commonly diagnosed cancer, with around 2.5 million new cases and an estimated 1.8 million deaths, making it the leading cause of cancer-related mortality globally (1). Lung cancer is generally categorized into non-small cell lung cancer (NSCLC), which accounts for approximately 90% of cases, and small-cell lung cancer (SCLC) (2). Over the past 10 years, the identification of predictive biomarkers has opened up new avenues for treatment, particularly through targeted therapies and immunotherapies (3,4). Molecular targeted therapy has significantly improved outcomes for patients with lung adenocarcinoma (LUAD) harboring oncogenic driver alterations. Recent advances have expanded the therapeutic arsenal beyond well-established targets, such as epidermal growth factor receptor (EGFR), anaplastic lymphoma kinase (ALK), and ROS proto-oncogene 1 (ROS1), to include KRAS G12C mutation, MET exon 14 skipping, RET rearrangements, BRAF V600E, and HER2 (ERBB2) mutations (5,6). In addition, immunotherapy is a relatively novel treatment approach applied in patients with NSCLC and negative for oncogenic drivers. The checkpoint inhibitors used in treating advanced NSCLC include the anti-programmed cell death 1 (PD-1) monoclonal antibodies (mAbs) pembrolizumab and nivolumab and the anti-PD-1 ligand 1 (PD-L1) mAb atezolizumab (7). Despite the considerable efforts exerted to achieve safe and effective therapy, the exposure of patients to unnecessary toxicity, distant metastases, local recurrence, and resistance remain key challenges. As the therapeutic responses and outcomes of patients with LUAD are highly variable, developing reliable predictors to evaluate clinical benefit remains a critical need.

Genome-wide association studies (GWASs) have revealed various genetic variants linked to the occupation and development of adenocarcinoma (8,9). Expression quantitative trait loci (eQTL) are loci that account for a portion of the genetic variation in gene expression phenotypes (10). Summary data-based Mendelian randomization (SMR), with cis-eQTL genetic variants serving as instrumental variable, can assess the association between the exposure (gene expression) and the outcome (phenotype) (11).

Machine learning can recognize relationships between molecular information and clinical symptoms from complex datasets and can also effectively predict the prognosis of patients with cancer (12). It has been extensively implemented in cancer research to predict patients’ clinical outcomes (13). In our study, we constructed a prognostic model for patients with LUAD using a machine learning-based integration program, with The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases being used to validate risk stratification ability. We additionally identified high- and low-risk groups based on the expression levels of identified genes. Furthermore, to assess the ability of the model to predict the prognosis of patients with LUAD, we examined several indicators, including immune checkpoints, tumor microenvironment (TME) score, tumor mutation burden (TMB) score, and abundance of tumor-infiltrating immune cells. We present this article in accordance with the TRIPOD reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2026-0318/rc).

Methods

Data collection

GWAS summary data for LUAD, comprising 1,590 cases and 314,193 controls, were obtained from the FinnGen study (https://www.finngen.fi/en) (14), and lung-related cis-eQTL data were obtained from the Genotype-Tissue Expression (GTEx) project (15). Single-cell RNA-sequencing (RNA-seq) data and spatial transcriptomic data were collected from the GEO database. The single-cell RNA-seq dataset was GSE189357, and the spatial transcriptomic dataset was GSE189487 (16). We used a total of six samples for analysis, with two cases each of adenocarcinoma in situ (AIS), minimally invasive adenocarcinoma (MIA), and invasive adenocarcinoma (IAC).

SMR analysis

To identify genes with a significant causal relationship with LUAD, we employed SMR analysis using SMR software v. 1.3.1 (17). Within the SMR framework, we used cis-eQTL genetic variants as instrumental variables for gene expression. We integrated the effect estimates through an inverse variance-weighted Mendelian Randomization approach. To ensure the observed gene expression-outcome associations were not due to linkage scenarios, we applied the heterogeneity in dependent instruments (HEIDI) test. For each probe, single-nucleotide polymorphisms (SNPs) with a linkage disequilibrium (LD) r² value between >0.90 or <0.05 were excluded. Additionally, we removed one SNP from each pair with an LD r² value >0.90 among the remaining SNPs. In the SMR analysis, a P value <0.05 was considered significant, while a P value >0.05 in the HEIDI test indicated that the observed causal relationship between the exposure and the outcome was not confounded by linkage disequilibrium.

Acquisition of LUAD patient datasets

We sourced RNA-seq datasets along with corresponding clinical information from the TCGA database. RNA-seq data were normalized to transcripts per kilobase million (TPM) values. To ensure the dataset’s integrity and reliability, we excluded cases with incomplete or missing data. Specifically, cases lacking overall survival (OS) data or with OS values under 30 days were removed to enhance the prognostic model’s accuracy. This yielded 493 cases with both gene expression and survival information. For external validation, we downloaded the GSE68465 dataset from the GEO database, which was sequenced with the HumanMethylation450 BeadChip (Illumina, San Diego, CA, USA) (18). After excluding cases with missing OS values or OS values below 30 days, we obtained 439 LUAD cases. To ensure seamless integration of the 316 SMR-identified genes with clinical cohorts, we extracted the corresponding expression profiles by matching gene symbols across the TCGA and GSE68465 datasets. To eliminate non-biological technical biases across the RNA-seq and Microarray platforms, we concatenated the expression matrices and implemented the ComBat function from the sva package. This empirical Bayes framework successfully adjusted for known batch effects, outputting a harmonized expression matrix for downstream machine learning training (19). The TMB, defined as the total count of nonsynonymous mutations across a patient’s exomes per million bases, was calculated.

Construction of a prognostic model via machine learning

After the SMR analysis, we applied univariate Cox regression to identify the genes associated with prognosis, with the significance level set at P<0.05. A prognostic model for patients with LUAD was then constructed through use of 10 distinct machine learning algorithms, including least absolute shrinkage and selection operator, ridge regression, elastic net, stepwise Cox regression, survival-support vector machine, Cox likelihood-based boosting, supervised principal component analysis (PCA), partial least squares regression for Cox, random survival forest (RSF), and generalized boosted regression modeling. These algorithms were evaluated via a leave-one-out cross-validation (LOOCV) approach across 101 model combinations based on the TCGA dataset and were further validated on the GEO dataset. Models with fewer than five genes were excluded. The concordance index (C-index) was computed for each model across all datasets. A risk score was assigned to each sample based on the expression levels of genes in the prognostic model. Through a method similar to that used in a previous study (20), the samples were subsequently categorized into low- and high-risk groups, with the median risk score serving as the threshold.

Validation of the prognostic model

Survival analysis was performed to validate the association between OS and the prognostic model. The model’s prognostic value was assessed through time-dependent receiver operating characteristic (ROC) curves and the corresponding area under the curve (AUC) via the survminer and timeROC R packages. To assess whether the risk score and clinical features were independent prognostic factors, we conducted both univariate and multivariate Cox regression analyses.

Establishment of the nomogram

A nomogram was developed to predict the 1-, 3-, and 5-year survival probabilities for patients with LUAD via the regplot and rms R packages. The nomogram incorporated factors such as risk score, age, sex, and tumor stage. Additionally, we evaluated the performance of the nomogram by creating a calibration curve based on the Hosmer-Lemeshow test.

Evaluation of the TME

We employed the Estimation of Stromal and Immune Cells in Malignant Tumor Tissues Using Expression data (ESTIMATE) algorithm to assess the components of the TME. This algorithm was employed to assess tumor purity, along with stromal, immune, and ESTIMATE scores. We examined antitumor immune responses and calculated pathway scores for all samples using single-sample gene set enrichment analysis (GSEA). Furthermore, we compared the expression of immune checkpoint genes between the low- and high-risk groups to evaluate the ability of the risk score to predict responses to immune checkpoint inhibitor (ICI) therapy.

GSEA and sensitivity to chemotherapy

We used GSEA to compare enriched hallmark functions and pathways between the low- and high-risk groups. Drug sensitivity data for human cancer cell lines were sourced from the Cancer Therapeutics Response Portal and the Profiling Relative Inhibition Simultaneously in Mixtures databases. The pRRophetic R package was applied to estimate the half-maximal inhibitory concentration (IC₅₀) of frequently used chemotherapy agents.

Single-cell transcriptome and spatial transcriptome analyses

The Seurat R package was used to analyze and process single-cell sequencing data and spatial transcriptomics data. For the scRNA-seq datasets, rigorous quality control was applied to screen out low-quality cells. We set the thresholds as follows: cells with >20% mitochondrial gene expression, >1% red blood cell gene expression or feature counts outside the range of 300 to 20,000 were excluded from the analysis. After data normalization, the top 3,000 highly variable genes (HVGs) for each dataset were identified using the “vst” selection method. Following data scaling, we performed dimensionality reduction using PCA. The first 18 principal components (PCs) were retained for neighbor graph construction (FindNeighbors), and similar cell groups were clustered using the Louvain algorithm with a resolution set to 0.5 (FindClusters). Non-linear dimensionality reduction was performed via uniform manifold approximation and projection (UMAP) for visualization. To identify differentially expressed genes (DEGs) defining each cluster, the FindAllMarkers function, based on the Wilcoxon rank-sum test, was employed with criteria set as an adjusted P value <0.05 and a log2(fold change) >0.25. By reviewing the literature, we annotated the cell clusters by comparing the DEGs significantly expressed in each cell group with commonly known marker genes of various cell types (21-23).

For the spatial transcriptomics data analysis, UMI counts and corresponding images were loaded using the Read10X_Image function. Normalization, scaling, and highly variable feature identification were concurrently processed using the SCTransform function to mitigate technical variations while preserving biological heterogeneity. For region division within the spatial context, we applied unbiased unsupervised clustering. PCA was conducted on the SCT-normalized data, and the first 30 PCs were selected to identify spatial neighborhoods and construct cell clusters, which were subsequently projected onto UMAP and spatial coordinates (SpatialDimPlot). Furthermore, the analytical logic for determining gene spatial expression relied on identifying spatially variable features. We utilized the FindSpatiallyVariableFeatures function, applying the Moran’s I spatial autocorrelation method, to quantitatively select top genes exhibiting significant spatial expression patterns across the tissue slices.

Statistical analysis

All statistical analyses and data visualizations were performed using R software (version 4.1.2) and SMR software (version 1.3.1). For the comparison of continuous variables between two groups, the Student’s t-test or the Wilcoxon rank-sum test was employed, depending on the normality of the data distribution. The Kaplan-Meier method, coupled with the log-rank test, was utilized to compare OS between the low- and high-risk groups. To determine the independent prognostic value of the risk score and other clinical characteristics, univariate and multivariate Cox proportional hazards regression analyses were conducted to estimate hazard ratios (HRs) and 95% confidence intervals (CIs). The predictive accuracy of the prognostic model was evaluated using time-dependent ROC curves and the AUC. For single-cell RNA-seq data analysis, the FindMarkers function based on the Wilcoxon rank-sum test was used, with significance defined as an adjusted P value <0.01 and an absolute log2(fold change) >0.25. In the SMR analysis, a P value <0.05 was considered statistically significant, while a P value >0.05 in the HEIDI test was used to rule out confounding by linkage disequilibrium. Unless otherwise specified, a two-sided P value <0.05 was considered to indicate statistical significance.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Results

SMR analysis

GWAS summary data of 20,191,338 SNPs were included from the finngen_R10_C3_NSCLC_ADENO dataset, and 914,625 SNPs were included from the GTEx_V8_cis_eqtl_summary_Lung.lite dataset. Following the evaluation of allele consistency for each SNP across paired datasets, a total of 763,748 SNPs were retained. Subsequently, 335 SNPs exhibiting allele frequency differences greater than 0.20 between any pair of datasets were excluded from further analysis. We conducted SMR analysis and HEIDI tests for 7027 probes to examine the relationships between gene expression and LUAD with genetic instruments in the cis region. Ultimately, 316 genes were selected based on a P value for SMR of less than 0.05 and a P value for HEIDI of greater than 0.05 (Figure 1A).

Figure 1 SMR analysis in lung adenocarcinoma and prognosis related genes. (A) Manhattan plots in SMR analysis. Red plots indicate the genes from the SMR findings with a P value <0.05, a false-discovery rate <0.05, and P_HEIDI value >0.05. (B) Univariate Cox regression analysis of 316 genes with causal associations. Chr, chromosome; HEIDI, heterogeneity in dependent instruments; SMR, summary data-based Mendelian randomization.

Identification of prognosis-related genes

Univariate Cox regression analysis was performed on 316 candidate gene expression values, and the HR and P value of each gene were computed. Among the genes, 20 were significantly associated (P<0.05) with poor prognosis and 13 with good prognosis (Figure 1B).

Machine learning

To establish a reliable and precise prognostic model, we employed a machine learning-based integration approach incorporating the 33 identified genes. In this process, we trained 101 predictive models within the LOOCV framework and evaluated the C-index (table available at https://cdn.amegroups.cn/static/public/tlcr-2026-0318-1.xlsx). The best-performing model, which integrated RSF and ridge regression, achieved the highest average C-index of 0.652 (Figure 2A). Gene signatures in this model included CERS4, IL22RA1, OPN3, PPIA, CYP4B1, PHF11, GALNT3, RNASET2, PYGB, SULT1A1, ATP1B2, YWHAE, PNP, CDCA3, ARHGEF4, RGS11, POLR3G, ARF3, RPA3, PCMTD2, FAIM2, EXOSC5, ZNF322, GOT2, SPATA5L1, and CLDN18 (table available at https://cdn.amegroups.cn/static/public/tlcr-2026-0318-1.xlsx). We then calculated the risk score from the GEO and TCGA databases and divided them into high- and low-risk groups. Patients with LUAD in the high-risk group had poor survival duration according to the prognostic model (Figure 2B,2C). The characteristics of the patients are shown in table available at https://cdn.amegroups.cn/static/public/tlcr-2026-0318-1.xlsx.

Figure 2 Machine learning-based prognostic model. (A) A total of 101 kinds of machine learning algorithms were used to obtain the optimal model, and the C-index of each model was calculated for all cohorts. (B) Kaplan-Meier curves of overall survival according to the risk score in TCGA-LUAD. (C) Kaplan-Meier curves of overall survival according to the risk score in GSE68465. GEO, Gene Expression Omnibus; TCGA-LUAD, The Cancer Genome Atlas-lung adenocarcinoma.

Evaluation of model performance

We constructed ROC curves incorporating risk score, age, gender, and stage and estimated their reliability by computing the AUC value. We found that risk score had a maximum AUC value of 0.757 (Figure 3A) Time-dependent ROC analysis indicated that the AUCs for predicting OS at 1, 3, and 5 years were 0.757, 0.707, and 0.686, respectively (Figure 3B). Furthermore, we conducted both univariate and multivariate Cox regression analyses to determine the prognostic value of age, gender, stage, and risk score. The analysis revealed that both stage (P<0.001) and risk score (P<0.001) were significantly correlated with prognosis and acted as independent risk factors both in the univariate (Figure 3C) and multivariate Cox regression analyses (Figure 3D). These results validated that our risk score-based model is a reliable and stable independent predictor of prognosis in LUAD patients.

Figure 3 Evaluation of the performance of model. (A) The 1-year ROC curves of the risk score and clinical characteristics. (B) The 1-, 3-, and 5-year ROC curves of the risk score. (C) Univariate Cox regression analysis of clinical characteristics and risk score in patients with lung adenocarcinoma. (D) Multivariate Cox regression analysis of clinical characteristics and risk score in patients with lung adenocarcinoma. AUC, area under the curve; CI, confidence interval; ROC, receiver operating characteristic.

Nomogram construction

Based on the risk score, age, gender, and grade, we developed a nomogram to predict the 1-, 3-, and 5-year OS rates for patients with LUAD (Figure 4A). To assess the nomogram’s accuracy, we generated a calibration plot (Figure 4B), and the results indicated that the nomogram effectively predicted patient risk with high accuracy.

Figure 4 Construction of the nomogram and GSEA. (A) Performance of the nomogram that included the risk score and tumor stage in predicting the probability of 1-, 3-, and 5-year overall survival. (B) Calibration curves for the nomogram. (C) Enriched pathways of the high-risk group according to GSEA. ***, P<0.001. CI, confidence interval; GSEA, gene set enrichment analysis; OS, overall survival.

GSEA of the low- and high-risk groups

We conducted GSEA to identify the potential key pathways in the high-risk group in order clarify the specific treatment response and clinical outcome of high-risk score group. GSEA showed that high-risk group was enriched to 10 pathways: Alzheimer disease, cell cycle, DNA replication, Huntington disease, oocyte meiosis, oxidative phosphorylation, Parkinson disease, pyrimidine metabolism, ribosome, and spliceosome (Figure 4C). The findings indicated that the poor prognosis in the high-risk group could be associated with disruptions in the cell cycle and genomic instability.

The TME and therapy response of the different groups

We analyzed the infiltration levels and functional roles of immune cells within the TME and found significant differences between the low- and high-risk groups for B cells, CD8+ T cells, checkpoint pathways, cytolytic activity, dendritic cells, human leukocyte antigens, mast cells, neutrophils, natural killer cells, T-cell coinhibition, T-cell costimulation, T helper cells, T follicular helper cells, tumor-infiltrating lymphocytes, regulatory T cells, and type II IFN response (Figure 5A). The analysis of immune checkpoint molecules revealed distinct expression profiles: the low-risk group had elevated expression of BTLA, CD40LG, TNFRSF8, CD200, CD200R1, CD244, CD48, CD160, TNFRSF14, NRP1, TNFRSF25, CD28, TNFSF14, CD80, LAIR1, TNFSF18, BTNL2 and IDO2, whereas the high-risk group had elevated expression of TNFSF9 and CD276 (Figure 5B).

Figure 5 Comparison of the TME and TMB between the high- and low-risk groups. (A) Difference in the immune pathways of the two groups. (B) Difference in the immune checkpoints between the two groups. (C) TME score for the two groups. (D) TMB for the two groups. (E) TMB was positively correlated with the risk score. *, P<0.05; **, P<0.01; ***, P<0.001. TMB, tumor mutation burden; TME, tumor microenvironment.

TME score has emerged as a valuable prognostic indicator, enhancing predictive accuracy beyond conventional clinicopathological metrics (24). Notably, in our study, patients with LUAD classified as low risk exhibited increased ESTIMATE, immune, and stromal scores relative to their high-risk counterparts, implying that elevated TME-related indices correlate inversely with the risk score (Figure 5C). Furthermore, in a number of malignancies, a higher TMB has been linked to better outcomes in individuals undergoing ICI therapy (25). Compared to the low-risk group, the high-risk group exhibited elevated TMB levels (Figure 5D). Moreover, TMB was positively linked to the risk score (Figure 5E), suggesting that patients with LUAD and higher risk scores may receive diminished benefit from immunotherapy.

We assessed the responsiveness of patient cohorts to commonly administered chemotherapeutic agents by determining their IC₅₀. Notably, the high-risk subgroup exhibited enhanced sensitivity to 5-fluorouracil, alisertib, alpelisib, cediranib, cytarabine, crizotinib, dasatinib, erlotinib, foretinib, and gallibiscoquinazole (Figure 6). Meanwhile, the low-risk subgroup was more sensitive to axitinib and doramapimod (Figure 6).

Figure 6 Sensitivity to chemotherapy of the low- and high-risk groups. The half-maximal inhibitory concentration was used to evaluate the sensitivity.

Single-cell transcriptome analysis

After the initial screening, we obtained 30,798 cells from the AIS group, 17,981 cells from the MIA group, and 28,212 cells from the IAC group for subsequent analysis. After identifying 2,000 HVGs and performing dimensionality reduction analysis, we divided the cells in the AIS group into 20 cell clusters, the cells in the MIA group into 21 cell clusters, and the cells in the IAC group into 19 cell clusters (Figure S1). By reviewing relevant marker genes, we annotated these cell clusters, and the distribution of cells in each group and the expression of key marker genes are shown in Figure 7. We further examined the 26 genes used to construct the prognostic model in terms of their expression in different cell types and pathological tissue types. As shown in Figure 8, CERS4, OPN3, PPIA, CYP4B1, PHF11, GALNT3, RNASET2, PYGB, SULT1A1, YWHAE, PNP, ARF3, RPA3, PCMTD2, ZNF322 and CLDN18 exhibited higher expression levels in the AIS, the MIA and the IAC groups. Among them, CLDN18 was significantly expressed only in epithelial cells, and its expression decreased as the disease progressed from AIS to MIA and then to IAC. CERS4 was found to be mainly expressed in ciliated cells, with its expression in epithelial cells decreasing with disease progression. CYP4B1 was primarily found to be expressed in both ciliated cells and epithelial cells, and similarly, its expression in epithelial cells decreased with disease progression. In the univariate Cox regression analysis of OS, the expression of CLDN18, CERS4, and CYP4B1 was associated with low-risk factors. The results of the single-cell transcriptome analysis corroborated those from the univariate Cox regression analysis. Additionally, regardless of disease progression, PPIA and YWHAE were highly expressed across multiple cell types.

Figure 7 Cell annotation and distribution. (A) Overview of cell distribution in the adenocarcinoma in situ group. (B) Overview of cell distribution in the minimally invasive adenocarcinoma group. (C) Overview of cell distribution in the invasive adenocarcinoma group. (D) Bubble plot for the expression of marker genes in the adenocarcinoma in situ group. (E) Bubble plot for the expression of marker genes in the minimally invasive adenocarcinoma group. (F) Bubble plot for the expression of marker genes in the invasive adenocarcinoma group.

Figure 8 Bubble plot for the expression of 26 genes in the prognostic model. (A-C) Results from the invasive adenocarcinoma group, minimally invasive adenocarcinoma group, and adenocarcinoma in situ group, respectively.

Spatial transcriptome analysis

Using the Seurat R package, we processed the spatial transcriptomics data corresponding to the six single-cell samples (Figure S2). After normalization, scaling, and data correction, we examined the spatial expression of 26 genes in tumor tissues for constructing the prognostic model. The expression data of each gene in the six samples can be found in Figures S3-S8. Interestingly, we found that CLDN18 was primarily expressed in the cells at the periphery of the tumor mass (Figure 9A). Regardless of whether it developed into invasive adenocarcinoma, PPIA was highly expressed in every cell type and at every spatial location within the tumor (Figure 9B). Finally, the expression of CERS4 and ARF3 was more evenly distributed, decreasing with disease progression to IAC (Figure 9C,9D).

Figure 9 Expression of genes in the spatial dimension of tumor tissue. (A-D) The expression distribution of CLDN18, PPIA, CERS4, and ARF3 in six samples, respectively. AIS, adenocarcinoma in situ; IAC, invasive adenocarcinoma; MIA, minimally invasive adenocarcinoma.

Discussion

Lung cancer is one of the most prevalent malignancies globally, accounting for a substantial proportion of cancer-related deaths worldwide. Despite continuous advancements being achieved in clinical treatment modalities, a substantial portion of patients undergoing therapy still experience distant metastases, local recurrence, and drug resistance. With the advancement of high-throughput molecular detection technologies, it has become possible to perform detailed stratification of patients for more precise and effective treatment by analyzing vast amounts of clinical data derived from patients. Meanwhile, machine learning has been applied to effectively determine the relationships between these data and aid in predicting the prognosis of patients with cancer (26). Given this, we sought to integrate machine learning with extensive clinical data to identify predictive features among different patients with LUAD. In our study, we first conducted SMR analysis to screen out genes with causal effects. Subsequently, we developed a novel computational framework incorporating 10 machine learning algorithms and 101 combinations of algorithm interactions to construct predictive models, aiming to enhance the prognostic classification of patients with LUAD.

First, 316 genes with a P value <0.05 in the SMR analysis and a P value >0.05 in HEIDI were identified. Univariate Cox analysis was employed to further screen these 316 genes, which ultimately identified 33 genes associated with patient prognosis. Among these, 20 genes, including peptidylprolyl isomerase A (PPIA), interleukin 22 receptor subunit alpha 1 (IL-22RA1), opsin 3 (OPN3), and cell division cycle-associated 3 (CDCA3), were determined to be significant risk factors, while 13 genes, including fas apoptotic inhibitory molecule 2 (FAIM2) and ceramide synthase 4 (CERS4), were found to be significant protective factors. PPIA can promote cancer progression in patients with NSCLC with nuclear factor erythroid 2-related factor 2 (NRF2) hyperactivation by modulating NRF2 activity (27). The IL-22-IL-22R1 axis can greatly contribute to cancer progression by inducing oxidative signaling downstream of growth, angiogenesis, metastasis, inflammation, chemoresistance, and signal transducer and activator of transcription 3 (STAT3), a highly procarcinogenic signal transducer (28). High expression of OPN3 promotes epithelial-mesenchymal transformation and cancer metastasis in LUAD cells (29). CDCA3, as part of the SKP1-Cullin-RING-F-box (SCF) ubiquitin ligase (E3) complex, can regulate cell cycle and promote LUAD progression by degrading the endogenous cell cycle inhibitor, WEE1 G2 checkpoint kinase (30). In most tumors, FAIM2 expression is associated with increased CD8⁺ T-cell infiltration and decreased myeloid-derived suppressor cell infiltration (31). CERS4, an important sphingolipid metabolizing enzyme, can enhance anti-PD-1 therapy for patients with LUAD (32). Collectively, these findings can provide novel perspectives for examining the mechanisms of ICI resistance in patients with LUAD.

Subsequently, 101 models were fitted via LOOCV in the machine learning process. We optimized these models and identified the combination of the RSF and ridge regression algorithms as the optimal model according to the highest average C-index. Using this model, we constructed time-dependent and clinically relevant ROC curves. The results demonstrated that this model possesses excellent ability to predict patient prognosis at 1, 3, and 5 years, outperforming conventional clinical indicators such as staging, gender, and age. Analysis of nomogram accuracy also indicated that risk score, age, gender, and grade had favorable predictive value (33).

To investigate the specific molecular mechanisms underlying the prognostic differences across patients, we employed GSEA to identify differentially expressed pathways. We found that the biological processes of cell cycle, DNA replication, oocyte meiosis, and oxidative phosphorylation pathway were enriched pathways. Recent studies indicate that increasing cell cycle activity in cancer cells inhibits antitumor immunity (34). These findings suggested that the poor prognosis of the high-risk group may be related to the abnormality of cell cycle and the instability of hereditary substance, which is closely linked to tumor proliferation and progression (35). Therefore, the inhibition of tumor cell cycle progression through cyclin-dependent kinase (CDK) inhibitors might be a beneficial therapy for this high-risk group. Benson et al. reported that combination treatment with miR-143 and miR-506 downregulated CDK4 and CDK1 levels and induced lymphocyte apoptosis (36).

The TME is critical to inducing cell proliferation, reducing cell death, promoting angiogenesis, facilitating cancer metastasis, increasing inflammatory responses, and counteracting immune responses (37). The TME consists of noncancerous host cells, including adaptive and innate immune cells, and their noncellular components (38). Specifically in LUAD, the crosstalk between cancer cells and stromal components, notably cancer-associated fibroblasts (CAFs), dictates tumor trajectory. CAFs exhibit profound phenotypic heterogeneity and drive disease progression via core regulatory mechanisms, including the secretion of transforming growth factor-beta (TGF-β) and CXCL12. These specific interaction modes facilitate extracellular matrix (ECM) cross-linking, establishing a physical barrier that restricts the infiltration of cytotoxic CD8⁺ T cells, thereby orchestrating an immune-excluded microenvironment (39). Therefore, targeting these core CAF-tumor regulatory nodes presents a robust strategy to overcome resistance to conventional ICIs. In the TME, the failure of immune cells is an important cause of poor prognosis, which suggests that restoring the immune function of tumor-infiltrating lymphocytes by targeting inhibitory receptors can be highly efficacious (40). ICIs inhibit immune cell suppression signaling and activate T cell-mediated cytotoxic responses to exert antitumor effects (41). Clinical oncologists are combining traditional therapies with immunotherapy to improve outcomes. To identify potential factors contributing to the variability in prognoses from the perspective of the TME and to assess the application of our model in clinical immunotherapy, we analyzed indicators associated with immunotherapeutic efficacy. These indicators included immune checkpoints, the infiltration levels of immune cells within the TME, immune functional pathways, TME scores, and TMB scores.

The high-risk group had a lower immune function score, lower expression of most immune checkpoints, and a lower TME score, which suggests high-risk populations benefits less from immunotherapy than do low-risk ones. The high expression of immune checkpoints is closely associated with better response to immunotherapy. As a specific type of ICP, CD276 can reflect the therapeutic effect after the administration of immunotherapy. Silencing CD276 inhibits lung cancer progression by modulating integrin signaling (42). We found a higher expression of CD276 in the high-risk group, which indicates benefits from the silencing of CD276 in suppressing lung cancer progression. In addition, higher TNFSF9 (CD137) expression was found in the high-risk group. Activated CD137 can provide costimulatory signals independent of CD28, enhancing T-cell activity and reducing ICI resistance (43). Geuijen et al. determined that MCLA-145 can act on both CD137 and PD-L1, thereby enhancing the antitumor activity of relevant immune cells, which might benefit high-risk patients (44).

The extent of immune cell infiltration in the TME substantially impacts the clinical benefit derived from ICI therapy in patients with advanced LUAD (45). In our study, the high-risk group generally exhibited lower levels of immune cell infiltration, suggesting that these patients may experience limited benefit from immunotherapy and can be regarded as inflammation-deficient “cold” tumors (46). Combining therapies that augment the priming of T-cell antitumor responses with approaches that remove immunosuppressive signals and provide immunostimulatory signals may be a promising strategy than can allow this patient population to benefit from immunotherapy (47). However, such integration may increase the incidence of side effects, and thus a careful assessment of efficacy and adverse effects is needed (47).

TMB is the amount of mutations per million bases in the coding region of the genome of a malignant cell and reflects the immunogenicity of the tumor. In most tumors, a higher TMB predicts better ICI efficacy (48). The US Food and Drug Administration approved pembrolizumab for patients with a high TMB. With the continuous advancement of related research, the detection methods and cutoff values for TMB are expected to become more standardized. Consequently, an increasing number of drugs and therapeutic approaches are anticipated to be available to patients with tumor whose TMB levels meet the treatment criteria (49). Interestingly, we found that the high-risk group in our study had a higher TMB score.

Chemotherapy is among the preferred treatments for LUAD, and so we further analyzed the IC₅₀ of common drugs in the high- and low-risk LUAD groups. The high-risk group may be more likely to benefit from treatment with 5-fluorouracil, alisertib, alpelisib, cediranib, cytarabine, crizotinib, dasatinib, erlotinib, foretinib, and gallibiscoquinazole. Meanwhile, the low-risk group may be more likely to benefit from treatment with axitinib and doramapimod. Our findings can help to inform individualized treatment regimens for patients and reduce drug resistance.

While previous studies have employed SMR and machine learning pipelines for tumor biomarker screening, our study is distinct in its multilayered omics integration. Conventional approaches often halt at bulk transcriptomic predictions, which fail to capture intra-tumoral heterogeneity. In contrast, our study mapped the macroscopic causal signals derived from large-scale GWAS to a highly localized microscopic context, revealing the spatial margins and specific epithelial origins of genes like CLDN18 and CERS4 along the continuum from AIS to IAC. Clinically, this prognostic model holds substantial translational potential. The identified 26-gene signature could be transitioned into a customized multiplex RT-qPCR or targeted RNA-seq panel utilized on routine surgical resections or biopsy specimens. By calculating the risk score, oncologists can identify high-risk patients who might require intensified adjuvant chemotherapy. More importantly, while existing biomarkers such as TMB and PD-L1 Tumor Proportion Score (TPS) provide baseline guidance for immunotherapy, they are imperfect predictors on their own. Integrating our risk score with these existing metrics could offer a multi-dimensional stratification tool, distinguishing patients who will genuinely benefit from ICIs from those who harbor an immune-excluded phenotype and may require CAF-targeted combinatorial therapies.

Admittedly, our study involved certain limitations that should be addressed. First, as we employed a retrospective design, a prospective cohort study is needed to test the robustness and validity of our method. Second, the sample used was insufficiently large to generalize the predictive model to a larger scale. Third, more accurate clinical validation is required before the molecular histology results can be applied clinically. Regarding the future of modeling in oncology, novel approaches should be developed to overcome the above-mentioned limitations.

Conclusions

We employed bioinformatics analyses in conjunction with multiple machine learning algorithms to construct a theoretical model capable of predicting the prognosis of patients with LUAD. Moreover, we identified low- and high-risk groups in this population, with different prognoses and immune characteristics. A higher risk score was associated with poorer outcomes from immunotherapy. Our prognostic model may be able to optimize decision-making for the individualized treatment of patients with LUAD.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2026-0318/rc

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2026-0318/prf

Funding: This study was supported by the Research Project Foundation of Hubei Provincial Health Commission (No. WJ2023M024).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2026-0318/coif). H.X. reports funding from the Research Project Foundation of Hubei Provincial Health Commission (No. WJ2023M024). The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Zhang Y, Vaccarella S, Morgan E, et al. Global variations in lung cancer incidence by histological subtype in 2020: a population-based study. Lancet Oncol 2023;24:1206-18. [Crossref] [PubMed]
3. Wang M, Herbst RS, Boshoff C. Toward personalized treatment approaches for non-small-cell lung cancer. Nat Med 2021;27:1345-56. [Crossref] [PubMed]
4. Herbst RS, Morgensztern D, Boshoff C. The biology and management of non-small cell lung cancer. Nature 2018;553:446-54. [Crossref] [PubMed]
5. Shin H, Hwang S, Jeong JH, et al. Targeting USP47 enhances the efficacy of KRAS inhibitor in KRAS(G12C) mutated non-small cell lung cancer by controlling deubiquitination of c-Myc. Pharmacol Res 2025;215:107722. [Crossref] [PubMed]
6. Sakamoto T, Matsubara T, Takahama T, et al. Biomarker Testing in Patients With Unresectable Advanced or Recurrent Non-Small Cell Lung Cancer. JAMA Netw Open 2023;6:e2347700. [Crossref] [PubMed]
7. Alexander M, Kim SY, Cheng H. Update 2020: Management of Non-Small Cell Lung Cancer. Lung 2020;198:897-907. [Crossref] [PubMed]
8. Dai J, Lv J, Zhu M, et al. Identification of risk loci and a polygenic risk score for lung cancer: a large-scale prospective cohort study in Chinese populations. Lancet Respir Med 2019;7:881-91. [Crossref] [PubMed]
9. Zhang R, Shen S, Wei Y, et al. A Large-Scale Genome-Wide Gene-Gene Interaction Study of Lung Cancer Susceptibility in Europeans With a Trans-Ethnic Validation in Asians. J Thorac Oncol 2022;17:974-90. [Crossref] [PubMed]
10. Shabalin AA. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 2012;28:1353-8. [Crossref] [PubMed]
11. Zhu Z, Zhang F, Hu H, et al. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nat Genet 2016;48:481-7. [Crossref] [PubMed]
12. Li R, Li L, Xu Y, et al. Machine learning meets omics: applications and perspectives. Brief Bioinform 2022;23:bbab460. [Crossref] [PubMed]
13. Xu Q, Ma L, Streuer A, et al. Machine learning-based in-silico analysis identifies signatures of lysyl oxidases for prognostic and therapeutic response prediction in cancer. Cell Commun Signal 2025;23:169. [Crossref] [PubMed]
14. Kurki MI, Karjalainen J, Palta P, et al. FinnGen provides genetic insights from a well-phenotyped isolated population. Nature 2023;613:508-18. [Crossref] [PubMed]
15. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 2020;369:1318-30.
16. Zhu J, Fan Y, Xiong Y, et al. Delineating the dynamic evolution from preneoplasia to invasive lung adenocarcinoma by integrating single-cell RNA sequencing and spatial transcriptomics. Exp Mol Med 2022;54:2060-76. [Crossref] [PubMed]
17. Wu Y, Zeng J, Zhang F, et al. Integrative analysis of omics summary data reveals putative mechanisms underlying complex traits. Nat Commun 2018;9:918. [Crossref] [PubMed]
18. de Goede OM, Razzaghian HR, Price EM, et al. Nucleated red blood cells impact DNA methylation and expression analyses of cord blood hematopoietic cells. Clin Epigenetics 2015;7:95. [Crossref] [PubMed]
19. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
20. Zhang Y, Wang Y, Chen J, et al. A programmed cell death-related model based on machine learning for predicting prognosis and immunotherapy responses in patients with lung adenocarcinoma. Front Immunol 2023;14:1183230. [Crossref] [PubMed]
21. Chen Z, Huang Y, Hu Z, et al. Landscape and dynamics of single tumor and immune cells in early and advanced-stage lung adenocarcinoma. Clin Transl Med 2021;11:e350. [Crossref] [PubMed]
22. Zhang J, Song C, Tian Y, et al. Single-Cell RNA Sequencing in Lung Cancer: Revealing Phenotype Shaping of Stromal Cells in the Microenvironment. Front Immunol 2021;12:802080. [Crossref] [PubMed]
23. Bu L, Baba H, Yoshida N, et al. Biological heterogeneity and versatility of cancer-associated fibroblasts in the tumor microenvironment. Oncogene 2019;38:4887-901. [Crossref] [PubMed]
24. Zhang H, Song J, Dong J, et al. Tumor Microenvironment Analysis Identified Subtypes Associated With the Prognosis and the Tumor Response to Immunotherapy in Bladder Cancer. Front Genet 2021;12:551605. [Crossref] [PubMed]
25. Samstein RM, Lee CH, Shoushtari AN, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet 2019;51:202-6. [Crossref] [PubMed]
26. Shi Z, Lei JT, Elizarraras JM, et al. Mapping the functional network of human cancer through machine learning and pan-cancer proteogenomics. Nat Cancer 2025;6:205-22. [Crossref] [PubMed]
27. Lu W, Cui J, Wang W, et al. PPIA dictates NRF2 stability to promote lung cancer progression. Nat Commun 2024;15:4703. [Crossref] [PubMed]
28. Lim C, Savan R. The role of the IL-22/IL-22R1 axis in cancer. Cytokine Growth Factor Rev 2014;25:257-71. [Crossref] [PubMed]
29. Zhang W, Feng J, Zeng W, et al. Integrated analysis of the prognostic and oncogenic roles of OPN3 in human cancers. BMC Cancer 2022;22:187. [Crossref] [PubMed]
30. Adams MN, Burgess JT, He Y, et al. Expression of CDCA3 Is a Prognostic Biomarker and Potential Therapeutic Target in Non-Small Cell Lung Cancer. J Thorac Oncol 2017;12:1071-84. [Crossref] [PubMed]
31. Cai J, Ye Z, Hu Y, et al. FAIM2 is a potential pan-cancer biomarker for prognosis and immune infiltration. Front Oncol 2022;12:998336. [Crossref] [PubMed]
32. Wang J, Li RZ, Wang WJ, et al. CERS4 predicts positive anti-PD-1 response and promotes immunomodulation through Rhob-mediated suppression of CD8(+)Tim3(+) exhausted T cells in non-small cell lung cancer. Pharmacol Res 2023;194:106850. [Crossref] [PubMed]
33. Lin JT, Chan TC, Li CF, et al. Downregulation of the cytochrome P450 4B1 protein confers a poor prognostic factor in patients with urothelial carcinomas of upper urinary tracts and urinary bladder. APMIS 2019;127:170-80. [Crossref] [PubMed]
34. Li J, Stanger BZ. Cell Cycle Regulation Meets Tumor Immunosuppression. Trends Immunol 2020;41:859-63. [Crossref] [PubMed]
35. Matthews HK, Bertoli C, de Bruin RAM. Cell cycle control in cancer. Nat Rev Mol Cell Biol 2022;23:74-88. [Crossref] [PubMed]
36. Hossian AKMN, Mackenzie GG, Mattheolabakis G. Combination of miR 143 and miR 506 reduces lung and pancreatic cancer cell growth through the downregulation of cyclin dependent kinases. Oncol Rep 2021;45:2. [Crossref] [PubMed]
37. Xiao Y, Yu D. Tumor microenvironment as a therapeutic target in cancer. Pharmacol Ther 2021;221:107753. [Crossref] [PubMed]
38. Xie Y, Xie J, Gupta P, et al. Editorial: Tumor microenvironment and cancer therapy. Front Cell Dev Biol 2023;11:1290456. [Crossref] [PubMed]
39. Altorki NK, Markowitz GJ, Gao D, et al. The lung microenvironment: an important regulator of tumour growth and metastasis. Nat Rev Cancer 2019;19:9-31. [Crossref] [PubMed]
40. Zhao Y, Shao Q, Peng G. Exhaustion and senescence: two crucial dysfunctional states of T cells in the tumor microenvironment. Cell Mol Immunol 2020;17:27-35. [Crossref] [PubMed]
41. Tang S, Qin C, Hu H, et al. Immune Checkpoint Inhibitors in Non-Small Cell Lung Cancer: Progress, Challenges, and Prospects. Cells 2022;11:320. [Crossref] [PubMed]
42. Genova C, Dellepiane C, Carrega P, et al. Therapeutic Implications of Tumor Microenvironment in Lung Cancer: Focus on Immune Checkpoint Blockade. Front Immunol 2021;12:799455. [Crossref] [PubMed]
43. Wu J, Wang Y. Role of TNFSF9 bidirectional signal transduction in antitumor immunotherapy. Eur J Pharmacol 2022;928:175097. [Crossref] [PubMed]
44. Geuijen C, Tacken P, Wang LC, et al. A human CD137×PD-L1 bispecific antibody promotes anti-tumor immunity via context-dependent T cell costimulation and checkpoint blockade. Nat Commun 2021;12:4445. [Crossref] [PubMed]
45. Rakaee M, Adib E, Ricciuti B, et al. Association of Machine Learning-Based Assessment of Tumor-Infiltrating Lymphocytes on Standard Histologic Images With Outcomes of Immunotherapy in Patients With NSCLC. JAMA Oncol 2023;9:51-60. [Crossref] [PubMed]
46. Bonaventura P, Shekarian T, Alcazer V, et al. Cold Tumors: A Therapeutic Challenge for Immunotherapy. Front Immunol 2019;10:168. [Crossref] [PubMed]
47. Whiteside TL, Demaria S, Rodriguez-Ruiz ME, et al. Emerging Opportunities and Challenges in Cancer Immunotherapy. Clin Cancer Res 2016;22:1845-55. [Crossref] [PubMed]
48. Sha D, Jin Z, Budczies J, et al. Tumor Mutational Burden as a Predictive Biomarker in Solid Tumors. Cancer Discov 2020;10:1808-25. [Crossref] [PubMed]
49. Bravaccini S, Bronte G, Ulivi P. TMB in NSCLC: A Broken Dream? Int J Mol Sci 2021;22:6536. [Crossref] [PubMed]

(English Language Editor: J. Gray)