While molecular subgrouping has revolutionized medulloblastoma classification, the extent of heterogeneity within subgroups is unknown. Similarity network fusion (SNF) applied to genome-wide DNA methylation and gene expression data across 763 primary samples identifies very homogeneous clusters of patients, supporting the presence of medulloblastoma subtypes. After integration of somatic copy-number alterations, and clinical features specific to each cluster, we identify 12 different subtypes of medulloblastoma. Integrative analysis using SNF further delineates group 3 from group 4 medulloblastoma, which is not as readily apparent through analyses of individual data types. Two clear subtypes of infants with Sonic Hedgehog medulloblastoma with disparate outcomes and biology are identified. Medulloblastoma subtypes identified through integrative clustering have important implications for stratification of future clinical trials.

While medulloblastoma is widely recognized as comprising four distinct subgroups, the degree of heterogeneity within the four subgroups, and the extent of overlap between the four subgroups is unknown. Applying similarity network fusion to integrate gene expression and DNA methylation profiling, we demonstrate that the degree of overlap between groups 3 and 4 is minimal after accounting for both expression and DNA methylation data. We identify medulloblastoma subtypes within each of the subgroups that have distinct somatic copy-number aberrations, differentially activated pathways, and disparate clinical outcomes. Integrated analysis has refined the boundaries between the four medulloblastoma subgroups, and identified clinically and biologically relevant subtypes, which will inform and improve preclinical modeling, as well as refine our current clinical classification.

Prior reports have recognized the existence of additional substructure within the four consensus subgroups, particularly within groups 3 and 4 (). Consequently, a medulloblastoma consensus conference established that subdivisions within the known subgroups would be defined as subtypes, and labeled α, β, γ, δ, ε, etc. (). In this study our goal was to resolve intra-subgroup heterogeneity and identify biologically distinct and clinically relevant medulloblastoma subtypes by studying a very large cohort of primary tumor samples.

Genome-wide transcriptional arrays and/or genome-wide methylation arrays are the current gold standard for medulloblastoma subgrouping (). These approaches have been used independently with the underlying assumption that they identify similar, perhaps even identical patient clusters. However, the subgroups identified using the two data types in isolation have not been compared head to head. More recently, methods of integrative clustering that analyze multiple data types in aggregate have been developed, including similarity network fusion (SNF) (). Integrative approaches using multiple data types have been suggested to provide superior results compared with the analysis of single data types in isolation. SNF creates a unified view of patients based on multiple heterogeneous data sources, as it can integrate both gene- and non-gene-based data. SNF avoids the bias of genes or features pre-selection, is robust to different types of noise, is highly scalable, and has been shown to outperform other approaches for data integration ().

WNT and SHH medulloblastomas are clearly identifiable and separable across the majority of transcriptional and methylation profiling studies, demonstrating minimal overlap with other subgroups (). Clear heterogeneity exists within the SHH subgroup, which includes infants, children, and adults, although the extent and nature of the substructure is not clearly defined (). The transcriptomes of group 3 and group 4 medulloblastoma are more similar to each other, and several cytogenetic features, such as isochromosome 17q (i17q), are found in both groups (). In response to this, the recent revision of WHO Classification of CNS Tumors has assigned groups 3 and 4 as provisional entities, and a recent consensus on high-risk medulloblastoma left this question unresolved (). Establishing the nature of the boundary between group 3 and group 4 is of clinical importance as outcomes differ, particularly in the setting of upfront metastatic dissemination ().

Genomics has substantially advanced our understanding of medulloblastoma (). While historically considered one entity, it is now clearly accepted that medulloblastoma comprises at least four distinct entities: WNT, SHH, group 3, and group 4; as reflected in the current revision of the WHO classification (). These four subgroups have distinct transcriptional profiles, copy-number aberrations, somatic mutations, and clinical outcomes (). Indeed, current clinical trials and risk stratification biomarkers incorporate the four molecular subgroups (), as do preclinical modeling and the development of novel therapeutics (). However, the extent to which there are additional layers of heterogeneity within the medulloblastoma subgroups is unknown, and a concerted global effort to analyze a very large cohort of tumors will be needed to resolve the question.

Group 3α tumors are enriched for photoreceptor, muscle contraction, and primary cilium-related genes ( Figure 7 B). Pathways involved in protein translation are enriched in groups 3β and 3γ, which are potentially actionable using modulators of protein synthesis such as proteasome inhibitors. Telomere maintenance is also more enriched in group 3γ, suggesting that telomerase inhibition may only be effective in one group. Several pathways are identified across group 4 subtypes, which, coupled with subtype-specific copy-number enrichment, further supports the existence of three group 4 subtypes ( Figure 7 C). Actionable pathways restricted to particular subtypes include MAPK and FGFR1 signaling in group 4β and PI3K-AKT signaling and ERBB4-mediated nuclear signaling in group 4γ. Cell migration pathways are more enriched in group 4α.

Pathway enrichment analysis was performed for each of the identified subtypes across all four subgroups using the top 10% of associated genes across each subtype. We observe several significantly enriched pathways for all identified subtypes (adj. p value < 0.05), supporting subtype-specific biological processes and transcriptional networks ( Figures 7 A–7D ). In particular, in SHH we observe several pathways enriched in SHH β and γ, with developmental pathways more enriched in SHH γ over β ( Figure 7 A). Genes involved in DNA repair and cell cycle are significantly enriched in SHH α. Several actionable pathways, as defined by the availability of approved drugs, are subtype specific. Specifically, sumoylation is enriched in SHH α, ion channels are enriched in SHH β and γ, and telomere maintenance is enriched in SHH α and δ. Receptor tyrosine kinase signaling is enriched in SHH γ and, to a lesser extent, in β. DNA repair pathways are enriched in SHH α, suggesting that strategies to inhibit the DNA damage response and increase replicative stress are more likely to be effective in this group.

(A–D) Enrichment maps representing biological processes and pathways enriched in subtype-specific upregulated genes for SHH subtypes (A), group 3 subtypes (B), group 4 subtypes (C), and WNT subtypes (D). Each node represents a process or pathway; nodes with many shared genes are grouped and labeled by biological theme. Processes and pathways connected at edges have genes in common. Nodes are colored according to the subtype(s) in which the process is enriched; processes enriched in more than one subtype have multiple colors. Nodes sizes are proportional to the number of genes in each process, in each subgroup. Enriched processes were determined with g:Profiler (FDR-corrected q value < 0.05) and visualized with the Enrichment Map app in Cytoscape. Connected nodes and unconnected but actionable nodes are shown.

iCluster was used successfully by TCGA to identify relevant subtypes (). When applying iCluster to our dataset, at k = 4, the four groups did not have the demographics and SCNA consistent with the four previously described groups. When comparing the four iCluster groups with those defined by SNF, WNT and group 3 do not separate, and SHH comprises two groups. When we analyze the iCluster results for five groups, we recover two SHH groups, plus WNT, group 3, and group 4, which in this case corresponds very well to the SNF subgroup (when considering the two SHH groups together). We then asked if we could recover similar subtypes to SNF using iCluster. As we could not recover the four main groups, subgroups defined by SNF were then individually analyzed using iCluster. We observe a near 80% concordance with the SNF subtypes. The childhood and the adult SHH subtypes as well as the group 4 subtypes are recapitulated (along with a single SHH infant group). However, we identified key differences particularly within the WNT, SHH, and group 3 subgroups. Only one WNT group is identified, the two infant SHH subtypes are not identified, and the two distinct group 3 subtypes with MYC amplifications and GFI1 activation are not observed. Clearly, the SNF method is superior at leveraging information of multiple datasets to identify meaningful groups of patients in a cancer cohort, specifically in a medulloblastoma cohort.

Two other integrative clustering methods have been employed by the The Cancer Genome Atlas (TCGA) consortium in previous studies of other cancer histologies. We applied both methods to our dataset; when applying the cluster of clusters (COCA) method used by TCGA in low-grade glioma and pan-cancer studies () we observed that the method was quite limited in the potential to leverage information from our two data types in the current manuscript. The COCA subgroups were driven by the samples that agree or disagree between the two data types clustered in isolation, which is the COCA input. COCA failed to identify one SHH infant subtype or group 3β.

At k = 2, we observe groups 4α and 4γ forming one group, and group 4β being largely preserved ( Figures S8 A, S8B, and S8E). At k = 4, group 4β continues to segregate from the other groups; however, no new groups emerge with any significant clinical or copy-number differences ( Figures S8 A, S8B, and S8F). Due to the enrichment of key SCNA at k = 3, we chose this as our preferred solution. Moreover, our classifier exhibits a decline in confidence at k = 4, suggesting these groups are not as robust as k = 3 ( Table S2 ).

Group 4 is the most prevalent subgroup comprising >40% of all medulloblastomas; previously described features include i17q, tandem duplications of SNCAIP, and high-level amplifications of MYCN and CDK6 (). We observe clear enrichment of key focal and arm-level SCNA at k = 3: group 4α (n = 98), group 4β (n = 109), and group 4γ (n = 119) ( Figures 6 A, S8 A, and S8B ). Clinically we observe group 4β have a slightly higher median age at diagnosis (8.22, 10, and 7 years for groups 4α, 4β, and 4γ; p = 1.34 × 10Pearson's chi-square test, Figure 6 B); however, there is no statistically significant difference in the overall survival ( Figure 6 C) or rate of metastatic dissemination at diagnosis (groups 4α 30/75, 4β 35/86, 4γ 36/94; p = 0.94 Pearson's chi-square test, Figure 6 D). Group 4α are enriched for MYCN amplifications (11/66, compared with none in group 4β and 4γ; p = 2.46 × 10Pearson's chi-square test, Figure S8 C; Table S4 ). Group 4α and 4γ are strongly enriched for 8p loss (group 4α 47/98, 4β 24/109, 4γ 87/119; p = 1.22 × 10Pearson's chi-square test) and 7q gain (group 4α 57/98, 4β 9/109, 4γ 62/119; p = 9.5 × 10, Pearson's chi-square test, Figure 6 E). Group 4β are strongly enriched for SNCAIP duplications (group 4α 4/66, 4β 11/74, 4γ 0/73; p = 0.0019 Pearson's chi-square test) and almost ubiquitous i17q (group 4α 40/98, 4β 87/109, 4γ 31/119; p = 9.75 × 10Pearson's chi-square test) with a paucity of other SCNA ( Figures 6 E and S8 C; Table S4 ). In addition, groups 4α and 4γ are enriched for focal CDK6 amplifications (group 4α 4/66, 4β 0/74, 4γ 6/73; p = 0.051 Pearson's chi-square test, Figure S8 C; Table S4 ). Previous studies have suggested GFI1 and GFI1B activation to be present in group 4, however we see GFI activation to be largely restricted to group 3β ( Figure S8 D).

We find less support for other solutions of group 3, specifically k = 2 and k = 4 ( Figures S6 C, S6D, S7 C, and S7D). At k = 2, we observe a group enriched for MYC amplification (c1 0/38, c2 7/48; p = 0.014 Pearson's chi-square test), and GFI1 family of oncogene activations cluster together (GFI1/1B activation: c1 1/71, c2 29/73; p = 1.14 × 10Pearson's chi-square test) without any meaningful clinical differences ( Figure S7 C). At k = 4, group 3α splits into two groups with minor contributions from the other two groups without any new meaningful clinical or copy-number enrichment ( Figures S6 D and S7 D). In addition the elastic net classifier performs strongly at k = 3 (89%–98.8% per-group accuracy), while at k = 4 one group is less reliably predicted (72% accuracy, Table S2 ).

Three very distinct subtypes of group 3 emerge from our analysis, each with characteristic copy-number and clinical variables: group 3α (n = 67), group 3β (n = 37), and group 3γ (n = 40) ( Figures 5 A, S6 C, and S6D ). A total of 60% of infants under the age of 3 years are in group 3α (age < 3: group 3α 14/63, 3β 4/36, 3γ 5/36; p = 0.021, Kruskal-Wallis test, Figure 5 B). Clinically, groups 3α and 3β have a more favorable prognosis compared with group 3γ ( Figure 5 C). Group 3β are slightly older (p = 0.021, Kruskal-Wallis test, Figure 5 B), and are infrequently metastatic (group 3α 23/53, 3β 5/25, 3γ 15/30; p = 0.058 Pearson's chi-square test, Figure 5 D). Group 3α and 3γ have a similar frequency of metastatic dissemination at diagnosis ( Figure 5 D). Chromosome 8q (MYC locus at 8q24) loss is more frequent in group 3α and gain more frequent in group 3γ (8q gain: group 3α 0/67, 3β 3/37, 3γ 22/40; p = 2.2 × 10Pearson's chi-square test, Figure 5 E), group 3β tumors have a higher frequency of activation of the GFI1 and GFI1B oncogenes, previously shown to be drivers of group 3 through a process termed enhancer hijacking via focal gains and losses on chromosomes 1 and 9, with a paucity of arm-level chromosomal gains and losses (GFI1 or GFI1B activation: group 3α 1/67, 3β 26/37, 3γ 3/40, p < 2.2 × 10Pearson's chi-square test, Figures S7 A and S7B) (). OTX2 amplifications are also enriched in group 3β, as are losses of DDX31 on chromosome 9; previously described to lead to activation of GFI1B through enhancer hijacking (OTX2: group 3α 0/35, 3β 6/28, 3γ 0/24; p = 0.0013; DDX31 deletion: group 3α 1/35, 3β 9/28, 3γ 0/24; p = 0.0031 Pearson's chi-square test, Figure S7 A; Table S4 ). Group 3γ have the worst prognosis (p = 0.036 log rank test, Figure 5 C), a trend to enrichment of i17q (group 3α 17/67, 3β 5/37, 3γ 10/40; p = 0.32 Pearson's chi-square test, Figure 5 E) and frequently harbor increased MYC copy number (group 3α 0/35, 3β 2/28, 3γ 5/24; p = 0.012, Figures 5 F and S7 A; Table S4 ), without other focal aberrations (). Group 3γ have a poor prognosis independent of MYC amplification, expanding the group of high-risk group 3 tumors beyond just MYC status (p = 0.026, log rank test, Figure 5 G).

We identify two WNT subtypes, WNT α (n = 49) and WNT β (n = 21) ( Figures 4 A, S6 A, and S6B); WNT α is comprised mainly of children ( Figure 4 B), has similar survival as WNT β (p = 0.5, log rank test, Figure 4 C), and has ubiquitous monosomy 6 (WNT α 48/49, β 6/21; p = 2.365 × 10Pearson's chi-square test, Figure 4 D). WNT β is enriched for older patients (p = 4.013 × 10, Kruskal-Wallis test, Figure 4 B) who are frequently diploid for chromosome 6 ( Figure 4 D). Monosomy 6 has previously been described as a defining WNT medulloblastoma feature; clearly, patients with WNT β will be misdiagnosed if this criterion is used alone. Prior reports suggesting that adult WNT medulloblastoma might have a different biology and worse prognosis than childhood WNT medulloblastoma, are supported by our current analysis (). At k = 3, we observe a new group, comprised primarily of WNT β without monosomy 6 ( Figures S6 A and S6B); however, in the absence of any other defining feature or clear clinical relevance, we chose k = 2 as our preferred solution.

To interrogate other possible solutions and to present the full results ( Figures S5 A–S5E), we also compared SHH subtypes when divided into three or five SNF groups. We refer to the clusters obtained by SNF for other numbers of groups (k = 3, k = 5 here) as c1, c2, c3, etc. (see Figures S4 A and S4B). When comparing k = 4 with k = 3, SHH α and δ correspond closely to c2 and c1, respectively, with c3 representing a group of infants comprising SHH β and γ ( Figures S4 A, S4B, and S5 A). SHH k = 5 reveals an additional group comprised primarily of a subset of SHH α patients with a group (c3) enriched for 9q loss with a good prognosis and a second group (c5) with a poor prognosis enriched for anaplasia ( Figures S4 B and S5 C–S5E). Several machine-learning classifiers using both data types suggest poor confidence (<80%) in predicting the c5 group. The machine-learning classifier with the best performance, elastic net (), is able to distinguish between four groups with >90% accuracy ( Table S2 ). The identification of two groups of infant medulloblastoma with distinct clinical behavior allows for more precise and rational planning of clinical trials for infants with SHH medulloblastoma ().

Interestingly, infant SHH tumors are mainly distributed across SHH β and SHH γ (age < 3: SHH α 5/65, β 23/35, γ 34/47, δ 0/76; p = 2.2 × 10Pearson's chi-square test, Figure 3 B), with disparate outcomes and copy-number profiles. SHH β tumors are frequently metastatic (33.3% versus 9.4% in SHH β and γ; p = 0.027 Pearson's chi-square test, Figure 3 G), harbor focal PTEN deletions (25% in SHH β versus none in γ), have multiple focal amplifications ( Figure S4 C; Table S4 ), and have a worse overall survival compared with SHH γ (HR of SHH β versus γ: 2.956 95% CI: 0.908–9.63; p = 0.059 Cox proportional hazards, Figure 3 C). The difference in outcomes between SHH β and γ is possibly related to the increased rate of metastatic dissemination in SHH β, as there is a clear trend toward metastases being a marker of poor outcome within SHH β (HR of SHH β metastatic versus non-metastatic: 3.621 95% CI: 0.798–16.44; p = 0.096 Cox proportional hazards). Conversely, SHH γ have a relatively quiet copy-number landscape, with no recurrent amplifications, only one low-level recurrent focal deletion, and no significant arm-level gains ( Figures 3 D and S4 C). Moreover, SHH γ are enriched for the MBEN (medulloblastoma with extensive nodularity) histology (20.9%; p = 2.34 × 10, Pearson's chi-square test, Figure 3 H), which is known to portend more indolent clinical behavior (). Although almost all SHH tumors with MBEN histology (n = 10) were assigned to it, only a minority of SHH γ tumors have MBEN histology, demonstrating that histology alone is an inadequate surrogate to identify SHH γ tumors. The survival difference of SHH γ patients is not statistically significant between MBEN and non-MBEN tumors, suggesting that subtype affiliation is a more powerful biomarker than histopathology in infants with SHH medulloblastoma (p = 0.268, log rank test, Figure 3 I). SHH δ are primarily composed of adults, have a favorable prognosis, and are strongly enriched for TERT promoter mutations (SHH α 6/34, β 2/22, γ 7/26, δ 38/42; p = 8.13 × 10, Pearson's chi-square test, Figure 3 J).

Applying SNF and spectral clustering on SHH subgroup samples at k = 4 identified four clinically and cytogenetically distinct groups: SHH α (n = 65), SHH β (n = 35), SHH γ (n = 47), and SHH δ (n = 76) ( Figures 3 A , S4 A, and S4B). SHH α tumors primarily affect children aged 3–16 years ( Figure 3 B), have the worst prognosis (p = 0.03, log rank test, Figure 3 C), and are enriched for MYCN amplifications (SHH α 8/37, β 3/23, γ 0/29, δ 1/48; p = 0.0034 Pearson's chi-square test), and GLI2 amplifications (SHH α 6/37, β 0/23, γ 0/29, δ 0/48; p = 0.0002 Pearson's chi-square test, Figure S4 C; Table S4 ). Specific CNAs including 9q loss (SHH α 42/65, β 8/35, γ 11/47, δ 17/76; p = 2.94 × 10Pearson's chi-square test), 10q loss (SHH α 29/65, β 6/35, γ 7/47, δ 6/76; p = 1.54 × 10Pearson's chi-square test), 17p loss (SHH α 24/65, β 5/35, γ 3/47, δ 8/76; p = 3.44 × 10Pearson's chi-square test, Figure 3 D), and YAP1 amplifications (SHH α 3/37, β 0/23, γ 0/29, δ 0/48; p = 0.04 Pearson's chi-square test, Figure S4 C; Table S4 ) are also enriched in SHH α. The recent WHO classification includes SHH-activated TP53 mutant tumors as a distinct category based on studies showing this group as being very high risk (). To further explore this association, TP53 was sequenced across 145 SHH samples. TP53 mutations are highly enriched in SHH α (SHH α 14/40, β 2/27, γ 2/31, δ 6/47; p = 0.0026 Pearson's chi-square test, Figure 3 E; Table S5 ). When survival is analyzed stratified by TP53 mutation and SHH α subtype, TP53 mutations are only prognostic in SHH α (HR TP53 mut versus WT: SHH α 6.006 [95% CI: 1.586–22.75; p = 0.00832] and non-SHH α 1.222 [95% CI: 0.2795–5.342; p = 0.79, Cox proportional hazards, Figure 3 F]).

To determine whether analysis of a single data type in isolation yielded similar results, we performed NMF clustering using gene expression or DNA methylation data individually. Using NMF clustering of the most variable expressed genes and methylated probes, we found that the two different types of data yield discordant subtypes as defined by both the cophenetic coefficient and silhouette value (>0.9) criteria ( Figures S3 A–S3D). In addition, the group memberships between the two modalities are divergent, indicating a lack of agreement between expression and methylation when analyzed in isolation ( Figures S3 E–S3H). When compared with the SNF subtypes, we found important differences, suggesting that both methylation and expression signatures contribute significantly and differently to define heterogeneity within the four subgroups; the data types provide distinct but complementary signals that improve over single-modality analyses. The subtypes identified by SNF are truly a combination of information present in both datasets, and therefore both data types are required to gauge the true intertumoral heterogeneity of medulloblastoma. For example, we observe that SHH α is mainly supported by the methylation data, but the defined group does not contain all SHH α samples (61%, Figure S3 B). SHH δ is strongly supported by both the expression and methylation data ( Figure S3 B). In addition, groups 3β and 3γ are mainly defined by the signatures found in the expression data and do not separate well using the methylation data alone ( Figure S3 C). Finally, group 4γ is very well supported by the methylation data, and corresponds to a group obtained with the expression data, but this latter group is missing 24.4% of group 4γ samples ( Figure S3 D). Group 4β is well supported by both data types ( Figure S3 D). We conclude that methylation and expression data are complimentary, and an integrated approach allows a unified view of the underlying groups that is very valuable in elucidating heterogeneity within subgroups.

For each subgroup, we identified the top associated genes and methylation probes that best support the final subtypes. Analysis of the top 1% of the associated genes and methylation probes for each subgroup demonstrates that the subgroups are supported by specific gene sets and methylation probes that vary substantially across subtypes ( Figures 2 A, 2B , and S2 A–S2D; Table S3 ). We evaluated the relationship between the associated genes and methylation probes in each subgroup. We first evaluated the number of associated genes that had methylation probes in their promoter region. Then we identified the subset of associated genes for which those probes were subgroup associated, and finally checked if we could detect an anti-correlation between the associated gene expression and the associated probe methylation levels. Only 3.7%, 8.3%, 6%, and 13% of WNT, SHH, group 3, and group 4 associated genes, respectively, follow all the criteria described above ( Figure 2 C). Therefore, only a small percentage of the associated genes are directly affected by DNA methylation. This is in support of both DNA methylation and gene expression contributing to the heterogeneity observed within each subgroup.

(C) Percentage of genes associated for each subgroup; (1) that have methylation probes in their promoter region, (2) for which those methylation probes are in the top 1% associated probes of the respective subgroup, and (3) for which an anti-correlation can be detected between the gene expression and methylation probes levels inside the subgroup. The numbers of genes in each category are indicated.

(A and B) Heatmap of the top 1% most associated genes (A) and the top 1% most associated methylation probes (B) for the subtypes inside each subgroup (left side color bar), respectively. Top color bars indicate the subgroup and subtype sample affiliation. Samples are ordered by subtype.

We applied SNF and spectral clustering within each of the four subgroups as defined by k = 4 across the entire cohort to determine the extent and nature of intra-subgroup heterogeneity. SNF and spectral clustering were selected to reduce the noise introduced by biased feature selection, and to leverage the full spectrum of our dataset. We identified clusters from k = 2 to k = 8 within each subgroup. In addition, we applied seven different machine-learning classifiers to predict the SNF subtypes. Cluster assignments from spectral clustering on the SNF fused similarity matrix was used as the “ground truth” subtype assignments. We split the dataset into a 70% training set and 30% testing set, trained the various classification models in 5-fold cross-validation on the training set and repeated the procedure 100 times ( Table S2 ). We then applied the following criteria a priori to select the optimal number of subtypes: (1) how similar are the SNF clusters on the sample-to-sample heatmap? (2) How subtype specific are the broad and focal SCNA? (3) How relevant are the clinical associations? (4) How robustly can these subtypes be predicted using supervised machine learning? Using these criteria, we identified 12 subtypes: two WNT, four SHH, three group 3, and three group 4. For each solution, we identified focal SCNA from SNP6 data and arm-level copy-number gains and losses using copy-number states inferred from the methylation arrays.

To determine the degree of overlap between groups 3 and 4, we undertook unsupervised clustering of 470 group 3 and 4 tumors using DNA methylation array data only, and then subsequently using transcriptional profiling data only. Both k-means and non-negative matrix factorization (NMF) consensus clustering revealed a small subset of tumors (2.9%–8.9%) that switched subgroup between k = 2 and k = 3 as determined through analysis of either transcriptional or methylation data ( Figures S1 G and S1H). Strikingly, the set of “ambiguous group 3–4 tumors” identified by gene expression profiling had very little overlap with those identified by DNA methylation profiling ( Figure 1 D) suggesting that the identification of the ambiguity may be a limitation of the particular type of measurement or data, rather than the identification of a truly distinct biological subtype. Examination of tumors within the “overlap” group does not reveal any demographic, clinical, or genetic commonalities, suggesting that it could be an artifact rather than a biologically discrete, clinically important group. Subsequent application of SNF and spectral clustering to this cohort of group 3 and 4 samples demonstrates that only 13/470 (2.8%) of samples change subgroup between k = 2 to k = 3, and of these 13 only 3 (0.64%) do not track back to their original subgroup when k > 3 ( Figures 1 E and S1 I). We conclude that group 3 and group 4 medulloblastomas are stable, mostly non-overlapping molecular subgroups, and that SNF followed by spectral clustering is a more robust method of delineating subgroups than using a single data type in isolation.

Groups 3 and 4 are more similar to each other than to SHH and WNT ( Figures 1 B and 1C). We tested the stability of these core subgroups, by counting samples that switch subgroup affiliation when the number of clusters increases ( Figure 1 A). Following each sample from k = 4 to k = 12, no sample changed affiliation between WNT and SHH, while a small minority of samples moved between groups 3 and 4.

To these samples, we applied SNF to integrate both gene expression and DNA methylation data, followed by spectral clustering ranging from 2 to 12 groups. At k = 4, four distinct subgroups are clearly identified. Those groups correspond clinically and structurally to the previously described consensus subgroups: WNT (n = 70), SHH (n = 223), group 3 (n = 144), and group 4 (n = 326) ( Figures 1 A–1C and S1 A–S1F) ().

(E) Tumor clusters obtained through spectral clustering on the SNF network fused data of group 3 and 4 samples (n = 470). A small minority of samples (n = 13, 2.8%) that were initially classified as group 3 samples at k = 2, subsequently move to group 4 at k = 3. Only 3/470 (0.64%) samples remain in group 4 after k = 5. These samples are tracked up to k = 8 (orange).

(B) Network representation of the relationships between tumors (k = 4). The shorter the edge between samples (nodes) is the more similar the samples are (only edges with a similarity value above the median value of all patient to patient similarity values are displayed).

Through the Medulloblastoma Advanced Genomics International Consortium, we assembled a cohort of 763 primary frozen medulloblastoma samples with high-quality DNA and RNA, and generated genome-wide methylation and expression profiles. Of these, 491 had DNA copy-number profiles generated by Affymetrix SNP6 microarrays (). Clinical data including age, tumor histology, metastatic status, and survival were available on 95.7%, 76.9%, 75.2%, and 82% of cases, respectively ( Table S1 ). Arm-level somatic copy-number aberrations (SCNA) were inferred from methylation arrays in 100% of cases.

Discussion

Figure 8 Graphical Summary of the 12 Medulloblastoma Subtypes Show full caption Schematic representation of key clinical data, copy-number events, and relationship between the subtypes inside each of the four medulloblastoma subgroups. The percentages of patients presenting with metastases and the 5-year survival percentages are presented. The age groups are: infant 0–3 years, child >3–10 years, adolescent >10–17 years, and adult >17 years. Our study identifies and delineates the intertumoral heterogeneity present within medulloblastoma subgroups. Leveraging a large cohort of medulloblastomas profiled by combined gene expression and DNA methylation, we have identified different subtypes within each of the four core subgroups. These subtypes have particular clinical and copy-number features, which allow for a refinement in our understanding of the genomic landscape of medulloblastoma ( Figure 8 ). Combining expression and methylation data using SNF adds further proof that groups 3 and 4 are largely different biological entities. The deeper we go in clustering medulloblastoma samples, the less consistent the groups become. This is exemplified by poor predictability of putative subtypes when a large number of subtypes is assumed. Defining clinical features and CNAs also tend to lose their distinctive profiles as we increase the number of clusters, suggesting that heterogeneity is bounded by a discrete number of optimal groups.

Comparison of SNF with consensus clustering of either gene expression or DNA methylation data analyzed in isolation clearly suggests that an integrated approach provides a much more refined and accurate classification. This is particularly striking when evaluating the boundary between groups 3 and 4, where samples that are deemed indeterminate using gene expression and DNA methylation in isolation are largely non-overlapping. Moreover, in elucidating the heterogeneity within subgroups, we observe significant disagreement between gene expression and DNA methylation in isolation, suggesting that each data type makes a unique and non-redundant contribution to defining the subtypes. The very low number of samples that change subgroup affiliation using SNF strongly advocates that definition of these two groups is largely enhanced using an integrative approach. A limitation of our approach is the bulk analysis of samples. At a subclonal level, a greater degree of overlap across groups 3 and 4 cannot be discounted. More detailed analysis at a cellular level, specifically applying single-cell methods, will help delineate the full subclonal structure, potentially uncovering subsets of group 3 and 4 samples with common mechanisms and cellular origins. Further studies integrating emerging technologies such as long non-coding RNA, proteomics, and histone modifications may allow an even more refined description of the medulloblastoma landscape; however, the large cohorts of frozen tissue required for these studies are presently not available.

Northcott et al., 2012b Northcott P.A.

Shih D.J.

Peacock J.

Garzia L.

Morrissy A.S.

Zichner T.

Stutz A.M.

Korshunov A.

Reimand J.

Schumacher S.E.

et al. Subgroup-specific structural variation across 1,000 medulloblastoma genomes. The identification of subtypes has significant biological and clinical implications. Several previously described copy-number alterations within medulloblastoma subgroups such as amplifications/gains of MYC, MYCN, OTX2, CDK6, SNCAIP, and ACVR1, as well as several arm-level events including i17q clearly segregate between subtypes (). Our identification of unique cytogenetic aberrations that occur in concert, as well as specific biological pathways enriched within specific subtypes, will serve to inform creation of rational preclinical models that closely mirror the human diseases. Several of these aberrations are actionable and largely restricted to subtypes, which will also allow for a more personalized treatment approach. Several subtypes, particularly in SHH and group 3, have clear and drastic clinical and prognostic differences, which will allow for more robust risk stratification in future clinical trials. Furthermore, a major hurdle to clinical trial design has been the overlap of groups 3 and 4 in current studies, which if applied today would make strata assignment difficult. The next generation of clinical trials for high-risk medulloblastoma will involve subgroup-specific therapies. The inability to stratify 10% of patients to either groups 3 or 4 has the potential to either deprive a patient of an innovative therapy or, of more concern, expose a child to an inappropriate escalation or de-escalation of therapy.

Ramaswamy et al., 2016a Ramaswamy V.

Remke M.

Bouffet E.

Bailey S.

Clifford S.C.

Doz F.

Kool M.

Dufour C.

Vassal G.

Milde T.

et al. Risk stratification of childhood medulloblastoma in the molecular era: the current consensus. Clinically, our observed groups have immediate implications. It has been shown that TP53 mutations are highly prognostic in SHH. We extend these findings whereby TP53 mutations are not only enriched in SHH α but also only prognostic in SHH α. This is highly relevant for clinical trial design, where TP53 mutant SHH has been identified as a very-high-risk group to be prioritized for novel therapies in both Europe and North America (); clearly, the observation that TP53 mutations are highly enriched and prognostic in SHH α has significant implications. A limitation of this is the absence of germline status, which, based on previous studies, are likely TP53 mutant enriched in SHH α.

The identification of two infant SHH groups has clear and immediate clinical significance. Currently, infant medulloblastomas are stratified by the presence or absence of desmoplastic morphology. However, several reports have suggested that infant SHH as a whole have a favorable prognosis independent of morphology. Our results suggest that clinical risk stratification can be refined by incorporation of integrated subtypes, whereby SHH γ are clearly a very low-risk group and could be spared the toxic effects of high-dose chemotherapy. Our observation that MBEN histology is almost exclusive to SHH γ, but represent a minority of cases within SHH γ, has significant implications for clinical trials. Current infant clinical trials stratify patients based on either classic or desmoplastic/MBEN histology. Indeed, the frequency of desmoplastic histology is similar across all four SHH subtypes, despite significant differences in survival between SHH subtypes. The most recent infant medulloblastoma study from the Children's Oncology Group ACNS1221 (NCT02017964) was closed prematurely due to an excess of relapses. This study selected infants with a “desmoplastic” morphology for treatment de-escalation, of which the vast majority are SHH. Indeed, our identification of two infant subtypes of SHH represents an example where more robust risk stratification has the potential to accurately select patients for de-escalation of therapy in future clinical trials. Overall, this further supports the idea that the incorporation of molecular stratification rather than subjective morphology alone has the potential for immediate clinical benefit.

Similarly, for group 3, we identify a high-risk group that is enriched for MYC amplification, but for which not all patients are MYC amplified. Interestingly, the majority of in vitro cell lines of medulloblastoma do not represent the clear intertumoral heterogeneity, but rather are MYC-amplified or MYC-activated models that actually represent only group 3γ. The identification of significant heterogeneity across group 3 underlies the urgent need to develop preclinical models that faithfully recapitulate the different subtypes within each subgroup. In group 4, there are currently no robust preclinical models, and the subgroups we describe, specifically the mutually exclusivity of MYCN amplifications and SNCAIP duplications, may help with future modeling.

Taken together, our results highlight the power of combining multiple data types compared with the use of single data types in isolation. This approach has identified that there may be a limit to the degree of substructure across medulloblastoma samples; however, only a study with a much larger cohort could fully assess the extent of intertumoral heterogeneity within the subgroups. We identify clinically important substructure within subgroups, which will allow further refinement of our biological and clinical risk stratification schemes. The identification of homogeneous subtypes may simplify the identification of targets for therapy, and could allow for therapies effective across subtypes. The development of reliable biomarkers to identify subtypes will provide much needed prognostic information for patient stratification, particularly in regard to SHH and group 3 medulloblastoma.