Transcriptional immunogenomic analysis reveals distinct immunological clusters in paediatric nervous system tumours

Abstract

Background

Cancer immunotherapies including immune checkpoint inhibitors and Chimeric Antigen Receptor (CAR) T-cell therapy have shown variable response rates in paediatric patients highlighting the need to establish robust biomarkers for patient selection. While the tumour microenvironment in adults has been widely studied to delineate determinants of immune response, the immune composition of paediatric solid tumours remains relatively uncharacterized calling for investigations to identify potential immune biomarkers.

Methods

To inform immunotherapy approaches in paediatric cancers with embryonal origin, we performed an immunogenomic analysis of RNA-seq data from 925 treatment-naïve paediatric nervous system tumours (pedNST) spanning 12 cancer types from three publicly available data sets.

Results

Within pedNST, we uncovered four broad immune clusters: Paediatric Inflamed (10%), Myeloid Predominant (30%), Immune Neutral (43%) and Immune Desert (17%). We validated these clusters using immunohistochemistry, methylation immune inference and segmentation analysis of tissue images. We report shared biology of these immune clusters within and across cancer types, and characterization of specific immune cell frequencies as well as T- and B-cell repertoires. We found no associations between immune infiltration levels and tumour mutational burden, although molecular cancer entities were enriched within specific immune clusters.

Conclusions

Given the heterogeneity of immune infiltration within pedNST, our findings suggest personalized immunogenomic profiling is needed to guide selection of immunotherapeutic strategies.

Background

Cancer immunotherapies have been clinically and experimentally investigated in paediatric oncology with a wide range of response rates. Objective responses to immune checkpoint inhibitors (ICI) as a monotherapy have been limited to 5–11% of paediatric cancers [1,2,3]. Addition of a monoclonal antibody targeting disialoganglioside GD2, granulocyte–macrophage colony-stimulating factor (GM-CSF) and interleukin-2 to standard therapy (isotretinoin) has been shown to improve overall survival of children with high-risk neuroblastoma treated with intensive multimodal therapy [4, 5]. More recently, anti-GD2 Chimeric Antigen Receptor (CAR) T-cell therapy has been investigated in diffuse midline gliomas in four patients, three of whom showed clinical improvement [6]. Experimentally, novel immunotherapies have been proposed and tested in preclinical models including CAR-T targeting immune checkpoint protein, B7-H3 [7], and targeting the myeloid compartment with anti-CSF1R [8] or anti-CD47 [9]. Considering the wide range of response rates for existing and emerging paediatric cancer immunotherapies, biomarkers for patient stratification are needed to identify potential candidates for clinical trials.

Immunogenomic analysis of tumours has been a major focus of biomarker discovery for immunotherapy. A prominent outcome of such studies is the FDA approval of tumour mutation burden (TMB) and microsatellite instability as the first immunotherapy-related biomarkers in adult cancers [10]. Several other biomarkers have been derived from large-scale genomic and transcriptomic datasets with a major focus on adult extracranial tumours [11,12,13,14,15]. Despite recent single-cell RNA sequencing (RNA-seq) studies in neuroblastoma [16], high-grade glioma [17], ependymoma [18, 19] and medulloblastoma [20], a comparison of the immune microenvironment across paediatric nervous system cancers and implications for informing immunotherapeutic interventions or patient selection have not been systematically analysed. In this study, we sought to characterize the immune microenvironment in 925 treatment-naïve paediatric central and peripheral nervous system cancers, as they share embryonal origins and similar levels of TMB [21]. We showed discrepancies across existing immune deconvolution tools when applied to paediatric nervous system tumours. To address this, we performed a consensus geneset analysis using publicly available datasets to identify specific immune genes with no tumour cell expression, followed by further characterization of pathway and gene expression, TMB and genomic alterations, T- and B-cell repertoire and cellular interactions.

Methods

Human subjects

We curated a total of 925 tumours from the Children’s Brain Tumour Network [22] (CBTN, n = 581), the Therapeutically Applicable Research To Generate Effective Treatments [23] (NCI TARGET, n = 149) and the International Cancer Genome Consortium [24] (ICGC, n = 195). We selected tumours from paediatric patients (median age 7 years) who had RNA-seq data generated from their primary tumours. We focused this study to major types of paediatric CNS tumours and neuroblastoma; ETMR (n = 9), neurofibroma (n = 11), choroid plexus tumours (n = 16), meningioma (n = 13), schwannoma (n = 14), craniopharyngioma (n = 27), ATRT (n = 31), ependymoma (n = 65), pedHGG (n = 83), neuroblastoma (n = 151), medulloblastoma (n = 208) and pedLGG (n = 298). As adult cancer comparator, we included TCGA participants with common types of adult cancers; glioblastoma multiforme (GBM, n = 153), low-grade glioma (LGG, n = 507), skin cutaneous melanoma (SKCM, n = 102), colorectal adenocarcinoma (COAD, n = 298), ovarian serous adenocarcinoma (OV, n = 373), prostate adenocarcinoma (PRAD, n = 497) and lung adenocarcinoma (LUAD, n = 522).

Tumour datasets

We used RNA-seq datasets from the TARGET [23] (https://ocg.cancer.gov/programs/target), CBTN [1] (https://cbtn.org/research/specimendata/), ICGC [24] (https://icgc.org/icgc/cgp/62/345/822) and TCGA [25] (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) in this study. For the ICGC dataset, we chose 195 primary CNS tumour samples with matched RNA-seq and WGS, and complete clinical data [21]. For the TARGET dataset, we chose 149 primary neuroblastoma samples with matched RNA-seq and WES data. For the CBTN dataset [22], we refined the original dataset (n = 996) and excluded data from 22 cell lines, 153 progression samples, 70 recurrences and 16 secondary malignancies. We further removed 98 samples from rare (defined as ≤ 10% of the cohort) or unannotated tumours. To enable comparison at the transcriptional level, we performed unsupervised clustering on primary samples and removed 51 samples that did not match their corresponding pathological annotation (Additional file 1: Fig. S1B-C). This sample curation resulted in 581 samples from the CBTN consortium (Additional file 1: Fig. S1A). These resulted in a cohort of 925 paediatric CNS tumours and neuroblastomas.

We accessed CBTN, TARGET and TCGA data through the Kids First Data Resource Centre [22] (https://portal.kidsfirstdrc.org/) and used the CAVATICA [22] platform for data processing and analysis (https://cavatica.sbgenomics.com/). In total, we analysed 925 cases of paediatric nervous system tumours (CNS and neural crest tumours) (CBTN, TARGET and ICGC). We included RNA-seq data from 79 PDX models (2 ATRT, 10 ependymoma, 10 HGG, 14 neuroblastoma, 43 medulloblastoma) obtained from ITCC-P4 as control for immune infiltration and immune cell-specific geneset analysis (https://www.itcc-consortium.org/). We used a neuroblastoma dataset as a validation set for immune checkpoint profiling that was available on CAVATICA pounder accession number phs001436.c1 [26]. This dataset was not included in the pedNST cohort, as tumour types (primary, relapsed, etc.) were not known.

Transcriptome data processing

RNA-seq reads from TCGA and TARGET were aligned to human genome 38 (hg38) as described on the GDC website (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline). We used RSEM output for the TCGA and TARGET from Toil [27]. The CBTN raw RNA-seq reads were aligned to hg38 using the STAR v2.5.2b [28] and quantified using RSEM v1.2.28 [29] (detailed workflow can be accessed under https://github.com/kids-first/kf-RNA-seq-workflow). The ICGC data were aligned and quantified using STAR [28] v2.7.6a and RSEM [29] v1.2.28, respectively. With the exception of the ICGC, all data processing was performed on the CAVATICA data analysis platform [22]. For immune inference, we used transcripts per million (TPM), as per instructions [30]. For gene set enrichment scores and consensus clustering, we used log₂ transformed TPM values corrected for batch effects (data source) using ComBat function from the sva [31] R package. To adjust for cancer type differences, we normalized log₂ transformed batch-corrected TPM values using median expression in each cancer type.

Variant and fusion calling

For the TCGA and TARGET, variant calls were obtained from whole exome sequencing (WES) while the ICGC and CBTN variant calls were from whole genome sequencing (WGS). For the TCGA and TARGET, we restricted our analyses to those with MuTect2 [32] calls available on the GDC (details available at https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline). For the ICGC, we used SNV and Indels calls as previously described [21]. CBTN WGS data were processed using MuTect2 [32] and Strelka2 [33] (detailed workflows at https://github.com/kids-first/kf-somatic-workflow.) and only overlapping calls were used for downstream analysis. For all datasets, we excluded somatic variants with less than 3% variant allele frequency. The Arriba [34] workflow was used for gene fusion calling on CBTN and ICGC datasets (https://github.com/suhrig/arriba). Detailed workflows are available at https://github.com/DKFZ-ODCF/RNAseqWorkflow [35] and https://github.com/kids-first/kf-rnaseq-workflow [36].

TMB and oncogenic pathways

To ensure all variants in our study were covered across WES and WGS datasets, we generated a common region list by intersecting.bed files used in exome experiments with 50 base pairs padding. The exome kits included Agilent Custom V2 Exome Bait (TCGA and TARGET), Agilent SureSelect All Exon 38 Mb V2 (TCGA), Agilent SureSelect All Exon 50 Mb (TCGA), SeqCap EZ Exome Probes v3.0 (TCGA), SeqCap EZ Human Exome Library v2.0 (TARGET) and SeqCap EZ HGSC VCRome 2.1 (TCGA). The final common.bed file consisted of 30,028,393 base pairs. We then used this region list to subset the WGS variant calls. We included non-synonymous coding SNV and Indels in the TMB calculations.

For molecular pathways and their associations with immune clusters, we identified samples with at least one alteration in genes involved in ten TCGA oncogenic pathways [37]. The Receptor Tyrosine Kinase (RTK) pathway contained 71 affected genes with a total of 419 alterations in 274 samples that encompassed 71% of pedLGG. Twenty six genes in the PI3-kinase pathway were altered in 74 pedNST samples, 32% of which were pedHGG. The Wnt pathway was disrupted in 70 samples, of which 36 and 20% were medulloblastoma and craniopharyngioma. We found 155 alterations in 51 genes involved in the Notch signalling pathway across 66 samples, primarily medulloblastoma (45%) and pedHGG (23%). Cancer entities with alterations in 32 genes involved in the HIPPO pathway included neuroblastoma (25%), pedHGG (22%) and medulloblastoma (20%). Five core genes involved in the TP53 pathway accounted for 56 mutations across 51 pedNST samples. At least one of 10 cell cycle genes was mutated in 16 pedNST samples, primarily in pedHGG (56%). Acknowledging that MYCN amplifications were not included in this analysis, we found 26 alterations in 10 genes related to MYC signalling in 15 samples. Nine samples showed alterations in the TGF-β pathway that consisted of 7 altered genes. NRF2 pathway contained three genes (KEAP1, CUL3 and NFE2L2) with 7 mutations across 6 samples. We studied the effects of somatic alterations in the Mismatch Repair (MMR) pathway on immune microenvironment by comparing samples with at least one somatic alteration in MLH1, MSH2, MSH6, PMS2, POLE or POLD1 to those without any alterations.

Immune infiltration and in silico simulations

We used ESTIMATE with default parameters to measure overall immune infiltration (https://bioinformatics.mdanderson.org/public-software/estimate/) [38]. We performed in silico simulations using bulk RNA-seq data from PDX models derived from medulloblastoma and RNA-seq data of four immune cell types (CD8 + T-, CD4 + T-, B-cell and NK cells) downloaded from ENCODE portal under identifiers ENCSR861QKF [39] (CD8 + T-cells), ENCSR463JBR [40] (CD4 + T-cells), ENCSR449GLL [41] (B-cells), ENCSR357XTU [42] (NK cells)). For each simulation, an equal percentage of reads from four immune cell types was randomly sampled using SAMtools [43] v1.9, converted to fastq files with BEDTools [44] v2.27.1 and then concatenated together with randomly sampled reads from PDX RNA-seq data to a total of 100 million reads. The generated pseudo-samples were then processed with Kallisto and Sleuth [45, 46], and TPM values were used as input for the ESTIMATE.

Immune microenvironment analysis

We used TIMER2 [30] web interface (http://timer.cistrome.org/) for comprehensive analysis of immune cell composition using six computational tools, CIBERSORT [47], EPIC [48], QUANTISEQ [49], MCPCOUNTER [50], TIMER [51] and XCELL [52]. We used the CRI-iAtlas Shiny app (https://isb-cgc.shinyapps.io/iatlas/) to cluster the pedNST samples with the immune subtype classifier previously published for adult cancers [14]. To estimate immune cell composition using ICGC methylation array data, we used EpiDISH (Epigenetic Dissection of Intra-Sample Heterogeneity) R package that uses cell-type-specific DNAse Hypersensitive Site (DHS) data and robust partial correlations to infer immune cell compositions [53]. We used a reference signature matrix previously constructed and used for immune deconvolution analysis of paediatric CNS tumours [54].

Development of immune cell-type gene sets and consensus clustering

To identify genes specific to the immune system, we leveraged four data sources: (1) We obtained median gene expression for 18 purified cell types, as previously reported [13]. We found 9897 genes were expressed > 75th percentile of normalized expression (median centred and Median Absolute Deviation (MAD) scaled) in at least one cell type within this dataset. We selected 1958 genes with ≥ 2 MAD difference between immune and non-immune cell populations. (2) We downloaded the human protein atlas data v20.1 (https://www.proteinatlas.org/about/download). The Human Protein Atlas consisted of single-cell data from 51 cell populations [55]. 15,302 genes were expressed > 75th percentile in at least one cell population. (3) In the Human Protein Atlas and across 37 tissues, 11,069 genes were annotated as specific to at least one tissue, determined as normalized expression ≥ 1 (median normalized expression 33.9). We found 7206 genes specific to only one tissue, 1257 of which were specific to blood, lymphoid tissue or bone marrow. (4) We included the ESTIMATE immune signature consisting of 141 genes [38]. In total, we found 3041 unique genes across four data sources with evidence of specificity to immune cell populations.

Next, we sought to identify genes that may be expressed in paediatric cancer and other non-immune cells. We compiled a paediatric cancer geneset using three data sources. (1) We used protein-coding gene expression data from 79 PDX models. We found 13,788 genes expressed > 75th percentile in at least one PDX. (2) We used data from 22 cell lines collected by CBTN. These consisted of one ependymoma, 18 pedHGG and 3 medulloblastoma. Across these cell lines, 11,195 protein-coding genes were expressed > 75th percentile in at least one cell line. (3) We leveraged published single-cell datasets from paediatric cancers, neuroblastoma [16] and ependymoma [18, 19] and used original cell annotations as reported by authors. As single-cell datasets are generally sparse, we focused our analysis to genes with maximum expression > 75th percentile in at least one cell population in each dataset. We then determined non-immune genes as those with ≥ 2 MAD difference between non-immune and immune cell populations in each dataset. This analysis resulted in 3036, 12,109 and 5040 non-immune genes from Jansky et al. [16], Gojo et al. [18] and Gillen et al. [19] datasets, respectively. Intersecting the immune and non-immune genes resulted in 791 immune-specific genes.

To accurately assign the identified genes to immune cell types, we compiled genesets from seven sources and conducted a consensus approach to assign genes to specific immune cell types. For this analysis, we aggregated immune cell subtypes (e.g. Tregs) into the following major immune cell types: T, B, NK cells, dendritic cells, monocytes, macrophages, granulocytes and myeloid cells. These sources included gene annotations from five immune deconvolution tools, MCPCOUNTER [50], QUANTISEQ [49], CIBERSORT [47], EPIC [48] and TIMER [51]. We used all 264 genesets that were used to develop XCELL [52]. We used two additional sources: (1) scaled average counts of 7172 genes specific to 28 purified immune cell types derived from 416 healthy donors and patients with immune diseases (Immunexut) [56]. We chose 3896 genes with scaled count > 75th percentile in at least one cell type, of which 2555 genes had ≥ 2 MAD difference between one cell type compared to all others. (2) The human protein atlas annotated 5934 genes as blood cell specific [57], determined as genes with normalized expression ≥ 1 in at least one cell type. We chose 4902 genes that were specific to one immune cell type. Our consensus analysis revealed 216 genes that were assigned to one specific immune cell type in at least two independent sources. We excluded myeloid cell geneset due to the low number of genes (n = 6) that may lead to inaccurate enrichment scores. To calculate enrichment scores, we applied single-sample geneset enrichment analysis (ssGSEA) from GSVA R package [58] to batch-corrected log₂-transformed TPM values across pedNST. We then applied consensus clustering to the scaled (median centred and MAD scaled) immune cell-type enrichment scores using the ConsensusClusterPlus R package [59]. We performed k-means clustering based on Euclidean distance with 200 subsamples using 80% of samples and 100% of features for 2 to 8 clusters. This analysis revealed four major immune clusters across pedNST.

For further delineate types of infiltrating T-cells, we derived enrichment scores of 40 T-cell subtypes using top 50 signature genes derived from recent pan-cancer analysis of tumour-infiltrating lymphocytes [60]. We focused this analysis on Paediatric Inflamed to avoid overestimation as a result of low overall immune infiltration in other clusters. We scaled enrichment scores within each T-cell subtype across Paediatric Inflamed and performed consensus clustering with parameters as previously described. For myeloid cell subtypes, we used 13 myeloid signatures derived from single-cell analysis of adult cancer patients [61]. We chose 9 gene signatures consisting of genes with ≥ 1.5 fold change in expression between any given myeloid cell-type and other clusters. We scaled enrichment scores from myeloid gene signatures within each myeloid cell subtype across Myeloid Predominant and performed consensus clustering with parameters as previously described. For microglia, we used signatures from two previous reports and derived enrichment scores in the pedCNS subset of Myeloid Predominant using ssGSEA [62, 63].

Immunohistochemistry

Formalin-fixed, paraffin-embedded TMA sections were analysed for CD4, CD8 and CD19 expression. The IHC staining was performed using the Ventana Discovery platform. IHC was optimized and performed with CD4 (Abcam Ab183685), CD8 (Leica NCL-L-CD8-4B11) and CD19 (e-Bioscience 14–0194) with dilutions of 1:500, 1:100 and 1:500 respectively. In brief, tissue sections were incubated in Tris EDTA buffer (cell conditioning 1; CC1 standard) at 95 °C for 1 h to retrieve antigenicity, followed by incubation with the respective primary antibody for 1 h. Bound primary antibodies were incubated with the respective secondary antibodies (Jackson Laboratories) with 1:500 dilution, followed by Ultramap HRP and Chromomap DAB detection. For staining optimization and to control for staining specificity, normal tonsil was used as control. Intensity scoring was done on a common 4-point scale. Descriptively, 0 represents no staining, 1 represents low but detectable degree of staining, 2 represents clearly positive staining and 3 represents strong expression. Expression was quantified as H-score, the product of staining intensity and percentage of stained cells. The TMAs used in the study were obtained from the Children’s Oncology Group (COG) and contains 9 ATRT, 13 pedHGG, 20 ependymoma, 64 medulloblastoma and 33 neuroblastoma cases, all of which are represented in duplicate cores.

To study protein levels of TIM3 and LAG3 in neuroblastoma, we purchased two serial TMA slides consisting of 26 neuroblastoma cases all in duplicate cores with 10 cores of normal peripheral nerve tissue (Biomax, NB642a). Briefly, TMA sections were dewaxed and rehydrated through an ethanol series to water and endogenous peroxidases were blocked using 3% H₂O₂ in PBS for 15 min at room temperature. Antigen retrieval was performed using 10 mM sodium citrate pH 6.0. Primary antibodies were incubated for 45 min at room temperature (TIM3 CST45208 1:150, LAG3 ab40466 1:200) followed by several washes with PBS. Secondary antibodies were incubated for 30 min at room temperature (BA-1000 and BA-9200 at 1:500 dilutions Vectors Labs). ABC kit (PK-6100, Vectors Labs) was applied for 25 min with DAB for 4 min (SK-4100) followed by PBS washes. Stained slides were scanned using a Nanozoomer 2.0HT (Hamamatsu Photonics) at 20 × or 40 × .

Computational tumour-infiltrating lymphocyte (TIL) analysis was performed on digitally scanned haematoxylin and eosin (H&E)-stained whole-slide images (WSI) from the Children’s Brain Tumour Network (CBTN) as previously described [64]. Briefly, a pre-trained Inception-V4 deep learning model was applied to each WSI, returning a model probability of TIL presence within 50 × 50 micron patches tiling across the entire WSI. Image patches that did not contain tissue were filtered out using a color standard deviation threshold of 18 (computed on an 8-bit [0,255] color scale), and a mean TIL score was computed as the average TIL probability on the remaining patches for each WSI. Upon manual review, WSI with poor scanning quality and/or with a high false positive rate of TIL detection were removed from the analysis, including medulloblastoma and other tumours arising in the cerebellum/posterior fossa. One outlier sample was removed (Grubbs test).

HLA-typing and neoantigen prediction

We used the Optitype [65] tool available on the CAVATICA platform to determine HLA-A, B and C types from the TARGET and CBTN RNA-seq datasets. To predict neoantigens from gene mutations, we used Mutant Peptide eXtractor and Informer (MuPeXI) tool using MuTect2 variant caller as input [66]. We predicted 8-11mer peptides for all HLA class I types. We determined strong or weak binding peptide-MHC complexes as previously described (percentile rank ≤ 0.5 for strong binding, > 0.5 and ≤ 2 for weak binding peptides) [67]. We included peptides predicted from mutations that overlapped between MuTect2 and Strelka2.

T- and B-cell repertoire analysis

To recover T- and B-cell clonotypes, we used MiXCR v2.1.12 [68] with default parameters for RNA-seq data processing. We applied the framework iNterpolation/EXTrapolation (iNEXT) to study immune repertoire diversity [69]. This method was primarily developed for diversity estimates in ecology and aimed to bring together asymptotic estimates and rarefaction/extrapolation methods of estimation for samples with different sizes. Three most popular diversity indices have already been translated to immune repertoire studies; richness, the total number of species, Shannon index (entropy), which puts moderate weight on abundance and Simpson index, which measures diversity of abundant species. The robustness of these indices is dependent on sample size and could be biased in under-sampled experiments, which is a known problem in biodiversity studies and attempts have been reported to reduce such empirical biases [70]. Chao et al. provided a more accurate estimation of diversity by integrating slopes obtained from accumulation curves into entropy formula followed by bootstrapping method for variance [70]. We applied similar principles on RNA-seq datasets to infer immune diversity. We found that the estimates were severely affected in samples with less than three clonotypes recovered from bulk RNA-seq or in samples with multiple clonotypes of the same clonal fractions and excluded them from our analysis. Note that these samples may be (A) samples with low number of immune cells and therefore low number of clonotypes that could not be reflected on shallow RNA-seq data, or (B) samples with sufficient immune infiltration, but highly uneven clonal distribution, therefore only highly abundant clonotypes were seen at shallow coverage, yet the true diversity could not be robustly estimated.

For B-cell repertoire, we used constant regions of immunoglobulin heavy chain (IGH C segments) to study the distribution of immunoglobulin isotypes. Due to somatic hypermutation in the B-cell repertoire, we removed sequences with ≤ 2 reads. As B-cell clones originating from a naïve B-cell may have mismatches in CDR3 sequences via somatic hypermutation events, we clustered IGH CDR3 octamers allowing for one mismatch, as described previously [71]. We considered sequences shorter than 8 amino acids as individual CDR3 sequences. Only samples with > 3 CDR3 sequences were included from this analysis and we used the gini index of inequality as a measure for uneven distribution of B-cell clusters.

TCR sequencing

To validate the T-cell diversity estimates from RNA-seq, we applied CapTCR-seq [72] hybrid capture protocol to RNA-seq libraries of adult nasopharyngeal carcinoma (n = 33), paediatric samples from the PRecision Oncology For Young peopLE (PROFYLE) programme (www.profyle.ca) (n = 10) and ten paediatric GBM samples (ICGC). We compared estimated values of TCRβ richness, Shannon and Simpson diversities from bulk RNA-seq data with observed richness, Shannon and Simpson diversities from CapTCR-seq experiment.

Statistics and visualization

Statistical tests, p-values and other details are noted in text and figure legends. For tumour subtypes, we used one-sided Fisher’s exact test (alternative = ‘greater’ in R) comparing each subtype to all others within each immune cluster. We used analysis of covariance (ANCOVA) to compare continuous variables such as TMB or age at diagnosis across immune clusters adjusting for cancer entities. For survival analysis, we defined overall survival (OS) as follows. For NCI TARGET, we used ‘Overall Survival Time in Days’ and ‘vital status’ reported in the metadata. For CBTN, we used ‘Last Known Status’ with ‘Deceased-due to disease’ recoded as ‘Event’. We calculated OS time using ‘Age at Collection’ of the first collected specimen (‘Initial CNS Tumor’) and ‘Age At Last Known Status’ of the last collected specimen, as reported in the metadata. For ICGC, we defined OS as the time between the date of first diagnosis and date of death or last follow-up. We defined progression-free survival (PFS) as follows. For NCI TARGET, we used ‘Event Free Survival Time in Days’ and ‘First Event’ as reported in the metadata with ‘Event’, ‘Death’, ‘Progression’, ‘Relapse’ and ‘Second Malignant Neoplasm’ recoded as ‘Event’. For CBTN cases with one surgical resection, we used ‘Last Known Status’ with ‘Deceased-due to disease’ recoded as ‘Event’. For these cases, we calculated PFS time using ‘Age at Collection’ and ‘Age At Last Known Status’ of the collected specimen, as reported in the metadata. For cases with multiple surgical resections, we calculated PFS time using ‘Age at Collection’ of the first collected specimen (‘Initial CNS Tumor’) and ‘Age at Collection’ of the second collected specimen (‘Progressive’ or ‘Recurrence’), as reported in the metadata. For ICGC, we calculated PFS as the time between date of first diagnosis and date of first relapse/progression or death (if death was the first event), or date of last follow-up (if no event).

We used log-rank tests for Kaplan–Meier analyses. We used Cox proportional hazards models to adjust for cancer entities and gender in multivariable analyses. We used Student’s t test to compare scaled values across groups and rank sum test to compare values with non-normal distribution. We applied the Cochran-Mantel–Haenszel (CMH) test to compare samples with mutations in oncogenic pathways across immune clusters controlling for the effects of cancer entities. In all boxplots, boxes show median and IQR and whiskers represent 1.5 times IQR. Boxplots are shown for groups with more than three datapoints. We annotated statistical significance levels as follows: *p < 0.05, **p < 0.01 and ***p < 0.001. We performed all analyses and visualizations in R v4.0. We used Adobe Illustrator v24.0.1 for aesthetic edits and figure alignments.

Results

Establishment of paediatric nervous system tumour cohort for immunogenomic analysis

To study immune attributes across paediatric nervous system tumours (pedNST), we compiled a cohort of bulk RNA-seq data from three consortia: the Children’s Brain Tumour Network (CBTN), the National Cancer Institute Therapeutically Applicable Research To Generate Effective Treatments initiative (NCI TARGET) [23] and the International Cancer Genome Consortium (ICGC) [24]. To better represent primary pedNST and to enable comparison at the transcriptional level, we performed unsupervised clustering on the CBTN RNA-seq dataset (n = 581, Additional file 1: Fig. S1A-C). We further annotated CBTN medulloblastoma and ATRT tumour subtypes using unsupervised clustering for downstream analysis (Additional file 1: Fig. S1D-E, ‘Methods’). We included 195 primary samples from the ICGC [21] with matched RNA-seq, whole genome sequencing (WGS) and DNA methylation data and 149 primary neuroblastomas from the NCI TARGET [23] with matched RNA-seq and whole exome sequencing (WES) data [23].

The final, aggregated non-overlapping paediatric dataset for immunogenomic analysis consisted of 925 samples with primary locations in the central nervous system (CNS) or peripheral nervous system: embryonal tumours with multilayered rosettes (ETMR, n = 9), neurofibroma (NFB, n = 11), choroid plexus tumours (CP, n = 16), meningioma (MNG, n = 13), schwannoma (SCHW, n = 14), craniopharyngioma (CPH, n = 27), atypical teratoid/rhabdoid tumours (ATRT, n = 30), ependymoma (EPN, n = 65), paediatric high-grade glioma (pedHGG, n = 83), neuroblastoma (NBL, n = 151), medulloblastoma (MB, n = 208), and paediatric low-grade glioma (pedLGG, n = 298) (Fig. 1A) (Additional file 2: Table S1). We included data from 79 patient-derived xenograft (PDX) models as part of the Innovative Therapies for Children with Cancer Pediatric Preclinical Proof-of-Concept Platform (ITCC-P4) project as negative controls lacking immune or stromal infiltration [73, 74] (Fig. 1B).

Fig. 1

Transcriptional analysis of 925 paediatric nervous system tumours (pedNST) reveals four distinct immune clusters. A Overview of cohorts and sample size in the present study (PDX: patient-derived xenografts, ETMR: embryonal tumour with multilayered rosettes, NFB: neurofibroma, MNG: meningioma, SCHW: schwannoma, CP: choroid plexus tumours, CPH: craniopharyngioma, ATRT: atypical teratoid/rhabdoid tumour, EPN: ependymoma, pedHGG: paediatric high-grade glioma, NBL: neuroblastoma, MB: medulloblastoma, pedLGG: paediatric low-grade glioma). B Distribution of percentage immune reads based on ESTIMATE immune score across paediatric and adult cancers. PDX samples serve as a negative control for immune infiltration. PRAD: prostate adenocarcinoma, LGG: low-grade glioma, OV: ovarian serous cystadenocarcinoma, SKCM: skin cutaneous melanoma, COAD: colorectal adenocarcinoma, GBM: glioblastoma multiforme, LUAD: lung adenocarcinoma. C Heatmap representing consensus clustering using enrichment scores derived from six immune-cell-specific genesets across the pedNST cohort. Cell-type scores correspond to normalized gene set enrichment scores. D Heatmap illustrating proportion of samples in each immune cluster across cancer entities. Barplots show the total number of samples for each row (cancer type) and column (immune cluster). E Heatmap showing fraction of samples in each CRI-iAtlas cluster across immune clusters in pedNST. Barplot shows the total number of samples in each CRI-iAtlas cluster. F Boxplot showing average tumour-infiltrating lymphocyte (TIL) scores as determined by segmentation analysis of pathological images across the paediatric CNS tumour samples (CBTN). Three sample images are shown representing 1% (bottom), 10% (middle) and 15% (top) TIL scores corresponding to lower quartile, mean and higher quartile, respectively. Boxes show median and interquartile range (IQR) and whiskers represent 1.5 times IQR. Two-sided rank sum test, *p < 0.05. G Barplots showing frequency (top barplot) and fraction (stacked barplot) of tumour subtypes across immune clusters (SHH: Sonic Hedgehog, TYR: Tyrosine, MYCN-NA: MYCN non-amplified, MYCN-A: MYCN amplified, WT: wildtype)

Immune infiltration analysis reveals high variability across and within cancer types

To determine the overall levels of immune infiltration across pedNST, we first assessed the performance of the deconvolution tool ESTIMATE [38] using in silico simulations and converted immune scores to ‘immune read percentage’ using a non-linear regression model (Additional file 1: Fig. S2B). To enable comparison with adult brain tumours and cancers in which ICI agents have been clinically studied [75,76,77,78,79], we performed a parallel analysis of 2452 primary adult tumour samples from The Cancer Genome Atlas (TCGA) [25, 80]: glioblastoma multiforme (GBM, n = 153), low-grade glioma (LGG, n = 507), skin cutaneous melanoma (SKCM, n = 102), colorectal adenocarcinoma (COAD, n = 298), ovarian serous adenocarcinoma (OV, n = 373), prostate adenocarcinoma (PRAD, n = 497), lung adenocarcinoma (LUAD, n = 522).

Immune read percentages varied considerably within and across pedNST (median 7.2%, range 1.9–47.9%) and in similar range compared to adult CNS tumours (median 7.8%, range 2.9–41.8%). Common non-CNS adult cancers showed higher overall levels of immune infiltration (median 10.3%, range 2.5–67.8%). In pedNST, the highest values were in extracranial entities including neurofibroma, schwannoma and meningioma (medians 22.6, 17.8 and 15.3%, respectively). Paediatric brain tumour entities ETMR and medulloblastoma had the lowest median immune infiltration of all cancers analysed (medians 2.7 and 2.8%) [81,82,83] (Fig. 1B). Choroid plexus tumours had consistently low variation of immune read percentage (range 3.1–13%, median 7.3%). In contrast, ATRT exhibited the widest distribution of immune infiltration (range 2.3–47.9%, median 5.5%), consistent with known associations between infiltration and ATRT subgroups [84, 85]. Similarly, a wide range of immune read percentages were observed in craniopharyngioma and neuroblastoma (range 3.9–40.9%, median 12.6% and 2.6–30.7%, median 8.8%). These data revealed notable examples of highly infiltrated and immune-excluded samples in paediatric nervous system tumours and identified general trends of immune infiltration.

To investigate whether the wide variability in immune infiltration could be recapitulated at the protein level, we applied immunohistochemistry (IHC) to a tissue microarray (TMA) from an independent cohort of 139 paediatric cancers, obtained from Children’s Oncology Group (COG). The median CD8 staining was greatest in ATRT followed by neuroblastoma and pedHGG, while medulloblastoma and ependymoma had median H-score (defined as product of staining intensity and percentage of stained cells) of zero (Additional file 1: Fig. S3A-B). We found CD4 staining in eight samples (H-score > 0, 3 ATRTs, 3 ependymoma, 2 medulloblastoma), while other samples showed no CD4 staining. CD19 staining was variable among cancer types ranging from 3% (medulloblastoma) to 51% (neuroblastoma) (Additional file 1: Fig. S3A-B). Summing H-scores across all three markers, we found highest staining in ATRT followed by neuroblastoma (medians 17 and 15) and medulloblastoma showing the lowest staining. These results confirm the high variability of immune cell-type infiltration within and across pedNST, as well as generalized trends inferred from gene expression.

Consensus gene set clustering reveals four distinct immune clusters in pedNST

In light of the variability of immune infiltrates within each tumour type, we sought to categorize immune microenvironments across pedNST. We detected unexpected immune signals inferred by various immune deconvolution algorithms when applied to PDX RNA-seq data (Additional file 1: Fig. S4A-G, Additional file 2: Table S2), which should not yield a (human) immune cell signature. This indicated wide discordance and non-specific signals across existing immune deconvolution tools. To address the non-specific immune signal that may originate from tumour cells, we sought to identify immune cell-specific genes that lack expression in paediatric nervous system cancer cells. We compiled 3041 immune-related genes by incorporating data from four sources: ESTIMATE signature [38], an immune cell compendium [13], the Human Protein Atlas (HPA) single cell [55] and tissue specificity datasets [86] (Methods). We excluded genes with evidence of expression in paediatric cancer cells using data from the PDX models (ITCC-P4), single-cell RNA-seq [16, 18, 19] and established cell lines (CBTN). This analysis identified that 791 of the 3041 immune genes are not expressed in paediatric nervous system tumour cells. We performed consensus analysis using gene sets from aforementioned immune deconvolution tools along with gene expression profiling of 28 immune cell types [56], and the HPA blood-cell-specific dataset [57]. Using this approach, 216 genes were assigned to a single cell-type in at least two data sources. These included genes specific for T-cells (n = 41), B-cells (n = 32), NK cells (n = 20), monocytes (n = 18), dendritic cells (n = 28) and granulocytes (n = 71) (Additional file 2: Table S3). Consensus clustering of normalized enrichment scores for these cell-specific genesets identified four distinct immune clusters (Fig. 1C,D) that we characterized further and designated Paediatric Inflamed, Myeloid Predominant, Immune Neutral, and Immune Desert based on differential expression of aforementioned genesets.

To place these clusters in context of prior immunogenomics work in adult tumours, we applied the CRI-iAtlas, a comprehensive classification method performed on TCGA dataset consisting of over 10,000 tumour samples spanning 33 adult tumour types [14], to pedNST. We found that a lower proportion of paediatric brain tumours were grouped in ‘Immunologically quiet’ or ‘Lymphocyte depleted’ clusters compared to adult counterparts (85%, 635/749 vs 98%, 654/668). In contrast, 16.5% of paediatric extracranial tumours (neuroblastoma, neurofibroma and schwannoma, 29/176) harboured cold immune microenvironment, compared to 10.5% of adult extracranial tumours (888/8458). Across pedNST, we found 19% of pedNST (n = 174) belonged to the ‘Inflammatory’ or ‘IFN-γ dominant’ clusters, while 7% (n = 64) were grouped in ‘Wound healing’. Twenty percent of pedLGG (60/298) were grouped in the ‘Inflammatory’ cluster, in contrast to adult LGG, 98% of which showed a cold immune microenvironment. Almost all pedHGG samples (82/83) were grouped in the ‘Immunologically quiet’ or ‘Lymphocyte depleted’ clusters, similar to adult GBM (151/154). Overall, we found 72% of pedNST (n = 664) were grouped in the ‘Immunologically quiet’ or ‘Lymphocyte depleted’ clusters indicating a generally cold immune microenvironment in pedNST (Fig. 1E).

The Paediatric Inflamed cluster (n = 90, 9.7%) had the highest immune read percentage across pedNST. 57.1% of schwannoma (n = 8) and 54.5% of neurofibroma (n = 6) were clustered in this group followed by 40.7% of craniopharyngioma (n = 11) and 30.8% of meningioma (n = 4) (Fig. 1D). Five of six cases grouped in the ‘TGF-β dominant’ cluster belonged in this cluster, followed by 42.3% of ‘IFN-γ dominant’ samples (n = 11) (Fig. 1E). This cluster was devoid of samples in the ‘Immunologically quiet’ or ‘Lymphocyte depleted’ cluster indicating that ~ 10% of pedNST samples harbour T-, B- and NK cells, monocytes, granulocytes and dendritic cells.

The Myeloid Predominant cluster (n = 279, 30.1% of pedNST) scored highly for monocyte, dendritic cell and granulocyte gene sets while harbouring lower levels of lymphoid cell types. 53.8% of meningioma (n = 7) and 53% of pedLGG (n = 158) cases clustered in Myeloid Predominant followed by 44.5% of craniopharyngioma (n = 12) (Fig. 1D). 50% of the samples within the CRI-iAtlas ‘IFN-γ dominant’ (n = 13) and 40% of the ‘Inflammatory’ samples (n = 59) were clustered in this group (Fig. 1E), suggesting the inflammatory component of this cluster is driven primarily by myeloid cells.

The Immune Neutral cluster (n = 393, 42.5%) had myeloid and lymphoid cell infiltration scores near the median of the entire pedNST cohort. This cluster contained a large fraction of the cancer types with low immune read percentage including ETMR (n = 5, 55.6%), medulloblastoma (n = 136, 65.4%) and ependymoma (n = 47, 72.3%) (Fig. 1D). This cluster included 56.4% of the ‘Immunologically quiet’ (n = 62) and 48.7% of the ‘Lymphocyte depleted’ (n = 270) tumours (Fig. 1E), indicating that intermediate immune infiltration across pedNST corresponds to cold immune microenvironment compared to adult cancers, which may contribute to intrinsic resistance to immune checkpoint inhibitors.

The Immune Desert cluster (n = 163, 17.6%) received low immune inference scores for all immune cell types (Fig. 1C). This cluster included 32.5% of neuroblastoma (n = 49) and 29.8% of medulloblastoma (n = 62) and none of the immune infiltrated cancer entities (Fig. 1D). Purity scores inferred from copy number alterations [87] in 156 ICGC samples revealed that samples in this cluster were of highest cancer cell content compared to Immune Neutral or Myeloid Predominant (Additional file 1: Fig. S5A-B).

To validate lymphocyte infiltration levels across the immune clusters, we leveraged a deep learning method of tissue image analysis for 355 pedCNS samples from CBTN with haematoxylin and eosin (H&E)-stained tissue slides and matched RNA-seq [64], as well as 195 ICGC samples with matched DNA methylation arrays and RNA-seq. Consistent with the RNA-seq immune inference analysis, choroid plexus tumours received the lowest tumour-infiltrating lymphocyte (TIL) scores determined by H&E tissue analysis (median 0.08), while meningioma and neurofibroma scored highest across cancer entities (medians 0.14 and 0.11, Additional file 1: Fig. S6A). TIL scores were variable across CNS tumour sites ranging from a median of 0.07 in cerebral hemispheres to a median of 0.18 in case of meninges (Additional file 1: Fig. S6B). Paediatric Inflamed showed significantly higher TIL score compared to Myeloid Predominant and Immune Desert (medians, 0.1, 0.09 and 0.08, respectively, two-sided rank sum test, p = 0.02 and 0.03, Fig. 1F, Additional file 2: Table S4). This difference was more pronounced in samples with immune read percentage > 5% (medians 0.1 and 0.03 in Paediatric Inflamed and Immune Desert) or in those collected from cerebral hemispheres (medians 0.09 and 0.05), suggesting that high immune infiltration and composition of tumour microenvironment may influence the H&E-based inference of TILs. These findings were recapitulated by immune inference analysis utilizing the ICGC DNA methylation array data [53, 54]. Paediatric Inflamed tumours displayed higher levels of T-cell infiltration (pairwise two-sided Student’s t test with Bonferroni correction, p < 0.01, Additional file 1: Fig. S6D). Although B-cell estimates were not significantly different between Paediatric Inflamed and Myeloid Predominant, lower levels of B-cells were estimated in Immune Neutral and Immune Desert compared to Myeloid Predominant (p < 0.001, Additional file 1: Fig. S6D). Similarly, Immune Neutral and Immune Desert showed depletion of DNA methylation-based estimates for NK cells and monocytes compared to Myeloid Predominant, in agreement with our RNA-seq clusters (p < 0.001, Additional file 1: Fig. S6D). We compared levels of CD8 + T-, CD4 + T- effector and regulatory T-cell (Treg) subtypes estimated from DNA methylation deconvolution across immune clusters. Paediatric Inflamed showed elevated levels of CD8 + T-cell and Treg estimates, while depleted in CD4 + T effector estimates. Conversely, Myeloid Predominant showed high levels of Treg and CD4 + T effector estimates with minimal CD8 + T-cells (Additional file 1: Fig. S6E). These results suggest presence of a lymphocyte population in Myeloid Predominant that is undetectable in bulk RNA-seq.

As pedNSTs consist of distinct molecular entities and subgroups [88, 89], we next investigated associations between cancer subtypes and immune clusters (Fig. 1G). MYC-like ATRTs appeared to be enriched in Paediatric Inflamed or Myeloid Predominant compared to non-MYC-like ATRT, and 8 out of 12 SHH-like ATRT tumours clustered into Immune Neutral, although it did not reach statistical significance (one-sided Fisher’s exact test with Bonferroni correction, p > 0.05) [84]. In neuroblastoma (n = 148), MYCN non-amplified neuroblastoma cases were enriched in Myeloid Predominant (n = 41, p = 0.007), while neuroblastoma MYCN amplified showed a trend toward Immune Desert (n = 16, p = 0.1) [90, 91]. SHH medulloblastoma showed a trend toward Myeloid Predominant and Immune Neutral (n = 4 and 28, p > 0.05), while 95% of WNT medulloblastoma were grouped into Immune Neutral and Immune Desert (19 out of 20). In pedLGG (n = 298), BRAF wildtype samples were enriched in Immune Desert (n = 24, p = 5 × 10⁻⁶), BRAF-KIAA1549 fusion-positive tumours were grouped into Myeloid Predominant (n = 82, p = 0.005) and BRAF p.V600E samples showed a trend toward Paediatric Inflamed (n = 7, p = 0.07). In summary, our findings reveal associations between tumour-intrinsic characteristics and the immune microenvironment, although it remains to be established whether the microenvironment is sculpted by the cancer cells or whether specific cancer subtypes thrive in specific microenvironments.

We next investigated whether there were associations between patient characteristics or outcome with the immune clusters. We found an association between male sex and pedNST, more prominently seen in medulloblastoma and ependymoma (Additional file 1: Fig. S7A). A logistic regression model adjusting for cancer type indicated a trend toward an association between male genetic sex and Immune Desert (Wald test, p = 0.05) (Fig. 2A). There were no significant associations between race or age and immune clusters independent of cancer type (Fig. 2B, C, Additional file 1: Fig. S7B-C). Kaplan–Meier analysis indicated significant differences in overall survival (OS) and progression-free survival (PFS) among the immune clusters (log-rank test, p = 0.001 and 0.02) (Fig. 2D, E). We observed poorer OS for Immune Neutral and Immune Desert compared to Paediatric Inflamed in a Cox proportional hazards model adjusting for cancer type and sex (hazard ratio (HR) = 1.56 and 1.97, confidence interval (CI) 1.13–2.14 and 1.39–2.79, p = 0.006 and 0.0001) (Additional file 2: Table S5). Myeloid Predominant cases showed a trend toward improved PFS compared to Paediatric Inflamed (HR = 0.68, CI 0.44–1.06, p = 0.08) (Additional file 2: Table S6). This trend may be a result of low but differential presence of T-cells in Myeloid Predominant compared to Paediatric Inflamed (Additional file 1: Fig. S6E). Our results suggest that independent of cancer type, infiltration of specific immune cell populations (Paediatric Inflamed and Myeloid Predominant) provides a survival advantage over low infiltration (Immune Neutral and Desert) that may indicate poor immune recognition.

Fig. 2

Associations of immune clusters with patients’ clinical parameters and survival. A, B barplots showing distributions of gender (A) and race (B) across immune clusters in pedNST. ANCOVA, not significant. C boxplot depicting age distribution across immune clusters. ANCOVA, not significant. Boxes show median and interquartile range (IQR) and whiskers represent 1.5 times IQR. D, E Kaplan–Meier curves and risk tables for overall survival (D) and progression-free survival (E) among immune clusters. Log-rank test p-values are denoted

Gene expression analysis reveals differential molecular pathways and immunoregulatory genes in immune clusters

To understand the cellular pathways underlying each immune cluster, we conducted quantitative set analysis for gene expression (QuSAGE) [92] comparing samples in one cluster against all others while adjusting for data source and cohort. Leveraging the Molecular Signatures Database (MSigDB) Hallmark genesets [93], we found that immune-related pathways such as IFN-α, IFN-γ and allograft rejection were significantly enriched in Paediatric Inflamed and Myeloid Predominant and depleted in Immune Neutral and Immune Desert (Fig. 3A, Additional file 2: Table S7). Mitotic spindle and G2M checkpoint pathways were significantly enriched in Immune Desert indicating a high proliferation rate with low amounts of non-malignant cells. To investigate whether these associations could be recapitulated at the protein level, we leveraged data from 147 samples with matched RNA-seq and proteomic profiles [62]. Compared to Paediatric Inflamed, Immune Neutral and Immune Desert showed a significantly lower average protein z-score in eight immune pathways identified using QuSAGE (pairwise two-sided Student’s t test with Bonferroni correction, p < 0.05, Fig. 3B). These findings reveal active immune pathways in Paediatric Inflamed and Myeloid Predominant consistent with the inflammatory and IFN-γ dominant adult clusters [14].

Fig. 3

Pathway and differential gene expression analysis confirm immune features across immune clusters. A Heatmap showing log-fold change in gene enrichment scores (derived from QuSAGE) in each immune cluster compared to all others in pedNST. Columns show MSigDB pathways with false discovery rate (FDR) < 0.1 in at least three immune clusters. Black box outlines pathways validated by protein in (B). B Ridge plots illustrating average z-scores of proteins [62] involved in immune pathways from A across immune clusters in 141 pedCNS samples (CBTN). Fraction indicates the number of protein readouts from the number of genes in each pathway. Pairwise two-sided Student’s t test with Bonferroni correction, *p < 0.05, **p < 0.01 and ***p < 0.001. Significance levels are shown for comparison to Paediatric Inflamed. C Volcano plots for differentially expressed genes (derived from DESeq2) in each immune cluster compared to other samples. Up- or downregulated genes with absolute log-fold change > 1.5 and FDR < 0.1 are shown in red or blue, respectively. Dashed line shows the log-fold change threshold. Top ten differentially expressed genes are annotated. D Heatmap showing median z-score expression of 59 genes with known immunoregulatory functions across immune clusters in pedNST

We further sought to determine key differential genes underlying each immune microenvironment. Using differential gene expression analysis [94], we found cell-type-specific genes upregulated in Paediatric Inflamed including B-cell (JCHAIN, MZB1 and CD79A), macrophages (MARCO) and granulocyte-specific genes (ADGRG3 and FPR2) (Fig. 3C). Consistently, differentially expressed genes comparing Paediatric Inflamed to other samples included genes specific to B-cells, plasma cells and macrophages (Additional file 1: Fig. S8A). In Myeloid Predominant, we found only one significantly upregulated gene, CHIT1, which encodes for chitotriosidase secreted by active macrophages [95] (Fig. 3C). The lack of other significantly upregulated myeloid-specific genes may suggest heterogeneity among myeloid cell types in different cancer entities. JCHAIN and CXCL10 were downregulated along with MMP9 in Immune Neutral (Fig. 3C). Immune-related genes such as JCHAIN, CD3 chains, IL6 and cytotoxicity genes GZMK and GZMA were significantly downregulated in the Immune Desert (Fig. 3C). Altogether, these results reveal core immunological pathways and genes driving immune clusters.

We sought to determine possible immunoregulatory players in pedNST. Comparison of expression levels of 59 genes with known regulatory functions [96] revealed differential expression across immune clusters after adjusting for cancer type and the total immune infiltrate (Analysis of Covariance (ANCOVA)). Members of CD28 superfamily receptor (PDCD1, CTLA4, BTLA and ICOS) were upregulated in Paediatric Inflamed (p = 6.61 × 10⁻⁹, 9.98 × 10⁻¹², 3.04 × 10⁻¹³ and 1.67 × 10⁻¹⁷, respectively, Fig. 3D). Genes encoding PD-L1 and PD-L2 (CD274 and PDCD1LG2) were upregulated in Paediatric Inflamed (p = 0.03 and 6.34 × 10⁻⁴). Genes involved in regulatory T-cell (Tregs) function and activation (IDO1, STAT3 and FOXP3) were upregulated in Paediatric Inflamed (p = 1.63 × 10⁻⁸, 0.002 and 2.31 × 10⁻¹¹, respectively). These findings reveal high expression levels of immune checkpoint genes as well as genes involved in immunosuppression in Paediatric Inflamed.

In Myeloid Predominant, TGFB3 encoding for TGF-β3 was significantly upregulated (ANCOVA, p = 0.007) (Fig. 3D). LGALS3 encoding for Galectin-3 with negative regulatory functions in macrophages [97] was upregulated in Myeloid Predominant suggesting the prominent macrophage influence in this cluster (p = 0.004). Other myeloid-specific genes with immunoregulatory functions were upregulated in both Paediatric Inflamed and Myeloid Predominant and included IL4I1 [98], CCL2 [99] and FGF2 [100] (Fig. 3D). These findings show preferential expression of myeloid-related regulatory genes in Myeloid Predominant.

In Immune Desert, we found that CD276 encoding immune checkpoint protein B7-H3 to be highly expressed (p = 9.61 × 10⁻⁴, Fig. 3D). Unexpectedly, among immune checkpoint genes, LAG3 was expressed at higher levels in Immune Desert (p = 0.004). Further inspection of immune checkpoint genes in PDX data confirmed high expression of LAG3 (Additional file 1: Fig. S8B). Across medulloblastoma samples in Immune Desert, WNT subgroup expressed LAG3 at significantly higher levels compared to SHH and there was a non-significant trend between WNT and Group3/4 (two-sided Student’s t test, p = 0.02 and 0.08, Additional file 1: Fig. S8C). These results suggest LAG3 expressed by tumour cells in pedNST calling for further validation of its function in the brain and identify B7-H3 as a possible mechanism of immune exclusion in subsets of pedNST.

In neuroblastoma, profiling expression of five immune checkpoint genes revealed ~ 10% of samples had exclusively elevated expression of HAVCR2 gene encoding TIM3 (Additional file 1: Fig. S8D). We confirmed this finding using a non-overlapping neuroblastoma dataset of 209 immune infiltrated tumours. In this cohort, we found 13 cases (6.2%) that showed elevated expression of HAVCR2 and low expression of LAG3 (Additional file 1: Fig. S8D). Immunohistochemical staining for TIM3 and LAG3 proteins using an independent TMA showed that TIM3 was detectable in 14 samples without LAG3 staining (Additional file 1: Fig. S8E). These findings identify an HAVCR2 expression in a subset of neuroblastomas suggesting a distinct mode of immunosuppression.

Immune microenvironment associations with tumour-intrinsic genomic alterations and tumour mutational burden

We next asked whether immune clusters are associated with tumour mutational burden (TMB). We did not find a statistically significant difference in the total number of non-synonymous single-nucleotide variants (SNV) per coding megabase (SNV/Mb) or SNV + Insertion/deletion (Indel) /Mb (SNV + Indel /Mb) across immune clusters (ANCOVA comparing to Paediatric Inflamed, Fig. 4A,B, Additional file 2: Tables S9, S10). However, when we looked at 63 pedHGG samples, we found that Myeloid Predominant and Immune Neutral harboured significantly higher TMB compared to Immune Desert (two-sided rank sum test, p = 0.02, Fig. 4C). Four of 6 cases with > 5 SNV + Indel/Mb and germline variations in MLH1, MSH2, MSH6, PMS2, POLE or POLD1 belonged to Myeloid Predominant, suggesting higher levels of immune infiltration in patients with biallelic Mismatch Repair Deficiency (bMMRD) syndrome. Our findings otherwise revealed no associations between the immune microenvironment and TMB in non-hypermutant pedNST, suggesting TMB, with the exception of MMR-associated tumours, may not be an appropriate biomarker for immune checkpoint inhibitors in this population.

Fig. 4

No associations between TMB or predicted neopeptides across immune clusters. A, B Boxplots showing TMB defined as SNV per megabase (Mb) (A) or SNV + Indel per Mb (B) across immune clusters in pedNST. ANCOVA, not significant. C Boxplot showing TMB for pedHGG samples across immune clusters. Two-sided rank sum test, *p < 0.05. D Boxplot showing number of predicted strong binding peptides (defined as binding affinity ≤ 0.5) for pedNST samples across immune clusters. ANCOVA, not significant. E Boxplot showing number of predicted strong binding peptides for pedHGG samples across immune clusters. Two-sided rank sum test, *p < 0.05. F Heatmap illustrating scaled number of samples (z-score) with at least one non-synonymous SNV/Indels in genes involved in oncogenic pathways, as defined by TCGA. Barplot shows the total number of samples with alterations in each pathway. Stacked barplot shows proportion of tumour types present in samples with alterations in each pathway. Numbers in brackets indicate the number of altered genes in each pathway. G Oncoprint illustrating the distribution of somatic mutations in the top 15 most commonly altered genes in pedNST across the four immune clusters. In all boxplots, boxes show median and IQR and whiskers represent 1.5 times IQR

To elucidate whether antigen presentation may be associated with immune clusters, we investigated the MHC presentation potential of somatic mutations in pedNST, we determined patients’ HLA class I types [65] (Additional file 1: Fig. S9) and performed the mutant peptide extractor and informer (MuPeXI) [66] analysis. We identified 7591 strong binding and 21,680 weak binding peptides, as defined previously [67]. The number of predicted strong binding peptides was highest in neuroblastoma and pedHGG (medians 19 and 17) and lowest in craniopharyngioma (median 2.5), consistent with its low TMB (median 0.07 SNV + Indel/Mb, Additional file 1: Fig. S10A-B). Differences in the number of strong or weak binding peptides did not reach statistical significance after adjusting for cancer type and data source (ANCOVA, Fig. 4D, Additional file 1: Fig. S10C, Additional file 2: Table S11). Consistent with our TMB observation, pedHGG showed a significantly higher number of strong or weak binding peptides in Myeloid Predominant compared to Immune Neutral or Immune Desert (Fig. 4E, Additional file 1: Fig. S10D). While these data indicate neopeptides are not independently associated with immune infiltration across cancer entities in the pedNST, higher numbers of neopeptides in pedHGG cases clustered in the Myeloid Predominant suggest associations with immune infiltration and potentially immune recognition in this cancer entity.

We hypothesized that disruptions in oncogenic pathways may be associated with distinct immune clusters. We identified samples with alterations (SNVs, Indels and fusions) in at least one gene in the ten TCGA oncogenic pathways [37] (Fig. 4F). Tumours with somatic alterations in members of the receptor tyrosine kinase (RTK) pathway were most frequently found in Myeloid Predominant (Cochran-Mantel–Haenszel (CMH) test, p = 6.1 × 10⁻⁵, Fig. 4F). Paediatric Inflamed showed a higher percentage of samples with mutations in the Mismatch Repair (MMR) pathway, although this difference did not reach statistical significance across cancer entities (p = 0.41, Fig. 4F). Among common driver mutations (Fig. 4G) [101], we found that BRAF-altered samples were differentially clustered among immune clusters, with BRAF-KIAA1549 fusion-positive samples grouped primarily in Myeloid Predominant (p = 3 × 10⁻⁶, Fig. 4G). While the causal relationship between tumour-intrinsic genomic alterations and immune microenvironment remains unclear, our results reveal associations of altered molecular pathways and microenvironmental features.

T- and B-cell repertoire analysis suggest antigen presentation and clonal outgrowth in immune infiltrated pedNST

Adaptive immune system consisting of T- and B-cells is a critical determinant of anti-tumour immunity with clonal expansion suggesting antigen exposure. To infer potential immune reactivity across pedNST, we sought to determine the extent of clonal diversity for T- and B-cells. Using an immune repertoire processing framework [68], we recovered a total of 23,842 complementarity-determining region (CDR3β) sequences from 582 pedNST samples. To validate the diversity estimates, we used the CapTCR-seq method [72] to enrich T-cell receptor (TCR) sequences in RNA-seq libraries from adult and paediatric cancer samples (n = 26). TCRβ estimated Shannon diversity inferred from bulk RNA-seq showed the highest correlation with the observed Shannon diversity obtained from CapTCR-seq (Fig. 5A, Additional file 1: Fig. S11A-B). Comparing the diversity estimates and total number of TCRβ reads showed a linear correlation (adjusted r² = 0.72, Fig. 5B). However, 12 samples displayed lower diversity relative to their number of TCRβ reads, suggesting a T-cell clonal expansion in these samples. Conversely, eight samples harboured outlier high diversity indicating several individual clones infiltrating these tumours (Fig. 5B). These findings suggest a method of T-cell diversity inference and identification of potential polyclonal and clonal repertoires using bulk RNA-seq data.

Fig. 5

T-cell repertoire analysis reveals associations of T-cell clonal expansion with B-cell repertoire and immune microenvironment in Paediatric Inflamed. A Scatterplot depicting a linear correlation between TCRβ Shannon diversity estimated from RNA-seq and measured by capturing all TCR sequences from the same RNA-seq libraries in adult and paediatric cancer samples using CapTCR-seq. B Scatterplot showing correlation between TCRβ estimated Shannon diversity and total number of TCRβ reads. Blue line shows fitted linear regression. Red and blue dots represent polyclonal and clonal T-cell repertoires defined as residuals greater than two absolute standard deviations. Circle plots to the right illustrate two examples of polyclonal (top) and clonal (bottom) T-cell repertoires. Each circle is one T-cell clone and circle diameters are proportional to TCRβ reads. C Boxplots showing differences in log-transformed estimated TCRβ Shannon diversity across immune clusters for neuroblastoma (left) and pedCNS (right). Two-sided Student’s t test with Bonferroni correction, *p < 0.05, **p < 0.01 and ***p < 0.001. D Boxplots comparing levels of T-cells, dendritic cells or monocytes, as determined in Fig. 1C, in samples with TCRβ residuals (obtained from the linear regression in B) ≤ 25th percentile and ≥ 75th percentile of Paediatric Inflamed (left) or Myeloid Predominant (right). Two-sided Student’s t test, *p < 0.05, NS: not significant. E Boxplot showing proportion of specific immunoglobulin isotypes in B-cell repertoires across immune clusters in pedNST. Two-sided Student’s t test, *p < 0.05 and **p < 0.01. F Boxplots comparing immunoglobulin clonality (gini index) across immune clusters for neuroblastoma (left) and pedCNS (right). Two-sided Student’s t test with Bonferroni correction, *p < 0.05 and **p < 0.01. G Boxplots comparing levels of gini index or tumour-type normalized expression of IGHG1 or IGHG3 in samples with TCRβ residuals (obtained from the linear regression in B) ≤ 25th percentile and ≥ 75th percentile of Paediatric Inflamed (left) or Myeloid Predominant (right). Two-sided Student’s t test, *p < 0.05, NS: not significant. In all boxplots, boxes show median and IQR and whiskers represent 1.5 times IQR

Using estimated Shannon diversity, we analysed 361 pedNST samples for TCRβ diversity (Additional file 1: Fig. S11C, ‘Methods’). Neuroblastoma samples in Paediatric Inflamed showed significantly higher diversity compared to neuroblastoma samples in Immune Neutral and Desert (two-sided Student’s t test with Bonferroni correction, p = 0.03 and 2.6 × 10⁻⁶, Fig. 5C). Within pedCNS, Paediatric Inflamed had higher diversity compared to Myeloid Predominant and Immune Neutral (p = 3.3 × 10⁻⁷ and < 2 × 10⁻¹⁶, Fig. 5C). To determine microenvironmental changes that co-vary with T-cell diversity independent of T-cell infiltration, we calculated residuals from a linear regression model for TCRβ diversity and reads (Fig. 5B). In Paediatric Inflamed, samples with residuals in lower quartile, corresponding to uneven T-cell repertoire, had significantly higher dendritic cell scores (p = 0.01, Fig. 5D). Conversely, within Myeloid Predominant, samples with residuals in lower quartile harboured significantly higher levels of monocytes (p = 0.03, Fig. 5D). These results reveal associations of myeloid compartment with clonal T-cell repertoire, suggesting the possibility of clonal outgrowth as a consequence of interactions with antigen-presenting cells.

We used a similar immune repertoire inference tool as for the T-cell repertoire analysis [68] and recovered 197,769 unique immunoglobulin heavy chain (IGH) CDR3 sequences. Across pedNST, we found the highest number of IGH isotypes in craniopharyngioma and neuroblastoma samples (Additional file 1: Fig. S11D). Consistent with findings in adult cancers [71], IGHG1 constituted the largest proportion of B-cell repertoire relative to total number of isotypes (Additional file 1: Fig. S11D). We found IGHG1 and IGHG3 were significantly enriched in Myeloid Predominant compared to Immune Neutral and Immune Desert (two-sided Student’s t test, p = 0.03 and 0.008, Fig. 5E), suggesting that potential antibody-producing cells infiltrated these tumours.

To determine the extent of B-cell clonal expansion in pedNST that may suggest antigen recognition [102], we used the gini index as a measure of uneven B-cell cluster distribution in each sample [103] (Additional file 1: Fig. S11E, ‘Methods’). Across neuroblastoma samples, Immune Desert had significantly lower gini index compared to Paediatric Inflamed (medians 0.29, and 0.53, two-sided Student’s t test with Bonferroni correction, p = 0.05, Fig. 5F). Within pedCNS, Paediatric Inflamed had the higher levels of gini index compared to Myeloid Predominant and Immune Neutral (medians 0.36, 0.26 and 0.23, respectively, p = 0.03 and 0.002, Fig. 5F). Within Paediatric Inflamed, samples with clonal T-cell repertoire (residuals ≤ 25th) also had clonal B-cell repertoire and expressed higher levels of IGHG1 and IGHG3 compared to samples with polyclonal T-cell repertoire (residuals ≥ 75th) (p = 0.02, 0.01 and 0.01, Fig. 5G). We did not find this association in Myeloid Predominant that harboured lower levels of T-cells (Fig. 5G), suggesting that both the extent and clonality of B-cell infiltration may contribute to T-cell clonal expansion.

Lymphoid and Myeloid subtypes provide insights into mechanisms of immune evasion

To further characterize T-cell subtypes in the Paediatric Inflamed cluster with high levels of lymphocyte infiltration, we obtained 40 gene signatures of CD4 + and CD8 + T-cell subpopulations from the single-cell profiling of tumour-infiltrating lymphocytes [60]. Consensus clustering revealed five distinct T-cell groups (TGs) within Paediatric Inflamed with notable patterns (Fig. 6A). Nine neuroblastoma samples were characterized by elevated enrichment of naïve T cells (TG2). TG3 consisted of 6 samples with enrichment of Naive and NME1 + T cell signatures while depletion of other T cell subtypes. TG4 (n = 30) showed intermediate enrichment for T-cell signatures, while TG5 had highest enrichment of T-cell signatures related to memory (Tm), effector memory (Tem), T follicular helper/helper (Tfh/h), Tregs and exhausted (Tex) T cells. These data identify subset of samples within Paediatric Inflamed that are infiltrated with naïve (TG2/TG3) or differentiated (TG4/TG5) T cells.

Fig. 6

Analysis of immune-cell subtypes in Paediatric Inflamed and Myeloid Predominant provides insights into mechanisms of immunosuppression. A Heatmap showing scaled gene set enrichment scores for 40 T-cell subtypes as annotated in [60] in the Paediatric Inflamed (Tn: Naïve T-cells, Th: T helper, ISG: Interferon-stimulated genes, Tm: memory T-cells, Trm: Tissue-resident memory T-cells, Tem: effector memory T-cells, Temra: effector memory re-expressing CD45RA T-cells, Treg: regulatory T-cells, Tfh: follicular helper T-cells, Th17: T helper 17, Tex: exhausted T-cells, KIR: Killer immunoglobulin-like receptor). B Boxplots showing log-transformed TCRβ estimated Shannon diversity (left) and TMB (right) across T-cell groups (TG). Pairwise two-sided Student’s t test with Bonferroni correction, *p < 0.05, **p < 0.01 and ***p < 0.001 (left). Two-sided rank sum test, *p < 0.05 (right). C Boxplots showing differences in immune-cell infiltration, as determined in Fig. 1C, across T-cell groups (TG). Pairwise two-sided Student’s t test with Bonferroni correction, *p < 0.05, **p < 0.01 and ***p < 0.001. D Heatmap illustrating median tumour-type normalized expression of selected genes in TG2 and TG5. E Heatmap depicting scaled gene set enrichment scores for 9 myeloid cell subtypes as annotated in [61, 104] (DC: dendritic cell, Mono: monocyte, Mac/Macro: macrophage, TAM: tumour-associated macrophage). F Heatmap showing scaled gene set enrichment scores for microglia [62, 63] in the pedCNS subset of Myeloid Predominant. G Heatmap illustrating median tumour-type normalized expression of selected genes in MG1, MG3 and MG5. In all boxplots, boxes show median and IQR and whiskers represent 1.5 times IQR

In Paediatric Inflamed, Cox proportional hazards model revealed that TG5, highly infiltrated with differentiated T cells, had worse PFS compared to TG1 (HR = 9, 95% CI 1.5–54, p = 0.02, Additional file 2: Table S12). TG5 also showed lower T-cell diversity compared to TG2 and TG3, which were enriched with T cell naïve and proliferating NME1 + T-cell signatures (Fig. 6A) (Student’s t test with Bonferroni correction, p = 0.002, Fig. 6B, left). While TG2 harboured significantly higher TMB compared to TG1 (Student’s t test, p = 0.02), TG5 did not show any significant difference in TMB compared to other groups (Fig. 6B, right). TG5 harboured higher infiltration of T, NK and myeloid cell infiltration suggesting antigen recognition by myeloid cells followed by clonal expansion and T-cell differentiation (Fig. 6C). TG5 showed higher expression of HLA class I genes (HLA-A, HLA-B, HLA-C), inhibitory cytokines and chemokines (CXCL8, CXCL1, IL6), genes encoding transcription factor AP-1 subunits (FOS, FOSB, JUN, JUNB, ATF3), markers of T-cell activation (CD86, CD69) and Toll-like receptors (TLR2, TLR4) (Fig. 6D). Conversely, TG2 showed higher levels of B-cell marker genes (CD19, CD22), immune checkpoint gene (BTLA) and Fc receptors (FCER2, FCRL2) (Fig. 6D) suggesting a B-cell-mediated mechanism of T-cell suppression. These findings suggest two immunosuppression mechanisms in Paediatric Inflamed: suppression of T-cell effector functions (TG5) and inhibition of T-cell activation and differentiation (TG2).

To delineate myeloid cell subtypes in Myeloid Predominant, we performed consensus clustering using single-cell signatures for tumour-infiltrating myeloid cells [61, 104] (‘Methods’). MG1 consisted of 22 samples with high enrichment of mast cells, conventional dendritic cells type 2 (cDC2), monocytes and macrophages (Fig. 6E). MG2 (n = 89) showed low enrichment of dendritic cells but high levels of PLTP-expressing macrophages. MG3 (n = 70) was characterized by elevated enrichment of plasmacytoid and conventional dendritic cells (pDC and cDC1), and tumour-associated macrophages (TAM). MG4 (n = 75) showed depletion of mast and macrophage signatures with moderate enrichment of monocytes and cDC1 signatures, while MG5 (n = 23) had the lowest enrichment of all myeloid signatures (Fig. 6E). Scoring microglia-specific gene signatures [62, 63] in CNS tumours showed pedLGG samples had the highest levels of microglia geneset enrichment, primarily in MG3 (Fig. 6F). These findings reveal two subgroups within Myeloid Predominant, characterized by high levels of mast cells, macrophages and monocytes (MG1) and high levels of pDC and TAMs (MG3).

Cox proportional hazards model adjusting for gender and cohort revealed that low Mast/Macro groups (MG3, MG4 and MG5) had a trend toward improved PFS compared to MG1 (HR = 0.4, 0.38 and 0.37, 95% CI 0.14–1.14, 0.14–1.02 and 0.12–1.12, p = 0.09, 0.05 and 0.08, respectively, Additional file 2: Table S13). In MG1, differential gene expression analysis showed genes associated with cytotoxic T cells (FGFBP2) and leukocyte trafficking (TSPAN8) [105] were expressed at higher levels compared to MG3 and MG5 (Fig. 6G). PENK, encoding an opioid receptor proenkephalin and a precursor for neuropeptides involved in immune-neural crosstalk [106] was expressed at higher levels (Fig. 6G). MG1 also showed high expression of ANGPTL7, a regulator of angiogenesis [107], and F2RL2, encoding proteinase-activated receptor 3 (PAR3), that may trigger myeloid infiltration [108]. IL34, a brain-specific ligand for CSF1R, was expressed at higher levels in MG1 suggesting a role in monocyte-macrophage differentiation in MG1 [109, 110]. In MG3, genes involved in endothelial regulation (ESM1, APLN, ALCAM) and ligands for αVβ3 integrin (EDIL3 and POSTN) were expressed at higher levels compared to MG1 and MG5 (Fig. 6G), the latter of which recruit tumour-associated macrophages (TAM) in adult glioblastoma [111, 112]. Notably, MG3 exhibited higher expression of Histamine receptor H1 (HRH1) and HNMT involved in histamine metabolism (Fig. 6G) suggesting immune evasion and T-cell dysfunction [113]. While further validation is warranted, these findings suggest cellular interactions influencing the myeloid compartment (periostin-integrins, IL34-CSFR1) and immune exclusion and evasion (histamine and collagens) in Myeloid Predominant.

Discussion

We report universal immune microenvironment groups across 12 cancer entities highlighting tumour-agnostic immunological similarities in paediatric primary nervous system tumours. Across our compendium, we found that ~ 72% of samples harboured generally cold immune microenvironments in higher frequency compared to adult tumour microenvironments [14]. Within our cohort, we found that 10% of samples showed an inflamed microenvironment harbouring high levels of lymphoid and myeloid cell infiltration and diverse T- and B-cell repertoire. In this cohort of non-hypermutant cancers, TMB was not independently associated with an inflamed microenvironment. This finding is in contrast to bMMRD cancers where samples with high CD8 + T-cell infiltration had higher SNV/Mb [114]. Although we could not directly study the associations between TMB and ICI response, our findings do not support TMB as a biomarker for use in non-hypermutant paediatric cancers. The dynamic range of cellular interactions in immune-infiltrated tumours, some of which we suggested in the presented association studies, reflects the existing challenge in defining a robust biomarker of anti-tumour immune response in paediatric tumours.

Analysis of immunoregulatory genes and immune cell subtypes provided insights into the ligand-receptor interactions with translational implications. In Paediatric Inflamed, a number of upregulated genes overlapped with those differentially expressed in post-treatment samples collected from melanoma patients responding to nivolumab [115]. These include immune checkpoint genes (PDCD1, TIGIT, CTLA4 and BTLA) as well as genes involved in T-cell cytotoxic functions (IL4I1, FASLG, TNFRSF9, TNFRSF4, CD244, CD27, CD80 and ICOS) and suggest samples in this cluster may be good candidates for ICI therapy. The Myeloid Predominant cluster showed elevated expression of two immunoregulatory genes with myeloid-specific functions: TGFB1 [116, 117] and SIRPA, encoding for SIRPα that negatively regulates phagocytosis [118]. Blocking this interaction with anti-CD47 promotes cell killing in preclinical models [9] and a phase I clinical study, the efficacy testing is ongoing (NCT02216409). Our results suggest disrupting cellular interactions involving the lymphoid and myeloid compartments may be beneficial among a subset of pedNST to reshape the tumour microenvironment and elicit anti-tumour immunity.

Although this study provides insights into immunological aspects of paediatric nervous system tumours, there are important limitations to note. We may have missed specific regions with significant immune infiltration due to biased sampling. The IHC analysis on tissue microarrays may not capture tumour heterogeneity accurately. Immune deconvolution analysis, based on bulk RNA-seq data, can be prone to overestimation due to overlapping gene expression in different immune cell types. Due to this limitation, we could not confirm the cell origin of LAG3 and TIM3 gene expression in the Immune Desert, which may have implications for immune modulation. Higher-resolution techniques like single-cell and spatial RNA and protein analysis can offer a more comprehensive understanding of the tumour microenvironment.

Conclusions

We identified distinct immune microenvironment clusters across paediatric nervous system tumours and proposed microenvironmental mechanisms of immune dysfunction and suppression. With immunotherapy becoming more widely available in the paediatric oncology armamentarium, our findings highlight the value of immunogenomic approaches to guide patient stratification and inform precision oncology programmes.