DIMEimmune: Robust estimation of infiltrating lymphocytes in CNS tumors from DNA methylation profiles.

The interaction of CNS tumors with infiltrating lymphocytes plays an important role in their initiation and progression and might be related to therapeutic responses. Gene expression-based methods have been successfully used to characterize the tumor microenvironment. However, methylation data are now increasingly used for molecular diagnostics and there are currently only few methods to infer information about the microenvironment from this data type. Using an approach based on differential methylation and principal component analysis, we developed DIMEimmune (Differential Methylation Analysis for Immune Cell Estimation) to estimate CD4+ and CD8+ T cell abundance as well as tumor-infiltrating lymphocytes (TILs) scores from bulk methylation data. Well-established approaches based on gene expression data and immunohistochemistry-based lymphocyte counts were used as benchmarks. The comparison of DIMEimmune to the previously published MethylCIBERSORT and MeTIL algorithms showed an improved correlation with both gene expression-based and immunohistological results across different brain tumor types. Further, we applied our method to large datasets of glioma, medulloblastoma, atypical teratoid/rhabdoid tumors (ATRTs) and ependymoma. High-grade gliomas showed higher scores of tumor-infiltrating lymphocytes than lower-grade gliomas. There were overall only few tumor-infiltrating lymphocytes in medulloblastoma subgroups. ATRTs were highly infiltrated by lymphocytes, most prominently in the MYC subgroup. DIMEimmune-based estimates of TILs were a significant prognostic factor in the overall cohort of gliomas and medulloblastomas, but not within methylation-based diagnostic subgroups. To conclude, DIMEimmune allows for robust estimates of TIL abundance and might contribute to establishing them as a prognostic or predictive factor in future studies of CNS tumors.

Introduction

The immune microenvironment is a key factor for tumor growth and progression in various cancers including central nervous system (CNS) tumors.[1,2] There has been an increasing interest in tumor immunology with the avenue of immunotherapies, as there is growing evidence that the tumor microenvironment influences the therapeutical outcome. Tumor-infiltrating lymphocytes (TILs) have been established as an important predictive and prognostic biomarker in several solid tumors, particularly in breast cancer (see[3-5] for review). In neuro-oncology, TILs have been primarily studied in glioblastoma and were shown to be associated with molecular alterations, such as NF1 and RB1 mutations, although the reported prognostic associations are still being controversially discussed.[6-8] Using transcriptomic approaches, immunological differences between molecular subgroups were identified in medulloblastoma (MB), ependymoma and ATRT.[9-11] Recently, profiling the immune microenvironment of over 6,000 primarily pediatric brain tumors (medulloblastomas, malignant rhabdoid tumors, and high-grade gliomas) using MethylCIBERSORT, Grabovska et al. showed associations of particular immune cells with molecular subgroups, mutations, as well as the overall survival in these entities.[12,13]

Different strategies have emerged for the quantification of tumor-infiltrating immune cells. TILs are frequently scored on hematoxylin and eosin (HE) stained slides, which is easily possible in e.g. breast cancer, lung cancer or glioblastoma.[3,4] Recently, machine learning methods have been applied to optimize morphological TIL quantification.[14] However, it becomes increasingly difficult in tumors with morphological similarities to immune cells, like small-round-blue-cells tumors, where immunohistochemical analyses are needed. Also, the high diversity of rare brain tumor entities, spatial heterogeneity, and potential interrater variations make it difficult to obtain solid TIL phenotyping from tissue sections. Utilizing molecular diagnostics such as methylation data to determine the degree of infiltration and subset phenotype of immune cells is therefore appreciated.

Various techniques have been developed to estimate tumor-infiltrating immune cells from bulk molecular profiling data. Deconvolution algorithms based on transcriptomic techniques are most widely used and include approaches based on specific gene signatures for microenvironment cell populations and the CIBERSORT method.[15-17] On the one hand, estimates based on specific immune cell signatures are usually defined as the average expression value of the cell-type-specific genes. However, most published signatures were not optimized for brain tumors resulting in potentially unspecific results. Therefore, we recently optimized these methods allowing for robust estimates of microenvironmental cells in medulloblastoma and pediatric high-grade glioma.[8,9]

On the other hand, CIBERSORT uses a support-vector regression-based deconvolution algorithm to quantify the relative amount of 22 immune subpopulations, which could be used to identify prognostic markers in various cancers.[18] The algorithm provides a p-value for each sample, which is a measure of confidence in the results. In a large gene expression-based meta-analysis of breast cancer, the overall amount of immune infiltration was negatively correlated with this p-value, indicating that lower amounts of tumor-infiltrating immune cells might result in less robust deconvolution results.[19]

These transcriptomic approaches are well established and validated and led to several new insights on the immunological landscape of solid tumors. However, transcriptomic data are only rarely used in routine diagnostics of brain tumors. Conversely, global DNA methylation analysis has been extensively used for molecular analyses of brain tumors,[20-24] and there is an increasing number of applications in the clinical routine for the diagnostic classification of brain tumors and other neoplasms.[25-28] Therefore, methods that can robustly quantify immune cells based on DNA methylation are needed and currently much more widely applicable in neuro-oncology than transcriptomic approaches.

Recently, DNA methylation-based algorithms have been established to study tumor-infiltrating immune cells. Jeschke et al. used differential methylation analysis between lymphocytes and breast cancer to establish a methylation signature (MeTIL) that recapitulates TIL evaluations and their prognostic value in breast cancer.[29] Adapting the CIBERSORT algorithm to DNA methylation data, Chakravarthy et al. introduced MethylCIBERSORT. Using methylation profiles from head and neck squamous cell carcinoma, they divided their data into immune hot tumors with better response and immune cold tumors with worse response to the therapy.[13] As previously described, this approach has been recently applied by Grabovska et al. to study a large cohort of brain tumors.[12] A prerequisite for the reliable immune cell estimation with MethylCIBERSORT is a tumor reference signature for the entity under consideration. Therefore, cell line data have been used, which is, however, not available for rare brain tumor entities. Grabovska et al. therefore used 25 cell lines from rhabdoid tumors and medulloblastoma only as reference.[12]

As nonspecific tumor signatures result in nonreliable deconvolution results, and as the reference data are unavailable for many rare brain tumor entities, we established DIMEimmune (Differential Methylation Analysis for Immune Cell Estimation), a robust DNA methylation-based method for the quantification of TILs (DIME-TIL), CD4+ (DIME-CD4) and CD8+ (DIME-CD8) lymphocytes in central nervous system tumors. The method is not relying on a reference signature of pure tumor cells from the studied tumor entity and can therefore be directly applied to any CNS tumor.

Materials and Methods

All data analyses were performed using the statistical programming language R version 3.6.0[30] with the packages minfi, lumi, missMethyl, ComplexHeatmap, MASS, beeswarm, TCGAbiolinks, MethylCIBERSORT, FlowSorted.Blood.450k and survival.[13,31-39]

Datasets and preprocessing of previously published methylation data

The methylation data used to train our method as well as the data used for comparison with gene expression-based deconvolution algorithms were previously published, analyzed on the Illumina Infinium Methylation 450K Bead Chip and are available from public data repositories. These raw data can be divided into training and validation data (Figure 1a). As reference for immune cells, publicly available methylation profiles of magnetic-activated cell sorted CD4+ and CD8+ lymphocytes from 101 cases were downloaded from Gene Expression Omnibus (GEO).[40,41] Not all samples contain data for both CD4+ and CD8+ lymphocytes, therefore profiles from patients with missing data either for CD4+ or CD8+ lymphocytes were excluded resulting in a dataset of 94 samples. The training cohort for the brain tumor classifier published by Capper et al. was used as a reference for CNS tumors.[25] As we computed differential methylation between bulk tumor methylation profiles and immune cells, samples from reference tissue and tumors with expected high immune infiltration were excluded (methylation classes: LYMPHO, PLASMA, CONTR INFLAM, CONTR REACT, MELAN, MELCYT) resulting in a dataset of 2706 samples from 85 diagnostic categories. As validation datasets, the validation cohort from Capper et al. was used (1104 samples), as well as 763 methylation profiles from medulloblastoma (Cavalli et al.[20]), 162 methylation profiles from atypical teratoid/rhabdoid tumor (ATRT),[23,42] 534 methylation profiles from lower-grade glioma (LGG) (TCGA[43]), 155 samples from glioblastoma (GBM) (TCGA[43]), and 557 samples from ependymoma (EPN) (Pajtler et al.).[22] Clinical annotations were extracted from the supplementary material of the corresponding publications cited above. All DNA methylation profiles were preprocessed with the R-package minfi[31] using single-sample normal-exponential out-of-band (noob) normalization, beta-scores were used for the final analysis.[38] CpG sites associated with single-nucleotide polymorphisms, sex chromosomes, and cross-reactive sites were excluded as previously reported, resulting in a dataset of 428799 CpGs.[25]

Figure 1.

Overview of the data and the DIMEimmune method. a. Data sources and number of samples; b. Overview of the method: the final lymphocyte-specific signature is the intersection of three sets of CpGs. The first two sets are obtained within the framework of differential methylation analysis and the third one by threshold analysis. The first group is selected based on the difference between brain tumor samples, and CD4+ and CD8+ T cells by applying an unpaired t-test. The second group is selected by applying a paired t-test on different immune cell types in the blood. Additionally, a threshold analysis is conducted, which gives the third group of CpGs. The final estimate for lymphocyte infiltration is obtained from dimensionality reduction using PCA; c. Validation: Methylation-based estimates are obtained by applying the learned rotation matrix of PCA on the validation data. Gene expression-based estimates are computed from specific signatures as previously described.[8] Both methods are compared on the same samples using RNAseq/gene expression and methylation-based values in the form of correlation and scatter diagrams

Matching gene expression data were obtained for medulloblastoma, ATRT, glioma and ependymoma from GEO or the European Genome-phenome Archive. Expression data for medulloblastoma (763 samples) were analyzed on the Affymetrix Human Gene 1.1 ST Array and preprocessed as previously described.[9] Expression data for ATRT (88 samples) were analyzed on the Illumina HT12 gene expression array and preprocessed with the R-package lumi.[23,39] Preprocessed RNAseq data (FPKM) for LGG (532 samples) and GBM (64 samples) were obtained from Rhaman et al.[44] Expression data for EPN (129 samples) were analyzed on the Affymetrix HG U133 Plus 2.0 microarray and preprocessed with the R-package affy using the custom chip definition file hgu133plus2hsentrezgcdf (v19.0.0).[45,46]

Immunohistological validation

Parallel methylation and immunohistological analysis were performed on 47 diagnostic cases including pituitary adenomas, gliomas, MBs, EPNs and ATRTs. Formalin-fixed paraffin embedded (FFPE) tissue was obtained from the archives of the Institute of Neuropathology, University Medical Center Hamburg-Eppendorf. Informed consent was obtained for all patients prior to the analysis.

For methylation analysis, DNA isolation was performed on FFPE tissue, 10 × 10 μm sections were cut and DNA isolated using the ReliaPrep™ FFPE gDNA Miniprep System (Promega) according to manufacturer’s instructions. About 100–500 ng DNA was used for bisulfite conversion by the EZ DNA Methylation Kit (Zymo Research). Afterward, the DNA Clean & Concentrator-5 (Zymo Research) and the Infinium HD FFPE DNA Restore Kit (Illumina) were employed to clean and restore the converted DNA. Finally, the Infinium MethylationEPIC BeadChip Kit (Illumina) was used to quantify the methylation status of 850,000 CpG sites on an iScan device (Illumina). Preprocessing was performed analogously to samples from public data repositories. FFPE human tumor samples were used for immunohistochemistry. Analyses were done on an automated Ventana system using anti-CD3 primary antibodies (Zytomed, M3974, 1:100). CD3+ cells were counted in three representative image regions of 2000 × 2000 pixels (at magnification x400) for each sample.

DNA methylation-based estimation of tumor-infiltrating lymphocytes

Three signatures were defined to estimate the amount of tumor-infiltrating CD8+ T cells (DIME-CD8), tumor-infiltrating CD4+ T cells (DIME-CD4), and a mixed signature of tumor-infiltrating lymphocytes (DIME-TIL). To this end, multiple steps of differential methylation analysis were applied to identify specific CpG sites. This was followed by dimensionality reduction using principal component analysis (PCA) resulting in an estimate for the studied cell populations. The details of the method are described in the following.

Differential methylation analysis between immune cells and CNS tumors

First, an unpaired one-sided Welch’s t-test was used to assess the significance of differential methylation between the profiles of 94 CD8+ T cells, 94 CD4+ T cells as well as the combined set of both T cells (188 profiles) and the methylation profiles assigned to the selected diagnostic categories in the training cohort of Capper et al.[25,36,47] Significance of hypomethylation and hypermethylation was computed separately, as we aimed for the identification of CpG sites consistently hypo-/hypermethylated between lymphocytes and tumors across all 85 diagnostic categories. The mean p-value for differential methylation over the 85 statistical tests was used to rank CpGs. Only those CpGs with a mean p-value < 0.05/428799 (Bonferroni correction) were retained for further analysis. This resulted in 373, 387, and 362 hypomethylated as well as 287, 217, and 243 hypermethylated CpGs between immune cells and tumor for CD8+, CD4+, and the mixed set of lymphocytes, respectively. As there was a large overlap between those CpG sites identified for CD4+ and CD8+ T cells, we next used immune cell profiles from the peripheral blood to identify specific CpG sites for each cell type.

Differential methylation analysis between immune cells from the peripheral blood

Methylation data from sorted blood cells using flow cytometric analyses of CD4+, CD8+ T cells, monocytes (CD14+), NK cells (CD56+), B cells (CD19+) as well as neutrophils, eosinophils, (mixed) granulocytes and (mixed) peripheral blood mononuclear cells (PBMCs) from six patients were obtained from Reinius et al. via the R-package flow.sorted.blood.450k.[36,47] In a first step, to obtain specific CpGs for CD8+ and CD4+ T cells, the mean p-value of a paired one-sided Welch’s t-test was computed between CD4+/CD8+ cells and the remaining populations of blood cells in a similar way as for the analysis of immune cells and tumor. Next, we selected those CpGs whose mean beta values were consistently higher/lower in CD4+/CD8+ than in all the remaining blood cell populations (threshold = 0.1). Differentially methylated sites were defined by the intersection of these sets resulting in 601 and 959 hypomethylated as well as 857 and 954 hypermethylated CpGs specific for CD8+ and CD4+ T cells, respectively. Specific CpG sites for lymphocytes were computed analogously comparing the CD4+ and CD8+ profiles with the CD14+, neutrophils, eosinophils and granulocytes using an unpaired t-test resulting in 19393 hypomethylated and 30604 hypermethylated sites.

Definition of specific immune signature and dimensionality reduction

To obtain the final signatures of specifically differentially methylated CpGs, the intersection of the hypo/hypermethylated CpGs from the comparison of immune cells and tumor was built with the hypo/hypermethylated CpG sites comparing the studied immune cell type and the immune cells from the peripheral blood (monocytes (CD14+), NK cells (CD56+), B cells (CD19+), neutrophils, and eosinophils). Next, hypo- and hypermethylated CpGs were combined. This resulted in 44, 30, and 197 specific CpGs for CD8+ T cells (DIME-CD8), CD4+ T cells (DIME-CD4) and tumor-infiltrating lymphocytes (DIME-TIL), respectively.

Based on the assumption that the main variance over the methylation of these specific CpGs is explained by the amount of infiltrating lymphocytes, we applied principal component analysis (PCA) on these CpGs on the brain tumor training dataset for dimensionality reduction and defined the first principal component as a marker for the immune cell population. The amount of immune infiltration on a new dataset can then be estimated by applying the learned PCA representation to the new dataset. The orientation of the first PCA component was selected to correlate with the negative of the mean of the hypomethylated sites. The R code for the computation of DIMEimmune estimates is available as supplementary material.

Statistical and computational analyses, data visualization

DNA methylation-based diagnostic classification of brain tumors was executed with the Heidelberg classifier v11b6.[25] Gene ontology enrichment analyses for significant CpGs were performed with the gometh function.[32] Gene expression-based estimates for T cells and CD8+ T cells were computed as previously described.[8] Briefly, well-established methods for TIL estimation[16,17] were adapted for the use with brain tumors by manually reviewing the signatures and excluding unspecific genes using correlation analyses. The MeTIL score was computed as described by Jeschke et al.[29] MethylCIBERSORT analysis was performed as implemented by Chakravarthy et al. using the brain tumor-specific signature matrix as published by Grabovska et al.[12,13] Gene expression-based/immunohistological TIL scores and methylation-based TIL scores were correlated with Spearman’s rank correlation coefficient and robust linear regression was applied as implemented in the R package MASS.[38] Differences in estimates of TILs between different tumor entities were computed with the Kruskal-Wallis test. Proportional hazards modeling was performed with the R-package survival introducing immune scores as continuous variables.[35] P-values < 0.05 were considered statistically significant.

Results

Definition of specific signatures for tumor-infiltrating lymphocytes

The major principle of our approach is differential DNA methylation analysis between tumor profiles and lymphocytes as well as between different immune cell populations in the peripheral blood (Figure 1). As reference for the different tumor types, 2,706 samples from 85 different entities published by Capper et al. were used.[25] For immune cells, we used large cohorts of profiles from CD4+ T cells and CD8+ T cells, published by Tserel et al.[40] Applying different t-test statistics on brain tumors and immune cells resulted in lymphocyte-specific CpG sites (Selected CpGs 1, Figure 1b, see Methods). As there was considerable overlap between specific sites identified for CD4+ T cells and CD8+ T cells, we next used differential methylation as well as a thresholding analysis on the average methylation value between immune cell populations from the peripheral blood to identify specific CpG sites for CD4+ and CD8+ T cells (selected CpGs 2 + 3, Figure 1b). Furthermore, a mixed TIL signature was computed by performing the differential methylation analysis for the combined set of CD4+ and CD8+ T cells analogously (see methods). The final signature was obtained by intersecting the three selected groups of CpGs. Based on the assumption that the main variance in tumor profiles for the selected CpGs is based on the number of infiltrating lymphocytes, PCA was applied for dimensionality reduction on the selected CpGs in the tumor cohort. The first principal component was defined as an estimate for the tumor-infiltrating immune cells (Figure 1b). To validate our method, we compared our estimated score for immune cells with lymphocyte counts based on immunohistological images and signatures obtained from gene expression analysis. To this end, the PCA obtained in the training phase is applied to the methylation values of the independent test samples for the CpGs of the signature. Hence, the first component is defined as an estimate for the infiltration of the corresponding lymphocyte population (Figure 1c).

Analysis of computed signatures

We defined specific methylation signatures for CD4+ (DIME-CD4) and CD8+ T cells (DIME-CD8) as well as a mixed signature for tumor-infiltrating lymphocytes (DIME-TIL) (Figure 2a-c). These can be used to estimate the amount of infiltration of the corresponding cell type from bulk methylation profiles. Both, specific hypomethylated and specific hypermethylated CpGs exist in the obtained signatures. DIME-TIL contains 105 hypermethylated and 92 hypomethylated CpGs. DIME-CD4 has 4 hypermethylated and 26 hypomethylated CpGs, and DIME-CD8 has 17 hypermethylated and 27 hypomethylated CpGs. To get an overview of all three signatures and compare them to MeTIL, they have been applied to the flow-sorted blood dataset[47] (Figure 2d).

Figure 2.

Visualization of the DIMEimmune signatures using heatmaps. Each colorful dot represents the methylation value of the corresponding CpG in the specified sample; a. DIME-TIL signature for tumor-infiltrating lymphocytes; b. DIME-CD4 signature of CD4+ T cells; c. DIME-CD8 signature for CD8+ T cells. d. Overview comparing estimates based on the three aforementioned signatures and MeTIL on the flow-sorted blood dataset. The DNA methylation values of samples of the same type have been averaged

The methylation values divide the heatmaps into two different vertical areas. One area corresponds to the studied immune cells, and the other to the rest of the cell populations. In DIME-CD4 (Figure 2b), CpGs are the most hypo/hypermethylated for the CD4+ T cells. This is also the case for DIME-CD8 and CD8+ T cells. (Figure 2c). In DIME-TIL, both T cell types are differentially methylated as expected (Figure 2a). Our scores were compared with those from the previously published method by Jeschke et al. (MeTIL), which is also shown in Figure 2d. In contrast to our signatures, MeTIL results in highest scores for the CD19+ B cells, followed by the other lymphocyte populations, showing a different weighting of the lymphocyte subpopulation during TIL estimation. All studied methods were able to differentiate between lymphocytes and CD14+ cells as well as granulocytes.

Next, gene ontology enrichment was performed to identify pathways related to the lymphocyte-specific CpG signatures (Table 1). The most significantly enriched pathways for the DIME-TIL were “regulation of innate immune response” (p = 3.51e-04, FDR = 0.64) and “T cell activation” (p = 5.86e-04, FDR = 0.64). The most significantly enriched pathways for the DIME-CD4 were “regulation of defense response to virus by virus” (p = 1.82e-05, FDR = 0.04) and “positive regulation of interleukin-2 biosynthetic process” (p = 1.19e-04, FDR = 0.13). Finally, the most significantly enriched pathways for the DIME-CD8 were “type I interferon signaling pathway” (p = 1.10e-04, FDR = 0.24) and “interferon-gamma-mediated signaling pathway” (p = 5.14e-04, FDR = 0.45). Most of the enrichment results are not significant after correction for multiple testing. Nonetheless, the most enriched categories are almost all related to the immune system, suggesting that the identified CpGs are related to genes belonging to immune system-related processes.

Table 1.

Enrichment analysis of immune signatures. Results are tabulated for: DIME-TIL signature of tumor-infiltrating lymphocytes; DIME-CD4 signature of CD4+ T cells; DIME-CD8 signature of CD8+ T cells

Cell Type	ID	TERM	DE	P.DE	FDR
TIL	GO:0045088	regulation of innate immune response	3	3.51e-04	0.64
	GO:0042110	T cell activation	4	5.86e-04	0.64
	GO:0060337	type I interferon signaling pathway	4	1.61e-03	1.00
	GO:0030217	T cell differentiation	3	4.00e-03	1.00
	GO:0031295	T cell costimulation	3	8.40e-03	1.00
	GO:1900017	positive regulation of cytokine production involved in inflammatory response	2	9.34e-03	1.00
	GO:0001816	cytokine production	2	9.54e-03	1.00
	GO:0030101	natural killer cell activation	2	9.86e-03	1.00
	GO:0043551	regulation of phosphatidylinositol 3-kinase activity	2	1.12e-02	1.00
	GO:0048535	lymph node development	2	1.16e-02	1.00
CD4	GO:0050690	regulation of defense response to virus by virus	3	1.82e-05	0.04
	GO:0045086	positive regulation of interleukin-2 biosynthetic process	2	1.19e-04	0.13
	GO:0006953	acute-phase response	2	4.76e-04	0.35
	GO:0050829	defense response to Gram-negative bacterium	2	1.20e-03	0.51
	GO:0033572	transferrin transport	2	1.27e-03	0.51
	GO:0034097	response to cytokine	2	1.56e-03	0.51
	GO:0042102	positive regulation of T cell proliferation	2	1.63e-03	0.51
	GO:0050731	positive regulation of peptidyl-tyrosine phosphorylation	2	7.97e-03	1.00
	GO:0033674	positive regulation of kinase activity	1	1.72e-02	1.00
	GO:0046039	GTP metabolic process	1	1.73e-02	1.00
CD8	GO:0060337	type I interferon signaling pathway	3	1.10e-04	0.24
	GO:0060333	interferon-gamma-mediated signaling pathway	3	5.14e-04	0.45
	GO:0043551	regulation of phosphatidylinositol 3-kinase activity	2	6.14e-04	0.45
	GO:0061512	protein localization to cilium	2	9.56e-04	0.52
	GO:0042110	T cell activation	2	3.09e-03	1.00
	GO:0048701	embryonic cranial skeleton morphogenesis	2	3.82e-03	1.00
	GO:0002479	antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-dependent	2	8.92e-03	1.00
	GO:0060071	Wnt signaling pathway, planar cell polarity pathway	2	1.59e-02	1.00
	GO:0019835	cytolysis	1	1.81e-02	1.00
	GO:0032897	negative regulation of viral transcription	1	1.97e-02	1.00

Validation and comparison of immune cell estimates with immunohistological results

As a first validation, methylation-based scores for infiltration of T cells (DIME-TIL, MeTIL, and the sum of the T cell signatures of MethylCIBERSORT) were compared to T cells counts from immunohistological images in a series of 47 diagnostic cases from our institution (Figure 3, Methods). Examples of cases with low/absent, intermediate, and high infiltration of CD3+ cells are shown In Figure 3a-c. The DIME-TIL score showed the strongest correlation with the number of CD3+ cells (R = 0.74, p = 2e-09), followed by MeTIL (R = 0.32, p = 0.029). The MethylCIBERSORT T cell signatures also showed a positive correlation with the number of CD3+ cells quantified using immunohistochemistry (R = 0.26), but the correlation was not significant.

Figure 3.

Immunohistological validation of methylation-based TIL estimates. Parallel immunohistological analyses for CD3 as well as methylation analyses were performed in a cohort of 47 brain tumors. Examples of tumors with low/absent CD3+ cells; a. anaplastic astrocytoma IDH-mutant, WHO grade III, intermediate CD3+ cells; b. rosette forming glioneuronal tumor, WHO grade I, and high numbers of CD3+ cells; c. glioblastoma IDH wild-type, WHO grade IV. Immunohistological counts of CD3+ cells are compared to DIME-TIL (d), MeTIL (e) and the sum of the T cells signatures analyzed by MethylCIBERSORT (f)

Validation and comparison of immune cell estimates with gene expression-based results

As a second benchmark, we used results from RNAseq and gene expression analysis. The validation data sets included 5 types of tumors: 763 medulloblastomas (MB), 129 ependymomas (EPN), 88 atypical teratoid/rhabdoid tumors (ATRT), 532 lower-grade gliomas (LGG), and 64 glioblastomas (GBM), for which both methylation and gene expression/RNAseq data were available. The remaining samples, with methylation data only, could not be used for this analysis. Our approach was also compared to the results for the two previously published methylation-based methods (MethylCIBERSORT and MeTIL). Whereas MeTIL provides only a global TIL score, MethylCIBERSORT can estimate several immune and stromal subpopulations. The comparative analysis for DIMEimmune and MeTIL is shown in Figure 4a.

Figure 4.

Comparison of DIMEimmune with previous studies for different tumor types (GBM, LGG, MB, ATRT, and EPN). Here, the gene expression/RNAseq-based results are used as the benchmark. In the diagrams, the dots are the samples, for which both gene expression/RNAseq and methylation data are available. The y-axis represents the gene expression/RNAseq-based estimates and the x-axis contains the methylation-based estimates; a. Comparison of DIME-TIL and MeTIL with gene expression/RNAseq-based estimation of T cells; b. Comparison of DIME-CD4 and CD4+ T cells estimates obtained from MethylCIBERSORT (MCS) with gene expression/RNAseq-based estimation of CD4+ T cells; c. Comparison of DIME-CD8 and CD8+ T cells estimates obtained from MethylCIBERSORT with gene expression/RNAseq-based estimation of CD8+ T cells

The correlation of DIME-TIL with the gene expression-based T cell signature (Figure 4a) is best for GBM (R = 0.72, p < 2.2e-16), followed by LGG (R = 0.56, p < 2.2e-16), ATRT (R = 0.53, p = 2.3e-6), EPN (R = 0.28, p = 0.0012) and MB (R = 0.13, p = 0.00032). For MeTIL, the correlation is also best for GBM (R = 0.6, p = 3.8e-7), followed by LGG (R = 0.5, p < 2.2e-16), ATRT (R = 0.42, p = 3e-04) and EPN (R = 0.21, p = 0.019), whereas it is not significant for MB (R = 0.042, p = 0.25). Overall, correlation is superior for all 5 tumor entities for DIMEimmune based estimates compared to MeTIL based estimates.

Next, DIME-CD4 and CD4 estimates from MethylCIBERSORT (“CD4_Eff” + “Treg”) were compared to the gene expression-based estimation of T cells (Figure 4b). For DIME-CD4, correlation is again best for GBM (R = 0.5, p = 3.3e-5), followed by LGG (R = 0.42, p < 2.2e-16) and ATRT (R = 0.41, p = 0.00047) and lowest for MB (R = 0.072, p = 0.045), whereas it is positive, but not significant for EPN (R = 0.092, p = 0.3). Correlation between MethylCIBERSORT CD4 and gene expression-based T cells estimates are best for ATRT (R = 0.25, p = 0.036) followed by LGG (R = 0.19, p = 7.8e-06) and not significant for GBM (R = 0.19, p = 0.14), MB (R = −0.0034, p = 0.93) and EPN (R = 0.084, p = 0.35).

Finally, methylation-based estimates for CD8 were compared with gene expression-based estimation of CD8+ T cells (Figure 4c). DIME-CD8 showed highest correlation in GBM (R = 0.52, p = 1.6e-5), followed by ATRT (R = 0.39, p = 0.00093), LGG (R = 0.33, p = 3.7e-15) and EPN (R = 0.25, p = 0.0039), whereas it was not significant for MB (R = −0.071, p = 0.051). Correlation for MethylCIBERSORT CD8 estimates were best in ATRT (R = 0.46, p = 6.4e-5). Correlation in LGG (R = 0.085, p = 0.049) was weak and correlation in MB (R = −0.08, p = 0.027) was negative. Correlation in GBM (R = 0.18, p = 0.15) and EPN (R = 0.15, p = 0.085) was not significant.

Overall, our method resulted in an increased correlation with well-established gene expression-based results compared to previously published algorithms for all brain tumor types except for the MethylCIBERSORT estimate of CD8+ T cells in ATRT.

Lymphocytic infiltration in brain tumors

Next, we used DIME-TIL estimation on the validation set of the Heidelberg Brain Tumor classifier (n = 1,104). Although larger, the training set was not used to avoid bias. As expected, the highest TIL score was found in CNS lymphoma (LYMPHO, Figure 5a). Further, there were high TIL scores in inflammatory and reactive control tissue (CONTR, REACT and CONTR, INFLAM). Other entities with high TIL scores included mesenchymal glioblastoma (GBM, MES), MYC ATRT, anaplastic pilocytic astrocytoma (ANA PA), melanoma (MELAN), chordoma (CHORDM) and pituitary adenomas (PITAD, TSH). Low TIL scores were found in posterior fossa B ependymoma (EPN, PF B), paraganglioma (PGG, nC) and group 4 medulloblastoma (MB, G4) (Figure 5a).

Figure 5.

Clinical application of immune cell estimation. The estimated score of tumor-infiltrating lymphocytes for different subgroups of the Capper et al. validation data set (a), the TCGA glioma samples (b), medulloblastoma samples (c), ATRT (d), and ependymoma samples (e). The p-values are calculated based on the Kruskal-Wallis test. f-i.: Survival analysis for the prognostic relevance of DIME-TIL estimates for the same data as in panel b-e in the form of forest plots. HR: hazard ratio, CI: confidence interval

To get a more detailed insight into the distribution of TILs in different brain tumor subgroups, we applied our method to the validation data of LGG, GBM, MB, ATRT and EPN. Methylation-based diagnoses were computed with the Heidelberg Brain Tumor classifier.[25] As WHO grade III and IV tumors were present in the LGG and the GBM dataset, both datasets were combined to a glioma dataset. The differences in TIL scores between the 6 glioma entities were highly significant (Figure 5b, p = 1.6e-36). Estimated infiltration was lowest in IDH (Isocitrate dehydrogenase) mutated oligodendroglioma and in IDH mutated astrocytoma and highest in mesenchymal glioblastoma. For medulloblastoma, although statistically significantly different between molecular subgroups (p = 1.7e-49), the TIL score was overall low. The lowest number was found in Group 4 medulloblastoma and the largest number in Group 3 (Figure 5c). ATRTs showed prominent lymphocytic infiltration. The ATRT MYC subgroup had significantly more TILs than TYR (intermediate) and SHH (lowest) ATRT (p = 0.00012, Figure 5d). The TIL scores were significantly different between ependymoma subgroups (p = 4.4e-0.8, Figure 5e) with the lowest scores in the YAP and the PF B subgroup and the highest scores in the PF A and the RELA subgroup.

Associations of tumor-infiltrating lymphocytes with survival

Finally, we investigated the associations between the TIL estimates and overall survival in glioma, MB, ATRT and EPN. In the overall glioma cohort, the DIME-TIL score was a strong negative prognostic factor (p = 6.06e-15, HR = 6.26). However, within methylation-based diagnostic subgroups, there was no significant association with survival (Figure 5f). Similarly, the DIME-TIL score was also a negative prognostic factor in medulloblastoma (p = 0.003, HR = 4.3), but not within methylation-defined subgroups (Figure 5g). There were no significant associations with survival in ATRT, likely due to low sample size compared to the other tumor types (Figure 5h). TIL scores were also not associated with survival in ependymoma and the studied subgroups, for which sufficient amount of survival data was available (Figure 5i).

Discussion

Although tumor-infiltrating lymphocytes in the CNS have been extensively studied using conventional immunohistological techniques (see Bienkowski and Preusser for review[5]) or gene expression-based analyses (e.g.[8,9,11,18,48]), reports using DNA methylation-based methods are still relatively rare. In contrast, DNA methylation has been central to brain tumor research in the last decade, and diagnostic tumor classification using DNA methylation data has been established in the diagnostic routine at several institutions.[49,50] This makes DNA methylation data widely available in the clinical routine in neuro-oncology. Although there are well-established techniques for immune cell estimation using, e.g., bulk gene expression data, we aimed for the development of a robust and user-friendly method, which can be used on DNA methylation data.

Specific reference signatures are key for quantification methods of immune cell signatures from DNA methylation data. Approaches based on MethylCIBERSORT do not only require specific signatures for immune cells, but also the tumor entity under consideration.[13] For the latter, DNA methylation data from cell lines are used.[12,13] While these data are widely available for common cranial and extra-cranial solid tumors, there are only few well-established cell lines for rare brain tumor entities. Further, ex vivo cell lines might have different methylation profiles from tumors in vivo. To allow for broader applicability of the method, we opted for differential methylation analysis using bulk data from brain tumors, which is far more widely available than cell lines.[25]

Overall, we established three methods to estimate CD8+ T cells (DIME-CD8), CD4+ T cells (DIME-CD4) and tumor-infiltrating lymphocytes (DIME-TIL). Several steps of analysis were performed to obtain specific signatures for the cell type under consideration. Although several CpGs were identified as lymphocyte-specific, there were only few showing specific hypo/hypermethylation in CD4+ and CD8+ cells (Figure 2), underlining the difficulty to precisely separate lymphocyte subpopulations based on methylation data.

First, we used the correlation between the TIL scores of methylation-based methods (DIMEimmune, MeTIL and MethylCIBERSORT) and lymphocyte counts based on immunohistological images to validate and compare the efficiency of our method. Second, we used correlation analysis of TIL estimates from gene expression and methylation data as validation. These datatypes were selected as paired gene expression and methylation profiling is available for large cohorts of brain tumors. Overall, the DIME-TIL approach showed improved correlations with gene expression-based signatures compared to MeTIL and the implementation of MethylCIBERSORT, which was used for a pan-central nervous system cancer analysis by Grabovska et al.[12]

The MeTIL algorithm has been established using a methodologically related approach for immune cell estimation in breast cancer.[29] While it shows robust results for breast cancer, our approach gives better results for brain tumors. The correlation analysis for medulloblastoma in Figure 4a shows a bimodal distribution of the immune infiltration score, which is likely due to a nonspecific signature in medulloblastoma. We observed similar results in a previous gene expression-based study of immune infiltrates in medulloblastoma, where published lymphocyte signatures developed using extracranial tumors as reference contained genes linked to embryonal brain development.[9] This highlights the need for a brain tumor-specific TIL quantification algorithm, as algorithms and reference signatures developed for extracranial tumors might result in unspecific results in the CNS.

The MethylCIBERSORT algorithm showed mixed results in the correlation analysis with gene expression data and was outperformed by DIME-CD8 and DIME-CD4 for all comparisons except for the estimation of CD8+ T cells in ATRT. For both CD4+ and CD8+ T cells the best performance of MethylCIBERSORT was in ATRT. The tumor reference signature used by Grabovska et al.[12] was composed solely of cell lines from medulloblastoma and rhabdoid tumors, but did not contain tumor cell lines from other brain tumors. This is likely to be the cause for better performance in ATRT than in gliomas, where the correlation is weak and highlights the need to include tumor-specific references. In their work, Grabovska et al. report immune cell contents from 75% to 100% for several brain tumors, in particular low-grade gliomas, which further underlines that tumor cells recognition by the algorithm is compromised for certain brain tumor entities due to unspecific tumor reference signatures. Further, the performance of CIBERSORT-based approaches might be compromised in tumors with low overall immune infiltration, as it has been previously shown for breast cancer.[19] Therefore, our approach might be more robust for the estimation of TIL abundance in CNS tumors, in particular for rare entities, for which there are only few or no cell line data available.

All algorithms showed only weak correlations with gene expression-based approaches in medulloblastoma, in particular for CD8+ T cells. Only the correlation reported for DIME-TIL (R = 0.13, p = 0.00032) and DIME-CD4 (R = 0.072, p = 0.045) were significant. Medulloblastoma has rather low numbers for tumor-infiltrating immune cells as seen here (Figure 5a) as well as in other studies.[9,12] This might result in an unfavorable signal-to-noise ratio, which might be particularly pronounced for methylation analysis due to bimodal distribution of the underlying data. Therefore, immunohistology, gene expression (in particular RNA-sequencing) or single cell sequencing-based methods might be more suitable to study tumors with very low amounts of infiltrating immune cells. Nonetheless, methylation-based analysis keeps the advantage of a much broader application field.

Overall, estimates of TIL abundance based on DIME-TIL were in line with previously reported findings. The high scores in CNS lymphoma as well as reactive and inflammatory tissue can be regarded as a validation of our method. Medulloblastoma showed rather low scores for tumor-infiltrating immune cells as previously reported.[9,51] However, the methylation-based analysis did not find larger amounts of T cells in the SHH subgroups as identified by gene expression analysis,[9,52] possibly due to a high noise to signal ratio.

Lower TIL scores in lower-grade gliomas as well as in IDH mutated compared to IDH wild-type gliomas are in line with previously reported results.[53-55] Further, the mesenchymal subgroup is well known to be particularly highly infiltrated by TILs.[6,8] The larger estimates of TILs in MYC ATRTs compared to SHH and TYR ATRT has also been previously reported.[10,56] PF A ependymoma are among the subgroups of ependymoma with the largest estimated number of TILs, which is well compatible with the inflammatory phenotype that has been shown to be a key feature of this subtype.[11,57]

Survival analysis showed a significant negative prognostic effect of the DIME-TIL score in the overall glioma and the medulloblastoma cohort, but not within methylation subgroups. In glioma, this effect may be caused by increasing numbers of TILs in more aggressive WHO grade IV forms of glioma.[55] Similarly, the difference in TIL estimates in medulloblastoma subgroups, which are highest in the poor prognosis Group 3 medulloblastoma might also act as a confounder. Overall, there are conflicting reports in the literature on the prognostic effects of tumor-infiltrating immune cells in glioma and medulloblastoma.[6-9,12,51,55,58] These results emphasize the need for studies on the prognostic role of TILs in brain tumors in large cohorts of well-defined diagnostic subgroups, as the histological or molecular subtype might otherwise confound the results. To this end, methylation data-based methods offer a unique opportunity, as such data are very widely available.

To conclude, we established a method for the estimation of tumor-infiltrating lymphocytes in CNS tumors from methylation data based on differential analysis. The method showed better performance than previous methods taking lymphocyte counts based on immunohistological images and gene expression-based TIL estimates as reference. Estimates are most robust in tumors with more pronounced lymphocytic infiltrates. As it can be applied to any brain tumor entity, it can contribute to the identification of TILs as a prognostic or predictive factor in oncoming studies involving methylation data from tumors of the central nervous system.

Click here for additional data file.