To investigate the TME of LUAD, we analyzed scRNA-seq data (23, 24) obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (GSE164789). In addition to analyzing gene expression patterns, the study combined clinical information and mutation data, which together provided a more comprehensive foundation for the analysis (25). Since the data were publicly accessible, ethical approval was not required.

The Seurat package (v4.4.0), along with R software (v4.3.3), was employed to process the raw gene expression data, thereby facilitating robust analysis (26). The DoubletFinder package (v2.0.3) was first employed to effectively eliminate possible doublet cells while also filtering out cells of low quality. To eliminate low-quality cells, we applied stringent filtering criteria: cells with nFeature (number of detected genes) outside the range of 300-6,000 or nCount (total number of counts) beyond 500-75,000 were excluded. Eliminated also were cells with mitochondrial gene expression more than 25% or red blood cell gene expression more than 5%.

The “Normalize Data” function in Seurat helped us to normalize the data, subsequently, the “Find Variable Features” tool found 2,000 extremely highly variable genes (27). The “ScaleData” tool helped to further standardize gene expression counts so suited for principal component analysis. Harmony R package (v1.2.0) lower sample batch effects (28). Moreover, the “Cell Cycle Scoring” tool indicated cell cycle phases to ensure appropriate characterizing of cellular states (29). For dimensionality reduction, clustering analyses were performed and gene expression patterns were visualized using UMAP (30, 31).

Cell subtypes were determined by performing clustering analysis in Seurat, utilizing the “FindNeighbors” and subsequently the “FindClusters” functions (32, 33). The diverse cell subtypes within the TME were accurately classified because clusters were annotated according to the average expression levels of established marker genes.

Stemness of cells was evaluated by computing gene set activity scores from scRNA-seq data using the AUCell method. Gene sets might be ranked based on degree of enrichment inside individual cells. This approach defined the variations in tumor cell subgroups and their stemness characteristics.

We used CytoTRACE, a computational tool for predicting developmental potential, to esitmate the differentiation status of tumor cell subtypes (34). Pseudotime trajectories were created from Monocle (v2.24.1), while lineage architecture was deduced from Slingshot (v2.8.0), which used a cluster-based minimal spanning tree (MST) (35). Smooth trajectory curves produced by the “getLineages” and “getCurves” let differentiation paths and branching events be seen.

Using “FindAllMarkers,” differentially expressed genes (DEGs) were found in Seurat. We investigated their biological relevance using Genomic Variant Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA). Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), and ClusterProfiler (v4.8.2) derived functional annotations (36, 37). Using an adjusted P-value threshold of 0.05, we detected a range of GO terms that showed significant enrichment, thereby helping to elucidate the regulatory processes associated with various tumor cell subtypes (38, 39).

The CellChat R package (v1.6.1) was used to analyze cellular communications within the TME (40). The “IdentifyCommunicationPatterns” function measured discrete communication patterns, and the variations in the intensity of interactions between cells were analyzed using the “netVisual diffInteraction” (41). This analysis revealed key signaling pathways and ligand-receptor pairs mediating cellular crosstalk in LUAD. To visually examine the signals that enter and exit every cell, we utilized scatter plots, heatmaps, and a range of other visualization techniques. Moreover, we considered associations among different cell types to be meaningful when the P value was less than 0.05.

We built clusters of TCs and single-cell regulatory networks using Python (v3.9.19) and the pySCENIC library (v0.12.1) (42). Using data on rankings of human gene motif from (https://resources.aertslab.org/cistarget/), we determined the top five TFs displaying the most notable expression changes.

The central goal of this investigation was to determine how effectively specific genes associated with various LUAD subtypes could predict patient survival outcomes. After identifying the most significant prognostic genes, we demonstrated that they served as strong predictors for constructing reliable prognostic models (43). This identification process was carried out through both univariate and multivariate Cox proportional hazards analyses (44). Thereafter, we employed a risk evaluation technique, in which the final risk value was obtained by adding up the multiplied values of gene expression amounts and their matching coefficients: Risk score = ∑ i n Xi × Yi . Furthermore, we derived optimal cut-off values using the “surv_cutpoint” function, which enabled a comparative analysis of prognostic differences among patient subgroups. Survival analysis was subsequently conducted by means of the “Survival” package in R (v4.3.3), and survival curves were also constructed using the “ggsurvplot” function so as to examine the predictive accuracy of our risk model (45). In addition, we assessed the model’s reliability by constructing ROC curves (46) with the package of “timeROC” (v0.4.0), thus offering a thorough assessment of both the model’s accuracy and its calibration.

To comprehensively analyze the immune landscape, we first used the CIBERSORT R software package (v0.1.0) to estimate immune cell scores for each patient (47, 48). Subsequently, we investigated the infiltration of immune cells in detail and also evaluated the differential expression of immune checkpoint-related genes (49). We further examined associations among risk scores, immune cell populations, and genes included in the model. Additionally, we made use of the Tumor Immune Dysfunction and Rejection (TIDE) platform (http://tide.dfci.harvard.edu) to predict patient responses to tumor immunotherapy.

The inferCNV R software package (v1.16.0, https://github.com/broadinstitute/inferCNV) served as our primary tool for CNV inference, enabling us to characterize the CNV patterns within multiple cell subtypes using the inferCNV algorithm. We assessed both relative gene expression and chromosome location data to infer CNV status across individual cells (50).

To strengthen the clinical significance of our results regarding drug applications, we evaluated the sensitivity of multiple agents. Specifically, we employed the “pRRophetic” package (v0.5) to estimate the half-maximal inhibitory concentration (IC50) for every individual case and subsequently compared drug sensitivities across the high-risk and low-risk groups.

The A549 cell strain was cultured in F-12K medium under standardized conditions, which included a temperature of 37 °C, 5% CO₂, and 95% humidity. The medium was fortified by incorporating 10% fetal calf serum along with 1% antimicrobials. Likewise, the NCI-H1975 cell strains were routinely preserved in RPMI-1640 medium with similar environmental conditions, and this medium also contained 10% fetal calf serum as well as 1% antimicrobials to ensure the cells attained their best growth rates.

A decrease in FOS expression was observed, which was partly a result of RNA synthesized by GenePharma (Suzhou, China). After introducing the FOS-targeting siRNAs (siFOS-1 and siFOS-2) combined with a control siRNA (si-NC), cells were distributed into 6- well chambers at a coverage level of 50%. Transfection was carried out in accordance with the instructions from the manufacturer, utilizing Lipofectamine 3000RNAiMAX (Invitrogen, USA). Each siRNA (RIbbio, China) was individually introduced into the cells.

Following achievement of 70% confluence in transfected cells, lysis was conducted with RIPA buffer. The lysates underwent centrifugation at 12,000 rpm for 15 minutes to clarify the extracts before being processed by SDS-PAGE. The separated proteins were then transferred onto a PVDF membrane, which was blocked at room temperature for 1.5 hours with 5% calf serum albumin. The membrane underwent incubation throughout the night at 4°C with an anti-FOS antibody, followed by a one-hour incubation with a secondary antibody. Finally, the presence of FOS protein was revealed using a chemiluminescent Western blot substrate.

To extract RNA, cell lysis was achieved with Trizol reagent. All water sources, laboratory instruments, and workspaces were maintained RNase-free throughout the protocol to ensure RNA stability. Once RNA was isolated, the PrimeScript™ reagent kit was used for reverse transcription. Finally, quantitative real-time PCR was conducted, using SYBR Green master mix as the fluorescent indicator for the detection of amplification.

To evaluate the survival of A549 and NCI-H1975 cells following transfection, the Cell Counting Kit-8 (CCK-8) method was applied (51). The cells were plated into 96-well chambers at 5 × 10³ cells per well and, following this, were left to adhere for 24 hours. Afterwards, every well gained 10 µL of CCK-8 reagent (A311-01, Vazyme), followed by an additional 2-hour incubation at 37°C in the absence of light. Light absorption at 450 nm was recorded each day from the first to the fourth day by means of a microtiter reader (A33978, Thermo). The collected light density data were averaged and graphed to illustrate changes in cell viability over the time course.

After the transfection procedure, A549 and NCI-H1975 cells were plated into 6-well chambers at 5 × 10³ cells per well. Cells were allowed to incubate for 24 hours at ambient temperature before adding the EdU working solution for another 2-hour period. Being rinsed two times with PBS, the cells were fixed in 4% paraformaldehyde for a duration of 15 minutes. After fixation, permeabilization was performed, and the cells were quenched for 15 minutes in a solution containing 0.5% Triton X-100 and 2 mg/ml glycine. Following this, 1 ml of Apollo solution and 1 ml of Hoechst stain were applied, and the cells were incubated for 30 minutes. Finally, fluorescence microscopy was employed at the end to image the cells and evaluate their proliferation.

Initially, cells were deprived of serum for 24 hours by incubating them in a serum-free medium. Following this, the cellular suspension was blended with Matrigel (BD Biosciences, USA) and the upper well of the Costar plate was loaded with the cell suspension, while the lower well received serum-containing medium to form a chemotactic gradient. Then cells were allowed to migrate and invade in a constant-temperature apparatus for a duration of 48 hours. Afterwards, 4% paraformaldehyde solution served to fix the cells, and crystal violet staining was employed to evaluate their invasion capacity visually.

Cells with stable transfection were cultured in 6- well chambers and incubated at 37°C within a moist environment supplemented with 5% CO₂ until they became confluent. Sterile 200 µL pipette tips were then used to create straight scratches in the cell monolayer. The wells were carefully washed with PBS to remove floating cells and debris. Following this, the cells were incubated with serum-free medium to facilitate migration. Photographs were captured at both the beginning (0 hours) and after 48 hours, and scratch widths were measured using Image-J for subsequent analysis.

The datasets were subjected to analysis using R (v4.3.3) and Python (v3.9.19). For comparisons between groups, Wilcoxon’s test and Pearson’s correlation coefficient were applied (52). We interpreted the significance at the statistical level following the criteria listed below: *P< 0.05, **P< 0.01, ***P< 0.001, and ****P< 0.0001, “ns” was used to denote results without statistical significance. These analytical methods and cutoffs were adopted to confirm the robustness of the findings.

To explore the ECs, we analyzed their CNVs using inferCNV with ECs as a reference ( Supplementary Figure 1 ). We examined the obtained dataset to explore the single-cell profile present in the LUAD environment. Figure 1 illustrated our procedure. By examining five localized adenocarcinomas and twenty-six infiltrating adenocarcinomas from GSE164789, we employed dimensionality reduction clustering with UMAP plots. This analysis initially exhibited how thirty-one individual samples were distributed and subsequently identified fourteen cell types: T and NK cells, epithelial cells (EPCs), macrophages, plasma cells, mesenchymal cells (MCs), endothelial cells (ECs), monocytes, fibroblasts, B cells, conventional dendritic cells type 2 (cDC2), proliferating cells, conventional dendritic cells type 1 (cDC1), myofibroblasts, and plasmacytoid dendritic cells (pDCs). Furthermore, we presented the distribution of cells within three groups (AIS, IAC, MIA) and across three different phases of the cell cycle (G1, G2/M, S) ( Figure 2A ). We examined the expression profiles of nCount RNA, nFeature RNA, cell stemness AUC, and pMT among the fourteen cell types and three different groups ( Figures 2B, C ). Additionally, Figure 2D illustrated the cellular distribution of the dataset. Volcano plots revealed the presence of different genes among these cell types ( Figure 2E ). In EPCs, TCP11L2, POM121C, HYAL2, HBP1, and SYNE2 were found to be upregulated, while BCOR, LHFPL6, IKZF2, SRPK1, and SRSF8 exhibited downregulation. In macrophages, EIF5AL1, ZNF468, AGO1, USP25, and TMEM165 showed upregulation, whereas UPRT, PDHX, MTFR1L, TRIOBP, and TMEM141 showed downregulation. We also provided the top five marker genes corresponding to the fourteen cell types, as depicted in Figure 2F . Through investigating the ratios of individual cell types in the three groups alongside cell cycle data, we found that each cell type distributed across separate phases within the cell cycle, with EPCs and macrophages constituting the main components of the IAC group, predominantly located in the differentiated G1 phase ( Figure 2G ).

We then conducted functional enrichment analysis. The analysis suggested that IAC group displayed increased activity in pathways connected to the regulation of apoptotic signaling, negative regulation of apoptotic signaling, entry into the host, viral life cycle, and biological processes involved in symbiotic interaction. EPCs indicated enrichment within pathways related to positive regulation of cell-cell adhesion and antigen processing and presentation, while macrophages were enriched in pathways associated with positive regulation of leukocyte cell-cell adhesion, antigen processing and presentation of peptide antigens via MHC class II, peptide antigen assembly with MHC class II protein complex, assembly of MHC class II protein complexes, and antigen processing and presentation of exogenous antigens ( Figure 2H ). Furthermore, EPCs demonstrated upregulation in pathways related to chemical homeostasis within tissues and surfactant homeostasis, whereas macrophages displayed significant upregulation in antigen processing and presentation of peptide antigens and antigen processing and presentation of peptide antigens via MHC class II ( Figure 2I ).

The involvement of the TME in tumorigenesis became widely recognized, as it harbored TCs that, through their interactions with other cells via the circulatory and lymphatic networks, played an essential part in both the emergence and advancement of cancer (53). We characterized the significance of TCs in the TME. Four subtypes of LUAD TCs were classified according to the expression of marker genes: C0 MAFF+ TCs, C1 GSTA1+ TCs, C2 GOLGA8B+ TC, and C3 CAV1+ TCs, and demonstrated the nCount RNA, nFeature RNA, cell stemness AUC, and pMT for each subtype ( Figure 3A ). Subsequently, we examined the primary distribution of TCs across different groups using UMAP plots ( Figure 3B ). Bubble plots illustrated the marker genes associated with various groups and tumor cell subtypes ( Figures 3C, D ). In MIA, we found that ATP5E, RPL11, SPINK5, MALAT1, and XIST were highly expressed. In AIS, MALAT1, WSB1, VMP1, EMP2, and NEAT1 exhibited high expression levels. In IAC, ATF3, AREG, C16orf89, ATP5E, and RPL11 showed elevated expression. The C0 subtype exhibited high expression of CXCL8, CXCL2, RRAD, ATF3, and AREG. The C1 subtype showed elevated levels of RPL11, ATP5E, GSTA1, C16orf89, and SPINK5. The C2 subtype had increased expression of NEAT1, MALAT1, VMP1, WSB1, and XIST. The C3 subtype displayed high expression of EMP2, CLIC5, CAV2, CAV1, and AGER. Figure 3E visualized these different genes. We found that JUN, ATF3, IRF1, IER5, and MAFF were upregulated in the C0 subtype. RPS1, RPL11, RPS14, RPS18, and RPL37 were upregulated in the C1 subtype. MALAT1, NEAT1, WSB1, MT-ND3, and VMP1 showed increased expression in the C2 subtype. CAV1, CAV2, AGER, EMP2, and NCKAP5 were upregulated in the C3 subtype. We subsequently analyzed how the top five marker genes were differentially expressed ( Figure 3F ).

Moreover, ratios within separate groups, phases of development, and specific cell subtypes were determined, discovering that the C0 subtype represented the highest proportion within the IAC, accounting for up to 47.2% ( Figures 3G-I ). Therefore, the heterogeneity of the IAC group might be related to the C0 subtype.

Next, we presented the results of nCount RNA, nFeature RNA, and cell stemness AUC for different subtypes using bar graphs ( Figure 3J ). The findings suggested that, in these aspects, C0 demonstrated elevated levels of expression in these metrics relative to the other subtypes. Consequently, we inferred that the C0 subtype likely corresponded to a higher degree of malignancy.

To explore the biological functions of TCs in LUAD, we conducted GSEA analysis across different subtypes. We observed that the C0 subtype was primarily enriched in the response to lipid ( Figure 3K ). Subsequently, we analyzed the biological processes associated with the four subtypes and the C0 subtype was primarily identified as being enriched in fat cell differentiation, the C1 subtype in cytoplasmic translation, the C2 subtype in regulation of RNA splicing, and the C3 subtype in cell-substrate adhesion ( Figure 3L ). With regard to the biological processes, we found that the C0 subtype was associated with differentiation ( Figure 3M ). While showing considerable negative regulation in pathways linked to mitochondrial translation and spliceosomal complex assembly, the GSEA results also revealed that the C0 subtype displayed a significant positive regulation in pathways related with response to lipopolysaccharide and response to molecules of bacterial origin ( Figure 3N ). C0 subtype linked responses to lipid and molecules of bacterial origin, fat cell differentiation, and hemopoiesis. These results greatly affected the viability and growth of TCs, so facilitating the advancement of LUAD.

In our investigation, we examined the expression of genes connected to stemness in the different tumor cell subtypes and determined that the expression levels of ATF3, CXCL2, RRAD, AREG, and CXCL8 were significantly higher in the C0 subtype ( Figure 4A ). UMAP visualization revealed how these genes were distributed ( Figure 4B ). The bar graphs indicated that the expression profiles of ATF3, CXCL2, RRAD, AREG, and CXCL8 in the C0 subtype were expressed at greater levels than in the remaining tumor cell subtypes ( Figure 4C ).

In seeking to understand how LUAD tumor cell subtypes arose and developed, we investigated both their lineage relationships and differentiation states. In Figures 4D-G , we analyzed the results of trajectory and CytoTRACE. CytoTRACE analysis indicated that C0 MAFF+ TCs exhibited a higher degree of stemness ( Figure 4E ). Stemness was closely related to cellular differentiation, higher levels of stemness correspond to lower degrees of differentiation, and TCs in the latter stages of differentiation often possessed greater stemness.

In order to define how the tumor cell subtypes differentiated, we utilized Slingshot analysis studied their developmental trajectories, and further presented the results through UMAP visualization. We first presented two main differentiation trajectories for the four tumor cell subtypes: lineage 1: C2→C1→C0 and lineage 2: C2→C1→C3. The differences between the two trajectories were primarily observed in later stages, where C2 GOLGA8B+ TCs were located at the initial stages of differentiation in both lineage 1 and lineage 2, and C0 MAFF+ TCs were found at the final stage of lineage 1 ( Figure 4H ). Additionally, we analyzed the dynamic trends of marker gene expression among the four subtypes, the high expression of MAFF+ was predominantly observed at later stages ( Figure 4I ), further validating our Monocle analysis findings.

We ultimately conducted analysis of GO-BP enrichment on the subtypes to verify the biological processes pertinent to both lineages ( Figure 4J ). The dynamic timing approach demonstrated how gene expression in TCs shifted along the two pseudotime trajectories.

In order to better comprehend the intricacies of cellular responses, we sought to investigate the networks and intercellular connections that underpin ligand-receptor signaling, which allowed us to depict the interactions among different cell types. By employing analysis of CellChat, we established a comprehensive network of cell-cell interactions covering the majority of cell types, such as C0 MAFF+ TCs, C1 GSTA1+ TCs, C2 GOLGA8B+ TCs, C3 CAV1+ TCs, T and NK cells, macrophages, plasma cells, MCs, ECs, monocytes, fibroblasts, B cells, cDC2, proliferating cells, cDC1, myofibroblasts and pDCs ( Figure 5A ). We found that in the outgoing signaling pattern, C0 MAFF+ TCs predominantly exhibited pattern 1, while JAM, CADM, GDF, OCLN and UGRP1 displayed high expression levels within this pattern. Correspondingly, in the incoming signaling pattern, C0 MAFF+ TCs also primarily represented pattern 1, where SCT, CDH, OCLN, DESMOSOME, and TWEAK showed elevated expression levels ( Figure 5B ). The heatmaps illustrated the communication patterns and signaling pathways of each tumor cell subtype ( Figure 5C ). We found that the degree as well as the count of connections involving C0 MAFF+ TCs and macrophages were more significant in both afferent and efferent signals, so it could be inferred that there was a strong interaction between C0 MAFF+ TCs and macrophages ( Figures 5D, E ).

To present the results more intuitively, we employed a hierarchical graph to depict the relationship between C0 MAFF+ TCs and macrophages. The results demonstrated that the primary interactions between C0 MAFF+ TCs and macrophages occurred through paracrine and autocrine signaling, leading to intense communication ( Figure 5F ). Subsequently, through network centrality analysis of the MIF signaling pathway, we examined the roles of C0 MAFF+ TCs and macrophages within this pathway. The findings revealed that C0 MAFF+ TCs acted as senders, mediators and influencers, whereas macrophages primarily functioned as influencers, senders and mediators, which could have been linked to the conversion of regular macrophages into tumor-associated macrophages (TAMs) ( Figure 5G ). We compared the receptor-ligand interactions between C0 MAFF+ TCs and other cell types, allowing us to infer that the MIF ligand from C0 MAFF+ TCs interacted with the CD74-CD44 receptors on macrophages ( Figures 5H, I ). Further supporting our conclusions was a circular diagram showing the interaction between C0 MAFF+ TCs and macrophages within the MIF-(CD74+CD44) signaling pathway ( Figure 5J ).

Overall, our work exposed the interactions between C0 MAFF+ TCs and macrophages in LUAD, stressing their possible relationship with the change of macrophages to TAMs, so fostering the development of LUAD.

TFs governed gene expression through their ability to bind specific consensus sequences located in the local chromatin setting (54). TFs were typically considered to comprise two main domains: a DNA-binding domain and a functional domain. Additionally, these domains were seen as independent and separable, with the DNA-binding domain responsible for gene targeting and the functional domain for facilitating transcriptional regulation (55).

We employed the SCENIC method to perform two-dimensional clustering of TCs from LUAD with reference to various subtypes and groups ( Figure 6A ). Afterwards, we generated a matrix for the connection specificity index based on the similarity of AUCell scores, which categorized LUAD TCs into four regulatory modules: M1, M2, M3, and M4 ( Figure 6B ). A visualization analysis of these modules was conducted, revealing that C0 MAFF+ TCs exhibited increased expression in the M1 module compared to other subtypes. Furthermore, UMAP visualizations were employed to examine how each tumor cell subtype was distributed among the various modules, indicating that C0 MAFF+ TCs were predominantly located in the M1 module ( Figures 6C, D ).

Afterwards, we evaluated the top five TFs for each distinct subtype of TCs, with a particular focus on their specificity scores across the various modules. Notably, in M1, C0 MAFF+ TCs displayed the highest regulatory activity score, supporting our earlier conclusions ( Figures 6E, F ). Furthermore, we ranked the TFs for each subtype and looked at the distribution of various subtypes ( Figures 6G, H ). In C0 MAFF+ TCs, we examined the distribution and expression levels of five important TFs (KLF6, ATF3, JUN, FOS, and FOSB) in various subtypes. The findings showed that FOS expression was considerably elevated in C0 MAFF+ TCs relative to those in other subtypes ( Figures 6I, J ). It is still unknown, though, exactly how FOS affects LUAD. Consequently, validating the role of FOS in LUAD cells through in vitro functional studies became a principal focus.

We explored the role of FOS in LUAD by performing in vitro assays using A549 and NCI-H1975 cell lines. First, we inhibited FOS and analyzed the related mRNA and expressions of protein at baseline and post-knockdown. Compared to the control group, both cell lines showed substantially decreased mRNA and protein expression ( Figure 7A ). Furthermore, cell viability was greatly diminished after FOS inhibition ( Figure 7B ). We subsequently examined the potential for TCs to proliferate. A significant decrease in cell counts was observed in the colony development assay as a result of the FOS suppression ( Figures 7C, D ). In the EDU experiment, we noted a decrease in colony density following the knockdown, data analysis indicated that siFOS impeded proliferation ( Figures 7E, F ). To evaluate cell migratory and invasive behaviors, we performed assays of transwell and wound healing. Transwell assays depicted that siFOS TCs showed fewer migrated cells and lower cell density compared to the si-NC group, demonstrating that the knockdown of FOS impaired their migration and invasion potential ( Figures 7G-I ). The wound healing assays illustrated that siFOS inhibited wound healing of TCs in both A549 and NCI-H1975 cell lines ( Figures 7J, K ).

Collectively, the results suggested that reducing FOS expression impaired cell activity and significantly decreased the migratory, invasive, and proliferative abilities of TCs. This underlined the important tumor-promoting action of FOS, which was essential in the evolution of TCs. Therefore, a main focus of LUAD treatment is FOS inhibition since it might increase patient survival and prognosis.

We utilized univariate Cox regression analysis to explore prognostic indicators in patients and identified twenty-five genes that exhibited prognosis-related associations. We could observe the HR<1 of SEMA4A, PDE4C and GDF15, while the HR values of the remaining genes were >1 ( Figure 8A ). To reduce potential multicollinearity across genes, we further refined the list of prognosis-related genes ( Figure 8B ). Following this, we performed multivariate Cox regression analysis on twelve selected genes to determine their respective genetic risk coefficients ( Figures 8C, D ). Additionally, our analysis of the curve chart and scatter plot demonstrated that the low MAFF+ tumor risk score (MTRS) group exhibited lower risk scores and better survival outcomes compared to the high MTRS group, which was linked to poorer prognosis ( Figure 8E ). Heatmap further illustrated the distinct prognostic gene expression patterns between these groups, with SEMA4A, PDE4C, and GDF15 more strongly expressed in the low STRS group, suggesting their potential positive impact on prognosis ( Figure 8F ). Moreover, the ROC curves and corresponding AUC values at 1-year, 3-years, and 5-years were 0.73, 0.68 and 0.63 respectively, highlighting the strong predictive power of the model ( Figure 8G ). Kaplan-Meier survival analysis further corroborated the worse prognosis observed in the high MTRS group ( Figure 8H ). In summary, our analysis revealed a negative association between the risk score and survival, supporting the observation that higher risk scores were linked to shorter survival periods ( Figure 8I ).

Further examination of the expression levels of twelve prognostic genes in the high MTRS group and low MTRS group, the scatter plots showed that GDF15 and SEMA4A were beneficial for patient prognosis, while DNAJA1, CEBPB, MAFF, SOX9, HSPA1A, ERRFI1, TM4SF1, KLF4, and IER3were associated with poor patient prognosis ( Figure 8J , Supplementary Figures 2A, B ).

The study compared gene expression levels and underlying biological processes across groups with high and low values, we made use of various visualization and enrichment analysis techniques. We first illustrated the estimated proportions of cell populations in both groups using a stacked bar graph ( Figures 9A, B ). Additionally, we utilized a bubble plot to illustrate the correlations among immune checkpoints, prognostic genes, overall survival, as well as the risk, revealing that KLF4 generally exhibited a positive correlation with immune checkpoint-associated genes, while HSPA1A showed a negative relationship ( Figure 9C ). Upon additional examination, we observed that higher risk scores were correlated with increased macrophage presence ( Figure 9D ). Moreover, the stromal and estimate appeared higher in the high MTRS group, while the immune score was lower compared to the low MTRS group ( Figure 9E ). Through box plot analysis, we discovered that the majority of immune checkpoint genes had considerably increase expression in the high MTRS group as opposed to the low group ( Figure 9F ).

To shed light on the differences between the high and low scoring cohorts, we examined genes that were differentially expressed ( Figure 9G ). The volcano plot further depicted the trends of gene upregulation and downregulation among these DEGs ( Figure 9H ). Following this step, enrichment analysis was undertaken to explore the biological processes associated with the identified genes. We initially performed GO enrichment analysis, which identified significantly enriched terms in biological processes (BP), cellular component (CC), and molecular functions (MF) ( Figure 9I ). In terms of GO-BP, the genes were prominently enriched in pathways involved the organization of intermediate filaments, the structuring of the intermediate filament cytoskeleton, processes based on intermediate filaments, and the differentiation of keratinocytes. The cornified envelope, keratin filaments, intermediate filament cytoskeleton, and intermediate filaments were identified as the main sites of enrichment in terms of GO-CC. GO-MF analysis revealed enrichment included structural constituent of the skin epidermis as well as serine-type endopeptidase inhibitor, endopeptidase inhibitor, and peptidase inhibitor activities.

Afterwards, the heatmap offered a clear visual comparison that further supported our earlier findings ( Figure 10A ). Additionally, we performed KEGG pathway analysis, which highlighted enrichment in pathways such as the estrogen signaling pathway, staphylococcus aureus infection, steroid hormone biosynthesis, and retinol metabolism ( Figure 10B ). To better elucidate the functionality and enrichment of pathway associated with LUAD, further investigations were carried out, we also conducted GSEA analysis ( Figure 10C ). These results demonstrated positive enrichment trends in keratinization, differentiation of keratinocytes, development of the skin, and the epidermis underwent development. Conversely, negative enrichment trends were observed for positive regulation of myoblast fusion, peptide antigen assembly with MHC protein complex, metal ion export, and DNA replication-dependent chromatin organization.

In conclusion, the results of our drug sensitivity analysis indicated that the high MTRS group showed heightened sensitivity to chemotherapeutic agents ( Figure 10D ), including GSK269962A, Midostaurin, SB.216763, WH.4.023, ZM.447439, AICAR, Pazopanib, and QS11. Furthermore, it was determined that the low MTRS group exhibited decreased IC50 values for VX.702 and Roscovitine compared to the high MTRS group. This finding suggested that VX.702 and Roscovitine may have more favorable therapeutic effects for low MTRS group patients.