Establishment of prognostic models for astrocytic and oligodendroglial brain tumors with standardized quantification of marker gene expression and clinical variables.

BACKGROUND: Prognosis models established using multiple molecular markers in cancer along with clinical variables should enable prediction of natural disease progression and residual risk faced by patients. In this study, multivariate Cox proportional hazards analyses were done based on overall survival (OS) of 100 glioblastoma multiformes (GBMs, 92 events), 49 anaplastic astrocytomas (AAs, 33 events), 45 gliomas with oligodendroglial features, including anaplastic oligodendroglioma (AO, 13 events) and oligodendraglioma (O, 9 events). The modeling included two clinical variables (patient age and recurrence at the time of sample collection) and the expression variables of 13 genes selected based on their proven biological and/or prognosis functions in gliomas (ABCG2, BMI1, MELK, MSI1, PROM1, CDK4, EGFR, MMP2, VEGFA, PAX6, PTEN, RPS9, and IGFBP2). Gene expression data was a log-transformed ratio of marker and reference (ACTB) mRNA levels quantified using absolute real-time qRT-PCR.
RESULTS: Age is positively associated with overall grade (4 for GBM, 3 for AA, 2_1 for AO_O), but lacks significant prognostic value in each grade. Recurrence is an unfavorable prognostic factor for AA, but lacks significant prognostic values for GBM and AO_O. Univariate models revealed opposing prognostic effects of ABCG2, MELK, BMI1, PROM1, IGFBP2, PAX6, RPS9, and MSI1 expressions for astrocytic (GBM and AA) and oligodendroglial tumors (AO_O). Multivariate models revealed independent prognostic values for the expressions of MSI1 (unfavorable) in GBM, CDK4 (unfavorable) and MMP2 (favorable) in AA, while IGFBP2 and MELK (unfavorable) in AO_O. With all 13 genes and 2 clinical variables, the model R(2) was 14.2% (P = 0.358) for GBM, 45.2% (P = 0.029) for AA, and 62.2% (P = 0.008) for AO_O.
CONCLUSION: The study signifies the challenge in establishing a significant prognosis model for GBM. Our success in establishing prognosis models for AA and AO_O was largely based on identification of a set of genes with independent prognostic values and application of standardized gene expression quantification to allow formation of a large cohort in analysis.

Background

Glioma is the major part of primary malignant brain and central nervous system tumors characterized by tumor cell components of astrocytic glial, oligodendroglial or mixtures of both features. The most malignant and common form of brain tumor is glioblastoma multiforme (GBM, WHO grade IV), comprising about 50% of glioma. The one-year survival rate for GBM patients is 34% and two-year survival is only 12%.1 Histopathological cellular morphology and tumor cytoarchitecture-based grading systems are used to classify gliomas, though providing useful prognostic diagnostic insights, these classification systems do not account for a significant proportion of variation in overall survival among individual glioma patients with the same histology, especially for non-GBM gliomas. Thus far, treatment options are limited for gliomas, for non-GBM gliomas in particular, to radiation and/or cytotoxic chemotherapy without stratification or triage. Prognosis models established using multiple molecular markers in cancer along with clinical variables should enable prediction of natural disease progression and residual risk faced by patients. For patients with glioma currently limited options exist in treatment with poor treatment efficacy. If we successfully achieve development and validation of prognostic models, we will have tools to enable patients to a starting point on the road to personalized treatment.

The power of a prognosis study relies on the number of events, and a meaningful analysis requires an average of 10 patient outcomes per variable used,2 and clinic utility of a prognosis model requires standardization of data. Recent microarray-based gene expression profiling has provided molecular sub-classifications of GBM and identification of genes with prognosis values.3,4 However, direct application of a large panel of gene signature data generated by microarray based studies to modeling prognosis is challenged by the two main issues cited above. In an endeavor of applying gene expression information from tumors to explain patient survival variation, we have taken an alternate approach in modeling glioma prognosis. Our strategy is to establish a prognosis model based on a small to medium size of gene expression variables and a large size of samples and events.5,6 The genes selected to be included in initial model establishment are based on their defined functions in cancer initiation and progression core pathways, and each prognosis value separately reported for glioma. In this study, we applied a standardized platform for real-time quantitative reverse transcription (qRT)-PCR to ensure data comparability.7

This study is a continuing exploration of modeling prognosis for gliomas with a select set of standardized gene expression data in a larger number of cDNA/tumor samples and patient’s follow-up information together with two clinical variables (age and recurrence at the time to tumor sample collection) from the University of Texas, M.D. Anderson Cancer Center (MDACC), the University of Arkansas for Medical Science (UAMS), and the University of California Irvine (UCI). In addition to the prior studied cancer pathway related genes (CDK4, EGFR, MMP2, VEGFA, PAX6, PTEN, RPS9, and IGFBP2),6 we included five cancer stem cell associated genes (ABCG2, BMI1, MELK, MSI1, PROM1) to explore their independent prognosis values in multivariate models for GBM (100 samples, 92 events), WHO grade III anaplastic astrocytoma (AA, 49 samples, 33 events), and oligodendroglial brain tumor, including a mix of WHO grade II oligodendroglioma (O, 18 samples, 9 events) and WHO grade III anaplastic oligodendroglioma (AO, 27 samples, 13 events).

The rationale for including the five stem cell associated genes in glioma prognosis modeling is on the emerging evidences on the existence of stem-like cells in brain tumors responsible for tumor resistance and recurrence.8–11 The most well studied glioma stem-cell associated antigen is CD13312 through expression of PROM1 gene. The expression of other neural stem cell associated genes in glioma have also been reported to have adverse prognostic effects for patients, including BMI1 in oligodendroglial tumors13 and MELK in GBM of younger age patients.14 Functionally involved in self-renewal of neural stem cells,15 BMI1 expression in glioma has been reported to determine tumor phenotype16,17 and control chemo-response via activation of NF-kappaB signaling.18 MELK expression has been shown to regulate the transition from GFAP-expressing progenitors to rapid amplifying progenitors in the postnatal brain19 and promote glioma cell proliferation.14

The expression pattern and prognostic effects of two early identified stem cell associated genes, ABCG2 and MSI1, have not been extensively studied in gliomas. ABCG2 is a marker of adult stem cells,20 and MSI1 marker of CNS stem cells and/or neural progenitor cells.21,22 The quantities of (neural) stem cell associated gene expressions may represent the percentage of NSLC within the brain neoplasm and hence have prognostic values. In this study, we included all of the above five gene expression variables (ABCG2, BMI1, MELK, MSI1, and PROM1) in a multivariate model with two clinical (age and recurrence) and/or eight previously studied neoplastic pathway related genes (CDK4, EGFR, MMP2, VEGFA, PAX6, PTEN, RPS9, and IGFBP2)6 to determine the overall contribution of each variable to glioma prognosis.

Methods

Patients and clinical data

Following informed consent, brain tumor specimens were collected from patients operated in M.D. Anderson Cancer Center (MDACC) at the University of Texas, the University of Arkansas for Medical Sciences (UAMS), and the University of California, Irvine (UCI) and included in this study. MDACC patients with AA, AO, and O were operated during 1987–1997. Majority of GBM patients (n = 59) of MDACC were operated during 1990–1997, while 7 in 2003. UAMS and UCI patients with GBM or AA, or AO were operated during 2003–2006. Following IRB approval, the clinical data for patient’s age, recurrence (at time of sample collection), last contact date, and survival status were provided by each institutes tumor registry and verified by each site investigator in this study. The OS data was calculated based on the time of sample collection and time of death or last contact date.

Study design

This study was designed to establish three main prognosis models with standardized gene expression variables for three distinct gliomas based on histology: GBM model with 100 cases, AA model with 49 cases, and a model for glioma with oligodendroglial components: 27 anaplastic oligodendrogliomas (AO) and 12 oligodendrogliomas (O), as shown in Table 1. We combined AO and O into the same model (AO_O model) to make a total of 45 cases in order to have a meaningful analysis for this type of glioma. The patient survival data was mature for GBM (92% dead event, not including those lost contact for over 5 years), almost mature for AA with death event of 67%, and a 49% death event for AO_O group. Genes selected for inclusion in this study were based on published pilot prognosis data with rationale for their involvement in glioma malignance and resistance as detailed in the Introduction section. We included two main clinical variables that correlated with survival: age and recurrence at the time of surgical removal of the tumor that were also included in our earlier prognostic study for GBM and AA.6 The 1p/19q deletion data for oligodendralgial tumors is not included in modeling due to lack of such information for large number of the studied subjects.

Table 1

Summary of 194 glioma samples and patient follow-up.

Grade1	Histology	CBTRUS 1973–2001	Tumor source			Patients		Survival time (years)			Age
			MDA	UAMS	UCI	Observed	Events2	Median	0.95 LCL	0.95 UCL	Median	Max
1	O	10 yrs at 53.2%	18	0	0	18	9	9.67	4.88	NA	34	59
2	AO	4 yrs at 48.0%	17	3	7	27	13	7.12	2.73	NA	45	60
3	AA	2 yrs at 44.2%	42	5	2	49	33	3.08	1.67	NA	38	53
4	GBM	1 yr at 29.1%	66	21	13	100	92	0.81	0.63	1.00	53	83

Notes: Histology treated as a four-level numeric variable;

Death events were updated to April 2010.

Abbreviations: O, oligodendroglioma; AO, anaplastic oligodendroglioma; AA, anaplastic astrocytoma; GBM, glioblastoma multiforme; LCL, lower 95% confidence limit; UCL, upper 95% confidence limit.

Tumor specimens and tumor cDNA samples

The cDNA samples from MDACC have been used by us and others in published studies, with indication of sample quality control, RNA, and details on cDNA sample processing.5,6,23 The RNA samples of gliomas from UAMS and UCI were processed from 2–4 mm3 snap-frozen tumor pieces using Ultraspec (Biotecx Laboratories, Houston) with an initial homogenization by passing the tissue through a 20-Gauge needle attached to a 1 ml syringe, then RNA extraction following manufacturer’s protocol. The integrity of RNA samples were examined by a RNA gel for the presence of 18S and 28S RNA bands. The cDNA was reverse transcribed from an aliquot of RNA (0.5–2 μg) in a 10 μl reaction using 100 Units of Supercript II reverse transcriptase (Invitrogen, Carlsbad), 10 pmol of poly(dT) 20VN primer, and other components provided in the kit following the manufacturer’s protocol. The cDNA synthesis reaction was diluted 30 times with 10 mM Tris.HCl (pH 7.5), and an aliquot of 4 μl diluted cDNA was quantified for ACTB. Based on ACTB quantity, we further diluted the cDNA, eg, 1000 copy number per 4 μl per quantification reaction, for efficient use of the tumor cDNA samples.

Gene expression data for prognosis study

We used Ziren® Human Real-Time AqRT-PCR Standard-1001 (Ziren Research LLC, Irvine) to quantify ABCG2, BMI1, MELK, MSI1, PROM1, and ACTB, and AqRT-PCR Standard-1020 to quantify PAX6, PTEN, VEGFA, and ACTB for the entire set of cDNA included in this study. IGFBP2 and ACTB were quantified with AqRT-PCR Standard-1009 for non-GBM samples. Quantification of CDK4, EGFR, and MMP2 was carried out in UAMS and UCI on two sets of cDNA samples based on the same single standard containing these three genes, but not the reference gene ACTB. For these three genes, the data was a relative ratio to ACTB. In order to combine the two sets of relative quantitative data in to a single prognosis model, the data in the 2nd set of samples (mainly composed of GBM) was adjusted by the mean of fold difference between the 2nd and the 1st quantification in a set of the same cDNA samples from the 1st set.

We used FAST-START DNA Master SYBR Green I mix (Roche, Indianapolis) in real-time PCR using a LightCycler 2.0 real-time instrument (Roche) or step one real-time PCR instrument (Applied Biosystems, Foster City). The primers are designed to amplify all transcription variants for genes, including EGFR, MMP2, PAX6, and VEGFA. Primer sequence information and PCR parameters are available upon request to Ziren Research (www.zirenresearch.com). The quantification for each gene was repeated 2–4 times for each cDNA sample, and the mean value was used for calculating the ratio of marker gene to ACTB. In our prior glioma prognosis study,5 we have shown that ACTB (called β-Actin there) is a fair reference gene to normalize the marker gene expression among different glioma samples, and we included some of these prior quantified data into this study dataset.

Statistical analyses

From our earlier two studies, we found that modeling glioma prognosis by dichotomizing patients based on recursive partitioning of raw gene ratios produced a model with a higher and biased R2 value than treating the gene expression data as continuous variables after log transformation.5,6 Thus, in this study, we applied multivariate Cox Proportional Hazard (PH) models to estimate the prognoses of GBM, AA, and AO_O using log-scaled ratios (logRatios) of marker vs. ACTB expression quantities. The application of log transformation also avoids outliers on the right and ensures a reliable result. To avoid taking log of zero we used an offset selected to be large enough (such as 0.01) to avoid outliers on the left.

We computed the univariate model coefficients for three glioma groups graded based on tumor malignancy (GBM as grade 4, AA as grade 3, AO_O as grade 1/2) from computer output for the models with grade-gene interactions. The model is represented as follows: log HR = b1*gene + b2*g1 + b3*g2 + b4*g1*gene + b5*g2*gene. b1, b2, b3, b4, and b5 are model coefficients; gene is the log(gene + 0.01) value; g1 is a dummy variable coded 1 for grade 3, 0 otherwise; and g2 is a dummy variable coded 1 for grade 4, 0 otherwise. We also assessed whether the log hazard ratios change with grade by a test for significance of interaction of genes with grade. Kaplan-Meier survival analysis was performed for a highly significant prognosis factor MSI1 dichotomized at the overall median expression in a total of 194 glioma samples including GBM, AA, and AO_O. We used Spearman Rank Correlation Test to analyze gene expression correlations, and Wilcoxon Rank Sum Test to analyze the difference on gene expressions among the three tumor grades (GBM, AA, and AO_O). All statistical analyses were performed using S-PLUS 2000 computer software (MathSoft Inc, Seattle, WA).

Results

The prognostic effect of age and recurrence in gliomas

As shown in Table 1, survival time for patients with GBM is much shorter and less varied than those with other types of gliomas. We performed Cox PH analyses for GBM, AA, AO_O to assess the prognostic value of age and recurrence. As shown in Model 1 in Tables 2 and 3, recurrence, not age, is a significant predictor of poor survival for AA, and none of the clinical variables are significant in GBM and AO_O.

Table 2

Cox PH glioma models based on clinical and gene expression variables.a

Histology	Model 1			Model 2			Model 3			Model 4			Model 5
	GBM	AA	AO_O	GBM	AA	AO_O	GBM	AA	AO_O	GBM	AA	AO_O	GBM	AA	AO_O
R squareb	0.5%	13.6%	5.3%	2.4%	32.7%	57.6%	11.2%	22.0%	35.6%	14.2%	45.2%	62.6%	11.1%	33.6%	58.1%
P valuec	0.774	0.0277	0.299	0.983	0.059	0.0016	0.106	0.095	0.0071	0.358	0.0285	0.0076	0.468	0.142	0.0054
Cases	100	49	45	100	45	34	100	49	45	100	45	34	100	45	34

Notes: An offset of 0.01 was chosen for log-transformed mRNA ratios to avoid outliers on the left;

as R2, an index of the Cox PH model showing the percentage of variation in survival explained by the model;

Likelihood ratio test compares each model to a null model (one with no covariates) to test whether all of the model coefficients are simultaneously equal to zero. The cutpoints used for significance is P < 0.05. Model 1: Two clinical variables (patient age, recurrence status); Model 2: Two clinical variables plus previous studied 7 (GBM) or 8 (non-GBM) genes (CDK4, EGFR, MMP2, VEGF, PAX6, PTEN, RPS9, and IGFBP2 for non-GBM models); Model 3: Two clinical variables plus 5 stem cell associated genes (ABCG2, BMI1, MELK, MSI1, PROM1); Model 4: Two clinical variables plus genes included in Models 2 and 3; Model 5: omitting clinical variables from Model 4.

Table 3

Estimated parameter values, their estimated standard errors and the P-values in a multivariate Cox models shown in Table 2.

To substantiate above finding, we combined all glioma grades in a Cox PH analysis with treatment of histology as a four-level numeric variable, as shown in Table 1. The result revealed that GBM histology has a significant power in prediction of poor survival versus non-GBM with Log (HR) = 1.23 (P < 0.0001), adjusted for other clinical variables (age and recurrence). There is a significant correlation between grade and age (R = 0.42, P < 0.00001), which is consistent with a bias in age with older population in GBM compared to non-GBM (see Table 1). The hazard ratio (HR) for a 20-year increase in age was 1.25 with P = 0.092 in all gliomas (1.41, 2.09, and 1.13 for AO_O, AA, and GBM, respectively), while HR for age did not vary significantly with grade (P = 0.27).

As shown in Tables 2 and 3, tumor recurrent status (recurrence) at time of sample collection, which is a binary variable versus non-recurrence, is an unfavorable prognostic factor in multivariable models for AA, but not for other glioma grades. Consistent with the finding, the HR for recurrence varied significantly with grade (P = 0.025): 1.83, 3.26, and 0.92 for AO_O, AA, and GBM, respectively.

Further analysis was carried out to exam the prognostic effect of grade-age and grade-recurrence interactions on patient’s OS. The data showed that there is a significant GBM–recurrence interaction with Log (HR) = − 1.2606, P = 0.0024, indicating that recurrence has a different function in predicting survival for GBM and non-GBM, which can be explained by the fact that GBM is the highest malignancy and reoperation is beneficial to OS, but recurrence in non-GBM is related to tumor progression into a higher malignancy thus a poor prognostic factor.

Prognostic effect of gene expression variable in univariate models of gliomas

We performed univariate Cox PH assay for each gene separately for GBM, AA, and AO_O. With consideration of sample size and based on their common histological features, we combined AO and O cases in this study, We plotted the hazard ratio (HR) vs. gene expression logRatio curves based on the univariate Cox PH model. As shown in Figure 1 and Table 4, from the 13 genes, only PTEN showed the same decreasing curve and consistently an unfavorable value for Log (HR) in all three glioma types. The other 12 genes showed different effects on prognosis in different glioma subsets classified based on histology.

Figure 1

Hazard ratio vs. gene expression logRatios curves for GBM, AA, and AO_O based on the univariate Cox PH model. The hazard ratio for a particular marker corresponding to each of the three grades was computed using a Cox PH model with 3 terms (grade, marker, and grade-marker interaction), for detail statistical analyses see Method. The hazard ratios are shown using zero as the comparator value. A decreasing curve indicates a favorable prognostic effect from the gene expression. In contrast, an increasing curve indicates an unfavorable prognostic effect, while a flat curve signifies no prognostic effect from the gene expression.

Table 4

Log-hazard ratios computed from the univariate model coefficients for each glioma group (AO_O as grade 1/2, AA as grade 3 and GBM as grade 4).

GENE	Log-hazard ratios			P-value
	AO_O	AA	GBM
PROM1	−0.35	1.14	−0.13	0.44
ABCG2	0.27	−0.20	−0.20	0.053
MELK	0.60	0.40	−0.17	0.097
BMI1	−0.49	−0.18	0.66	0.74
MSI1	−2.67	0.21	1.15	0.0011
PAX6	0.19	−0.15	−0.11	0.44
PTEN	−0.11	−0.30	−0.10	0.91
VEGFA	0.06	0.59	−0.06	0.19
CDK4	0.75	0.37	0.00	0.12
EGFR	−1.14	0.05	0.01	0.091
MMP2	−0.53	−0.59	0.03	0.24
IGFBP2	1.58	0.23	−0.19	0.0046
RPS9	0.35	−0.19	−0.10	0.14

Note: P values are from a test for significance of interaction of genes with grade.

In GBM, the HR curves for five pathway associated genes (CDK4, VEGFA, EGFR, and MMP2) are in general not altered by gene expression levels, with the univariant coefficient Log (HR) being around zero. In contrast, all these genes have either favorable (MMP2, EGFR) or unfavorable (CDK4, VEGFA) prognostic values in univariate models of AO_O and AA (except EGFR). Although over-expressed in GBM, IGFBP2 expression showed a favorable effect on prognosis for GBM. In contrast, IGFBP2 expression in non-GBM gliomas has an unfavorable prognostic value. In accordance with these findings, there is significant interaction of IGFBP2 with glioma grades (Table 4).

The HR-gene expression curve and Log (HR) values revealed consistently that BMI1 expression had an unfavorable prognostic effect for GBM, but favorable for AA and AO_O, and MELK expression a favorable for GBM, but unfavorable for AA and AO_O. PROM1 expression was shown to be unfavorable prognostic factor for AA, but favorable for AO_O and GBM. In contrast, the expressions of ABCG2, PAX6 and RPS9 are favorable prognostic factors for both AA and GBM, but unfavorable for AO_O.

Based on Kruskal-Wallis Rank Sum test, MSI1 is one of the two genes (the other being CDK4) that is not differentially expressed among the three types of gliomas with P > 0.05. However, MSI1 expression shows a strong opposing effect on prognosis for GBM and AO_O that is out of displaying range of HR-gene expression plot, with a significant interaction to grades shown in Table 4. Further analysis of MSI1 prognosis function was carried out in univariate models as continuous and dichotomized variables shown below.

Establishment of multivariate models for GBM, AA, and AO_O

The overall analyses revealed different prognosis effects of the same set of genes in different glioma histopathology classifications, stressing the need of establishing multivariate prognosis models based on histology-classified glioma groups. Using a single data set with logRatios of 13 gene expressions to ACTB, two clinical variables (numeric data for patient age, and binary data for recurrence), and the patients OS time, we performed Cox PH regression analyses for GBM, AA and AO_O. We analyzed different combinations of variables by generating sub-models with two clinical variables (Model 1), with addition of the 8 pathway related genes from our previous study (Model 2), or the 5 stem cell associated genes (Model 3), and addition of all 13 genes (Model 4). We also examined a model only with the 13 gene expression variables (Model 5) to assess the prognostic significance of the genes independent of patient’s age and recurrence of the tumor. Table 2 summarizes the R2 and P value from a likelihood ratio test for the models. Table 3 shows each variable’s log hazard ratio and the statistical significance of each model. The individual models were generated to compare the effect of adding the 5 stem cell associated markers to the original model with 8 cancer pathway associated genes, in order to gain more biological insights on the function of cancer stem cell associated gene expression to residual risk of glioma. We report below a summary of the results from each model for each glioma type.

GBM model

In a multivariate model including 100 GBM patients with 92 events on OS, the two clinical variables (age and recurrence) did not show significant prognostic value and failed to produce a significant prognosis model (R2 = 0.5%, P = 0.77) (GBM model 1). We have shown in our prior glioma prognosis study6 that IGFBP2 expression is significantly correlated with GBM histology and four of the 8 pathway associated genes (MMP2, VEGFA, RPS9, and PAX6) and lacks a significant prognostic value in a multivariate model with these variables, thus it is not included in GBM prognosis modeling in this study in an attempt to increase the ratio for events to variables number. We analyzed prognostic model for GBM by including expression variables of CDK4, EGFR, MMP2, VEGFA, PAX6, PTEN, and RPS9, which made no significant improvement with a model R2 = 2.4% (P = 0.983) (Table 2) and none of the variables showed a significant prognostic value (Table 3). The three genes PTEN, RPS9, CDK4 that showed significant prognostic value in a GBM_AA mixed model with 41 GBM and 43 AA cases in our earlier study6 failed to show prognostic significance in the model with 100 GBM cases.

In the GBM model 3 with 5 stem cell associated gene expressions, a marginally significant improvement was seen in the prognosis model (R2 = 11.2%, P = 0.106). The expression of the genes PROM1 (= 0.41) and ABCG2 (R = − 0.31) showed a significant correlation (P < 0.00001) with grade, R = 0.41, − 0.31, respectively. Both genes showed significant prognostic value in univariate models based on all 194 gliomas—PROM1, Log (HR) = 0.188 (P = 0.015), and ABCG2, Log (HR) = − 0.274 (P < 0.0001). However, when adjusted with dichotomous GBM/non-GBM variables, both genes lost their prognostic value for GBM, suggesting that ABCG2 and PROM1 prognosis values are confounded by glioma grade. In contrast, with a lack of correlation to glioma histology, MSI1 expression was found to be a statistically significant negative prognostic factor for GBM. The unfavorable prognostic value of MSI1 to GBM was seen with or without the inclusion of clinical variables and/or the 7 pathway related genes.

AA model

Recurrence was found to be an unfavorable prognostic factor for the AA group with 2/3 mature patient’s survival information, and together with age variable showed prognostic significance in the AA model 1 (R2 = 13.6%, P = 0.0277). The addition of the 8 cancer pathway related genes (CDK4, EGFR, MMP2, VEGFA, PAX6, PTEN, RPS9, and IGFBP2) or the 5 stem cell associated genes, improved the R2 of the AA model (AA Model 2, R2 = 32.7% and AA model 3, R2 = 22.0%) with marginal significance on the likelihood ratio test. The AA Model 4 with two clinical and 13 gene expression variables together achieved likelihood significance with an explanation of 45.2% of the variation in OS of the patients. In the AA Model 5 excluding the two clinical variables and with the remaining 13 gene expression variables the R2 dropped to 33.6% with a lack of significance in the likelihood ratio test (P = 0.142), consistent with the fact that recurrence is a strong unfavorable prognostic factor for AA.

In AA, the expression of CDK4 and MMP2 showed independent significant prognostic values; CDK4, as an unfavorable prognostic factor and MMP2 as a favorable prognostic factor, based on the Log (HR) and P values shown in Table 3. None of the stem cell associated genes showed prognostic significance in the multivariate models of AA. Although the R2 increased from 13.6% in the model with two clinical variables to 22.0% by adding the 5 stem cell associated gene variables, the model was not significant in the likelihood ratio test. There are multiple pairwise correlations among these genes, which explain the lack of individual significance but an overall improvement of the model in explaining the variation in survival.

AO_O Model

In the multivariate model of AO_O with 48%–50% events in each grade, neither of the two clinical variables showed prognostic significance and failed to produce a significant prognosis model (R2 = 5.3%, P = 0.299). However the addition of 8 pathway related genes (AO_O Model 2) improved the model R2 value from 5.3% to 57.6%, and the model was significant based on a likelihood ratio test (P = 0.0016). Including the 5 stem cell associated gene expression variables also greatly improved the model R2 values from 5.3% to 35.6%, with a significant P value 0.007 in likelihood ratio test of the model. There was a further increase of model R2 to 62.6% (AO_O Model 4) with the addition of all 13 gene expressions. The exclusion of clinical variables only decreased the R2 by 4.5%, consistent with a R2 of 5.3% in a model with only the clinical variables (AO_O Model 1). Thus the 13 gene expression variables we included in this study have significant prognostic value for oligodendroglial tumors.

In the AO_O prognosis model, IGFBP2 expression showed an independent significant unfavorable prognosis effect for oligodendraglial tumors, with Log (HR) = 2.1 (P = 0.022) in AO_O model 2, in which the number of cases dropped from 45 to 34 due to lack of cDNA samples for 9 cases. In a total of 45 AO_O patients, MSI1 expression showed a favorable prognostic value with Log (HR) = −3.1 (P = 0.020), while the expression of MELK an unfavorable prognostic value, with Log (HR) = 1.1 (P = 0.036), independent of clinical variables and the other 4 stem cell associated gene expression variables.

Comparison of models with different sets of variables

We compared Model 4 to Model 1 to assess the contribution of the whole set of genes to the model (with just the two clinical variables. The log-ratio P-values for comparing these models are: <0.0001 for AO, 0.0006 for AA, and 0.25 for GBM, indicating that in aggregate the whole 13 genes contribute significantly for AA and AO_O models, but not for GBM. We compared Model 4 to Model 2 to assess the contribution of the 5 stem cell markers over the cancer pathway associated genes. The log-ratio P-values for comparing these models are: 0.64 for AO, 0.093 for AA, and 0.025 for GBM. Thus, in aggregate the 5 stem cell markers contribute significantly only for GBM. The contribution of the 5 stem cell markers improved predictive accuracy for AA and AO_O, but not to a level of statistical significance. So for AA and AO_O, the models with the new (5 stem cell markers) + previously included (the 8 cancer pathway associated markers) genes were significantly better than models with no genes but not significantly better than models with just the previously included genes.

Differential prognostic impact of MSI1 expression in GBM, AA, and AO_O

The results from this study revealed an interesting finding about the prognostic impact of MSI1 expression in glioma; an unfavorable prognostic value for GBM, no effect for AA, and a favorable prognostic value for AO_O patients. Consistent with the results from the multivariate Cox PH analysis, MSI1 expression as a continuous variable has a significant unfavorable prognostic value in the univariate model for GBM [Log (HR) = 1.23, R2 = 5.2%, P = 0.02], significant favorable prognostic value for AO_O [Log (HR) = −2.83, R2 = 19.2%, P = 0.002], and a lack of prognosis value for AA [Log (HR) = 0.21, R2 = 0.4%, P = 0.65]. We further explored if patient’s survival can be distinguished by dichotomizing based on a biologically relevant threshold of MSI1 expression in gliomas. We set the threshold at an overall median of 0.0012 in the raw (absolute ratio of MSI1 to ACTB) data from all 194 gliomas, and analyzed the hazard ratio (HR) and the survival variation by Kaplan-Meier survival curves. As shown in Figure 2, in contrast to the results from treating MSI1 expression as a log-scaled continuous variable, dichotomizing raw MSI1 expression ratios failed to separate GBM patients with a significant difference in survival. In agreement with the data in log-scale as a continuous variable, higher MSI1 expression showed an unfavorable prognostic value with HR = 1.3, 95% CI = (0.9, 2.0).

Figure 2

Opposing effect of MSI1 in prognosis for GBM and oligodendroglia tumors. Upper panel shows Kaplan-Meier survival curves for GBM, AA, and AO_O based on absolute ratio of MSI1 to ACTB dichotomized at the overall median of 0.0012 for all 194 gliomas. Bottom panel shows log scaled MSI1 univariate models for GBM, AA and AO_O.

Dichotomizing raw MSI1 expression ratios was able to separate AO_O patients with a significant difference in survival time (P = 0.0045) and showed a favorable prognostic value with HR = 0.3, 95% CI = (0.1, 0.7). The effect appears to be independent of AO and O histology, as the distribution of these two groups is not skewed with MSI1 expression. In agreement with non-significant positive prognosis value of MSI1 expression as a continuous variable, dichotomized MSI1 variable in AA also showed a non-significant positive effect on prognosis (P = 0.11) with HR = 0.6, 95% CI = (0.3, 1.2).

Different molecular signatures in GBM, AA, and AO_O based on gene expression correlation analyses

Using a Spearman rank correlation test, we analyzed the pair-wise correlation on gene expressions in different glioma types graded based on histology and survival rate: 4 for GBM, 3 for AA, 2 for AO, and 1 for O. As shown in Table 5, among the 5 stem cell associated genes, there is a lack of significant correlation in GBM. A significant positive correlation between BMI1–ABCG2 (R = 0.35, P = 0.013) and an unfavorable correlation between MSI1–MELK (R = −0.36, P = 0.011) was seen in AA. The BMI1–ABCG2 positive correlation is greater in AO_O (R = 0.50, P = 0.0005), but the other genes lack a significant correlation in AO_O.

Table 5

Spearman rank correlation coefficient matrix with R (upper part) and P (lower part) values for GBM, AA, and AO_O. (correlations with P < 0.02 are in bold face, P ≤ 0.001 are in box).

We further analyzed the correlation between the stem cell associated genes and the 8 previously studied cancer pathway associated genes. As shown in Table 5, three stem cell associated genes (ABCG2, BMI1, and MELK) are positively correlated with PAX6 or PTEN in GBM, AA and AO_O, which have been shown to have a decreased expression in GBM compared to AA or surrounding normal tissues5 and play tumor suppression functions in GBM-derived cell lines.24–27 There is a significant positive correlation between the pro-angiogenic gene VEGFA and the expression of different stem cell associated genes in different tumor grades; PROM1–VEGFA in GBM (R = 0.46, P < 0.0001), MSI1–VEGFA in AA (R = 0.36, P = 0.0117); MELK–VEGFA in AO_O (R = 0.39, P = 0.0074). There is also a significant positive correlation between ABCG2–PAX6 and ABCG2–PTEN across all glioma histologies. VEGFA expression is also significantly positively correlated with PAX6 (R = 0.48, P = 0.0008) and PTEN (R = 0.42, P = 0.004) in AO_O, but not seen in AA and at a reduced level in GBM.

In GBM, there is a general lack of significant correlation between stem cell associated genes and those directly in control of signaling pathways related to tumor aggressiveness, such as CDK4, EGFR, MMP2, and IGFBP2, all were shown to have an increased expression in GBM compared to AA or surrounding normal tissues.5,6 In contrast, there is a strong significant positive correlation between pro-proliferation gene CDK4 and several stem cell associated genes (ABCG2, BMI1, and MELK) in AA and AO_O. MSI1 expression has a significant positive correlation with EGFR (R = 0.59, P < 0.0001) and VEGFA (R = 0.36, P = 0.012) in AA. In AO_O, MSI1 is also positively correlated with EGFR (R = 0.36, P = 0.0141), and MMP2 (R = 0.43, P = 0.0034), but not with VEGFA.

In agreement with results of an independent study on gene expressions in AA and GBM,6 there are significant positive correlations for PAX6–PTEN in either combined (R = 0.53, P < 0.0001) or separate glioma grades with R = 0.52 (P < 0.0001) in GBM and R = 0.51 (P < 0.0001) in AA. The negative correlations between PAX6–IGFBP2 (R = −0.35, P = 0.0004) and PTEN–IGFBP2 (R = −0.43, P < 0.0001) occur only in GBM, while a positive correlation between MMP2–IGFBP2 occurs in both GBM and AA, consistent with IGFBP2 up-regulation of MMP2 expression28 and PAX6 and PTEN suppression of malignant behaviors of glioma cells.24,29 Apart from the above indicated gene expression correlations in AA and GBM, positive correlations between PAX6–VEGFA (R = 0.48, P = 0.0008) and PAX6–RPS9 (R = 0.60, P = 0.0002) are seen in AO_O, and specific correlations between PTEN–VEGFA (R = 0.42, P = 0.004) and PTEN–EGFR (R = 0.50, P = 0.0004) occur only in AO_O.

Overall data from gene correlation analyses revealed different molecular signatures in GBM, AA, and AO_O, which support above analysis of the prognostic effects of these gene expressions separately in each of these tumor types.

Discussion

This is a study of prognosis based on multiple gene expression information in the tumor specimens for histology classified glioma patients without discriminating their difference on treatments received, because the outlook for patients with malignant gliomas has improved very little since the first randomized prospective clinical trials for malignant astrocytomas published in 1978.30,31 Only a modest 2.5 month median survival increase has been achieved by adding concomitant temozolamide to radiotherapy after surgery according to a report published in 2005.32 Our endeavor in modeling glioma prognosis follows the principle of generating reliable models by satisfying the statistical criterion for the ratio between the number of events and variables. Standardized gene quantification ensures data comparability,7 thereby allowing combination of different cohort data to from a large training set for prognosis study. We selected candidate genes for improvement of R2 of the prognosis model, based on their individual prognostic value in univariate models on pilot sets of gliomas from published studies, as well as their functional roles in tumor suppression and progression. Our evidence-based selection of genes for modeling prognosis has produced statistically significant prognosis models for AA and oligodendroglial tumors.

Gene expression information in relation to OS of glioblastoma multiforme

In this prognosis study, in order to achieve a large sample size, we combined GBM samples from patients operated between 1990–1997 from one institute (MDACC) and between 2003–2007 from three institutes (MDACC, UAMS, and UCI), regardless of whether patients were subjected to intensive or less intensive treatments. GBM patients, especially the recurrent ones, have been subjected to different concurrent chemo/radiation clinical trials over the last 40 years without much improvement.30–32 Based on result of this and our prior studies on glioma prognosis, the overall multivariate prognosis model for GBM with 13 gene expressions and 2 clinical variables lacks the power to explain a significant portion of variation on OS. This is unlikely due to sample size issue, given the analysis was carried out based on 100 tumors and 92 events.

Although most of the molecular or clinical variables have shown prognostic values in univariate Cox PH models for a mixed glioma set in our prior studies, their prognostic values were lost when adjusted for GBM in this study. It indicates that GBM diagnosis is in itself a strong prognosis factor, thus genes functionally associated with GBM diagnosis hallmarks (high proliferation index, anaplastic, micro-vascular amplification and/or necrosis) lack independent prognosis values for GBM. Although in general we failed to identify molecular markers for prognosis of GBM, we have identified the MSI1 gene expression as an unfavorable prognostic factor for GBM. Its effect in prognosis was shown by treating it as continuous variable, as we failed to find a cut point for MSI1 expression to dichotomize GBM with a significant difference on OS.

Gene expression information in relation to os of anaplastic astrocytomas

The study of prognosis for AA has been challenged by a lower incidence (7.5% tumors of neuroepithelial tissue) and a longer survival (2 yrs at 44.0%) compared to GBM. AA is often combined with GBM to increase sample size similar to two of our previous prognosis studies.5,6 This study provides for the first time a multivariate prognosis model for AA, with main effects from 13 gene expressions (in log scale) and 2 clinical variables, to explain 45.2% of the survival variation with statistical significance. In contrast to GBM and AO_O, recurrence at operation was found to be a significant unfavorable prognostic factor for AA, independent of the 13 gene expression variables included in the multivariable model.

By modeling prognosis in the multivariable model, we identified CDK4 expression to have an independent significant unfavorable prognosis value for AA, which is consistent with its function in promoting cancer cell proliferation. The other gene with significant prognostic value for AA is MMP2, the matrix metallopeptidase gene overexpressed in glioma and functions in promoting glioma cell invasion.33–35 The data from this study show for the first time that MMP2 expression is an independent significant favorable prognostic factor for AA. Consistent with the idea that angiogenesis drives tumor progression, VEGFA expression has an unfavorable prognostic effect in an univariable model for AA, but not in multivariate models, suggesting VEGFA prognostic function is confounded by other prognostic factors in the model.

Gene expression information in relation to os of oligodendraglial tumors

The same issues that challenge AA prognosis apply to modeling the prognosis of oligodendroglial tumors, which comprise about 8.8% of the overall tumors of neuroepithelial tissue with patient outcomes better than astrocytic gliomas of the same WHO grade. We included mainly those patients operated during 1987–1997 to ensure 50% cases with mature survival information. Based on current follow-up information with the 13 genes and 2 clinical variables, we generated a prognosis model that explains 62.6% of survival variation with statistical significance based on a likelihood ratio test (P = 0.0076). The main contributions come from the 8 cancer pathway related genes that markedly improve the model based on 2 clinical variables (R2 = 5.3%, P = 0.299) to a model with R2 = 57.6% (P = 0.0016). IGFBP2 expression has an independent unfavorable prognostic value for AO_O. There are multiple pairwise correlations among the remaining 7 genes in AO_O, which probably explains the lack of individual significance but an overall improvement of the model. Based on results from univariate analysis, MMP2 expression has a favorable effect for prognosis of AO_O, as seen in AA. Different from that in AA or GBM and in contrast to its usual oncogenic role, EGFR has a favorable prognostic effect in AO_O. Although the P values show a lack of significance, the negative log (HR) values for MMP2 and EGFR are consistent with their favorable prognostic effect in AO_O.

Adding the 5 stem cell associated gene expression variables also greatly improved the model R2 of 5.3% with the 2 clinical variable to a R2 of 35.6% (P = 0.0071). This improvement apparently comes from independent prognostic values of MSI1 (favorable) and MELK (unfavorable). The other three stem cell associated genes have prognostic values in univariate models of AO_O, but are not significant in the multivariate model, which may be explained by their expression correlations, such as the one between ABCG2–BMI1 (R = 50%, P < 0.001).

Both MELK and IGFBP2 have been shown to promote cell proliferation14 and invasion28 in glioma and thus their unfavorable prognostic values are related to differential activation of these two pathways in AO_O. Our finding of a favorable prognosis effect of MSI1 expression in oligodendroglial tumors is in contrast to its unfavorable prognosis effect in GBM. In contrast to GBM in which we were unable to dichotomize patients based on MSI1 expression, we were able to set up a threshold, a median level of MSI1 expression in three glioma sets in combination, to dichotomize patients with AO_O to show a significant difference on OS, regardless of tumor grades (AO or O). This data supports our combining of AO and O into a single set in a prognosis study.

Association of prognostic effects of stem cell associated genes with glioma histology

Results from univariate Cox PH analysis in this study revealed interesting opposing effects from expression of the stem cell associated genes to the prognosis of glioma with different histopathology characteristics. The expressions of ABCG2, MELK, and the neural stem/progenitor cell-associated PAX6 showed unfavorable prognostic effects for AO_O, but favorable prognostic effect for AA and GBM. The expressions of other three stem cell associated genes (MSI1, BMI1, and PROM1) showed a favorable prognostic effect for AO_O, but MSI1 and BMI1 are unfavorable factors for GBM. In contrast to results on CD133 immunostaining and microarray expression data showing that PROM1 is an unfavorable prognostic factor for patients with GBM and oligodendroglial tumors36–41 our data from real-time qRT-PCR quantification in this study set revealed PROM1 as a favorable prognostic factor for AO_O, lack of prognostic value for GBM, and unfavorable for AA. This discrepancy needs to be further investigated for difference in relation to the detection methods as well as source samples.

Concerns of sample size for aa and AO_O prognostic models

This study generated statistically significant prognosis models that are able to explain the variations on OS for 45% and 63% of patients with AA and AO_O. However, based on the statistical criterion for prognosis modeling, approximately 10 patient outcomes per variable,2 the model 4 for AA and AO_O with 49 and 34 patients, respectively, needs to be reassessed in a model with proportionate sample size. Although the P values from the likelihood ratio test showed significance for both models, there is a need for validating the AA and AO_O prognosis models using an independent test set with increase of sample size.