{'fraction_non_alphanumeric': 0.03770833717861164, 'fraction_numerical': 0.022038059026560107, 'mean_word_length': 5.057015092230296, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrative Modeling Reveals Annexin A2-mediated Epigenetic Control of Mesenchymal Glioblastoma

Glioblastomas are characterized by transcriptionally distinct subtypes, but despite possible clinical relevance, their regulation remains poorly understood. The commonly used molecular classification systems for GBM all identify a subtype with high expression of mesenchymal marker transcripts, strongly associated with invasive growth. We used a comprehensive data-driven network modeling technique (augmented sparse inverse covariance selection, aSICS) to define separate genomic, epigenetic, and transcriptional regulators of glioblastoma subtypes. Our model identified Annexin A2 (ANXA2) as a novel methylation-controlled positive regulator of the mesenchymal subtype. Subsequent evaluation in two independent cohorts established ANXA2 expression as a prognostic factor that is dependent on ANXA2 promoter methylation. ANXA2 knockdown in primary glioblastoma stem cell-like cultures suppressed known mesenchymal master regulators, and abrogated cell proliferation and invasion. Our results place ANXA2 at the apex of a regulatory cascade that determines glioblastoma mesenchymal transformation and validate aSICS as a general methodology to uncover regulators of cancer subtypes.


Introduction
Glioblastoma multiforme (GBM) is the most frequent primary brain tumor in adults. Despite intense efforts at defining novel targeted therapies, the prognosis of GBM patients remains poor, with a median survival of approximately 12 months. Molecular profiling studies of GBM surgical samples have revealed distinct subtypes of GBM (Brennan et al., 2009;Phillips et al., 2006;Verhaak et al., 2010). Although the proposed subtyping systems differ in the details, a shared feature of current proposals is a subtype with high expression of mesenchymal marker transcripts. Mesenchymal subtype GBM exhibits decreased expression of the Neurofibromin 1 gene (NF1) and high levels of necrosis-associated genes (Brennan et al., 2009;Phillips et al., 2006;Verhaak et al., 2010).
While there is a strong association between mesenchymal gene expression and tumor invasiveness, in both GBM and other tumors (Balbous et al., 2014;Carro et al., 2010;Kalluri and Weinberg, 2009;Thiery et al., 2009), the cellular networks that drive mesenchymal transformation are far from understood.
Genes that modulate mesenchymal transformation in GBM have previously been found by computational analysis of mRNA profiling data. These studies have revealed the transcription factors CEBPB and STAT3 (Carro et al., 2010) and WWTR1 (a.k.a TAZ) (Bhat et al., 2011) as master regulators of mesenchymal GBMs. Moreover, Bhat and colleagues showed that NF-κB can induce the mesenchymal signature through induction of CEBPB, STAT3 and WWTR1 (Bhat et al., 2013). Amplified and overexpressed RHPN2 also promotes mesenchymal transformation (Danussi et al., 2013) and microRNAs miR-128a and miR-504 have been shown to negatively correlate with mesenchymal gene expression (Ma et al., 2012). Despite the progress in mapping transcriptional networks of GBM phenotypes, our understanding of how these networks are modulated by genetic, epigenetic factors or miRNAs remains incomplete. For such analyses to be possible, new approaches are warranted that take several types of evidence into account.
To conduct such an open and data-driven search for regulators of all GBM subtypes, we here adapt the recently described Augmented Sparse Inverse Covariance Selection (aSICS) method for integrative cancer data analysis (Kling et al., 2015), which has two suitable features for the problem at hand. First, unlike previous tools for detection of subtype regulators (Cantini et al., 2015;Carro et al., 2010), it combines multiple types of data into a single network model, thus allowing us to efficiently screen for associations across several classes of molecular events along with sample data on mRNA and miRNA expression, DNA copy number aberrations (CNA) of genes and miRNAs, DNA methylation, loss of heterozygosity (LOH), point mutations, and clinical information. Secondly, the procedure is based on robust estimation of partial correlations to detect the coupling between two variables after correction for all the other variables. Previously, aSICS has been shown to exhibit robust properties compared to existing network construction methods on cancer data (Kling et al., 2015).
We first demonstrate that, in addition to including multiple layers of data, aSICS outperforms a standard method in terms of consistency between two glioma data sets. Next, we use aSICS to build a multi-layered model of regulation of glioma subtypes. Our model identifies both known and novel regulators of mesenchymal, proneural and classical subtypes. Among the predicted regulators, Annexin A2 (ANXA2) stands out as an epigenetically controlled master regulator of mesenchymal transformation in glioma, associated with patient survival. To validate the functional relevance of this model prediction, we abrogate ANXA2 expression in patient-derived, IDH1 wildtype glioma cancer stem cells, leading to loss of mesenchymal signature genes. Interestingly, ANXA2 loss led to reduced phosphorylation of the previously described key regulator of mesenchymal transformation, STAT3 (Carro et al., 2010). Indepth analysis of two independent patient cohorts further supports that ANXA2 suppression in lower grade glioma is likely explained by methylation induced by IDH1 mutation. However, ANXA2 also retains IDH1-independent prognostic power in IDH1 wildtype higher grade glioma. Analysis of heterogeneous glioma materials, both surgical samples from different brain regions and single cell data, is used to show that ANXA2 correlates regionally with mesenchymal transformation.
Our analysis and results thus support a new functional link between epigenetic regulation and mesenchymal transformation, and identifies ANXA2 as a possible therapeutic target against mesenchymal glioma. Application of aSICS to Breast, Ovarian and Colorectal cancer identified multiple genetic, epigenetic and transcriptional regulators of subtypes in these cancers, illustrating its generality. The aSICS method is made available as a free-of-charge Matlab package.

Applying aSICS to Define Regulators of Cancer Subtypes
The below Methods description focuses on our adaptation of aSICS for the purpose of uncovering subtype regulators. For a full technical description of aSICS (e.g. fitting algorithm, FDR estimation, parameter values, distance metrics, and different data types), we refer to (Kling et al., 2015).
Step 1: Data integration. All datasets except somatic mutations were downloaded as Level 3 (gene level) data from TCGA (The Cancer Genome Atlas, http://cancergenome.nih.gov) database. We chose the platform for each data type to maximize the number of patients in that dataset. RNA and miRNA measurements were used as provided by TCGA without further normalization. DNA copy number aberrations were mapped to mRNA and miRNA transcripts by mapping the genomic coordinates of each transcript to the log2 relative DNA copy number segmentation map provided by TCGA. If more than one segment was mapped, copy number was estimated as a weighted average across the involved segments, calculating the weights based on the relative overlaps. Genes or miRNAs with a CNA value, but lacking an expression measurement were discarded from the analysis. Methylation beta values (defined as the ratio of the methylated probe intensity and the sum of methylated and unmethylated probe intensities) were filtered to keep probes with standard deviation across the patients N0.05. Somatic mutation state was assigned as 0 or 1, the state 1 assigned if a gene had at least one missense mutation called by TCGA (level 2 data) in that sample. Silent mutations were not considered. Genes with at least one missense mutation in fewer than 5 patients were not included in the analysis. Glioblastoma subtype assignments were assigned by GSEA (gene set enrichment analysis) (Mootha et al., 2003;Subramanian et al., 2005) using the signatures provided in (Verhaak et al., 2010). The four subtypes are classical (n = 159), neural (n = 93), mesenchymal (n = 164) and proneural (n = 113). The collected data were subsequently organized into a block correlation matrix: where each block Sab contains the Pearson linear correlation across patients for all pairs of variables in data types a, b ∈ {1,⋯, k}, respectively (i.e. in this notation, we consider 1 = point mutations, 2 = copy number aberrations etc) (the combination of continuous (like expression) and binary (like mutations) is discussed and validated in (Kling et al., 2015)).
Step 2: Network construction. Given the correlations S as input, aSICS uses an iterative algorithm estimate a network of partial correlations, Θ by solving the following optimization problem: The penalty function P (Θ, p, τ) is, in turn, defined as: where α is the elastic net parameter (Zou and Hastie, 2005), selected as in (Kling et al., 2015). λ 1,ij = τ ν ij is a link specific penalty for the pair of network nodes i and j. τ regulates the overall stringency, i.e. the sparsity of the network. ν ij take on three possible values: 1, p (b1), or ∞, and makes it possible to emphasize specific network features. ν ij = ∞ is used on CNA and mRNA interactions from differing locus, based on the assumptions that CNA affect transcription locally, and for within datatype interactions for CNAs and methylations, to avoid that the network model is strongly dominated by links between CNAs or methylation variables of close genetic proximity. ν ij = p is used for interactions 1) between miRNAs with their predicted mRNA targets, as defined by miRanda (John et al., 2004) prediction (MicroCosm Targets Version 5 (Griffiths-Jones et al., 2008, http://www.ebi.ac.uk/enright-srv/microcosm/htdocs/targets/v5/), 2) between cis localized methylations probes with their corresponding mRNA, as defined by associations between genes and methylation probes provided in the TCGA level 3 data, and 3) between all interactions involving a point mutation. The implications of this prior is further discussed in (Kling et al., 2015). To increase the confidence of network links we apply bootstrapping and aggregate the results, as in (Kling et al., 2015). Specifically, we randomly chose 90% of the tumor samples and estimate the network 100 times. The bootstrap results are aggregated into a final network Z = z ij , in which z ij = +1 if the link between variables i and j is present in at least 80% of bootstrap runs and the correlation sij is positive. Similarly, z ij = −1 when the link is found in at least 80% of bootstraps and sij is negative (s ij denotes the pairwise correlation between variables i and j). We have previously shown that high levels of network stringency tend to enrich for functional interactions (Kling et al., 2015).
Here, we used a variant approach by which the network stringency (τ) was selected by maximizing the fold enrichment (FE) of genes from our in-house database of genes that have been functionally validated in mouse GBM models (c.f. Supplementary Table 3). FE is defined as the odds ratio P (in database | selected as regulator)/P (in database|not selected as regulator).
Given the aggregated network Z = {z ij }, the net connectivity, c, between a subtype r and a regulator s, is given by the sum where RNA the set of indices for variables representing mRNAs (as a concrete example, cg08081036 is negatively linked to the mesenchymal subtype in Table 1, since it is the negative regulator of ANXA2 which, in turn is a positive regulator of mesenchymal subtype, i.e. c = (− 1)(+ 1) = − 1). The number of regulators (c ≠ 0) and the magnitude of c is controlled by the stringency parameter τ, discussed above.

Exploration of Additional Cancers
In our exploration of additional cancers (Supplement), we selected three cancers for which there is clear evidence of molecular subtypes: ovarian adenocarcinoma (TCGA 'OV' cohort) Network et al. (2011), a recent reanalysis of colorectal cancer from the (TCGA 'COAD' cohort) (Cantini et al., 2015) and breast cancer (TCGA 'BRCA' cohort) (Network et al., 2012), using the provided subtyping data in each respective publication. aSICS was run as above using the same settings as for GBM.

Validation in Two Independent Cohorts
We used the Phillips et al. GBM cohort of 100 cases (Phillips et al., 2006) and the 529 expression profiles from the TCGA cohort to confirm that similar sets of TF regulators were detected. The analysis was repeated for a range of values of τ to control the number of predictions. Master regulator was applied as described (Carro et al., 2010), using the ARACNE I0 parameter to control the number of predictions. For each method, we recorded the regulators that did not overlap between data sets (Fig. 1d, x axis) and the ones that did (Fig. 1d, y axis). The slope for each method and the null expectation, respectively were assessed and tested by linear regression.

Functional Annotation of Networks
For Supplementary Fig. 1, hierarchical clustering of the network was performed using the adjacency matrix indicating existence or absence of links. The pairwise topological overlap (Ravasz et al., 2002) between nodes was used for the distance matrix and clustering was performed for each connected component of the network separately. The number of clusters for each component was chosen as the minimum number for which adding another cluster did not improve the average silhouette width. Gene set enrichment was computed for each cluster using the Fisher's exact test comparing overlap of node gene assignments with v5.0 of the MSigDB database (Subramanian et al., 2005), including only the c2, c3, c5, c6 and c7 gene set collections. p-Values were adjusted using the Benjamini-Hochberg method and considered significant if less than 0.05.

Cell Lines and Culture Conditions
Glioblastoma cell lines and HEK 293T cells were routinely grown in DMEM with 10% FBS. Primary glioblastoma-derived BTSCs were prepared from tumor specimens collected at the University of Freiburg (BTSC 168 and BTSC 161) and University of Uppsala (U3047MG) according to a study approved by the University ethic committees and upon patient approval. Tumor tissues were collected in the operation room and immediately transferred to a sterile hood for cell dissociation. The tissue was initially washed in sterile PBS and subsequently cut into small pieces with a scalpel. After PBS removal by low-speed centrifugation (700 rpm for 5 min), the tissue was incubated with Cell Dissociation Solution (Sigma) for 5 min. During this time the cells were also mechanically dissociated by pipetting with a 200 l pipette tip. The Cell Dissociation Solution was then removed by centrifugation and the tissue resuspended in 5 ml of ACK Lysis Buffer (Life Technologies) and incubated for 5 min at room temperature. After centrifugation, the cells were resuspended in 5 ml of neurobasal medium and filtered through a 10 m Nylon-Filter (BD Falcon), collected in a new tube, and transferred to a 75 cm 2 cell culture flask (low adherent, Nunclon). The established BTSCs were then cultured as neurospheres in Neurobasal medium (Invitrogen) containing B27 supplement (Invitrogen), FGF (20 ng/ml, R&D Systems), EGF (20 ng/ml, R&D Systems), LIF (10 ng/ml, Genaxxon biosciences), Heparin (2 μg/ml, Sigma), and glutamax (Invitrogen). U3047MG BTSC were grown on Laminin.

Vectors and Viral Infection
Knockdown of ANXA2 was obtained with a shRNA lentiviral vector (pLKO, Sigma; ); ANXA2, originally provided in a 1st generation lentiviral vector (kindly gifted by L. Zheng, Johns Hopkins University, Baltimore, MD), was overexpressed by cloning it into a pCHMWS lentiviral vector (kindly provided by V. Baekelandt, Katholieke Universiteit Leuven) in frame with a FLAG tag. To infect BTSCs, supernatant from transfected 293T cells was ultracentrifuged to produce high-titer, serum-free virus; Polybrene 4 μg/ml (Sigma) was added to the lentiviral particles during the infection. Selection of infected cells was conducted with 2 μg/ml of Puromycin (Sigma).

Quantitative Real-time PCR
RNA was prepared using the RNAeasy kit or the All Prep DNA/RNA Protein Mini Kit (Qiagen) and used for first strand cDNA synthesis using random primers and SuperscriptIII Reverse Transcriptase (Invitrogen). Quantitative real-time PCR (qRT-PCR) was performed using the following pre-validated TaqMan Assays (Applied Biosystems); ANXA2: Hs00743063_s1, 18srRNA: Hs99999901. Quantitative RT-PCR with SYBR green (Applied Biosystems), and primers listed in Supplementary Table 3, were used to validate a panel of mesenchymal genes.

Invasion Assays
The Matrigel invasion assay was done using BioCoat Matrigel Invasion Chambers (BD Biosciences) according to the manufacturer's instructions. 2.5 × 10 4 cells infected with lentivirus were seeded in the upper compartment. PDGF-BB (20 ng/ml, R&D) was used as a chemoattractant. Images were collected using a wide-field microscope (Axiovert, Zeiss).

EdU Cell Proliferation Assay
Cell proliferation was assessed using the EdU-Click594 Cell Proliferation Imaging Kit (Baseclick GmbH) according to the manufacturer's instructions. 2.0 × 10 4 cells were seeded on laminin-coated glass coverslips in a 24-well cell culture plate. Pictures were acquired using an Axiovert Microscope (Zeiss).

Osteogenesis Differentiation Assay
The osteogenesis differentiation assay was performed using the StemPro Osteogenesis Differentiation Kit (Life Technologies) according to the manufacturer's instructions. Briefly, 5 × 103 cells/cm 2 were seeded on laminin-coated glass coverslips in a 24-well cell culture plate. Cells were incubated in MSC Growth Medium at 37C, 5% CO 2 for 21 days, replacing the medium every 4 days. Cells were then fixed with 4% formaldehyde, stained with Alizarin Red S solution (pH 4.2) and mounted on microscope slides. Pictures were acquired using an Axiovert Microscope (Zeiss).

Gene Expression Array and Gene Set Enrichment Analysis (GSEA)
For gene expression profiling total RNA was prepared using the RNeasy kit or the all Prep DNA/RNA/Protein mini kit (Qiagen) and quantified using 2100 Bioanalyzer (Agilent). 1.5 μg of total RNA for each sample was sent to the genomic facility of the German Cancer Research Center (DKFZ) in Heidelberg (Germany) where hybridization and data normalization were performed. Hybridization was carried out on Illumina HumanHT-12v3 expression BeadChip. Microarray data were analyzed using the GSEA software (http://www.broadinstitute.org/ gsea).

Classification of Brain Tumor Samples
The classification of brain tumor samples was performed by using 510 genes out of the 840 classifier genes used by Verhaak et al. to classify 260 glioblastoma samples (Verhaak et al., 2010), and 529 glioblastoma tissue samples from TCGA with assigned subtypes as reference (Network, 2008). The 510 genes were selected such that they classify well the extended set of 529 TCGA samples and were represented on the Illumina HumanHT-12v3 expression BeadChip arrays. The expression levels for these genes on the Illumina arrays and in the TCGA data set were converted into z-scores and the combined matrix was used to classify each glioma sample based on a k-nearest neighbors (k = 10) and voting procedure, in which a subtype was assigned based on the majority subtype among the 10 TCGA samples with highest correlation coefficients for these genes with respect to the tumor sample. All data manipulations were performed in R (https://www.r-project.org/) and MATLAB (The MathWorks, Inc., Natick, MA, United States).

Pyrosequencing
Genomic DNA was extracted from tissues and cell lines using all Prep DNA/RNA/Protein mini kit (Qiagen) and quantified using NanoDrop2000c (Thermo Scientific). For pyrosequencing, genomic DNA was bisulfite modified using the EZ DNA methylation gold kit (The Epigenetics company) according to the manufacturer's instructions. Pyrosequencing analysis was performed using a PyroMark Q96 instrument (Qiagen), following the manufacturer's protocol. The following primers were used: ANXA2-pyr-for GGGTAGGGGTGAGTTATTTTTGATTT, ANXA2-pyr-rev biotin-ACAAAACCCAC-TAACCTAAATAAAACTTTTATACC, ANXA2-pyr-seq GGGTAGGGGTGAGTTATTTTTGATTT. The results were analyzed by PyroMark CpG software (Qiagen). The methylation index for each sample was calculated as the average value of CpG methylation in the CpG examined.

Methylation Array and Analysis
The methylation array was performed using the Illumina Infinium HumanMethylation450 chip according to the manufacturer's instructions. Data analysis was performed by R software and RnBeads software package.

MRI-localized GBM Biopsy Collection
Before surgery, standard gadolinium-enhanced T1-weighted MRI 3D datasets were acquired. According to the MRI data, three tumor regions were identified for sampling: i) a non-enhancing, peritumoral area referred as edema zone; ii) a contrast-enhancing area within the tumor; iii) a non-enhancing intratumoral area, referred to as core. During surgical resection of the tumor, pre-originated MRI coordinates were used to accurately track the sampling positions. All biopsy specimens were collected upon patients approval by an experienced neurosurgeon, subsequently snap-frozen in liquid nitrogen, and stored at −80°C. This study was approved by the University of Freiburg ethic committee.

A Comprehensive Pipeline to Uncover Modulators of Cancer Subtypes
For aSICS analysis of tumor subtypes, we first developed a novel pipeline with three steps. In step 1 (data integration), we combine different genomic data types with binary subtype information, e.g. according to Verhaak et al.'s system (Verhaak et al., 2010) (Fig. 1a). We subsequently compute all pairwise correlations between variables in the joint data set. Such correlations can be informative, but will also contain multiple correlations that are due to indirect or spurious effects in the data. In step 2 (aSICS network construction), we process the correlation matrix to remove such indirect correlations. This is done by an efficient statistical optimization algorithm to select a set of robust pairwise partial correlations (Fig. 1b). Partial correlations with strong evidence in the data are interpreted as network links, which can either be positive (+) or negative (−). To ensure that the links have strong support in data, we aggregate the results from 100 bootstrap runs, keeping only links with consistent evidence support in at least 80% of the runs. The procedure is adjusted by only two parameters: the stringency τ (which controls the sensitivity/specificity) and the prior p (which controls the weight of a network, used to model e.g. miRNA-mRNA interactions and promoter structure). For an in-depth discussion of the model, bootstrap procedure and the prior, we refer to (Kling et al., 2015). In step 3 (subtype regulators) we analyze the network to predict direct and indirect regulators of subtypes. As direct regulators, we consider transcription factors or miRNAs that are directly linked to at least one subtype (i.e. they are found to contribute independently to subtype variability). As indirect regulators, we consider mutations, CNA loci, epigenetic changes or miRNAs connected to at least one subtype variable via at least one mRNA (Fig. 1c). The sign and connectivity of each regulator is summarized by its net connectivity, defined as the sum of positive and negative influences between a regulator and a subtype (Fig. 1c). To evaluate this pipeline as a tool to uncover subtype regulators, we performed cross-validation using two independent transcriptional data sets: the TCGA data with 529 GBM cases (Brennan et al., 2013) and the study by Phillips et al. of 100 glioma cases (Phillips et al., 2006). When compared to the Master regulator algorithm (MRA, a technique to predict TF (Carro et al., 2010) and miRNA (Cantini et al., 2015) regulators of cancer), aSICS was approximately twice as sensitive at any fixed level of false detections (p ≤ 0.001) (Fig. 1d). The number of predicted regulators will depend on the stringency τ (Fig. 1e). To set the stringency cross-validation can be used (Jornsten et al., 2011), or an empirical control. Here, we chose the latter strategy, and adjusted stringency to maximize the relative overlap between the network and our internal database of published candidate genes in GBM extracted by text mining and manual curation (Supplement) (Fig. 1f). In all, the proposed pipeline presents a comprehensive framework to predict several types of transcriptional and genomic determinants of cancer subtype, and their net effect.

New Regulators of Glioblastoma Subtypes
At the preferred stringency setting we found 12 direct and 38 indirect regulators of GBM subtypes, from an aSICS model with 3400 links (  Supplementary Fig. 1). Three trends emerged. Firstly, we confirmed known key subtype selective mutations, including NF1 (mesenchymal) and IDH1 (proneural), as well as amplifications of PDGFRA (proneural) and EGFR amplification (classical). The model confirmed multiple transcriptional regulators with key roles in GBM, in particular CEBPB, WWTR1, ASCL1, SOX4 and SOX10 (Carro et al., 2010;Glasgow et al., 2014;Su et al., 2014;Zhang et al., 2014), as well as the miRNA regulators miR-222 and miR-21 (Stewart et al., 2013) (Table 1). The previously reported mesenchymal driver STAT3, however, was not detectable in TCGA data. In contrast to the Phillips et al. material, from which it was first detected, STAT3 mRNA does not correlate with mesenchymal subtype in TCGA material (r = 0.13); this may reflect either a technical difference or the inclusion of lower grade glioma in the Phillips et al. cohort. Secondly, while multiple regulators were detected for the proneural and mesenchymal subtypes (Table 1), there were relatively few predictions for the classical subtype and not a single one for the neural subtype. The neural subtype has limited prognostic relevance and was recently found to be rare among primary cell cultures from GBM (Xie et al., 2015). While more research is warranted, one possible interpretation is that the neural subtype might be better explained by the presence of brain parenchyma in the TCGA samples where it was detected (Verhaak et al., 2010). The prediction that stood out for the classical subtype was methylation of the EVI2A gene promoter, in close proximity to the NF1 locus. This methylation event (cg23352695) may warrant further investigation as a possible determinant of classical subtype GBM. Thirdly, we predicted a notable number of epigenetically controlled regulators of the mesenchymal sub-type, proneural subtype, or both (with opposing signs). The methylation events included the recently reported suppression of podoplanin, PDPN (Peterziel et al., 2012), the aSICS detected promoter regulation of brevican (BCAN) (Lu et al., 2012), ANXA2, SOX10 and MYT (Table 1). Such detections are not a consequence of IDH1 mutation, since this mutational event was identified as an independent node in the model, and was thus accounted for. We further concluded that despite the tendency for DNA copy number aberrations and DNA methylations to appear across entire genomic regions, aSICS was generally successful in restricting the prediction to relevant CNAs (EGFR, PDGFRA) and cis-acting promoter methylations, and the detected regulators were not significantly clustered on any particular chromosome. This is due to the combined effect of rigorous variable selection and a network prior (this is further analyzed in (Kling et al., 2015)).

Evidence of ANXA2 as a Methylation Controlled Marker of Glioma Grade
The high number of predicted epigenetic regulators of the mesenchymal subtype might suggest a general principle for how GBM cells adopt more mesenchymal and invasive properties. We chose one such regulator, ANXA2, for further investigation. According to the aSICS analysis, ANXA2 was linked to the mesenchymal subtype and promoter methylation (cg08081036), positive and negative interaction sign, respectively (Fig. 1g, Supplementary Fig. 1). A member of the annexin family of Ca 2+ and membrane-binding proteins (Rescher and Gerke, 2004), it has been suggested as a biomarker for glioma grade stratification (Gao et al., 2013) and found to regulate cell migration and proliferation (Tatenhorst et al., 2006;Zhai et al., 2011). Interestingly, inspection of the regulators in Table 1 highlighted a number of functional partners and targets of ANXA2: S100A10 and PLAUR (controlled by miR-21 and miR-222, respectively) are involved in plasminogen to plasmin conversion, a reaction promoted by the presence of ANXA2 in the cleaving complex (Stewart et al., 2013), while the mesenchymal gene PDPN whose promoter methylation is positively connected to the proneural signature and which plays a role in tumor invasion, is controlled by the PI3K/AKT/AP1 pathway that can be activated by the binding of ANXA2 to tenascin C (Chung and Erickson, 1994;Gong et al., 2010;Peterziel et al., 2012). This observation further strengthened our Table 1 aSICS derived regulators of GBM subtypes. CL = classical subtype, MES = mesenchymal subtype, NL = neural subtype, PN = proneural subtype. Numbers in boxes are net connectivity between each regulator and the corresponding subtype. Color indicates positive (red) and negative (blue) net connectivity between regulator and subtype, respectively (Methods). hypothesis for a pivotal role of ANXA2 in mesenchymal tumors. To explore ANXA2 as a possible methylation-controlled regulator of mesenchymal transformation in GBM, we first analyzed ANXA2 expression and promoter methylation in a set of independent samples. Like in the TCGA cases (Fig. 2a,b), expression of ANXA2 was significantly higher in mesenchymal compared to non-mesenchymal GBM in an independently collected set of samples from the Freiburg University Medical Center (Fig. 2c,d). Consistent with a role of DNA methylation in ANXA2 gene expression regulation, ANXA2 promoter methylation in TCGA samples showed a clear correlation with grade, primary GBM exhibiting a lower degree of methylation compared to both LGG grad II and III (Fig. 2e), which resulted in a negative correlation between ANXA2 expression and its promoter methylation (Fig. 2f). Finally, we looked at the methylation status of the ANXA2 promoter (cg08081036) in a cohort of in-house high and low grade gliomas analyzed by methylation array and gene expression profile ( Supplementary  Fig. 2). Tumor sample classification based on expression profiling data allowed to characterize the low-grade tumors as proneural/neural and the high-grade ones as mesenchymal. We decided to consolidate the neural group with the proneural one as a classification distinct from the mesenchymal subtype since it is still controversial whether the neural subtype represents a distinct subclass (Verhaak et al., 2010). As expected, the proneural/neural low grade tumors showed higher global methylation compared to the high-grade/mesenchymal samples, consistent with previous findings (Noushmehr et al., 2010) (Supplementary Fig. 2a). Examination of the ANXA2 promoter (cg08081036) confirmed its methylated status in low-grade/proneural/neural samples ( Supplementary Fig. 2b). Analysis of ANXA2 expression by qRT-PCR confirmed its association with high-grade/mesenchymal gliomas (Supplementary Fig. 2c). Overall, our data indicate that ANXA2 presents a low promoter methylation and elevated expression in high-grade mesenchymal gliomas.

ANXA2 is Suppressed by Promoter Methylation and Correlates With Patient Prognosis
Since the network model and subsequent validation of ANXA2 promoter methylation in TCGA samples was based on a single probe (cg08081036), we performed pyrosequencing analysis on several low and high-grade tumor samples and brain tumor stem cells-like (BTSCs) from the Freiburg Medical Center. The analyzed region covered 7 CpG sites, including the one corresponding to probe cg08081036 (   Fig. 1. A novel pipeline to discover determinants of cancer subtypes, applied to glioblastoma. Pipeline principle. (a) In step 1, different types of data are integrated with molecular subtype information. (b) In step 2, aSICS is applied to obtain a robust network of sparse partial correlations. The user sets the stringency of the solution and adjusts the strength of prior information (e.g. miRNA-mRNA links) with two tuning parameters (p and τ). (c) In step 3, the best derived network is parsed to define direct (transcription factors and miRNAs) and indirect (miRNAs, DNA methylation events, DNA copy number aberrations and point mutations) determinants of subtypes. The indirect determinants are assigned a net signed connectivity to each subtype (asterisk: methylation event 'X' with connectivity score: −3 + 1 = −2). (d) Reproducibility of TF regulator detection between TCGA and Phillips et al. data sets using aSICS and the Master Regulator Algorithm (MRA), respectively. (e) Application to GBM data, showing increasing number of detected regulators, depending on stringency. (f) Enrichment of genes in a database of experimentally studied GBM genes is optimal for 30-50 detected regulators, suggesting that this is an optimal network size. (g) Examples of detected subtype regulators (c.f. Table 1). 3a and Supplementary Fig. 3). As shown in Fig. 3b, ANXA2 methylation was higher in a fraction of low grade and in secondary GBM, which were defined according to their likely evolution from lower grade tumors, and was generally consistent upon recurrency (Fig. 3b). In accordance with these findings, immunostaining of an independent set of 6 GBM and 6 LGG tissues sections showed a considerably higher expression of ANXA2 protein in GBM compared to LGG (Fig. 3c-d). Analysis of IDH1 status revealed that samples with low ANXA2 methylation were mostly IDH1 wild-type (19/22) (Fig. 3b). Survival analysis of the same group of Freiburg patients with known ANXA2 methylation status (Fig. 3b) showed a significant association of ANXA2 methylation with survival (p = 0.013) (Fig. 3e). Although limited in number, secondary GBMs seemed to be associated with a more favorable prognosis. Anecdotally, the only secondary GBM with low ANXA2 methylation had a survival time of only 15 months, further reinforcing the idea that ANXA2 methylation could be a prognostic factor (Fig. 3e). In accordance with this, and consistent with our network, we observed a clear association of ANXA2 methylation with survival in TCGA patients (Fig. 3f). Similarly to what was observed in the Freiburg sample set (Fig. 3b), ANXA2 methylation was significantly associated with IDH1 mutation in TCGA samples (Fig. 3g). To further confirm a role of DNA methylation in regulating ANXA2 expression, we treated the cell line BTSC23, characterized by a methylated ANXA2 promoter ( Supplementary Fig.  3), with the hypomethylating agent 5′-Aza-2′-deoxycytidine (5-Aza-2-dC) alone or in combination with the histone deacetylase inhibitor Trichostatin A (TSA), because histone deacetylation has been shown to cooperate with DNA methylation to epigenetically repress transcription (Nan et al., 1998). Our results showed increased expression of ANXA2 upon 5-Aza-2-dC treatment and a further increment with the combination of 5-Aza-2-dC and TSA, indicating that ANXA2 expression is epigenetically regulated (Fig. 3h). Since tumors often contain cells with distinct transcriptional signatures and such heterogeneity could have important implication on the effect of therapeutic treatments (Patel et al., 2014), we performed a regional analysis of ANXA2 expression within a cohort of 13 GBMs in which 3 samples per tumor were collected through localized MRI-guided sampling (Fig.3i). The analysis showed a significant degree of co-expression with mesenchymal marker CHI3L1 across biopsies (Fig.3j). There was, however, no gradient from central to peripheral tumor regions, suggesting that its expression is dependent on the GBM molecular signature rather than associated with any particular region of the tumor (Fig.3j). We further noted that, in the single cell material collected and profiled by Patel et al. (Patel et al., 2014), ANXA2 was predictive of mesenchymal signature (p b 10 −8 , Supplementary Fig. 4)  The edema zone is highlighted in light blue, the contrast-enhancing zone in red, and the core in black. (j) Box-plot graph correlating ANXA2 expression with GBM regions and CHI3L1 expression. The graph shows the smallest and largest observations (upper and lower whiskers, respectively), the interquartile range (box), and the median (black line); data points more than 1.5 times the interquartile range lower than the first quartile or 1.5 times the interquartile range higher than the third quartile were considered to be outliers. For each GBM region n = 13.
a b d f i c e g h j expression is associated with the higher GBM tumor grade and the mesenchymal subgroup and that ANXA2 is repressed in LGGs and secondary GBMs at least partially via DNA methylation. In addition, ANXA2 covaries with mesenchymal transformation within individual tumors, even down to a resolution of individual cells.

ANXA2 Regulates Mesenchymal Transformation of Brain Tumor Stem Cells (BTSCs)
To further investigate the role of ANXA2 in GBM, we selected two BTSC lines (BTSC161 and 168) in which ANXA2 was expressed to perform knockdown experiments (Fig. 4a,b, Supplementary Fig. 5a,b). Silencing of ANXA2 led to suppression of mesenchymal marker transcripts (Phillips et al., 2006;Verhaak et al., 2010), as measured by GSEA (Fig. 4c,d and Supplementary Fig. 5c). We decided to validate our results in BTSC168 since the phenotype appeared stronger in this cell line. A panel of mesenchymal genes which showed downregulation upon ANXA2 knockdown in the gene expression array, was then validated by qRT-PCR (Fig. 4e). Downregulation of the mesenchymal gene CHI3L1 upon ANXA2 silencing was also validated at the protein level ( Fig. 4f-h). Furthermore, ANXA2 overexpression in BTSCs lacking ANXA2 (3047) induced both mesenchymal genes ( Supplementary Fig.  6a,b) and, interestingly, ANXA2 itself, suggesting a positive feedback mechanism ( Supplementary Fig. 6c,d). A comparatively more favorable outcome for GBM has been associated to a subgroup of tumors, G-CIMP, characterized by a higher DNA methylation level (Noushmehr et al., 2010). This was also confirmed by our network model ( Supplementary  Fig. 1) in which patients survival was linked to methylation events, which in turn tended to be in the promoter region of G-CIMP reported genes (Fisher test p-values 2.5 · 10 −25 and 2.5 · 10 −61 for genes directly or one network step away from the survival node, respectively). Our data indicate that ANXA2 is methylated in the group of tumors (lower grade, mostly belonging to the proneural and neural subclass) which show high methylation of similar nature to the previously identified G-CIMP ( Supplementary Figs. 2,3). This is also in agreement with the study of Noushmehr and colleagues showing that ANXA2 is a G-CIMP gene (Noushmehr et al., 2010). Since we observed that ANXA2 silencing causes a loss of the mesenchymal gene signature, and knowing that most of the non-G-CIMP tumors are mesenchymal, we tested whether ANXA2 also affects the methylation and expression of G-CIMP genes. Although no significant change was observed in global methylation upon ANXA2 knockdown in two BTSC lines (data not shown and Supplementary Fig. 7a), GSEA showed a gene expression profile shift towards the G-CIMP signature; in particular, genes methylated ( Fig. 4i) or downregulated (Fig. 4j) in G-CIMP tumors were downregulated, whereas the expression of genes upregulated in G-CIMP tumors was enriched ( Supplementary Fig. 7b). Altogether, these results suggest that ANXA2 plays a role in the maintenance of the mesenchymal signature and that its removal from the cellular system leads to a gene expression profile similar to that of the less aggressive G-CIMP subgroup of GBM.

ANXA2 Acts Upstream of Known Master Regulators of Glioma Transformation
The expression of mesenchymal genes in GBM has shown to be controlled by a handful of transcription factors that act as master regulators of the gene expression profile (Bhat et al., 2011;Carro et al., 2010). Since our results implicated ANXA2 in the maintenance of the mesenchymal signature, we investigated whether it plays a role in the regulation of the mesenchymal master regulators. ANXA2 knockdown in BTSC168 (Fig. 3a,b) led to a downregulation of 4 out of 7 tested master regulators, although STAT3 and CEBPB which have been shown to drive the expression of the other regulators, did not significantly change (Fig. 5a). Then, we decided to look at STAT3 phosphorylation, which has been demonstrated to be required for its activation in GBM (Bhat et al., 2011). Here, two new primary cell lines (BTSC380 and BTSC407) were used, since the previous cells (BTSC161 and 168) were no longer available due to a loss of their previous phenotype and to a decrease in cell growth rate. Interestingly, silencing of ANXA2 in BTSC380, BTSC407, and in the GBM cell line SNB19 resulted in decreased STAT3 phosphorylation (Fig.  5b,c). The halt to the activation of STAT3 pathway upon ANXA2 knockdown was then confirmed by GSEA, showing a loss of expression of genes upregulated by STAT3 and conversely, an enrichment of genes downregulated by STAT3 (Fig. 5d). These data suggest that ANXA2 is controlling the maintenance of the mesenchymal signature by regulating the expression and/or the activation of a subset of mesenchymal master regulators.

ANXA2 Suppression Abrogates Mesenchymal Characteristics of Brain Tumor Stem Cells
Upon ANXA2 silencing, two primary BTSC lines (BTSC168 and BTSC407) lost their characteristic elongated shape and their dispersed distribution, tending instead to cluster together to form large spheres (Fig. 6a, Supplementary Fig. 8a), a phenotype reminiscent of mesenchymal-to-epithelial transition described in several different cancer types (Baum et al., 2008). To support this observation, stainings for the stemness and EMT marker nestin (NES), and for the EMT markers MMP2 and vimentin (VIM) showed an attenuated protein abundance upon ANXA2 silencing in the mesenchymal primary cell line BTSC407 (Fig. 6b, Supplementary Fig. 8b-d). Moreover, ANXA2 modulated both cellular proliferation, as shown by EdU incorporation assay and Ki67 staining ( Fig. 6c-d, Supplementary Fig. 8e-g), and migration, as shown by a matrigel invasion assay (Fig. 6e,f). Finally, ANXA2 silencing rendered BTSCs incapable of mesenchymal differentiation into osteocytes when exposed to differentiating conditions (Fig. 6g). Thus, ANXA2 knockdown in BTSCs seems to alter the gene expression profile of the cells, inducing drastic changes of important cellular properties such as proliferation, invasion, and cell fate.

Discussion
The progression of a GBM depends on invasive growth marked by mesenchymal cell features. It therefore remains a critical priority to expand our knowledge of how the mesenchymal features are modulated. We have employed aSICS as an integrative modeling method to identify regulators of cancer subtypes. Applied to GBM, our model implied that a number of methylation events, independently of IDH1 mutation, contribute to a mesenchymal/proneural axis in GBM. We confirmed ANXA2 as a key regulator of mesenchymal transformation and demonstrated its importance for viability, invasiveness and maintaining a mesenchymal gene signature. The results thus show that integrative modeling can uncover a new mesenchymal modulator. Additional predictions provide a rich source for future investigation (Table 1). Our results warrants further investigation of aSICS as a general tool to uncover cancer subtypes. Compared to existing approaches, aSICS includes broader spectrum of data types, and benchmarking confirmed reproducible performance between independent cohorts. Explorative analyses of three additional cancers reveal that TF, miRNA, mutations and DNA methylation regulators are also detected in breast, ovarian and colorectal cancers, suggesting that new regulators can be predicted (Supplementary Fig. 9). A second important venue of investigation will be to explore the impact of alternate classification schemes, or using the model to improve subclassifications. For instance, the absence of regulators of the neural subtype is notable and may indicate a lack of mechanistic support for this subtype. Furthermore, there are currently inconsistent observations regarding the prognostic value of mesenchymal signatures from individual biopsies (Noushmehr et al., 2010;Ozawa et al., 2014;Phillips et al., 2006;Verhaak et al., 2010). These observations motivate extended method development and broader applications of the aSICS framework, also taking into account tumor heterogeneity, reserved for future work. Furthermore, the bootstrapping framework can be generalized, e.g. to compute confidence intervals of connectivity scores, reserved for future work. The software itself is currently available in Matlab and can be obtained from the authors upon request. In addition to the finding observation that ANXA2 is a key regulator of mesenchymal targets, our analysis adds to the characterization of ANXA2 as a possible biomarker in solid tumors (Liu et al., 2015;Zhang et al., 2012). Based on our results, ANXA2 promoter methylation shows promise as a prognostic marker, and we also find that ANXA2 expression levels to be predictive of patient survival. When extending the analysis to lower grade glioma (LGG) and secondary GBM, which tend to be IDH1 mutant and hypermethylated, we noted elevated ANXA2 promoter methylation and suppressed expression. In the G-CIMP signature described by Noushmehr et al. (2010), ANXA2 is a methylated gene. Together, these observations imply that ANXA2 can be suppressed by IDH1 mutation in LGG but also that IDH1-independent mechanisms can modulate ANXA2 via methylation of its promoter in GBM cells. Although we did not observe any significant effect of ANXA2 on global DNA methylation, several G-CIMP target genes were among the genes affected by ANXA2 knockdown. Thus, it is possible that ANXA2 knockdown leads to a general shift in gene expression that effectively mimics the expression modulation normally achieved through DNA methylation in G-CIMP tumors. Since G-CIMP tumors are associated with the most favorable patients' prognosis (Noushmehr et al., 2010), we speculate that the acquisition of a G-CIMP-like gene expression signature upon ANXA2 knockdown in BTSCs might explain the connection between ANXA2 methylation and survival in our two cohorts. IDH1 mutation correlates significantly with both ANXA2 methylation or expression. However, as estimated by a correlation analysis on TCGA cases, IDH1 status only explains circa 62% of the variation of ANXA2 promoter (cg08081036) methylation and 26% of the expression variation. Thus, removing IDH1 mutant cases from the analysis, a significant degree of correlation between ANXA2 expression and methylation is retained (r = −0.21, p = 0.03). We further find that after removal of IDH1 mutant cases, ANXA2 expression is predictive of longer survival ( Supplementary Fig. 10A). In addition, in the Freiburg patient survival analysis (Fig. 4e) ANXA2 methylation was not always associated with IDH1 status. In fact, in the group of secondary GBM (shown in Fig. 3b and in the survival analysis (Fig. 3e, indicted by a "x") two patients with ANXA2 methylation but wildtype IDH1 status showed nonetheless a more favorable prognosis. Combination of ANXA2 expression with recently described prognostic indicators 1p19q co-deletion and TERT promoter mutation (Eckel-Passow et al., 2015) also suggested that ANXA2 expression adds independent prognostic power ( Supplementary Fig.  10B), motivating combined analysis of ANXA2 and other indicators in larger cohorts. Together, our results therefore imply that while IDH1 mutation is a key mechanism behind ANXA2 suppression, particularly in LGG and secondary GBM, it is likely regulated by additional factors, and may have clinical promise as a prognostic complementary to IDH1 mutation, or as an indirect marker of G-CIMP cases. Our MRI-localized regional analysis of ANXA2 expression in GBM showed a variation across the tumor, with elevated expression in regions expressing a mesenchymal markers. This finding, in turn, was consistent with a re-analysis of the data by Patel et al, in which ANXA2 expression predicts mesenchymal transformation in individual cells (c.f. Supplementary Fig. 4). Reduced ANXA2 expression in the edema zone compared to contrast-enhancing zone and core of mesenchymal tumors might be related to the particularly heterogeneous nature of this region in which ANXA2-expressing neoplastic cells are intermixed with non-neoplastic cells without ANXA2 expression (Gill et al., 2014). GBM tumors often contain cells with distinct transcriptional signatures; such heterogeneity might reflect the evolution process of tumor cells by the accumulation of temporal mutation events (Ozawa et al., 2014;Patel et al., 2014;Sottoriva et al., 2013). The different level of ANXA2 expression in the tumor may thus suggest that it is an interesting target to re-program mesenchymal subsets of GBM cells.
Functional experiments in BTSCs showed that ANXA2 is a causative regulator of mesenchymal transformation, which both regulates and co-depends on previously described regulators. Firstly, ANXA2 knockdown and overexpression in BTSCs suppressed and induced mesenchymal genes, respectively. In BTSC168 ANXA2 knockdown resulted in a more pronounced loss of the mesenchymal signature by Phillips and Verhaak compared to the other tested cell line, BTSC161 (Fig. 4c, consistent with previous indication that glioblastoma tend to shift to a mesenchymal phenotype upon recurrence (Phillips et al., 2006). Exogenous expression of ANXA2 appears to induce a positive feedback loop that induces endogenous ANXA2 in BTSCs. This is interesting in the light of ANXA2 presence in microvesicles secreted by GBM cells which are able to deliver material into neighboring cells (Bronisz et al., 2014). Potentially, such a positive feedback loop might represent a mechanism of paracrine interaction between tumor cells to coordinate cellular behavior across the cell population. Effects of ANXA2 on cellular properties in GBM has been previously reported ( Tatenhorst et al., 2006;Zhai et al., 2011), although these studies only used human GBM cell lines or murine models. Here, we investigated the ANXA2-dependent mesenchymal signature modulation effect on cell behavior in BTSCs. As previously described for GBM cell lines (Tatenhorst et al., 2006;Zhai et al., 2011), ANXA2 expression levels also determine cellular proliferation and invasiveness in BTSCs. In BTSCs, loss of the mesenchymal signature prevented the cells from differentiating into mesenchymal progeny. Thus perturbation of ANXA2 translates into an actual change in cell phenotype by affecting key properties such as proliferation, invasion, and notably, cell fate. Consistent with the phenotypic response to ANXA suppression (Fig.  6), suppression of ANXA2 led to an alteration of not only mesenchymal signatures but also cell cycle related pathways, such as cell cycle checkpoint genes (Gene Set Enrichment Analysis p = 8.6 × 10 −7 , q = 8.6 × 10 −6 , Supp Table 1). Intriguingly, additional pathways are also detected as significantly affected by ANXA2 perturbation, including e.g. Tumor Necrosis Factor Signaling (p = 1.2 × 10 −8 , q = 2.08 × 10 − 7 ) and MHC class II antigen presentation (p = 1.02 × 10 − 4 , q = 1.38 × 10 −2 ). Thus further studies will be needed to elucidate the full cellular impact of ANXA2 knockdown. In order, an extended study of ANXA2 knockdown in more cell lines is motivated, e.g. to determine differences in vulnerability to ANXA2 suppression. Such a systematic study appears quite possible, since ANXA2 is broadly expressed in glioma cells and does not co-vary significantly with cell culture passage (Supp Fig 11).
Our data provides new evidence of a functional interaction between ANXA2 and transcription factors (TFs) known to be involved in GBM mesenchymal programs. Previously, STAT3 phosphorylation has been linked to regulation of three TFs that regulate mesenchymal genes: RUNX1, FOSL2, and BHLHB2 (Carro et al., 2010). It is notable that ANXA2 knockdown in BTSCs suppresses both STAT3 phosphorylation and all three target transcription factors. This would place ANXA2 upstream in a regulatory hierarchy. Interestingly, ANXA2 is not categorized as a transcription factor (TF), the class of protein most commonly associated with the gene expression master regulator function (Bhat et al., 2011;Carro et al., 2010). Danussi et al. identified the RhoA-binding protein, RHPN2, as an example of mesenchymal signature driver that does not belong to the TF category (Danussi et al., 2013). Similarly, ANXA2 might represent a mechanism to control cellular programs alternative to the canonical TF drive. Further studies would be required to fully understand the ANXA2 downstream activation mechanism. Since ANXA2 cellular localization is more accessible than that of TFs, which generally pose as difficult drug targets, our identification of ANXA2 as a regulator of the mesenchymal signature might represent an opportunity to therapeutically target this molecule in order to more efficiently treat brain tumors and improve patients' survival.

Author Contributions
TK, RJ, PJ, and SN were responsible for the computational part of the study, and RF, DOA, DHH, FD, IV and MSC were responsible for the experimental part. AW provided surgical samples and supervised the collection of MRI-localized GBM biopsies. TK, RF, MSC, and SN planned the work and wrote the manuscript. All authors participated in interpreting the results and refining the manuscript.