A detailed description of the computational algorithms used can be found in the Supplementary Methods section. The Multi-Reg algorithm was developed for this research, and the software is available at: http://www.c2b2.columbia.edu/danapeerlab/html/software.html.

SF188, SNB19 and 293T cell lines were grown in DMEM plus 10% FBS (Gibco/BRL). Primary human astrocytes (Lonza Inc.) were grown in Astrocyte Medium (SciencCell). Mouse NSCs (clone C17.2)(14) were cultured in DMEM plus 10% heat-inactivated FBS (Gibco/BRL), 5% horse serum (Gibco/BRL) and 1% l-glutamine (Gibco/BRL). Neuronal differentiation of mouse NSCs was induced by growing cells in DMEM supplemented with 0.5% horse serum.

Lentiviral expression vectors pLOC RHPN2 and pLOC VEC; and lentiviral vectors carrying RHPN2 shRNAs were purchased from Thermo Scientific Open Biosystems. To generate lentiviral particles, each expression plasmid was co-transfected with pCMV-dR8.91 and pCMV-MD2.G vectors into human embryonic kidney 293T cells using Fugene 6 (Roche). Lentiviral infections were performed as previously described (15).

RNA was prepared with RiboPure kit (Ambion), and used for first-strand cDNA synthesis using random primers and SuperScript II Reverse Transcriptase (Invitrogen). qRT–PCR was performed using Power SYBR Green PCR Master Mix (Applied Biosystems). Primers are listed in Supplementary Data, Table S1. qRT–PCR results were analysed by the ΔΔC_T method (16) using 18S as a housekeeping gene.

RNA amplification for mouse array analysis was performed with Illumina TotalPrep RNA Amplification Kit (Ambion). One-and-a-half micrograms of amplified RNA was hybridized on Illumina Mouse ref8 v2 expression BeadChip according to the manufacturer’s instructions. Hybridization data was obtained with an iScan BeadArray scanner (Illumina) and pre-processed by variance stabilization and robust spline normalization implemented in the lumi package under the R-system (17).

Cells were grown on polylysine (Sigma) treated glass cover slips and fixed with PBS 4% PFA for 15 min. Then, they were permeabilized (with PBS, 1% BSA, 0.1% Triton X-100, and 2% FCS) for 5 min and saturated with the blocking buffer (PBS, 1% BSA, and 2% goat serum) for 30 min. The primary and the following secondary antibodies were incubated at room temperature for 1 h. These were: SMA (mouse monoclonal, Sigma), βIII-tubulin (mouse monoclonal, Promega), fibronectin and Paxillin (mouse monoclonal, BD Biosciences), phospho-Cofilin (rabbit monoclonal, Cell Signaling), goat anti-mouse and anti-rabbit Cy3 conjugated (Life Technologies). Actin cytoskeleton was stained with Alexa-Fluor 568 Phalloidin (Life Technologies) and nuclei were visualized with DAPI (Invitrogen). Images were acquired with Nikon A1R MP confocal microscope. Quantification of the fibronectin intensity staining in mouse NSCs was performed using NIH Image J software (http://rsb.info.nih.gov/ij/).

Exponentially growing cells were serum starved for 24 h, detached with Accutase solution (Innovative Cell Technologies Inc.) and adhered to FN-coated dishes (10 µg/ml; Sigma) for 30 min. Then, cells were lysed and 300 µg of protein was tested in a rhotekin-RBD bead pulldown assay (Rho Activation Assay Biochem Kit, Cytoskeleton Inc.) for 1 h at 4°C. After thorough washes, the samples were boiled for 5 min in Laemmli buffer to detach active GTP-bound Rho and then loaded on 4–20% SDS-PAGE gels (Invitrogen) and immunoblotted using an anti-RhoA antibody.

An emerging trend in cancer treatment is drugs that target genes and signaling pathways that are only activated in specific cancer cells (18, 19). However, genomics has revealed incredible heterogeneity in cancer, which makes identification of specific genes contributing to cancer progression (driver genes) difficult. Targeting novel drivers is especially important in GBM, since the median survival with conventional therapy is 12–15 months (5).

A key challenge in identifying driver genes from DNA copy number is that amplification and deletions frequently involve large regions of DNA, each consisting of multiple genes. To pinpoint the driver genes within such genomic regions, we previously developed CONEXIC (20), a computational algorithm that integrates copy number and gene expression data, to identify driver genes and connect these to their expression signatures. A key limitation to our previous approach is that CONEXIC can only identify the one dominant driver for each expression signature. However, multiple drivers can contribute to the same effect, sometimes acting in parallel. For example, almost all GBM patients have activated Receptor Tyrosine Kinase signaling and disrupted p53/RB signaling, but each patient could have a different combination of deletions, amplifications and mutations in some of the many genes known to influence these signaling pathways (21).

We therefore developed Multi-Reg, based on a similar framework as CONEXIC, but seeking multiple regulators for each phenotype. This improvement was achieved by a change in the main statistical model. While CONEXIC and other methods begin from gene expression signatures and attempt to find a driver for each one (22), Multi-Reg begins from candidate drivers and then finds a signature associated with each driver.

Multi-Reg begins from regions that are significantly altered in copy number, either amplified or deleted (Fig. 1A). It identifies all genes in each region as candidate driver genes (Fig. 1B). Then for each candidate driver it generates its gene expression signatures. i.e. the list of candidate target genes associated with this driver. Comparing the expression signatures between drivers from the same region allows us to focus on significant drivers (Fig. 1C). The final step of the analysis involves assigning the expression signature of each predicted driver to a distinct subtype of GBM (MES, Proliferative, Proneural). This leads to the testable hypothesis that some genes altered in copy number are drivers of distinct biological functions and the same region may contain multiple drivers that effect distinct GBM subtypes (Fig. 1D). More details on the algorithm can be found in the Supplementary Methods.

We analyzed gene expression and copy number measurements from 136 samples (corrected for batch effects) of primary GBM from the TCGA cohort (21), see Supplementary Methods for details. We identified 238 regions that are recurrently altered in copy number using a modified version of GISTIC called JISTIC (23) and generated a list of 747 candidate drivers contained within these regions. Next, we applied Multi-Reg to integrate copy number with gene expression and mutation data to pinpoint the top ranking candidate drivers and their targets. Applying Multi-Reg resulted in the identification of 83 high-scoring drivers, which associate with a total of 12125 targets.

Multi-Reg identified many of the well-known oncogenes and tumor suppressors in GBM, including EGFR, NF1, CKDN2B, p53, PIK3CA, RB1, PTEN (21), and more (see Table S2 for complete list). The successful identification of known key brain tumor drivers increases our confidence in the novel predictions discussed below. To determine the functional properties of the identified drivers, we compared the gene expression signatures previously associated with MES, Proneural and Proliferative GBM subtypes (12) with the gene modules associated with each driver identified by Multi-Reg, using hypergeometric enrichment.

In addition to well-known drivers of GBM, we identified new drivers that regulate the GBM signatures. Overall, out of the 83 drivers we identified, 23 drivers were associated with the MES, 14 with the Proliferative and 11 with the Proneural signature (Table S2 for all drivers and see Fig. 2A and B for a schematic representation of these 48 drivers and the chromosomal location of the drivers associated with the MES signature). Thus our algorithm is able to identify drivers for each of the three crucial GBM subtypes.

We focused on the genes implicated in the Mesenchymal subtype, since this subtype has the worst clinical prognosis. We sorted the 23 genes associated with the MES signature according to q-value. The top three genes, ERBB2, COL1A1 and ITGB3, were mutated, but they did not harbor copy number changes. The next two top genes, C5orf32 and RHPN2, were both amplified and overexpressed. The RHPN2 gene is located on Chr19q12-q13 and, according to data from TCGA and Rembrandt databases, it is amplified and overexpressed in >30% of gliomas (Fig. S2). This gene codes for a RhoA-binding protein, called Rhophilin 2, which is not well characterized. We focused further experiments on RHPN2, given its potential biological function as regulator of Rho GTPases (24), which are key factors for cell migration and invasion, two hallmarks of the MES phenotype.

To further investigate the genes associated with this driver, we used Gene Ontology enrichment to analyze the function of predicted up- or down-regulated target genes identified by Multi-Reg for RHPN2 (Fig. 3). We found that RHPN2 up-regulated genes were significantly enriched for plasma membrane, extracellular region, transmembrane transporter activity, regulation of protein kinase cascade and cell adhesion (Fig. 3A and B), indicating that RHPN2 regulates the expression of genes involved in cell-ECM interactions, which is consistent with the induction of a MES phenotype. Genes down-regulated by RHPN2 were significantly enriched for the Gene Ontology terms such as development, neuron projection morphogenesis, regulation of gene expression and glial cell differentiation (Fig. 3C and D) and for the Proneural signature (q-value = 5.33 × 10⁻¹³). Thus, the Gene Ontology results match the predictions based on GBM signatures, RHPN2 induces a MES phenotype while repressing neuronal and glial cell differentiation.

A survival study performed in the Rembrandt database supported the importance of RHPN2 in GBM and encouraged us to investigate the role of this protein. According to this analysis, RHPN2 amplification and overexpression predicts a markedly poor clinical outcome of glioma patients (Fig. 3E). Among the top ranked drivers selected by Multi-Reg and associated with the MES genes, RHPN2 displayed the strongest statistical association with patient survival (data not shown).

To experimentally validate the new driver for GBM aggressiveness inferred by the Multi-Reg algorithm, we performed biological assays, overexpressing and/or silencing RHPN2 in different cell models. We started by testing the expression of specific genes and proteins that are representative of GBM phenotypes, performed functional assays and finally explored in detail the mechanism of action of RHPN2.

We used recombinant lentiviruses to express RHPN2 in C17.2 (Fig. S3A), a mouse immortalized neural stem cell (NSC) line previously used to investigate MES transformation of high-grade glioma (12). Consistent with the computational predictions, microarray gene expression analyzed by GSEA showed that expression of RHPN2 in C17.2 led to significant enrichment of MES genes (Fig. 4A; p-value = 0.0003). This result was validated by quantitative-RT-PCR of a representative panel of MES genes. Acta1, Acta2, Ctgf, Tnc, SerpinE1, Itga7, Osmr and C1rl were up-regulated after ectopic expression of RHPN2 (Fig. 4B). In addition, RHPN2 triggered the expression of MES proteins as shown by immunostaining for the MES markers fibronectin (FN) and smooth muscle alpha actin (SMA, encoded by the Acta2 gene; Fig. 4C–F; p-value < 1×10⁻⁴).

Next, we asked whether RHPN2 altered the default neuronal differentiation pathway of NSCs. Normally, C17.2 NSCs undergo neuronal differentiation upon mitogen removal (14). However, ectopic expression of RHPN2 conferred a fibroblast-like morphology (Fig. 4G) and inhibited neuronal differentiation, as evidenced by a decrease in βIII-tubulin immunostaining (Fig. 4H–I, p-value < 1×10⁻⁴). These findings suggest that RHPN2 blocks neuronal differentiation by reprogramming NSCs towards an aberrant MES lineage.

To elucidate the full scope of the biological phenotypes triggered by RHPN2, we used both gain and loss of function experiments. Based on the Multi-Reg predictions about the MES phenotype, we focused on invasion and migration assays. In addition to mouse NSCs we used primary Human Astrocytes (HA), candidate cells-of-origin for GBM, and two human glioma cell lines, SF188 and SNB19, which display the lowest and highest RHPN2 expression levels among several glioma cell lines respectively (Fig. S3B). One of the most distinguishing features of MES transformation is increased invasiveness (12). Indeed, expression of RHPN2 in C17.2 cells promoted invasion through the extracellular matrix in a Matrigel invasion assay (Fig. 5A) and enhanced cell migration in a wound assay (Fig. 5B). Interestingly, ectopic expression of RHPN2 (Fig. S3C) resulted in dramatic changes of the morphology of HA. Compared to control vector infected cells, HA expressing RHPN2 acquired an elongated and spindle-shaped morphology (Fig. S3D). Notably, these changes were associated with increased ability to invade through Matrigel (Fig. 5C; p-value < 0.001). A significant gain of invasion was observed also when RHPN2 was introduced in SF188 cells (Fig. S3E and 5D, p-value < 1×10⁻⁵). Conversely, silencing of RHPN2 in the human glioma cell line SNB19 by 4 different shRNA sequences (Fig. S3F) significantly decreased the invasive capacity through Matrigel (Fig. 5E; p-value < 0.005). Neither the ectopic expression of RHPN2 in HA and SF188 or its silencing in SNB19 affected cell proliferation (Fig. 5F–H).

Altogether these data indicate that RHPN2 overexpression is sufficient and necessary to drive transformation of glioma cells along the MES lineage.

We sought to identify the molecular mechanisms by which RHPN2 amplification/overexpression generates MES transformation of glial cells. To address this question, we first asked whether RHPN2 affects RhoA activity in primary human astrocytes (HA). We selected astrocytes as appropriate cellular models for the next series of experiments, since they are candidate cells-of-origin for GBM and display low rates of invasion. Ectopic expression of RHPN2 (Fig. S3C) resulted in more than 4-fold increase of RhoA activity (measured as RhoA-GTP), when compared with control cells (Fig. 6A).

Next, we asked whether the increase in RhoA activity by RHPN2 affected organization of the actin cytoskeleton, a key determinant of cell migration and invasion. After 30 min of adhesion on fibronectin, vector-transduced HA displayed a well-organized network of actin stress fibers that were distributed in the whole cell body (Fig. 6B, upper panel). RHPN2-expressing cells manifested a profound reorganization of the actin cytoskeleton, characterized by lack of stress fibers and accumulation of actin at the cell border in a ring-shaped manner, resembling a non-polarized lamellipodium (Fig. 6B, lower panel). Immunostaining for Paxillin, a protein implicated in cell migration (25), showed the expected localization of this protein at the focal contacts in control cells (Fig. 6B, upper panel). However, in HA transduced with the RHPN2-expressing lentivirus Paxillin formed concentric circles at the cell edge, co-localizing with actin (Fig. 6B, lower panel). Moreover, the aberrant activation of the RhoA pathway by RHPN2 resulted in the accumulation of membrane speckles of phospho-Cofilin (pCofilin; Fig. 6B), a downstream effector of RhoA (26).

Finally, we asked whether the increased RhoA activity induced by RHPN2 was responsible for the enhanced cell invasion. Treatment of vector and RHPN2-expressing cells with a specific RhoA inhibitor (named RhoA inhibitor I or C3 exoenzyme) completely abrogated RhoA activity after 6 hours (Fig. 6C). An invasion assay revealed that inhibition of RhoA did not affect the basal invasion capacity of HA transduced with the control lentivirus, but it completely reversed the RHPN2-induced invasion (Fig. 6D; p-value < 3×10⁻⁵). Together, these results indicate that activation of RhoA is specifically recruited to enhance invasion of HA following ectopic expression of RHPN2.