{'fraction_non_alphanumeric': 0.03546373529249296, 'fraction_numerical': 0.020659559189084683, 'mean_word_length': 5.297964863454514, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Resolving medulloblastoma cellular architecture by single-cell genomics

Summary Medulloblastoma is a malignant childhood cerebellar tumour comprised of distinct molecular subgroups. Whereas genomic characteristics of these subgroups are well defined, the extent to which cellular diversity underlies their divergent biology and clinical behaviour remains largely unexplored. We used single-cell transcriptomics to investigate intra- and inter-tumoural heterogeneity in twenty-five medulloblastomas spanning all molecular subgroups. WNT, SHH, and Group 3 tumours comprised subgroup-specific undifferentiated and differentiated neuronal-like malignant populations, whereas Group 4 tumours were exclusively comprised of differentiated neuronal-like neoplastic cells. SHH tumours closely resembled granule neurons of varying differentiation states that correlated with patient age. Group 3 and Group 4 tumours exhibited a developmental trajectory from primitive progenitor-like to more mature neuronal-like cells, whose relative proportions distinguished these subgroups. Cross-species transcriptomics defined distinct glutamatergic populations as putative cells-of-origin for SHH and Group 4 subtypes. Collectively, these data provide novel insights into the cellular and developmental states underlying subtype-specific medulloblastoma biology.


Introduction
Medulloblastoma (MB) comprises a series of molecularly and clinically diverse malignant childhood cerebellar tumours 1 . While advances in treatment have improved survival, many patients suffer from neurological sequelae or still succumb to their disease. Genomic studies of bulk patient cohorts have defined four consensus molecular subgroups (WNT, SHH, Group 3, and Group 4) 2 , each characterized by discrete genomic landscapes, patient demographics, and clinical phenotypes [3][4][5][6][7] . The association between genotypes, transcriptional profiles, and patient age at diagnosis suggests that distinct MB subgroups arise from the transformation of different cell types in precise spatio-temporal patterns. Such genotype to cell type associations have been partially investigated for WNT-and SHH-MBs, which are thought to originate from cells in the extracerebellar lower rhombic lip 8 and from cerebellar granule neuron progenitors (GNPs) 9,10 , respectively. In contrast, cellular origins of Group 3-and Group 4-MB remain unconfirmed. Overlapping transcriptional and epigenetic signatures observed in bulk profiling studies have consistently hampered definitive classification of Group 3 and Group 4 tumours and suggest they may share common developmental origins 3,11 . Thus, a better understanding of MB cellular composition and substructure according to subgroup is a critical goal, especially for poorly characterized Group 3 and Group 4.
Single-cell RNA-seq (scRNA-seq) has emerged as a powerful method to comprehensively characterize cellular states within healthy and diseased tissues 12 . While in central nervous system malignancies scRNA-seq has been applied to decipher adult and pediatric gliomas [13][14][15][16] , such approaches have yet to be deployed across MB subgroups. Here, we applied full-length scRNA-seq across consensus MB subgroups to infer cellular trajectories, deconvolute bulk MB expression cohorts, and nominate developmental origins. We find that WNT, SHH, and Group 3 tumours exhibit subgroup-specific cellular trajectories that consist of malignant undifferentiated and differentiated neuronal-like populations, whereas Group 4 tumours recapitulate more differentiated populations of known lineage. Collectively, these data provide new insights into the molecular and cellular architecture of MB across all subgroups, with the potential to inform future studies aimed at improving patient outcomes.
To classify single cells into malignant and non-malignant subsets, we used two complementary strategies. First, we inferred genome-wide CNVs from the scRNA-seq data as previously described 13 (see Methods). This analysis identified large-scale genomic gains and losses in most patient samples (21/25), including hallmark alterations such as monosomy 6 (WNT) and isochromosome 17q (Groups 3 and 4; ED Figure 2a-e). Few cells (n=36) from these patients lacked discernable CNVs (see Methods). Second, we clustered single cells across all samples according to their transcriptional profiles. A minority of single cells in our cohort clustered with reference immune cells (n=6) or oligodendrocytes (n=22). (ED Figure 3a,b). All cells that lacked CNVs and/or clustered with normal reference populations were deemed as non-malignant and excluded from further analysis (n=43). Across individual tumours, 96-100% of cells were classified as malignant, consistent with prior estimates of high MB tumour cell fractions based on genome sequencing 20 . We further validated these assignments by quantifying genetic mutations identified by bulk tumour DNA sequencing in our scRNA-seq data (1,937 mutant and 1,952 wild-type transcripts detected; see Methods; ED Figure 3c-f).
To relate MB single-cell profiles to normal developmental hierarchies, we leveraged recently generated scRNA-seq data for murine cerebellar development spanning 13 21 . Canonical correlation analysis (CCA; see Methods) facilitated cross-species comparisons between our murine cerebellar single-cell, human MB single-cell, and bulk 3 expression datasets. SHH-MB was highly correlated with GNP populations (cosine distance=0.54), consistent with literature 9,10 supporting GNPs as the cell-of-origin for this subgroup ( Figure  1d; ED Figure 4f). Notably, Group 4-MB was highly correlated with unipolar brush cells (UBC; cosine distance=0.50) and glutamatergic cerebellar nuclei (GluCN; cosine distance=0.49). In contrast, we did not detect high-confidence correlations between any cerebellar populations and either WNT or Group 3 subgroups.

Malignant trajectories within WNT-MB
WNT subgroup patients account for ~10% of childhood MBs and have an excellent prognosis 22 . Somatic CTNNB1 mutations or germline APC mutations, both of which drive constitutive WNT signaling, are found in virtually all WNT-MBs 3,23 . Five WNT tumours were included in our dataset. Pairwise correlation analysis revealed multiple distinct transcriptional states that were consistently identified within these tumours ( Figure 2a). Inferring CNVs from our scRNA-seq data identified four cases with monosomy 6 -a stereotypic genomic feature of this subgroup (Figure 2b; ED Figure 2a). The fifth case (BCH807) exhibited chromosome 19 gain and was negative for nuclear beta-catenin by immunohistochemistry (IHC; data not shown), both of which are atypical characteristics for this subgroup despite high-confidence molecular classification as WNT-MB (Figure 1b). SJ99 exhibited striking heterogeneity at both a transcriptional and genetic level, with evidence for two distinct subclones. Subclone SJ99-A exhibited monosomy 6 and chromosome 17p loss, whereas subclone SJ99-B exhibited broad gains and losses affecting nearly every chromosome. Investigation of genetically supported SNVs confirmed expression of mutant transcripts in 57.2% of cells (including key WNT-MB driver genes CTNNB1, DDX3X, and TP53; Figure 2c). Non-negative matrix factorization (NMF) was applied to define underlying transcriptional programs specific to each tumour (ED Figure 5a,b; Supplementary Table 2a; see Methods). This analysis uncovered remarkably similar programs in all five WNT-MBs, which we grouped accordingly into four meta-programs (WNT-A/B/C/D). To interpret the characteristics of each meta-program, we evaluated their underlying gene signatures. WNT-A contained numerous markers of cell cycle activity (e.g. TOP2A, CDK1, and RRM2; P<0.001, Fisher's exact test; Supplementary Table 2b). WNT-C was characterized by markers of neurogenesis/neuronal differentiation (STMN2, KIF5C, SYT11; P<0.001; Figure  2d). WNT-B consisted of ribosomal and metabolic genes (NME2, HK2, PGM5), while WNT-D contained select WNT pathway genes (LRP4, APCDD1) and immediate early response genes (JUNB, EGR1; Figure 2d). Cells scoring highest for both WNT-B and WNT-D expressed elevated levels of additional canonical WNT pathway genes (DKK2, AXIN2, WIF1) and MYC (ED Figure 5c; Supplementary Table 2c). We thus interpreted that these meta-programs reflect cell cycle activity (WNT-A), neuronal-like differentiation (WNT-C), and two WNT-driven states (WNT-B and WNT-D), with WNT-B characterized by elevated protein biosynthesis and metabolism ( Figure 2d). RNA in situ hybridization (ISH) performed on the same tumours validated the expression of specific meta-program marker genes in subpopulations of cells (ED Figure 5d). Moreover, scoring each cell in our cohort for these four meta-programs defined a putative developmental trajectory for WNT-MB, with cell cycle activity restricted to cells that were both WNT-B high and WNT-C/D low (Fisher's exact test, P<0.001), suggesting that this subpopulation fuels WNT-MB growth ( Figure 2e). Importantly, each meta-program was identified in at least four samples (ED Figure 5a), suggesting that they reflect shared features of WNT-MB.

Developmental trajectories within SHH-MB
As the dominant subgroup in both infants (≤3yrs) and adults (≥18yrs) 22 Figure 3b). SHH-C was defined by markers of neuronal differentiation (STMN2, MAP1B, TUBB2B, SEMA6A; P<0.001; Figure 3b). We interpret that these programs reflect cell cycle activity (SHH-A), undifferentiated progenitors (SHH-B), and more differentiated neuronal-like programs (SHH-C). Scoring each SHH-MB cell for these programs defined a putative developmental trajectory, with proliferating cells restricted to undifferentiated progenitors ( Figure 3c). These respective programs were partially recapitulated in SHH subgroup PDX models (ED Figure 6c,d).
To investigate the developmental significance of these findings, we used CCA to compare SHH-MB meta-programs to murine cerebellar populations. SHH-B correlated with undifferentiated UBC/GNP and GNP populations, while SHH-C correlated with UBC/GN

Malignant programs within Group 3/4-MB
Group 3 and Group 4 tumours account for ~60% of MB diagnoses and remain the least understood with respect to disease biology and developmental origins 7 . Group 3 tumours are frequently metastatic at diagnosis and typified by genomic amplification/over-expression of MYC, which is associated with unfavorable outcomes 11,28 . Group 4 tumours are metastatic at diagnosis in approximately one third of patients and harbor recurrent chromatin modifier alterations 28,29 . Recent bulk profiling studies have demonstrated marked molecular and clinical heterogeneity in Group 3/4, with a subset of tumours exhibiting overlapping molecular signatures that confound robust classification 3,30,31 . We next investigated whether recently described 3,32 DNA methylation-based subtypes of Group 3/4-MB were related to the meta-programs inferred from scRNA-seq. We found that DNA methylation subtypes I and V, both of which contain a mixture of Group 3-and Group 4-MBs, were significantly enriched for tumours with intermediate expression patterns (P<0.001, Fisher's exact test; Figure 5c; ED Figure 9c). These results suggest that a continuum of cellular states accounts for the molecular substructure seen in Group 3/4 that complicates accurate consensus classification.

Lineage-specific correlates of Group 4-MB
We next sought to compare and inter-relate the different subgroup-specific meta-programs. To this effect, we applied all observed meta-programs (n=10) to all 7,745 malignant cells in our dataset. Pairwise correlation of expression scores confirmed high similarity among cell cycle programs (WNT-A, SHH-A, Group 3/4-A; average r=0.99) (Figure 6a). The undifferentiated progenitor-like programs (WNT-B, SHH-B, Group 3/4-B) exhibited low correlations (average r=0.23), in agreement with their distinct underlying biology. In contrast, the neuronal-like differentiation programs (WNT-C, SHH-C, and Group 3/4-C) were highly correlated (average r=0.77; Figure 6a,b; ED Figure 9d), consistent with shared capacity for neuronal differentiation across subgroups. We reasoned that the neuronal-like differentiation programs defined in each subgroup consist of general neuronal differentiation markers, potentially masking markers of specific lineages. To elucidate markers that might inform developmental origins, we compared genes specific to neuronal-like cells in the different subgroups (n=260; relative to undifferentiated cell populations; see Methods). Half of these genes (52%) were shared between at least two subgroups and included general markers of neuronal differentiation (e.g. ENO2, SYT11, TUBB3, MAP2), while the remainder were exclusive to individual subgroups (13-20%; Figure

Discussion
Despite extensive characterization of MB genomic landscapes, effective subgroup-specific therapies have yet to emerge, suggesting that a deeper understanding of the biological and cellular basis of MB is essential. This is particularly urgent for Group 3-and Group 4-MB, which often bear inferior outcomes. As a first challenge, these subgroups have proven difficult to accurately classify, confounded by transcriptional and epigenetic ambiguity. Our combined single-cell analysis of Group 3/4-MBs confirmed that prototypic Group 3-MBs are dominated by undifferentiated progenitor-like cells, whereas prototypic Group 4-MBs are comprised almost exclusively of more differentiated neuronal-like cells. Importantly, we identified a subset of 'intermediate' tumours characterized by varying proportions of both undifferentiated and more differentiated populations (ED Figure 10g). These findings offer a novel molecular and cellular explanation for the challenges associated with Group 3 and Group 4 sub-classification and provide a framework for future classifications that incorporate population heterogeneity.
Cellular origins for WNT-and SHH-MB have been mostly informed from genetically faithful mouse models 8,9 . Cross-species transcriptional analyses performed here confirmed significant correlations between SHH-MB and GNPs of variable differentiation states that were associated with patient age. Moreover, our analyses identified UBCs and GluCN as cellular correlates of Group 4-MB subtypes, building on prior studies that have implicated glutamatergic cellular origins for Group 4 4 . For WNT-MB, we failed to identify significant correlation between malignant single-cell programs and cerebellar populations, compatible with an extracerebellar origin for this subgroup 8 . No significant correlations were detected between Group 3-MB and our cerebellar dataset. This observation may be attributed to transformation and cellular reprogramming driven by specific oncogenes (i.e. MYC) or may imply that Group 3-MBs have an extracerebellar origin. It is also plausible that our murine reference atlas was incomplete and lacked populations pertinent to either WNT-or Group 3-MB origins. Technical limitations of comparing single-cell datasets between species should not be underestimated, warranting future studies of the cellular correlates between human cerebellar and MB single cells.
In conclusion, our work provides the first cellular atlas of MB across all subgroups and a cross-species comparison to cerebellar development, shedding light on putative subgroupspecific origins. Our analyses also define the cellular states underlying each MB subgroup, disentangling determinants of intra-and inter-tumoural heterogeneity. These findings will enable future studies to assess translational opportunities and to evaluate the impact of therapeutic approaches on the spectrum of cellular states that drive MB.

Tissue sample collection and dissection
Human primary tumours-Patients and their parents at Boston Children's Hospital, the Medical University of Vienna and St. Jude Children's Research Hospital gave consent preoperatively according to Institutional Review Boards. Fresh tumours were collected at the time of surgery and processed directly. Tumour samples from Boston Children's Hospital and the Medical University of Vienna were mechanically and enzymatically dissociated using a papain-based brain tumour dissociation kit (Miltenyi Biotec). Tumour samples from St. Jude Children's Research Hospital were precut and dissociated for 30 minutes at 37°C in papain solution (10 units/ml, Worthington, LS003126) containing N-acetyl-L-cysteine (160 ug/ml, Sigma-Aldrich, A9165) and DNase I (12 ug/ml, Sigma-Aldrich, DN25), rinsed in Neurobasal medium and filtered using a 40 μm strainer.
Murine cerebellum-Murine cerebellar tissue from Crl:CD1 (ICR) mice at 13 distinct developmental time points was previously isolated 21 . Embryonic time points include each day between e10 to e18 and postnatal time points include P0, P4, P7, and P10. Two biological replicates were included at each time point and three at e14. Cerebella were dissociated as previously described.

Generation and processing of DNA methylation data
All single-cell patient and PDX samples were analysed using Illumina Infinium Methylation EPIC BeadChip arrays according to the manufacturer's instructions. Data was generated from both freshly frozen and formalin-fixed paraffin-embedded (FFPE) tissue samples. Medulloblastoma subgroup predictions were obtained from a web-platform for DNA methylation-based classification of central nervous system tumours (www.molecularneuropathology.org, version 11b4 35 ). Resulting assignment of samples to WNT, SHH, Group 3, and Group 4 subgroups were used for all downstream analyses. A similar classification system was used for predicting medulloblastoma subtypes 32 . Copynumber variation (CNV) analysis from EPIC methylation array data was performed using the conumee Bioconductor package. Identified CNVs were compared to those predicted from the single-cell data (shown in EDF2).

Generation of whole-exome and whole-genome sequencing data
Human genomic whole-exome sequencing libraries were generated using the SureSelect XT kit specific for the Illumina HiSeq instrument (Agilent Technologies), followed by exome enrichment using the SureSelect XT Human All Exon V5 without UTRs bait set. The resulting exome enriched libraries of tumour and normal samples were then sequenced using paired end 100 cycle sequencing on a NovaSeq 6000 (Illumina) according to manufacturer's instructions. Whole-genome sequencing libraries were constructed using the TruSeq DNA PCR-free sample preparation kit according to the manufacturer's instructions. Tumour and normal samples were sequenced on an Illumina HiSeq 2500 instrument as previously described 36 . Somatic SNVs and INDELs were determined via the Mutect2 algorithm as implemented in GATK version 4.0. Coding and splice-related variants were subsequently annotated utilizing the Medal Ceremony annotation pipeline. Additionally, all reported somatic variants were manually curated in IGV.

Human single-cell RNA-seq data generation and processing
Whole transcriptome amplification, library construction and sequencing were performed as previously described following the Smart-seq2 modified protocol 13 . Expression levels were quantified as E i,j =log 2 (TPM i,j /10+1), where TPM i,j refers to transcript-per-million for gene i in sample j, as calculated by RSEM 37 . TPM values were divided by 10 since we estimate the complexity of single cell libraries in the order of 100,000 transcripts and would like to avoid counting each transcript ~10 times, as would be the case with TPM, which may inflate the difference between the expression level of a gene in cells in which the gene is detected and those in which it is not detected.
To filter out poor quality cells, we first removed cells for which less than 2,500 genes were detected. For each processed 96-well plate, we then determined the average number of genes detected per cell minus two times its standard deviation. We then additionally filtered out the cells that were below that threshold. For the remaining cells, we calculated the aggregate expression of each gene as E a (i)=log 2 (average(TPM i,1…n )+1), and excluded genes with E a <4. In each subgroup (WNT, SHH, and Group 3/4), we defined relative expression by centering the expression levels, Er i,j =E i,j -average[E i,1…n ] for the remaining cells and genes.
On average, we detected ~4,500 genes per cell. Gene expression values were uploaded to the Gene Expression Omnibus (accession number GSE119926).
Pearson correlation coefficients between expression profiles of cells that passed quality filtering was calculated using centered gene expression levels (for each subgroup separately, shown in Figures 2a, 3a, 4a). Cells were ordered by hierarchical clustering using 1correlation coefficient as the distance and Ward's linkage, within each sample or genetic subclone (for samples SJ99 and BCH825).

Identification of copy number variations in single-cell data
Copy number variations were estimated as previously described 13 by applying a moving average to the relative expression, with a sliding window of 100 genes within each chromosome after sorting genes by their chromosomal location (shown in EDF 2). Nonmalignant tumour cells were determined by unsupervised clustering of all single-cell-derived copy-number profiles for each sample with 190 copy-number profiles derived from two nonmalignant cell types (tumour-associated oligodendrocytes and immune cells 14 ). For the majority of tumours (21/25), most of the cells did not cluster with the non-malignant cells but formed their own cluster(s) and showed clear evidence of CNVs. A small fraction of tumor cells clustered with the non-malignant cells (<4%). Given the high percentage of malignant cells in these tumors, we decided to classify all cells from the remaining four tumours (MUV41, SJ577, MUV34, and SJ625) as malignant. For two samples (SJ99 and BCH825) genetic subclones were identified based on their CNV profiles.

Identification of single-nucleotide variants in single-cell data
To detect mutant transcripts in our full-length scRNA-seq expression data (shown in Figure  2c and ED Figure 3c-f), sequencing reads were first aligned to the human genome build hg19 using STAR version 2.5.1b. RefSeq gene annotations were supplied to guide alignment. Variants were then quantified in each single cell at the genomic position in which they were detected in the WGS/WES data using samtools mpileup version 1.3. For some genes multiple variants were detected (e.g. four different variants were detected for CTNNB1 in WNT-MB) and quantified separately. To detect mutant and wild-type transcripts we required one or more supporting reads. We then filtered variants that were detected as mutant in less than three cells, or that were considered erroneously called as they were detected at elevated frequency in samples in which they were not detected in the genome sequencing data. A total of 82 variants remained after this filtering step. Mutations in highly expressed transcripts were detected in the majority of cells from the respective sample (e.g. OTX2 Q103R mutation in MUV39). Mutations in less highly expressed transcripts were detected less frequently. Only a small number of mutant transcripts were detected in samples in which they were not detected by genome sequencing (e.g. only 10 mutant transcripts were detected for the respective other CTNNB1 variants in the WNT-MB single cells; Figure 2c), illustrating the high specificity of our approach.

Identification of intra-tumour NMF programs and cellular hierarchies
Transcriptional programs were determined as previously described 14 by applying Nonnegative Matrix Factorization (NMF) to the centered expression data 38 . Negative values were converted to zero. Analysis was performed for each sample and subclone individually (excluding samples for which less than 100 cells were profiled), using only the malignant cells and setting the number of factors to 4 for WNT and 3 for SHH and Group 3/4 tumours. For each of the resulting factors we considered the 30 genes with the highest NMF scores to be characteristic of that factor (provided in Supplementary Table 2a). All single cells within the WNT, SHH, or Group 3/4 subgroups where then scored for these NMF programs (as described below, shown in ED Figures 5a, 6a, and 8a). Hierarchical clustering, with one minus Pearson correlation as the distance metric and Ward's linkage, of the scores for each program revealed four (WNT subgroup) or three (SHH and Group 3/4 subgroups) main correlated sets of programs. The 30 genes with the highest average NMF score within each correlated program set (excluding ribosomal protein genes) were then used to define a total of 10 subgroup-specific meta-programs (provided in Supplementary Table 2b).
To interpret the characteristics of each meta-program, we manually inspected their underlying gene signatures. Additionally, we tested for enrichment of described gene sets (GO biological processes cell cycle and neuron differentiation, KEGG hedgehog signaling pathway, and manually curated ribosomal proteins and translational initiation/elongation factors) in each meta-program using Fisher's exact test.

Generation of single-cell program expression scores
Single-cell expression scores were generated in a similar way as described previously 13 . Given a set of genes (G j ) reflecting a NMF program or meta-program, we calculate for each cell i, a score, SC j (i), quantifying the relative expression of G j in cell i, as the average relative expression (Er) of the genes in G j , compared to the average relative expression of a control gene-set (G j cont ): SC j (i)=average[Er(G j ,i)]-average[Er(G j cont ,i)]. The control geneset contains 100 genes with the most similar aggregate expression level. In this way, the control gene-set has a comparable distribution of expression levels to that of the considered gene-set, and the control gene set is 100-fold larger, such that its average expression is analogous to averaging over 100 randomly-selected gene-sets of the same size as the considered gene-set.
Single cells were assigned to different cell populations based on the maximum expression score for their respective subgroup-specific meta-programs, excluding the cycling programs. The fraction of cells per tumour sample assigned to each cell population is provided in Supplementary Table 2d. Scores for the cycling programs were binarized into cycling and non-cycling (larger and smaller than 1, respectively). For illustration of the cellular hierarchies in SHH-MB, scores were normalized by minimizing the average minimum difference of all cells per sample to −1 or 1 (shown in Figure 3c). For the pan-subgroup analysis of all malignant medulloblastoma cells we re-centered expression values across the dataset and calculated expression scores for each of the 10 meta-programs. The pairwise correlation of expression scores is shown in Figure 6a.

Determination of cell population-specific genes
For comparison of cell populations in WNT-tumours, we calculated the average expression level of all cells per population (log 2 -transformed, un-centered expression levels were used). For this analysis we excluded cells from BCH807, as it was very different from the other four WNT-MBs and represents an atypical case (without cells scoring highest for the WNT-B meta-program, highly proliferative, negative staining for nuclear beta-catenin and lacking monosomy of chromosome 6). We then determined all genes with a difference smaller or larger than 1 between the average log 2 -transformed expression levels when comparing the undifferentiated proliferative population (highest for meta-program WNT-B) against the neuron-like population (WNT-C) or undifferentiated post-mitotic population (WNT-D). A total of 511 genes were identified in this way (provided in Supplementary Table 2c; ED Figure 5c).
For comparison of neuron-like cell populations between medulloblastoma subgroups, we first determined genes that were specific to any of the neuron-like populations. For every subgroup (WNT, SHH, and Group3/4), the average expression level of all neuron-like cells was compared to the average expression level of all undifferentiated cells from each subgroup, determining genes with a difference larger than 1 between the average log 2tranformed expression levels. This way we determined a total of 260 genes that were specific to the neuron-like cell population of at least one subgroup (provided in Supplementary Table 3). Genes that were specific to two or three subgroups were grouped as shared genes (shown in Figure 6c and ED Figure 9e).

RNA in situ hybridization
Paraffin-embedded tissue sections from two WNT-MB tumours of the single-cell cohort (SJ99 and SJ129) were obtained from St. Jude Children's Research Hospital. Sections were mounted on glass slides and stored at −80°C. Slides were stained using the RNAscope 2.5 HD Duplex Detection Kit (Advanced Cell Diagnostics (ACD), 322430). Slides were baked for 1 hour at 60°C, deparaffinized and dehydrated with xylene and ethanol. The tissue was pretreated with RNAscope Hydrogen Peroxide (ACD, 322335) for 10 minutes at room temperature and RNAscope Target Retrieval Reagent (ACD, 322000) for 15 minutes at 98°C. RNAscope Protease Plus (ACD, 322331) was then applied to the tissue for 30 minutes at 40°C. Hybridization probes were prepared by diluting the C2 probe (red) 1:50 into the C1 probe (green). Advanced Cell Technologies RNAscope Target Probes used included Hs-MKI67 (ACD, 591771 and 591771-C2), Hs-DKK2 (ACD, 531131-C2), Hs-STMN2 (ACD, 525211-C2), Hs-ZFP36 (ACD, 427351) and Hs-EGR1 (ACD, 457671). Probes were added to the tissue and hybridized for 2 hours at 40°C. A series of 10 amplification steps were performed using instructions and reagents provided in the RNAscope 2.5 HD Duplex Detection Kit. Tissue was counterstained with Gill's hematoxylin for 25 seconds at room temperature followed by mounting with VectaMount mounting media (Vector Laboratories).

Murine single-cell RNA-seq data generation and processing
Single cells from developing murine cerebellar tissue were processed using the microfluidics-based 10x Chromium protocol, as previously described 21 . In brief, single cells were prepared using the Chromium v1 Single Cell 3' Library and Gel Bead Kit according to the manufacturers specifications. Quantification and quality checks for the library were performed using an Agilent Technologies DNA 1000 chip. Libraries were sequenced on an Illumina HiSeq 2500 machine. Raw sequencing data has been uploaded to the European Nucleotide Archive (accession PRJEB23051).
Murine developing cerebellum cells were filtered and normalized using the scanpy python package 39 . Genes expressed in less than 50 cells and cells expressing less than 200 genes were removed. Additionally, cells with less than 524 and greater than 3,206 total counts (±3 MADs) were removed. Furthermore, those cells with greater than 5% of their total counts mapping to mitochondrial genes were removed. Gene expression values were then divided by the total number of transcripts and multiplied by 10,000. Normalized values were calculated by natural-log transforming these values. We calculated scaled expression (zscores for each gene) for downstream analysis.

Identification of cell types in developing murine cerebellum
The scanpy package implemented in python was applied to identify cell types among 82,228 cells expressing a total of 16,475 genes. After two rounds of clustering (using the Louvain method), populations predicted to be of non-cerebellar origins were excluded. Removed populations were enriched for haemoglobin, oligodendrocycte, and/or immune associated genes. The remaining 78,156 cells were visualized by t-SNE (using the first 100 principle components as input) and clustered a third time. We then merged clusters if mantel spearman correlation between gene distance matrices (using Manhattan distance) was greater than 0.9 (shown in Figure 1c and ED Figure 4c). Resulting clusters, in conjunction with marker genes, were used to identify major cell types in the developing cerebellum.

Integrated analysis of murine and human data sets
Gene expression matrices for human and mouse data sets were restricted to the 16,919 highconfidence homologous genes with gene order conservation (GOC) and whole genome alignment scores (WGA) greater than 75 percent, as defined by Ensembl. We removed genes without expression in at least 200 cells and filtered out those with gene dispersion across cells/samples less than equal to zero in each dataset. We also regressed out individualspecific effects in the single-cell data.
For canonical correlation analysis (CCA) the first 30 canonical correlation vectors were calculated to project each expression matrix into the maximally correlated subspace, as similarly described previously 40 . Briefly, CCA is implemented as singular value decomposition, by implicitly restarted Lanczos bidiagonalization algorithm, of a distance matrix between two gene expression matrices.
We adopted a correlation of differential expression approach to measure similarity between biological groups in two different studies. Such a procedure has previously been shown to be effective in implicating cellular origins for WNT and SHH medulloblastoma subgroups 41 .
Gene expression for each cell, or centroids for each cluster when at the cluster level, is subtracted by mean gene expression of all other cluster centroids to determine differential expression. Cosine distance is then used to calculate correlations between differential expression vectors between studies as a metric for similarity (shown in Figures 1d, 3d, 3e, and ED Figures 4f, 7j, 10a, 10f). Significance is assessed by 10,000 permutations, followed by FDR correction, for cluster labels of interest. Genes driving CCA differential correlations between human and murine datasets were investigated by identifying genes both differentially expressed in the cell type of interest (Mann-Whitney U test) and correlated with the CCA correlation (Pearson correlation). Significant genes were those predicted to drive CCA differential correlations.
Non-negative Matrix Factorization applied to the centered murine expression data, with negative values assigned to zero and rank set to two, determined an undifferentiated and differentiated program. Both programs were projected onto a centered dataset of interest, scaled to a range of zero and one, then differentiated programs were subtracted from undifferentiated programs to calculate differentiation scores (shown in ED Figure 4b 42 ). detected in a single cell, only the mutant transcript is shown. Variants were initially detected by genome sequencing and subsequently quantified in the scRNA-seq data. Sample BCH807 was not subjected to genome sequencing, and the CTNNB1 variant was manually detected by examining scRNA-seq alignments. Mutations are detected almost exclusively in single cells from samples in which they were detected by genome sequencing, illustrating the high specificity of single-cell variant detection. e, Heatmap shows mutant and wild-type transcripts for 15 variants in each cell (n=1,135, rows) of the SHH-MB dataset. Sample SJ454 was not subjected to genome sequencing, and the TP53 mutation was manually identified by examining scRNA-seq alignments. f, Heatmap shows mutant and wild-type transcripts for 28 variants in each cell (n=3,172, rows) of the Group 3/4-MB samples that were subjected to genome sequencing. during CB development (images obtained from the Allen Brain Atlas). j, Radar plot showing the CCA cosine correlation coefficients between each murine CB cell type and the MB single-cell cohort from cells scoring highest for meta-programs SHH-B (GNP-like cells) and SHH-C (Granule neuron-like cells).
shows that the same programs and cell populations that are identified by the NMF analysis are also supported by PCA and t-SNE clustering. Furthermore, no additional programs and cell populations are identified (starting from PC4 components are less informative