{'fraction_non_alphanumeric': 0.03425762910798122, 'fraction_numerical': 0.014646909233176839, 'mean_word_length': 5.453684708852769, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Single-Cell RNA Sequencing Unravels Upregulation of Immune Cell Crosstalk in Relapsed Pediatric Ependymoma

Ependymoma (EPN) is a malignant glial tumor occurring throughout central nervous system, which commonly presents in children. Although recent studies have characterized EPN samples at both the bulk and single-cell level, intratumoral heterogeneity across subclones remains a confounding factor that impedes understanding of EPN biology. In this study, we generated a high-resolution single-cell dataset of pediatric ependymoma with a particular focus on the comparison of subclone differences within tumors and showed upregulation of cilium-associated genes in more highly differentiated subclone populations. As a proxy to traditional pseudotime analysis, we applied a novel trajectory scoring method to reveal cellular compositions associated with poor survival outcomes across primary and relapsed patients. Furthermore, we identified putative cell–cell communication features between relapsed and primary samples and showed upregulation of pathways associated with immune cell crosstalk. Our results revealed both inter- and intratumoral heterogeneity in EPN and provided a framework for studying transcriptomic signatures of individual subclones at single-cell resolution.


INTRODUCTION
Ependymomas (EPNs) are primary tumors of the central nervous system that commonly present in childhood. Although the diagnosis and stratification of EPN patients have been facilitated by identification of nine EPN molecular groups from genome-wide DNA methylation studies (1,2), EPN patients display a high prevalence of relapse, and recurrence typically results in much poorer outcomes (3). Between molecular groups, posterior fossa group A (PFA) EPN and supratentorial (ST) EPN with C11orf95-RELA fusions have been reported to show worse prognoses than PF group B, ST-EPN with YAP1fusions, and spinal EPNs (1,4). With the advent of highthroughput single-cell RNA sequencing technologies, recent studies have provided resources for understanding the molecular landscape of EPN, revealing a cellular hierarchy in these tumor cells characterized by an undifferentiated progenitor population, which transitions into distinct cell lineages including neuronal precursor-like, glial progenitor-like, and ependymallike cells (5,6). However, these studies focused on characterizing EPN from different major molecular groups and anatomical locations but did not reveal the differences in molecular signatures between subclones within individual tumors.
To date, cellular heterogeneity has typically been viewed as a consequence of hyperproliferation and genomic instability that can give rise to intratumoral subclones during tumor progression (7). In the context of EPN, previous studies have revealed potential signaling pathways involved in driving the expansion of therapy-resistant EPN subclones, which contribute to tumor relapse and disease progression (8). These findings suggest that unravelling subclone-to-subclone variability at single-cell resolution could shed light on the molecular mechanisms underpinning EPN pathogenesis.
In particular, mutant genotypes can grant selective advantage on specific cellular subclones, leading to their outgrowth and allowing them to establish dominance in different types of tissue environments (7). To date, cellular heterogeneity has been viewed as a consequence of hyperproliferation and genomic instability that can give rise to intratumoral subclones during tumor progression (7). For example, genomic instability associated with intratumoral heterogeneity can manifest in the form of extensive subclonal evolution demonstrated to be correlated with a higher risk of recurrence or death in nonsmall-cell lung cancer (9). Indeed, the emergence of subclonal diversity is a fundamental characteristic of intratumoral heterogeneity and has been found to be significantly associated with patient survival across diverse cancer types in a pan-cancer study, including lower-grade glioma and glioblastoma multiforme (10). Interestingly, subclones of glioblastoma were shown to display remarkable heterogeneity of drug resistance wherein characteristics of coexisting subclones could be linked to distinct drug sensitivity profiles, hinting at the therapeutic potential of targeted treatments for tumor subclones associated with differential survival outcomes (11).
To interrogate subclonal heterogeneity within tumor populations in pediatric EPN, we performed single-cell RNAseq on EPN samples across PF-A and ST regions. Using a deconvolution approach provided by inferCNV to compare subclones in a single PF-A sample as a proof-of-concept, we identified a subclone-specific cilia-associated program within an individual PF-A EPN sample. We further incorporated a trajectory score analysis to predict correlations between survival outcomes in EPN and molecular characteristics, and primary and recurrent tumor populations. Finally, we identified cell-cell communication features between relapsed and primary EPN samples and show upregulation of pathways associated with immune cell crosstalk. Our results reveal gene expression profiles associated with subclonal variability, providing a framework for studies on transcriptomic signatures of brain tumor subclones.

Single-Cell Transcriptomic Profiling Reveals Stemness Signature Differences Between Intratumoral Subclones in PF-EPN
To characterize intratumoral heterogeneity in human EPN, we performed single-cell RNA sequencing on four EPN patients with the 10× Genomics platform and profiled the transcriptome of 35,102 qualified cells with an average of 3,472 genes per cell (Supplementary Table S1). The cells with high percentage (>12%) of mitochondrial genes and low number (<1,500) of genes and those regarded as doublets were removed from the dataset for subsequent analysis (Supplementary Figure S1). We first performed copy number variation (CNV) analysis to distinguish neoplastic cells from non-malignant (NM) cells and identified putative subpopulations of malignant tumor cells with high CNV in each tumor samples ( Figure 1A). Moreover, we applied the approach to score the differentiation state of malignant tumor cells using a panel of differentiationassociated genes, which revealed that CNV-inferred neoplastic cell populations were in a less differentiated state ( Figure 1B).
After subsetting malignant cells from our single-cell dataset, we further examined the presence of different subclones within individual EPN samples. Notably, the inference of subclonal CNV events uncovered the presence of two putative subclones within the PFA-EPN sample GTE009, based on the hierarchical clustering of inferCNV matrix (Figures 1C, D; Supplementary Figure S2). Subsequently, the differentially expressed genes (DEGs) of various cell types included neural stem cell (NSC), neuron (NEU), radial glial cell (RGC), oligodendrocyte precursor cell (OPC), oligodendrocyte (OD), astrocyte (AS), ependymocyte (EpC), endothelial cells (EC), microglia (Mic), and T cells (TC) from human and rodent embryonic and postnatal cortex scRNAseq data (Supplementary Table S2) and were enriched as signatures to classify cell types in our tumor samples (19)(20)(21)(22)(23)(24) (see Materials and Methods). We further applied signature enrichment (SE) analysis and reversed-SE (rSE) for the accuracy of cell-type classification (Supplementary Figures S3,  S4, S5F), which was supported by correlation analysis (Supplementary Figure S5). Similar to previously published EPN single-cell datasets (5,6), our analysis identified NSC-, EpC-, NEU-, RGC-, OPC-, OD-, and AS-like cells in malignant populations, and endogenous NEUs, ECs, Mic, and other cells ( Figures 1E-F; Supplementary Figure S4). Although both subclones in GTE009 sample encompassed the same malignant cell types, cell-type composition analysis of these subclones interestingly showed a lower proportion of NSC-like cells in subclone 1 ( Figure 1G) with a corresponding increase in EpClike cells. Moreover, Gene Ontology (GO) analysis revealed enrichment of cilium-related terms based on the DEGs in cells from subclone 1 (Figures 1H, I) consistent with studies highlighting the role of cilium-related genes in disease processes linked to tumorigenesis (25,26), suggesting that molecular characteristics correlated with EPN pathogenesis may be associated with intratumoral subclonal heterogeneity.

Upregulation of Cilium-Associated Genes Is Associated With CNV Amplification in Highly Differentiated EPN Subpopulations
In spite of transcriptomic and CNV differences between intratumoral subclones, RNA velocity analysis revealed a classic molecular trajectory originating from NSC-like cells to EpC-like cells similar to previously published single-cell EPN datasets (5,6) (Figures 2A-C; Supplementary Figure S6). Given the marked differences in cell states between populations from separate intratumoral subclones, we further examined the expression level of differentially expressed cilium-related genes and their corresponding CNV score (Supplementary Table S3). We found that cilium-related genes possessed both higher expression and genome amplification in cells from the more highly differentiated subclone 1 compared to subclone 2, suggesting a correlation between CNV amplification and genes associated with more differentiated cell states; for example, the expression level of the cilium-related gene DYNC2H1, which has  (12). (F) Heatmap of DEGs calculated by cell types and pathogenic sites from scRNA-seq data in this study and bulk-DEGs. Aforementioned bulk-DEGs were calculated by pathogenic sites from online bulk-seq data [Gene Expression Omnibus (13,14); Seq: GSE89448 (15); Array: GSE64415 (1,16,17); aligned to human reference genome GRCh38(hg38)] through DESeq2 (18). (G) Histogram of cell types in PF-EPN sample GTE009 colored by cell types in percentage and outlined by subclone annotation showing significant difference (p value = 7.975e−05) in cell-type proportions using asymptotic two-sample Fisher-Pitman permutation test. (H) Workflow of gene ontology enrichment analysis comparison between PF-EPN sample GTE009 subclones 1 and 2. (I) Gene ontology analysis of upregulated genes in the PF-EPN sample GTE009 subclone 1 compared to the subclone 2 ordered by adjusted p-value. The "*" means significant difference (p-value): * p < 0.05; ** P < 0.01; *** p < 0.001; **** p < 0.0001. been implicated in the formation of hypothalamic hamartoma (28), was both upregulated and amplified in chromosome 11 of subclone 1 ( Figures 2D-F). Indeed, other genes associated with cilium-related terms in GO analysis (27) were found to be more highly expressed in subclone 1 compared to subclone 2 ( Figure 2G). The marker genes of EpC-like cells in subclone 1 compared to that of subclone 2 (EpC-Sub1) may be indicative of more mature cellular populations, which was consistent with an inverse correlation with the undifferentiated score (r = −0.65) in EpC-like cells ( Figures 2H ,I). These findings suggest that CNV amplification of EpC-related genes is correlated with differentiation of malignant cells, manifesting as alterations in cell composition within intratumoral subclones while maintaining cardinal features of EPN tumorigenesis. This may be relevant in the context of diagnosing malignant EPN samples, as previous findings have shown direct evidence for the roles of CNV-amplified genes in preventing differentiation, inhibiting cell death, and promoting tumor growth, which were in turn correlated with poor patient outcomes (29).

Trajectory Score Analysis Identifies Cellular Compositions Associated With Worse Survival Outcomes in EPN
To further analyze the cellular subpopulations in subclones 1 and 2, we performed GO analysis on the DEGs between subclones 1 and 2 in NSC-like and EpC-like cells, respectively. We detected high expression of cell-cycle-related genes in NSC-like cells from subclone 2, and enrichment of cilium-related genes in EpC-like cells in subclone 1 ( Figure 3A). Given that previous studies have demonstrated that high expression of NSC-like cells and low expression of EpC-like cells are correlated with poor patient survival and vice versa (5, 6), we hypothesized that performing a combinatorial analysis integrating information from both subclone and cell types may help in predicting EPN patient outcomes.
To provide a numerical representation of molecular trajectory information at the single-cell level, we developed a trajectory score for downstream analyses ( Figure 3B Supplementary Table S2) in PF-EPN sample GTE009. (I) Pearson correlation between undifferentiated score and normalized average expression of EpC-Sub1 (p < 0.0001) in the PF-EPN sample GTE009. The "*" means significant difference (p-value): * p < 0.05; ** P < 0.01; *** p < 0.001; **** p < 0.0001. the subclone 2; NSC-Sub2) markers was subtracted by differentiated-EpC (EpC-Sub1) markers. The trajectory score was presented in the tSNE plot ( Figure 3C), which resembled the trajectory analysis (see Figure 2B). The trajectory score allows for the identification of differentiation state at the single-cell level, which complements the molecular trajectory analysis in merged samples. Notably, well-characterized stemness-associated markers (10, 31) such as FTL, LGALS1, MEG3, MEST, TUBB, TMSB4X, and STMN1 were found in the DEGs of trajectory-high group compared to trajectory-low group ( Supplementary Table S4), while cilium-related terms were enriched in the trajectory-low group compared to the trajectory-high group (Supplementary Figure S7A).
Based on the hypothesis that trajectory score was directly correlated with survival outcome and could be easily used to predict the prognosis of EPNs by its simple calculation method, we applied this scoring methodology to two published EPN scRNA-seq datasets (33 EPN patients in total) and showed that our results are consistent with the respective survival outcomes reported in these studies based the respective cell-type composition of individual samples (5, 6) ( Figure 3D). On the contrary, the application of the aforementioned undifferentiated score led to relatively more inconsistent results (Supplementary Figure S7B), suggesting that trajectory score analysis could be a useful tool to investigate EPN prognosis. Indeed, samples with a high trajectory score were found to have correspondingly poorer survival outcomes in published PF-EPN samples and PF/ST-EPN (5, 6) (Supplementary Figure S7C), although the comparison did not reach significant difference due to the small sample size. Likewise, a higher percentage of recurrent patients compared to primary  (5,6) of EPN using previously defined cell-type annotations (p < 0.0001; Kruskal-Wallis test). The "*" means significant difference (p-value): * p < 0.05; ** P < 0.01; *** p < 0.001; **** p < 0.0001. patients had higher trajectory score (Supplementary Figure S7D), supporting the association between EPN relapse state and overall survival outcomes of this disease.

Cell Compositions Correlated With Poor Prognosis in EPN Recurrent Patients Are Revealed by Trajectory Score Analysis
Relapse rates for EPN can be as high as one-third of the patients (3), and relapse is known to lead to significantly worse survival outcomes based on published data (5, 6) ( Figure 4A). To determine transcriptomic signatures at the single-cell level associated with poor survival outcomes in EPN relapse, we compared the cell-type composition in 36 EPN patients (5,6) and revealed a significant difference with higher percentage of NSC-like cells in relapse patients ( Figure 4B). Similarly, we found a significantly higher trajectory score in recurrent NSC-like cells compared to primary NSC-like cells, and a similar trend was observed for EpC-like cells ( Figure 4C). Given the association between high stemness and poor survival outcome found in relapse patients based on the trajectory score, we further performed GO analysis on the DEGs between NSC-like cells in primary and recurrent patients and revealed enrichment of cilium-and immune-related terms ( Figure 4D). This suggests that NSC-like cells from recurrent patients were not only in a more immature cell state but also that there is a likelihood of extensive immune cell crosstalk within these cellular populations, consistent with previous findings reporting the association between cell-cell communication in immune cells and tumor progression (32,33).

Relapsed EPN Show Upregulation of Distinct Signaling Pathways Associated With Immune Cell Crosstalk
Based on studies demonstrating the role of tumor-infiltrating NM cells such as Mic in brain tumor (32,33), we performed crosstalk analysis to investigate cell-cell interactions between the different cell-types profiled. Although there has been increasing interest in examining communication patterns between cell populations using scRNA-seq, crosstalk analysis on intracranial EPN samples has not been extensively studied given the relatively lower number of tumor-infiltrating NM cells profiled in previous datasets (5,6). To investigate cell-cell interactions in EPN, we first examined the expression of ligand-receptor pairs in different cell-types across four EPN samples ( Figure 5A Table S5). For example, we observed strong outcoming events from NSC-like cells towards other cells (Supplementary Figure S8B). Moreover, we uncovered the overlap in simulated spatial 3D position between NSC-like cells and Mic by CSOMAP (34) ( Figure 5C; Supplementary Figure S8C). Interestingly, these events that had higher expression in recurrent samples than that in primary samples in 36 EPN patients ( Figure 5D) showed interaction between NSC-like cells and Mic (see Supplementary  Figures S8D, E), including the MK pathway that promoted brain tumor growth (32,33) and the EGFR pathway that inhibited glioblastoma invasion via pharmacological inhibition of EGFR (35). For example, MDK (ligand) and NCL (receptor) were highly expressed in the tumor microenvironment (TME) of recurrent samples ( Figure 5E), consistent with previous studies implicating the roles of this ligand-receptor pair in tumorigenesis (36,37) and the MK-deficiency-reduced tissue infiltration of microglia (38). To further elucidate signaling pathways involved in crosstalk between normal and malignant cells, the inferred gene regulatory networks also revealed multiple pathways shared in crosstalk (Supplementary Figure  S9

DISCUSSION
The increasing accessibility of scRNA-seq technologies has accelerated our understanding of cellular function in health and disease. Here, we generated a high-resolution EPN singlecell dataset with a particular focus on the comparison of subclone differences within tumor populations. Our analysis on four EPN samples profiled 35,102 single-cell transcriptomes and uncovered 17 major cell types including NSC-like, EpC-like, and microglia populations that are present across different EPN groups. We further revealed differences in cell proportions within highly differentiated populations within tumor subclone cells by integrating CNV pattern analysis with singletranscriptome data in this study and previously published datasets. Additionally, differential gene expression analysis also identified gene programs associated with tumor subclones and survival outcomes.
Treatment of heterogeneous tumors such as EPN can favor selection of resistant subclones given that different subclones respond differently to intrinsic and extrinsic signaling cues. EPN relapses after surgical resection and treatment of EPN are common and have poor outcomes, and recent findings have demonstrated that administration of radiation and chemotherapy can lead to a significant increase in EPN mutational burden in conjunction with changes to the tumor subclonal architecture, without eliminating the original founding clone (8). In this study, we report the presence of subclones within a single EPN tumor sample characterized by molecular signatures reflecting different stages of cellular differentiation. This suggests that in addition to stemness signature gradients between tumors (5, 6), intratumoral heterogeneity can also be uncovered within individual samples containing multiple subclones. Importantly, further interrogation of the subclones identified from CNV analysis also showed that EPN subpopulations that were more differentiated exhibited an increase in cilium-associated genes. For example, in the sample GTE009, 6% of 494 DEGs from EpC-like cells compared to the other malignant tumor cells were found to be overlapped with genes related to cilium assembly, organization, and movement such as PIFO, ZMYND10, DNAH9, TEKT1, DYNLL1, SPEF1, MNS1, DNAAF1, IQCG, SPAG16, and FOXJ1. Indeed, ciliary signaling is known to be a mediator of paracellular signals controlling cancer metastatic processes and responses to therapy, and mutations leading to defects or structural abnormalities in cilia have been shown to be directly correlated with cancer pathogenesis (25). Given that ependymal cells in EPN are multiciliated cells (26), it is plausible that changes in ciliation of EPN subpopulations and/or cells of the tumor TME during EPN development can contribute to disparities in outcomes within tumors of the same molecular group.
Notably, our analysis on the sample GTE009 showed that cells from different subclones on the tSNE plot appear to cluster together without displaying pronounced segregation. This is largely due to the clear difference between the selected transcriptome for tSNE reduction and transcriptome-based genome for subclonal classification. For example, the inferCNV analysis selects the expression (RNA) of a panel of genes located consecutively along a chromosome to establish the CNV situation (DNA) of those genes; however, these genes are not all selected as variable features for tSNE reduction. Hence, this results in a difference between the visualization of the transcriptome-based tSNE reduction and the inferred subclones, meaning that subclones may not polarize within the large cluster. Consistent with our analysis, previous findings on glioblastoma intratumoral heterogeneity revealed that most subclones within glioblastoma samples appear to show an overlapping distribution, even though a small proportion of subclonal populations could also exhibit distinct cell states (10). Additionally, since our subclone definition in this study was based on only CNV analysis due to inherent limitations of scRNA-seq data, future studies could also encompass the analysis of possible subclones differentiated by single nucleotide alterations to decipher the impacts of such changes on EPN metastasis.
To complement classical pseudotime molecular trajectory methodologies (39-42), we applied a curated trajectory score to our single-cell dataset to integrate EPN subclone and cell-type information in our analysis and found that higher trajectory scores were found in patients with EPN samples exhibiting elevated stemness signatures [FTL, LGALS1, MEG3, MEST, TUBB, TMSB4X, and STMN1 (10, 31)], which is also associated with worse prognoses. This trajectory score analysis allows for a numerical representation of EPN trajectory stages at the single-cell level to facilitate comparison with datasets of interest. Moreover, our analysis takes into consideration both cell types and stemness signatures and can be easily applied to other transcriptomic datasets of interest to quantify the relative correlation degree between survival outcomes and stemness signatures across different samples. As a proof-of-concept, we performed validation of our dataset with previously published results on EPN samples (derived from 36 patients in total) and indeed demonstrated that this analysis reveals consistent trends in EPN survival outcomes. In addition to using trajectory score analysis to examine the association between EPN stemness signatures and survival outcomes, we further applied this method to compare differences in survival outcomes between primary and recurrent EPN samples. Previous findings have identified an enrichment of undifferentiated programs (NSC-like) in recurrent PF-EPN relative to primary PF-EPN samples from comparing three matched samples at the single-cell level (5). Here, we find that this is consistent across the single-cell transcriptome of 36 published EPN samples, as trajectory score analysis shows a clear significant difference in cell-type composition between recurrent and primary EPN samples with a higher percentage of NSC-like cells in recurrent EPN. Further analysis of subpopulations within recurrent and primary samples revealed that recurrent samples also show higher trajectory score in NSC-and in EpC-like cells compared to the corresponding cell types in primary samples. These findings suggest that trajectory score analysis can uncover multiple types of association in EPN samples in a quantifiable form, such as the correlation between stemness and tumor occurrence with patient mortality.
In EPN and other brain cancers, there is increasing evidence that the brain TME functions as a key regulator of cancer progression in brain malignancies (43). Hence, we performed cell-cell communication analysis on our EPN single-cell transcriptomes to assess the crosstalk between different cell types in EPN, given that the TME contains non-cancerous cell types such as pericytes, endothelial cells, and immune cells in addition to cancer cells. We revealed putative interactions between malignant cells and tumor-infiltrating NM cells, such as enrichment in interactions between NSC-like cells and microglia. For example, the MK pathway, which has been implicated in brain tumor pathogenesis (32,33), was not only found to be significantly upregulated in this study's dataset but also showed higher expression in recurrent samples compared to primary samples based on 36 published EPN single-cell transcriptomes. Indeed, MK deficiency has been shown to reduce tissue infiltration of microglia, leading to reduced neuroinflammation and apoptosis (38). Given that inflammatory crosstalk with immune cells has previously been shown to play a key role in driving tumor growth in the EPN microenvironment (44), these findings suggest that immune cell crosstalk analysis may serve as a useful resource for the identification of candidate genes for future in vitro and in vivo validation studies.
In summary, we report a curated EPN atlas focusing on the comparison of intratumoral heterogeneity in this study. We used an integrative analysis approach to show both changes in celltype composition and cell-type-specific gene expression associated with different tumor groups and subclones. Moreover, we also applied a novel trajectory scoring method as a parallel tool to traditional molecular trajectory analysis and demonstrated its robustness in recapitulating survival outcomes within individual EPN samples and across primary and recurrent tumors. This approach will complement existing published datasets and provide valuable insights into cell-type-specific properties of EPN, laying the foundation for therapeutic treatments of this disease.

EPN Sample Preparation for scRNA-Seq
Fresh tumor samples were processed as previously described with minor modifications (45). A fresh EPN tissue was excised by

Data Processing of scRNA-seq
The sequencing data were processed by CellRanger v3.1.0 with the reference genome hg19-3.0.0 to generate filtered expression matrixes, which were analyzed using Seurat v3.2.0 (46). Doublet Finder (47) was first applied to erase doublets with default settings. Genes detected in at least 10 cells were used for analysis, and cells that possessed transcription numbers fewer than 1,500 or cells with mitochondrial genes taking up more than 12% of reads were removed. After normalizing the data, we used 5,000 highly variable features for downstream analysis, and cell cycle variation was regressed out as previously described (30). tSNE analysis was performed with the top 50 significant principal components from principal component analysis, and cells were clustered using "FindClusters" function based on tSNE reduction. CNV analysis was performed using inferCNV of the Trinity CTAT Project (48) (https://github.com/broadinstitute/ inferCNV) as described in the following section. Calculation of the undifferentiated score was performed using CytoTRACE (49) as previously described (45). The EPN single-cell datasets covering malignant cells and tumor-infiltrating NM cells were analyzed by "iCytoTRACE" function with the quantification of the number of expressed genes as an indicator of differentiation potential; genes' associated differentiation (deduced by CytoTRACE) and their relative expression levels were used to perform this calculation. Kyoto Encyclopedia of Genes and Genomes (KEGG)/GO analysis was performed based on the DEGs between malignant tumor cells and NM cells of the merged dataset using "FindAllMarkers" function in Seurat and clusterProfiler (50). The correlation of malignant tumor cells and NM cells among samples was in favor of malignancy separation (Supplementary Figure S1E) utilizing "cor" function of stats package in R v3.6.3.

Estimating CNVs in scRNA-seq Data
The initial CNVs of single cells were estimated from their wholegenome wide expression level by inferCNV (30,48). To perform comparison across samples, 300 cells of OPCs were sampled from the GTE009 as a common reference, and all the nonimmune cells of each sample were tested against it, with the parameters of "min_max_counts_per_cell = c(5e2, 6e6); cutoff = 0.1; min_cells_per_gene = 5" and other default parameters in inferCNV. Then, the estimated CNV values were re-scaled to 0-2, with 1 as the normal copy, to compare among samples. The CNV cluster of each sample was deduced by the hierarchical clustering (ward.D2) of inferCNV matrix. The CNV level of each cell was also calculated as previously described (51). The estimated CNV values were re-standardized as −1 to 1, and the CNV level of each cell was then calculated as the quadratic sum of all the expressed genes.

Estimating CNVs in WES Data
Sample DNA was extracted and sequenced via Agilent SureSelect Human All Exon v6 and Illumina platform. Whole-exome sequencing reads were aligned to human reference genome (b37), using BWA (52), followed by marking of duplications via Picard (http://broadinstitute.github.io/picard/). CNVkit (53) was used to call CNVs from targeted regions of exons in each sample, following the default workflow, with the bin size of 1 kb. A flat reference was made to run CNVkit with each sample, and during segmentation, the "CBS" method was applied, with 1e−4 as the significance threshold and parameters of "-drop-lowcoverage -drop-outliers 3."

Statistical Analysis
CNV score, gene expression, and undifferentiated score comparison between subclones was analyzed by D'Agostino and Pearson normality test, Shapiro-Wilk test, and Mann-Whitney test in GraphPad v8.3.0. Asymptotic two-sample Fisher-Pitman permutation test of cell-type composition between subclones of the sample GTE009 and between recurrent and primary samples was performed by "oneway_test" function of coin package in R v3.6.3. Correlation analysis was performed by R package "corrplot" from Taiyun Wei and Viliam Simko (2021) [Visualization of a Correlation Matrix (Version 0.90)] and cor() function in R v3.6.3.
The survival analysis was performed on published data (5, 6) by R package "survminer" from Alboukadel Kassambara, Marcin Kosinski, and Przemyslaw Biecek (2021) [Drawing Survival Curves using "ggplot2" (0.4.9)] and R package "survival" (54). The groups were separated by the relapse situation (recurrent or primary annotated by original authors (5, 6)) or the mean of trajectory scores of samples. First, a survival object was created by "Surv" function; second, the survival curves were created by "survfit" function based on a tabulation of the number at risk and at death time of events from the supplementary files of these published studies (5,6); and third, "ggsurvplot" function was used for the visualization of these curves.

Cell-Type Annotation
For cell malignancy analysis, we combined the following approaches to achieve a combinational separation of nonmalignant cells and malignant tumor cells. First, all cells were sorted by t-distributed stochastic neighbor embedding (t-SNE) projections in each patient and colored based on the cell clusters identified by Seurat (46), as malignant cells were often comprised of multiple clusters and were contiguous in t-SNE projection (55). Second, the malignancy of cells was explored by copy number variation (CNV) scores through whole-exon sequence by CNVkit (53) (Supplementary Figure S2) and through modified inferCNV (https://github.com/broadinstitute/ inferCNV) (30, 48) ( Figure 1A) for each sample. Third, the malignancy was further supported by the high undifferentiated score calculated by CytoTRACE (49) ( Figure 1B). Combining the result of cell-cycle stages (Supplementary Figure S1D), we classified each sample into malignant tumor cells and nonmalignant cells. For better exploration of the intratumoral genome, the pattern revealed by inferCNV was used for separation of malignant tumor cells into different subclones of each sample (Figures 1C, D). To support the genetic information inferred by transcriptome, we applied whole-exome sequencing analyzed by CNVkit (53) (Supplementary Figure S2) and obtained similar results in chromosomal level comparing to that from inferCNV. Fourth, the high correlation between malignant tumor cells and between non-malignant cells among samples supported malignancy separation (Supplementary Figure S1E). As validation, we performed KEGG analysis (50) on the differentially expressed genes (DEGs) of malignant tumor cells compared to non-malignant cells of merged samples, which showed enrichment on cell cycles and cancer-related terms: breast cancer, hepatocellular carcinoma, and proteoglycans in cancer (Supplementary Figure S1F). To sum up, we obtained non-malignant cells and malignant tumor cells and explored the genome of each sample.
To achieve cell-type classification, we applied signature enrichment analysis (Supplementary Figure S3). First, the DEGs of various cell types were calculated from each published data of the corresponding area (Supplementary Table S2) and were used as signatures to distinguish cell types in our tumor samples (19)(20)(21)(22)(23)(24) (Supplementary Figure S5F), which included RGC, AS, EC, EpC, NEU, NSC, OD, OPC, Mic, and T cells, from human and rodent embryonic and postnatal cortex scRNA-seq data. The enrichment of these gene signatures was calculated using the "AddModuleScore" function by subtracting the aggregated expression of control genes from the average expression levels of gene signatures (30). Then, the signature enrichment of those gene sets in cellular level was summarized in cluster level on average and was utilized as the standard to ascertain the cell type. Conclusively, the cell types of clusters were determined by the highest signature score.
Studies on scRNA-seq have also employed a reversed approach for cell-type classification (rSE; extracting signatures from unknown clusters and enriching them on published data of identified cell types; for details, see Materials and Methods; Supplementary Figure S3) (5,6,56,57). We then compared the similar results of these two methods in our data and obtained a high degree of correlation (Supplementary Figure S3). Hence, the following research applied the result of the first mentioned method, namely, SE analysis. This method and calculated signatures (Supplementary Table S2) were further tested on other published data, which displayed high similarity in cell-type classification compared to the cell types determined by original authors: human cortex (58) (Supplementary Figure S5B

Trajectory Analysis
Unbiased visualization of cell types without cells stack on each other was performed by SCUBI (12). For developmental trajectory analysis, the BAM files from Cell Ranger were processed by Velocyto (39) to obtain loom files containing spliced and unspliced transcript counts, which was used as input for scVelo (40) and Velocyto (39). Pseudo-time reconstruction and evaluation was performed by TSCAN (41) and by Slingshot (42) with modified default settings. Monocle (59,60) was applied on 2,000 variable features detected within more than 5% of cells to obtain reduced coordinates by "DDRTree." The trajectory score was inferred by "AddModuleScore." The undifferentiated trajectory score was calculated based on the significantly highly expressed genes in NSC-like cells in the subclone 2 of the sample GTE009, while the differentiated trajectory score was calculated based on the significantly highly expressed genes in EpC-like cells in the subclone 1 of the sample GTE009. The final trajectory score was then calculated by subtracting the differentiated trajectory score from the undifferentiated trajectory score of each cell.

Crosstalk and Gene Regulatory Network Analysis
Crosstalk analysis on EPN was performed on the integrated dataset of four EPN samples using CellChat (61), and the simulated 3D spatial structure of different cell types was calculated by CSOMAP (34). In details, CellChat preprocessed the expression data of our integrated four ependymomas for cellcell communication analysis; the cell-cell communication network was inferred by computation of the communication probability at a signaling pathway level and the calculation of the aggregated data frame. The elevated crosstalk pathways were validated in previously published EPN single-cell datasets (5,6).
Similarly, CSOMAP (34) computed the network of ligandreceptor interaction and then calculated optimized 3D coordinates.
Regulons for individual cell types were computed using the single-cell regulatory network inference and clustering (SCENIC) pipeline (62) on our integrated four EPN samples and validated by previously published EPN single-cell datasets (5,6). A log-normalized expression matrix of the four integrated EPN samples was used as an input into the pySCENIC workflow with default settings to infer regulon activity scores. To examine relevant networks using cell-cell communication analysis, we identified genes involved in crosstalk (ligands or receptors expressed in NSC-like cells) or gene regulatory networks (regulons with significantly high activity in NSC-like cells) of interest and then used genes classified under the same enriched terms in GO/KEGG analysis. Genes that had significantly higher expression in recurrent EPN than that in primary EPN were enriched in GO and KEGG analysis to highlight key terms in crosstalk and regulatory networks.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI-GEO with accession GSE189939 [https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE189939].

ETHICS STATEMENT
Fresh EPN tissue was excised by physicians with signed informed consent documents as approved by the ethics committee of Beijing Tiantan Hospital of Capital Medical University. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
Q-FW, WJ, SW, and HW conceived the project. JX and HW carried out the single-cell transcriptome amplification and library construction. YT and ZL did the sample collection. HW, RF, Y-HZ, and Z-HC performed the bioinformatics analysis. HW, SW, and Q-FW interpreted the results and wrote the manuscript with help from all the authors. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to express our gratitude to patients who donated samples for research purposes.  Table 2) and used as signatures to distinguish cell types. All unknown cells are then clustered at high resolution to obtain multiple clusters, and the highest signature enrichment score in each cluster is designated as the cell type identity for these clusters. Bottom: Schematic for cell type classification by reversed signature enrichment (rSE). In the reciprocal analysis pipeline, unknown cells are first clustered at high resolution to obtain multiple clusters and the DEGs of all clusters are calculated and used as signatures to distinguish cell types. The signatures are then compared with published transcriptomic datasets of human and rodent embryonic and adult cortex (details see Supplementary Table 2), and the highest signature enrichment score is assigned as the name for the unknown cluster. The correlation result of SE-and rSE-determined cell types by 'cor' function of stats package in R v3.6.3 confirms high correlation across the two analysis pipelines.