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

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. Supplementary Information The online version contains supplementary material available at 10.1186/s13073-023-01219-x.


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.

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.bedfiles 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.bedfile 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.

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 nonimmune cell populations.(2) We downloaded the human protein atlas data v20.1 (https:// www.prote inatl as.org/ about/ downl oad).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 nonimmune 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 nonimmune 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 2 -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 2 O 2 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 CAVAT-ICA 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.profy le.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.

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 RNAseq and whole exome sequencing (WES) data [23].

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]: (See figure on next page.)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.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.
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, twosided 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.
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.14and 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.

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].
We further sought to determine key differential genes underlying each immune microenvironment.Using differential gene expression analysis [94], we found celltype-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.
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 −4 , 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.
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 −5 , 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 −6 , 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 RNAseq showed the highest correlation with the observed Shannon diversity obtained from CapTCR-seq (Fig. 5A, Additional file 1: Fig. S11A-B).Comparing the diversity 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 estimates and total number of TCRβ reads showed a linear correlation (adjusted r 2 = 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.
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 −6 , Fig. 5C).Within pedCNS, Paediatric Inflamed had higher diversity compared to Myeloid Predominant and Immune Neutral (p = 3.3 × 10 −7 and < 2 × 10 −16 , 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.
To delineate myeloid cell subtypes in Myeloid Predominant, we performed consensus clustering using singlecell 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).

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 Tand 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 immuneinfiltrated 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.

Fig. 2
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

Fig. 3
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

(
See figure on next page.)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 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.