A novel approach to the discovery of survival biomarkers in glioblastoma using a joint analysis of DNA methylation and gene expression

Abstract

Glioblastoma multiforme (GBM) is the most aggressive of all brain tumors, with a median survival of less than 1.5 years. Recently, epigenetic alterations were found to play key roles in both glioma genesis and clinical outcome, demonstrating the need to integrate genetic and epigenetic data in predictive models. To enhance current models through discovery of novel predictive biomarkers, we employed a genome-wide, agnostic strategy to specifically capture both methylation-directed changes in gene expression and alternative associations of DNA methylation with disease survival in glioma. Human GBM-associated DNA methylation, gene expression, IDH1 mutation status, and survival data were obtained from The Cancer Genome Atlas. DNA methylation loci and expression probes were paired by gene, and their subsequent association with survival was determined by applying an accelerated failure time model to previously published alternative and expression-based association equations. Significant associations were seen in 27 unique methylation/expression pairs with expression-based, alternative, and combinatorial associations observed (10, 13, and 4 pairs, respectively). The majority of the predictive DNA methylation loci were located within CpG islands, and all but three of the locus pairs were negatively correlated with survival. This finding suggests that for most loci, methylation/expression pairs are inversely related, consistent with methylation-associated gene regulatory action. Our results indicate that changes in DNA methylation are associated with altered survival outcome through both coordinated changes in gene expression and alternative mechanisms. Furthermore, our approach offers an alternative method of biomarker discovery using a priori gene pairing and precise targeting to identify novel sites for locus-specific therapeutic intervention.

Keywords: :

• glioma
• DNA methylation
• gene expression
• biomarker
• mediation analysis

Introduction

Glioblastoma multiforme (GBM) is the most aggressive of all brain tumors and accounts for approximately 70% of all malignant gliomas.¹ Despite current treatments, patients with GBM have a median survival of only 12–15 mo.¹ This disease is thought to result from the outgrowth of clonal populations that harbor a combination of complex somatic gene alterations.¹ Genetic alterations include dysregulation of many angiogenic and proliferative pathways, including amplification of EGFR and overexpression of VEGF.¹ In addition, dysregulation of many members of the PI(3)K /Akt/RAS signaling pathway have also been implicated in the disease.¹ In 2006, Phillips et al. used these genetic alterations, as well as copy number variations (CNV), to distinguish subclasses of GBM with prognostic implications.² These analyses were further supported by several studies that assessed known, prevalent mutations in GBMs (EGFR, PTEN, IDH1, TP53, and NF1), copy number alterations, and expression changes in an integrative approach in order to more precisely define GBM subtypes important for survival prediction. These data and approaches strongly support the hypothesis that GBMs harbor a complex combination of somatic alterations that determine their phenotype.^3,4

Recently, Frattini et al. (2013) used a novel statistical approach to identify drivers of gliomagenesis through integration of somatic mutations and CNV.⁵ They classified three types of GBM: (1) GBM having deletions at sites containing mutations, (2) GBM having amplifications at sites containing mutations, and (3) GBM with recurrent mutations and no alteration in CNV.⁵ They also identified fusion products involving the EGFR-SEPT14 loci. Their integrative analysis further added to the genetic understanding of GBM pathogenesis and marked specific targets for possible therapeutic intervention.⁵

Epigenetics (particularly DNA methylation) plays an important role in gliomagenesis and glioma survival. Gene promoter DNA methylation has long been associated with gene silencing. Research has now identified a role for methylation in selecting alternate transcripts and gene promoters, giving rise to somatic events that can affect disease survival.^6-10 Our group and others have reported an association between isocitrate dehydrogenase 1 and 2 (IDH1/2) mutations and a hypermethylator phenotype in gliomas that is associated with early age of onset and increased patient survival, specifically in lower-grade gliomas and secondary GBM.^6,11 Our data, which made use of The Cancer Genome Atlas (TCGA) population (a population independent of our original data set), also demonstrated an association between TP53 and G-CIMP, a lack of association between EGFR and G-CIMP, and an overall increase in methylation genome-wide.⁶

DNA methylation does not act solely through the mediation of gene expression (the mechanism that we designate as an expression-based association). DNA methylation has also been found to associate with chromosomal instability, the induction of splice variants, alterations in enhancer regions, changes in microRNA binding regions and expression control regions, and mutations. These somatic changes (which we designate as an alternative association) could also greatly influence survival but are much less well studied.^6-10

These reports have highlighted the crosstalk between various types of carcinogenic somatic alterations and the need for a better understanding of the complex nature of somatic gene inactivation patterns involving genetic and epigenetic alterations that impact both the genesis and survival rates of glioma. Although there has been a call for integrative biomarkers that can sharpen predictive tools, most research has focused solely on the integration of genetic alterations (e.g., mutations) and their association with survival.^5,12,13 Here, we have made use of TCGA data sets to test our bioinformatics-based approach for identifying novel biomarkers of phenotypically-important relationships between DNA methylation, gene expression, and survival in GBM.

Results

DNA methylation and gene expression are significantly associated in GBM samples

After removal of all IDH1 mutant samples and replicates to prevent survival bias, the final phase 1 and phase 2 data sets contained n = 73 and n = 168 samples, respectively. Patient demographic data for all 241 GBM samples are presented in Table 1. Expression and methylation loci were paired by gene symbol for all 241 samples, resulting in a total of 66 202 unique methylation and expression pairs, which were used for the following analysis. In order to ensure functionality of methylation loci in the following analysis, an initial screen was conducted to determine the association of methylation and expression within each gene. To identify the methylation loci that regulate gene expression level, a linear model, as specified in Equation 2 (see Materials and Methods), was performed using the combined phase 1 and phase 2 data sets (n = 241). Pairs were designated as significant if they had a q-value < 0.05. Out of all 66 202 corresponding loci for both expression and methylation, 9821 were found to be significantly associated with each other (84.3% negatively correlated, 15.7% positively correlated). Samples were then separated back into the original phase 1 (n = 73) and phase 2 (n = 168) sets for survival analysis.

Table 1. Patient demographic and tumor* characteristics

		Data Sets
Characteristic	Phase 1 (n = 73)	Phase 2 (n = 168)	Pooled (n = 241)
Age, years
Median	56	60	59
Range	18–86	10–86	10–86
Sex, n (%)
Female	31 (42.5)	69 (41.1)	100 (41.5)
Male	42 (57.5)	99 (58.9)	141(58.5)
**Survival (months)
Median	12.58	10.6	11.3
Range	1.37–60.0	1.08–60.0	1.08–60.0

* All tumor data obtained from The Cancer Genome Atlas (TCGA). **Censored at 60 mo (5 y).

DNA methylation and gene expression pairs are significantly associated with patient survival in GBM samples

To determine which DNA methylation and gene expression pairs are not only significantly associated with each other, but also significantly associated with survival, a Cox proportional hazards model was run on phase 1, phase 2, and pooled data sets. We used the Cox model to investigate the effect of gene expression, DNA methylation, and their interaction term on survival, adjusting for age, gender, and study phase (phase 1 vs. phase 2). “Study” was included as a model variable as a precautionary measure due to the inherent difference in how the presence of IDH1 mutation was determined for each of the two data sets. Tumors with a G-CIMP phenotype or IDH mutation were removed from this analysis due to their association with increased survival in GBM patients. Analysis of the phase 1 data set (n = 73) yielded 878 pairs (from the original 9821) that were significantly associated with survival (P < 0.05). Those 878 pairs were re-run using the phase 2 data set (n = 168) using the same model, which revealed 100 pairs with P < 0.05. Finally, we assessed the effects of the 100 pairs on overall survival using the pooled data set (n = 241) (Supplemental Material, Table S1). Pairs that significantly correlated with survival were chosen based on the q-value (BH) of the pooled model (cutoff: q < 0.10). Thirty-six unique methylation/expression pairs from 29 genes were significantly associated with survival. Of these 36 unique pairs, CpG locus cg23134520 was found to contain a SNP (rs6032566) and was removed from further analysis. This yielded 35 unique methylation/expression pairs from 28 different genes, which were used for the final mediation analysis (Table 2).

Table 2. Final 35 DNA methylation/gene expression pairs that are significantly associated with survival

TargetID_Reporter.REF	SYMBOL
cg17942096_A_23_P165180	RFXANK
cg18345635_A_23_P147345	SLC16A3
cg23943801_A_23_P128166	RAB21
cg27626424_A_23_P34449	LOR
cg05743054_A_23_P419947	MLF1
cg18345635_A_23_P158725	SLC16A3
cg18345635_A_23_P147349	SLC16A3
cg11558474_A_23_P94552	TMEM2
cg01781266_NM_018222_2_3793	PARVA
cg05845503_A_24_P141275	GRHPR
cg05845503_A_23_P60225	GRHPR
cg04551925_A_23_P19894	AQP1
cg00973286_A_23_P139715	TNFRSF1A
cg16773028_A_32_P40593	KCNA2
cg03138091_A_24_P388860	OSMR
cg26475085_A_24_P388860	OSMR
cg24812523_A_23_P14346	AKAP6
cg24302095_A_24_P235266	GRB10
cg24302095_A_24_P235268	GRB10
cg22166290_A_24_P402580	BCL11A
cg03764161_A_23_P203330	FAM111A
cg17726022_A_24_P261734	SLC38A1
cg17726022_A_23_P326510	SLC38A1
cg07663789_A_23_P327451	NPR3
cg04006554_A_23_P214244	ENPP5
cg04006554_A_23_P214240	ENPP5
cg04006554_NM_021572_2_2378	ENPP5
cg05788437_A_23_P80826	FYTTD1
cg06038049_A_23_P35029	CPSF3L
cg20089715_A_23_P405754	CACNB1
cg24219058_A_23_P310921	PCDH7
cg20091959_A_23_P210445	L3MBTL
cg18138552_A_23_P67464	PSMD8
cg20161089_A_24_P270460	IFI27
cg18320336_A_24_P406335	STEAP1

Association of methylated loci with survival can be decomposed into (1) those whose action is mediated through expression and (2) those whose association with survival is not mediated in this fashion

We first estimated the association of DNA methylation with survival mediated through its presumptive effect on gene expression (expression-based association) and then assessed the association not directly mediated through gene expression (alternative association). The expression-based and alternative associations of paired loci with survival were estimated for the top 35 unique methylation/expression pairs (chosen from the linear model and Cox proportional hazards model) by using an accelerated failure time (AFT) model (see Supplemental Material, Table S2). This analysis yielded 10 unique methylation/expression pairs in which expression-mediated methylation was associated with survival outcome (or significant expression-based associations) (Fig. 1A), 13 methylation/expression pairs where methylation did not work through expression of the same gene to affect survival (significant alternative association) (Fig. 1B), and 4 methylation/expression pairs where methylation exerted its effect on survival outcome directly and through gene expression modulation (both significant alternative and expression-based associations) (Fig. 1C). Of the 27 significant methylation and expression pairs, 22 DNA methylation loci were located within a CpG island. In general, pairs within the same gene had similar effects on survival (Fig. 1A–C). In addition, all but three of the locus pairs (associated with CACNB1, RFXANK, and RAB21) had negative correlations, suggesting that the majority of the methylation/expression pairs were inversely related (see Supplemental Material, Fig. S1). Exon locations of methylation loci from significant pairs can be seen in the supplementary material (Fig. S2).

Figure 1. Significant expression-based and alternative associations of DNA methylation on gene expression and survival. The 35 unique DNA methylation/gene expression pairs were subjected to an Accelerated Failure Time (AFT) survival model and applied to alternative and expression-based equations (2–5 in methods). This yielded a total of 27 significant methylation/expression pairs, 10 had significant expression-based associations (A), 13 had significant alternative associations (B), and 4 had both significant expression-based and alternative associations (C). Grey lines indicate alternative associations, black lines indicate expression-based associations, gray circles indicate that the methylation locus for that gene pair was found in a CpG island, and black circles indicate that the methylation locus for that gene pair was not found in a CpG island. The y-axis indicates the change in survival time per 5% increase in methylation; therefore, effects that fall above the line are associated with an increase in survival, and effects that fall below the line are associated with a decrease in survival.

Discussion

The association of alterations in DNA methylation and gene expression in GBM with disease survival has been a major focus of recent studies, as it is apparent that outcome is not solely driven by somatic mutation. These previous studies generally identified loci whose methylation was inversely correlated with expression and examined the impact of those loci on patient outcome. Our study uniquely focused on methylation and attempted to classify the effects of methylation on survival into those mediated by expression and those not mediated by expression, thereby expanding the potential biomarker pool.

In 2013, Wang et al. used an integrative Bayesian analysis (iBAG) approach to analyze the association of DNA methylation with changes in gene expression and subsequently evaluated the association of changes in gene expression on GBM survival.¹⁴ This linear approach resulted in the identification of several genes significantly associated with gene expression modulated by methylation. Consistent with these data, several genes that we showed were significantly modulated by DNA methylation, including OSMR, STEAP1, and GRB10, were also reported by Wang et al.¹⁴ However, methylation not only exerts its effects on survival through expression of its associated gene, but can also operate through a variety of other mechanisms, including chromosomal fragility/instability, splicing variants, enhancer regions, and dysregulation of microRNA.^6-10 Etcheverry et al. (2010) investigated the impact of DNA methylation on gene expression and outcome in GBM.¹⁵ Their analysis focused on the relationship between DNA methylation and gene expression and the association of methylation with survival. They identified 421 CpG sites that were significantly inversely correlated between methylation and expression, 291 of which matched what we found to be correlated in our analysis. They also identified 13 genes that appeared to have consistently differential methylation and expression (between GBM and control brain) but were negatively correlated, suggesting that the regulation of these genes may be epigenetically modulated.¹⁵ However, Wang et al. did not consider the joint effect of methylation and expression on outcome. In addition, IDH1 mutant-associated samples were removed from our study to ensure that the final results would not reflect a bias toward the IDH1 hypermethylator phenotype due to its association with increased survival.⁶

Our final model focuses not only on how methylation acts through expression to affect survival but also assesses how methylation can associate with survival directly or as a proxy for alternative mechanisms (Fig. 2). The final 27 significant methylation/expression pairs contain genes associated with invasion, angiogenesis, and metabolism, and many have been previously linked to brain/glioma (Table 3). Of the 20 genes that contained the significant pairs, none appear to be associated with common amplifications or deletions found in GBM.¹⁶ Ten pairs (from seven genes) had a significant expression-based association with survival, suggesting that DNA methylation in these genes affects survival outcome via expression of the associated gene. Interestingly, two genes contained multiple significant methylation/expression pairs. One of these genes, oncostatin M receptor (OSMR), contained two significant pairs, both with the same gene expression probe, but paired with different DNA methylation loci. The DNA methylation loci for these pairs fall in a CpG island within 550 bp of the transcription start site of the OSMR gene, and the pairs had a negative correlation of methylation and expression, suggesting that methylation of these loci could inhibit gene expression. The locus pairs (cg03138091_A_24_P388860 and cg26475085_A_24_P388860) were associated with a significant expression-based association for each CpG. It is known that OSMR β associates with Interleukin 31 Receptor α (IL31RA) to form the Interleukin 31 receptor (IL31) complex which activates the signal transducer and activator of transcription 3 (STAT3).¹⁷ Priester et al. (2013) recently demonstrated that silencing of STAT3 inhibits glioma single cell infiltration and tumor growth, suggesting that STAT3 plays an important role in the invasiveness of gliomas.¹⁸ If OSMR is silenced via DNA methylation of its promoter, this could lead to a decrease in OSMR gene expression and its association with IL31RA, inhibiting the subsequent activation of STAT3. Without activated STAT3, GBM growth and infiltration could be attenuated, potentially causing increased survival. This proposed mechanism supports the expression-based association of OSMR methylation on survival found in the present study.

Figure 2. Model for mediation analysis. First a linear model adjusted for study was used to determine significantly correlated methylation/expression pairs. Next, a Cox proportional hazards model was used to find significant association between survival and expression, methylation, and their interaction term (adjusting for age, gender, and study phase). An Accelerated failure time model (AFT) was used to estimate the association between survival and expression, methylation, and their interaction term (adjusting for age, gender, and study phase), and a mediation analysis was performed to estimate the alternative and expression-based associations on glioma survival.

Table 3. Functions of significant genes and potential mechanisms in glioma

Symbol	Name	Function (GeneCards®)	Potential expression-based role in glioma survival	Potential alternative role in glioma survival	Ref. #
CACNB1	calcium channel, voltage-dependent, beta 1 subunit	Involved in modulating G protein inhibition	It has been proposed that CACNB1 can protect neurons from Ca(2+)-induced cell death by modulating Ca(2+) channels; therefore, methylation-induced inhibition of CACNB1 could lead to loss of their neuroprotective activities (Ruan B et al 2008).		41
IFI27	interferon, alpha-inducible protein 27	Promotes cell death through mediation of IFN-alpha	?
MLF1	myeloid leukemia factor 1	Oncoprotein that may be involved in lineage commitment	MLF1 and MLF1-like protein were found to co-localize and be over expressed in GBM tumors suggesting they play a role in glioma pathogenesis and survival. (Hanissian SH et al 2005). Dysregulation in expression of MLH1 via methylation could lead to differential survival outcomes.		42
OSMR	oncostatin M receptor	Member of the type 1 cytokine receptor family which heterodimerizes with interleukin 31, which as a complex can induce signaling events	Dysregulation of STAT3 activation via epigenetic induced silencing (Chattopadhyay et al 2007; Priester et al 2013).		17-18
RFXANK	regulatory factor X-associated ankyrin-containing protein	Forms a complex with regulatory factor X-associated protein and regulatory factor 5, which can then bind X box motif regions of some major histocompatibility (MHC) class II molecules, leading to activation	Methylation-induced decrease in RFXANK could inhibit MHC class II activation, which is associated with glioma tumor invasion (Zagzag D et al, 2005).		43
SLC16A3	solute carrier family 16, member 3 (monocarboxylic acid transporter 4)	Part of a family of monocarboxylate transporters that catalyze lactic acid and pyruvate transport across plasma membranes	Differential SLC16A3 expression causing dysregulation of glycolytic metabolism via MCTs (Halestrap AP et al 2004 and 2013; Miranda-Gonçalves V et al 2013; Colen CB et al 2011).		44-46, 53
TNFRSF1A	tumor necrosis factor receptor superfamily, member 1A	This receptor can activate NF-kappaB, mediate apoptosis, and function as a regulator of inflammation	Methylation induced changes in gene expression can dysregulate NF-kappaB pathway, which has been previously associated with glioma tumorigenesis and could be a possible therapeutic target of this disease (Atkinson GP et al 2010).		47
AQP1	aquaporin 1 (Colton blood group)	Molecular water channel protein		Methylation-mediated dysregulation of microRNA mir-320a binding region (Papadopoulos MC et al 2013; Bonomini F et al 2010; Wolburg H et al 2012; El Hindy Ner et al 2013; Saadoun S et al 2005; Sepramaniam S et al 2010).	19-24
ENPP5	ectonucleotide pyrophosphatase/ phosphodiesterase 5 (putative)	It may play a role in neuronal cell communication		Possible dysregulation in angiogenic signaling (Smith SJ et al 2012).	48
FYTTD1	forty-two-three domain containing 1	Required for mRNA export from the nucleus to the cytoplasm		?
KCNA2	potassium voltage-gated channel, shaker-related subfamily, member 2	Voltage-gated ion channel that has a multitude of different functions ranging from regulation of neurotransmitter release, heart rate, insulin secretion, and neuronal excitability		Contains an alternatively spliced product in glioma cells which could contribute to the inactivation rate of the k(+) current Akhtar S et al 1999)	49
L3MBTL	Lethal (3) Malignant Brain Tumor-like 1 (Drosophila)	Polycomb group gene which functions to regulate gene activity via chromatin modifications		?
NPR3	natriuretic peptide receptor C/guanylate cyclase C (atrionatriuretic peptide receptor C)	Natriuretic peptide receptor that regulates blood volume/pressure, pulmonary hypertension, cardiac function and some metabolic/growth processes		?
PSMD8	proteasome (prosome, macropain) 26S subunit, non-ATPase, 8	Regulatory subunit of the 26S multicatalytic proteinase complex, which is involved in the ATP-dependent degradation of ubiquitinated proteins		?
RAB21	RAB21, member RAS oncogene family	GTP-binding protein involved in integrin internalization and recycling		Rab21 expression has been found to attenuate Epidermal growth factor (EFG)-mediated mitogen-activated protein kinase (MAPK) by inducing EGF-receptor degradation (Yang X et al 2012).	50
STEAP1	six transmembrane epithelial antigen of the prostate 1	Found to be upregulated in multiple cancer cells lines and may be a potential metalloreductase		?
TMEM2	transmembrane protein 2	Involved in coordination of myocardial and endocardial morphogenesis (Totong R et al 2011, Smith KA et al 2011)		?	51-52
CPSF3L	cleavage and polyadenylation specific factor 3-like	Catalytic subunit of the integrator complex, which mediates the 3-prime end processing of small nuclear RNAs U1 and U3	?	?
GRB10	growth factor receptor-bound protein 10	Growth receptor binding protein that interacts with insulin and insulin-like growth-factor receptors	Methylation induced loss of imprinting (Blagitko N et al 2009; Monk D et al 2009; Yu Y et al 2011 ;Nord H et al 2009).	Methylation changes in splice variants, leading to expression of alternatively functioning isoforms (Blagitko N et al 2009; Monk D et al 2009; Yu Y et al 2011 ;Nord H et al 2009).	25-28
GRHPR	glyoxylate reductase/hydroxypyruvate reductase	Enzyme that plays a role in metabolism and reduces hydroxypyruvate to D-glycerate and glyoxylate to glycolate and oxidizes D-glycerate to hydroxypyruvate	?	?

?, possible mechanisms relating to glioma and signficant expression-based or alternative assocaiation are unknown.

In addition to the 10 pairs with significant expression-based associations, there were also 14 methylation/expression pairs (in 12 genes) with significant alternative associations. In these genes, DNA methylation is associated with survival either directly or through mechanisms other than direct changes in gene expression. For instance, aquaporin 1 (AQP1) contained one methylation/expression pair, which is located in a CpG island within 300 bp of the transcription start site of the AQP1 gene. The pair showed a negative correlation, suggesting that methylation of this locus could inhibit gene expression. The major function of aquaporins (AQPs) is transportation of water across cell membranes, the disruption of which has been shown to disturb the blood-brain barrier and lead to cerebral edema.^19-21 AQP1 and AQP4 are most abundantly expressed in the nervous system; the expression of both has been observed in GBM and found to correlate with malignancy, specifically cytotoxic cerebral edema, angiogenesis, and migration/invasion.^19,22,23 Recently, it has been shown that both AQP1 and AQP4 are direct targets of several microRNAs including microRNA 320a (miR-320a); furthermore, increased miR-320a is associated with a reduction in AQP1/4 expression.²⁴ Therefore, a possible mechanistic explanation for the alternative association we observe involves methylation of the microRNA target region on AQP1, inhibiting the binding of regulatory mircroRNAs such as miR-320a and ultimately allowing transcription of AQP1.

Interestingly, four methylation/expression pairs (three genes) had both significant alternative and expression-based associations. Of interest is the gene growth factor receptor-bound protein 10 (GRB10), which contained two significant pairs, both with the same DNA methylation locus but paired with different gene expression probes. The DNA methylation locus for these pairs falls in a CpG island of the GRB10 gene, and the pairs showed a negative correlation. The locus pairs (cg24302095_A_24_P235266 and cg24302095_A_24_P235268) have significant alternative associations that suggest a decrease in survival may be observed with a 5% increase in methylation; however, the pairs also have significant expression-based associations. GRB10 is an imprinted gene that is differentially expressed from two promoters. In the brain, it is paternally expressed.²⁵ GRB10 interacts with receptor tyrosine kinases and signaling molecules, most commonly insulin receptors and insulin-like growth factor receptors.^25,26 In addition, monoallelic expression appears to be limited to fetal brain, skeletal muscle, and, most recently, placenta.^25,26 Not only is expression of GRB10 tissue-specific, but it is also isoform specific.²⁵ Currently, 13 different splice variants of GRB10 have been identified, with all but one being expressed in the brain.²⁶ Overexpression of some isoforms has been shown to suppress growth.²⁵ Yu et al. (2011) found decreased expression of GRB10 in many human tumor types, including gliomas, compared with corresponding normal tissue.²⁷ These tumor samples demonstrated a negative correlation between GRB10 and PTEN expression. Furthermore, in a murine cell line, stabilization of Grb10 due to mTORC1-mediated phosphorylation resulted in inhibition of PI3K and ERK-MAPK pathways, suggesting a role for Grb10 as a tumor suppressor.²⁷ Conversely, Nord et al. (2009), using a 32K bacterial artificial chromosome array, found human GRB10 to be a putative novel oncogene in glioblastoma.²⁸ Mechanistic differences may be attributed to inherent imprinting differences in GRB10 between mice and humans. Nonetheless, DNA methylation of this CpG locus has the potential to cause alternative splice sites and may be responsible for the different isoforms of GRB10. Therefore, it is plausible that both the alternative and expression-based associations of this gene have a significant effect on survival. Further potential mechanisms for genes containing significant pairs can be found in Table 3.

There were several limitations to our work. First, we relied upon publically available data, which did not have complete IDH1 mutation data or survival data. We used a previously validated approach ^6,37 to control for this, but this remains a limitation. To address the issue of missing survival data, we used an accelerated failure time model to predict the survival time of censored values. In order to ensure functionality of methylation loci in our analysis, an initial screen was conducted, and only methylation and expression pairs that were significantly correlated within the same gene were used. It should be noted that promoter methylation of MGMT, which has been found in approximately 35–45% of GBM,^29,30 was significantly correlated with MGMT gene expression (data not shown), but was not observed in our final list of significant pairs. This may be attributable to the relatively large number of subjects required to detect an association between treatment and methylation at this locus. Furthermore, MGMT expression has been found to be very low (no more than 15 000 molecules per cell)³¹ and low sensitivity expression arrays have difficulty detecting lowly expressed genes. It has been demonstrated that even in cases where MGMT promoter methylation has been associated with survival outcome, it may not be simultaneously associated with MGMT expression.³² Other mechanisms, such as polycomb repression, may suppress MGMT expression without associated DNA methylation.³³ Therefore, it is not surprising that our study did not detect a significant correlation between MGMT methylation and expression. Furthermore, literature has demonstrated several cases where MGMT methylation has not been associated with survival in temozolomide-treated GBM patients.^34,35,36

Additionally, there was one pediatric patient out of the 241 samples (age 10) that was not removed from the study prior to analysis. Finally, our approach focused on methylation that regulates expression of the same gene, as mentioned above, but would miss methylation loci that do not regulate gene expression and are associated with survival through alternative mechanisms. When establishing significant loci with no gene expression associations, difficulties such as distinguishing null findings arising due to severe multiple comparisons from those with true biology will be an issue.

Overall, our findings are consistent with the accepted concept that DNA methylation can associate with survival outcome via alterations in gene expression (e.g., OSMR). Our findings also suggest that methylation can associate with survival outcome through mechanisms other than dysregulation of gene transcription, including disruption of microRNA function, as is suggested in the case of AQP1. Additionally, some methylation/expression pairs have both significant alternative and expression-based associations, suggesting that different tumors are using discrete mechanisms and yielding different survival outcomes, as described for the proposed alternative and expression-based associations of GRB10. Importantly, our data suggest that this approach might be profitably applied to cancers other than GBM. Our method also brings to light pathways for future study into potential mechanisms in the pathogenesis of glioma. Though additional validation is needed, our work supports the concept that DNA methylation can function both through gene expression and alternative mechanisms to modulate survival outcomes among glioblastoma patients.

Materials and Methods

External data sets

Methylation, expression, and mutation data for glioblastoma multiforme (GBM) were downloaded from The Cancer Genome Atlas (TCGA) for two different sample sets. Level 1 HumanMethylation27 (Illumina) DNA methylation data and level 2 AgilentG4502A_07_1 and 2 gene expression data were downloaded for all available GBM batches. GBM batches 1, 2, 3, and 10 were used as the phase 1 set and GBM batches 16, 20, 26, 38, and 62 were used as the phase 2 data set. Patient samples lacking covariate data were removed; samples were further restricted to patients diagnosed with glioblastoma who were alive 30 d after their date of diagnosis. Data sets were not combined in further analyses due to the fact that phase 2 data did not have definitive IDH mutation status. Since IDH mutations are associated with survival, we were hesitant to combine the two data sets as mis-identification of IDH mutations could grossly affect findings.

Recursively partitioned mixture model to determine IDH1 mutation status

Patient survival, DNA methylation, gene expression, and IDH1 mutation data (phase 1 set only), was obtained for primary glioblastoma multiforme (GBM) samples. It has been widely acknowledged that IDH1 mutants are almost exclusively associated with a hypermethylater (G-CIMP) phenotype, and this phenotype is associated with increased survival in glioma.^10,11 Therefore, we wanted to remove IDH mutant samples from our study so results would not be biased due to increased survival associated with this mutation. Since IDH mutation data was not available for the phase 2 sample set, we employed a recursively partitioned mixture model (RPMM) as described by Houseman et al.³⁷ and used in Christensen and Smith et al.⁶ The RPMM successfully divided the phase 1 set into seven classes (see Supplemental Material, Fig. S2), and the samples in the top two most highly methylated classes, along with the samples having IDH mutations in the phase 1 set, were removed (TCGA.14.1458, TCGA.16.1460, TCGA.19.1788, TCGA.14.1456, TCGA.28.1756, TCGA.14.4157, TCGA.32.4208).

Methylation data

Methylation β values were extracted from raw idat files using GenomeStudio software (Illumina), which calculates β values using M/(M + U + 100), where M is the methylated signal, U is the unmethylated signal, and 100 is an arbitrary offset. Replicates that did not correlate were removed (TCGA.06.0137, TCGA.06.0145). For methylation loci, all loci that contained a detection P > 0.05 for any sample were removed from further analysis. Since approximately 25% of the survival data are censored, censored survival times were estimated using an accelerated failure time (AFT) model based on the Equation 1 below.

log(T) = b₀ + b₁Age + b₂Gender + b₃Study + b₄(Age + Study) + b₅(Study + Gender) + με

where T follows a Weibull distribution³⁸ (μ is a scale parameter, and ε follows an extreme value distribution). Next, methylation values were normalized for bead chip to control for potential batch effect using the ComBat method³⁹ with adjustment for age, gender, survival, censored data, and survival-censored interaction.

Expression data

TCGA expression and methylation subject identification numbers were matched; all non-matching samples were removed from the data sets. Replicates in expression samples were either averaged or chosen based on the closest mean and standard deviation to the methylation distribution across all samples. The final data sets consist of a phase 1 data set (n = 73) and a phase 2 data set (n = 168) that contain complete data on overall survival, DNA methylation, and gene expression, with samples considered G-CIMP removed.

Final methylation/expression locus pairs

Methylation and expression loci were merged based on gene of origin. Annotation files for both platforms (HumanMethylation27 and AgilentG4502A_07_1 and 2) were downloaded from TCGA and matched by gene symbol, (using the manufacturer’s annotation) yielding 66 202 methylation/expression pairs. It should be noted that there are usually several methylation loci and/or expression probes found within each gene, so while each pair is unique upon merging, an individual methylation or expression locus may be repeated among several pairs.

Statistical analysis

To choose statistically significant methylation and expression pairs, expression was regressed on methylation in the pooled (n = 241) data set. The associated p-values were adjusted for false discovery rate (FDR) using the Benjamini-Hochberg (BH) procedure.⁴⁰ All methylation/expression pairs that had a q-value < 0.05 were identified as being significantly associated with each other (n = 9821 pairs).

To further siphon out statistically significant pairings, pairs were then assessed using a Cox proportional hazards model for the effect of expression, methylation, and their interaction on survival, controlling for age, gender, and study phase (phase 1 and phase 2, when applicable). A three degree-of-freedom (DF) Chi-square test was performed to test for significance of expression, methylation, and their cross-product interaction. The three-DF models were repeated for both phase 1 (n = 73) and phase 2 data sets (n = 168) separately and the pooled data set (n = 241). In order to reduce false positives, final statistically significant pairs were selected for having P < 0.05 in both phase 1 and phase 2 data sets and q-values of <0.1 in the pooled data set.

The associations of methylation and expression on survival were determined by a mediation analysis adopted from VanderWeele³⁸ using the following equations for the expression-based and alternative associations of methylation on survival:

Equation 2. E[E|M, c] = β₀ + β₁M + β₂c

Equation 3. log(T) = θ₀ + θ₁M + θ₂E + θ₃EM + θ₄c + νε

Equation 4. △_M→E→T = (θ₂β₁ + θ₃β₁m)(m − m*)

Equation 5. △_M→T = (θ₁ + θ₃[β₀ + β₁m* + β₂c + θ₂σ²])(m − m*) + 0.5θ₃²σ²(m² − m*²),

where T is survival time, E is expression, M is methylation, c is study, σ² is the variance of the error term in Equation 2, ε is a random error in Equation 3 following the extreme value distribution, and ν is a scale parameter. For our purposes, m* is median methylation and (m-m*) is the change in methylation we are interested in observing. For example, we would set m-m* to 0.05 if we wanted to look at the change in survival for a 5% increase in methylation. Equation 2 represents the linear model for the association between expression and methylation, and Equation 3 represents the accelerated failure time model with interaction between methylation and expression. β₀-β₂ are the regression parameters for the linear model, and θ₀-θ₄ are the regression parameters for the accelerated failure time model. We used a stepwise mediation analysis that considers the relationships between methylation and expression (Eqn. 2) and their joint effect on survival (Eqn. 3). In our case, an alternative association is the effect that methylation alone (or as a proxy for alternative mechanisms) has on survival, and expression-based association is the effect of methylation on survival mediated through gene expression. Equation 4 represents the expression-based association, and Equation 5 represents the alternative association of methylation on survival,³⁸ both of which can be estimated by fitting the models in Equations 2 and 3. We used bootstrap to find the variances and confidence intervals of the expression-based and alternative associations.

To determine directionality of the association of methylation on expression, we looked at the coefficient in the linear model regressing expression on methylation (Eqn. 2). A negative coefficient suggests that methylation and expression are inversely related (i.e., increased methylation is associated with decreased expression and vise versa). A positive correlation demonstrates that methylation and expression are directly related (i.e., increased methylation is associated with increased expression).

Abbreviations:
GBM	=	glioblastoma multiforme
CNV	=	copy number variation
G-CIMP	=	Glioma CpG Island Methylator Phenotype
AFT	=	accelerated failure time
iBag	=	integrative Bayesian analysis
BH	=	Benjamini-Hochberg
FDR	=	false discovery rate
DF	=	degrees-of-freedom

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

The authors would like to thank all individuals involved in the TCGA, particularly the patients who donated samples for use in this research.

Funding

This work was funded by grants from the National Cancer Institute (CA100679, K.T.K.); and the National Institutes of Health (CA126831, J.K.W).

10.4161/epi.28571