A 4-Gene Signature Associated with Clinical Outcome in High-Grade Gliomas

High-grade gliomas (HGG) are brain tumors associated with high morbidity and mortality. They are classified as either grade III or grade IV on the basis of histopathologic features established by the World Health Organization (WHO) (1). In combination with other clinical parameters, the grade has long provided important prognostic information (2). Recently, molecular biomarkers have been shown to be strongly associated with the prognostic of these tumors. O(6)-methylguanine-DNA-methyltransferase (MGMT) promoter hypo-methylation is involved in glioblastoma resistance to temozolomide chemotherapy (3) and mutations of the isocitrate dehydrogenase 1 (IDH1) gene are associated with better outcome of patients (4).

Recent studies have demonstrated that molecular and genetic analysis of gliomas could help in their classification and in the design of treatment protocols (5, 6). Microarray expression profiling has characterized molecular subtypes of brain tumors associated with tumor grade, progression, and prognosis (6–11), though only a few genes have been consistently identified (12). To overcome such a lack of reproducibility, the best approach is to analyze multiple data set simultaneously to combine the results from relevant studies. Such analysis applied to microarray data has been shown to be a powerful tool to identify candidate biomarkers and biological pathways (13).

The 2 most comprehensive glioma microarray classifications schemes published to date (6, 9) are based on unsupervised analysis, and they clearly show a strong association between the tumor grading and the defined glioma subtypes. These 2 classifications proposed by Phillips et al. (8) and Li et al. (6) show that HGG patients with better-than-expected survival could be classified in an enriched grade III subtype designated proneural or oligodendroglioma-rich, respectively.

Here we considered a supervised approach to account for the current WHO grading when deriving gene biomarkers associated with clinical outcomes (7, 10, 11). Following this rationale, we performed a combined analysis to identify such biomarkers from a robust signature related to tumor aggressiveness. We identified and validated a risk-score model based on the weighted expression of 4 genes: CHAF1B, PDLIM4, EDNRB, and HJURP. Two independent validations were performed: the first, on a public microarray study, and the second, on a new set of HGGs by quantitative reverse transcription-PCR (RT-qPCR). We also compared the performances of our risk-score model with the prognostic value of currently admitted clinical and molecular risk factors (MGMT promoter methylation status and IDH1 mutational status).

The local cohort totalized 194 patients with newly diagnosed and untreated HGGs admitted to the University hospitals involved in the French Cancéropôle Grand-Ouest Glioma Project (see Supplementary References). Patients were selected retrospectively in the period 1998 to 2008 with a follow-up time of minimum 2 years. Tumor samples were collected in accordance with the French regulations and the Declaration of Helsinki. All patients gave their informed consent before inclusion. Initial histology was confirmed by a central review involving at least 2 neuropathologists according to the WHO classification of central nervous system tumors (1). Patient characteristics are summarized in Table 1. Total DNAs and RNAs were isolated from frozen samples of primary brain tumors stored (−80°C) in the Cancéropôle Biological Resource Centers (AC-2008-77). Quality of DNA samples was assessed on 1% agarose gel and RNA integrity was confirmed using the Agilent 2100 Bioanalyzer (Agilent Technologies).

RT-qPCR reactions were performed as described previously (14) with B2M (β-2 microglobulin) and HPRT1 (hypoxanthine phosphoribosyltransferase) as internal controls. All the primer sequences used in this study can be found in the Supplementary Experimental Procedures.

Exon 4 of the IDH1 gene was amplified with the use of a PCR assay and sequenced in DNA from the tumor from each patient, as described previously (15). Patients were screened for somatic mutations affecting the R132 residue of IDH1.

The pyrosequencing methylation assay was performed with the PyroMark Q96 CpG MGMT kit (Qiagen), according to the manufacturer's protocol. Samples were considered methylated if they had average CpG methylation of 9% or more and unmethylated if they had average methylation less than 9%, in duplicate reactions (16).

External microarray data for 326 patients were collected from 4 Gene Expression Omnibus HGGs data sets (7, 8, 10, 17). There were 22,215 common probe sets in the 3 data sets. Base-10 log-transformed intensities were centered using the scale function of the R base package. Data sets characteristics and analysis workflow are presented in Fig. 1.

Combined analysis was performed on 267 patients (GDS1962, GSE4271, and GSE4412) using the Bioconductor RankProd package (18). This package utilizes the rank product nonparametric method to identify up- and downregulated genes between anaplastic astrocytomas and glioblastomas (19). The RankProd package was chosen for its ability to easily combine data sets from different origins (laboratories and environments) into a single analysis. It was also shown that this nonparametric method outperforms other meta-analysis methods in terms of sensitivity and specificity (18). Individual analyses were also performed for each study (2-sided Student's t test) and results were combined. Genes were considered to be differentially expressed for a corrected P value (false discovery rate, FDR) less than 0.05 and a fold change greater than 2 in at least 1 of the 2 approaches. Functional annotation analyses were assessed using the Database for Annotation, Visualisation, and Integrated Discovery (http://david.abcc.ncifcrf.gov/) and unsupervised PCA with integration of biological knowledge (20). We used Benjamini corrected P values for multiple testing (P < 0.05).

We performed a cross-study analysis of genes that can assist in the prognostication of survival by univariate Cox regression analysis. Gene expression was used as a predictor and survival time (in months) as the response. To select the significant genes, we controlled the FDR with the Benjamini–Hochberg correction and set the P value threshold at 0.01. To build an optimal survival model, we selected survival-associated genes with the rbsurv R package. Briefly, this package allows a sequential selection of genes based on the Cox proportional hazard model and on maximization of log-likelihood. To increase robustness, this package also selects survival-associated genes by repetition (1,000 times) of a separation between the training and validation sets of samples. Regression coefficients of the optimal survival model were estimated after adjustment on the study factor. Risk scores were determined using classical Cox model risk formulae with a linear combination of the gene expression values weighted by the estimated regression coefficients. Time-dependent receiver-operating characteristic (ROC) curve analyses were used to select the optimal risk cutoffs for the stratification of patients. The Kaplan–Meier method was used to estimate the survival distributions. Log-rank tests were used to test the difference between survival groups. Analyses were carried out with the survival and survivalROC R packages.

We constructed a model including clinical factors—age, treatment, histologic grade, and risk classes as defined by the Radiation Therapy Oncology Group (RTOG) by recursive partitioning analysis (RPA; ref. 21)—along with MGMT methylation status and IDH1 mutational status. We evaluated the discriminatory capability of the model with the gene expression risk score as compared with the model without the gene expression risk score using C statistics. Differences in discrimination were evaluated using a nonparametric approach (22). We assessed model calibration using the Hosmer–Lemeshow χ² test (23). Analyses were performed using the Hmisc and Design R packages.

Data sets characteristics and analysis workflow are summarized in Fig. 1.

Combined analysis and individual study approaches were performed to define a consensus gene expression signature in HGGs that could be used to find biomarkers associated with clinical outcomes. This signature was composed of 438 gene probe sets with 65 identified by both approaches (Supplementary Fig. S2; Supplementary Table S1). Associated enriched GO processes (Supplementary Tables S2 and S3) were related to invasion (Supplementary Fig. S5), angiogenesis (Supplementary Fig. S3), response to stress (Supplementary Fig. S6), and morphogenesis (Supplementary Fig. S4). Among the consensus genes strongly associated with grading, we selected and validated 9 genes (CHI3L1, ADAM12, S100A4, TIMP1, NDRG2, NTRS2, LUZP2, ALDH5A1, and RASL10A) by RT-qPCR (P < 0.001) on a subset of 90 HGG samples (Supplementary Fig. S7).

To assess the survival prognosis capabilities of the 438 selected probe sets, we performed univariate Cox analyses of the expression data for these genes, with overall survival (OS) as a dependent variable. We ranked the genes on the basis of their predictive power (univariate z score). We then selected the genes having a highly significant association with survival and identified 40 genes with high predictive power (Supplementary Fig. S8). According to the univariate z score, 26 were risk genes and 14 were protective genes. Risk genes were related to gene ontology (GO) biological process cell cycle (CDC25A, ASPM, CHAF1B, CENPE, CEP55, CDC20, NCAPG, AURKA) and to ECM-receptor interaction and Focal Adhesion KEGG pathways (HMMR, COL1A2, COL4A2, COL1A1, COL4A1, MET). Interestingly, 5 of the protective genes were related to GO biological process nervous system development (EDNRB, ABLIM1, ALDH5A1, NDRG2, FGF12).

We performed multivariate Cox regression analyses to create an optimal gene-based survival model. We used the 40 selected genes to sequentially construct survival models. The model best associated with survival (P < 0.001) and with good discrimination ability (C statistic, 0.843; 95% CI, 0.647–0.827) was based on the expression of 4 genes: CHAF1B, PDLIM4, EDNRB, and HJURP. The relative contributions of each of the 4 genes in the multivariate analysis are summarized in the portion of the Cox risk equation that captures the individual risk profile: (0.587 × CHAF1B) + (0.326 × PDLIM4) + (−0.470 × EDNRB) + (0.532 × HJURP). Patients were ranked according to their risk score. The optimal risk cutoff was assessed and used for the stratification of patients into 2 groups: low risk of death and high risk of death. Patients with a low risk of death (25 anaplastic astrocytomas and 36 glioblastomas) had a median OS of 46.6 months (95% CI, 28.7–73.9), which was significantly longer than 11.7 months (95% CI, 9.0–13.5) for patients with a high risk of death (4 anaplastic astrocytomas and 79 glioblastomas), P < 0.001 by the log-rank test (Fig. 2A).

During this work, the MD Anderson group published a 9-gene panel (AQP1, CHI3L1, EMP3, GPNMB, IGFBP2, LGALS3, OLIG2, PDPN, and RTN1) to predict outcome in glioblastoma (24). Six of these genes were also found in our consensus gene selection. We compared our 4-gene panel with the MD Anderson group 9-gene predictor. Both models were highly significant (P = 1e-08 and P = 3e-05, respectively). The discrimination of the 4-gene model was significantly higher than the discrimination of the 9-gene model [C statistic, 0.80 (95% CI, 0.72–0.86) versus 0.76 (95% CI, 0.64–0.81), P < 0.001, respectively; Supplementary Fig. S9], showing the relevance of the 4-gene panel.

We performed an external validation of the 4-gene survival model with an independent microarray data set comprising 56 HGGs with survival data reported by Petalidis et al. (7). Patients were divided into 2 groups on the basis of the 4-gene model (low or high risk of death). The low-risk group was composed of 12 anaplastic astrocytomas and 5 glioblastomas and the high-risk group of 5 anaplastic astrocytomas and 34 glioblastomas. The OS was higher in low-risk HGGs than in high-risk HGGs [17.8 months (95% CI, 9.6–47.9) versus 9.3 months (95% CI, 7.2–13.9), respectively; P < 0.001; Fig. 2B]. The discrimination was as good as in the original data (C statistic, 0.852; 95% CI, 0.673–0.933).

Model validation was also performed to determine whether the 4-gene expression data contained survival-predictive information that was distinct from the prediction embedded within histologic grade. In the whole anaplastic astrocytoma set, the OS was higher in low-risk patients (n = 9) than in high-risk patients [n = 37; 69.4 months (95% CI, 41.8 to not reached) versus 19.7 months (95% CI, 13.7 to not reached), respectively; P < 0.05; Fig. 2C]. In the whole glioblastoma set, low-risk patients (n = 34) had a much higher OS (30.07 months; 95% CI, 17.7–54.2) than high-risk patients (n = 120; 9.3 months; 95% CI, 7.6–11.7; P < 0.001; Fig. 2D).

A cohort of 194 patients with extensive bioclinical parameters was used to validate the performances of the 4-gene classifier (Table 1). Univariate analyses showed that the gene expression risk score, the DNA methylation status of the MGMT promoter, and the IDH1 mutational status were significantly associated with the OS in this cohort. In the whole cohort, patients were divided into 2 groups on the basis of the risk-score model with log2-transformed data issued from RT-qPCR analysis. The OS was clearly higher for low-risk patients (55.8 months; 95% CI, 26.0 to not reached) than for high-risk patients (14.5 months; 95% CI, 12.5–16.0; P < 0.001; Fig. 3A). In this population, MGMT-methylated tumors, compared with unmethylated tumors, had a significantly better OS [19.5 months (95% CI, 16.7–29.4) versus 14.5 months (95% CI, 11.4–16.2), respectively; P < 0.001; Fig. 3B]. Similarly, in this group, IDH1-mutated tumors had a much higher OS (median survival not reached; 95% CI, 42.5 to not reached) than IDH1-nonmutated tumors (14.9 months; 95% CI, 13.7–16.5; P < 0.001; Fig. 3C).

Two multivariate models were built, both including age, treatment, grade, RTOG RPA classes, MGMT methylation status, and IDH1 mutational status; 1 with and 1 without the 4-gene expression risk score. These models were used to estimate the prognostic value of the gene expression risk score (i) for the 176 of 194 patients with complete data for all variables and (ii) for a subset of patients treated with temozolomide chemoradiation (n = 105). Results are provided in Table 2. In both cases, the gene expression risk score was strongly associated with survival (hazard ratio = 0.49; 95% CI, 0.30–0.81; P = 0.005; and hazard ratio = 0.37; 95% CI, 0.18–0.77; P = 0.008; respectively) and all models showed excellent discrimination, with C statistics more than 0.80. In the whole cohort and for the patients treated with temozolomide chemotherapy, the C statistic improved significantly with the addition of the gene expression risk score in the model (0.816 versus 0.846, P < 0.001 and 0.792 versus 0.822, P < 0.001, respectively), showing that the 4-gene risk score added beyond standard clinical parameters and beyond both the MGMT methylation status and the IDH1 mutational status.

The performances of the gene expression risk score were also evaluated on a subset of 98 patients with glioblastoma who underwent tumor resection and who were treated with radiotherapy plus concomitant and adjuvant temozolomide. After adjustment for RTOG RPA classes and MGMT promoter methylation status, multivariate analysis confirmed that the 4-gene expression risk score was an independent marker robustly associated with outcome for glioblastoma patients treated with standard protocol (hazard ratio = 0.386; 95% CI, 0.164–0.910; P = 0.03).

Molecular studies of HGGs have highlighted the heterogeneity of these tumors, and have linked molecular signatures to their natural history and to differences in survival rates. It is likely that the ability to identify such molecular subtypes of tumors will be essential for guiding therapeutic advances. In this study we report a risk-score model based on the expression of 4 genes for the stratification of patients with HGGs. This risk calculation is based on a consensus gene expression signature and is strongly associated with survival independently from current clinical risk factors, IDH1 mutational status, and MGMT promoter methylation status.

The initial step of our study consisted in a discovery phase for the identification of biomarkers repeatedly correlated with both tumor aggressiveness and patient outcome. It should be noticed that information regarding the therapeutic regimens was not incorporated in the meta-analysis of microarray data sets. Although this could have weakened this discovery phase, combining multiple and independent data sets was also an asset to identify robust biomarkers. Moreover, the RT-qPCR validation of the 4-gene signature in an external cohort of patients showed that the 2 risk groups had significant differences in OS independently from treatment. These results suggest that the 4 genes are relevant molecular markers in HGGs.

One explanation for the association between the 4-gene signature and clinical outcome could be that it may detect the molecular fingerprints inherent to glioma aggressiveness. The proposed multimarker panel is based on the expression of EDNRB, CHAF1B, PDLIM4, and HJURP. In this model, the overexpression of EDNRB correlates with better prognosis. EDNRB encodes the endothelin receptor type B implicated in tumor proliferation, survival, invasion, angiogenesis, and metastasis (25). Freije et al. (10) have reported EDNRB as a member of the neurogenesis-related genes group that portends the longest survival. The 3 other genes of our model (CHAF1B, PDLIM4, HJURP) are correlated with a higher risk of death. CHAF1B encodes the p60 subunit of the chromatin assembly factor I, which plays a major role in chromatin assembly after replication and DNA repair. It has been proposed as a specific marker of actively proliferating cells (26) and as a predictor of poor outcome in squamous cell carcinoma of the tongue (27). PDLIM4, a LIM domain gene also known as RIL, is suspected to have tumor suppressor functions in myeloid diseases (28) and prostate cancer (29) by loss of heterozygosity (LOH), deletion, or hypermethylation. However, its extreme upregulation by integrin-promoted demethylation has been recently reported (30) in breast carcinoma cells together with other genes also validated in our study (S100A4, NCAPG), suggesting a potential oncogenic function of PDLIM4. The Holliday Junction Recognition Protein (HJURP) was recently shown to be an indispensable factor for cell-cycle regulation of centromeric chromatin assembly (31, 32) and for chromosomal stability in immortalized cancer cells (33). It has also recently been suggested that HJURP could be implicated in glioma malignancy (34). These studies and our findings suggest that these 4 genes are important molecular components of astrocytic tumors aggressiveness.

The 2 risk groups defined by the 4-gene classifier are also characterized by the expression change of genes related to cancer malignancy or survival of gliomas (Supplementary Table S4). Genes highly expressed in high-risk HGGs are remarkably related to cell cycle and cytokinesis, in accordance with the fact that aggressive tumors exhibit a high percentage of cycling cells. This was also reported for the Proliferative subgroup of HGGs identified by Philips et al. (8). Most of the genes highly expressed in low-risk HGGs are related to the development of the nervous system. Other authors (8, 10, 35) also described a correlation between neuronal markers and the favorable subclasses of HGGs. These findings underline that the 2 risk groups have distinct molecular phenotypes and suggest that they may respond differently to therapeutic regimens.

Multivariate analysis confirmed that both the mutations of IDH1 and the presence of MGMT promoter methylation were associated with a survival benefit in the whole cohort of HGGs and in the subgroup of patients with glioblastoma treated similarly with temozolomide chemoradiation. This analysis also showed that the 4-gene expression risk score was strongly associated with outcome, independently from clinical and molecular risk factors. The performance evaluation indicated that the 4 genes added beyond the prognostic capabilities of all these factors. These results suggest that the 4-gene status, along with the existing clinical and other molecular markers, could be used to optimize patient stratification. As an illustration, when combined with the IDH1 mutational status, the 4-gene risk score allowed the identification of 3 groups of HGGs (good-, intermediate-, and poor-outcome groups) with significant differences in OS (P < 0.001; Fig. 4). The group of HGGs with intermediate-outcome (nonmutated/low-risk or mutated/high-risk) was characterized by a median survival of 20.6 months (95% CI, 16.5–72.1), as compared with 14 months (95% CI, 12.3–15.2) for the poor-outcome group (nonmutated/high-risk) and with a median survival not reached (95% CI, 83.2 to not reached) for the good-outcome group (mutated/low-risk). For this intermediate-outcome group (representing 24% of the whole cohort), the MGMT methylation status did not provide any predictive information (P = 0.5) and the median survival time was similar to that of patients with methylated MGMT promoter. These results suggest the importance of taking into account the 4-gene signature as a stratification factor for the design of future comparative therapeutic trials, though it needs to be further investigated in a prospective clinical study.

No potential conflicts of interest were disclosed.

This work was supported by grants from the Fonds européen de développement régional (FEDER), the Cancéropôle Grand-Ouest Glioma Project, Institut national du cancer (INCa), Ligue contre le Cancer foundation (LNCC), Région Bretagne (PRIR), and CRITT Santé Bretagne.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.