The transcriptome dataset, along with relevant clinical data, was obtained from 59 control tissue samples and 513 LUAD tumor tissue samples in The Cancer Genome Atlas (TCGA)-LUAD dataset from TCGA database (https://cancergenome.nih.gov/) (26). Additionally, the GSE30219 (GPL570) dataset included 85 LUAD tumor tissue samples and 14 control tissue samples retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30219) (27). The GSE31210 dataset (GPL570) comprised 226 LUAD tumor tissue samples, while the GSE131907 (GPL16791) single-cell dataset included 11 LUAD tumor tissue samples and 11 control tissue samples, also sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131907) (27). A total of 59 PMRGs were retrieved from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Genome-wide association summary statistics for differentially expressed (DE)-PMRGs and 65,864 patients with LUAD were collected from the Genome Wide Association Studies (GWAS) database (https://www.ebi.ac.uk/gwas/). The GWAS data for LUAD (ieu-a-984) included 11,245 patients with LUAD and 54,619 controls with 10,345,176 SNPs, which were used as the summary association statistics for the outcome.

Principal component analysis (PCA) clustering analysis was performed on TCGA-LUAD dataset to investigate clustering patterns between LUAD and normal samples using the plotPCA function from the ‘DEseq2’ (28) R package (V1.38.0; Posit Software, PBC). The primary aim of the present study was to identify differentially expressed genes (DEGs) in LUAD and correlate them with PMRG functional modules. The ‘DEseq2’ R package was utilized to identify DEGs in TCGA-LUAD tumor tissues compared with control tissues (adjusted P<0.05 and |log₂ fold change (FC)| >2) (28). The ‘ggplot2’ package (V3.4.4) (29) and the ‘pheatmap’ package (V1.0.12) (30) were used to generate volcano plots and heatmaps, respectively. To control the false discovery rate (FDR) arising from multiple comparisons, P-values were adjusted using the Benjamini-Hochberg method. DEGs were identified based on the criteria of adjusted P<0.05 and log₂FC >2.

To explore the core gene modules associated with PMRG functions, the single sample gene set enrichment analysis (ssGSEA) method (31) was applied to determine the enrichment fraction of PMRGs (PMRG scores) within individual samples. The PMRG scores were treated as traits in WGCNA performed using the ‘WGCNA’ package in R (V1.71) (32) to identify module genes highly associated with PMRG scores. Initially, hierarchical clustering was applied to the samples to identify and remove any outliers. A soft threshold (β) with connectivity close to 0 was chosen by setting R²>0.85. A scale-free network was constructed based on the selected soft threshold, which divided all genes into multiple modules visually distinguished by different colors (minModuleSize=200; deepSplit=3; mergeCutHeight=0.25). Spearman correlations between these modules and PMRG scores were computed and the genes within the modules that exhibited the highest positive and negative correlations were selected as the module genes associated with PMRG scores.

The core genes (DE-PMRGs) that were differentially expressed and associated with the function of PMRGs were identified through intersection analysis, forming the basis for subsequent causal verification and model construction. DEGs and PMRG module genes were intersected using the ‘VennDiagram’ R package (V1.7.3) (33) to obtain DE-PMRGs. Enrichment analysis of DE-PMRGs was performed using the ‘clusterProfiler’ R package (V4.7.1.003) (34), including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, with a significance level of P<0.05 considered statistically significant. To investigate protein-protein interactions (PPIs) within the DE-PMRG group, the STRING online network tool (https://string-db.org/) was used to examine protein associations, with a significance threshold of P<0.05 (interaction score >0.7). A network diagram was generated and visualized using Cytoscape software (V3.9.1) (35).

To verify the causal relationship between DE-PMRGs and LUAD and eliminate confounding factors, two-sample MR and multivariate MR analyses were conducted using the ‘extract_instruments’ function from the ‘TwoSampleMR’ R package (V0.5.7) (36). In this analysis, DE-PMRGs were the exposure of interest, LUAD was the outcome and SNPs served as instrumental variables (IVs). The MR method relied on the following assumptions: i) IVs are strongly associated with DE-PMRG risk; ii) IVs influence LUAD risk only through their effect on DE-PMRG risk; and iii) IVs are independent of confounders.

GWAS data for expression quantitative trait loci of DE-PMRGs were used as the summary association statistics for the exposure. The exposure factors were assessed, and IVs were screened using the ‘TwoSampleMR’ R package (V0.5.7) (36) function ‘extract_instruments’ to identify IVs with significant associations with the exposure factors (P=5×10⁻⁸). IVs for linkage disequilibrium were removed (clump=TRUE; r²=0.001; kb=100). The ‘extract_outcome_data’ function was used to retrieve the outcome data and eliminate IVs significantly associated with the outcome. Simultaneously, IVs with an F-value <10 were excluded (F=β²_exposure/se²_exposure).

A total of five MR methods, including MR Egger (37), weighted median (38), inverse variance weighted (IVW) (39), simple mode (40) and weighted mode (41), were employed to perform robust causal analyses. A significant causal effect of DE-PMRGs on LUAD was determined using P<0.05 from the IVW method among the five MR approaches. The ‘mr_heterogeneity’ function within the ‘TwoSampleMR’ R package (42) was used to perform a heterogeneity test (Cochran's Q test; P>0.05). The ‘mr_pleiotropy’ function (43) was used to evaluate horizontal pleiotropy (P>0.05). Finally, the ‘directionality_test’ function (44) was used for Steiger analysis, with correct_causal_direction=TRUE and P<0.05 as the criteria for establishing the causal relationship.

To identify genes associated with the prognosis of LUAD, univariate Cox regression analysis was performed on the MR acquisition genes using the ‘survival’ package (V3.8-3), based on LUAD samples with survival data from TCGA-LUAD dataset. This aimed to initially screen genes related to LUAD prognosis (criteria, HR≠1 and adjusted P<0.05). The proportional hazards (PH) assumption test was performed for the genes obtained from univariate Cox regression (P>0.05). In the univariate Cox regression analysis, FDR correction was applied to the P-values to account for multiple testing. Genes with an adjusted P<0.05 and HR=1 were considered candidate prognostic markers for further least absolute shrinkage and selection operator (Lasso) regression analysis. Subsequently, the ‘glmnet’ package (V4.1.4) was used to conduct Lasso analysis on the candidate prognostic genes. At the optimal λ value, the Lasso model minimized the error rate and successfully identified prognostic genes with non-zero regression coefficients, analyzed through 10-fold cross-validation. The risk model was constructed by multiplying the coefficients derived from Lasso regression by the gene expressions: Risk score=Σ[coef(gene_i) × expr(gene_i)] for i=1 to n. Based on their risk scores, individuals were divided into high-risk and low-risk groups, using the optimal cutoff value of −0.5605245. Survival differences between these groups were analyzed using Kaplan-Meier curves via the ‘survminer’ package (V0.4.9) (https://CRAN.R-project.org/package=survminer). Additionally, reciever operating characteristic (ROC) curves were generated using the ‘survivalROC’ package (V1.0.3) (https://CRAN.R-project.org/package=survivalROC). The validity of the risk model was further assessed by replicating the analysis in the GSE30219 dataset. The differential expression of prognostic genes between LUAD tumor and control tissues was examined using Wilcoxon tests (P<0.05) on TCGA-LUAD, GSE30219 and GSE31210 datasets, visualized with box plots generated using the ‘ggplot2’ package (V3.4.4) (29).

To enhance the clinical applicability of the model, a nomogram was constructed by integrating clinical features and risk scores. Univariate and multivariate Cox regression analyses were performed to assess the potential of risk score and various clinical characteristics [including sex, age, stage, tumor (T), lymph node (N) and metastasis (M) staging] as independent predictors in patients with LUAD (criteria, HR≠1 and P<0.05). The ‘rms’ R package (V6.5-0) (45) was used to generate the nomogram, predicting 3-, 5- and 7-year overall survival (OS) based on clinical characteristics and risk score. Calibration curves [‘rms’ R package V6.5-0 (45)] and decision curves (‘ggDCA’ R package V1.6.0) (46) were generated to assess the performance of the nomogram.

To explore the relationship between clinical characteristics and risk scores, the Wilcoxon test was used to examine differences in risk scores across various subgroups based on clinical features within TCGA-LUAD dataset (P<0.05).

To examine the molecular localization of prognostic genes and provide insights into their functional mechanisms, the distribution of prognostic genes on chromosomes was analyzed using the ‘Circos’ R package (V1.2.2) (47). Subcellular localization was also assessed to identify the specific cellular regions where the prognostic genes are located. mRNA sequences (FASTA format files) of prognostic genes were obtained from the National Center for Biotechnology Information's GENE database (https://www.ncbi.nlm.nih.gov/gene/?term=) and analyzed for subcellular localization using the mRNALocater online database (http://bio-bigdata.cn/). The results were visualized by plotting histograms using the ‘ggplot2’ R package (V3.4.4) (29).

To expand the functional association network of prognostic genes and identify potential synergistic genes, GeneMANIA (http://www.genemania.org/) was employed to analyze genes related to the functions of prognostic genes. GeneMANIA was used to predict genes associated with the functional roles of prognostic genes and their corresponding functions.

The functional similarity of prognostic genes was assessed using GO term-based Friends analysis with the ‘GOSemSim’ R package (V2.24.0) (48). Similarity scores were averaged, ranked from high to low and the results were visualized.

To explore the upstream regulatory mechanisms of prognostic genes, a multi-level regulatory network was constructed. The miRDB (http://mirdb.org) and TargetScan (https://targetscan.org) databases, accessed through the ‘multiMiR’ R package (V1.20.0) (49), were used to predict upstream microRNAs (miRNAs/miR) associated with the prognostic genes. The intersection of miRNAs identified from both databases was determined using the ‘VennDiagram’ R package (V1.2.3) (33), resulting in a target miRNA set. The Starbase database (https://rnasysu.com/encori/) was then used to predict upstream long non-coding RNAs (lncRNAs) for the target miRNAs. Additionally, the hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget) was utilized to predict transcription factors (TFs) regulating the prognostic genes. The lncRNA-miRNA-prognostic gene and TF-prognostic gene-miRNA regulatory networks were subsequently constructed using Cytoscape software.

To examine the functional differences between high- and low-risk groups and uncover the potential biological mechanisms underlying the model, GSEA was performed for GO and KEGG pathways in TCGA-LUAD dataset, comparing samples with high and low risk scores. Pathways were considered enriched if they met the criteria of normalized enrichment score >1 and P<0.05, with the top five most significantly enriched pathways ranked by their P-values.

To explore the association between risk models and the tumor immune microenvironment and provide a basis for immunotherapy, the present study utilized the ssGSEA algorithm in the training set and performed the Wilcoxon test to compare differences between high- and low-risk groups (P<0.05). Additionally, a correlation analysis was conducted in TCGA-LUAD dataset using Spearman's rank correlation coefficient to examine the relationship between prognostic genes and differential immune cell populations, with a significance threshold of Ρ>0.3 and P<0.05.

To identify individualized therapeutic drugs based on risk models and enhance their clinical applicability, the present study screened for candidate drugs in the Genomics of Drug Sensitivity in Cancer 2 database (https://www.cancerrxgene.org/), prioritizing those associated with risk scores. Half-maximal inhibitory concentration (IC₅₀) values for various drugs were calculated for both risk groups using the ‘oncoPredict’ tool (V1.2.0) (50). Spearman's analysis, conducted using the R package ‘psych’ (V2.3.6) (https://CRAN.R-project.org/package=psych), was used to assess the correlation between drug IC₅₀ values and risk scores. The Wilcoxon test was then employed to compare IC₅₀ expression for common chemotherapeutic agents between high- and low-risk groups, with statistical significance set at P<0.05.

To analyze the cellular expression characteristics and key cell populations of prognostic genes at the single-cell level, the ‘seurat’ R package (V5.0.1) was used to process single-cell data from GSE131907, which were converted into Seurat objects (51). Quality control was applied to select high-quality cells, characterized by a gene count between 300 and 10,000 per cell and a mitochondrial proportion <10%. The LogNormalize method (52) was applied to normalize feature expression measurements of each cell against total expression. These normalized values were scaled by a factor of 10,000 and logarithmically transformed. The FindVariableFeatures function was used to identify highly variable genes (HVGs). PCA was performed to reduce the dimensionality of these genes, and relevant principal components (PCs) were selected for cluster analysis using the Uniform Manifold Approximation and Projection (UMAP) method (53). Significant marker genes within clusters were identified using the Find All Markers function, which facilitated cell type identification based on marker expression (Table SI). These marker genes were compared with those reported in the literature (27) to identify distinct cell clusters. To pinpoint key cell clusters, differential expression of prognostic genes between annotated LUAD and control tissue cells was compared using the Wilcoxon test (P<0.05), with significant genes defining the key clusters.

To elucidate the biological pathways through which key cell clusters function, functional enrichment analysis was performed on key cell clusters from LUAD samples using the R package ‘ReactomeGSA’ (V4.3.0) (54) to explore signaling pathways associated with LUAD progression.

Cell communication among annotated cells was analyzed in CellPhoneDB (http://www.cellchat.org/), using the CellChatDB.human database as a reference. The ‘CellChat’ package (V1.6.1) (55) was utilized to investigate cell-to-cell interactions (P<0.05). The single-cell data, pre-processed and annotated via Seurat, served as the input. Based on the built-in human ligand-receptor database of CellChat, significantly overexpressed ligands or receptors were first identified across different cell populations. Communication probabilities between cell types were then calculated, and significant interactions were determined through 1,000 permutation tests. Subsequent analyses involved signaling pathway mapping, identification of communication patterns and network centrality evaluation to pinpoint key interaction axes. The core results were visualized for further interpretation. To uncover the developmental trajectories and interactions of key cell clusters, pseudotime analysis was conducted using the Monocle algorithm (V2.30.0) (56). This approach aimed to investigate the differentiation direction and explore the effect of prognostic gene expression on the differentiation degree of key cell clusters.

To validate the expression patterns of prognostic genes, clinical sample experiments were conducted to enhance the reliability of the research. Differential expression between LUAD and control samples was analyzed using the Wilcoxon test (P<0.05) in TCGA and GSE30219 datasets. To further validate these findings, 10 clinical samples (5 LUAD and 5 control samples) were prospectively recruited between December 2022 and June 2023 from patients undergoing radical surgery at Hunan Cancer Hospital (Changsha, China). All participants provided written informed consent prior to sample collection. The study protocol was approved by the Ethics Committee of Hunan Cancer Hospital (approval no. SBQLL-2022-127). All clinical procedures were performed in strict accordance with the approved guidelines and the Declaration of Helsinki. The inclusion criteria for patients with LUAD in the present study were: i) Pathologically confirmed LUAD; ii) patients who underwent radical surgery; and iii) no prior history of radiotherapy or chemotherapy before the operation. The exclusion criteria included: i) Presence of other primary malignancies; ii) secondary lung cancer; and iii) samples with poor RNA integrity or microbial contamination. Regarding the control group, the control samples were obtained from adjacent healthy lung tissues (located ≥5 cm away from the tumor margin) from the same patients with LUAD, rather than from separate healthy individuals.

RT-qPCR was used to assess gene expression. Total RNA was extracted using the Trizol kit (Ambion; Thermo Fisher Scientific, Inc.; cat. no. 15596-018CN), and cDNA was synthesized via reverse transcription using the SureScript First Strand cDNA Synthesis Kit (Wuhan Servicebio Technology Co., Ltd). Each RT-qPCR reaction contained 3 µl cDNA, 5 µl 2X Universal Blue SYBR Green qPCR Master Mix (Wuhan Servicebio Technology Co., Ltd) and 1 µl each forward and reverse primer (10 µm; Table SII). The amplification protocol began with 95°C for 1 min, followed by 40 cycles of denaturation at 95°C for 20 sec, annealing at 55°C for 20 sec and extension at 72°C for 30 sec. The results were analyzed using the 2^−ΔΔCq method (57) with GAPDH as the internal reference gene (Table SIII). Gene expression differences between LUAD and control samples were compared using GraphPad Prism 5 (Dotmatics; P<0.05).

Statistical analyses were performed using R (V4.2.2; Posit Software, PBC), with group comparisons conducted using the Wilcoxon test. Statistical significance was defined as P<0.05, unless otherwise specified. For high-throughput data analyses, including differential expression and initial prognostic gene screening, multiple testing corrections (FDR/Benjamini-Hochberg method) were systematically applied to minimize false-positive results. A P-value or adjusted P-value of <0.05 was considered to indicate a statistically significant difference.

PCA was performed on TCGA-LUAD dataset to characterize the distribution patterns of LUAD samples vs. normal samples. The analysis revealed a noticeable separation between the two groups in the PCA space (Fig. S1). To identify polyamine-related module genes in LUAD, differential analysis of PMRG scores was conducted between LUAD and control groups, showing significant differences (P<0.05; Fig. 1A). Following WGCNA (β=9; R²=0.85) and module-gene screening (Fig. 1B and C), the MEblue module showed the lowest correlation with PMRG scores (R=−0.72; P<0.05), while the MEmagenta module exhibited the highest correlation (R=0.70; P<0.05). These two modules identified 2,882 module genes (Fig. 1D and E). Using TCGA-LUAD data, 2,253 DEGs were found between LUAD and control samples, with 1,684 upregulated and 569 downregulated (Fig. 2A and B). By intersecting DEGs with module genes, 470 DE-PMRGs were identified (Fig. 2C).

The DE-PMRGs were involved in 1,399 GO terms, including 1,158 biological process (BP), 114 cellular component (CC) and 127 molecular function (MF) terms, as well as 45 KEGG pathways (Fig. 2D and E). Specifically, these were associated with ‘ameboidal-type cell migration’ and ‘vascular processes in the circulatory system in BP’, ‘collagen-containing extracellular matrix’ and ‘membrane rafts’ in CC, ‘signaling receptor activator activity’ and ‘receptor-ligand activity’ in MF and ‘vascular smooth muscle contraction’ in KEGG pathways. The PPI network analysis revealed extensive protein interactions among most DE-PMRGs, suggesting a synergistic regulation of LUAD malignant phenotypes (Fig. 2F).

MR analysis primarily based on the results from the IVW algorithm identified 30 DE-PMRGs significantly causally associated with LUAD (P<0.05; Fig. 3; Table SIV). Among these, 17 genes, including activin A receptor like type 1 and adenoreceptor β-1 (ADRB1), were protective factors for LUAD occurrence [odds ratio (OR) <1], indicating their potential to inhibit tumorigenesis and progression via polyamine metabolism regulation. A total of 13 genes, including carbonic anhydrase 4 (CA4) and caveolin 2, were identified as risk factors (OR >1), with their abnormal expression potentially promoting LUAD malignant transformation by disrupting polyamine metabolism. Correlation scatterplots of SNP-exposure factor effects and SNP-outcome effects confirmed the high consistency between the present study results and the IVW findings (Fig. S2). The forest plots further corroborated the consistency of the IVW results (Fig. S3, Fig. S4, Fig. S5, Fig. S6, Fig. S7). MR randomness testing indicated that the 30 candidate genes adhered to Mendel's second law (Fig. S8). The heterogeneity test (Q-P-value >0.05) indicated no heterogeneity among the samples (Table SV). The horizontal pleiotropy test showed P>0.05, suggesting no significant confounding effects (Table SVI). SNP elimination tests further confirmed the reliability and stability of the MR results (Fig. S9, Fig. S10, Fig. S11, Fig. S12, Fig. S13). Steiger filtering analysis of the SNPs from the 30 selected candidate genes revealed that all SNPs had a TRUE direction, excluding reverse causality interference. This preliminary analysis supports the causal association between these DE-PMRGs and LUAD (Table SVII).

Univariate Cox regression and the PH assumption test identified four genes [calcium voltage-gated channel auxiliary subunit α2δ2 (CACNA2D2), ADRB1, immunoglobulin superfamily member 10 (IGSF10) and CA4] as prognostic markers (Fig. 4A; Table I). These four genes were selected at λ=0.01362608 [log(λ)=−4.29577] under the condition of the lowest evaluated error rate, indicating their retention in the model (Fig. 4B and C). A risk score model was constructed using gene expression and regression coefficients (Table II). LUAD samples from TCGA dataset were categorized into high-risk and low-risk groups based on the optimal risk score cutoff (Fig. 4D and E). Kaplan-Meier survival analysis revealed a significantly lower survival rate in the high-risk group compared with the low-risk group (P<0.0001; Fig. 4F). ROC curve analysis demonstrated that the model's area under the curve (AUC) values for predicting 3-, 5- and 7-year survival rates of patients with LUAD were all >0.6, suggesting modest but consistent predictive performance (Fig. 4G). The prognostic gene heatmap showed significant downregulation of IGSF10, ADRB1, CA4 and CACNA2D2 in LUAD samples (Fig. 4H).

Using the same methodology, LUAD samples in the GSE30219 and GSE31210 datasets were similarly divided into high-risk and low-risk groups (Fig. 5A and B). The survival time of patients in the high-risk group was significantly shorter than that in the low-risk group (Fig. 5C). ROC curve analysis also showed that the AUC values for predicting 3-, 5- and 7-year survival rates were >0.6, further supporting the stability of the model (Fig. 5D). These results demonstrate that the prognostic model constructed in the present study exhibits high accuracy and stability across both datasets, affirming its clinical value for predicting the outcomes of patients with LUAD.

Univariate Cox regression identified seven factors (including risk score, age and stage) that correlated with prognosis (P<0.05; Fig. 6A). The PH assumption test (P>0.05) and multivariate analysis confirmed risk score, T stage, M stage and N stage as independent prognostic factors (P<0.05; Fig. 6B, Table III). A nomogram was developed incorporating these independent prognostic factors (Fig. 6C), allowing for the prediction of survival outcomes for patients with LUAD at 3-, 5- and 7 years. The prognostic model demonstrated certain diagnostic potential in predicting the probabilities of 3-, 5- and 7-year survival (Fig. 6D). Decision curve analysis indicated that the net benefit curve corresponding to the nomogram remained a significant distance from the extreme curves (‘All’ and ‘None’) for the 3-, 5- and 7-year assessments. Furthermore, the net benefit of the nonogram surpassed that of individual variables such as T stage, N stage and M stage across the entire threshold probability range, suggesting that the model provides superior net benefits for clinical decision-making in the individualized prognostic assessment of patients with LUAD (Fig. 6E and F). Variations in risk assessment within subgroups defined by clinical features in TCGA-LUAD revealed significant disparities in risk scores across subgroups based on age, T stage, N stage and overall clinical stage (P<0.05; Fig. 6G).

The GeneMANIA database identified 20 genes associated with the functions of prognostic genes, including solute carrier family 4 (SLC4) member 4, SLC4 member 1 and alkaline ceramidase 2. These prognostic genes were predicted to be involved in 105 functions, such as peptide hormone secretion, insulin secretion, regulation of hormone secretion and voltage-gated cation channel activity, among others (Fig. 7A). The functional similarity of prognostic genes was ranked from highest to lowest based on the average similarity scores. Results indicated that IGSF10, ADRB1, CA4 and CACNA2D2 exhibited low mean functional similarity (threshold >0.5; Fig. 7B).

Chromosomal localization revealed that ACNA2D2 and IGSF10 were primarily located on chromosome 2, ADRB1 on chromosome 10 and CA4 on chromosome 17 (Fig. 7C). Subcellular localization analysis showed that ADRB1, IGSF10 and CA4 were predominantly located in the cytoplasm, whereas CACNA2D2 was primarily localized in the endoplasmic reticulum (Fig. 7D).

The miRDB and TargetScan databases predicted upstream miRNAs for the prognostic genes, identifying 97 and 111 miRNAs, respectively. A total of 26 target miRNAs were common between the two databases. The upstream lncRNAs of these target miRNAs were predicted using the Starbase database, resulting in the identification of 73 lncRNAs. The lncRNA-miRNA-mRNA regulatory network was then constructed (Fig. 7E). Notably, lncRNAs XIST and NEAT1 co-regulated OLFML2B via hsa-miR-18a-5p. Additionally, 37 TFs (such as FOXA2, histone deacetylase 2 and TATA-box binding protein associated factor 3) were predicted to regulate IGSF10, ADRB1 and CACNA2D2 in the hTFtarget database. Both IGSF10 and ADRB1 were found to be simultaneously regulated by CEBPB in the miRNA-mRNA-TF network (Fig. 7F).

At a significance threshold of P<0.05, the high-risk group was significantly enriched for 756 GO pathways, including structural components of ribosomes, ribosomal subunits, complex-containing mitochondrial proteins and oxidative phosphorylation (Fig. 8A). In KEGG pathway analysis, 29 pathways were significantly enriched in the high-risk group, such as proteasome, ribosome and pentose phosphate pathways, while the low-risk group was enriched in 11 KEGG pathways, including ‘α-linolenic acid metabolism’ and ‘taste transduction’ (Fig. 8B).

Immune infiltration analysis revealed significant variations in 15 immune cell types between high- and low-risk groups, including CD56⁺ natural killer cells, memory B cells and T follicular helper cells (Fig. 8C). Additionally, eosinophils showed a strong positive correlation with all prognostic genes (Fig. 8D).

A total of 33 drugs were identified with significant correlations between their IC₅₀ values and risk scores (|cor| >0.3; P<0.05), including CCT007093, Parthenolide, CMK, RO-3306 and A-443654 (Fig. 8E). Among these, the IC₅₀ values of 15 drugs showed significant differences between the high- and low-risk groups, with drugs such as Parthenolide, CMK and RO-3306 exhibiting lower IC₅₀ values in the high-risk group (P<0.05; Fig. 8F).

After filtering, the GSE131907 dataset retained 83,883 cells and 25,498 genes (Fig. S14A). The data were standardized, and 2,000 HVGs were extracted (Fig. S14B). Clustering analysis of the top 20 PCs identified 21 distinct cell clusters (Fig. S14C and E). A total of eight cell types were annotated based on marker genes, including epithelial cells, fibroblasts, endothelial cells, T lymphocytes, natural killer cells, B lymphocytes, myeloid cells and mast cells (Fig. 9A). The edited genes are highly expressed in all cell types (Fig. S14D). Cell proportion analysis in both LUAD and control groups showed that T lymphocytes and myeloid cells represented relatively high proportions in both groups (Fig. 9B). Significant differential expression of prognostic genes was observed in epithelial cells between LUAD and control tissues, leading to the selection of epithelial cells as a key cell cluster for subsequent analysis (P<0.01; Fig. 9C). Enriched signaling pathways in epithelial cells included TWIK releasing acid-sensitive K⁺ channels (TASK), ATP-sensitive potassium channels and the binding of phenylacetate to glutamine (Fig. 9D). Pseudotime analysis revealed that epithelial cells differentiated across three distinct stages and 19 clusters. Cluster 19 was predominantly localized in stage 2, while cluster 12 was primarily found in stages 1 and 3, with cluster 9 predominantly observed in stage 3 (Fig. 9E). In pseudotime, the expression of ADRB1 and IGSF10 in epithelial cells tended to increase, while the expression of CACNA2D2 initially increased and then decreased (Fig. 9F). Cell communication analysis highlighted frequent and intense interactions between epithelial and myeloid cells, suggesting a synergistic role for these two cell types in tumor microenvironment remodeling and immune regulation (Fig. 9G).

Wilcoxon test analysis revealed notable differences in the expression levels of the four prognostic genes between LUAD and control samples, with notable downregulation in LUAD samples (Fig. 10A and B). Consistent with these findings, RT-qPCR results (Fig. 10C-F) corroborated the dataset outcomes, indicating that these four genes could potentially improve the prognosis of patients with LUAD. To validate the expression levels of the four prognostic genes, RT-qPCR was performed on clinical samples. As shown in Fig. 10C-F, the mRNA expression levels of CACNA2D2, ADRB1, IGSF10 and CA4 were significantly lower in LUAD tissues compared with control tissues (all P<0.01), consistent with the bioinformatic predictions.