High-Resolution Cartography of the Transcriptome and Methylome Landscapes of Diffuse Gliomas

Abstract

Molecular mechanisms of lower-grade (II-III) diffuse gliomas (LGG) are still poorly understood, mainly because of their heterogeneity. They split into astrocytoma- (IDH-A) and oligodendroglioma-like (IDH-O) tumors both carrying mutations(s) at the isocitrate dehydrogenase (IDH) gene and into IDH wild type (IDH-wt) gliomas of glioblastoma resemblance. We generated detailed maps of the transcriptomes and DNA methylomes, revealing that cell functions divided into three major archetypic hallmarks: (i) increased proliferation in IDH-wt and, to a lesser degree, IDH-O; (ii) increased inflammation in IDH-A and IDH-wt; and (iii) the loss of synaptic transmission in all subtypes. Immunogenic properties of IDH-A are diverse, partly resembling signatures observed in grade IV mesenchymal glioblastomas or in grade I pilocytic astrocytomas. We analyzed details of coregulation between gene expression and DNA methylation and of the immunogenic micro-environment presumably driving tumor development and treatment resistance. Our transcriptome and methylome maps support personalized, case-by-case views to decipher the heterogeneity of glioma states in terms of data portraits. Thereby, molecular cartography provides a graphical coordinate system that links gene-level information with glioma subtypes, their phenotypes, and clinical context.

Keywords: DNA methylation; gene expression; grade II–IV gliomas; integrative bioinformatics; molecular subtypes; self-organizing maps machine learning; tumor evolution; tumor heterogeneity.

Figure 2

Cartography of ‘all-glioma’ transcriptomes considering pilocytic astrocytomas (PA) and IDH-wt, IDH-A, IDH-O, and neuronal (NL) gliomas of WHO-grade II- IV (see text): (A) The pairwise similarity heatmap reveals two major clusters formed by IDH-wt and IDH-mut (and NL) gliomas, respectively. The IDH-mut subgroup E3 reveals similarities with PA and, to a lesser degree, also with CL and MES IDH-wt (GBM-like) gliomas. (B) The similarity net visualizes mutual similarities between glioma specimen (dots). It separates virtually all subgroups (PA, CL, MES, and PN for IDH-wt and E2–E8 for IDH-mut). In the right part, samples of each group are separately colored in black. (C) Heatmap of expression profiles of selected gene sets indicate up- and downregulation in a subtype-specific fashion (see Supplementary File 1: Figures S1–S4 for details). (D) The SOM maps of the different expression groups ‘portrait’ their expression patterns in terms of up- (in red) and downregulated (in blue) gene clusters. Note that each group is characterized by a unique ‘fingerprint’ expression portrait. (E) The overexpression summary map provides an overview about the spots upregulated in any of the groups, as illustrated by selected expression profiles. They indicate upregulation of the respective genes in a subtype-specific fashion. For example, ‘E3_UP’ assigns a spot that specifically upregulates in group E3 and partly in PA (see also the group portraits of E3 and PA in part D). (F) Box plots of the group-related expression of selected gene sets reveal characteristic effects. While PRC2 targets (and healthy brain functions) loose expression in gliomas compared with the NL subtype, cell cycle activity, immune response, and hypoxia (and EMT, epithelial mesenchymal transition) functionalities gain in expression, especially in IDH-wt gliomas. Interestingly, PA and E3 concertedly change in virtually all situations (see arrows). The gene set maps indicate accumulation of the genes of the respective set (shown by dots) in distinct areas of the map. Gene sets were taken from the literature [41,42,43,57]. (G) The prognostic map colors areas in which upregulation of the respective genes associates with high (red) or low (blue) hazard ratio (HR). The barplots show the composition of groups expressing the respective genes together with the respective overall survival curves (see Figure S7 for details).

Figure 9

Schematic overview: (A) Major characteristics, cell functions, and methylation modes of the LGG types. (B) Cell-level view suggesting that aberrant DNA methylation shapes IDH mut gliomas into a developmental hierarchy while IDH-wt gliomas consist of multiple cellular states. The scheme was partly adapted from [98] on the basis of recent single-cell transcriptomic studies [66,67,99,100] and our mapping of single cell glioma signatures (Figure S11). See text.

Figure 1

Overview about glioma strata and of their analyses: ‘All glioma’ analysis comprises expression data of pilocytic astrocytomas and of grade II–IV gliomas (LGG and GBM), which were stratified as shown in the figure. Then, combined analysis of expression and DNA methylation data of LGG was performed in order to study the effect of DNA methylation on glioma biology. About 85% of LGG carry mutations in the IDH gene, which causes aberrant methylation of DNA and histone side chains via repression of demethylating enzymes by the onco-metabolite 2 hydroxy-glutamate (2HG). E-groups and M-groups of LGG were stratified strictly on the basis of gene expression and DNA methylation data, respectively (see [20]), which, in consequence, gives rise to only partly matched E- and M-groups as well as genetic WHO-classes based on IDH mutation and chromosome 1p and 19q co-deletion status. Abbreviations and color coding of the groups are used throughout the paper.

Figure 3

Cartography of the LGG transcriptome. (A) Mean SOM expression portraits of the E-groups were specific for each class and revealed correspondence with the expression portraits of the M-groups. (B) The LGG transcriptomes divide into 10 major modules of co-expressed genes with (C) characteristic profiles and mostly negative correlation between group averaged expression and gene promoter methylation (the scatter plots show mean expression versus mean methylation averaged over all genes included in the spot and all gliomas per group). (D) The gene set maps plot genes of selected functional context into the expression landscape. Please note their accumulation in different areas, which were assigned to the spots introduced above. (E) The spot-number distribution estimates the degree of the heterogeneity of sample portraits in terms of expressed spot modules and the spot implication maps join spots frequently appearing together (confidence > 0.5) in each of the subtypes by lines.

Figure 4

Cartography of the LGG methylome: (A–E) subgroup portraits, spot summary map, spot profiles, and functions. See legend of Figure 3 for details. The heatmap in part C provides an overview about the methylation spot modules. They were named as follows: anti-GCIMP (spot A’), resembling the CpG-HypOmethylation module of IDH-mutated tumors (CHOP, [29]); GPCR-module (B’); KER-module (C’, keratinization); GCIMP module (D’); GCIMP with specific hypermethylation of IDH-O (GCIMP-O, E’); and hypermethylation of IDH-wt, particularly of the RTK II type (GCIMP-wt, F’) [9]. The spots enrich selected GBM methylation signatures according to [9], which were taken from [29].

Figure 5

Alterations of transcriptional programs and of DNA methylation patterns between subtypes and upon recurrence: (A) differential SOM-portraits of the E- and M-subgroups with respect to the NL subtype revealed alterations of expression and methylation patterns upon glioma development and their functional context. Stratification of portraits of E7 and M6 with respect to the IDH mutation status showed almost identical patterns, meaning that E7 and M6 are suited as reference state, which is dominated by healthy brain characteristics. Differential portraits were also calculated between WHO grades III and II (see also Figure S19) and between primary and recurrent tumors [24]. Cell cycle (CC), oxphos (OX), and/or inflammatory (G) spots gain in most comparisons. (B) Group-averaged expression (GSZ-score) of PRC2 targets as a function of cell cycle activity. Overall, PRC2 targets negatively correlate with cell cycle where effect is largest for IDH-wt (CL and MES subtypes). For IDH-A tumors, one finds a nearly linear decay as indicated by the red line. (C) The summary scheme locates the subtypes in triangular coordinates spanned by ‘archetypic’ cellular functions. It shows selected methylation profiles indicating reduced levels in of the IDH-A subtypes. (D) Mutational load (log number of single nucleotide polymorphisms) per tumor is largest in IDH-wt (E1, M1) and smallest in NL (E7, E8, M1). It increases upon IDH-A progression in the order M4, M3, M2. SNV numbers were taken from TCGA-matched LGG [20].

Figure 6

Relations between expression and methylation: (A) The covariance maps highlight genes co-regulated by mostly anti-correlated expression and methylation values. The frames mark selected spot areas where red and blue color assign opposite alterations of expression and methylation levels as indicated in the figure. (B) Mutual mapping of E- and M-spot genes (red frames) into the M- and E-SOM (arrows) indicates spot ‘melting’ and suggests divergence of genes related to distinct molecular mechanisms with respect to concerted transcriptional activation and methylation changes. Note that capital letters in the SOM-portraits assign spot modules.

Figure 7

Combined view of the expression and methylation landscapes: (A) The network between expression and methylation spots (see Figure S20 for details) was obtained on the basis of gene overlap between the spots and mutual correlations between the spot profiles. (B) The correlation map between expression and methylation spot values revealed negatively correlated M–E spot pairs with considerable gene overlap (see also Figure S20). (C) The network in part A divided into four main coupled regulatory modules in the four consensus subtypes.

Figure 8

Molecular and phenotype maps of the transcriptome (left part of the figure) and methylome (right part) landscapes of LGG. Phenotype maps show associations between gene expression and promoter methylation with the hazard ratio (HR), patient’s age at first diagnosis, sex of the patients, and telomere length ratio between tumor and leukocytes (TLR). Phenotype maps were generated as described in Figure S7 for prognostic maps. The ‘female difference score’ estimates the deviation from the mean percentage of female patients in units of ‘percent-of-percent’. Overall, phenotype maps enable the comparison of expression and methylation levels in the different subtypes with the respective phenotype features.