Epitranscriptomic analysis reveals clinical and molecular signatures in glioblastoma

Abstract

This study characterizes the glioblastoma (GB) epitranscriptomic landscape in patient who evolve to progressive disease (PD) or pseudo-progressive disease (psPD). Novel differences in N6-Methyladenosine (m6A) RNA methylation patterns between these groups are identified in the first biopsy. Retrospective data of patients that were eventually deemed to have progressive disease or pseudoprogressive disease was captured from the electronic health record, and RNA from the first resection specimen was utilized to evaluate N6-methyladenosine (m6A) biomarkers from FFPE samples. Molecular analysis of m6A methylation modified RNA employed ACA-based RNase MazF digestion. After Quantitative Normalization with ComBat to mitigate batch effects, we identifed differentially methylated transcripts and gene expression analyses, co-expression networks analyses with WGCNA, and subsequently performed gene set GO and KEGG enrichment analyses. Enrichments for metabolic biological processes and pathways were identified in our differential methylated transcripts and select module eigengene networks highlighted key co-expressed genes intricately tied to distinct phenotypes/traits in patients that would ultimately be deemed PD or psPD. Our study identified key genes and pathways modified by m6A RNA methylation associated with cell metabolism alterations, highlighting the importance of understanding m6A mechanisms leading to the oncometabolite accumulation governing PD versus psPD patients. Furthermore, these data indicate that epitranscriptomal differences between PD versus psPD are detected early in the disease course.

Keywords: Epitranscriptome; Glioblastoma; Novel enhancement; Progression disease; Pseudo-progression.

Fig. 1

Methodology workflow and Kaplan-Meier survival analysis comparing glioblastoma patients. (A) The study design illustrates the clinical and imaging diagnostic process and patient monitoring, highlighting the specific points of biopsy collection and treatment initiation. This shows the time points at which our two study cohorts were sampled. Additionally, it delineates the blood testing and molecular marker assessment points utilized in this study, along with the clinical data. (B) The methodology workflow to access the m6A RNA methylation patters of this study., begins with retrospective Glioblastoma patient data collection, categorizing patients into progression disease (PD) and pseudo progression disease (psPD). FFPE tissue specimens was used for molecular analysis. RNA extraction, quality assessment, and m6A RNA methylation site identification was performed using the m6A Single Nucleotide Array Service using the MazF enzyme to digest the unmethylated (ACA) sites on genes and labeling with Cy5 fluorophore (red) in the microarray analysis, while the methylated sites was undigested (m6ACA) by MazF enzyme are labeling with Cy3 fluorophore (green). The data of Cohort B and Cohort A cohorts pass underwent normalization, batch adjustment, and analysis of m6A variance. Further analysis involved differential gene expression, methylation, and exploration of site-specific m6A methylation relationships. Co-expression networks were constructed using WGCNA. R software and packages was employed to analysis and visualization throughout the entire process and statistical analysis. (C) showcases the overall comparison between Cohort 1 and Cohort 2 survival in both Groups. The impact of (D) EGFR, (E) KI-67, (F) MGMT, and (G) P53 markers on survival outcomes. Kaplan-Meier survival curves and log-rank test for statistical analysis for assessing significance (p < 0.05)

Fig. 2

Validating the purity of hierarchical cluster labeling (A) Workflow outlines the sequential steps involved in the comparison of label purity in clinical and histopathological groups. (B) Heatmap showcasing the comprehensive array of Diagnosis variable among the PD/psPD patients with glioblastoma. (C) Mosaic plot compares cluster formation across the Diagnosis variable within the studied PD/psPD patient cohort by Fisher’s Exact Test. (D) Workflow of the performance and validation of the HC. (E) Boxplot showing silhouette and ARI metric. (F) Barplot showing RMSE, ARMSE, MAPE and CVIM metrics comparing Cohort A as a training and Cohort B as test and the inverse. (G) Distribution plot of the three clusters, showing the two most relevant principal component

Fig. 3

Enrichment representation of top biological pathways. (A) workflow exploring the variance within methylated transcripts and delineate their impact on biological pathways. (B) density plot showcasing the distribution of variance among methylated transcripts. (C) A bar graph is presented illustrating the most highly enriched GO and KEGG pathways, ranked according to their p-values (p < 0.05). (D) A specific subnetwork is depicted, comprising of 23 critical genes, placing emphasis on the two primary GO pathways and the most prominent KEGG pathway((D(

Fig. 4

Analytical workflow employed in the WGCNA. (A) Workflow and Cluster Dendrogram derived from WGCNA, revealing the subdivision of genetic data into 23 distinct modules represented by unique colors. (B) Heatmap detailing the correlation patterns between these modules and various clinical and histopathological traits of the patients. (C) Boxplot compares Diagnostic PD versus psPD by each module colors. (D, H, and L) illustrate the differences in each Cohort A or Cohort B batch analyses of diagnostic data compared to modules colored in magenta, purple, and light cyan, respectively. (E, I, and M) highlight genes displaying differential methylation patterns within the magenta, purple, and light cyan modules, respectively. (F, J, and N) represent the enriched gene-gene associations and showcase the top 5 GO processes and KEGG pathways within augmented networks specific to the magenta, purple, and light cyan modules, respectively. All pathways present p-value < 0.05. Initial comparisons were screened by ANOVA and t test were used for specific groupings and significant effects are indicated by *p < 0.05; **p < 0.01; ***p < 0.001; **** p < 0.0001

Fig. 5

Comprehensive insight into differential methylation transcripts. Figure (A) showcases a heatmap depicting the distribution of 87 DMTs among patients. (B) the volcano plot delineates 38 hypomethylated and 49 hypermethylated DMTs. (C) Demonstrates a network diagram featuring key genes (depicted as green circles) with hyper- or hypomethylated status, their gene co-expressions (indicated by red lines), or protein-protein interactions (indicated by purple lines). These genes are associated with the top eight GO pathways (depicted as pink circles) and KEGG pathways (depicted as blue circles). All pathways present p-value < 0.05. Venn diagram showing the hub genes between DMT and (D) Magenta, (E) Purple, and (F) Lightcyan module

Fig. 6

Integrated analysis of gene expression and DNA methylation. (A) Venn diagram illustrates the overlap of differentially expressed genes with hyper- or hypomethylated status. (B) Depicts a heatmap displaying differentially methylation patterns of the 14 genes that are overlapping DMT and DEG. Matrix represents the m6A percentage 0 = 0% to 5 = 1200%. (C) Demonstrates a network diagram featuring key genes (depicted as green circles) with hyper- or hypomethylated status, their gene co-expressions (indicated by red lines), or protein-protein interactions (indicated by purple lines). These genes are associated with the top eight GO pathways (depicted as pink circles) and KEGG pathways (depicted as blue circles). Additionally, a node representing another gene (depicted as a green circle with a red border) not included in our primary analysis is shown, which typically exhibits co-expression with PPIF. All pathways present p-value < 0.05. Venn diagram showing the hub genes between DMT and (D) Magenta, (E) Purple, and (F) Lightcyan module