We collected gene expression profiling data from past LUAD patients that were included in the TCGA database so that 513 lung adenocarcinoma and 56 normal tissue samples were utilized in our study (https://portal.gdc.cancer.gov/). Corresponding patient data, such as age, gender, living status, tumor status, TNM stage, radio- and chemotherapy and survival data, and were available from TCGA. During data preprocessing we transformed gene names into official gene symbols using Perl language, and genes needed to have a non-zero expression value in a minimum of half of the sample type to be kept.

Differentially expressed analysis was conducted by the package “limma”, and the cut-off criteria was set as |log (Fold Change)| > 1 and p-value < 0.05 (Pencina and D’agostino, 2004). The R package “pheatmap” was used to draw heatmap and volcano plot. The box plots of the key genes for validation were plotted by R, using the package “ggpubr”.

We constructed co-expression modules using the WGCNA R package based off the DEGs we identified between normal and tumor tissues in LUAD patients. WGCNA was able to group highly correlated genes and identify important modules or genes that are associated with clinical characteristics of interest. We also used Pearson correlation coefficients between each identified gene module to construct a matrix than can establish module-trait relationships between DEG expression and the associated mRNAsi with respect to the β value (soft threshold value). After the most important module was identified we then calculated the gene significance (GS) and module membership (MM). We defined key genes as those with the GS > 0.7 and MM > 0.7 for this module.

We used Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis with the “clusterProfiler” R package to explore potential gene pathways and biological functions in the significant module. Statistical significance was defined as a P. adjust value < 0.05. The R package “org.Hs.eg.db” was used to map important genes with the Ensemble ID. Bubble plots were assembled by R to visualize the top 10 results.

The protein-protein interaction (PPI) networks were constructed based on genes in the most mRNAsi-related module, using the STRING database (version 11.0) of known and predicted protein-protein interactions (Ritchie et al., 2015), which now covers 24,584,628 proteins from 5,090 organisms. The users only need to submit a list of gene symbols and species, and the website provides interaction relationships among submitted proteins. These interactions include direct (physical) and indirect (functional) associations.

To study the function of hub genes in the key module, we separated 513 LUAD patients into different subgroups by the approach of “ConsensusClusterPlus”, an algorithm for determining clusters using an unsupervised analysis based on gene expression (Szklarczyk et al., 2019). The consensus clustering tool provides users with quantitative and qualitative evidence that allows the estimation of unsupervised class counts in a dataset. The maximum evaluated k (max K) was set to 9 and other parameters of ConsensusClusterPlus were set to default. The R package (R v3.5.1) of PCA analysis was adopted to explore the gene expression patterns in subgroups of LUAD.

To explore the prognostic value of hub genes in hub module resulting from WGCNA, we used a least absolute shrinkage and selection operator (LASSO) regression to narrow the range of target genes via the “glmnet” R package because the predictor variable was much larger than the sample content in the gene expression data. The risk score formula for predicting the prognosis of LUAD patients was: risk score = the sum of the multivariate LASSO regression coefficient ratio of each mRNA multiplied by the expression level of each mRNA. We divided the LUAD patients into two groups, high- and low-risk, based on the median value of risk score. In the two groups, the clinicopathological characteristics of each patient, including age, gender, living status, tumor status, TNM stage, survival status, chemotherapy, radiotherapy, and gene expression profile, were presented via the “pheatmap” and “survival” R packages. In addition, we conducted time-dependent receiver operating characteristic (tdROC) curve and Kaplan-Meier survival curve analyses to validate the signature in both the training set and testing set. Log-rank test was applied to calculate the difference of overall survival rate between the high-risk and low-risk groups. “p < 0.05” was considered statistically significant. In the validation phase, we verified the nomogram in the GEO by using another LUAD cohort, GSE17536 and GSE17537.

To validate the independence of the risk model, we conducted univariate and multivariate Cox regression analyses to evaluate the predictive efficacy of the model. A nomogram encompassing the risk score model and clinicopathological factors was plotted by the “rms” R package. The accuracy of the nomogram was examined using the consistency between the actual and the predicted outcomes. Next, we submitted these outcomes to the calibration curve to visualize the performance of the nomogram. The 45° line represented the best prediction (Wilkerson and Hayes, 2010). Based on the different clinicopathological characteristics and the risk score of each patient, we calculated the total score to predict 1, 3, and 5-years prognosis of LUAD patients. To determine the prognostic value of the nomogram, we use the Kaplan–Meier survival curve to compare the predictive value of nomogram for OS. The predictive efficiency of the nomogram for 1/3/5-years survival was assessed using the tdROC curve analysis.

Differences between variables were assessed with independent t-tests. Kaplan-Meier curves and log-rank tests were used to analyze the survival data, and univariate Cox regression analysis was used to identify independent prognostic factors. Time-dependent ROC analysis was used to evaluate the accuracy of the prognostic predictive model. The area under the ROC curve (AUC) > 0.60 was regarded as acceptable for predictions, and AUC >0.75 was deemed to have great predictive value. R software was used to perform all statistical analyses, and p < 0.05 was considered statistically significant.

The flowchart shown in Figure 1 summarizes the overall bioinformatics analysis of our study. The WGCNA R package was applied to build a co-expression network, and 4,430 differentially expressed genes (DEGs) was listed in Supplementary Table S1. In this study, we set the soft threshold to the power of β = 4 to guarantee a scale-free network (Figures 2A,B). To investigate the gene sets that are associated with mRNAsi of LUAD, we applied WGCNA because it defines transcriptional modules using Pearson correlation and establishes a relationship between different colored modules and clinical traits. As a result, we identified 9 distinct co-expression modules (Figure 2C).

To detect whether modules were significantly associated with specific clinical traits, we looked for strong associations between eigengenes and external traits. The correlation between module and trait was then visualized as a heatmap (Figure 2D). These results showed that patient disease state is significantly correlated with seven modules. Based on the correlation coefficients, genes clustered in black, pink, green, brown, blue, and yellow modules are downregulated in LUAD tissues, while genes in the turquoise module are highly expressed in LUAD tissues. Genes clustered in turquoise modules have the strongest positive correlation (Cor = 0.7, p = 1e-76) with patients’ disease status. This data indicates that genes in the turquoise module are significantly associated with the mRNAsi of LUAD patients. The scatterplot below illustrates the strength of the link between the mRNAsi signature and the module membership for each gene found in blue, brown, and turquoise modules (Figures 2E–G).

Next, the turquoise module underwent GO and KEGG pathway analysis and we identified the top 10 GO and top 7 KEGG pathway enrichment terms, as shown in Figures 3A,B. The GO analysis showed that this module is enriched in organelle fission, chromosomal, and ATPase activity. The KEGG analysis indicated that the genes are mainly involved in cell cycle, cellular senescence, oocyte meiosis et al., which are pivotal in the regulation of immune responses.

As indicated in Figure 4A, we found a notable difference between the mRNAsi of LUAD and the mRNAsi of normal tissues. The mRNAsi of LUAD tissues was found to be higher compared to that of normal tissues. We also observed significant differences in N stage (Figure 4B), T stage (Figure 4C), and AJCC stage (Figure 4D). The PPI network, consisting of the top 22 hub genes in turquoise module, was constructed using the STRING database. In total, 22 nodes and 460 edges were included in this PPI network (Figure 4E), with an average node degree of 20.9 and strong correlations. The expression levels of the top 22 hub genes were higher in tumor tissue than in normal tissue (Figure 4F).

ConcensusCluster analysis was utilized to classify the tumor samples. According to the expressional similarity of the gene expression in the turquoise module identified above as the powerful independent prognostic factors, k = 3 could be the optimal choice with clustering increasing from k = 2–9 (Figures 5A,B). However, we found that the interference found between subgroups only measured as minimal when k = 2. As such, we used k = 2 for consensus clustering analysis and, though this, we identified two subgroups named cluster 1 and cluster 2 (Figures 5C,D). Patients in cluster 1 were found to have a significantly shorter overall survival (OS, Figure 5E) and recurrence-free survival (RFS, Figure 5F) when compared with patients in cluster 2.

To filter out an mRNAsi-based classifier to predict the prognosis of LUAD patients, a LASSO regression model was carried out using the 76 miRNAsi-related genes from the turquoise module. Furthermore, the nine mRNAsi survival-related DEG prognostic model (CENPW, MCM2, STIL, RACGAP1, ASPM, KIF14, ANLN, CDCA8, and PLK1) was constructed with LASSO regression to improve the predicted accuracy for overall survival in LUAD when log(lambda) was between -3.0 and −4.0 (Figures 6A,B). Risk scores were based on gene expression levels multiplied by its corresponding regression coefficient (Table 1). The formula was shown as: risk score = (ANLN*0.1574) (ASPM*0.0163) −(CDCA8*0.1315) +(CENPE*0.0241) +(KIF14*0.0598) −(MCM2*0.0431) +(PLK1*0.0907) −(RACGAP1*0.0916) −(STIL*0.0478). We then deeply analyzed the relationship among these 9 genes, and the results were shown in Figure 6C. Expression of the 9 genes and clinicopathological features in low- and high-risk patients from the TCGA dataset was demonstrated in the heatmap (Figure 6D). We found significant differences between the high- and low-risk groups associated with tumor status, stage_N, stage_T, AJCC stage, recurrence, and living status. We then used the median risk score to classify 513 LUAD samples as either low risk (n = 256) or high risk (n = 257) (Figures 6E,F). K-M survival curve indicated that patients in the high-risk group showed markedly poorer overall survival (OS) than those in the low-risk group (p = 7.852e−09; Figure 6G). The area under the ROC curve (AUC) for OS was 0.716 (Figure 6H), suggesting that this prognostic model exhibited a great sensitivity and specificity.

To determine whether the risk signature could be used as an independent risk factor for LUAD patients, we performed univariate and multivariate Cox analyses (Figure 7A). The covariates included age, gender, T stages, M stages, N stages, chemotherapy, radiotherapy, tumor status, and risk model. Univariate Cox regression analysis showed that age, stage_T, stage_N, chemotherapy, radiotherapy, and risk model correlated with the prognosis of LUAD (p < 0.05). However, subsequent multivariate Cox analysis showed that the stage_N (HR = 1.378, 95% CI = 1.150–1.652, and p < 0.001), tumor status (HR = 6.324, 95% CI = 4.290–9.323, and p < 0.001), and risk model (HR = 2.373, 95% CI = 1.659–3.394, and p < 0.001) were independent risk factors for overall survival. We constructed nomogram maps to predict 1-, 3-, and 5-years overall survival in LUAD patients with stage_N, tumor status, and risk score (Figure 7B). The corresponding score of each factor was listed as Table 2. We validated the nomogram by calibration curve, and the 45° line represented the best prediction. Calibration plots indicated that the nomogram performed well in predicting 1-, 3-, and 5-years survival (Figure 7C). Patients were then divided into 3 subgroups, equally, according to the total score resulting from the nomogram. Kaplan-Meier survival curve revealed that patients in the low score subgroup had a better clinical outcome than those in the moderate and high score subgroups (Figure 7D). ROC curve analysis exhibited that the risk score AUC values of the model were 0.797, 0.821, and 0.836 for 1-, 3-, and 5-years survival, respectively (Figure 7E). These results reveal that the nomogram constructed by the 9-gene signature has a high accuracy in predicting the overall survival of LUAD patients.

The efficacy of the classifier was further evaluated using another outcome of different type of survivals, and the results were similar to what was seen in the TCGA cohort. To validate our signature, we first calculated the risk score for each patient according to the coefficient value of the 9 genes. Patients were divided into high-risk and low-risk groups with the median risk score utilized as the cut-off value. In the GSE17536 dataset, patients in the high-risk group had poorer outcomes in disease-free survival (DFS, Figure 8A), disease specific survival (DSS, Figure 8B), and overall survival (OS, Figure 8C). Similarly, in the GSE17537 dataset, patients in the high-risk group also had a worsened prognosis in DFS (Figure 8D) and OS (Figure 8E) when compared with those in the low-risk group. Altogether, this external validation indicated that risk signature is stable and highly precise in its prediction of LUAD patient prognosis.