Microarray expression data and clinical information for LUAD were obtained from the GEO and TCGA databases. There were two cohorts in the treatment group, GSE63459 and GSE176348, with 89 specimens (including 45 tumor samples and 44 normal samples). In addition, external validation using the TCGA-LUAD cohort with 598 samples (including 539 tumor samples and 59 normal samples) was performed. All sequencing information for normal samples comes from adjacent tissues.

We used the data normalization and probe annotation from the R software (version 4.2.1) “limma” and “GEOquery” packages for the data of GSE63459 and GSE176348, with P-value < 0.05 and |log fold change (FC) | > 1 for the DEGS screening criteria (13, 14).

We used the WGCNA to process expression profile data from GSE63459 and GSE176348 datasets to establish a weighted co-expression network. Then we investigated the genes that deviate from the top 25% of the median (10). The data integrity is checked by the ‘Good SampleGenes’ function. We chose a suitable soft threshold value (β) and validated the ability of the soft threshold value. The matrix data was transformed into an adjacency matrix by us, followed by clustering to identify modules based on the topological overlap. Then, the module feature element (ME) is calculated, the similarly expressed modules are combined into the cluster tree according to the ME, and we draw the hierarchical cluster tree graph. Then, the module and phenotype data are combined, and then the gene significant (GS) and module significant (MS) are calculated; the calculation results are used to evaluate whether the gene and clinical information are essential and to analyze the correlation between the module and the model. In addition, We calculated the module membership (MM) for each gene to analyze the GS values of each module.

The gene with the highest inter-module connectivity was selected as the candidate HUB gene. The GS values for biologically significant genes are also generally higher. Therefore, we chose candidate genes with an absolute GS value> 0.20 and an MM value > 0.80. We then intersected these genes with DEGS using the “glmnet” package in the R software package and used the LASSO analysis to determine the final HUB genes (11). We used box plots to probe the HUB gene expression levels in LUAD samples and healthy samples. With the help of ‘MedCalc’ software (version 2.0.1), we draw the receiver operating characteristic curves (ROC) to determine these HUB genes’ diagnostic specificity and accuracy. A dataset (TCGA-LUAD) is also available for external verification of the HUB gene’s expression level and diagnostic value.

With the help of the “Survival” and “SurvMiner” packages in the “R” software, we divided the samples in TCGA-LUAD into two groups (high and low expression groups) using the median expression of the HUB gene. Lastly, survival curves for HUB and LUAD genes were plotted using the Kaplan-Meier method with the aid of the software package “ggplot2”.

Results of immunohistochemical staining of the HUB gene in normal lung tissue and lung cancer specimens from The Human Protein Atlas(www.proteinatlas.org).

We quantified the infiltration of 28 immunocytes in the GSE63459 and GSE176348 datasets using the ssGSEA algorithm (15). The box plots we established indicate the differences in the expression levels of these immune cells. We also calculated the Spearman correlation of these immune-infiltrating cells with the candidate HUB genes and visualized the calculated results with the ‘ggplot2’ program package.

We performed GO analysis of DEGs, KEGG analysis, and GSEA analysis through the ‘clusterProfiler’ and ‘enrichplot’ package of the R software package (16). We used the immunological signature genomes from the Molecular Signature Database (MsigDB) as the reference, and the significantly enriched genomes had to meet the P <0.05 and the false discovery rate (FDR) q-value <0.05.

Several HUB genes (DPYSL2, OCIAD2, and FABP4) were selected for study to verify the HUB genes’ role further. Normal human lung epithelial cells (BEAS-2B) and lung adenocarcinoma cells (A549 and H1299) were collected for culture, moreover, extracted RNA from the three cells using the Trizol reagent. For cDNA synthesis, the synthesis was performed using the reverse transcription reagent VAZYME R232. The final PCR reaction was performed on a quantitative real-time PCR instrument. The reaction parameters included the denaturation process (30s at 95°C), followed by 40 PCR cycles (5s at 95°C and 34s at 60°C). The melting curve of the PCR product was established, and after the amplification reaction, the temperature was slowly heated from 60°C at (95°C, 15s, 60°C, 60s, 95°C, 15 s) to 99°C (instrument automatic ramp rate 0.05°C/s). Quantitative real-time PCR reactions were performed for target and housekeeping genes for each sample. We calculated the relative expression levels of the three genes using the 2 ^ (-δδ ct) method. Since the experimental results of FABP 4 and DPYSL2 were fit to a normal distribution, the analysis was performed using the one-way ANOVA test. We used the Kruskal-Wallis test for statistical analysis for OCIAD2 genes whose results do not conform to the normal distribution ( Supplementary file 1 ). The sequence of the primers is as follows:

DPYSL2, 5’-CCCTGCAGAACATCAACGAC-3’ (forward) and 5’-GGCATCTGGAAACGAGTGTG-3’ (reverse); OCIAD2, 5’-GTCTGCTCGTGGAAACCAAG-3’ (forward) and 5’-CAAGAGACCAGCAAGTGCAAC-3’ (reverse); FABP4, 5-GATGACAGGAAAGTCAAGAGCAC-3’ (forward) and 5’-GACGCCTTTCATGACGCATTC-3’ (reverse); and GADPH,5’-TCTGACTTCAACAGCGACACC-3’ (forward) and 5’-CTGTTGCTGTAGCCAAATTCGT-3’ (reverse).

The absolute deviations in the median top 25% of genes were selected for constructing WGCNA, and missing values and outliers in the samples were subsequently removed by cluster analysis. To maintain consistency with the scale-free network, we set the soft threshold β to 5 (scale-free R2 = 0.91; slope =-1.67) ( Figures 1A, B ). We also analyzed the gene expression in the normal and LUAD samples and plotted the results as a heatmap( Figure 1C ). We built a co-expression matrix and obtained seven modules with the help of dynamic hybridization shear ( Figure 2A ). The relationship between these seven modules and the LUAD and healthy samples is shown in the heatmap. The turquoise module has the highest correlation (cor) (cor = 0.89; P=1e-31) ( Figures 2B, C ). Moreover, after our correlation analysis, we found that in the turquoise module, GS and MM are also well correlated(cor=1.51; P=1.3e-08) ( Figure 2D ).

For the DEGs, our filtering criteria were P <0.05 and | logFC |> 1, including 354 differential expressed genes and displaying the results on the volcano plot ( Figure 3A ). We further selected 87 genes with higher connectivity in the turquoise module using |GS|> 0.20 and |MM|> 0.80 as screening criteria. The results of these two screens were compared, and their intersection was selected as candidate HUB genes, and 85 genes were combined ( Figure 3B ). Ultimately, after a further screening of these genes using the LASSO analysis, we were able to obtain 12 genes (ADAMTS8, CD36, DPYSL2, FABP4, FGFR4, HBA2, OCIAD2, PARP1, PLEKHH2, STX11, TCF21, TNNC1) ( Figures 3C, D ).

Next, we investigated the function of the DEGs screened above; We performed GO analysis on 354 genes. According to the results, we know that DEGs mainly focus on the regulation of genes or pathways (e.g., transforming growth factor receptor signaling pathway, positive regulation of gene expression, negative regulation of transcription by RNA polymerase II promoter), vascular development (e. g., angiogenesis, angiogenesis, vascular development), immune response and inflammatory response (e. g. inflammation, cellular response of interleukin-1, positive regulation of interleukin-6 production), and even play an essential role in alveolar development ( Figure 4 ). According to the KEGG analysis, We can also learn that DEGs are mainly enriched in the following pathways, pathways of immunological and inflammation-related diseases (AGE-RAGE signaling pathway in diabetic complications, fluid shear stress, and atherosclerosis, rheumatoid arthritis), there are also some immune-related pathways (IL-17 signaling pathway, TNF signal channel, the TGF signaling pathway, PPAR signaling pathway, and other pathways)( Figure 5 ). The GO analysis and KEGG analysis showed that there are many biological processes and signaling pathways related to the immune-inflammatory response in the development and development of LUAD.

IHC staining results were paired as shown in Figure 6 . On the left of each pair of images is the gene staining in normal lung tissue, and on the right is the staining in lung cancer samples. We can estimate the expression level of HUB genes, and it is clear that two HUB genes, OCIAD2 and PARP 1, have higher expression in lung cancer samples ( Figures 6A, B ), while the remaining HUB genes have more expression in normal samples ( Figures 6C–K ). Unfortunately, we were unable to find the PLEKHH2 staining results, and we will continue to follow up on the results in follow-up studies.

We assessed the expression levels of the 12 HUB genes by box plots. The results indicated that only OCIAD2, PARP 1 were significantly increased in the control group ( Figures 7G, H ) (P <0.001), while the other ten genes, ADAMTS8, CD36, DPYSL2, FABP4, FGFR4, HBA2, PLEKHH2, STX11, TCF21, TNNC1 were higher in the control group ( Figures 7A–F, I–L ). Next, we also externally verified the expression levels of these 12 HUB genes using the TCGA-LUAD dataset, and validation results were in agreement with our experimental group ( Figures 7M–X ). In the ROC curve analysis of the 12 HUB genes, the area under the curve (AUC) of the HUB gene represents the sensitivity and specificity of the HUB gene for the diagnosis of LUAD. From the ROC curve, we can know that the AUC values of all 12 HUB genes were> 0.93, indicating the high value of HUB genes for the diagnosis of LUAD ( Figure 8A ). While in the TCGA-LUAD cohort, the AUC values, except for PAPR 1 and PLEKHH2, were 0.884 and 0.839, respectively. The AUC values for the remaining HUB genes were all> 0.95, which coincident with the findings from our cohort study above ( Figure 8B ).

We partitioned LUAD samples into two groups (high and low expression groups) based on the TCGA-LUAD database. Kaplan-Meier curves were performed for the HUB gene in order to analyze its prognostic relationship to LUAD patients. According to the results, the high expression of OCIAD2 and PARP1 is linked to poor prognosis in LUAD patients ( Supplementary Figures 1G, H ). However, high expression of ADAMTS8, CD36, DPYSL2, PLEKHH2, STX11, and TCF21 tends to lead to better prognosis ( Supplementary Figures 1A–F, I–L ), which coincides with the difference in expression of these HUB genes in normal samples and lung cancer samples.

We compared and analyzed the immune cell infiltration of LUAD samples and the control group with the ssGSEA algorithm. The graph shows the distribution of 28 immunocytes in two datasets, GSE63459 and GSE17634 ( Figure 9A ). Shown according to its results, the CD4+T cells, CD8+T cells, natural killer (NK) cells, and Bcell in LUAD samples were higher than that in normal samples, and this result indicates that these cells play a significant role in the progression of LUAD ( Figure 9B ). According to the correlation analysis of HUB genes and immune cell infiltration, most of these HUB genes showed a positive correlation with immune infiltrating cells, such as macrophage, CD4 + T cell, CD8 + T cell, and NK cell. CD8 + T cell exerts antitumor effects in a wide range of cancers. It should be noted that OCIAD2 and PARP1 genes are negatively associated with numerous immune cells. This fits with their results leading to the poor prognosis associated with patients with LUAD ( Figure 9C ).

To make out that immunogenetics may exist during the progression of LUAD, we have used the immunologically signature gene set in the MsigDB database as a reference standard for DEGS GSEA. A total of 819 gene sets were significantly enriched (|normalized enrichment score (NES)|> 1; FDR Q value <0.05). These genes were mainly concentrated in CD8 T cells, NK Cells, CD4+T cells, monocytes, and regulatory T cells ( Supplementary Table S1 ). Based on the above findings, it appears that many immune genes play an essential role in LUAD progression ( Figure 10 ).

After statistical analysis of the PCR results, we found that DPYSL2 (p <0.001), OCIAD2 (p <0.05), and FABP4 (p <0.001) had higher expression in A549 compared to BEAS-2B, showing a statistically significant difference. The expression content of DPYSL2 was lower in H1299 cells compared to t BEAS-2B (p <0.01). Although the expression difference of FABP 4 and OCIAD2 genes in H1299 was not statistically significant, they both showed a trend to increase ( Figure 11 ). These experimental results can better support our study. Nevertheless, the results may require further study with a larger sample size.