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 YAP1-fusions, and spinal EPNs (1, 4). With the advent of high-throughput 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 ependymal-like 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 non-small-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 RNA-seq 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.

Results

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 differentiation-associated genes, which revealed that CNV-inferred neoplastic cell populations were in a less differentiated state (Figure 1B).

FIGURE 1

Figure 1 scRNA-seq analysis reveals intratumoral subclone heterogeneity in PF-EPN. (A) CNV score calculated by modified inferCNV of PF-EPN sample GTE009 presented on tSNE reduction. (B) Undifferentiated score calculated by CytoTRACE of PF-EPN sample GTE009 presented on tSNE reduction. (C) CNV heatmap (rows represent cells, and columns represent CNV score of genes) of malignant tumor cells from four EPN samples labeled by genetic subclone information for each sample. (D) Subclonal populations in malignant cells and NM cells of PF-EPN sample GTE009 classified by CNV pattern presented on tSNE reduction. (E) tSNE plot of all clusters in PF-EPN sample GTE009 color coded by cell types with unbiased visualization by SCUBI (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.

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 scRNA-seq data (Supplementary Table S2) and were enriched as signatures to classify cell types in our tumor samples (19–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 EpC-like 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 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).

FIGURE 2

Figure 2 Highly differentiated cells in PF-EPN subclonal populations show CNV amplification and enrichment of cilium-associated genes. (A) Schematic of RNA splicing analysis and cell differentiation using RNA velocity and trajectory deduction methodologies. (B) RNA velocity inferred by Velocyto and scVelo of malignant tumor cells presented on tSNE reduction and colored by cell types in PF-EPN sample GTE009. (C) Differentiation trajectory inferred by Monocle of malignant tumor cells in PF-EPN sample GTE009. (D) Volcano plot showing genes with differentially expressed CNV values highlighting DYNC2H1 in the PF-EPN sample GTE009 subclone 1 compared to the subclone 2. (E) Heatmap of chromosome 11 showing inferCNV scores colored by cell types designated in Figure 2B and subclone annotation in PF-EPN sample GTE009. DYNC2H1 is highlighted by black vertical bar. (F) Violin plot showing significant difference (p < 0.0001; Mann Whitney–test) in gene expression (RNA) and CNV level of DYNC2H1 between subclones in PF-EPN sample GTE009. (G) Heatmap showing relative expression of identified genes from cilium-related terms in GO analysis (27) colored by subclones in the PF-EPN sample GTE009. (H) Correlation analysis of undifferentiated score in EpC-like cells, normalized average expression of markers in EpC-like cells in subclone 1 (EpC-Sub1) and subclone 2 (EpC-Sub2), and normalized average expression of Mic (Control; see 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.

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.

FIGURE 3

Figure 3 Trajectory score analysis can predict EPN cell compositions associated with poor survival outcomes. (A) Gene ontology analysis of differentially expressed genes from EpC- and NSC-like cells between subclones classified by annotation of subclones and cell types in PF-EPN sample GTE009. (B) Workflow for calculating trajectory score based on published computation method (30). E, expression; TPM_i,j, transcript per million (TPM) for gene i in sample j; Er, relative expression. (C) tSNE plot of trajectory score in combined subclone 1 and 2 datasets with EpC- and NSC-like cell populations labeled in the PF-EPN sample GTE009. (D) Validation of trajectory score on published scRNA-seq data (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.

To provide a numerical representation of molecular trajectory information at the single-cell level, we developed a trajectory score for downstream analyses (Figure 3B): the normalized average expression of undifferentiated NSC (NSC-like cells in 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 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).

FIGURE 4

Figure 4 Cellular populations in recurrent EPN with poor prognosis are associated with higher trajectory score. (A) Survival plot of primary and recurrent EPN patients. The solid line refers to overall survival (OS; p-value = 0.0018), and the dotted line refers to progression-free survival (PFS; p = 0.00026), which are colored by relapse situations on published scRNA-seq data (5, 6). (B) Histogram of cell types in primary and recurrent EPN colored by cell types and outlined by primary/recurrent conditions showing significant difference (p < 2.2e−16) between cell types using asymptotic two-sample Fisher–Pitman permutation test. (C) Trajectory score analysis comparison between primary and recurrent samples in NSC- and EpC-like cells using Kruskal–Wallis test (all p < 0.0001). (D) Gene ontology analysis of differentially expressed genes in NSC-like cells between primary and relapse conditions. The "*" means significant difference (p-value): * p < 0.05; ** P < 0.01; *** p < 0.001; **** p < 0.0001.

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) and identified the presence of numerous crosstalk events (Figure 5B; Supplementary Figure S8A; Supplementary 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). Taken together, crosstalk analysis on 35,102 individual cells in conjunction with validation using 36 EPN patients revealed elevated cell–cell interactions between malignant cells and tumor-infiltrating NM cells, such as between NSC-like cells and Mic, consistent with studies demonstrating a key role of the central nervous system TME in the pathogenesis of EPN.

FIGURE 5

Figure 5 Crosstalk analysis reveals cell–cell interactions implicating immune cell populations in recurrent EPN. (A) Cell numbers of all cell types in four EPN samples. (B) Crosstalk net analyzed by CellChat. Individual lines represent the crosstalk from source to target cells, highlighting interactions from NSC-like cells and Mic to other cell types. Related to Supplementary Figures S8A, B. (C) Simulated 2D spatial structure showing overlap of Mic and NSC-like cell populations by CSOMAP (34). Related to Supplementary Figure S8C. (D) Heatmap of ligands or receptors with significantly higher expression in recurrent samples compared to primary samples, colored by cell type and gene class (ligands or receptors) using published single-cell transcriptomes of 36 EPN samples (5, 6). (E) Expression of MDK (ligand) and NCL (receptor) between recurrent and primary samples of 36 EPN patients (5, 6) (p < 0.0001; Mann–Whitney test).The "*" means significant difference (p-value): * p < 0.05; ** P < 0.01; *** p < 0.001; **** p < 0.0001.

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 single-cell 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 single-transcriptome 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 cell-type 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.

Materials and Methods

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 physicians, with signed informed consent documents and was approved by the Ethics Committee of Beijing Tiantan Hospital of Capital Medical University. Samples were delivered on ice to the Institute of Genetics and Developmental Biology of Chinese Academy of Sciences immediately. Microdissected tissues were transferred to a 24-well cell culture plate and digested by buffer comprising 20 U/ml Papain (LK003178; Worthington, Lakewood, USA), 100 U/ml DNaseI (LK003172; Worthington), 10 U/L chondroitinase ABC (C3667; Sigma-Aldrich, St. Louis, USA), 0.07% hyaluronidase (R006687; Rhawn, Canton, China), 1× Glutamax (35050061; Life Technologies, Waltham, USA), 0.05 mM (2R)-amino-5-phosphonovaleric acid (APV; 010510; Thermo Fisher Scientific, Waltham, USA), 0.01 mM Y27632 dihydrochloride (T9531; Sigma), and 0.2× B27 supplement (17504044; Thermo Fisher Biosciences) in Hibernate-E media (A1247601, Life Technologies) for 1–2 h at 37°C, and then pooled with Hibernate-E buffer containing 1× Glutamax, 0.05 mM APV, 0.2× B27, and 0.01 mM Y27632 dihydrochloride. Tissues were gently triturated through Pasteur pipettes with finely polished tips of 600, 300, and 200 μm diameters in order, and washed once with Hibernate-E buffer to generate single-cell suspensions. After filtration through a 40-μm strainer (130-101-812; Thermo Fisher Scientific), 1× red blood cell lysis solution (130-094-183; Miltenyi Biotec, Bergisch Gladbach, Germany) was added to remove blood contamination during surgery followed by 1,800 μl debris removal solution (130-109-398; Miltenyi Biotec). Subsequently, the dissociated cells were stained with 4′,6-diamidino-2-phenylindole (DAPI) (0.2 µg/ml) to identify dead cells. scRNA-seq libraries were constructed following the manufacturer’s instructions provided by 10× Genomics accompanying single-cell 3′ Library and Gel Bead Kit V3 (1000075; 10× Genomics, Pleasanton, USA). A Chromium Single Cell Controller (10× Genomics) was loaded in cell suspensions (300–600 living cells per microliter determined by Count Star) to generate single-cell gel beads in the emulsion (GEM). Quality control was performed on the generated cDNA library using Agilent 4200, and scRNA-seq was performed on a Illumina Novaseq6000 sequencer.

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 whole-genome 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 non-immune 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-low-coverage –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 non-malignant 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 non-malignant 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–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), rodent cortex (24) (Supplementary Figure S5C), ependymoma (5) (Supplementary Figure S5D), and childhood ependymoma (6) (Supplementary Figure S5E).

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 cell–cell 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 ligand–receptor 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.

Funding

The work was supported by the National Key R&D Program of China (2019YFA0801900 and 2018YFA0801104), the National Natural Science Foundation of China (81891002, 31921002, 32070972, and 31771131), the Strategic Priority Research Program of Chinese Academy of Sciences (XDB32020000), and the Beijing Municipal Science & Technology Commission (Z210010 and Z181100001518001).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to express our gratitude to patients who donated samples for research purposes.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.903246/full#supplementary-material

Supplementary Figure 1 | Quality Control of scRNA-Seq Analysis of Human EPN. (A) Doublets identified by Doublet Finder presented on tSNE reduction. High doublet cells are filtered out. (B) Density plot of mitochondria percentage of cells in sample GTE009. Cells with more than 5% of mitochondria percentage are filtered out. (C) Density plot of captured gene numbers of cells in sample GTE009. Cells with less than 1500 transcripts are filtered. (D) Cell cycle phases of cells in sample GTE009 re-calculated after quality control presented on tSNE reduction. (E) Correlation plot of transcriptome of malignant tumor cells (Mal) and non-malignant cells (NM) among samples (GTE001, GTE002, GTE009, and GTE012). (F) KEGG analysis on DEGs of malignant tumor cells compared to non-malignant cells of four merged samples (GTE001, GTE002, GTE009, and GTE012). (G) Enrichment of RGC signatures in malignant tumor cells compared to other cells types using single-cell transcriptomes from 1) this study, 2) EPN (5) and 3) childhood EPN (6) (one-way ANOVA analysis; p value < 0.0001).

Supplementary Figure 2 | CNV Analysis of Whole-exome Sequencing. CNV heatmap of whole-exon sequence data labeled by samples.

Supplementary Figure 3 | Workflow of Cell Type Classification. Top: Schematic for cell type classification by signature enrichment (SE). The DEGs of different cell populations are obtained from published transcriptomic datasets of human and rodent embryonic and adult cortex (details see Supplementary 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.

Supplementary Figure 4 | Additional scRNA-Seq Analysis of Human EPN samples. (A) CNV score calculated by modified inferCNV of samples (GTE001, GTE002, and GTE012) presented on tSNE reduction. (B) Cell cycle phases in cellular level of samples (GTE001, GTE002, and GTE012) presented on tSNE reduction. (C) Undifferentiated score calculated by CytoTRACE of samples (GTE001, GTE002, and GTE012) presented on tSNE reduction. (D) Classified non-malignant cells and malignant tumor cells of samples (GTE001, GTE002, and GTE012) presented on tSNE reduction. (E) Annotated clusters of samples (GTE001, GTE002, and GTE012) presented on tSNE reduction with unbiased visualization by SCUBI (12). (F) Enrichment of signatures in malignant tumor cells compared to other cell types in samples (GTE001, GTE002, and GTE012; one-way ANOVA analysis; p value < 0.0001). See also Supplementary Table 1.

Supplementary Figure 5 | EPN cell type classification validation. (A) Heatmap of DEGs in annotated clusters of the sample GTE009. (B–E), Correlation of cell types classified by method of signature enrichment (rows) and that of original cell types determined by original authors [B: human cortex (58), C: rodent cortex (63), D: ependymoma (5), and E: childhood ependymoma (6)]. (F) Workflow of cell-type classification. Signature markers genes are obtained from public transcriptome databases of human and rodent cortex and used for cell-type assignment and cluster annotation.

Supplementary Figure 6 | Trajectory and pseudotime analysis. The columns are samples from patients (from left to right: GTE001, GTE002, GTE009, and GTE012). (A) RNA velocity inferred by Velocyto and scVelo of malignant tumor cells presented on tSNE reduction and colored by cell types. (B, C) Pseudo-time reconstruction and evaluation by TSCAN (41) (B) and by Slingshot (42) (C). (D) Differentiation trajectory inferred by Monocle of malignant tumor cells.

Supplementary Figure 7 | Additional trajectory score analysis of EPN. (A) Gene ontology analysis of upregulated genes in patients with low trajectory score compared to patients with high trajectory score in published ependymoma scRNA-seq data (5, 6). (B) Undifferentiated score analysis on published scRNA-seq datasets of ependymoma (5, 6) (Kruskal-Wallis test). (C) Overall survival analysis on the trajectory score. A p value of 0.39 and 0.72 were obtained from the difference between two groups of trajectory score of patients (separated by the mean of trajectory score), with a trend of worse survival in patients with high trajectory score. The high and low group were separated by the mean of trajectory score on published scRNA-seq data (5, 6). (D) Histogram showing percentage of cells with high and low trajectory score and outlined by subclone annotation in samples from primary and recurrent patients. Permutation test shown significant compositional difference between primary and recurrent samples (p value = 0.01241; asymptotic two-sample Fisher-Pitman permutation test). *p < 0.05; **P < 0.01; ***p < 0.001; ****p < 0.0001

Supplementary Figure 8 | Additional crosstalk analysis of EPN. (A) Interaction strength between cell types profiled in EPN samples inferred by CellChat. (B) Crosstalk net analyzed by CellChat. Individual lines represent the crosstalk from source to target cells. (C) Simulated 3D spatial structure of all cells and 2D angle of the simulated spatial structure by CSOMAP (34) of Mic, NSC-like cells, and NEU cells respectively, colored by pink (cell type of interest) and blue (other cell types). (D, E) Circle plots of ligands and receptors with higher expression in recurrent samples than that in primary samples. Lines represent the crosstalk between specific ligands and colors represent the cell type origin for each interaction

Supplementary Figure 9 | Additional regulon analysis of EPN. (A) Gene regulatory networks were inferred by SCENIC and were clustered by cell types (bottom) and regulons (right). (B) List of significantly upregulated genes in recurrent samples from crosstalk and gene regulatory network analysis of NSC-like cells which share the same enriched terms in GO/KEGG analysis. (C, D) Visualization of genes using GO and KEGG enrichment analysis.

Supplementary Table 1 | Clinical data from EPN patients. Clinical data of four sequenced EPN samples from patients in this study.

Supplementary Table 2 | Signatures for cell type classification. For recognition of cell type in SE algorithms, Seurat was used to calculate the DEGs (latter utilized as signatures) of each cell type from each reference data 12-16,62 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.

Supplementary Table 3 | Differentially expressed genes between subclones of the sample GTE009. Differentially expressed genes calculated from subclone 1 compared to the subclone 2 of the sample GTE009.

Supplementary Table 4 | Analysis on cilium-related genes and CNV score. DEGs calculated from subclone 1 compared to the subclone 2 of the sample GTE009 marked by GO terms, annotated with results of statistical analysis from CNV score.

Supplementary Table 5 | Crosstalk analysis output. List of ligand-receptor interactions annotated by cell type source and target.

References

1. Pajtler KW, Witt H, Sill M, Jones DT, Hovestadt V, Kratochwil F, et al. Molecular Classification of Ependymal Tumors Across All CNS Compartments, Histopathological Grades, and Age Groups. Cancer Cell (2015) 27(5):728–43. doi: 10.1016/j.ccell.2015.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ramaswamy V, Hielscher T, Mack SC, Lassaletta A, Lin T, Pajtler KW, et al. Therapeutic Impact of Cytoreductive Surgery and Irradiation of Posterior Fossa Ependymoma in the Molecular Era: A Retrospective Multicohort Analysis. J Clin Oncol (2016) 34(21):2468–77. doi: 10.1200/jco.2015.65.7825

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Merchant TE, Li C, Xiong X, Kun LE, Boop FA, Sanford RA. Conformal Radiotherapy After Surgery for Paediatric Ependymoma: A Prospective Study. Lancet Oncol (2009) 10(3):258–66. doi: 10.1016/S1470-2045(08)70342-5

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Panwalkar P, Clark J, Ramaswamy V, Hawes D, Yang F, Dunham C, et al. Immunohistochemical Analysis of H3K27me3 Demonstrates Global Reduction in Group-A Childhood Posterior Fossa Ependymoma and is a Powerful Predictor of Outcome. Acta Neuropathol (2017) 134(5):705–14. doi: 10.1007/s00401-017-1752-4

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gojo J, Englinger B, Jiang L, Hübner JM, Shaw ML, Hack OA, et al. Single-Cell RNA-Seq Reveals Cellular Hierarchies and Impaired Developmental Trajectories in Pediatric Ependymoma. Cancer Cell (2020) 38(1):44–59.e9. doi: 10.1016/j.ccell.2020.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gillen AE, Riemondy KA, Amani V, Griesinger AM, Gilani A, Venkataraman S, et al. Single-Cell RNA Sequencing of Childhood Ependymoma Reveals Neoplastic Cell Subpopulations That Impact Molecular Classification and Etiology. Cell Rep (2020) 32(6):108023. doi: 10.1016/j.celrep.2020.108023

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Miller CA, Dahiya S, Li T, Fulton RS, Smyth MD, Dunn GP, et al. Resistance-Promoting Effects of Ependymoma Treatment Revealed Through Genomic Analysis of Multiple Recurrences in a Single Patient. Cold Spring Harb Mol Case Stud (2018) 4(2):a002444. doi: 10.1101/mcs.a002444

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Jamal-Hanjani M, Wilson GA, McGranahan N, Birkbak NJ, Watkins TBK, Veeriah S, et al. Tracking the Evolution of Non-Small-Cell Lung Cancer. N Engl J Med (2017) 376(22):2109–21. doi: 10.1056/NEJMoa1616288

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell (2019) 178(4):835–49.e821. doi: 10.1016/j.cell.2019.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Reinartz R, Wang S, Kebir S, Silver DJ, Wieland A, Zheng T, et al. Functional Subclone Profiling for Prediction of Treatment-Induced Intratumor Population Shifts and Discovery of Rational Drug Combinations in Human Glioblastoma. Clin Cancer Res (2017) 23(2):562–74. doi: 10.1158/1078-0432.Ccr-15-2089

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hou W, Ji Z. Unbiased Visualization of Single-Cell Genomic Data With SCUBI. Cell Rep Methods (2022) 2(1):100135. doi: 10.1016/j.crmeth.2021.100135

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI Gene Expression and Hybridization Array Data Repository. Nucleic Acids Res (2002) 30:207–10. doi: 10.1093/nar/30.1.207

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: Archive for Functional Genomics Data Sets–Update. Nucleic Acids Res (2013) 41(Database issue):D991–5. doi: 10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bayliss J, Mukherjee P, Lu C, Jain SU, Chung C, Martinez D, et al. Lowered H3K27me3 and DNA Hypomethylation Define Poorly Prognostic Pediatric Posterior Fossa Ependymomas. Sci Trans Med (2016) 8(366):366ra161. doi: 10.1126/scitranslmed.aah6904

CrossRef Full Text | Google Scholar

16. Pajtler KW, Wen J, Sill M, Lin T, Orisme W, Tang B, et al. Molecular Heterogeneity and CXorf67 Alterations in Posterior Fossa Group A (PFA) Ependymomas. Acta Neuropathol (2018) 136(2):211–26. doi: 10.1007/s00401-018-1877-0

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Vladoiu MC, El-Hamamy I, Donovan LK, Farooq H, Holgado BL, Sundaravadanam Y, et al. Childhood Cerebellar Tumours Mirror Conserved Fetal Transcriptional Programs. Nature (2019) 572(7767):67–73. doi: 10.1038/s41586-019-1158-7

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Love MI, Huber W, Anders S. Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Polioudakis D, de la Torre-Ubieta L, Langerman J, Elkins AG, Shi X, Stein JL, et al. A Single-Cell Transcriptomic Atlas of Human Neocortical Development During Mid-Gestation. Neuron (2019) 103(5):785–801.e788. doi: 10.1016/j.neuron.2019.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Oetjen KA, Lindblad KE, Goswami M, Gui G, Dagur PK, Lai C, et al. Human Bone Marrow Assessment by Single-Cell RNA Sequencing, Mass Cytometry, and Flow Cytometry. JCI Insight (2018) 3(23):e124928. doi: 10.1172/jci.insight.124928

CrossRef Full Text | Google Scholar

21. Velmeshev D, Schirmer L, Jung D, Haeussler M, Perez Y, Mayer S, et al. Single-Cell Genomics Identifies Cell Type-Specific Molecular Changes in Autism. Science (2019) 364(6441):685–9. doi: 10.1126/science.aav8130

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zheng GX, Terry JM, Belgrader P, Ryvkin P, Bent ZW, Wilson R, et al. Massively Parallel Digital Transcriptional Profiling of Single Cells. Nat Commun (2017) 8:14049. doi: 10.1038/ncomms14049

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Franzén O, Gan LM, Björkegren JLM. PanglaoDB: A Web Server for Exploration of Mouse and Human Single-Cell RNA Sequencing Data. Database (2019) 2019:baz046. doi: 10.1093/database/baz046

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zeisel A, Hochgerner H, Lönnerberg P, Johnsson A, Memic F, van derZwan J, et al. Molecular Architecture of the Mouse Nervous System. Cell (2018) 174(4):999–1014.e22. doi: 10.1016/j.cell.2018.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Liu H, Kiseleva AA, Golemis EA. Ciliary Signalling in Cancer. Nat Rev Cancer (2018) 18:511–24. doi: 10.1038/s41568-018-0023-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhang J, Chandrasekaran G, Li W, Kim D-Y, Jeong IY, Lee SH, et al. Wnt-PLC-IP-Connexin-Ca Axis Maintains Ependymal Motile Cilia in Zebrafish Spinal Cord. Nat Commun (2020) 11(1):1860. doi: 10.1038/s41467-020-15248-2

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: Tool for the Unification of Biology. The Gene Ontology Consortium. Nat Genet (2000) 25(1):25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Fujita A, Higashijima T, Shirozu H, Masuda H, Sonoda M, Tohyama J, et al. Pathogenic Variants of DYNC2H1, KIAA0556, and PTPN11 Associated With Hypothalamic Hamartoma. Neurology (2019) 93:e237–51. doi: 10.1212/WNL.0000000000007774

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Venkatesan AM, Vyas R, Gramann AK, Dresser K, Gujja S, Bhatnagar S, et al. Ligand-Activated BMP Signaling Inhibits Cell Differentiation and Death to Promote Melanoma. J Clin Invest (2018) 128:294–308. doi: 10.1172/jci92513

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, et al. Dissecting the Multicellular Ecosystem of Metastatic Melanoma by Single-Cell RNA-Seq. Science (2016) 352:189–96. doi: 10.1126/science.aad0501

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Shah PT, Stratton JA, Stykel MG, Abbasi S, Sharma S, Mayr KA, et al. Single-Cell Transcriptomics and Fate Mapping of Ependymal Cells Reveals an Absence of Neural Stem Cell Function. Cell (2018) 173:1045–57.e1049. doi: 10.1016/j.cell.2018.03.063

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hambardzumyan D, Gutmann DH, Kettenmann H. The Role of Microglia and Macrophages in Glioma Maintenance and Progression. Nat Neurosci (2016) 19:20–7. doi: 10.1038/nn.4185

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Wesolowska A, Kwiatkowska A, Slomnicki L, Dembinski M, Master A, Sliwa M, et al. Microglia-Derived TGF-Beta as an Important Regulator of Glioblastoma Invasion–an Inhibition of TGF-Beta-Dependent Effects by shRNA Against Human TGF-Beta Type II Receptor. Oncogene (2008) 27:918–30. doi: 10.1038/sj.onc.1210683

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ren X, Zhong G, Zhang Q, Zhang L, Sun Y, Zhang Z. Reconstruction of Cell Spatial Organization From Single-Cell RNA Sequencing Data Based on Ligand-Receptor Mediated Self-Assembly. Cell Res (2020) 30:763–78. doi: 10.1038/s41422-020-0353-2

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Coniglio SJ, Eugenin E, Dobrenis K, Stanley ER, West BL, Symons MH, et al. Microglial Stimulation of Glioblastoma Invasion Involves Epidermal Growth Factor Receptor (EGFR) and Colony Stimulating Factor 1 Receptor (CSF-1R) Signaling. Mol Med (2012) 18:519–27. doi: 10.2119/molmed.2011.00217

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Jia W, Yao Z, Zhao J, Guan Q, Gao L. New Perspectives of Physiological and Pathological Functions of Nucleolin (NCL). Life Sci (2017) 186:1–10. doi: 10.1016/j.lfs.2017.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Filippou PS, Karagiannis GS, Constantinidou A. Midkine (MDK) Growth Factor: A Key Player in Cancer Progression and a Promising Therapeutic Target. Oncogene (2020) 39:2040–54. doi: 10.1038/s41388-019-1124-8

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Takada S, Sakakima H, Matsuyama T, Otsuka S, Nakanishi K, Norimatsu K, et al. Disruption of Midkine Gene Reduces Traumatic Brain Injury Through the Modulation of Neuroinflammation. J Neuroinflamm (2020) 17:40. doi: 10.1186/s12974-020-1709-8

CrossRef Full Text | Google Scholar

39. La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, et al. RNA Velocity of Single Cells. Nature (2018) 560:494–8. doi: 10.1038/s41586-018-0414-6

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA Velocity to Transient Cell States Through Dynamical Modeling. Nat Biotechnol (2020) 38:1408–14. doi: 10.1038/s41587-020-0591-3

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ji Z, Ji H. TSCAN: Pseudo-Time Reconstruction and Evaluation in Single-Cell RNA-Seq Analysis. Nucleic Acids Res (2016) 44:e117. doi: 10.1093/nar/gkw430

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: Cell Lineage and Pseudotime Inference for Single-Cell Transcriptomics. BMC Genomics (2018) 19:477. doi: 10.1186/s12864-018-4772-0

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Quail DF, Joyce JA. The Microenvironmental Landscape of Brain Tumors. Cancer Cell (2017) 31:326–41. doi: 10.1016/j.ccell.2017.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Griesinger AM, Josephson RJ, Donson AM, Mulcahy Levy JM, Amani V, Birks DK, et al. Interleukin-6/STAT3 Pathway Signaling Drives an Inflammatory Phenotype in Group A Ependymoma. Cancer Immunol Res (2015) 3:1165–74. doi: 10.1158/2326-6066.CIR-15-0061

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zhang YH, Xu M, Shi X, Sun XL, Mu W, Wu H, et al. Cascade Diversification Directs Generation of Neuronal Diversity in the Hypothalamus. Cell Stem Cell (2021) 28(8):1483–1499.e8. doi: 10.1016/j.stem.2021.03.020

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, et al. Comprehensive Integration of Single-Cell Data. Cell (2019) 177:1888–902.e1821. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

47. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst (2019) 8:329–37.e324. doi: 10.1016/j.cels.2019.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, et al. Single-Cell RNA-Seq Highlights Intratumoral Heterogeneity in Primary Glioblastoma. Science (2014) 344:1396–401. doi: 10.1126/science.1254257

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Gulati GS, Sikandar SS, Wesche DJ, Manjunath A, Bharadwaj A, Berger MJ, et al. Single-Cell Transcriptional Diversity is a Hallmark of Developmental Potential. Science (2020) 367:405–11. doi: 10.1126/science.aax0249

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Peng J, Sun BF, Chen CY, Zhou JY, Chen YS, Chen H, et al. Single-Cell RNA-Seq Highlights Intra-Tumoral Heterogeneity and Malignant Progression in Pancreatic Ductal Adenocarcinoma. Cell Res (2019) 29:725–38. doi: 10.1038/s41422-019-0195-y

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Li H, Durbin R. Fast and Accurate Long-Read Alignment With Burrows-Wheeler Transform. Bioinformatics (2010) 26:589–95. doi: 10.1093/bioinformatics/btp698

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Talevich E, Shain AH, Botton T, Bastian BC. CNVkit: Genome-Wide Copy Number Detection and Visualization From Targeted DNA Sequencing. PLoS Comput Biol (2016) 12:e1004873. doi: 10.1371/journal.pcbi.1004873

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. New York: Springer (2000).

Google Scholar

55. Yuan J, Levitin HM, Frattini V, Bush EC, Boyett DM, Samanamud J, et al. Single-Cell Transcriptome Analysis of Lineage Diversity in High-Grade Glioma. Genome Med (2018) 10:57. doi: 10.1186/s13073-018-0567-9

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, et al. Single-Cell RNA-Seq Supports a Developmental Hierarchy in Human Oligodendroglioma. Nature (2016) 539:309–13. doi: 10.1038/nature20123

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Filbin MG, Tirosh I, Hovestadt V, Shaw ML, Escalante LE, Mathewson ND, et al. Developmental and Oncogenic Programs in H3K27M Gliomas Dissected by Single-Cell RNA-Seq. Science (2018) 360:331–5. doi: 10.1126/science.aao4750

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Bakken TE, Jorstad NL, Hu Q, Lake BB, Tian W, Kalmbach BE, et al. Comparative Cellular Analysis of Motor Cortex in Human, Marmoset and Mouse. Nature (2021) 598:111–9. doi: 10.1038/s41586-021-03465-8

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The Dynamics and Regulators of Cell Fate Decisions are Revealed by Pseudotemporal Ordering of Single Cells. Nat Biotechnol (2014) 32:381–6. doi: 10.1038/nbt.2859

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed Graph Embedding Resolves Complex Single-Cell Trajectories. Nat Methods (2017) 14:979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan C-H, et al. Inference and Analysis of Cell-Cell Communication Using CellChat. Nat Commun (2021) 12:1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: Single-Cell Regulatory Network Inference and Clustering. Nat Methods (2017) 14:1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Zhang J, Chandrasekaran G, Li W, Kim DY, Jeong IY, Lee S-H, et al. Wnt-PLC-IP(3)-Connexin-Ca(2+) Axis Maintains Ependymal Motile Cilia in Zebrafish Spinal Cord. Nat Commun (2020) 11:1860. doi: 10.1038/s41467-020-15248-2

PubMed Abstract | CrossRef Full Text | Google Scholar