The microarray dataset GSE32863 of human LUAD with patient clinical information was extracted from Gene Expression Omnibus (GEO) database (Selamat et al., 2012). GSE32863 dataset was obtained using Illumina HumanWG-6 v3.0 expression beadchip. It includes 58 LUAD samples, which were used to construct a co-expression network followed by extraction of hub genes. R package was used to annotate raw data, generate an expression matrix, and match probes with official gene symbols. The median absolute deviations (MADs) were ranked from the largest to smallest, and the expressions of the top 25% genes with the largest differences in the samples were extracted for in-depth analysis.

We used the R package “WGCNA” to perform WGCNA on selected genes to find out the expression patterns between different genes (Langfelder and Horvath, 2008). Genes with similar expression patterns were regarded as a specific module and marked with a unique color. Subsequently, the correlation between these modules was calculated, and a heatmap was created to show the independence between each module. Then, a correlation analysis was performed to find modules related to clinical information. Modules closely related to tumor recurrence were chosen for further analysis.

In order to understand the biological functions and signal pathways involved in the genes from the selected modules, differentially expressed genes (DEGs) in these modules were analyzed through the R package “clusterProfiler” (Yu et al., 2012). Then, another R package “ggplot2” was used to perform the top 10 enrichment terms of Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Enriched GO terms and KEGG pathways were identified base on the cutoff criterion of p < 0.05.

For hub gene screening, genes from the selected modules were uploaded to the STRING database to build a protein–protein interaction (PPI) network (Szklarczyk et al., 2019). The interaction score >0.4 was defined as the threshold of the key genes in the PPI network. Then, a Cytoscape plug-in cytoHubba was used to extract the top 1% targets in these modules based on the degree method (Shannon et al., 2003).

Based on the cytoHubba analysis, the top 15 genes in the selected modules and the genes enriched in KEGG pathway were selected as candidate hub targets for further analysis. In this study, the Kaplan–Meier plotter (http://kmplot.com/analysis/) was used to draw the survival analyses of the candidate hub targets, while the Gene Expression Profiling Interactive Analysis (GEPIA) webserver (Tang et al., 2017) (http://gepia.cancer-pku.cn/) was used to verify the outcomes of survival analysis so as to screen out real hub genes. With the use of the Human Protein Atlas (http://www.proteinatlas.org) database, the real hub genes were further validated by immunohistochemistry (IHC). Moreover, Single Cell Portal (https://singlecell.broadinstitute.org/single_cell) was used to detect the expression of key targets in various cell types of lung cancer. In addition, the cBioPortal tool (http://www.cbioportal.org/; version: 2.2.0) was used to compare the genetic variations of the real hub genes in LUAD. Finally, the cBioPortal tool was also used to present the co-expression analysis of real hub genes.

In the validation set GSE116959, samples of LUAD were divided into low and high groups according to the expression level of the genes. To explore the potential function of real hub genes in LUAD, gene set enrichment analysis (GSEA) was conducted and mapped into KEGG pathway enrichment database. The terms with p-value <0.05 were chosen as the cutoff criteria.

A regulating network on hub genes and LUAD was constructed by Cytoscape software. Then, the plugin iRegulon of Cytoscape was applied to forecast transcription factor (TF) regulation networks. In addition, Comparative Toxicogenomics Database (http://ctdbase.org/), a robust and publicly available database that aims to reveal how environmental exposures affect human health, was used to search for drugs for the hub genes. Chemicals supported by at least one database were selected as the potential drugs. Potential drug–key target interaction was constructed by Cytoscape. The final list only involves drugs that may interact with at least two key targets and have been reported to have an anticancer pharmacological activity.

The flowchart of strategy in this study is summarized in Figure 1. GSE32863, a GEO dataset containing gene expression profiling data of 58 matched LUAD and non-tumor lung samples, was analyzed using the R package WGCNA. Clinical characteristics of tumor pathological stage, gender, recurrence, egfr, kras, lkb1, and smoking status of LUAD patients were denoted. After screening by MADs arranged from large to small, the expression of the top 25% genes (6,360 genes) with the greatest differences in samples were selected for further analyses.

After classification of DEGs with similar expression patterns into modules by average linkage clustering, 14 modules were finally identified by merging similar modules when the MedissThres was set at 0.25 (Figure 2A). Among them, module eigengenes of three modules (green, purple, and brown) were found to be more related to the tumor recurrence as others. Scatter plots further indicated positive correlations between members of these three modules and gene significance for LUAD recurrence (Figures 2B–D). Herein, the green, purple, and brown modules were chosen as key modules associated with LUAD recurrence for further investigation.

To understand the biological functions of genes in the selected modules, GO and KEGG pathway analyses were firstly applied. As depicted in Figure 3, genes in three key modules were expected to exert their functions in terms of cell adhesion molecule binding, endosome membrane, ribosomal subunit and autophagy, etc. Meanwhile, these genes could only be enriched in one single cellular pathway called Tight Junction signaling pathway (Figure 4 and Table 1).

Next, genes in three key modules were extracted to establish a PPI network through STRING database and Cytoscape software. As a result, a total of 926 DEGs of the three key modules were mapped into the PPI network, including 1,426 nodes and 9,522 edges. In addition, the top 15 key targets within PPI network were selected with the help of cytoHubba plug-in in Cytoscape software based on rank of degree (Figure 5).

According to the results above, 40 genes (25 genes could be mapped into the KEGG pathway plus 15 key targets in PPI network) remained the potentials to be the hub genes for LUAD recurrence. For the purpose of scope reduction and effect evaluation, overall survival analysis was firstly performed on these candidate genes using the survival information of Kaplan–Meier plotter. After the first round of screening, only genes with significant impact on LUAD patient survival were further validated by survival analysis with GEPIA. As a result, a total of eight candidate hub genes (ACTR3, ARPC5, RAB13, HNRNPK, PA2G4, WDR12, SRSF1, and NOP58) were identified to be significantly correlated with LUAD patient survival with both survival analysis methods (Figure 6). Remarkably, all eight candidate genes were found to be correlated with poor LUAD patient survival through GEPIA, five of which (ACTR3, ARPC5, RAB13, HNRNPK, and SRSF1) were inversely revealed to be associated with favorable outcome using Kaplan–Meier plotter.

We next sought to investigate the expression of these candidate genes in LUAD. Firstly, comprehensive single-cell RNA sequencing data mapping myeloid populations in non-small cell lung tumor and peripheral blood was utilized to detect the expression of our candidate genes in different types of immune cells (Zilionis et al., 2019). As presented in Figure 7, candidate genes ACTR3, ARPC5, and HNRNPK were highly expressed in almost all types of lung cancer immune cells, while RAB13, PA2G4, and WDR12 were found to be less abundant in the majority of myeloid populations in lung cancer. As compared with the rest of the hub genes, ACTR3, ARPC5, and HNRNPK are all characterized in regulation of actin polymerization (Yoo et al., 2006; Kopitar et al., 2019). Thus, the hyperactivation of these three genes in lung cancer cells may reflect the significance of actin polymerization during the development of cancer. Next, the distribution and expression of these eight hub gene candidates were checked using IHC staining data of both LUAD and adjacent non-tumor samples from The Human Protein Atlas database. Representative IHC images revealed elevated protein levels of HNRNPK, PA2G4, WDR12, SRSF1, and NOP58 in LUAD samples (Figure 8), suggesting the dysregulation of these genes during progression of LUAD.

The genetic alterations of these candidate genes were then determined using cBioPortal. Interestingly, seven out of the eight genes were found to be amplificated in LUAD, while only half of them were reported to be partially mutated (Figures 9A,B). Gene co-expression analysis of the eight candidate genes was also performed using cBioPortal database. As shown in Figure 9C, highly positive co-expression was observed in pairs of PA2G4–NOP58, HNRNPK–SRSF1, and NOP58–SRSF1, while significantly negative co-expression was revealed in pairs of HNRNPK–ARPC5 and SRSF1–ARPC5.

To investigate the biological characteristics of right candidate genes associated with LUAD recurrence, the GSEA assay was further applied. The results in Figure 10 showed that upregulation of these candidate genes was enriched in multiple cellular processed such as Fc gamma R-mediated phagocytosis, regulation of actin cytoskeleton, spliceosome, and tight junction.

In total, all the above findings indicated that all eight candidate genes may serve as the hub genes associated with tumor recurrence in LUAD.

To enhance the significance of our study, we next sought to establish the transcriptional regulatory network of hub genes and TFs by a Cytoscape plug-in iRegulon. As revealed in Figure 11, a total of 57 TFs and eight hub genes were involved in this network, such as models of MYC–NOP58, MYC–WDR12, MYC–SRSF1, MYC–PA2G4, MYC–HNRNPK, PML–HNRNPK, TAF1–HNRNPK, YY1–HNRNPK, and TAF1–HNRNPK. To facilitate the future targeted drug screening, Comparative Toxicogenomics Database was used to search for drugs specialized in targeting the eight hub genes. Drugs that may interact with at least two hub genes and have been reported to have an anticancer pharmacological activity were selected. Finally, 11 drugs including arsenic trioxide, cisplatin, copper, ICG 001, Jinfukang, tretinoin, doxorubicin, sodium selenite, quercetin, sunitinib, and epigallocatechin gallate were discovered (Figure 12).