Available GBM datasets with DNA methylation data measured from brain tissue were systematically collected from main public repositories [The Cancer Genome Atlas (TCGA), NCBI Gene Expression Omnibus (GEO) repository, GEOmetadb and ArrayExpress], more than 600 samples with clinical information were obtained. In order to eliminate platform heterogeneity, only those data derived from Illumina Human Methylation 450 array (GPL13534) were kept in this study. Patients with an overall survival (OS) time ≥0.5 month were kept for survival analysis. At last we achieved two satisfactory datasets as TCGA_450k (n = 138, median age is 60.77 years) and GSE60274_450k (n = 62, median age is 51.5 years). Furthermore, 97 samples in the TCGA_450k cohort have clinical indicator information about MGMT, G_CIMP, IDH1, and Treatment strategy (Supplementary Table S1), based on which the Cox proportional hazards regression analysis was carried out for the newly identified methyl-probe signature.

Illumina Human Methylation 450 array contains 485,577 probes, which are derived from different genomic regions of 9,988 genes. Methylation level of each probe was represented as β value, ranging from 0 (completely unmethylated) to 1 (completely methylated). Preprocessing was performed on the TCGA_450k cohort to eliminate probes with null value, and 382,452 probes with methylation level (0 < β < 1) were kept for downstream analysis. Prognostic signatures were constructed via a series of processes. Firstly, univariate Cox proportional hazards regression analysis was carried out for each probe in the TCGA_450k cohort. Significant methylation probes were chosen with the adjusted P-value cutoff (P < 0.05, adjusted by false discovery rate). 35,708 probes without null values were kept. Secondly, 62 known genes associated with glioma were found in both tumor gene list of Catalogue of Somatic Mutations in Cancer (COSMIC) and glioma gene list of Online Mendelian Inheritance in Man (OMIM) (Supplementary Table S1 – Sheet S3). 4,366 out of 5,390 probes for the 62 genes were without null value in TCGA_450k cohort and used. Thirdly, the overlapping of significant methylation-and 62 gene-related probes were selected, and dispersion statistic with boxplot method was created based on the hazard ratio (HR), outliers that lie beyond the extremes of the whiskers were excluded. Finally, the HRs were used to build the percent weighted coefficients, and probes with weight coefficients over 0.001 were kept to construct survival risk classification model, which was the methyl-probe signature with 111 probes. The median of the methyl-probe score from the training dataset, TCGA_450k, was set as the cutpoint for stratifying high- and low-risk GBM tumors. Validation phase was performed in both TCGA_450k cohort and an independent cohort of GSE60274_450k.

All statistics were performed in R (version 3.4.4). Survival package (version 2.42) (Therneau and Grambsch, 2000) was used for survival risk comparison. Chi-square test was used to evaluate the relationship between methyl-probe and Gender, MGMT, G_CIMP, IDH1, and Treatment. Two-sided t-test was performed for age, PFS_day and os_month. Univariate and multivariate Cox proportional hazards regression were tested on Age, Gender, MGMT status, G-CIMP status, IDH1 mutation, Treatment, Three-CpG panel, Six-CpG panel, and Methyl probe signature. Descriptive statistics with boxplot method was carried out for eliminating outliers of methylation probes. Discovery of significant differential methylation probes was based on the β value of Illumina Human Methylation 450 array by Empirical Bayes (Limma package, version 3.34) between high- and low-risk GBM tumor groups. 1,768 genes of the 3,000 most significant methylation probes (Supplementary Table S3) were determined by the Illumina Human Methylation 450 platform information, and GO analysis was carried out by MetaScape¹.

A total of 62 glioma-tumor related genes were identified and their methylation probes were used for clustering 138 GBM patients. The group of 46 patients with shorter OS, had a non-significant tendency in lower promoter methylation of MGMT or G_CIMP, and a trend but not significant in fewer mutations of IDH1 (Figure 1A). This indicated that methylation loci had the potential to act as effective biomarkers to evaluate the prognosis for GBMs.

To generate the optimal representative methylation probes, genome-wide methylation arrays were used. 449,869 probes with none significant hazard (P > 0.05) were filtrated after univariate Cox analysis, and 35,708 prognostic methylation probes were overlapped with 4,366 ones associated with those 62 glioma-tumor genes, the remained probes were used for further odd probe elimination (Figure 1B). The per centum of hazard value for each remained probes of TCGA_450k was calculated after boxplot filtration. Only the weight value >0.001 were left as prognostic coefficient, and 111 methyl-probes constituted the joint-score formula with the optimal cutoff of 0.4386 (the median risk value from the discovery cohorts) for stratifying low-risk and high-risk patients (Figure 1C, Supplementary Figure S1, and Supplementary Table S2). The methyl-probes covered 60 glioma-tumor genes, 34 of which were involved in the cancer pathway and 20 genes participated in the process of Cell fate or TP53/MYC regulation signaling (Figure 1D and Supplementary Table S2).

The prognostic efficiency of the 111 methyl-probes was validated for GBM risk stratification. Patients were divided into high-risk group (with low joint-score, n = 69) and low-risk group (with high joint-score, n = 69) in the discovery cohort of TCGA_450k, where high-risk patients were consistently associated with shorter OS and progression-free survival (PFS) than low-risk ones (Figures 2A,B and Table 1). Besides, well-known molecular marks, G_CIMP and IDH1, showed certain interrelationships with the methyl probe signature. And the unmethylated G_CIMP and wild IDH1 tumors contributed to low methyl score and relative high prognostic risk. Compared to these molecular markers in NCCN guideline (G_CIMP and IDH1) and previous reported prognostic signatures, the 111 methyl-probes showed more efficiency in glioblastomas risk clarification (Supplementary Figure S2). The joint signature was further tested in an independent cohort of heterogeneous population, GSE60274_450k, and it accurately stratified high-risk and low-risk patients with chemotherapy or radiation or combined treatments (Figure 2C), indicated its generality of predicting power for GBM-specific survival progression.

Patient characteristics of stratified cohorts supported that the methyl-probe signature was not only a predictive indicator for part of the clinical index, but also a prognostic factor for GBMs (Table 1). Within TCGA_450k samples, univariate Cox regression model revealed that age, IDH mutation status, Three-CpG panel, and the methyl-probe signature were significantly correlated with OS (Table 2). MGMT promoter and G-CIMP methylation status showed the risk effect but not significant. Meanwhile, the six-CpG panel failed to show risk relation to the TCGA_450k GBM tumors. Multivariate Cox model further demonstrated that the methyl-probe signature was an independent and protective prognostic indicator (Table 2).

All the GBM patients were divided into high- and low-score groups according to the median of methyl-probe joint-score, and differential methylation analysis was carried out between these two groups. Results demonstrated that 8,361 out of total 9,988 genes covered by Illumina Human Methylation 450 array were in differential methylation status, and 8,153 differential genes were in hypermethylation status in the high-score signature group (Figures 3A,B). In the 111 methyl probes, 45 genes black colored in Supplementary Table S2, were in the significant hypermethylation status, which implied that the methyl-probe signature was greatly relevant to the major subgroup of methylation differential genes. The low-score signature group with high hazard risk significantly had overall low level of methylation in most methylation differential genes, which implied that hypomethylation of the methyl-probe signature accompanied with relatively high risk in GBMs. Gene ontology analyses performed on 1,768 genes from 3,000 methylation differential probes between high and low risk GBM groups. Among these genes, NGFR, PRKN, STAT3, IDH2, etc. were associated with the neuronal system and cancer pathway as the significant terms (P < 0.001) (Figure 3C and Supplementary Table S3). This indicated that methyl-probe signatures have great potential to play as new biomarkers in clinical diagnosis and prognostic evaluations for GBM.