A systematic search of the GEO database for LUAD gene expression data was conducted. Our analysis included four LUAD expression profile cohorts: GSE10072, GSE32863, GSE40791, and GSE68465 (Landi et al., 2008; Shedden et al., 2008; Selamat et al., 2012; Zhang et al., 2012). We downloaded the matrix files for each GEO cohort for further analysis. All acquired RNA expression profiles were normalized using log2 (TPM +1) transformation. In total, we included 654 LUAD tumor tissue samples and 226 normal samples. To correct for batch effects caused by non-biological technical variation, we used the ComBat algorithm from the svaR package.

Differentially expressed genes (DEGs) between normal tissues and LUAD samples were identified using empirical Bayes methods in the limma R package (Ritchie et al., 2015). Significant DEGs were identified based on the cutoff of P < 0.05 and a dynamic threshold of |logFC| ≥ 0.828 (calculated as |logFC| ≥ [mean(|logFC|) + 2sd(|logFC|)]) (Rong et al., 2020). Additionally, we performed a protein-protein interaction (PPI) analysis using the STRING database (https://cn.string-db.org) to explore the relationships among the DEGs. The PPI network was constructed with Cytoscape software, and key hub genes were identified by calculating Degree scores using the CytoHubba plugin (Shannon et al., 2003). Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (GO/KEGG) pathway enrichment analyses were deployed to identify the potential functions of DEGs.

By utilizing the ConsensusClusterPlus (Wilkerson and Hayes, 2010), an unsupervised consensus clustering analysis was performed to investigate the expression profile data of prognostic genes. The optimal number of clusters was determined based on the cumulative distribution curve, and the process was repeated 1,000 times to ensure the stability of the results. We have uploaded the complete code and parameter settings of ComBat and ConsensusClusterPlus to GitHub. For more details, please visit: https://github.com/mikelu1997/consensus.

Multiple algorithms, including ESTIMATE (Liu et al., 2020) and TIMER (Zhou et al., 2020), were employed to analyze the tumor immune microenvironment. The ESTIMATE algorithm was used to evaluate ESTIMATE scores, immune scores, and stromal scores, while the TIMER algorithms was utilized to assess the abundance of immune cell infiltration across different categories.

The TIDE score is an online tool (http://tide.dfci.harvard.edu/) designed to evaluate the efficacy of immunotherapy in different risk groups and assess the likelihood of tumor immune evasion (Wang et al., 2022). A higher TIDE score indicates poorer response to ICIs. In this study, we identified LCN2, MUC4, CDH17, COMP, GREM1, COL11A1, MMP9, THBS2, COL1A2, COL3A1, PLAU, CNTF, CXCL13, CXCL9, and CCL19 genes were associated with response to ICIs (Chen D. et al., 2021). We further extracted the expression values of these 15 genes to analyze the expression patterns of immune checkpoint-related genes in different groups.

We developed a RiskScore signature to comprehensively evaluate the role of IRGs in patient prognosis and response to immunotherapy. The glmnet R package was used to perform LASSO-Cox regression analysis with 10-fold cross-validation. Ultimately, a linear equation for immunoregulation-related gene index (IRGI) was constructed to predict overall survival (OS) of early-stage LUAD patients: RiskScore = [coef (Siegel et al., 2023) × GeneExp (Siegel et al., 2023)] + [coef (Zeng H. et al., 2024) × GeneExp (Zeng H. et al., 2024)] + …+ [coef(i) × GeneExp(i)] (Zhou et al., 2019). Kaplan-Meier curves were generated using the survival and survminer R packages to conduct prognostic analysis and to evaluate 2-year, 3-year, and 4-year survival rates in the test cohort. To validate the effect of IRGI risk model in predicting patient prognosis and therapeutic response to immunotherapy, external datasets, including GSE72094, the Cho cohort (Cho et al., 2020), and VanAllen cohort (Van Allen et al., 2015), were utilized.

The GSE229353 (Hui et al., 2023) data set includes single-cell RNA-sequencing (scRNA-seq) data from six LUAD samples treated with immunotherapy. Quality control procedures were applied to filter single cells based on the following criteria: 1) Each gene was expressed in more than 200 genes and in more than 3 cells; 2) The number of genes expressed in a cell ranged from 500 to 50,000; 3) The percentage of mitochondrial genes in a single cell was less than 15% (Yates et al., 2025). After applying these filters, a total of 15,293 cells were retained. The SEURAT R package was used to process the scRNA-seq data. To eliminate batch effects across the six samples, the “FindIntegrationAnchors” and “IntegrateData” functions were applied. Next, the “ScaleData” function was used to scale scRNA-seq data, and principal component analysis (PCA) was performed to reduce dimensionality. Moreover, “FindClusters” and “FindNeighbors” functions were executed with parameters set to dim = 9 and resolution = 0.1 to cluster the single cells into different subgroups. We obtained cell markers for different cell types from CellMarker 2.0 (Hu et al., 2023) and used these markers to annotate the single cells. Finally, single cells were visualized using t-SNE plot generated by the “RunTSNE” function.

The XGBoost model-related packages were loaded in R, and data was divided into training (70%) and validation (30%) sets using a stratified random sampling method. The XGBoost machine learning model was trained on the training set, and the predictive performance was evaluated using area under the curve (AUC) values. To achieve a deeper understanding of the contribution of each feature in the model prediction, we applied SHAP to interpret the XGBoost model. The shap package was used to calculate the contribution values of each feature for individual samples. SHAP value analysis was conducted to rank the importance of each feature and identify its directional impact within the model. Furthermore, we generated feature importance plots, dependence plots, and force plots to visually represent the SHAP analysis results.

Total RNA was extracted from normal human lung epithelial cells (BEAS-2B) and A549 lung adenocarcinoma cell lines using Trizol (Thermo Fisher Scientific, Sweden) (all cell lines were maintained according to the supplier’s recommendations). Reverse transcription was performed according to the instructions of the Takara Kit (Takara, Maebashi, Japan). SYBR Green Premix Ex Taq kit (Takara) was used to quantitatively detect PLAU and GRME1 in normal cells and lung adenocarcinoma cells. QRT-PCR was performed using the Roche Applied Science Light Cycler 480 instrument. The cycle threshold (CT) (2–△△CT) method was employed to calculate the data. The normalized expression levels were compared with those of β-actin using the comparative CT method. The primers used in this study are presented in Supplementary Table S3.

All data calculations and statistical analyses were performed using R programming (version 4.2.1). The Kaplan-Meier method was applied for survival analysis, and the predictive performance of the risk model was evaluated using time-dependent receiver operating characteristic (ROC) curves via the timeROC package. For comparisons between two groups of continuous variables, independent Student’s t-tests were used to analyze normally distributed variables, while the Mann-Whitney U test was employed for non-normally distributed variables. All two-sided statistical p-values less than 0.05 were considered statistically significant.

The flow chart of this research is presented in Figure 1. First, to explore DEGs and their functional roles in regulating TME, we integrated four publicly available datasets in LUAD: GSE10072, GSE32863, GSE40791, and GSE68465. A comprehensive comparison of gene expression between normal and tumor samples determined 484 DEGs, consisting of 276 downregulated and 208 upregulated genes (Figure 2A). To further identify key genes that regulate TME, we constructed a PPI network using the STRING database (network type: physical subnetwork, minimum required interaction score: 0.4) and Cytoscape software (Cytohubba plugin, Degree >15). This analysis identified a core PPI network comprising 178 genes (Figure 2B). These genes were hypothesized to be involved in the regulation of TME and LUAD progression. To assess whether these crucial genes could distinguish between LUAD and normal tissue, we performed PCA to reduce the dimensionality of the gene expression data. The PCA results revealed two completely separate clusters, with normal samples and LUAD tumor samples forming distinct groups, suggesting significant differences in the expression patterns of key genes between the two groups (Figure 2C). A heatmap of these genes further emphasized their differential expression patterns across the samples (Figure 2D). Moreover, we conducted GO and KEGG enrichment analyses on the identified genes to uncover their biological functions. The enriched results revealed that the DEGs were significantly associated with immune response, metabolic processes, cell signaling, and cell proliferation and regulation (Figures 2E,F). These findings suggest that the identified DEGs may play pivotal roles in immune regulation and tumor progression, thereby highlighting their potential as therapeutic targets or biomarkers in LUAD.

We performed WGCNA to further identify co-expression modules and determine immune-regulatory key genes related to LUAD progression. A co-expression network was constructed based on mRNA expression data, which allowed to uncover tumor-associated modules and genes (Figure 3A). To construct the scale-free network, the optimal soft-thresholding power was determined to be β = 3, as evidenced by the analysis of network topology (Figures 3B,C). The co-expression network revealed a total of 10 distinct modules. Notably, the black module exhibited the strongest correlation with immune clusters (r = 0.38, P = 7e-31, Figure 3D). This module contained 1,181 genes that significantly correlated with tumor biology (Supplementary Table S1). To further visualize the gene relationships within the black module, we generated a Topological Overlap Matrix (TOM) heatmap, which illustrated the degree of association and co-expression between genes (Figure 3E). Next, we intersected the DEGs, WGCNA tumor-related modules, and immune-related genes from the ImmPort database (https://www.immport.org/shared/) (Xin et al., 2021), leading to the identification of 15 candidate IRGs (Figure 3F). These genes were significant upregulation in LUAD tumor samples compared to normal tissues (Figure 3G). Finally, we investigated the mutation landscape of these IRGs. Among the 517 LUAD samples, 35.0% of the patients exhibited mutations in at least one of the key genes. COL11A1 was the most frequently occurred mutation, followed by COL3A1, THBS2,COL1A2, CDH17, MUC4, MMP9, and COMP (Figure 3H). These findings highlight the pivotal role of these IRGs in LUAD, including their prognostic value as immune regulatory factors in cancer biology, and their predictive effect for immunotherapy efficacy.

In order to investigate the relationship between immune-related prognostic genes and LUAD subtypes, we performed consensus clustering analysis. Based on the cumulative distribution function (CDF) values, we classified LUAD patients into two distinct clusters (k = 2, Figures 4A–D). Cluster 1 (n = 438) exhibited high expression levels of these immune-related genes and was defined as the IRG-high pattern, while Cluster 2 (n = 215) displayed low expression levels, categorized as the IRG-low pattern (Figure 4E). Notably, survival analysis revealed a significant survival difference between the two patterns. The IRG-high pattern was associated with a more favorable prognosis, while the IRG-low pattern was linked to poorer survival outcomes (Figure 4F). These results suggest that the expression levels of immune-related genes play a crucial role in influencing the prognosis of LUAD patients and could serve as a potential biomarker for patient stratification and therapeutictargeting.

We utilized the ESTIMATE algorithm to compare the immune microenvironment between the two molecular patterns. The results demonstrated that the IRG-high model exhibited higher ESTIMATE, immune, and stromal scores compared to the IRG-low model (Figure 4G). Additionally, according to TIMER scores, the IRG-high pattern was associated with increased abundance of immune cells, including B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, and dendritic cells (DC), as shown in Figure 4H. TIDE scores were also calculated to assess the potential clinical efficacy of immune therapy in different risk groups (Fu et al., 2020). TIDE is a computational tool that reflects the potential for tumor immune evasion. Higher TIDE scores are generally associated with poorer efficacy of ICIs. In this study, we revealed that the IRG-low pattern had significantly higher TIDE scores compared to the IRG-high pattern, suggesting that patients in the IRG-high group might benefit more from ICI therapy (Figures 4I,J). These findings underscore the differences in the tumor immune microenvironment between the IRG-low and high groups, as well as their impact in determining the prognosis of LUAD patients and their potential response to immunotherapy. The differential TIDE scores further highlight the potential for personalized treatment strategies based on immune-related molecular patterns.

Using Spearman correlation analysis, we observed that most of the IRGs exhibited positive correlations with each other (Figure 5A). Furthermore, expression of some molecules, such as GREM1, MMP9, PLAU, CXCL13, CXCL9, and CCL19, was positively correlated with immune checkpoint molecules, especially programmed cell death-ligand 1 (PD-L1) and cytotoxic T-lymphocyte associated protein-4 (CTLA-4), with PLAU and GREM1 showing particularly significant correlations (Figures 5B–D). Based on the results from LASSO-Cox regression analysis (Figures 5E,F), we selected PLAU and GREM1 as the optimal indicators for constructing the prediction model. GO/KEGG analysis indicated that both PLAU and GREM1 were enriched in the pathways of immune response regulation, immune cell infiltration, and the remodeling of the extracellular matrix. PLAU was involved in angiogenesis, cell migration and invasion, and GREM1 was associated with the transforming growth factor (TGF)-beta signaling pathway, which was important for tumor progression and immune evasion (Supplementary Figure S1). We also explored single nucleotide variants and copy number variations correlation with the targeted PLAU and GREM1 gene expression, and both high PLAU and GREM1 expression was significantly associated with more TP53, CSMD3, and LRP1B mutations (Supplementary Figure S2). Upon reviewing the Human Protein Atlas (HPA) database, we analyzed immunohistochemistry (IHC) data for these key molecules and qualitatively observed obvious expression differences between normal and LUAD samples (Figure 5G).

Furthermore, we integrated four databases, including hTFtarget, ENCODE, GTRD and ChIP-Atlas, and identified five upstream regulatory transcription factors, namely CTCF, EP300, MAX, RAD21, and KDM4A, which are common to both PLAU and GREM1(Supplementary Figure S3). The risk score was calculated as follows: RiskScore = (−0.095) × (GREM1 expression) + (−0.002) × (PLAU expression). These results were then incorporated into the IRGI risk model. We calculated the IRGI risk score for each patient in the LUAD dataset based on the expression levels and risk coefficients of two key IRGs (PLAU and GREM1; Figure 5H). Low-risk patients had significantly longer OS than high-risk patients (P < 0.001, hazard ratio [HR] = 0.49 [0.36–0.66]; Figure 5I). The ROC curve indicated the accuracy of this model, yielding an AUC value of 0.794 (95% confidence interval [CI]: 0.743–0.844) (Figure 5J). Furthermore, in the TCGA validation cohort, we confirmed OS of low-risk patients was significantly higher than that of high-risk patients (p = 0.032, HR = 1.52 [1.04–2.24]) (Figure 5K). To further evaluate the efficacy of this model, we utilized time-dependent ROC curves, which demonstrated good predictive ability of this IRGI risk model over a 4-year period (Figure 5L). We also performed chemotherapeutic drug predictions for patients with LUAD in the high-risk and low-risk groups to provide optimal treatment regimens. Cisplatin, Cyclophosphamide, Gemcitabine, Paclitaxel, Vinorelbine and Doceaxel were more effective for the high-risk group (Supplementary Figure S4).

To validate the predictive value of the IRGI risk model in ICI therapy, we applied independent GSE and external datasets. In the GSE72094 validation cohort, patients with high expression of PLAU (P = 0.006; Figure 6A) and GREM1 (P < 0.001; Figure 6B) showed significantly worse OS compared to patients with low expression of these genes. However, when evaluating ICI therapy response in the Cho and VanAllen cohorts, high expression of PLAU and GREM1 was significantly associated with better survival outcomes, suggesting that these genes might be predictive of improved efficacy of anti-programmed death 1 (PD-1)/PD-L1 and anti-CTLA-4 therapy (Figures 6C–F). Further validation was carried out using the scRNA-seq data from GSE229353, which included 6 NSCLC patients who underwent neoadjuvant chemotherapy or combination immunotherapy with chemotherapy. A total of 26,930 single cells were initially obtained, with 15,293 cells remaining after stringent quality control procedures for subsequent analysis. The gene expression levels were normalized, and the single cells were clustered into 20 distinct clusters (Figure 6G). Moreover, based on markers from CellMarker2.0, the clusters were categorized into 7 cell types (Figures 6H–J; Supplementary Table S2). The majority of the cells were identified as monocytes/macrophages and T-cells (both CD4⁺ and CD8⁺ T-cells), which are key participants in the immune response. These findings support that ICI therapy leads to increased T-cell infiltration and a better immune response, further corroborating the potential of the IRGI risk model in predicting the efficacy of immunotherapy.

Sankey diagram showing the distribution of patients across different characteristics, and most of patients with the IRG-high pattern had IRGI low-risk score (Figure 7A). The combined model integrating IRG signature patterns with IRGI risk score also confirmed that patients with IRG-high pattern and IRGI low-risk score exhibited the significantly best survival, while patients with IRG-low pattern and IRGI high-risk score had the worst survival outcomes (P < 0.001; Figure 7B). To further improve the predictive effect, we developed XGBoost machine learning predictive model using pre-identified IRGs, and AUC value of the XGBoost model was strikingly improved to 0.975 (95% CI, 0.962–0.989; Figure 7C). SHAP analysis indicated that GREM1, CCL19, COMP genes ranked in the top three important contributions to the XGBoost model predictions (Figure 7D). The impact of each IRG feature on the prediction is illustrated in SHAP waterfall and force plots (Figures 7E,F). With the basal prediction of −0.172 in the XGBoost model, COL11A1, MMP9, CCL19, and PLAU increased the overall prediction by +1.36, +1.35, +0.802, and +0.605, respectively (Figures 7E,F). SHAP importance scatter plot demonstrated the SHAP values and feature values for a specific IRG (Figure 7G). Each dot represents a patient, and the color indicates the feature value. GREM1 had the highest median SHAP values, indicating the greatest impact on the overall prediction.

The RT-qPCR results showed the mRNA expression levels of the two genes (GREM1 and PLAU) involved in our IRGI risk model in LUAD cells and normal lung epithelial cells. Specifically, compared with BEAS-2B cells, the mRNA expressions of GREM1 and PLAU in A549 cells were significantly upregulated, consistent with our database analysis findings (Figure 8).