Identification and Validation of a Seizure-Free-Related Gene Signature for Predicting Poor Prognosis in Lower-Grade Gliomas

Introduction

Lower-grade gliomas (LGGs), comprising grade II and grade III diffuse gliomas, represent an important class of malignant primary brain tumors with high morbidity and mortality worldwide.¹ The poor prognosis has predominantly been determined by lack of early symptoms and high rate of relapse. Epileptic seizures, often referred to simply as seizures, are common in patients with LGGs.² LGGs are at the highest risk for developing seizures among the primary brain tumors. The incidence of seizures is approximately 65–90% in LGGs.³ In most cases, seizures represent as an initial symptom and prompt diagnosis. Interestingly, many studies have shown that seizures predict more favorable outcomes in patients of LGGs. Two recent studies revealed that seizures at presentation had better survival outcomes than absence of seizures in LGGs.^4,5 A meta-analysis confirmed that LGGs patients presented free of seizures had a worse overall survival (OS) than those presented seizures.⁶ The reasons for this survival gap between seizures and seizures-free patients are considered to be multi-factorial and still poorly understood. A known explanation is that LGGs patients presenting with seizures can prompt early detection of this disease.⁷ In contrast, several studies have suggested that one appealing mechanism behind the survival gap might be isocitrate dehydrogenase 1 (IDH1) mutations or synaptic and electrical integration in neural circuits rather than early detection prompted by seizures.^8–10 Besides, it was reported that gliomas located supratentorially and cortically had a greater risk of seizures.¹¹ A meta-analysis suggested that frontal gliomas were associated with a higher risk of preoperative seizures, while gliomas in the occipital lobe were associated with a lower seizure risk.¹²

Over the past decades, the development of molecular biology has prompted our ability of early diagnosis and prognosis assessment of LGGs. A number of gene expression signatures have been identified to predict the survival of LGGs patients.^1,13,14 But it is still paramount to find more powerful signatures as all the existing prognostic signatures have been shown to offer only a marginal clinical utility.¹⁵ Recently, an increasing number of studies have showed that the combination of differentially expressed genes (DEGs) and clinicopathological parameters result in more effective prognostic evaluation of cancers.¹⁶ However, no prognostic signatures have been established in LGGs patients based on DEGs associated with seizures-free yet. In this study, we constructed a novel signature associated with seizures-free in survival risk stratification of LGGs patients. With this risk model, high-risk patients could be picked out for more effective and personalized treatment strategies.

Materials and Methods

Databases

The level 3 RNA sequencing and clinical data for LGGs patients in The Cancer Genome Atlas (TCGA) (http://tcga-data.nci.nih.gov/tcga/) database were obtained using R software and TCGAbiolinks package.¹⁷ Subsequently, we evaluated the eligibility of clinical data. We included patients aged 18 years or older who had a lower-grade (WHO grade II and III) glioma (astrocytoma, oligoastrocytoma, or oligodendroglioma). Also, only individuals with complete clinical information, such as age, gender, tumor histological grade, OS, radiation and themotherapy were included in the study. Overall, 272 individuals from TCGA database were included in the study. All samples obtained from the TCGA platform were used as the training dataset. A RNA-seq expression dataset including 157 patients with LGGs was downloaded from the Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn/) database and used as the validation dataset. Due to the open access to RNA-Seq expression and clinical data in TCGA and CGGA, additional approval from local ethics committee was not needed for this study.

Identification of DEGs Between Seizures and Seizures-Free

The publicly accessible TCGA RNA-Seq (level 3) data is preprocessed in the RESM normalization method, so that the initial normalization is not needed before differential expression analysis. EdgeR package of the R platform was applied to identify DEGs between seizures and seizures-free group of LGGs patients. In accordance with previous studies, the differential expression thresholds of DEGs were set as |log₂ fold change (FC)| >0.5 and P-value <0.05.^18,19 The same thresholds were also used to identify DEGs for breast cancer and ovarian cancer.^20,21 All DEGs were then enrolled for the construction of gene signature.

Construction of the Gene Signature

Least absolute shrinkage and selection operator (Lasso) and multivariate regression analyses are commonly applied to select prognostic genes in various conditions. Lasso regression model selects candidate prognostic genes by estimating the likelihood deviance and the ideal coefficient. In R platform, Lasso regression method is implemented by glmnet package. First, all DEGs were subjected to Lasso approach for preliminary selection of prognostic genes. Next, multivariate regression was employed to select the most prognostic genes. The same approach was used to identify gene signatures for LGGs and endometrial carcinoma.^22,23 And genes were considered significant with a cut-off of P-value <0.05. The prognostic gene signature was computed as follows: risk score = (βgene₁ × expression value of gene₁) +(βgene₂ × expression value of gene₂) +(βgene₃ × expression value of gene₃)+ … +(βgene_i × expression value of gene_i). βgene_i stood for regression coefficient of gene_i.

Assessment and Validation the Gene Signature

The prognostic performance of this gene signature in LGGs OS was evaluated by time-dependent receiver operating characteristic (ROC) through the timeROC package. According to the median risk score, 272 LGGs patients were dichotomized into high- and low-risk groups. Survival analysis of high- and low-risk groups was carried out using the survival and survivalROC packages. The prognostic power of this risk model was subsequently validated in an external dataset from CGGA.

Independent Prognostic Value of the Gene Signature

Moreover, Cox proportional hazards regression analysis was conducted to estimate hazard ratios (HRs) and their 95% confidence intervals (CIs) for potential prognostic factors such as age, gender, radiation therapy, risk score, and so forth. Univariate Cox proportional regression analysis was used to screen out the relevant factors affecting the OS. Multivariate Cox proportional hazards regression model was used to analyze the independent prognostic factors influencing the OS of patients. Only factors with statistical significance in univariate regression analysis will be further subjected to multivariate Cox proportional hazard regression model. Cox proportional hazards regression models were tested using the coxph() function of the R survival package. P-value <0.05 was considered as statistically significant.

Bioinformatics Analysis

DEGs between high- and low-risk patients were identified using R edgeR package. The differential expression thresholds of DEGs were set as |log₂ fold change (FC)| >0.5 and P-value <0.05. To further analyze biological processes of DEGs, bioinformatics analysis including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was carried out via the Database for Annotation, Visualization, and Integration Discovery (DAVID, https://david.ncifcrf.gov). P-value <0.05 and enriched gene count >2 were considered statistically enriched in GO terms and the KEGG pathways.

Statistical Analysis

All statistical analyses were conducted by R software and SPSS v.23.0 software (IBM Corp.). Student’s t-test for independent samples was performed to assess the significant differences between the two groups. Spearman correlation analysis was used for the correlation analysis.

Results

Identification of DEGs Between Seizures and Seizures-Free LGGs Patients

Two hundred and seventy-two LGGs patients were divided into the seizures (n = 223) and seizures-free (n = 49) groups according to their first presenting symptom. Table 1 summarizes patient demographics and clinical characteristics in TGGA cohort. For covariates, age, gender, tumor grade, tumor location, loco-regional surgery, molecular therapy, and radiation therapy all yielded no statistical differences between LGGs patients with and without seizures. In total, 253 DEGs were achieved between LGGs patients presented with seizures and free of seizures (Supplementary Table 1). Among them, 223 were up-regulated and 30 were down-regulated between two groups.

Table 1 Clinical Characteristics and Therapy Information in TCGA Cohort

Construction of the Prognostic Model

A total of 253 DEGs were incorporated into Lasso approach for preliminary selection of prognostic genes. As a result, 24 candidate prognostic genes were obtained through Lasso regression model (Figure 1A and B) and the name and coeff-value of each gene are shown in Figure 1C. As shown in Figure 1D, the area under the ROC curve (AUC) is 0.8343, which means that the prognostic effect of the whole 24 genes is 83.43%. Following proportional regression model, five genes (HIST1H4F, HORMAD2, LILRA3, PRSS33, and TBX20 gene) were ultimately identified as prognostic signature. The linear prediction model was characterized on the basis of the regression coefficient of each prognostic gene in the multivariate analysis. The risk score was imputed as follows: risk score = (−3.6638 × expression value of HIST1H4F) + (0.5897 × expression value of HORMAD2) + (0.5533 × expression value of LILRA3) + (0.2498 × expression value of PRSS33) + (−1.55280 × expression value of TBX20).

Figure 1 The preliminary selection of candidate prognostic genes by Lasso methods. (A) Ten-time cross-validation for tuning parameter selection in the Lasso model. (B) Partial likelihood deviation map with Lasso approach. (C) Coeff-values of the preliminary selected 24 genes. (D) ROC curve for the preliminary selected 24 genes. Abbreviations: Lasso, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Assessment of the Prognostic Model

According to the median risk score, 272 LGGs patients were divided into high- (n = 136) and low-risk (n = 136) groups. Kaplan–Meier analysis revealed that high-risk patients were linked to worse OS outcome compared with low-risk ones (P < 0.0001; Figure 2). To evaluate the prognostic signature in different clinical parameters, we adopted stratified survival analyses. Similarly, the results on OS were in favor of low-risk patients for age (younger, P < 0.0001; older, P = 0.036), gender (male, P < 0.0001; female, P = 0.0046), tumor grade (G2, P < 0.0001; G3, P < 0.0001), and seizures (with, P < 0.0001; without, P < 0.0001). Kaplan–Meier curves of all above stratified survival analyses are provided in Figure 3. ROC curves of the 1-, 3- and 5-year OS was constructed to reveal the predictive performance of the 5-gene signature. The AUC of the signature for 1-, 3- and 5-year OS is 0.885, 0.811, and 0.798, respectively (Figure 4).

Figure 2 Kaplan–Meier curves of the five-gene signature for high- (n = 136) and low-risk (n = 136) patients with lower-grade gliomas in TCGA database. Abbreviation: TCGA, The Cancer Genome Atlas.

Figure 3 Kaplan–Meier curves of the five-gene signature for high- (n = 136) and low-risk (n = 136) patients with lower-grade gliomas stratified by age, gender, tumor grade, and seizures status in TCGA database. Abbreviation: TCGA, The Cancer Genome Atlas.

Figure 4 ROC of the five-gene signature for predicting the 1-, 3- and 5-year survival of 272 patients with lower-grade gliomas in TCGA database. Abbreviations: TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

Validation of the Prognostic Model in CGGA Database

To verify the validity and reliability of the 5-gene signature impact on the prognosis of patients, we downloaded and analyzed an external dataset with a cohort of 157 LGGs cases from CGGA database. In this dataset, Kaplan–Meier survival showed that patients in the high-risk group had a significantly lower survival probability than those in the low-risk group (P < 0.001) (Figure 5). To further validate the accuracy of the risk prediction model, ROC curves of the 1-, 3- and 5-year OS were applied to reveal the predictive performance. Similar to that observed in TCGA database, the results showed that the AUC of the signature for 1-, 3- and 5-year OS is 0.778, 0.703, and 0.669, respectively (Figure 6).

Figure 5 Kaplan–Meier curves of the five-gene signature for high- (n = 136) and low-risk (n = 136) patients with lower-grade gliomas stratified by age, gender, tumor grade, and seizures status in CGGA database. Abbreviation: CGGA, the Chinese Glioma Genome Atlas.

Figure 6 ROC of the five-gene signature for predicting the 1-, 3- and 5-year survival of 157 patients with lower-grade gliomas in CGGA database. Abbreviations: CGGA, the Chinese Glioma Genome Atlas; ROC, receiver operating characteristic.

Independent Prognostic Value of the Signature

To determine the independent prognostic value of the five-gene signature for LGGs patients, univariate and multivariate regression analyses were conducted in TCGA database. Univariate proportional risk analysis demonstrated that the 5-gene signature (HR = 1.132, 95% CI 1.087 −1.179, P < 0.001), age (HR = 1.055, 95% CI 1.034 −1.076, P < 0.001), tumor grade (HR = 3.009, 95% CI 1.772 −5.107, P < 0.001), and radiation therapy (HR = 0.452, 95% CI 0.252 −0.813, P = 0.008) were all significantly correlated with the OS of LGGs patients. Differentiation remained statistically significant when entered in multivariate proportional risk regression with the five-gene signature (HR = 1.085, 95% CI 1.034 −1.134, P < 0.001), except radiation therapy (Table 2).

Table 2 Cox Regression Analysis of Clinical Pathologic Features for OS in TGGA Cohort

Bioinformatics Analysis

To further analyze biological processes of DEGs between high- and low-risk patients, GO and KEGG pathway enrichment analyses were performed via the online biological tool DAVID. As shown in Figure 7, these DEGs were mainly enriched in biological processes related to positive regulation of transcription from RNA polymerase II promoter, G-protein coupled receptor signaling pathway, and pathways of cytokine–cytokine receptor interaction, chemokine signaling pathway, and so forth.

Figure 7 Functional enrichment analysis of DEGs between different seizures presentations in lower-grade gliomas. (A) GO enrichment analysis of DEGs between different seizures presentations in lower-grade gliomas. (B) KEGG pathways of DEGs between different seizures presentations in lower-grade gliomas. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Discussion

Seizures refer to symptom related to abnormal excessive electrical activity within the brain. Seizures are frequently the first symptoms of a brain tumor, especially of LGGs. Interestingly, there is a seizures survival advantage in patients with LGGs. Scott and Gibberd first reported the prognostic value of seizures at presentation in gliomas in 1980.²⁴ Seizures have received increasing attention ever since. A meta-analysis with 2088 individuals confirmed the association between seizures at presentation and better survival outcomes in patients with gliomas.⁶ The reasons for this survival gap between seizures and seizures-free patients are considered to be multi-factorial and still poorly understood. Tumor location, early diagnosis, and IDH mutations might all contribute to the better survival outcomes of seizures patients with LGGs.^5,7,12 In addition, exploring the pathogenesis of LGGs and development of more effective prognostic markers are still hot topics in the branch of medicine. However, prognostic signatures based on different seizures presentations have not been identified in patients with LGGs till now.

To our knowledge, this study first constructed prognostic gene signature in LGGs based on different seizures presentations. Considering the combination of DEGs between different seizures presentations result in more effective prognostic signature in LGGs, we identified a total of 253 DEGs in patients with LGGs. Moreover, to gain a deep understanding of biological functions of DEGs between high- and low-risk patients, bioinformatics analysis of GO and KEGG pathways was employed. It revealed that these DEGs were mainly enriched in biological processes related to positive regulation of transcription from RNA polymerase II promoter, G-protein coupled receptor signaling pathway, and pathways of cytokine–cytokine receptor interaction, chemokine signaling pathway, and so forth. Coincidentally, Xue et al conducted a bioinformatics analysis of differentially methylated genes (DMGs) in gliomas and found that DMGs were also enriched in biological processes of positive regulation of transcription from RNA polymerase II promoter and immune response.²⁵ Chemokine and its receptor expression by tumor cells contribute to tumor growth and angiogenesis.²⁶ These biological processes and signaling pathways may contribute to the high risk of gliomas patients.

Lasso and multivariate Cox model is a popular tool to generate prognostic genes in the context of high-dimensional data. Compared to a single-gene biomarker, signatures composed of multiple genes are believed to be more robust and more accurate. With the help of R software and related packages, the prognostic signature was constructed in our study which consisted of five genes, including HIST1H4F, HORMAD2, LILRA3, PRSS33, and TBX20. As expected, Kaplan–Meier analysis found that the five-gene signature model could differentiate the survival outcome of low- and high-risk group of patients. The risk score was confirmed to be effective in different age groups, gender, tumor grade, and seizures status (with vs without). Furthermore, the AUC value of 1-, 3- and 5-year was greater than 0.79, indicating good performance of this risk model in predicting the survival. In addition, a multivariate Cox regression model analysis involving the risk score, age, gender, seizures status, tumor location, loco-regional surgery, radiation and molecular therapy demonstrated that the five-gene signature was an independent prognostic factor. In the TCGA and CGGA cohorts, there was also a prognostic difference between high- and low-risk patients.

The five-gene signature also provided new potential molecular treatment targets of LGGs. Histones are major essential components of chromatin and conserved in eukaryotic cells.²⁷ HIST1H4F, as one of the histone genes, was randomly distributed in human genome. Recent evidence has demonstrated that HIST1H4F was universally hypermethylated in patients with different cancer types, which may act as promising biomarker in screening for early cancer diagnosis in clinical applications.²⁸ For HORMAD2, it was a conserved meiotic chromosomal protein-coding gene in many organisms.²⁹ HORMAD2 was expressed at high levels in the lung adenocarcinoma tissue compared to the adjacent normal lung tissues and detected more frequently in early-stage lung adenocarcinoma, suggesting that HORMAD2 may be a candidate lung adenocarcinoma risk gene.³⁰ To some extent, this is consistent with our finding that high expression of HORMAD2 increased the risk of LGGs. There is a paucity of studies investigating roles of LILRA3, PRSS33, and TBX20 in LGGs type. Further research is needed to elucidate their role in this disorder.

Inevitably, there were several innate limitations towards the present study. First, though the consistence of our results with published data has suggested the validity of our results, it is necessary to be validated in a prospective clinical trial. Second, due to the IDH information was not available in the TCGA-LGGs dataset, the prognostic value of seizures in IDH mutations LGGs was not evaluated in our study. Further prospective studies are necessary to investigate if there is a seizures survival advantage in IDH mutations LGGs. Despite these drawbacks, our findings demonstrated that the five-gene signature could serve as reliable prognostic predictor of LGGs patients.

Conclusion

A five-gene signature was established as a novel biomarker for prognosis and provides guidance for the therapeutic strategy for these patients.