The transcriptome profiles and corresponding clinical data of LADC patients were downloaded from TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases. TCGA-LADC comprised a total of 594 (535 tumor sample and 59 normal samples) adenocarcinoma cases. The main characteristics of the analysis included the following: age, sex, and pathologic stage; details of patient clinical information are described in Table 1 . GSE31210 comprised a total of 226 primary LADC of pathological stage I-II. The median age was 67 years and the range was 30-76 years, and there were 105 male and 121 female patients. GSE72094 comprised a total of 442 LADC cases. The median age was 70 years and the range was 38-89 years, and there were 202 male and 240 female patients. The samples from TCGA database were defined as the training set, the samples from the GSE31210 database were defined as the validation set, whereas LADC cases from GSE72094 were set as testing set. LADC patients with missing survival values or follow-up time < 1 days were excluded. A total of 1,114 samples (490 from TCGA, 226 from GSE31210, and 398 from GSE72094) were used in our study.

The protein expression data of the MRGs of LADC patients were evaluated using the Human Protein Atlas (https://www.proteinatlas.org/), which is derived from antibody-based protein profiling using immunohistochemistry.

MRGs were obtained from two MRG datasets (188 from GSE5099_CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_UP and 194 from GSE5099_CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_DN) (15) from the gene set enrichment analysis (GSEA) website (http://www.gsea-msigdb.org/gsea/msigdb/). Finally, a total of 382 MRGs were utilized in this study ( Supplementary Table 1 ).

The univariate Cox regression analysis was used to screen out the genes significantly correlated with OS based on the 382 MRGs (for P-values < 0.05). Next, the overlapped prognosis-related MRGs from TCGA and GEO databases were selected for the least absolute shrinkage and selection operator (Lasso) with ten-fold cross-validation which was subsequently applied using “glmnet” and “survival” packages. Afterwards, a multivariate Cox regression was applied out to select candidate OS-related MRGs and determine a prognostic signature. The risk score was calculated as follows: Risk score = β1 × (expression of RNA1) + β2 × (expression of RNA2) + ··· + βn × (expression of RNAn). Median MRG risk scores were used to differentiate high-risk subgroups from LADC patients. The regression coefficient (β) was obtained from the multivariate Cox regression analysis. Additionally, Kaplan-Meier survival analysis was conducted to assess the predictive performance of the prognostic signature. The signature was also externally validated with the GEO dataset using the same formula. All analyses were carried out using R language, version 4.0.5 (www.r-project.org).

Gene set enrichment analysis (GSEA) was conducted to investigate various molecular pathways differentially activated between high- and low-risk subgroups. False discovery rate q-values < 0.05 and |NES| > 1 were defined as statistically significant difference.

We utilized CIBERSORT algorithm to estimate the proportion of 22 immune cells between low- and high-risk patients. The sum of ratio of 22 immune cell types in each sample is 1.

To provide a more individualized predictive model, a nomogram combining the MRG signature and other clinical variables was constructed using the training cohort. The discrimination ability of the nomogram was assessed using the calibration curves and receiver operating characteristic (ROC) curves in the training and validation subgroups. Next, decision curve analysis (DCA) was applied to evaluate the clinical usefulness of the nomogram in the training and testing sets.

Continuous variables were presented as means ± SD, whereas categorical variables were displayed as percentages. The statistical significance of the differences in survival rate was measured using the log-rank test with a threshold of P-value < 0.05. Kaplan-Meier plots were applied to display the differences in survival duration. All statistical analyses were conducted using the software R (version 3.5.2) with corresponding packages.

To limit the candidate prognosis-related MRGs, the OS-related MRGs that were overlapping in the data from TCGA and GEO databases (23 MRGs) were identified ( Figure 1A ). Next, the 23 MRGs were used to the Lasso-Cox proportional hazards regression and ten-fold cross-validation to construct the best gene signature, and 14 candidate MRGs were ultimately identified ( Figures 1B, C ). Further, a multivariate Cox regression was used, and results exhibited that CTSL, ECT2, HCFC2, HNRNPK, LRIG1, OSBPL5, P4HA1, and TUBA4A were the independent prognostic MRGs ( Figure 1D ). We created a risk score according to the expression of the eight MRGs as follows: Risk score = CTSL × 0.001326639 + ECT2 × 0.023009173 - HCFC2 × 0.257179317 + HNRNPK × 0.010298027 - LRIG1 × 0.024832171 + OSBPL5 × 0.071303241 + P4HA1 × 0.007389189 + TUBA4A × 0.008003706.

We investigated the protein levels of these genes, detected using immunohistochemistry and obtained from the HPA database. The immumohistochemical staining of MRGs were based on the normal alveolar and tumor tissues. We discovered that the protein levels of HNRNPK, P4HA1, and TUBA4A were significantly upregulated, while CTSL, HNRNPK, and OSBPL5 were significantly downregulated in the tumor tissues compared to those of normal tissues ( Figure 2A ). The quantitative analysis results for each immunohistochemistry were show in Supplementary Material . We also investigated the expression level of the identified MRGs for normal and tumor samples using RNA-Seq data from the training set. The results showed roughly the same trend as the one observed for the protein expression ( Figure 2B ).

To identify the MRG signature suitable for LADC survival prediction, the LADC patients were separated into low-risk (N = 245) and high-risk groups (N = 245) based on the median risk score. Kaplan-Meier curve analysis depicted that high-risk patients were associated with poorer OS compared to the low-risk patients (P < 0.001, Figure 3A ). Furthermore, the ROC curve analysis demonstrated that the area under the ROC curve (AUC) of the prognostic MRG model at 1, 3, and 5 years were 0.707, 0.707, and 0.65 in the training set, respectively ( Figure 3B ). The distribution survival status and time for each patient from the training set were plotted with a division line indicating the risk score cutoffs ( Figures 3C, D ). Next, we conducted the univariate and multivariate Cox regression analyses to analyze the signature and clinicopathological independent indices predicting survival. The results showed that the MRG-based signature was able to be an independent prognostic indicator ( Figures 3E, F ). Its prediction capacity was also evaluated through calculating C‐index in the training set. The results showed that the C-index for the prediction of OS of the identified MRG signature was 0.72 (95% CI = 0.65–0.76).

We verified the prediction performance of this signature using LADC cases from the GSE31210 dataset. The risk score of each patient was calculated based on the indicated formula and separated into low-risk (N = 113) and high-risk groups (N=113) according to the median risk score. Kaplan-Meier analysis demonstrated that high-risk patients were associated with poorer OS and relapse-free survival (RFS) compared to the low-risk patients ( Figures 4A, B ). ROC curve analysis revealed that the AUC of the prognostic MRG model for predicting OS at 1, 3, and 5 years were 0.704, 0.625, and 0.677, respectively ( Figure 4C ). Furthermore, the ROC curve analysis revealed that the AUC of the prognostic MRG model for predicting RFS at 1, 3, and 5 years was 0.661, 0.619, and 0.647, respectively ( Figure 4D ). The survival status and time distribution for each patient from the validation set were plotted with a division line representing risk score cutoffs ( Figures 4E–G ). The expression profiles of the eight prognostic MRGs are illustrated in Figure 4H .

In addition, we further tested the prediction performance of this signature using LADC cases from the GSE72094 dataset. the LADC patients were separated into low-risk (N = 199) and high-risk groups (N = 199) based on the median risk score. Kaplan-Meier curve analysis depicted that high-risk patients were associated with poorer OS compared to the low-risk patients (P < 0.001, Figure 5A ). Furthermore, the ROC curve analysis demonstrated that the area under the ROC curve (AUC) of the prognostic MRG model at 1, 3, and 5 years were 0.621, 0.659, and 0.707 in the test set, respectively ( Figure 5B ). The distribution survival status and time for each patient from the testing set were plotted with a division line indicating the risk score cutoffs ( Figures 5C, D ).

The prognostic significance of the signature was further assessed using subgroups with different demographics and clinical characteristics from the training set, including age, sex, TNM stage, and pathological stage ( Supplementary Figure 2 ). We discovered that the MRG signature was useful for most subgroups ( Table 1 ). For the validation set, the model can also accurately predict the OS and RFS of low- and high-risk groups in these subgroups ( Supplementary Table 2 ).

We further investigated the association between the MRGs and clinicopathological characteristics such as age, gender, pathological stage, and TNM stage for patients in the training cohort. We observed the differential expression of CTSL, ECT2, HCFC2, HNRNPK, LRIG1, and TUBA4A ( Figures 6A–F ).

Additionally, we explored the differentially signaling pathways between high- and low-risk LADC patient through GSEA. In the high-risk group, the top five enriched GO terms included cadherin binding, cellular response to heat, chromosomal region, chromosome segregation, and mitotic nuclear division ( Figure 7A ). The top five enriched KEGG pathways were basal transcription factors, cell cycle, oocyte meiosis, p53 signaling pathway, and ubiquitin mediated proteolysis ( Figure 7B ).

We further investigate the tumor-infiltrating immune cells from the high- and low-risk patients using CIBERSORT. The results displayed that the tumors of high-risk patients exhibited a higher proportion of plasma cells, resting CD4 T memory cells, monocytes, resting dendritic cells, resting mast cells, and eosinophils; while activated CD4 T memory cells, M0 macrophages, M1 macrophages, and activated mast cells were higher in low-risk group ( Figures 8A, B ).

To accurately predict a certain clinical outcome, a nomogram was established by integrating stage, age, gender, and the eight MRGs using a Cox model ( Figure 9A ). For the training set, the AUCs of the nomogram at 1-, 3-, and 5-year OS were 0.769, 0.765, and 0.75, respectively ( Figure 9B ). In the validation set, the AUCs of the nomogram at 1-, 3-, and 5- year OS were 0.896, 0.779, and 0.738, respectively ( Figure 9C ). For convenient clinical application and visualization of the prognostic model, we established an easy-to-use web-based calculator (https://emergency.shinyapps.io/LADC/) for predicting the overall survival of LADC ( Figures 9D, E ).

We evaluate the predictive ability and clinical usefulness of the nomogram using calibration curves and DCA. The calibration plots displayed that the nomogram could accurately predict OS ( Supplementary Figures 3A, B ). In addition, we used DCA to assess the clinical usefulness of the nomogram; the results showed good clinical usefulness of this model both in the training and validation sets ( Supplementary Figures 3C, D ).