Co-expression network Construction
Through conjoint analysis of co-expression network and clinical features, modules with biological significance can be confirmed in this step [15]. When the value of was chosen as 6, the scale-free topology fit index reached 0.98, meeting the standard of approximate scale-free topology (Fig. S1). LncRNAs with similar expression patterns were divided into the same module by cluster dendrogram trees and eight modules were obtained here. The results of the investigation into the relationships between modules and glioma WHO grading were presented in Fig. 1A. By setting the threshold values as Pearson correlation coefficient > |0.5| and p < 0.01 to select the significant modules, and by inquiring into mean gene significance across all genes in one module, the brown, green and yellow module were considered to be closely related to the WHO grade of glioma (Fig. 1B). Meanwhile, the module membership (MM) vs. gene significance (GS) analysis of the brown, green and yellow modules showed that the three were endowed with higher correlation between MM and GS (Fig. 1C-E). Within these three modules, we obtained 14 hub lncRNAs in total, by selecting lncRNAs with |GS|>0.85 and |MM|>0.65.
Survival analysis of risk score and clinical features
Survival analysis had been firstly performed on the training set of 181 glioma patients. From the univariate Cox analysis, all the 14 hub lncRNAs detected in the step of WGCNA were considered to be statistically connected with the clinical outcomes of gliomas (Fig. 2). Subsequently, the risk score system based on the multivariate Cox analysis was constructed, and the following formula was adopted here: risk score = CYTOR * 0.3447 + MIR155HG * -0.8509 + LINC00641 * -0.7135 + AC120036.4 * -0.5505 + PWAR6 * -0.5748. The median risk score of samples was calculated in the training set, and set as a cut-off in risk stratification for glioma patients (Fig. 3A, 3C). The same analysis was conducted for the validation dataset in order to validate the training set results (Fig. 3B, 3D). And the expression profiles of the five lncRNAs were visualized by heat map in the training (Fig. 3E) and validation (Fig. 3F) datasets. Moreover, we used the K-M curve evaluated whether there were significant survival time differences for glioma patients after risk stratification. The results showed the mortality of patients in the high-risk group was significantly higher than that in the low-risk group (Fig. 4A, 4B). The values of AUC of the signature at 1-, 3-, and 5-year OS were 0.84, 0.92 and 0.9 in the training set and 0.8, 0.85 and 0.77 in the validation cohort, demonstrating the great reliability of the prognosis signature (Fig. 4C, 4D). In the training set, univariate Cox regression suggested that risk score, age, 1p19q status and WHO grading had prognostic values (p < 0.05), while gender not. Then, in multivariate Cox regression analysis of risk score and clinicopathological variables, risk score was still an independent and powerful prognosis-predicting factor (HR = 2.002, 95%CI [1.584-2.530], p <0.001). The similar results were obtained in the validation set (HR = 1.243, 95%CI [1.063-1.469], p = 0.007) (Table 1).
Nomogram construction and Accuracy Assessment
To facilitate the clinical prognosis assessment for glioma patients, we established a nomogram to perform the prediction of the overall survival probability at 1-, 3-, and 5-year in the training cohort (Fig. 5A). The observed and predicted probabilities for the specific observations decrease along the diagonal line in calibration plot (Fig. 5B-D). To measure the predictive accuracy of the merged nomogram and the clinical risk factors at 1-, 3- and 5-year OS, the values of the area under the ROC curve were calculated. The 1-year AUC was 0.87 for nomogram, and 0.37 for age, 0.80 for glioma WHO grading, 0.38 for 1p19q status (Fig. 5E). In assessing the predicting efficacy at 3-year OS, AUC was 0.94 for nomograph, and 0.38 for age, 0.89 for glioma WHO grading, 0.37 for 1p19q status (Fig. 5F). Furthermore, the 5-year AUC was 0.93 for alignment chart, and 0.38 for age, 0.86 for glioma WHO grading, 0.33 for 1p19q status (Fig. 5G).
Differential expression analysis of hub lncRNAs
We confirmed that the 5 hub lncRNAs were expression-dysregulated between LGG and GBM tissues by t-test in the training and validation sets. In comparison with expression levels of lncRNAs in low grade glioma, lnc-CYTOR and MIR155HG were significantly over-expressed in GBM tissues, while three lncRNAs (LINC00641, AC120036.4, PWAR6) were remarkably down-regulated in GBM samples (Fig. S2).
Construction of ceRNA Regulation Network and Enrichment Analysis
A competing endogenous RNA (ceRNA) regulation network was constructed and visualized by Cytoscape software (Fig. 6A). And in the network, there were 238 nodes and 512 edges. Moreover, we explored the possible biological mechanisms of the hub lncRNAs related to glioma. By applying R package clusterProfiler, mRNAs in the network were subject to function annotation analysis. It indicated that the mRNAs regulated by the hub lncRNAs were mainly involved in regulation of cell cycle process, regulation of cellular senescence and in cell−matrix adhesion. Moreover, pathways analysis of target mRNAs showed a statistically significant association with glioma and p53 signaling pathway, confirming potential roles in cancer development of the hub lncRNAs (Fig. 6B, 6C).