Identification of Co-Expression Gene Modules
We used the WGCNA package on R to construct gene co-expression networks from the TCGA-LUAD and GSE32863 datasets. Totally, 15 modules in the TCGA-LUAD and 7 modules in the GSE32863 (Figure 2A and Figure 3A) were identified, and a gray module was excluded, because it was assigned into no cluster. Then, the heatmap of module-trait relationships was plotted to evaluate the relevance between the two clinical traits and each module. The results in Figure 2B and 3B, revealed that the turquoise module in the TCGA-LUAD and brown module in the GSE32863 were proved to be highest associated with normal lung tissue (turquoise module: r = 0.79, p = 5e−127, and brown module: r = 0.97, p = 9e−71). Moreover, compared with normal tissues, the expression of the genes in the turquoise module and brown module was all downregulated in LUAD tissues.
Extraction of Overlapping Genes between the DEGs and Co-expression Networks
We discovered that 3,583 DEGs in the TCGA-LUAD dataset (Figure 4A) and 957 DEGs in the GSE32863 dataset (Figure 4B) were dysregulated in tumor tissues using the limma package (cut-off criteria of |logFC| ≥ 1.0, adj. p < 0.05). A total of 5313 and 1699 co-expression genes were enrolled in the turquoise module of TCGA-LUAD dataset and the brown module in GSE32863, respectively. Totally, as shown in Figure 4C, 358 overlapping genes were extracted to verify the genes among co-expression networks.
Functional Enrichment Analysis
In order to further explore the potential functions of the 358 overlapping genes, GO and KEGG enrichment analyses were performed by the cluster profiler package. Several GO-enriched gene sets were displayed in Figure 5. The results of the biological process (BP) revealed that these genes are mainly enriched in extracellular structure organization, cell−substrate adhesion and response to acid chemical. The cellular component (CC) of these genes were mainly involved in cell−cell junction and collagen−containing extracellular matrix. Most importantly, molecular function (MF) analysis proved that DNA−binding transcription activator activity, glycosaminoglycan binding and growth factor binding were suggested to be associated with these genes. In the KEGG analysis, the main pathways were enriched by these genes such as Drug metabolism − cytochrome P450, Ether lipid metabolism and Arachidonic acid metabolism (Figure 6).
Construction of PPI Network and Identification of Hub Genes
We applied the STRING database to establish the PPI network among the overlapped genes (Figure 7A). Figure 7B gives information that hub genes were extracted from the PPI network with a connectivity degree ranking. Finally, the top 12 highest-scored genes, including GNG11, ADRB2, ADCY4, TGFBR2, IL6, GPC3, VIPR1, GRK5, CAV1, RAMP3, RAMP2, and CALCRL were selected as hub genes.
Clinical relevance of hub gene expression in LUAD patients
After the 12 hub genes were screened out, based on the dataset in TCGA, we validated that ADCY4 expression had a strong positive correlation with tumor individual cancer stages by applying package cliCor in R (p<0.05; Figure 8A). Then with VIPR1, stage I had a significant expression difference with stage IV (p=0.033, Figure 8E). Finally, TGFBR2 expression of stage I and stage III had significant differences (p=0.024, Figure 8I).
Prognostic Values and Protein Expression of 12 Hub Genes
We applied the R survival package and the GEPIA2 database to perform OS (Figure 9) and DFS (Figure 10) analyses of the twelve hub genes by Kaplan–Meier plotter, respectively. The Kaplan–Meier analysis revealed that the lower expression level of ADCY4, VIPR1, and TGFBR2 was obviously related to the worse OS of the LUAD patients (ADCY4: p=0.003, VIPR1: p=0.019, and TGFBR2: p=0.002) (Figure 9A, 9E, and 9I). No significant difference was found in LUAD patients with the expression level of ADCY4, VIPR1 and TGFBR2 in DFS (p>0.05) (Figure 10A,10E, and 10I). Furthermore, based on the HPA database, the protein level of the ADCY4 and VIPR1 gene was extremely lower in LUAD tissues, compared with normal tissues (Figure 11). The expression level of TGFBR2 is not discovered in the HPA database.