The detailed research workflow is illustrated in Figure 1 .

Gene expression data for LUAD were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) using the R software package “TCGAbiolinks”.

Clinical data were directly acquired from TCGA, comprising 535 LUAD patients and 59 healthy controls. After rigorous quality control and preprocessing, which included removal of duplicate and incomplete records, we established our training set with 55 LUAD patients and 38 healthy controls. Concurrently, we constructed a TCGA validation set by randomly selecting 30 samples with complete grouping information and unbiased gene expression profiles, excluding those used in the training set.

From January 2020 to April 2022, we prospectively collected medical records and tissue specimens from 10 patients undergoing surgical resection for primary LUAD at the First Affiliated Hospital of Wenzhou Medical University. These 10 LUAD specimens and their matched adjacent normal tissues (n = 10) were subsequently processed for gene sequencing analysis, serving as an independent external validation cohort. This study protocol was approved by the Institutional Review Board of the First Affiliated Hospital of Wenzhou Medical University (Approval No. YS2023-706).

Paired-end sequencing (150 bp) was conducted on the Illumina NovaSeq platform (Illumina, California, USA) following standard protocols. Raw sequencing data were processed on the BMK Cloud bioinformatics platform (https://www.biocloud.net) for downstream analysis.

In the R software environment, the “limma” package was used to normalize and annotate probe data from a training set of LUAD patients (55 cases) and control population (38 cases) in public databases. Differential gene expression analysis was performed to identify differentially expressed genes (DEGs) between the lung cancer group and the control group. To visualize the expression patterns of these DEGs, hierarchical clustering was conducted using the heatmap package, while the ggplot2 package was employed to generate volcano plots, providing an intuitive representation of gene expression differences.

To elucidate the interactions among genes in LUAD, we employed WGCNA on expression profile data from the TCGA training set. Data integrity was assessed using the “Good Samples Genes” function to ensure robust downstream analysis.

Following the removal of duplicate entries, IRGs were retrieved from the Immunology Database and Analysis Portal (ImmPort, https://www.immport.org), yielding a total of 1,793 unique genes (10). These genes were subsequently cross-referenced with those identified through gene clustering analysis using Venn diagrams, enabling the precise selection of immune-related genes.

To investigate the association between DEGs and patient overall survival, univariable Cox regression analysis was performed. Hazard ratios (HRs) were calculated to classify genes as either protective (HR < 1) or risk factors (HR > 1) for prognosis. For visual representation of the results, a volcano plot was generated to illustrate the statistical significance and risk ratios of DEGs. Additionally, a forest plot was employed to display detailed univariable Cox regression results, including HR values and their 95% confidence intervals (CIs).

Gene Ontology (GO) term (https://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (https://www.genome.jp/kegg/) enrichment analyses were conducted using the R package “cluster Profiler” (11).

In this study, we adopted three distinct machine learning algorithms—LASSO, SVM-RFE, and RF—to screen out HUB genes from specific modules. By independently applying these three algorithms and intersecting their selected genes, we identified robust diagnostic biomarkers.

To evaluate the differential expression of HUB genes between tumor and normal tissues, we generated box plots using the “ggplot2” package. Statistical significance was assessed using Student’s t-tests. To ensure the robustness of our results, we further validated these findings in an independent external dataset.

For the key genes identified during screening and validation, we utilized RNA-Seq data from LUAD, and adjacent non-tumor tissues obtained from the TCGA database. To assess the diagnostic potential of individual genes, we performed logistic regression analysis, treating gene expression levels as the independent variable and LUAD status (tumor vs. non-tumor) as the dependent variable. The diagnostic performance of each gene was quantified by calculating the area under the receiver operating characteristic (ROC) curve (AUC).

We employed Pearson correlation analysis to examine pairwise relationships between key genes, followed by visualization using heatmaps. This approach enables the identification of co-expression patterns and potential functional interactions among genes, providing insights into their regulatory networks.

To investigate the potential associations between gene expression levels and clinicopathological characteristics, we performed chi-square tests in this study.

We constructed a prognosis model for lung adenocarcinoma patients using univariate and multivariate Cox regression analysis methods. First, univariate Cox regression analysis was used to identify individual variables significantly associated with patient survival time. Subsequently, multivariate Cox regression analysis further evaluated the independent effects of these variables on patient prognosis while considering other potential confounding factors. The risk score was calculated as a linear combination of gene expression levels and corresponding Cox regression coefficients. To facilitate clinical use, we created a nomogram using the “rms” package in R.

Model performance was assessed through: 1. Calibration curves to compare predicted vs. actual survival, validated on an independent dataset; 2. ROC curves via the “pROC” package, with AUC values measuring predictive accuracy across thresholds.

We used the “neuralnet” and “neuralnettools” packages in R to construct an ANN diagnostic model based on hub gene expression. The model architecture included an input layer (hub genes), two hidden layers (with 8 and 3 neurons, respectively), and an output layer representing “normal” and “tumor” classes, with softmax activation. To avoid overfitting, we applied 10-fold cross-validation and used gene weight information during training.

Model performance was assessed using accuracy, precision, recall, F1-score, and AUC in both the training and validation sets. ROC curves were also used to evaluate the discriminative ability of both individual hub genes and the ANN model across the TCGA training set, internal validation set, and external cohort. ROC analysis assessed diagnostic accuracy, with AUC values reflecting its classification ability.

To identify potential molecular subgroups in LUAD samples, we applied the ConsensusClusterPlus algorithm (12), which determines the optimal number of clusters (k) by resampling the dataset and evaluating clustering stability.

To validate the robustness of the classification, we used two high-dimensional visualization techniques: UMAP and t-SNE. Additionally, Kaplan-Meier survival analysis was performed to assess the prognostic relevance of the identified subgroups by comparing survival curves across groups.

To assess immune infiltration patterns in LUAD, we leveraged the TIMER database (https://cistrome.shinyapps.io/timer) (13). Here, we utilized TIMER to examine correlations between the expression of four candidate genes (CBLC, GDF10, LTBP4, and FABP4) and Tumor immune infiltrating cells (TII Cs) populations, providing insights into their potential immunomodulatory roles in LUAD.

For a more comprehensive evaluation of the tumor immune microenvironment, we applied single-sample gene set enrichment analysis (ssGSEA) (14). Using the GSVA and GSEA Base R packages, we quantified immune infiltration levels in each LUAD sample based on 28 immune-related gene sets curated from the TISIDB (http://cis.hku.hk/TISIDB/) database.

To investigate the transcriptomic characteristics of LUAD tumor samples and identify key biological pathways potentially involved in disease pathogenesis, we performed Gene Set Enrichment Analysis (GSEA). Pathways with a nominal P-value (NOM P) < 0.05 and a false discovery rate (FDR Q-value) < 0.05 were considered statistically significant.

The Human Protein Atlas (HPA; https://www.proteinatlas.org/) is a comprehensive online repository that maps the expression and localization of human proteins across various cells, tissues, and organs using diverse biological techniques. We analyzed the immunohistochemical (IHC) staining patterns of selected hub genes in normal lung tissue and lung cancer specimens using data and analytical tools available on the Human Protein Atlas platform.

Significant differences between groups were assessed using the log-rank test and univariate COX regression analysis, with P-values and HR with 95% CI calculated. Multivariate COX regression analysis and stratified analysis were performed to evaluate the independence of the risk score model. The performance of gene prognostic signatures was assessed using ROC curves, with the area under the curve (AUC) as the metric. Statistical significance was defined as P < 0.05. All statistical analyses were conducted using R language version 4.3.0.

The original mRNA expression data for LUAD were obtained from TCGA database, comprising 535 LUAD tumor samples and 59 matched adjacent normal tissue samples. Following quality control and data filtering, our final study cohort consisted of 55 LUAD samples as the experimental group and 38 adjacent normal tissue samples as controls. The complete clinical characteristics of the 55 LUAD cases are summarized in Table 1 .

There is a significant difference between the tumor group and the control group. Differential analysis identified a total of 1822 DEGs (log2FC > 2, P < 0.001), including 719 downregulated genes and 1103 upregulated genes. These differences are intuitively displayed in the volcano plot ( Figure 2A ) and heat map ( Figure 2B ).

We adopted the WGCNA method to identify gene modules associated with LUAD. First, by evaluating the scale-free fit index and average connectivity under different soft thresholds, the optimal soft threshold was determined to be β=17 ( Figure 3A ). Subsequently, sample clustering analysis was performed ( Figure 3B ), a topological overlap matrix was constructed, and hierarchical clustering was applied to identify modules. By using the dynamic tree cut algorithm and merging modules with similarity greater than 0.75, a final module clustering map was obtained ( Figure 3C ). Then, correlation analysis between modules and clinical features was conducted, and a heatmap was drawn to show the correlation coefficients and P-values ( Figure 3D ). A total of 4 modules were identified, among which the blue module with the highest correlation with the tumor group (r=0.96, P=4e-50) was selected, containing 604 genes.

A set of data consisting of 1793 IRGs was obtained from IMMPORT database (15). To further explore the immune genes associated with LUAD, we compared these immune genes with 604 genes screened by gene clustering methods. A total of 102 immune genes closely associated with LUAD were identified, and these findings were visualized in a Venn diagram ( Figure 4A ).

After a one-way regression analysis of 102 differential immune genes related to LUAD, the results showed that 68 differential immune genes were significantly correlated with the overall survival rate of LUAD patients (P < 0.05) ( Figures 4B,C ). These significantly related differential immune genes may affect tumor progression and patient survival prognosis by affecting the immune microenvironment of LUAD.

GO and KEGG analyses revealed the functional localization of 68 differentially expressed immune genes significantly associated with overall survival in lung adenocarcinoma. GO analysis indicated that these genes are mainly involved in signal pathway regulation (such as cell surface receptor signaling, G protein-coupled receptor pathway), cell proliferation, substance metabolism (such as protein phosphorylation), and immune inflammatory response (such as chemotaxis, MAPK cascade regulation) ( Figure 5A ). KEGG analysis further found that these genes are enriched in immune-related pathways (such as PI3K-Akt and Jak-STAT pathways) and tumor proliferation-related pathways (such as Ras, Rap1, and MAPK signaling pathways) ( Figure 5B ). These results help to better understand the role of differentially expressed immune genes in the occurrence and development of LUAD.

Machine learning identified 68 candidate genes. LASSO, SVM-RFE, and RF algorithms further narrowed them to 16, 18, and 16 genes, respectively. By overlapping the three groups, four key diagnostic biomarkers were finally identified ( Figure 6 ). The information related to 4 genes is shown in Table 2 .

In lung cancer research, by comparing the lung cancer samples in the TCGA cohort with the normal control group, it was found that the expression level of CBLC genes was significantly increased, while the expression levels of the three genes FABP4, GDF10, and LTBP4 were significantly decreased ( Figure 7A ). This result was validated in an independent externally sequencing sample cohort ( Figure 7B ), showing gene expression patterns consistent with the TCGA cohort, enhancing the reliability and universality of this finding.

To evaluate the predictive utility of these genes as biomarkers of lung cancer, the study performed a ROC curve analysis. In the TCGA cohort, the AUC values ​​of CBLC, FABP4, GDF10, and LTBP4 genes were 0.982, 0.997, 0.980, and 0.983, respectively ( Figure 7C ), indicating that these genes have a high diagnostic efficiency in distinguishing lung cancer samples from normal control groups. In an independent external sequencing sample cohort, these genes also performed impressively, with the AUC value for CBLC at 0.800, while the AUC values for FABP4, GDF10, and LTBP4 all reached 1.000 ( Figure 7D ), further validating their potential as diagnostic markers for lung cancer.

To gain a deeper understanding the role of the four genes CBLC, FABP4, GDF10, and LTBP4, and their genetic interactions in LUAD, Pearson correlation coefficients were calculated. The calculation results ( Figure 8A ) show that, except for CBLC, there is a good correlation between the three genes FABP4, GDF10, and LTBP4, indicating that they may play similar or synergistic roles in the development of LUAD.

In addition, through statistical analysis of the TCGA dataset, the study further explored the relationship between these gene expression levels and the clinicopathological status of patients with LUAD ( Figures 8B–E ). The analysis results revealed, the high expression of CBLC in LUAD is positively correlated with poor clinical and pathological stages in patients, which means that the high expression of CBLC may be an indicator of poor progression and prognosis of lung adenocarcinoma. In contrast, high expression of FABP4, GDF10, and LTBP4 is associated with better clinical and pathological stages, indicating that the increased expression levels of these genes may be related to a better prognosis of LUAD.

A multi-factor COX regression analysis was used to construct an immune prognostic model for predicting the 1-year, 3-year, and 5-year survival rates of LUAD patients, which was visualized as a nomogram ( Figure 9 ). The analysis results of the calibration curve ( Figures 10A, C ) showed that there is a good consistency between the model predictions and the ideal model. In addition, the AUCvalues of the model training set were 0.80 (95% CI: 0.91–0.69) for 1-year survival rate, 0.85 (95% CI: 0.95–0.76) for 3-year survival rate, and 0.83 (95% CI: 0.95–0.70)for 5-year survival rate, respectively ( Figure 10B ), while in the TCGA validation set, the AUC values were 0.76 (95% CI: 0.97–0.55), 0.76 (95% CI: 0.88–0.65), and 0.77 (95% CI: 0.93–0.61), respectively ( Figure 10D ).

Furthermore, based on four HUB genes, an ANN model was developed ( Figure 11A ) to enhance prediction accuracy. The TCGA dataset was used as the training set, while independent sequencing results served as the test set. The diagnostic performance of the ANN model was evaluated using ROC curves. In the training set, the model achieved an AUC of 0.929 (95% CI: 0.878–0.980) ( Figure 11B ), while in the TCGA validation set and independent external validation set, the AUC values were 0.876 (95% CI: 0.812–0.940) and 0.906 (95% CI: 0.849–0.963), respectively ( Figures 11C, D ).

When analyzing the gene expression profiles of 55 tumor samples in the TCGA training set cohort, we first performed cluster analysis using the consensus clustering method. The analysis determined that when k=2, the classification of the samples was both highly reliable and stable ( Figures 12A–C ). To further validate the differences between these two subgroups, two dimensionality reduction techniques, Uniform Manifold Approximation and Projection (UMAP) and t-Distributed Stochastic Neighbor Embedding (t-SNE), were employed. The results confirmed the significant differences between the two subgroups ( Figure 12D ). Based on this classification, the samples were divided into two clusters: Cluster 1 (C1, N=29) and Cluster 2 (C2, N=26). Next, survival analysis was performed to explore the prognostic differences between these two subgroups. The analysis results showed ( Figure 12E ) that there were significant survival differences between the two subgroups (P < 0.001).

This study utilized the ssGSEA algorithm to analyze the immune cell infiltration and immune gene correlation in LUAD samples and the control group. The results show ( Figure 13B ): CBLC is negatively correlated with various immune cells (such as activated CD8 T cells, dendritic cells, regulatory T cells, etc.); FABP4 is positively correlated with various immune cells (such as memory T cells, macrophages, neutrophils, etc.); GDF10 is positively correlated with memory CD4 T cells, eosinophils, plasma cell-like dendritic cells, etc.; LTBP4 is positively correlated with memory T cells, eosinophils, and dendritic cells. Further analysis ( Figure 13C ) found that CBLC is mainly negatively correlated with dendritic cells, macrophages, and CD8 T cells; FABP4 is positively correlated with dendritic cells, macrophages, and neutrophils; GDF10 and LTBP4 are mainly positively correlated with macrophages and CD4 T cells, and LTBP4 is also positively correlated with dendritic cells. These results indicate that different genes play different roles in immune cell infiltration, reflecting the complexity and diversity of the LUAD tumor microenvironment.

Further research compared the differences in 28 immune cell types between two immune subgroups and found significant differences and heterogeneity in immune cell infiltration between the high-risk and low-risk groups ( Figure 13A ). Compared to the C2 group, the C1 group had significantly elevated levels of central memory CD4 T cells, follicular helper T cells, regulatory T cells, natural killer cells, myeloid-derived suppressor cells, plasmacytoid dendritic cells, and mast cells ( Figure 13D ). In further studies, we found significant increases in CD4 T cells, CD8 T cells, macrophages, and dendritic cells in the C1 group ( Figure 13E ).

According to the immune staining intensity comparisons provided by the HPA database ( Figure 14 ), the expression levels of FABP4, GDF10, and LTBP4 were higher in normal samples, while the expression level of CBLC was higher in lung adenocarcinoma samples. This result is consistent with previous studies on gene expression differences.

The GSEA analysis of four HUB genes in the progression of lung adenocarcinoma (LUAD) indicates that these genes are mainly enriched in the following biological processes and pathways ( Figure 15 ): 1. Allogeneic rejection response 2. Asthma 3. Amino acid biosynthesis 4. DNA replication.