Coexpression Network Analysis Identifies Transcriptional Modules Related to Proastrocytic Differentiation and Sprouty Signaling in Glioma

Gliomas are the most common primary brain tumors in adults (1). Despite advances in their surgery, radiotherapy, and chemotherapy, the majority of gliomas remain lethal within 1 to 3 years from the time of diagnosis (2). A key challenge for diagnosing and treatment of glial tumors is their high interindividual heterogeneity (3). Histologically, gliomas are classified into 2 major subtypes: astrocytoma and oligodendroglioma (1), which are further divided into grades ranging from I (benign tumors) to IV (malignant glioblastoma) according to presence of anaplastic features such as necrosis, mitotic figures, vascular proliferation and others (1). However, tumors falling within the same category often substantially differ in their clinical behavior which reflects their diversity at the molecular level (4). Therefore, molecular approaches for investigation of the tumors heterogeneity and pathogenesis are required (3).

Microarray gene expression profiling offers an opportunity for genome-scale, quantitative evaluation of glioma biology by simultaneously measuring expression levels for thousands of genes (5). Molecular subtypes of glioma determined according to tumor expression profiles were shown to correlate with survival more accurately than does conventional histology (6, 7). Due to this higher accuracy in capturing the tumor phenotype, expression-based classifications offer an opportunity to improve clinical diagnostics of gliomas (6–8). Gene expression analyses have also yielded insights into molecular aspects of tumor progression (9), recurrence (8), and differentiation (10). However, a key challenge remains the high intrinsic complexity of the glioma transcriptome, which encompasses expression profiles of thousands of genes driven by a wide range of factors. This complexity hinders development of a standard expression-based classification of gliomas, different studies suggesting from 3 to 6 transcriptional glioma subtypes (6, 8, 11, 12), and leads to the need for an integrated view of expression variation in the tumors. We hypothesized that a more structured approach may facilitate understanding of the glioma transcriptome and provide insight into molecular diversity and pathogenesis of glial tumors.

Analysis of gene coexpression relationships is a powerful method for exploring intrinsic organization of a transcriptome (13). In contrast to the conventional differential expression analysis which searches for expression changes of individual genes between given conditions, coexpression networks explore gene-to-gene relationships and reveal intrinsic modules of coordinately expressed genes in an unsupervised way (14). This delineates expression changes that are driven by distinct cellular processes and provides a functional framework for understanding gene expression patterns of the tumors (13, 15). Several studies based on this approach yielded insights into oncogenic signaling pathways (15), chemotherapy resistance (16), and molecular heterogeneity of glioma (17, 18). Although resulting in important findings, these studies were limited in their ability to provide an integrated view of the structure of the tumor transcriptome due to methodologic constraints: 1) incomplete genome coverage (17, 18), 2) restriction of the analysis to only glioblastoma tumors (15, 16, 18), 3) use of algorithms that produce “tight” coexpression modules with a limited gene content (16, 17).

Here, we apply improved methods based on weighted gene coexpression network analysis (WGCNA; ref. 14) to publicly available microarray data covering multiple histologic types of glioma and the majority of genes in the genome. We identified genome-scale modules of coexpressed genes with clear functional interpretations and obtained an overview of the higher order organization of the transcriptome based on module-to-module relationships. The analysis revealed a novel prognostic signature which was associated with a proastrocytic pattern of tumor cell differentiation and provides potential gene markers for tumor subtype diagnostics. Our results also suggest a functional link between Sprouty proteins and oncogenic epidermal growth factor receptor (EGFR) signaling in glioma.

We analyzed 5 microarray data sets from 4 gene expression studies in glioma (6, 8, 11, 12). The data sets were originally generated using Affymetrix arrays (U133 Plus 2.0 and U133A platforms). We labeled the data sets according to the names of the first authors in the studies: “Gravendeel” (6), “Li_Train” (12), “Li_Test” (12), “Phillips” (8), and “Freije” (11). The primary analysis was focused on the Gravendeel data set as it contained the largest number of samples. Normalization of the raw data sets was performed using MAS5 algorithm (19) followed by quantile normalization. See Supplementary Methods for further information on the data sets and their normalization.

WGCNA (14) was performed in the Gravendeel data set. A coexpression network was constructed on the basis of 4,000 genes with highest expression variance (20). For all possible pairs of the variable genes, Pearson correlation coefficients were calculated across all samples. The correlations matrix was raised to the power 9, thus producing a weighted network (weighted correlation = correlation⁹, Supplementary Methods; ref. 14). The weighted network was transformed into a network of topological overlap (TO)—an advanced coexpression measure that considers not only the correlation of 2 genes with each other, but also the extent of their shared correlations across the weighted network (14). Genes were hierarchically clustered on the basis of their TO. Modules were identified on the dendrogram using the Dynamic Tree Cut algorithm (21). The resulting modules were tested for sensitivity to random removal of samples, and the stable modules were retained (see Supplementary Methods for details).

For each gene, we determined its connectivity within its module of residence by summing up the TOs of the gene with all the other genes in the module. By definition, highly connected (hub) genes display expression profiles highly characteristic for their module of residence (22). To obtain a condensed representative expression profile of each module (module eigengene, ME; ref. 22), we summarized expression levels of the top 10 hub genes in the module (Supplementary Methods; ref. 14). The resulting MEs were used to extend the modules' gene composition from the 4,000 network genes up to the genome scale (∼19,000 genes) using a newly developed procedure (Supplementary Methods).

Functional annotation of the modules was performed on the basis of analysis of their gene composition. We used DAVID (http://david.abcc.ncifcrf.gov/) to test each module for enrichment in genes with particular GO, SwissProt, and InterPro terms compared with the background list of all genes on the array (23). Association of modules with genomic rearrangements was detected using DAVID on the basis of overrepresentation of genes encoded at neighboring chromosomal locations (23). Relation of modules to differentiation programs and cell types was detected on the basis of enrichment in genes preferentially expressed in purified populations of neurons, astrocytes or oligodendrocytes (Fisher's exact test; ref. 24). In addition, we related modules to processes in glioma on the basis of the literature data.

For each gene, data were normalized into Z-scores before clustering. In the Gravendeel data set, we performed a k-means clustering of samples into 5 groups on the basis of the top 25 hub genes from the consensus modules M8, M11, M15, M20 (a total of 100 genes, Euclidean distance). For each cluster, a centroid was obtained that consisted of average expression values of the 100 genes. Tumors in the 5 independent data sets were assigned to classes according to the nearest centroid as measured by Euclidean distance between the centroids and expression levels of the 100 marker genes.

To compare our subtypes of glioma with those proposed by Phillips and colleagues (8), we reproduced the centroids (Mes, Prolif, and PN) described in the original study by calculating mean expression levels of the published centroid genes (8) in the published groups of samples (8) within the Phillips data set. Samples in each of the 5 microarray datasets were assigned to a class on the basis of the nearest centroid (Euclidean distance).

Gene set enrichment analysis (GSEA; ref. 25) was used to compare tumors in the “cluster 17” (38 samples) versus all the other tumors (238 samples) in the Gravendeel data set. The analyzed gene set (424 genes) was compiled of genes that showed at least a 5-fold upregulation in astrocytes compared to neurons and oligodendrocytes (24). The analysis was performed by randomly permuting phenotype labels 1,000 times (25).

From each data set, we selected differentiated tumors, that is those where the proastrocytic (M15) and/or neurogenesis (M20) modules were upregulated (Z-score normalized expression values of the proastrocytic ME and/or proneural ME higher than 1.0). Pearson correlation coefficients were calculated between the proastrocytic ME and the neurogenesis ME across these tumors in each data set.

To investigate functional organization of the glioma tran scriptome, we analyzed the largest-to-date public microarray glioma data set generated by Gravendeel and colleagues (Gravendeel data set; ref. 6). The data set consisted of 276 tumor samples including glioblastoma, anaplastic astrocytoma, anaplastic oligodendroglioma, low-grade astrocytoma, and low-grade oligodendroglioma (Supplementary Table S1). We applied WGCNA algorithm to explore transcriptional relationships between genes in this data set (14). A weighted gene coexpression network was constructed on the basis of 4,000 highly informative genes (14, 20). Because WGCNA evaluates coexpression for a pair of genes based on their topological overlap (TO), which considers not only the correlation of the 2 genes with each other, but also the degree of their shared neighbors across the entire network, the analysis highlighted consistent patterns of gene coexpression in the transcriptome (see Methods). Hierarchical clustering of the TO network followed by a screening for cluster stability revealed 20 stable modules (Fig. 1; modules M1–M20), each of them containing coordinately expressed genes potentially involved in shared cellular processes.

To determine gene composition of the modules at a genome scale, we first obtained a condensed expression profile for each module (i.e., ME) by summarizing expression patterns of the module hub genes into a single characteristic expression profile (14, 22). The resulting MEs captured the major expression trends of the modules and served to represent each module in a summarized form. We used the MEs in a newly developed procedure to extend the modules from the 4,000 seed genes to the entire data set of about 19,000 genes (Supplementary Methods), which increased the modules' size on average 3-fold. Thus, the analysis resulted in determination of 20 genome-scale coexpression modules in the Gravendeel data set.

Microarrays are notorious for differences in gene expression profiles across data sets (26). To test whether the identified modules could be found in independent patient cohorts, we analyzed 4 additional public expression data sets from distinct institutions (see Supplementary Methods for details on the datasets). We first projected MEs from the Gravendeel data set onto each of the 4 additional patient cohorts and next used the ME projections to determine gene composition of the modules in these data sets (Supplementary Methods). A striking cross–data set similarity of the resulting modules was observed (Fig. 2), with P values ranging from 10⁻¹¹ to 10⁻⁶¹⁹ (see Supplementary Table S2). Furthermore, projections of a given module were always more similar to each other than to projections of the other modules. This confirmed that the modules are common across the data sets generated at distinct medical institutions and represent a consistent feature of the glioma transcriptome organization.

Despite the high overall consistency, we found some genes to be data set–specific. To determine the “core” gene composition of the modules, we compiled consensus modules by including genes that belonged to the same module in at least 3 of the 5 patient cohorts (FDR, false discovery rate, < 2%; Supplementary Methods). This resulted in 20 modules with a refined gene composition that are provided in Supplementary Table S3. Although similar meta-analysis approaches have been used in other cancer types (27, 28), this is, to our knowledge, the first meta-analysis of gene coexpression modules in glioma. The consensus modules represent a collection of high-quality gene expression signatures (Supplementary Table S3) and transcriptionally characterize the key processes inherent to glial tumors (see in the following text).

To determine which processes and structures underlie the observed gene coexpression, we statistically analyzed gene composition of the consensus modules. The modules were found to be associated with a wide range of functions that can be grouped into several categories: processes commonly involved in cancer (e.g., proliferation, hypoxia, angiogenesis, and immune response; ref. 29), glioma differentiation phenotypes (mesenchymal, proneural, and proastrocytic; refs. 1, 8), genetic alterations (e.g., 12q13 amplification and 10q deletion; refs. 3, 30), normal brain functions (infiltrated neurons and myelination), and others (Table 1). A detailed functional characterization of the modules based on their gene composition is provided in the Supplementary Table S4.

Distinct processes, such as angiogenesis and hypoxia, are known to interact in cancer (29). We hypothesized that individual modules may, therefore, be organized into a higher order structure which describes functional relationships between them. To test this, we clustered modules and patients in the Gravendeel data set. The modules were clustered on the basis of Pearson correlations between their MEs (Fig. 3, dendrogram of modules on the left), whereas the patients were clustered on the basis of expression levels of the top 25 hub genes from each module (Fig. 3, dendrogram of patients at the top). The clustering was robust to the choice of the number of hub genes (Supplementary Fig. S1) and reproducible across the data sets (Supplementary Fig. S2). To cross-reference the expression patterns with clinical behavior of the tumors, we displayed clinical characteristics, such as tumor grade, histopathologic diagnosis, and survival below the clustered heatmap (Fig. 3, panels at the bottom). We also analyzed the clinical data quantitatively on the basis of correlations of these parameters with the MEs (Table 2).

Notable in this structured view was the presence of 2 large module coexpression groups (groups A and E) that were expressed in a mutually exclusive way (Fig. 3). These groups consisted of the modules associated with poor and favorable prognoses, respectively (Table 2). Specifically, group A contained the modules associated with cancer progression and was linked to glioblastoma, the most aggressive glioma type (Fig. 3). In contrast, group E included modules related to neurogenesis and cell proneural differentiation and was activated in oligodendroglioma which is less aggressive (Fig. 3). The large scope and strong intensity of this mutually exclusive expression pattern of the 2 module groups underscores biological divergence between glioblastoma and oligodendroglioma. In line with this, a recent study of molecular heterogeneity in brain tumors suggested that gliomas fall into 2 major molecular classes: G type (glioblastoma-like) and O type (oligodendroglioma-like; ref. 12).

The structured view of the transcriptome also illuminates cell differentiation programs in glioma. As shown in the Figure 3, the functionally related modules M19 and M20 associated with proneural differentiation and neurogenesis, respectively, showed a notable distinction in their expression patterns: M20 was upregulated in oligodendroglioma only, whereas M19 was activated both in oligodendroglioma and in the samples containing a large amount of normal neurons (Table 2, Fig. 3). Together with the module M13 which was expressed selectively in the neuron-containing samples, the trio of these modules (M13, M19, M20) thus describes the shared and unique transcriptional features of normal neurons and proneural oligodendroglioma. These results are consistent with the previous reports claiming that although tumor cells in oligodendroglioma tend to morphologically resemble oligodendrocytes (1, 3), they are typically marked by a proneural pattern of gene expression including activation of multiple genes characteristic for mature neurons (M19, e.g., INA, L1CAM, MYT1L, NRXN1, SEPT3; refs. 10, 31). At the same time, oligodendroglioma is clearly distinguished from nonneoplastic brain tissue due to low expression levels of other neuronal genes (M13, e.g., GABRA1, GRIN1, KCNA1, SYN1) and elevated expression of genes characteristic for immature progenitors of brain cells (M20, e.g., ASCL1, BMP2, DLL3, HES6, NLGN1; Fig. 3; ref. 10). Taken together, individual modules are organized into a higher order structure that correlates with clinical characteristics and provides insight into biology of glioma.

Gliomas are notorious for their molecular and clinical heterogeneity (1, 3). To determine molecular subtypes of glioma predictive of its clinical behavior, we examined the patients' clustering in the Gravendeel data set (Fig. 3). Patients with distinct diagnoses (glioblastoma, lower grade astrocytoma, and oligodendroglioma) tended to segregate into separate clusters. Specifically, the glioblastoma-enriched cluster was marked by upregulation of the modules related to tumor mesenchymal differentiation (M8) and proliferation (M11), whereas the clusters enriched in astrocytoma and oligodendroglioma showed an upregulation of the proastrocytic (M15) and neurogenesis (M20) modules, respectively (Fig. 3).

The mesenchymal proliferation and neurogenesis modules are known to define prognostic subtypes of glioma (8, 11). In contrast, the proastrocytic module represents a novel signature, as only a few individual expression markers of glioma proastrocytic differentiation have been reported so far [e.g., APOE (32), ID4 (33), PEA15 (34)]. To evaluate prognostic relevance of the proastrocytic signature in the context of the known glioma subtypes, we reclustered tumors on the basis of hub genes from the mesenchymal (M8), proliferation (M11), neurogenesis (M20), and proastrocytic (M15) modules only. The k-means clustering defined the 4 tumor subclasses (mesenchymal, proliferative, proneural, and proastrocytic) each marked by upregulation of the respective module, as well as the “other” subclass in which all of these modules were expressed at a low level (Supplementary Fig. S3). For each subclass, we calculated a centroid by averaging expression values of the marker genes (Supplementary Table S5) and the centroids were used to classify tumors in the 5 independent data sets (see Methods section). The resulting stratification of samples into classes and the underlying expression patterns are presented in Supplementary Table S6 and Figure S5, respectively.

We compared survival between the obtained patient groups (Fig. 4). In agreement with the previous studies (8, 11), the mesenchymal and proliferative subtypes were associated with poor, whereas the proneural subtype with favorable prognosis (Fig. 4). The “other” subtype was characterized by short survival. Notably, in all the datasets, patients with proastrocytic tumors survived on average several fold longer than those with tumors of the mesenchymal, proliferative, and “other” subtypes (P = 0.005 for the least significant data set, log-rank test; Fig. 4). Survival rates in the proastrocytic and proneural subtypes were highly similar (P > 0.1). Thus, the proastrocytic subtype is strongly and reproducibly associated with favorable prognosis.

We compared our subtypes with those from the earlier studies (6, 8, 11, 12), by calculating proportions of patients shared between them in the respective data sets (Supplementary Methods). Because we observed a substantial overall agreement of our results with those from the previous studies (Supplementary Table S7), we further focused on a detailed comparison of the proastrocytic subtype.

The proastrocytic subtype was found similar to “cluster 17” that was identified by Gravendeel and colleagues as a biologically uncharacterized group of tumors with favorable prognosis (ref. 6; P < 3 × 10⁻⁷, Fisher's exact test). Using GSEA (25), we found that this cluster is characterized by upregulation of known astrocyte-lineage markers (24) compared with the rest of the tumors (P < 0.001, Supplementary Methods). This provides a functional annotation to “cluster 17” and suggests that tumors marked by upregulation of proastrocytic genes have previously been observed as a distinct cluster (6); however, it was not recognized as related to proastrocytic differentiation. These results, thus, for the first time, implicate proastrocytic differentiation as a feature that substantially contributes to expression-based heterogeneity of glioma.

To accurately compare our subtypes with the classification introduced by Phillips and colleagues (8), we classified tumors in each of the data sets into the Mes, Prolif, and PN classes using the respective centroids from the original publication (ref. 8; Supplementary Methods). We found that, in all the data sets, proastrocytic tumors represented strictly a subgroup of the PN class (Supplementary Table S6). In line with this, proastrocytic and proneural genes were jointly upregulated in many tumors which indicates their mixed differentiation phenotype (Supplementary Fig. S4). However, a strong and reproducible negative correlation between expression levels of the proastrocytic and proneural modules in the differentiated gliomas was observed (Supplementary Fig. S5). This inverse relationship between the 2 differentiation programs suggests that the proastrocytic and proneural tumors represent biologically distinct subtypes of glioma (see Supplementary Discussion).

Because proastrocytic and proneural tumors have similar survival rates, addition of the proastrocytic subtype into the classification scheme does not improve the overall strength of association between gene expression and survival. However, as a distinct group of markers, the proastrocytic genes may be important for confirmation of favorable survival prognosis in ambiguous cases during clinical diagnostics. Further research may also identify a distinction in clinical behavior between tumors with proastrocytic and proneural differentiation. Within the proastrocytic module M15 (Supplementary Table S3), genes with the highest mean value of kME (see Supplementary Methods; e.g., DAAM2, NRTK2, AGXT2L1, ALDH2, ADD3, and others) are expected to be the most robust markers of the proastrocytic subtype for potential use in diagnostics.

Taken together, the previously unrecognized proastrocytic module defines a molecular subtype of glioma associated with favorable prognosis. This module complements the known prognostic signatures (“mesenchymal”, “proliferative”, and “proneural”) and offers potential novel markers for diagnostics of the glioma subtypes.

Molecular heterogeneity of gliomas represents a challenge for development of therapies targeted at signaling pathways with aberrant activity in the tumors (2). The EGFR signaling pathway is one of those known to be constitutively activated in a large fraction of highly aggressive gliomas, often via EGFR gene amplification (35). Although this oncogenic pathway is therapeutically promising, molecular mechanisms by which it is regulated are still insufficiently understood (36). We, therefore, analyzed how expression levels of the transcriptomic modules are correlated with the amplification status of the EGFR gene (Supplementary Table S8). The module of protein kinase regulation M3 was found to be, by far, the most significantly and strongly upregulated module upon EGFR amplification (P = 4.6 × 10⁻¹², Wilcoxon test). Furthermore, this module was also the one whose expression profile was most strongly correlated with the expression level of EGFR (Pearson correlation = 0.37, P = 1.3 × 10⁻¹⁰; Supplementary Table S8).

We analyzed gene composition of the protein kinase regulation module (M3) in more detail and found it to contain a number of genes known to play significant roles in glioma such as stem cell marker NES, platelet-derived growth factor PDGFA and a poor prognosis marker FABP7 (Supplementary Table S3). The module also contained genes related to cell adhesion, such as cadherins (CDH2, CDH4, CDH6) and integrins (ITGA7, ITGB8), which suggests its potential upregulation in more invasive tumors (37, 38). Meanwhile, a notable feature of the module M3 was overrepresentation of genes involved in regulation of the mitogen-activated protein kinase signaling cascade (MAPK cascade regulation, P = 1.8 × 10⁻⁴) which is known to act downstream of EGFR (39). Specifically, the module included dual specificity phosphatases (DUSP4, DUSP6), Sprouty genes (SPRY1, SPRY2, SPRY4) and genes of the Sprouty-related family (SPRED1, SPRED2; Supplementary Table S3).

Sprouty proteins were shown to regulate EGFR signaling in normal human cells (40, 41), but it is unknown whether they also play a role in glioma. The observed coexpression of Sprouty and Sprouty-related genes (P = 2.8 × 10⁻⁸) suggests that they are coordinately activated in glial tumors with amplified EGFR and may regulate EGFR signaling in these tumors. Although Sprouty were originally discovered as negative regulators of MAPK signaling (42, 43), SPRY2 has been shown to inhibit ubiquitination and endocytosis of EGFR and consequently enhance EGFR signaling in non–brain human cells (40, 41). This raises the possibility that Sprouty genes that are activated in tumors with amplified EGFR act to sustain constitutively increased EGFR signaling, thus functioning as oncogenes in glioma. While establishing the role of Sprouty in pathogenesis of glioma requires further research, our data indicate that these genes may play an important role in biology of glial tumors and may represent potential targets for antiglioma therapy.

In summary, we explored the structure of glioma transcriptome and found it to be highly similar across multiple tumor sets from distinct medical institutions. The microarray data meta-analysis allowed determination of refined expression signatures for multiple processes and provided a systems-level view of gene expression variation in glioma. The structured view of the transcriptome uncovered a proastrocytic gene expression signature associated with favorable prognosis and revealed that Sprouty proteins which modify EGF signaling likely participate in glioma pathogenesis.

No potential conflicts of interest were disclosed.

We thank Dr. P.J. French for providing us the Gravendeel data set before it became accessible through the GEO database.

This research was supported by a grant from the Centre for Medical Systems Biology within the framework of the Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO) and grant RFBR-10-04-01385-аa from the Russian Foundation for Basic Research. AEI, MGS, Russian Foundation for Basic Research (grant RFBR-10-04-01385-аa); PA 't Hoen, grant from the Centre for Medical Systems Biology within the framework of the Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO).

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