Integrative Analysis of Single-Cell and Bulk RNA Sequencing Reveals Prognostic Characteristics of Macrophage Polarization-Related Genes in Lung Adenocarcinoma

Introduction

Lung cancer is a global health concern, being the most diagnosed cancer and the leading cause of cancer-related deaths worldwide. Among the various histological subtypes of lung cancer, lung adenocarcinoma (LUAD) is the most frequent.^1,2 According to the revised classification of lung cancer by the World Health Organization in 2021, LUAD is categorized into adenomatous atypical hyperplasia (AAH), adenocarcinoma in situ (AIS), minimally invasive adenocarcinoma (MIA), and invasive adenocarcinoma (AI).³ Research has identified smoking as the major risk factor for LUAD, followed by air pollution, occupational exposure, and genetic predisposition.^4–7 While existing therapies such as surgical resection, chemotherapy, radiotherapy, and immunotherapy have been used to treat advanced LUAD, their efficacy in improving the prognosis of patients is limited.⁸ Therefore, identifying novel characteristic genes and developing prognostic models may contribute to understanding the pathogenesis and molecular mechanisms of LUAD, as well as provide potential targets for clinical diagnosis and treatment.

Infiltrating immune and stromal cells are essential components for the function of tumor microenvironment (TME). Increasing data suggest that they have significant roles in tumor growth, recurrence, metastasis, and chemotherapy resistance. The significant role of infiltrating immune cells provides a novelty insight for revealing the tumor biology and will be an essential prognostic factor. Tumor-associated macrophages (TAMs) are a type of macrophage that infiltrates tumor tissues and constitute the major component of infiltrating leukocytes in the TME.^9,10 LUAD exhibits inter-patient and intra-tumoral heterogeneity in both tumor cells and TME.¹¹ The TME components vary markedly among different tumors and play crucial roles in tumor initiation, progression, and metastasis.^12,13 TAMs can be polarized into three main phenotypes: resting-state macrophages M0 (precursors of polarized macrophages), classically activated macrophages M1 (induced by Toll-like receptor and IFN-γ) and alternatively activated macrophages M2 (induced by IL-4 and IL-13).^14–16 In tumor progression, M1 macrophages typically play anti-tumor roles, while M2 macrophages exert pro-tumor functions.¹⁷ M1 macrophages kill tumor cells through direct-mediated cytotoxicity and antibody-dependent cell-mediated cytotoxicity (ADCC).¹⁸ In contrast, M2 macrophages can secrete tumor-stimulating factors, inhibit T cell-mediated immune response, and promote proliferation, invasion and angiogenesis.^19–22 Previous studies have shown that macrophages are involved in the development of lung adenocarcinoma (LUAD). For example, M2 macrophages enhance the expression of vascular endothelial growth factor (VEGF) family members A and C, which promote tumor angiogenesis and lymphangiogenesis. Additionally, M2 macrophages can promote LUAD cell epithelial-mesenchymal transition (EMT) by activating the ERK1/2/Fra-1/slug signaling pathway.^23,24 Zheng et al found that the lower density of M1 at the tumor center and higher proximity of tumor cells to M2 at the invasive margin predicted poor prognosis in lung cancer. It is worth mentioning that macrophages are highly plastic and heterogeneous.¹⁶ Affected by the microenvironment changes, TAMs can be converted from one polarized phenotype to another. Inducing the phenotypic switch from M2 to M1 macrophages has been regarded as a novel therapeutic idea for malignant tumors. The study of macrophage polarization-related genes (MPRGs) helps to understand the underlying mechanism of macrophage polarization and shed new light on polarization control.

Integrating single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing has significant potential in identifying biomarkers. scRNA-seq provides a more detailed and accurate understanding of the heterogeneity of cell populations by analyzing gene expression at the single-cell level.^25,26 In contrast, bulk RNA sequencing is more cost-effective, has higher throughput, and provides greater coverage.^27,28 By integrating scRNA-seq data to identify specific cell types or subpopulations with bulk RNA sequencing data to identify differentially expressed genes or pathways, more accurate and reliable biomarkers can be identified.^29–32 These biomarkers can be used for diagnosis, prognosis, and personalized treatment of diseases. Therefore, integrating scRNA-seq and bulk RNA sequencing can provide a more comprehensive understanding of the molecular mechanisms underlying disease and improve patient outcomes.

In this study, the researchers utilized a combination of bulk RNA sequencing and single-cell RNA sequencing to investigate the role of macrophage polarization-related genes in LUAD. By integrating these techniques with bioinformatics methods, the team was able to identify novel biomarkers for malignant tumors and construct a new prognostic model for LUAD. This approach enabled a more comprehensive understanding of the heterogeneity of a particular cell subset in the tumor microenvironment and provided valuable insights into the molecular mechanisms underlying LUAD progression.

Materials and Methods

Data Sources

Bulk RNA-seq data and clinical information of 59 normal tissue samples and 513 LUAD tissue samples with survival information in the TCGA-LUAD dataset were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The scRNA-seq data of 11 LUAD tissue samples in the GSE131907 dataset, the bulk RNA-seq data, and survival information of 398 LUAD tissue samples in the GSE72094 dataset were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/). Referring to the reported studies,³³ 35 macrophage polarization genes (MPGs) were obtained (Supplementary Table 1). The flow chart of this research was shown in Figure 1.

Figure 1 The flow chart of the research. Abbreviations are defined as follows. Abbreviations: DEGs, Differential expressed genes; K-M, Kaplan-Meier; DE-MPRGs, Differentially expressed macrophage polarization-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, Protein-protein interaction (PPI) network; LASSO, Least absolute shrinkage and the selection operator.

Analyzing of scRNA-Seq Data of LUAD

The “Seurat” package³⁴ was utilized to obtain the core cells, and several quality criteria were applied, including the exclusion of genes detected in only 3 or fewer cells, low-quality cells with less than 200 detected genes, and cells with mitochondrial expression genes ≥5%. The top 2000 highly variable genes were identified using the FindVariableFeatures and NormalizeData functions, and Principal Component Analysis (PCA) was performed on all samples, with the top 50 principal components selected for further analysis. UMAP algorithm³⁵ was used to perform overall dimensionality reduction on the samples based on the top 50 principal components. Marker genes were identified using the FindAllMarkers function in the Seurat package, and the ‘SingleR’ package³⁶ was used to annotate and visualize different cell types and extract high-margin genes in macrophage clusters. Furthermore, cellular communication between different cell types was analyzed using the ‘cellchat’ package.³⁷ These rigorous computational analyses helped to identify and characterize different cell types, extract key marker genes, and provide insights into cellular communication and molecular mechanisms underlying LUAD progression.³⁷

Acquisition of Macrophage Polarization-Related Subtypes in the TCGA-LUAD Dataset

The R package “ConsensusClusterPlus”³⁸ was employed to identify macrophage polarization-related subtypes in the TCGA-LUAD dataset based on the expression of 35 MPGs. The overall survival (OS) among different subtypes was explored using the “Survival” package. Additionally, differential expression analysis was performed using “DESeq2”³⁹ to obtain the differentially expressed genes (DEGs) between the two subtypes, which were defined as MPRGs. The screening criteria for MPRGs were adj.P-value <0.05 and |log2FoldChange|>0.5. Similarly, the DEGs between normal and LUAD tissue samples in the TCGA-LUAD dataset were identified using the DESeq2 package. Volcano plots and heatmaps were generated using the R packages “ggplot2” and “pheatmap”, respectively, to visualize DEGs and MPRGs.

Identification and Biological Functional Analysis of Differentially Expressed and Highly Mutated MPRGs in TCGA-LUAD

The DEGs, MPRGs, and hypervariable genes of macrophages were intersected to identify the differentially expressed and highly mutated MPRGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of these intersecting genes were performed using the “clusterProfiler” package,³⁸ and the results were visualized through bubble plots created using the “ggplot2” package in R. Additionally, the protein-protein interaction (PPI) network of the intersecting genes was constructed using the STRING database (https://string-db.org).

Construction and Validation of the Prognostic Risk Model

The intersection genes were subjected to univariate Cox analysis in the TCGA-LUAD dataset to identify prognosis-related genes with a significance level of P < 0.05. Then, using the least absolute shrinkage and selection operator (LASSO)⁴⁰ analysis, the characteristic genes were identified and sequentially eliminated. The resulting genes were used to construct a prognostic risk model, and the LUAD patients were divided into two groups based on the median risk score: high-risk and low-risk groups. The Kaplan-Meier (KM) curves were used to display the difference in overall survival (OS) between the two groups. To assess the prognostic capability of the model, the ROC curves were displayed using the “survROC” package.⁴¹ Finally, the external GSE72094 dataset was used to validate the stability of the model. The “rms” package was used to construct a nomogram that predicts survival probability based on the characteristic genes. The calibration curve was used to validate whether the nomogram can be used as an optimal model for clinical decision-making.

Assessment of the Prognostic Risk Model

The distribution and differences in clinical characteristics between the high and low-risk groups were analyzed using the Wilcoxon test after combining clinical characteristics such as age, sex, T-stage, N-stage, M-stage, and survival status with risk scores. To determine whether clinical characteristics and risk scores were independent predictive factors for LUAD patients, univariate and multifactorial Cox analyses were performed. The “rms” package was used to construct a nomogram that predicts survival probability based on independent prognostic criteria. The calibration curve and decision curve analysis (DCA) were both used to validate whether the nomogram can be used as an optimal model for clinical decision-making. The “maftools” package was utilized to analyze the mutation frequencies of genes between the two groups, while the Spearman correlation coefficient was used to analyze the correlation between characteristic genes and significantly mutated genes. Additionally, all genes in the high/low-risk groups were analyzed for correlation with characteristic genes, and the functional pathways of characteristic genes were analyzed using gene set enrichment analysis (GSEA).

Immune Infiltration Analysis in LUAD

The CIBERSORT algorithm⁴² was used to assess the infiltration of immune cells in all samples between the high and low-risk groups. The differences in immune cell abundance between the two groups were examined using the Wilcoxon test. Spearman correlation analysis was performed to investigate the correlations between the 22 immune cell types.

Analysis of Characteristic Genes in Key Cell Clusters

To identify macrophage subtypes, we extracted the macrophage clusters in the GSE131907 dataset annotated in Analyzing of scRNA-Seq Data of LUAD and re-annotated them by dimensionality reduction clustering and subgroup annotation to subdivide the macrophage clusters. The resulting clusters were further analyzed using the Monocle package to perform trajectory analysis.⁴³ From these clusters, key cell clusters of macrophage types were selected for further analysis. The expression levels of characteristic genes were then analyzed among subtypes using consensus clustering and the Wilcoxon test. Additionally, GSEA was performed to understand the potential functions of the characteristic genes among subtypes.

Identification of Transcription Factors (TFs) Affecting the Expressions of Characteristic Genes

To explore the potential regulatory mechanisms affecting the characteristic genes in different subpopulations, we adopted the single-cell regulatory network inference and clustering (SCENIC) approach⁴⁴ to predict specific transcription factors (TFs). The Wilcoxon test was then utilized to analyze differences in the expression of TFs among subtypes.

Collection of Tissues and Analysis of Characteristic Genes Expression

In this study, we investigated the expression of 7-MPRGs in LUAD tissues. A total of 5 paired LUAD and adjacent normal lung tissue samples were collected from the Second Affiliated Hospital of Xi’an Jiaotong University. And this study was approved by the Ethics Committee of the Second Affiliated Hospital of Xi’an Jiaotong University. Written informed consent was obtained from all patients. Total RNA was extracted from tissue samples using the FAST1000 kit (Pioneer, China), and RNA reverse transcription was performed using the PrimeScript™ RT-PCR Kit (Takara, Japan). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed using the TB Green^® Premix Ex Taq™II (Takara, Japan) kit. All primers used in this study were synthesized by Sangon Biotech Co. Ltd. (Shanghai, China) (Supplementary Table 2).

Results

Integration and Clustering of scRNA-Seq Data

To screen for highly variable genes in LUAD scRNA-seq, we first performed quality control and screened a total of 29,747 core cells (Supplementary Figure 1). Following data normalization, 2000 highly variable genes were identified for subsequent analysis (Figure 2A). PCA analysis was performed on the single-cell samples, and the sample dispersion distribution was logical (Figure 2B). The top 50 principal components were selected for subsequent analysis (Figure 2C and D). The core cells were then divided into 23 independent cell clusters using the UMAP algorithm (Figure 2E), and the marker genes in each cluster are shown in Figure 2F. The “singleR” software package was used to annotate the different cell clusters, resulting in nine cell clusters (Figure 1G). Among them, the macrophage cluster had 424 highly variable genes (Supplementary Table 3). The number of ligand-receptor interactions is shown in Figure 1H, and macrophages could interact with other cells.

Figure 2 Highly variable gene screening and cell clustering of LUAD scRNA-seq data from GSE131907. (A) 2000 highly variable genes. (B) PCA results. (C-D) P-value of first 50 principal components (PCs). (E) Cluster diagram of the core cells. (F) Heatmap showing the top 5 marker genes in each cell cluster. (G) 23 cell clusters were annotated into 9 cell types. (H) CellChat diagram of each cell type.

Identification of Macrophage Polarization-Related Subtypes

Based on the expression of 35 MPGs, we classified 513 LUAD patients into two subtypes (Figure 3A-C). It was observed that cluster2 had a worse prognosis than cluster1 (Figure 3D). This indicates that the degree of expression of MPGs had an impact on the prognosis of LUAD patients, adding rationality to the subsequent analysis. We found 3920 DEGs in the two subtypes, which were named macrophage polarization-related genes (MPRGs) (Figure 3E and F, Supplementary Table 4). In the TCGA-LUAD dataset, we observed a total of 9561 DEGs between normal/LUAD tissue samples (Figure 3G and H, Supplementary Table 5).

Figure 3 Identification of macrophage polarization-related subtypes of TCGA-LUAD samples. (A) Cumulative distribution function (CDF) curve of K=2-9. (B) The relative change in area under the CDF curve of K=2-9. (C) The consensus matrix heatmap of K=2. (D) Kaplan-Meier survival curves of the two subtypes. (E-F) Volcano plot and heat map representing DEGs (|log2FC|>0.5, adj.P-value<0.05) in the two subtypes. (G-H) Volcano plot and heat map of DEGs (|log2FC|>0.5, adj.P-value<0.05) between TCGA normal and LUAD samples.

The Function Analysis of Differentially Expressed and Highly Mutated MPRGs

Firstly, considering that the degree of expression of MPGs had an impact on the prognosis of LUAD patients, we took the intersection of the two datasets of MPGs and DEGs. However, the intersection of MPGs and DEGs may not be expressed in macrophages. To ensure that the selected differential genes were expressed in macrophages, it was changed to take intersecting genes in DEGs, MPRGs, and highly variable genes, defined these 28 genes as differentially expressed MPRGs (DE-MPRG) (Figure 4A, Supplementary Table 6). These 28 genes are the genes whose expression in macrophages of LUAD tumor tissue is different from that in normal tissue. Biological process (BP) analysis revealed that these genes are involved in the activation of cysteine-type endopeptidase activity, apoptosis process, and leukocyte aggregation. In terms of cellular components (CC), these genes were associated with phagocytosis-related CC, including secretory granule lumen, cytoplasmic vesicle lumen, and vesicle lumen. Regarding molecular function (MF), these genes were enriched for immunoglobulin-binding functions (Figure 4B). The intersecting genes were also enriched in pathways such as IL-17 and Apoptosis (Figure 4C). The final PPI network of the intersecting genes containing 21 points with 26 interacting edges was obtained (Figure 4D). ALDH2 and PLB01 were found to only interact with each other in this network.

Figure 4 KEGG and GO analysis of differentially expressed and highly mutated MPRGs. (A) Venn diagram showing the 28 genes obtained by intersecting. (B) GO terms on BP, CC and MF levels enriched for the crossover genes. (C) KEGG pathways enriched for the 28 genes obtained. (D) PPI network of 26 intersecting genes.

Construction and Validation of a Seven-Characteristic Gene-Based Prognostic Risk Model

In this study, univariate Cox regression analysis was performed using the TCGA-LUAD dataset, and 13 genes were found to be significantly associated with OS (Figure 5A). Next, a LASSO algorithm was used to screen for genes to construct a model (Figure 5B). The results showed that seven characteristic genes were screened at the lowest cross-validation error. The risk score was calculated using the following formula: Risk score = (−0.15120133 × RGS13) + (−0.10521700 × ADRB2) + (0.10214897 × DDIT4) + (−0.01788635 × MS4A2) + (−0.07556347 × ALDH2) + (−0.09135983 × CTSH) + (0.19599794 × PKM). Based on the median risk score (0.9332402), patients were divided into high-risk and low-risk groups. The K-M analysis revealed that the prognosis was better for low-risk patients, and the ROC curves indicated that the model was feasible (Figure 5C). Furthermore, the prognostic risk model demonstrated strong predictive power in the GSE72094 dataset (Figure 5D). To construct a reliable prognostic model, a nomogram was developed based on these 7 characteristic genes (Figure 5E). The nomogram indicated that the overall survival rate of patients decreased as the overall score increased. The calibration curve of the nomogram was close to the diagonal, indicating that the prediction of the model was accurate and reliable (Figure 4F).

Figure 5 Construction and validation of differentially expressed and highly mutated MPRGs associated with LUAD. (A) Univariate Cox regression analysis of the TCGA-LUAD screened the differentially expressed and highly mutated MPRGs associated with LUAD. (B) 7 genes were further screened by LASSO regression analysis. The trajectory of each independent variable and the confidence interval under each lambda were showed. (C) Predictive value of 7 characteristic genes for LUAD in TCGA training set. (D) Predictive value of 7 characteristic genes for LUAD in GSE72094 validation set.

Assessment of Prognostic Risk Models

Next, we will construct a nomogram model based on the seven characterized genes of LUAD. Risk scores differed between age, gender, T.stage, and N.stage (Supplementary Figure 2). Further analysis using univariate and multifactorial Cox regression revealed that the risk score, T.stage, N.stage, and M.stage were all independent prognostic factors for LUAD (Figure 6A and B). To construct a reliable prognostic model, a nomogram was developed based on these independent prognostic factors (Figure 6C). The nomogram indicated that the overall survival rate of patients decreased as the overall score increased. The calibration curve of the nomogram was close to the diagonal, indicating that the prediction of the model was accurate and reliable (Figure 6D). Therefore, the nomogram was considered the optimal model for predicting the prognosis of LUAD patients (Figure 6E).

Figure 6 Construction of a nomogram model based on the 7 characteristic genes for LUAD. (A-B) Univariate and multifactorial Cox analysis of the TCGA-LUAD cohort. (C) Probabilistic nomogram for predicting 1-, 2-, and 3-year survival probability of patients with LUAD. (D) Calibration plot for conformance testing between 1–3 years survival predictions and actual outcomes. (E) Time-dependent ROC curve showing the risk score had significant prognostic value.

GSEA Analysis of the High/Low-Risk Groups

The top 5 mutated genes in the high/low-risk groups were TP53, TTN, CSMD3, MUC16, and RYR2 (Supplemental Figure 3). In addition, ZNF536 and DNAH9 were mutated in the low-risk group, while COL11A1 was mutated in the high-risk group, and these genes were highly correlated with the characteristic genes (Figure 7A). These findings suggested that these mutated genes may be involved in the regulation of characteristic genes, which in turn affects the prognosis of LUAD patients. Further gene set enrichment analysis (GSEA) showed that genes associated with characteristic genes were mainly enriched in the cell cycle pathway, which is closely related to tumor development and is an important cause of tumorigenesis and malignant growth due to impaired cycle regulation (Figure 7B).

Figure 7 Mutated genes analysis and immune infiltration analysis in the high/low-risk groups. (A) Matrix heatmaps of the correlation analysis results. (B) GSEA enrichment analysis of high/low-risk groups. (C) Violin plots depicting differential analysis of immune cells in different risk groups. ****p<0.0001, ***p<0.001, *p<0.05.

Correlation Analysis of Risk Scores and Immune Microenvironment

Supplemental Figure 4A shows the results of analyzing the percentage of 22 immune cell components in the TCGA-LUAD dataset after extracting features using data from scRNA-seq. Differential immune cells were identified and 11 of them were significantly different between the high- and low-risk groups (Figure 7C). Furthermore, the Spearman analysis revealed the interactions between immune cells (Supplemental Figure 4B), and Monocyte, T cell CD8+, B cell plasma, and Mast cell activated were found to be highly correlated with other immune cells.

Role of Characteristic Genes in Macrophage Types Cell Clusters

To further subdivide macrophages, the macrophage clusters were further annotated and downscaled, resulting in the identification of five distinct classes of cell clusters (Figure 8A). Analysis of gene expression revealed that ALDH2, PKM, and CTSH were present in macrophages, while the remaining four characteristic genes were present in NK cells (Figure 8B). These findings suggested a potential association between macrophages and NK cells. The differentiation of macrophage clusters followed four directions (Figure 8C). Notably, NK cells: CD56hiCD62L+ differentiated earlier and were mainly present on wo branches, while macrophages were present on two branches and differentiated later (Figure 8D).

Figure 8 Reannotation and analysis of the macrophage clusters. (A) Five clusters were identified after further annotation. (B) Expression of characteristic genes in different cell clusters. (C) Pseudotime analysis of the macrophage clusters. (D) Cell annotation for the pseudotime analysis.

In examining the expression of genes in the various cell populations, three of the five cell populations annotated as subpopulations directly belonging to the macrophage type were selected as key cell populations for further follow-up analysis, given that our analysis was based on macrophages as a whole. The three cell clusters were Macrophage:Alveolar: B._anthacis_spore, Macrophage: Alveolar, Macrophage: monocyte-derived: M-CSF. Based on these three macrophage clusters, genes were extracted from the cells. The key cells were divided into two subpopulations based on the expression of characteristic genes (Figure 9A). Five characteristic genes, namely ADRB2, ALDH2, CTSH, DDIT4, and PKM, were found to be differentially expressed between the two subtypes. All of them were highly expressed in cluster2, except for PKM in cluster1 (Figure 9B). GSEA results between the two subpopulations suggested that these cells may be involved in the apoptotic signaling pathway and the ERBB signaling pathway, playing a role in LUAD (Figure 9C and D).

Figure 9 Expression and function of characteristic genes in macrophage types cell clusters. (A) Key cells from three macrophage-type-cell clusters were divided into two subpopulations. (B) Expression of characteristic genes in different cell subpopulations. (C-D) GSEA enrichment analysis of two subpopulations. ***p<0.001, **p<0.01, NS≥0.05.

BCLAF1 and MAX May Be Transcription Factors (TFs) Affecting the Expression of Characteristic Genes

To identify transcription factors for characterized genes, the 10 transcription factors with the highest regulon specificity score (R-DSS) were shown as the subgroup-specific TFs in Figure 10A. To visualize the expression of the top 10 R-DSS TFs in different cell clusters for the two subgroups, UMAP plots were performed (Figure 10B). Among the top 10 R-DSS TFs, only BCLAF1 and MAX were found to be differentially expressed in the two subpopulations (Figure 10C). These results suggest that these two TFs may be associated with the expression of the characteristic genes.

Figure 10 Identification of transcription factors of characteristic genes. (A) The top 10 TFs with the highest regulon specificity score in each subgroup. (B) Expression of the top 10 TFs in each cell cluster. (C) Expression of BCLAF1 and MAX in two subgroups. ****p<0.0001, **p<0.01.

Expression of Seven Characteristic Genes in LUAD Tumor Tissues and TCGA-Database

We conducted a differential expression analysis on the data obtained from TCGA, which indicated that RGS13, ADRB2, MS4A2, ALDH2, and CTSH were significantly downregulated in LUAD tumor tissues, while DDIT4 and PKM were upregulated (Figure 11A). To validate our findings, we further analyzed the expression of the 7-MPRGs in 5 pairs of LUAD and adjacent non-tumor lung tissues using qRT-PCR. Our results showed that the expression levels of MS4A2 and CTSH were significantly decreased in the LUAD tumor tissues when compared to the adjacent non-tumor lung tissues (Figure 11B).

Figure 11 Relative mRNA expression levels of 7-MPRGs in TCGA-database and clinical samples. (A) Differential expression in TCGA-LUAD-database. (B) Differential expression in five paired LUAD tissues and adjacent non-tumor lung tissues. ***p<0.001, *p<0.05.

Discussion

In recent years, there have been significant improvements in the diagnosis and treatment of LUAD. However, the prognosis for many patients is still poor, especially those diagnosed at an advanced stage of the disease.⁴⁵ It is now widely recognized that the tumor microenvironment plays a critical role in tumor progression, and the interactions between immune cells and cancer cells are key factors in this process.^46,47

In this study, we conducted a comprehensive analysis of gene expression profiles associated with macrophage polarization in LUAD tissues. By identifying and analyzing a panel of macrophage polarization-related genes, we were able to develop a risk score model that can be used to predict individual LUAD patient prognosis with high accuracy. This study provides valuable insights into the underlying mechanisms of LUAD progression and highlights the importance of macrophage polarization in tumor microenvironment. It also underscores the potential of using prognostic gene expression profiles to develop personalized treatment strategies for LUAD patients.

In recent years, several studies have highlighted the importance of post-transcriptional modifications such as polyadenylation and methylation alterations in the development of lung cancer, including LUAD. Differentially expressed genes, miRNAs, methylation sites, and APA-related genes associated with NSCLC (LUAD and LUSC) were identified in large datasets for the prediction of metastasis in cancer tissues.^48,49 Based on cell cluster analysis and microarray data, seven key genes associated with LUAD prognosis were identified in a July paper.⁵⁰ Considering the effect of macrophages, our study identified seven MPRGs using scRNA-seq and bulk RNA-seq data (RGS13, ADRB2, DDIT4, MS4A2, ALDH2, CTSH and PKM) as characterized genes and modeled the prognostic risk of LUAD. This comprehensive model, which includes both clinical and gene factors, can provide individualized assessment of LUAD prognosis and has better prediction ability for overall survival (OS) than single factors. While RGS13 has been associated with immune-related diseases, such as asthma and B lymphoma, its role in solid tumors and LUAD remains unknown and needs further validation.^51–54 ADRB2, on the other hand, has been shown to promote tumorigenesis and development in several malignancies, but its expression was significantly decreased in LUAD, and low expression was associated with shorter OS and disease-specific survival.^55–59 DDIT4 was confirmed to be an independent predictor of OS for LUAD patients and was shown to enhance the migration and invasion abilities of LUAD cells by activating the mTORC2/AKT pathway.^60,61 MS4A2, a favorable prognostic indicator, was identified as an independent prognostic marker in early-stage LUAD patients.⁶² ALDH2, a major enzyme for detoxification of endogenous acetaldehyde, was shown to inhibit malignant features of LUAD cells when overexpressed.^63–65 CTSH, a lysosomal cysteine protease, had a controversial role in tumors, with high levels observed in gliomas and melanoma but lower levels associated with better prognosis in colorectal cancer and thyroid carcinoma.^66–69 Interestingly, its levels were decreased in tumor tissues but elevated in sera of lung cancer patients.⁷⁰ PKM, a well-known promoter of LUAD progression, was found to be consistent with our risk model.^71–74 In conclusion, the role played by the RGS13 and CTSH genes in the progression of the tumors is unclear. Elevated expression levels of the DDIT4 and PKM genes may indicate a worsening of the tumors, whereas elevated expression levels of the ADRB2, MS4A2, and ALDH2 genes may indicate an improvement in the tumor condition. Overall, our study provides new insights into the molecular mechanisms underlying LUAD and offers a potential tool for personalized clinical management.

GSEA revealed that several pathways were enriched in the high-risk group, including cell cycle, ribosome, and spliceosome pathways. The deregulation of cell cycle can lead to abnormal cell proliferation, which is a basic mechanism of cancer. NSCLCs, including LUAD, often have cell cycle abnormalities that have a negative effect on prognosis.⁷⁵ Ribosomes and spliceosomes are essential for protein synthesis, sustaining tumor cell growth and proliferation. Dysregulation of the splicing process and ribosomal modifications are both closely related to tumorigenesis.^76–80 These findings provide insight into the mechanisms underlying the high-risk group of LUAD patients and could potentially lead to the development of targeted therapies to improve prognosis for this group.

In our study, we found that TP53, TTN, CSMD3, MUC16, and RYR2 were among the top 5 genes that were mutated in both groups. These somatic mutations may impact LUAD by affecting the expression of characteristic genes. For instance, DDIT4 is a transcriptional target of p53 following DNA damage.⁸¹ Additionally, a previous analysis involving 12 types of cancer found a strong association between TP53 mutations and PKM expression.⁸² Moreover, TTN expression has been shown to be related to CTSH,⁸³ while Hu et al demonstrated that mice with ALDH2 mutations had higher levels of phosphorylated RYR2 protein than wild-type mice.⁸⁴

To better understand how characteristic genes affect immune function, a CIBERSORT analysis was performed to assess the differences in immune infiltration between high- and low-risk groups. The results showed significant differences in 11 types of immune cell infiltration between the two groups, with certain immune cells being more abundant in the high-risk group. Specifically, macrophage M0, CD4+ memory activated T cells, and resting NK cells were found to be enriched in the high-risk group. M0 macrophages are thought to play a role in tumor initiation and progression by promoting angiogenesis, remodeling the extracellular matrix, and suppressing antitumor immune responses.^85–87 M0 macrophages can promote tumor growth and invasion by secreting cytokines and growth factors such as IL-6, IL-10, VEGF, and MMPs.⁸⁸ In LUAD, M0 macrophages have been shown to play a role in tumor progression and metastasis. Studies have found that M0 macrophages are increased in both tumor tissues and peripheral blood of LUAD patients, especially those with lymph node metastasis.⁸⁹ CD4+ memory-activated T cells have a crucial role in the immune response against tumors, including lung adenocarcinoma. However, their role in cancer is complex and may depend on various factors. In some instances, CD4+ memory-activated T cells may promote cancer growth and progression by suppressing the immune system or promoting inflammation. Regulatory T cells, which are a type of CD4+ T cell that suppresses immune responses, have been found to be enriched in some cancers, including lung cancer, and can contribute to tumor immune evasion and progression.^90–92 Current evidence suggests that the role of resting NK cells in predicting the prognosis of lung cancer is not entirely clear. Some studies have found that the density of resting NK cells in lung adenocarcinoma tissues is associated with a worse prognosis, possibly due to their immunosuppressive effects. Resting NK cells have been shown to produce cytokines such as TGF-β and IL-10, which can inhibit the function of other immune cells such as T cells and dendritic cells. Additionally, resting NK cells can promote the accumulation of regulatory T cells (Tregs), which can further suppress the function of other immune cells.⁹³ These immunosuppressive effects can contribute to tumor immune evasion and progression, ultimately leading to a worse prognosis for the patient. When elevated numbers of macrophage M0, CD4+ memory activated T cells, and resting NK cells are found in vivo, it should be realized that this may indicate a worsening of the LUAD tumor situation. Adopting appropriate means to inhibit the function of these immune cells may be a new therapeutic idea to inhibit tumor progression. In addition to these types of cells, the level of monocyte, T-cell CD8+, and mast cell infiltration may also have a correlative link to the prognostic level of patients with LUAD, and appears useful in predicting survival in advanced LUAD.⁹⁴

From the onset to the end of developmental differentiation, the major components of the cell clusters transition from NK cells to three types of macrophages, all of the three cell clusters are associated with the progression of LUAD. After dividing the three key clusters into two subgroups based on the expression levels of characteristic genes, GSEA analysis was performed, which revealed an enrichment of the ERBB signaling pathway. The ERBB tyrosine kinase family includes important members such as EGFR, HER2, HER3, and HER4. Oncogenic alterations of these genes have been known to activate tyrosine kinase aberrantly, thereby driving the development and progression of tumors and immune escape in LUAD.^95,96 These findings highlight the potential significance of ERBB signaling in the development of LUAD and could potentially lead to the development of new therapeutic targets for this disease.

The study identified two potential transcription factors, BCLAF1 and MAX, that may be associated with the differential expression of characteristic genes, but conclusive evidence for their relationship with the 7-MPRGs is currently lacking.

Finally, the researchers analyzed the differential expression of the seven model genes in LUAD RNA sequencing data from the TCGA database. The results were consistent with the risk model developed in the study. Furthermore, the reliability of most of the conclusions was verified using qRT-PCR. However, we acknowledged limitations such as an insufficient sample size and the potential impact of individual differences and other confounding factors. Therefore, it is important to consider these limitations when interpreting the results obtained from existing databases. Overall, these findings provide important insights into the potential roles of the identified genes in the development and progression of LUAD.

In summary, this study identified seven MPRGs (RGS13, ADRB2, DDIT4, MS4A2, ALDH2, CTSH, and PKM) as prognostic genes for LUAD through integrative analysis of scRNA-seq and bulk RNA-seq data. The developed model showed improved prediction compared to single clinical factors. However, the molecular mechanisms of these genes in LUAD and macrophage polarization remain unclear and require further experimental validation. Moreover, the application of the predictive model requires expansion of clinical samples and data support. The researchers will continue to focus on follow-up studies of these genes to gain more insights into their potential roles in LUAD.