The Cancer Genome Atlas (TCGA) database, jointly established by the U.S. National Cancer Institute (NCI) and National Human Genome Research Institute (NHGRI), provides multi-omics data including transcriptomic expression profiles, genomic variation data, and clinical annotations. Clinical datasets were retrieved from the TCGA portal (http://portal.gdc.cancer.gov/), and 562 samples meeting inclusion criteria (complete clinical metadata) were retained for subsequent analyses. All datasets utilized in this study were derived from public repositories, thus exempting the requirement for ethics committee approval.

Differentially expressed genes (DEGs) in lung adenocarcinoma were identified using the DESeq2 package in R (version 4.3.0), with thresholds set at an adjusted p-value (Benjamini-Hochberg false discovery rate [FDR]) < 0.05 and absolute log2-fold change (|log2FC|) > 1. Visualization of DEGs was performed using ggplot2 (volcano plots) and heatmaps (hierarchical clustering heatmaps). Subsequently, gene expression profiles were integrated with survival data (overall survival status and time). Prognostically significant genes were preliminarily screened through univariate Cox regression analysis (threshold: p < 0.05). Functional annotation of these prognosis-related genes was conducted using clusterProfiler and enrichplot packages, including: Gene Ontology (GO) enrichment analysis (biological processes, molecular functions, cellular components), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Results were visualized via dot plots, enrichment maps, and circular dendrograms following best practices for omics data visualization.

The GeneMANIA platform (http://genemania.org/), a widely used online tool for predicting protein-protein interactions and functional associations, was employed to construct a PPI network. The PPI network of these mitochondrial-related prognostic genes was subsequently constructed using the GeneMANIA database.

Differentially expressed genes (DEGs) were integrated with LUAD patient survival data, and samples with incomplete clinical information were excluded. The Glmnet package in R was utilized to perform LASSO (Least Absolute Shrinkage and Selection Operator) regression, followed by Cox proportional hazards modeling for survival analysis. Forest plots were generated using the forestplot package. To mitigate overfitting in the prognostic model, genes identified through univariate Cox regression were subjected to LASSO regression via the Glmnet package, excluding genes with regression coefficients of zero. Subsequently, multivariate Cox regression analysis was performed to identify genes significantly associated with prognosis. A risk score was calculated for each patient based on the final gene set.

The final prognostic risk score model was constructed using the following parameters:

A total of 503 LUAD samples were stratified into high and low risk groups based on the optimal cut-off value determined by receiver operating characteristic (ROC) curve analysis. The pheatmap package was utilized to generate: a survival status heatmap illustrated the distribution of high and low risk groups. An expression heatmap of prognosis-associated genes were used in the model. Survival curves were plotted using the survival package, and the independent prognostic value of the risk score was validated through Cox regression analysis. The predictive performance of the model was further evaluated by constructing ROC curves using an internal validation dataset.

Samples were collected from patients with NSCLC in the Department of Thoracic Surgery, the First Affiliated Hospital of Ningbo University. All patients provided written informed consent. NSCLC tissues and adjacent tissues were collected for research purposes. The samples were sectioned and stored frozen in liquid nitrogen at -80 °C. All procedures involving human participants were conducted in accordance with the 1964 Declaration of Helsinki and its subsequent amendments.

Proteins were extracted from tumor tissues and separated by polyacrylamide gel electrophoresis (Beyotime, Shanghai, China). The separated proteins were transferred to a polyvinylidene fluoride (PVDF) membrane (Invitrogen, USA) using an electro-transfer device (Servicebio, Wuhan, China), which was placed in an ice bath. Blocking was performed with 5% bovine serum albumin (BSA; Beyotime, Shanghai, China), and the membrane was placed on a shaker at a slow speed for 1 hour. The primary antibody was then added for incubation to allow specific binding, followed by incubation with horseradish peroxidase (HRP)-labeled secondary antibody to form an HRP-primary antibody conjugate (1:10000, Abbkine, China). The details of the primary antibodies are as follows: GREB1L (1:1000, Novus, USA), MFI2 (1:1000, CST, USA), SPRR1B (1:1000, Zen-bio, China), GPR115 (1:1000, Baijia, China), LIPK (1:1000, Novus, USA), and β-actin (1:1000, Proteintech, USA). Enhanced chemiluminescence (ECL; Millipore, USA) was used for signal development, which was stopped once clear bands appeared. The film was scanned, and the gray values of the target bands were analyzed using ImageJ software.

Total RNA was extracted from NSCLC cells using TRIzol reagent (Invitrogen, USA), and the RNA concentration was measured. qRT-PCR was performed using SYBR Green (Tiangen, China) on an ABI Illumina instrument (Foster, USA). The information of primers was listed in Table 1 .

All statistical analyses and presentations were performed using the R 4.3.0 software package. The differential expression of immune-related genes in LUAD and normal tissues was compared using Wilcox test. Multivariate Cox regression analysis was used to identify immune-related genes associated with poor prognosis in LUAD. T test was used to evaluate the correlation between prognostic genes and transcription factors. Survival analysis was performed using Kaplan-Meier curves. P<0.05, P<0.01 and P<0.001 was considered a statistically significant difference.

Analysis of LUAD patient data from the TCGA database identified 1,251 upregulated genes and 285 downregulated genes ( Figures 1A, B ). Using stringent thresholds (adjusted p-value [FDR] < 0.05, |log2 fold change (FC)| > 1) and according to high and low expression of CD274 (PD-L1), we detected DEGs comprising 873 upregulated and 1010 downregulated genes. A heatmap depicting CD274 (PD-L1) expression stratification was shown in Figure 1C . Subsequently, Volcano plot demonstrated top 20 differential CD274-related genes (CRGs, Figure 1D ). Upregulated genes included CD274, CXCL11, TBX21, CXCL10, GBP5, GBP1, WARS, SAMD9L, PDCD1LG2 and LILRB2. Downregulated genes were included CBR1, NPAS3, MYCN, CALCB, UGT2B4, PCSK2, KLK12, PGC, INSM1 and PCSK1. To elucidate the relationship between CD274 and LUAD-associated genes, Venn diagram analysis was performed ( Figures 1E, F ). 873 genes were overlapped with 1,251 LUAD-upregulated genes and the intersected gene was 110. 1,010 CD274-downregulated genes were overlapped with 285 LUAD-downregulated genes and the intersected gene was 9.

In order to further screen out valid genes, survival analysis was performed to 110 gene and p < 0.01 was screening criteria. Furthermore, Kaplan-Meier survival curves demonstrated that patients with high-expression of GPR115, MF12, GREB1L, SPRR1B, and LIPK exhibited significantly poorer prognosis compared to low-expression ( Figures 2A-F ). What’s more, comparative analysis revealed GPR115, MF12, GREB1L, SPRR1B, and LIPK was upregulated in LUAD samples ( Figures 2G ), which confirmed by hierarchical clustering heatmaps ( Figure 2H ). Intriguingly, heatmaps revealed a progressive increased in expression levels of GPR115, MFI2, GREB1L, SPRR1B, and LIPK in CD274-high expression ( Figure 2H ), suggesting potential co-regulatory mechanisms between these genes and immune checkpoint pathways.

The malignant progression and prognosis of LUAD were related to a variety of genes and factors. Therefore, this section would like to use Univariate Cox regression to demonstrate whether the five genes were associated with the prognosis of LUAD. Univariate Cox regression analysis of SPRR1B, MFI2, LIPK, GREB1L and GPR115 revealed hazard ratios (HR) > 1 ( Figure 3A ), indicating that high expression of these genes was associated with increased mortality risk. Furthermore, to elucidate functional relationships among these genes and their interactome, network analysis was performed ( Figures 3B ), revealing distinct interaction patterns. Co-expression and co-localization of SPRR1B was associated with keratin family genes KRT14, KRT6B and lymphocyte antigen LY6D. MFI2 was physically interacted with endocytic trafficking regulators. GREB1L was shared protein domains with GREB1 (retinoic acid receptor coactivator). Genetic interactions of GREB1L were with transcription factor TF and TMEM67. GPR115 was functionally associations with oncogenic regulators JUN, integrin ITGA9, and epithelial marker KRT5. Finally, we performed KEGG enrichment analysis on DRGs between CD274-high and low expression samples ( Figures 3C ). The results of KEGG analysis showed that the genes with high expression of CD274 in LUAD samples were mostly concentrated in Cytokine-cytokine receptor interaction, IL-17 signaling pathway and other pathways. The genes with low expression of CD274 were concentrated in the pathways of Neuroactive ligand-receptor interaction, Metabolism of xenobiotics by cytochrome P450 and Arachidonic acid metabolism.

We performed GO functional enrichment analysis ( Figures 4A ) and GSEA analysis ( Figures 4B ) to differential genes between CD274-high and low samples. GO analysis showed that most of the differentially expressed genes were mainly enriched in leukocyte-mediated immune response and immune regulation, receptor ligand activity and other processes. These processes played an indispensable role in the occurrence and development of tumors. GSEA analysis showed that the genes with high expression of CD274 and high expression in tumor patients were enriched in Necroptosis, JAK-STAT signaling pathway and Cell adhesion molecules. The genes with low expression of CD274 and low expression in tumor patients were enriched in Neuroactive ligand-receptor interaction, Hormone signaling and cAMP signaling pathway.

Risk scoring models were a core tool for bridging clinical data and practice by quantifying risk, simplifying decision-making, and facilitating precision medicine. 503 tumor samples were divided into high and low risk groups based on risk scores ( Figures 5A, B ). Survival curves indicated that the high-risk group had a significantly shorter overall survival (OS) than the low-risk group (P = 0.0015). Univariate Cox regression analysis was performed for different influencing factors ( Figure 5C ). Further multivariate Cox regression analysis identified five genes significantly associated with LUAD, which were used to build a prognostic risk scoring model. They were GREB1L (HR = 0.2124, coef = 1.218e - 04), MFI2 (HR < 0.001, coef = 8.705e - 05), SPRR1B (HR = 0.6224, coef = 1.773e - 05), GPR115 (HR = 0.7805, coef = - 3.601e - 05), and LIPK (HR = 0.395, coef = 1.453e - 03). For GREB1L, MFI2, RSPR1B, GPR115 and LIPK, a hazard ratio (HR) > 1 means high expression was linked to high risk ( Figure 5D ). Further multivariate Cox regression analysis was done by combining with other clinical parameters (such as age, gender, and tumor stage) ( Figure 5E ). Results showed that age (HR = 1.011, 95% CI = 0.9954 - 1.027, P = 0.1708), gender (male HR = 3, female HR = 0.992, 95% CI = 0.7353 - 1.339, P = 0.9585), and tumor stage (i stage HR = 0.748, 95% CI = 0.0658 - 8.497, P = 0.8147; ia stage HR = 0.971, 95% CI = 0.2298 - 4.104, P = 0.9682; ib stage HR = 1.161, 95% CI = 0.2780 - 4.848, P = 0.8378; ii stage HR = 5.894, 95% CI = 0.5262 - 66.031, P = 0.1501; iia stage HR = 3.046, 95% CI = 0.7065 - 13.130, P = 0.1352; iib stage HR = 2.157, 95% CI = 0.5094 - 9.134, P = 0.2965; iiia stage HR = 3.440, 95% CI = 0.8207 - 14.421, P = 0.0911; iiib stage HR = 2.430, 95% CI = 0.4852 - 12.169, P = 0.2801; iv stage HR = 3.848, 95% CI = 0.8722 - 16.975, P = 0.0752) were not associated with prognosis in lung adenocarcinoma patients. However, the risk model factor (HR = 0.661, 95% CI = 0.4878 - 0.895, P = 0.0074) was associated with prognosis in lung adenocarcinoma patients. These results confirm the predictive role of the risk score in lung adenocarcinoma prognosis. Importantly, a nomogram was built by combining risk scores with clinical characteristics, with the risk score making a significant contribution to the predictive model ( Figure 5F ).

The risk scores of each patient were calculated according to the risk scoring formula. By drawing a heat map ( Figure 6A ), it was found that the risk scores of above five genes (GPR115, MFI2, GREB1L, SPRR1B and LIPK) were high when they were highly expressed, indicating that they were more likely to have poor prognosis when they were highly expressed. Heatmap showed the differential genes between high-risk and low-risk groups ( Figure 6A ). The volcanic map ( Figure 6B ) also showed that top highest genes such as UPK1B, LHX1, GREB1L, PADI1, A2ML1, TGM5, IL1A, CD109, CCNE1 and MFI2 were found. Top lowest genes were CBR1, CACNA2D2, TMED6, WIF1, SLC38A8, GKN2, SLC14A2, PCSK2, PGC and CALCA. We performed KEGG enrichment analysis on DEGs ( Figure 6C ). Staphylococcus aureus infection, Estrogen signaling pathway, Cytokine-cytokine receptor interaction and Phtotransduction so on were involved positively in high-risk groups. Neuroactive ligand-receptor interaction, Metabolism of xenobiotics by cytochrome P450, and Drug metabolism - cytochrome P450 so on were involved negatively in high-risk groups. GO analysis showed these genes were mainly enriched in channel activity, passive transmembrane transporter activity, monoatomic ion channel activity and so on ( Figure 7A ). GSEA indicated they were predominantly involved in cell cycle, NOD-like receptor signaling pathway, cellular senescence and so on ( Figure 7B ).

Infiltrating immune cells in tumors are able to directly kill tumor cells and are associated with a good prognosis. Natural killer cells (NK cells) induce apoptosis of tumor cells by releasing perforin and granzyme. However, regulatory T cells (Tregs) are able to suppress the immune response and promote tumor immune escape. M2 tumor-associated macrophages secrete proangiogenic factors to support tumor growth and metastasis (43). To explore the correlations among immune cells, we performed pairwise correlation analysis of immune cell gene expression levels and visualized the results with a heatmap using the “ggplot2” package ( Figure 8A ). There were extremely strong correlations between Effector memory CD4 T cells and Type 2 T helper cells, between Effector memory CD8 T cells and MDSCs/regulatory T cells, and between natural killer T cells and Type 1 helper cells. To further evaluate the infiltration levels of immune cells between the high- and low-risk groups, we generated a heatmap ( Figure 8B ). To further assess the infiltration levels of immune cells between high- and low- risk groups, we used the CIBERSORT algorithm to calculate the relative proportions of 28 immune cell subtypes in each lung adenocarcinoma sample. Then, we compared the immune cells with significant differences between the two groups ( Figure 8C ). The high group had higher proportions of Activated CD4 T cells, Effector memory CD8 T cells, Regulatory T cells, and Natural killer T cells, while the low-risk group had higher proportions of Eosinophils and Mast cells.

LASSO regression solves the core problem of high-dimensional data modeling through automatic feature selection and regularization in the prognostic model, making the model more concise, explainable and generalizable. Its combination with survival assays such as Cox LASSO further expands its application value in medical research (44). Therefore, to optimize the model, we used LASSO regression to identify the best genes ( Figure 9A ). With the internal validation cohort, the time-dependent ROC curves revealed 1–3-year AUC values of 0.37-0.687, indicating good predictive performance of the model (Figure 9B).

To confirm the expression of PD-L1 related genes, PCR, Western blotting and immunohistochemistry were used to verify the expression levels of GREB1L, MFI2, SPRR1B, GPR115 and LIPK in LUAD tumor tissues. PCR results showed elevated expression of GREB1L, MFI2, SPRR1B, GPR115 and LIPK mRNA in LUAD tumor tissues ( Figure 10A ). Western blotting results showed that the protein expression levels of GREB1L, MFI2, SPRR1B and LIPK were increased, but the expression levels of GPR115 did not change significantly ( Figure 10B ). Immunohistochemistry confirmed that the expression levels of GREB1L, MFI2, SPRR1B, and LIPK proteins were elevated in LUAD tissues, and GPR115 did not change significantly ( Figure 10C, D ). These results suggest that GREB1L, MFI2, SPRR1B, GPR115 and LIPK might be involved in the pathological process of LUAD.