Identification and validation of a three-gene signature as a candidate prognostic biomarker for lower grade glioma

Dataset of patients with LGG

The LGGs RNA sequencing (RNAseq) data and corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA) hub by the University of California, Santa Cruz, Xena browser (https://xenabrowser.net/hub/) and CGGA data portal (http://www.cgga.org.cn/) respectively. The TCGA RNAseq data (level 3) shows the gene-level transcription estimates, as in log2(x + 1) transformed RSEM normalized count. The CGGA data displays the gene expression level as fragments per kilobase transcriptome per million fragments (FPKM), which has been standardized. Expressed gene defined only if its expressed level is larger than zero at half of samples. Only patients with a clear information of survival and detailed history of radiotherapy and chemotherapy/molecular therapy were included in the study. Finally, 456 cases from TCGA dataset and 159 cases from the CGGA dataset were included in the training set and validation set respectively. Table 1 summarized the clinical characteristics and therapy information of the training set and validation set. The workflow presentation of this study is shown in Fig. 1.

Identification of survival-related genes and construction of the prognostic model

By using the rbsurv package in R, a robust likelihood-based survival model was conducted to identify survival-related genes (Cho et al., 2009). The rbsurv package is a software program, which selects survival-associated genes based on the partial likelihood of the Cox model and adopts a cross-validation approach for robustness. According to the description of the rbsurv package, prior gene selection such as univariate survival modeling can be performed if necessary and the univariate survival modeling can be performed in this software program. Compared to the survival modeling without an adjustment of risk factors, the robust likelihood-based survival model can improve the ability to discover truly survival-associated genes by modeling genes after adjusting for certain risk factor. Thus, we directly conduct a robust likelihood-based survival model to screen for the prognostic genes. The robustness test was performed on 20,530 genes and 456 samples. After 10 iterations, 29 prognostic related genes were selected. With the help of glmnet and survival package in R, least absolute shrinkage and selection operator (LASSO) regression and the multivariable Cox proportional hazard regression method were used to further identify the survival-related prognostic model. The same approach was used to identify gene signatures for endometrial carcinoma (Ouyang et al., 2019). At last, three prognostic survival-related genes that were independent survival predictors and their regression coefficients were obtained at a threshold of p < 0.05. Based on the median expression value of each survival-related gene, we dichotomized 456 LGGs patients into low and high expression groups and compare the survival rate between the two groups by Kaplan–Meier plots and Log-rank test. According to the estimated regression coefficients, a prognostic risk score for each patients was then calculated (Wang et al., 2019). The risk score = (0.4470 × expression level of WEE1) + (−0.1530 × expression level of CRTAC1) + (−0.3723 × expression level of SEMA4G). With the three-gene signature, 456 LGGs patients were divided into high-risk and low-risk groups with the median risk score as the cut-off value. Kaplan–Meier curves were performed to estimate and compare the survival for TCGA LGGs patients with a high score or a low score. The receiver operating characteristic (ROC) curve and area under the curve (AUC) were applied to evaluate the prediction accuracy of the risk score model. Furthermore stratified survival analysis was performed in patients with different age group (younger, old), gender (male, female) and pathologic grade (G2, G3).

Univariate and multivariate Cox hazard regression analysis were conducted for the potential prognostic factors such as age group (younger vs. old), gender (male vs. female), pathologic grade (G2 vs. G3), radiation therapy (Yes vs. No), molecular therapy (Yes vs. No) and risk score (High vs. Low).

Validation of the prognostic model in the CGGA

The prognostic model was validated in the CGGA mRNAseq_325 cohort. Only patients with a clear information of survival, detailed history of radiotherapy and chemotherapy were included in the study. Finally, 159 cases from the CGGA mRNAseq_325 cohort were included in the validation set.

Exploring co-expression genes by WGCNA

To explore the regulatory network of the three genes, Gene Co-expression Network Analysis (WGCNA) was performed in training set by the R package WGCNA (Langfelder & Horvath, 2008). The top 50% variance of genes were selected for WGCNA. In other words, WGCNA based on 456 samples and 10,256 genes. First, RNAseq data were filtered to reduce outliers. Using the absolute value of the correlation between the expression levels of transcripts, a co-expression similarity matrix was constructed. Then, the co-expression similarity matrix was transformed to the adjacency matrix by choosing nine as a soft threshold. Co-expression gene module was established by the topological overlap measure. In order to identify the significance of each module, gene significance (GS) was calculated to estimate the correlation between genes and sample traits. Module significance (MS) was defined as the average GS within modules and was calculated to measure the correlation between modules and sample traits (vital status). Finally, the “vital status” related modules that contain the three genes as members and genes belong to such modules were identified. Genes interacted with those three genes were screened and the co-expression network was constructed by Cytoscape software (Shannon et al., 2003).

Functional enrichment analysis

Using Enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were conducted via the clusterProfiler package in R language (Yu et al., 2012) for those genes that belong to the “vital status” related modules associated with the three genes. Benjamini–Hochberg (BH)-adjusted p-value < 0.05 were considered significant.

Three prognostic genes were identified in TCGA dataset and validated in CGGA dataset

A total of 456 patients and 20,530 genes were included in the TCGA-LGG to train the prognostic model. The robust likelihood-based survival model found 29 survival-related genes, 13 genes were obtained through LASSO Cox method (Fig. 2). We further reduced the dimensionality of these high-dimensional data by multivariate Cox proportional hazard regression model. Finally, three genes that were independent survival predictors were identified as survival prediction signature. Those three genes included in the model were WEE1, SEMA4G, CRTAC1. It has been reported that WEE1 is closely related to the growth, invasion and migration of glioma (Wu et al., 2019). Currently, there is no study revealing the role of SEMA4G and CRTAC1 in gliomas. After calculating the risk score, patients were divided into a high- and low-risk group based on the median cut-off point of the risk score. The three-gene signature risk score distribution is shown in Fig. 3A. Besides, the relationship between risk score and the status of the LGGs was calculated (Fig. 3B). As shown in the heat map of the Fig. 3C, a remarkable high expression was noted for WEE1 in the high-risk group, while a lower expression was observed for the other genes in the high-risk group (Fig. 3C). Patients in the high-risk group were significantly worse off the overall survival time compared to the low-risk group (p < 0.0001) (Fig. 4A). The area under ROC curve of the signature for 1-, 3- and 5-year overall survival was 0.908, 0.878 and 0.827, respectively, in training set. (Fig. 4B). A similar result can be noted in the validation dataset (Fig. 4C). The area under ROC curve of the signature for 1-, 3- and 5-year overall survival was 0.781, 0.806 and 0.807, respectively, in validation set. (Fig. 4D). Moreover, the predicting power of the risk score model was not decreased in subgroup analysis for age group (younger, p = 0.00012; old, p < 0.0001), gender (male, p < 0.0001; female, p < 0.0001), and pathologic grade (G2, p = 0.00013; G3, p < 0.0001) in the training set (Figs. 5A–5F). The same trend can be observed in the validation dataset (Figs. 6A–6F). For the WEE1, the member of high expression group had significantly shorter survival than those in low expression group (p < 0.0001) (Fig. 7A). For the SEMA4G and CRTAC1, the member of high expression group had significantly longer survival than those in low expression group (p < 0.0001) (Figs. 7B and 7C). The expression level of WEE1 was significantly higher in grade III compared to grade II (p < 0.0001), while the other are opposite (Fig. 8A). These results can also be verified in the validation dataset (Figs. 7D–7F and 8B).

Multivariate Cox proportional hazard regression demonstrated that age group (HR = 0.274, p = 2.21E−09), pathologic grade (HR = 2.49, p = 0.00011) and risk score (HR = 0.198, p < 0.000000000168) were independent prognostic factors in the training set, while pathologic grade (HR = 3.799, p = 0.00000151), 1p19q status (HR = 4.566, p = 0.0000388), radiation therapy (HR = 0.524, p = 0.046), and risk score (HR = 0.415, p = 0.000653) were independent prognostic factors in validation dataset (Table 2).

Calculation of module-trait correlation in LGGs and module visualization of the network connections

Using the R package WGCNA, gene modules were identified based on the top 50% variance of genes. To analyze the relationship between gene modules and sample clinical information, we used the module eigengene (ME) as the overall gene expression level of the corresponding modules and calculated correlations with clinical phenotypes, for example, vital status. We obtained 16 gene modules (Figs. S1A–S1D) with size ranging from 31 to 1,501 genes. We assigned each co-expression module an arbitrary color for reference: black, blue, brown, cyan, green, greenyellow, lightcyan, magenta, midnightblue, pink, purple, red, salmon, tan, turquoise and yellow. These modules contained 449, 1,352, 850, 46, 519, 91, 31, 201, 43, 336, 135, 462, 51, 90, 1,501 and 845 genes, respectively. As a single group, the non-co-expressed group designated as “grey” based on the WGCNA developer’s rationale. Vital status related modules, such as yellow, green, black modules that contain the three genes as members and genes belong to such modules were screened (Fig. S1D). Finally, 32 genes were discovered to be co-expressed with CRTAC1, 181 genes were co-expressed with WEE1, six genes with SEMA4G. We exported the screened genes and three prognostic survival-related genes into Cytoscape and constructed the co-expression network (Fig. 9).

GO and KEGG analysis of screened genes interacted with three-gene signature

For the “biological processes” (BP), chromosome segregation, nuclear division, mitotic nuclear division, organelle fission, mitotic sister chromatid segregation were the commonly enriched categories (Fig. 10A). For the “cellular component” (CC), the enriched categories were correlated with condensed chromosome, chromosome/centromeric region, chromosomal region, kinetochore, condensed chromosome/centromeric region (Fig. 10B). For the “molecular function” (MF), these screened genes mainly enriched in microtubule binding, tubulin binding, histone kinase activity, DNA-dependent ATPase activity, protein serine/threonine kinase activity (Fig. 10D). KEGG pathway enrichment analysis suggested that cell cycle was the most important pathway for these selected genes. The following pathway also involved many screened genes, including, oocyte meiosis, progesterone-mediated oocyte maturation, fanconi anemia pathway, homologous recombination (Fig. 10C). Additionally, For the GO analysis, these three co-expression gene modules (yellow, green, black) enriched results can been seen in Figs. S2 and S3.

A comparison between our and other models

Recently, Zeng et al. (2018) reported a model containing three genes (EMP3, GSX2, EMILIN3) based on integrative analysis of DNA methylation and gene expression in TCGA dataset. Zhang et al. (2019a) also reported a 4-gene (EMP3, GNG12, KIF2C, IFI44) prognostic signature based on genes encodes by chr1p/19q. To compare the prognostic values of our prognostic signature and their model, we performed time-dependent ROC curve analysis in our model and other models based on the risk score calculated by the regression coefficients which obtained by themselves and the expression level of members in their signature were shown in both TCGA and CGGA dataset. The results exhibited that our model displayed a satisfactory predictive value in predicting 1-, 3- and 5-year overall survival compared to other models (Figs. 11A–11F). In other words, our 3-gene model had a satisfactory efficiency in predicting both short- and long-term prognosis.