Twelve single-cell RNA sequencing samples from five lung adenocarcinoma patients in the GSE127465 database were included in this study. Bulk sequencing data, mutation data, and clinical information for lung adenocarcinoma patients were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. Microarray sequencing data and clinical information of the GSE31210 database were downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) as an external independent validation set for the prognostic signature. Furthermore, the GSE126044 database was used as the immunotherapy response validation cohort (anti-PD-1 treatment). Detailed clinical information for TCGA and GSE31210 database is listed in Supplementary Table S1.

The available cell clustering and cell type annotation information of GSE127465 was used in the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/) (Sun et al., 2021), and single-cell analysis was performed using the R package ‘Seurat’ (Hao et al., 2021). Based on the annotation results, M2 macrophages were extracted for further analysis. In this study, the number of hypervariable genes was set to 2000, and the resolution for cell clustering to 0.6. Principal component analysis (PCA) was conducted based on 2000 hypervariable genes. In addition, dimensionality reduction of single-cell data was used by the t-distributed stochastic neighbor embedding (tSNE) method (Kobak and Berens, 2019), and the ‘FindAllMarkers’ algorithm was performed to search for characteristic differentially expressed genes among different cell clusters. R package ‘monocle2’ was used for differentiation trajectory and pseudotime analysis of M2 macrophages (Zhou et al., 2022). Subsequently, the ‘BEAM’ (branched expression analysis modeling) statistical algorithm was used to identify heterogeneous genes that play a key role in the differentiation of M2 macrophages (Wang et al., 2022). R package ‘GSVA’ and ‘scMetabolism’ determined the enrichment of signaling pathways at the single-cell level (Hanzelmann et al., 2013; Wu et al., 2022). Finally, cell-to-cell communication analysis was performed using the R package ‘iTALK’ (Wang Y. et al., 2019).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment assays were performed using the R package ‘clusterprofiler’ (Wu T. et al., 2021; Kanehisa et al., 2021). GO analyses include three parts: biological process (BP), cell composition (CC), and molecular function (MF). In addition, the R package ‘limma’ was used to identify differentially expressed genes in the prognostic signature between high- and low-risk groups (Ritchie et al., 2015). Gene Set Variation Analysis (GSVA), an unsupervised algorithm, was performed to calculate enrichment scores of hallmark gene sets (Molecular Signatures Database (MSigDB), http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp).

Univariate Cox analysis of heterogeneous genes associated with M2 macrophage differentiation was performed using the R package ‘survival’ to screen for prognostic genes associated with lung adenocarcinoma among them. In addition, the TCGA datasets were randomly grouped on a 3: 2 scale by the “sample_frac” function in the R package “dplyr” to obtain the training and testing datasets. Based on these prognostic genes, a Lasso-Cox regression analysis was performed using the R package ‘glmnet’ (Friedman et al., 2010; Wang H. et al., 2019). Next, a multivariate Cox stepwise regression approach was performed to construct a prognostic model related to the differentiation trajectory of M2 macrophages. The formula for the signature was: risk score = [Coef (gene 1) x gene Exp (1)] + [Coef (gene 2) x gene Exp (2)] +. . . . . . + [Coef (gene 1) x gene Exp (i)]. R packages ‘survival’ and ‘survminer’ were used for Kaplan-Meier prognostic analysis. The R package ‘timeROC’ was used to assess AUC values for time-dependent ROC curves. To further improve the prediction efficiency of the risk score, the R package ‘rms’ was used to combine the pathological stage with the risk score to construct a more accurate nomogram (Balachandran et al., 2015).

Estimate, EPIC, MCPcounter, TIMER, and ssGSEA algorithms were used to calculate immune infiltration abundance in lung adenocarcinoma patients with different risk scores. Among these five algorithms, the Estimate algorithm calculated the Estimate score, tumor purity, immune score, and stromal score (Yoshihara et al., 2013). The EPIC algorithm calculated the abundance of seven immune cell types (Racle et al., 2017). The MCPcounter algorithm calculated the abundance of 10 immune cell types (Becht et al., 2016). The TIMER algorithm calculated the abundance of six immune cell types (Li et al., 2020). Subsequently, the ssGSEA algorithm was utilized to calculate the enrichment score of 24 immune cell gene sets (Bindea et al., 2013). These algorithms revealed differences in immune cell infiltration abundance between high- and low-risk groups.

Based on the TCGA mutation data (maf format), mutations in the high- and low-risk groups were analyzed using the R package ‘Maftools’ and mapped waterfall plots (Mayakonda et al., 2018). In addition, the TIDE algorithm was utilized to analyze the sensitivity of high- and low-risk groups to immune checkpoint inhibitors (ICIs) (Jiang et al., 2018). Based on the signature formula constructed above, the risk score of lung adenocarcinoma patients was calculated in GSE126044 to assess differences in immunotherapy efficacy. Subsequently, the chemotherapeutic drug sensitivity (IC50) of patients in the high- and low-risk groups was analyzed using the R package ‘pRRophetic’ (Geeleher et al., 2014). These studies helped provide personalized treatment strategies.

Statistical analyses were performed using R software (v 4.1.3), and the results were visualized using the R packages. For non-normally distributed data, Wilcoxon rank-sum test, as a non-parametric test method, was used to examine the differences between the two groups of continuous variables, while for three and more groups we used the Kruskal–Wallis test for statistical testing. Using the Cox regression method, Kaplan-Meier prognostic analysis calculated the hazard ratio (HR). A two-sided p < 0.05 was considered statistically significant. Spearman method was applied for correlation analysis (*p < 0.05, **p < 0.01).

To better understand the heterogeneity of M2 macrophages in the TME of lung adenocarcinoma and its potential value for prognosis and drug treatment screening, we extracted and analyzed lung adenocarcinoma samples at the single-cell level. Based on the meta-information and cell type annotation information from the GSE127465 database on the TISCH website, we extracted 12 lung adenocarcinoma samples that had been quality-controlled and standardized. To overcome technical noise in scRNA-seq data, we performed Principal Component Analysis (PCA), and each principal component (PC) was considered a “meta-feature” (Supplementary Figure S1A). We identified the most suitable number of PCs (24 PCs) for downstream analysis by calculating the standard deviation of each principal component (Supplementary Figure S1B). In addition, we used tSNE, a nonlinear dimensionality reduction algorithm, to demonstrate the distribution of single-cell data from 12 lung adenocarcinoma samples (Figure 1A). We also examined the cell distribution of lung adenocarcinoma patients of different ages and clinical stages (Supplementary Figure S1C, D). Subsequently, we used the R package ‘Seurat’ to classify the cells in 12 samples into 24 clusters (Figure 1B). These 24 clusters can also be categorized into 13 cell types: B cells, CD4⁺ Tn cells, CD8⁺ Tex cells, endothelial cells, fibroblasts, M2 macrophages, malignant cells, mast cells, monocyte cells, neutrophils cells, NK cells, and plasma cells (Figure 1C). We counted the frequency of these 13 types of immune cells and found a higher proportion of M2 macrophages in each of the 12 samples (Figure 1D). We then extracted the M2 macrophages and re-clustered them using the ‘Seurat’ package (Figure 1E). The result suggested that M2 macrophages can be classified into 9 clusters (0–8) based on different molecular markers (Figure 1F). M2 macrophages of lung adenocarcinoma patients showed considerable heterogeneity in the different clusters. Patient seven had a more significant proportion of cluster0 and cluster4 in the M2 macrophages (Figure 1G). Patient six and patient four had a high percentage of cluster1. Patient five had a larger ratio of cluster5 and cluster6. However, patient three had a greater portion of cluster3 and cluster6. Based on the above clustering results, we analyzed the signature genes of the 9 clusters of M2 macrophages using the ‘FindAllMarkers’ algorithm and visualized the results with scatter plots and heat maps (Figures 1H,I). We observed that the genes had distinct expression differences between different clusters.

The heterogeneity among M2 macrophages was intriguing, and to further investigate the biological functions of essential genes in the differentiation of M2 macrophages, we performed differentiation trajectory analyses. We found that M2 macrophages can be divided into five differentiation states (Figure 2A). Meanwhile, we found that subpopulations of M2-type macrophages were differentially distributed on differentiation trajectory (Supplementary Figure S1E). Subsequently, we performed pseudotime analysis on M2 macrophages (Figures 2B,C). Purple indicated the initial state of cell differentiation, and yellow indicated the terminal state. Since state2 had a smaller number of cells and a high overlap with state4 in differentiation trajectory and pseudotime, we merged state2 with state4 as a whole. State1 accounted for 32%, state2&4 accounted for 11%, state3 accounted for 31%, and state5 accounted for 26% of all M2 macrophages (Figure 2D). In addition, we found that heterogeneous genes associated with the differentiation trajectory of M2 macrophages illustrated four expression patterns (Figure 2E).

We performed KEGG and GO analysis based on these statistically significant heterogeneous genes (p < 0.0001). KEGG enrichment analysis indicated that these heterogeneous genes were involved in activating numerous signaling pathways (Figure 2F). Examples include the Chemokine, IL-17, HIF-1, B cell receptor, and PI3K-Akt signaling pathway. The diverse activation levels of these pathways suggested that the different states of M2 macrophages might play distinct roles in the progression of lung adenocarcinoma. GO analysis phenotyped the heterogeneous genes in biological processes, cellular components, and molecular functions (Figures 2G–I, Supplementary Table S2). The results revealed that these genes activate multiple immune cells in the TME, suggesting crosstalk between M2 macrophages and immune cells.

To further examine the functional differences between the different classes of M2 macrophages we distinguished, we performed a GSVA enrichment analysis. First, we verified the potential differences in molecular mechanisms among the 9 clusters of M2 macrophages. Although M2 macrophages exhibited pro-oncogenic activity, the different clusters of M2 macrophages showed significant differences in the activation levels of the 50 gene sets contained in Hallmark (Figure 3A). For instance, cluster7 and cluster8 have lower activation levels in numerous pathways than the other seven types of clusters. This suggested that cluster7 and cluster8 might be relatively weak in oncogenic activities. To better understand this heterogeneity within M2 macrophages, we also performed an enrichment analysis of the different states in the differentiation trajectory (Figure 3B). The results confirmed substantial heterogeneity among the different states. It was found that state1 had significantly higher enrichment levels in the epithelial-mesenchymal transition (EMT). In comparison, states2&4 had significantly higher activation levels in the reactive oxygen species (ROS) pathway, mitotic spindle, and interferon-gamma response. State3 was significantly activated on the apical junction and notch signaling and had the lowest activation on E2F targets, G2M checkpoint, and the genes upregulated by ultraviolet (UV) radiation. State5 was notably enriched in angiogenesis, hypoxia, and genes upregulated by KRAS signaling. Subsequently, we analyzed the differences in the metabolic activities of M2 macrophages in these states using the R package ‘scMetabolism’ (Figure 3C). The red color in the heat map corresponds to the higher activation level, and we can identify that the different states of M2 macrophages have distinct metabolic levels. This heterogeneity in metabolic levels may reveal differences in the functional levels of the M2 macrophages in different states. Finally, we analyzed the cellular communication between M2 macrophages and the other cells in the TME (immune checkpoints, cytokines, and growth factors). Regarding immune checkpoints (Figures 3D,E), we found that the M2 macrophages of state1 mainly communicated with NK cells, neutrophils, monocyte cells, and cDC2 cells. M2 macrophages in state2&4 communicated with malignant cells, CD4⁺ Tn cells, and CD8⁺ Tex cells in addition to the above 4 cells. M2 macrophages in state3 had extensive communication with plasma cells and malignant cells. In contrast, M2 macrophages in state5 communicated predominantly with monocyte cells. In terms of cytokines, the different states of M2 macrophages also have diverse levels of cellular communication (Figures 3F,G). However, in terms of cell growth factors, there was no significant difference in the level of communication between these states of M2 macrophages and immune cells (Figures 3H,I). Together, the above results demonstrated heterogeneity in the level of cellular communication among different states of M2 macrophages.

After extracting heterogeneous genes associated with the differentiation trajectory of M2 macrophages, we performed a univariate Cox analysis in the TCGA database and obtained 289 prognosis-related genes of lung adenocarcinoma. Utilizing Lasso-Cox with multivariate stepwise regression, we constructed a prognostic model for lung adenocarcinoma in the training set (Figures 4A,B). In addition, Figure 4C shows the coefficients of the seven genes incorporated into the formula. These seven genes were: CCL20, BIRC3, CRYL1, SLC46A3, MAP3K8, TMED10, and CCR2. The formula of the model was: risk score = [0.144635490549132 * CCL20 Exp] + [0.343940732365844 * BIRC3 Exp] + [-0.291003194691483 * CRYL1 Exp] + [-0.241116620394011 * SLC46A3 Exp] + [-0.317179567218481 * MAP3K8 Exp] + [0.813174718146032 * TMED10 Exp] + [-0.353159490585916 * CCR2 Exp]. Using this formula, we calculated the risk score values in the testing set, the entire TCGA database, and the external validation dataset GSE31210. Lung adenocarcinoma patients with high-risk scores had higher deaths in these four datasets (Figures 4D–G). The heat map results also indicated that the expression of the above seven genes had significant differences in the two risk groups. We then performed Kaplan-Meier prognostic analysis to explore the potential value of our constructed model for patients with lung adenocarcinoma. Patients with high-risk scores had a worse prognosis in the training set (HR = 3.38, p < 0.001), the testing set (HR = 2.48, p = 0.001), the entire TCGA set (HR = 2.89, p < 0.001), and the GSE31210 database (HR = 5.34, p < 0.001, Figures 4H–K). We also performed a time-dependent ROC curve analysis on these four databases to judge our model’s accuracy in predicting prognosis. The AUC values of our model in the training set for 1-, 3-, and 5-years overall survivals were 0.742, 0.762, and 0.758, respectively (Figure 4L–O). In the testing set, the AUC values of our model for 1-, 3-, and 5-years survival were 0.741, 0.739, and 0.741, respectively. In the entire TCGA dataset, the AUC values of our model for 1-, 3-, and 5-years survival were 0.741, 0.753, and 0.741, respectively. While, in the GSE31210 dataset, the AUC values of our model for 1-, 3-, and 5-years survival were 0.784, 0.658, and 0.708, respectively.

Furthermore, we analyzed the relationship between numerous clinical factors and risk scores (Figures 5A–D). It was found that there was no statistical difference in the risk score between the two groups of patients aged≥65 years and those aged<65 years. The difference between the two risk groups of patients with different gender and smoking history was also not statistically different. However, the differences among patients with different pathologic stages were statistically significant. Patients with high pathological stages tended to have higher risk scores. In addition, we performed Kaplan-Meier prognostic analysis of lung adenocarcinoma patients with different clinical characteristics separately (Figure 5E–L). The results showed that patients with high-risk scores had a poor prognosis in all age groups (≥65, <65), all gender groups (Male, Female), and all pathological stage groups (Pathological stage IandII, Pathological stage IIIandIV). There was a significant prognostic difference between the high- and low-risk groups among the smoking group, with patients in the high-risk group having a poor prognosis (HR = 2.83, p < 0.001), however, there was no significant prognostic difference between the high- and low-risk groups in the non-smoking group. We also performed the univariate and multivariate Cox regression analyses regarding the risk scores. Our results showed that the risk score is an independent prognostic factor for lung adenocarcinoma and can be used as a clinical parameter to determine the prognosis of patients (Figure 5M, N).

As demonstrated in the above study, lung adenocarcinoma patients had significantly different prognoses between the high- and low-risk groups. Differential expression analysis was performed for the high-risk versus low-risk groups to investigate the mechanisms involved. Firstly, we used the R package ‘limma’ to identify differentially expressed genes (|FC|>1.5, FDR<0.05) and mapped the volcano (Figure 6A), Then, we performed enrichment analysis using GO and KEGG (Figure 6B, Supplementary Table S3). KEGG results showed significant enrichment in the cell cycle and IL-17 signaling pathway (Figure 6C). In addition, we performed GSVA enrichment analysis and plotted heat map and histogram for lung adenocarcinoma patients in high- and low-risk groups. Fifty gene sets from Hallmark were selected for GSVA analysis. By comparing the enrichment scores of the two groups in these 50 gene sets, 25 gene sets showed a statistically significant difference (Figures 6D,E). These results might explain the underlying mechanism for the difference in the prognosis of lung adenocarcinoma patients with different risk scores.

To understand the differences in TME among patients with different risk scores, we also analyzed the abundance of immune infiltration in high- and low-risk groups. Analysis using the Estimate algorithm showed that the high-risk group had higher tumor purity but lower immune and stromal scores (Figures 7A–D). This suggested that patients with high-risk scores exhibited a state that promoted tumor escape due to the lack of anti-tumor immune cells in the TME. Then, to investigate the immune status of the two risk groups in more detail, we used the ssGSEA algorithm to calculate the infiltration abundance of 24 immune cell types (Figure 7E). The results revealed that the high-risk group generally had a lower infiltration abundance of immune cells, including B cells, T cells, CD8⁺ T cells, NK cells, DC cells, and mast cells, compared to the low-risk group. Spearman correlation analysis also showed that the abundance of immune cells was negatively correlated with the risk scores for almost all immune cell types except Th2 cells (Figure 7F). We also used the MCPcounter, TIMER, and EPIC algorithms to confirm these results. By plotting the heat map, we visualized that lung adenocarcinoma patients in the high-risk group had lower levels of immune infiltration (Figure 7G).

Since gene mutation status significantly affects tumor formation and progression, we performed a mutation analysis of lung adenocarcinoma patients in high- and low-risk groups. We further analyzed the differences between the top 20 genes with the highest mutation frequencies in the high and low-risk groups respectively (Figures 8A,B). The frequency of TP53 mutations in the high-risk group was 55.1%, while it was 48.2% in the low-risk group. The frequency of KRAS mutations was 31.8% in the high-risk group versus 27.3% in the low-risk group. The difference in gene mutation frequency might be the reason for the poorer prognosis in the high-risk group. Based on the above analyses, we explored the treatment strategies for different risk scores in depth. Since immunotherapy is commonly used in lung adenocarcinoma, we first calculated the response to immunotherapy in high- and low-risk groups using the TIDE algorithm (Figure 8C). The results suggested that the low-risk group had a better treatment response upon immunotherapy (OR = 1.731, p = 0.004). We further validated this result using the GSE126044 database. Results showed that patients with CR/PR after immunotherapy had a significantly lower risk score than SD/PD (p = 0.038), suggesting that patients with a low-risk score are susceptible to benefit from immunotherapy (anti-PD-1 treatment) (Figure 8D). However, the high-risk group had a worse prognosis, so we performed a prediction of response to chemotherapeutic agents for patients in the high- and low-risk groups. We utilized the R package ‘pRRophetic’ to calculate the IC₅₀ of chemotherapeutic drugs (Figure 8E). We found that the high-risk group had better sensitivity to A-443654, BIBW-2992, Docetaxel, Paclitaxel, Embelin, and RO-3306. Taken together, we provided a personalized treatment option for clinical reference by predicting the treatment effect of patients in different subgroups and compensated for the poorer effect of immunotherapy in the high-risk group.

The above study indicated that the risk score could act as an independent prognostic factor that can be used to determine the prognosis of patients. Therefore, to further improve our signature’s predictive efficiency, we constructed a nomogram based on the TCGA database, incorporating factors such as pathological stage and risk scores (Figure 9A). We could visualize the risk assessment by calculating the score and assessing the outcome probability for each patient. In addition, we tested the predictive efficacy of the nomogram using a calibration plot (Figure 9B). The nomogram’s 1-, 3-, and 5-years survival predictions were more accurate than the theoretical values. To further compare the predictive efficacy of the nomogram score with other clinical factors, we plotted the ROC curve. The nomogram score had the highest AUC value, and the predictive efficiency was further improved based on the risk score (Figure 9C). By performing the time-dependent ROC curve based on the nomogram score, we found that the AUC values of 1-, 3-, and 5-years overall survival for lung adenocarcinoma patients were 0.786, 0.789, and 0.776, respectively. To verify the accuracy of the nomogram, we constructed the nomogram again based on the external dataset GSE31210. Time-dependent ROC analysis showed that the AUC values of 1-, 3-, and 5-years overall survival for lung adenocarcinoma patients were 0.918, 0.777, and 0.744, respectively (Figures 9D,E). These results suggested that the nomogram we constructed had good accuracy. We also performed Kaplan-Meier prognostic analyses using the TCGA and external database GSE31210. The results indicated that the high nomogram score group had a significantly worse prognosis than the low score group (Figures 9F,G). In summary, the nomogram significantly improved the accuracy of determining the survival status of lung adenocarcinoma patients.