Integrative analysis of mRNA and miRNA sequencing data for gliomas of various grades

The purpose of this study was to characterize subtype-specific patterns of mRNA and miRNA expression of gliomas using The Cancer Genome Atlas (TCGA) data to search for genetic determinants that predict prognosis in terms of overall survival and to create interaction networks for grade 2 and 3 (G2 and G3) astrocytomas, oligodendrogliomas and grade 4 (G4) glioblastoma multiforme. Based on open-access TCGA data, 5 groups were formed: astrocytoma G2 (n = 58), astrocytoma G3 (n = 128), oligodendroglioma G2 (n = 102), oligodendroglioma G3 (n = 72) and glioblastoma G4 (n = 564); normal samples of brain tissue were also analysed (n = 15). Data of patient age, sex, survival and expression patterns of mRNA and miRNA were extracted for each sample. After stratification of the data into groups, a differential analysis of expression was carried out, genes and miRNAs that affect overall survival were identified and gene set enrichment analysis (GSEA) and interaction analysis were performed. A total of 939 samples of glial tumours were analysed, for which subtype-specific expression profiles of genes and miRNAs were identified and networks of mRNA-miRNA interactions were constructed. Genes whose aberrant expression level was associated with survival were determined, and pairwise correlations between differential gene expression (DEG) and differential miRNA expression (DE miRNA) were calculated. The developed panel of genes and miRNAs allowed us to differentiate glioma subtypes and evaluate prognosis in terms of the overall survival of patients. The regulatory miRNA-mRNA pairs unique to the five glioma subtypes identified in this study can stimulate the development of new therapeutic approaches based on subtype-specific mechanisms of oncogenesis.


Background
Gliomas are the most common type of tumours of the central nervous system and are characterized by a very poor prognosis, mainly due to the complexity of the treatment and the lack of targets suitable for developing effective therapeutic approaches. However, recent progress in molecular studies of gliomas is reflected in their classification and potential treatment goals. The main criteria for the modern classification of glial tumours (2016 edition of the World Health Organization classification) are the degree of malignancy (grades 1-4) and molecular markers [1,2]. However, the heterogeneity of gliomas is the main obstacle to their accurate molecular classification, and therefore, single markers are not enough to subtype and select the optimal treatment regimen.
With the development of next-generation sequencing technologies, researchers are now able to determine the most detailed molecular genetic profile for each sample. This influx of big data generated a need to create such large consortia as The Cancer Genome Atlas (TCGA), which allow accumulation of huge amounts of data for collaborative research of various types of neoplasia. Such studies have led to the construction of interactive mRNA-miRNA networks demonstrating regulatory mechanisms that are specific to certain nosologies [3][4][5][6]. The study of the expression profiles of mRNA and miRNA of various subtypes of gliomas and the construction of RNA-miRNA interaction networks can serve as a prerequisite for the development of highly accurate diagnostic tests and the development of personalized therapeutic approaches that will increase the life expectancy and improve the patient's physical condition.
The aims of our study were to research subtypespecific patterns of mRNA and miRNA expression of diffuse gliomas using TCGA data, search for genetic determinants that determine the prognosis in terms of overall survival and construct interaction networks for glial tumours of various grades.

Sample and data collection
Based on open-access TCGA data, 5 main groups of samples were formed: diffuse astrocytoma grade 2 (G2), anaplastic astrocytoma G3, diffuse oligodendroglioma G2, anaplastic oligodendroglioma G3 and glioblastoma (GBM), and these groups of samples were subjected to RNA and miRNA sequencing analysis (TCGA-GBM and TCGA-LGG datasets), as well as microarray analysis of miRNA data (for GBM) (Fig. 1). Normal brain tissues (n = 15) were used as comparison samples. Data collection and processing were performed using R 3.6.2 software ("TCGABiolinks" [7] and "GDCRNATools" packages [8]). The miRNA names were updated to miRBase version 21 using the "miRBaseConverter" package [9].
The clinical and demographic characteristics of the patients are shown in Table 1.

Analysis of differential gene and miRNA expression
For the differentially expressed gene (DEG) and differentially expressed microRNA (DE miRNA) analyses, the "limma" package [10] and the glmLRT test (generalized linear models likelihood ratio test) were used. The cutoff points were FDR < 0.01 and | logFC | > 1. The results are presented in the form of a scatter plot (volcano plot).

Gene set enrichment analysis
The genes selected in the analysis of differential expression were divided into groups depending on their genetic ontology (molecular functions, biological processes, cell localization and signalling pathways) using the "TCGAbiolinks" package [7] and the DAVID (Database for Annotation, Visualization and Integrated Discovery) database [11].

Survival analysis
The survival analysis involved several steps. First, genes were selected using the Kaplan-Mayer method in the TCGABiolinks package, the expression of which had a significant effect on the duration of overall survival. The expression of each gene with the quantile cut-off points (upper-0.6 and lower-0.33) was divided into low, moderate and high. The effect of the genes with high and low expression on survival was then evaluated using Cox multivariate regression analysis. Genes selected in the regression analysis were clustered based on the protein-protein interactions of their products (via the STRING database) [12]. We visualized the dependence of overall survival on the expression of the genes selected on the basis of the previous two analyses and interactive analysis using the "UCSCXenaTools" [13], "survival" and "survminer" packages.

Interactome analysis
Correlation analysis of gene and miRNA expression was performed for screening (via the "Hmisc" package). Gene-miRNA pairs with a Pearson correlation coefficient greater than − 0.5 were excluded. The p value threshold used for statistical significance was equal to 0.05. The remaining pairs were verified using ten databases ("multiMiR" package) [14]: miRecords, miRTarBase, TarBase, DIANA-microT-CDS, ElMMo, MicroCosm, miRanda, miRDB, PicTar and PITA. The final network of miRNA-target interactions included validated pairs (databases showed a status of "validated"). Predicted pairs (databases showed a status of "predicted") were selected based on the overall survival result. The interaction network was built using the "igraph" package.

Analysis of differential gene expression
A total of 939 glial tumour samples were analysed. Of the studied glioma subtypes, the DEGs (Table 2) were especially obvious for diffuse astrocytoma (Fig. 2).
When compared with non-tumour samples, diffuse astrocytoma samples showed 3392 DEGs, and anaplastic astrocytomas showed 3387 DEGs. The aberrant tumour expression profile was dominated by cases of increased transcriptional activity in astrocytoma G2 samples (1988 upregulated vs. 1404 downregulated) and, conversely, decreased expression in astrocytoma G3 samples (1400 upregulated vs. 1987 downregulated). In the GBM group, 2324 genes had differentially increased expression, and 2491 genes had differentially reduced expression. A slightly smaller number of differentially expressed genes was described for oligodendrogliomas (Table 2): 2723 DEGs in  When assessing the genetic ontology of the DEGs of astrocytomas, genes associated with EIF2 signalling, axonal guidance and neuropathic pain signalling in dorsal horn neurons were generally identified (Table 3). In addition, DEGs associated with G protein-coupled receptors and cAMP-mediated signalling were determined for diffuse astrocytomas, and role the DEGs of diffuse anaplastic astrocytomas were enriched in processes related to NFAT in cardiac hypertrophy and synaptic long-term potentiation.
Gene ontology analysis (Table 3) made it possible to identify common significant signalling pathways for oligodendrogliomas, such as the EIF2 pathway, mTOR signalling and dopamine-DARPP32 feedback in cAMP signalling, as well as subtype-specific pathways. The top signalling pathways of diffuse oligodendroglioma samples included axonal guidance signalling and neuropathic pain signalling in dorsal horn neurons, while the top pathways included G protein-coupled receptor signalling and synaptic long-term potentiation for anaplastic oligodendroglioma samples. It is worth noting that expression related to synaptic long-term potentiation signalling was aberrant in both astrocytoma and oligodendroglioma G3 samples, suggesting that it may be related to advanced stage ( Table 3).
The top five signalling pathways for GBM samples were axonal guidance signalling, Rho family GTPases pathway, breast cancer regulation by stathmin1, molecular mechanisms of cancer and IL-8 signalling. In this case, the last four of the listed pathways were able to differentiate GBM into a separate group. Figure 3 illustrates the distribution of genes that are differentially expressed either in all the tumours analysed, in some subtypes or a specific subtype. As a result, 1036 upregulated and 517 downregulated universal DEGs were revealed in samples of glial tumours. Of the 1553 nonspecific DEGs, the most represented were participants in the following signalling pathways: calcium, MAPK, oxytocin and neuroactive ligand receptor pathways. These signalling pathways can mediate neural network activity and cell proliferation. A more pronounced overlap of aberrant DEGs was observed between the GBM, astrocytoma G2 and astrocytoma G3 groups, which can be explained by the astrocytic origin of GBM.
Gene ontology analysis of the 210 unique DEGs in diffuse astrocytomas, shown as an example in Fig. 4,    identified changes in transcriptional activity of genes encoding passive transmembrane transport proteins, which are located mainly on the plasma membrane of cells (Fig. 4).
The 119 DEGs specific to astrocytoma G2 samples were identified to be involved in the biological processes of translation and translational elongation. Most of the DEGs perform the molecular function of transmembrane transport or are structural components of ribosomes. Sixty-three DEGs specific to astrocytoma G3 samples were found to generally be involved in translational elongation and ion transport, as were the 147 DEGs specific to oligodendroglioma G2 samples (Online Resource 6). GBM was obviously different from the other subtypes, demonstrating 1716 unique DEGs, some of which were found to be involved in biological processes such as sensory perception of chemical stimuli, cognitive and systemic neurological processes and performing transmembrane transport functions.

Analysis of differential miRNA expression
The results of evaluating the differential expression of miRNAs in different types of glial tumours are shown in Table 2. As an example, the results for diffuse astrocytomas are illustrated in Fig. 5. Of note, the number of identified DE miRNAs was an order of magnitude less than the number of identified DEGs. Moreover, there were~2 times more upregulated DE miRNAs than downregulated DE miRNAs for astrocytomas and oligodendrogliomas. The number of GBM-specific miRNAs was an order of magnitude less than the number of GBM-specific DEGs, and their distribution according to the degree of altered activity was approximately equal ( Table 2).
Venn diagrams (Fig. 6) were constructed to show the differential expression of miRNAs depending on the glioma subtype.  As a result, 35 microRNAs were identified, and the degree of change for these microRNAs was similar for all types of glial tumours, including GBM. It should be noted that 207 common DE miRNAs were identified in astrocytoma and oligodendroglioma G2 and G3 samples, while GBM samples had only 69 DE miRNAs common to other subtype samples. The largest number of specific small noncoding RNAs was determined for G3 tumours: 65 (vs. 21) miRNAs discriminated astrocytoma G3 samples; and 51 (vs. 10) miRNAs discriminated oligodendroglioma G3 samples.
Association of overall survival with the expression of markers differentiating the studied subtypes of glial tumours Based on the initial screening data, a pool of genes was formed, the differential expression of which not only was unique for each subtype but also correlated with the overall survival of patients (Online Resource 20). The resulting panel included 26 DEGs unique to astrocytoma G2; 15 DEGs unique to astrocytoma G3; 17 DEGs unique to oligodendroglioma G2; 7 DEGs unique to oligodendroglioma G3; and 151 DEGs unique to GBM. Table 4 summarizes the unique survival-related DEGs that were not previously mentioned in glial tumour studies in PubMed and Cochrane.

Interactome analysis of miRNA-mRNA interactions
For each tumour subtype, an analysis of DEG associations with overall survival was performed. A total of 360 and 835 DEGs were identified in astrocytoma G2 and G3 samples, respectively, and 227 and 468 DEGs were identified in oligodendroglioma G2 and G3 samples, respectively. Genes that correlated with survival were clustered depending on the significance level of PPIs (protein-protein interactions according to STRING analysis). For example, 6 clusters were revealed for astrocytoma G2 samples (Fig. 7), 10 clusters were revealed for astrocytoma G3 samples and 1, 8 and 8 clusters were revealed for oligodendroglioma G2, oligodendroglioma G3 and GBM samples, respectively (Online Resource 13,14,15,16).
Thereafter, Pearson correlation analysis was performed on the expression levels of selected genes and miRNAs, and pairs with correlation values of − 0.5 or lower were considered significant. Figure 8 presents an example of the results of the interactome analysis for diffuse astrocytoma, and mRNA-miRNA pairs, which were verified against databases of validated/predicted miRNA targets, are visualized.
Thus, 173 DE miRNA and DEG mRNA pairs were associated with survival in the samples of diffuse astrocytoma, 76 DE miRNA and DEG mRNA pairs were associated with survival in diffuse oligodendroglioma samples and 199 and 372 DE miRNA and DEG mRNA pairs were associated with survival in anaplastic astrocytoma and anaplastic oligodendroglioma, respectively. For GBM, no correlation was found between DEGs and DE miRNAs.

Discussion
In our study, for each of the five subtypes of glioma, specific gene and microRNA expression profiles were described ( Table 2). As a result, patterns of specifically expressed markers were determined: 211 unique DEGs and 21 DE miRNAs were identified in astrocytoma G2, 63 and 65 were identified in astrocytoma G3, 119 and 10 were identified in oligodendroglioma G2 and 147 and 51 were identified in oligodendroglioma G3.
The identified expression patterns were generally related to signalling pathways associated with neuronal activity, proliferation and cancer progression. In glial tumour samples at earlier stages (G2 and G3 samples), a similarity in the set of signalling pathways with aberrant activity was demonstrated (Table 3). At the same time, signalling pathways that were not shared by the astrocytomas or oligodendrogliomas were identified for the GBM group.
Based on the initial screening data, a panel of genes was developed, the differential expression of which not only was unique for each subtype but also correlated with the overall survival of patients (Online Resource 20). The panel includes 26 DEGs unique to astrocytoma G2, 15 DEGs unique to astrocytoma G3, 17 DEGs unique to oligodendroglioma G2, 7 DEGs unique to oligodendroglioma G3 and 151 DEGs unique to GBM.
Currently, the role of individual miRNAs in glial tumours is not well understood. MiRNAs are small noncoding regulatory RNAs that decrease the stability and/ or inhibit the translation of target mRNAs with which they have full or partial complementarity. These regulatory molecules, playing an important role in the pathogenesis, development and progression of cancer and the response to therapy, are especially promising for diagnostic strategies, including minimally invasive methods (measuring the level in plasma, serum or cerebrospinal fluid).
Primary data analysis via TCGA allowed us to differentiate diffuse astrocytomas and oligodendrogliomas from other subtypes of gliomas by the expression of 21 and 10 miRNAs, respectively; anaplastic astrocytoma   and oligodendroglioma could be differentiated from other gliomas based on the expression of 65 and 51 miRNAs, respectively. For the GBM group, 13 unique DE miRNAs were identified. After the Pearson correlation analysis and the validation of the correlation pairs in the databases, unique miRNAs that were found to target the DEGs were identified for some subtypes, and these miRNAs were associated with survival. Below, we briefly consider their characteristics.
In diffuse astrocytomas, a decrease in the levels of hsalet-7c-5p and hsa-miR-125b-5p was detected. In parallel, it was noted that the expression of hsa-let-7c-5p was associated with the level of transcriptional activity of the BEGAIN gene; the link between hsa-miR-125b-5p and  the mRNA level of the DUSP3 and RAB3D genes was similarly found. Earlier, in an experimental study on cells of acute erythroleukaemia, the antiproliferative effect of let-7c-5p was demonstrated to be associated with the targeting of the PBX2 oncogene [24]. In colorectal cancer cells after exposure to 5-FU or starvation, Hou N et al. identified a miRNA pool including hsa-let-7c-5p [25]. The identified correlation between the expression of hsa-let-7c-5p and BEGAIN is consistent with the presence of one binding site (target score-90) within the 3′-UTR of the mRNA of this gene. The BEGAIN protein is involved in the regulation of the activity of postsynaptic neurotransmitter receptors (GO: 0098962); the contribution of this locus to oncogenesis has not been confirmed to date. The antitumour effect of miR-125b has been described in several studies of lung cancer, hepatocellular carcinoma and breast cancer [26][27][28][29]. For example, the study of Yuan M et al. demonstrated the antiproliferative effect of increased expression of hsa-miR-125b-5p in cell lines of diffuse astrocytoma in children [30]. In contrast, a number of studies of endometrial cancer, prostate cancer and GBM indicate the oncogenic effects of hsa-miR-125b-5p, which targets the known tumour suppressors TP53 and GJA1 [31][32][33]. The target of hsa-miR-125b-5p in diffuse astrocytomas is probably the mRNA of the DUSP3 gene, a member of a subfamily of protein phosphatases that negatively regulate the members of the mitogenactivated protein kinase superfamily (MAPK/ERK, SAPK/JNK and p38) and in such a way control proliferation and differentiation. The RAB3D protein (which has three 3′-UTR binding sites for miR-125b-5p) is described as an oncogenic small GTPase that promotes both cell migration in glial tumours and metastasis in breast cancer via the intracellular AKT/GSK3β signalling pathway [34,35]. In a single tumour, hsa-miR-125b-5p exhibits both oncosuppressive and oncogenic effects.
Anaplastic astrocytomas were characterized by a change in the level of three unique miRNAs: hsa-miR-190b-3p, hsa-miR-196a-3p and hsa-miR-218-2-3p. Correlation analysis allowed us to associate the level of hsa-miR-190b-3p with the expression of the EFNB1 gene, which is involved in cell adhesion and the formation of neural synapses. In this case, the expression of the hsa-miR-190b-3p target was increased, which can be explained by the presence of additional mechanisms regulating the transcription of EFBN1. However, hsa-miR-190b-3p may mitigate the increased expression of EFBN1. An increase in the level of hsa-miR-190b-3p has also been shown in other types of malignant neoplasms. For example, validation of an analysis of 1083 breast cancer samples (from TCGA) confirmed an increased level of expression of hsa-miR-190b-3p in both surgical biopsies and breast cancer cell lines, and this overexpression significantly contributed to cell proliferation and migration [36]. The expression of hsa-miR-190b-3p correlates with the expression of the C8orf34 gene, which is involved in the regulation of gene expression and the cell cycle. Mutations at this genetic locus can mediate severe toxicity during chemotherapy [37]. The increased activity of hsa-miR-196a-3p in colorectal cancer was demonstrated by the study of Wang YN [38]. Similarly, Chen Y. et al. showed that a decrease in hsa-miR-196a-3p expression was associated with a high risk of metastasis in breast cancer. The authors found that transforming growth factor β1 negatively regulates the expression of hsa-miR-196a-3p and activates neuropilin-2, which causes the migration and invasion of tumour cells [39].
Our study demonstrated a decrease in hsa-miR-218-2-3p expression. As a result of the correlation analysis and subsequent database validation, the correlation of hsa-miR-218-2-3p and LMO7 expression was revealed. In contrast, earlier in the Feng Z study, overexpression of hsa-miR-218-2-3p was detected both in cell lines and in surgical biopsy specimens of glial tumours of various grades. In the study, an increase in the activity of has-miR-218-2-3p was associated with the invasive and proliferative potential of glial tumours [40]. As a possible target of hsa-218-2-3p, LMO7 is described as a transcription factor and stabilizer of adhesion compounds, a change in the expression of which is associated with some oncological nosologies [41,42]. Decreased LMO7 expression in mouse models leads to the development of lung adenocarcinoma [43], and low expression in human lung adenocarcinoma was associated with a negative prognosis [42]. However, contrasting results have also been shown, in which high expression of LMO7 was identified as a negative prognostic factor in non-smallcell lung cancer associated with LRIG1 expression [44].
Despite the numerous miRNA-mRNA correlation pairs (n = 75) identified in diffuse oligodendrogliomas, there were no subtype-specific miRNAs among them.
Earlier, in a study of colorectal cancer (CRC), the association of hsa-miR-20b-5p with survival, as well as a negative correlation of this miRNA with one of the E2F5 cell cycle regulators, was noted [45]. In addition, increased levels of this microRNA in the blood were associated with a better prognosis in patients with metastatic CRC receiving chemotherapy based on bevacizumab [46]. Our analysis demonstrated a negative correlation between hsa-miR-20b-5p expression and PRKACB. This genetic locus has been identified as a significant oncogene involved in the progression of endocrine cancer by modulating the signalling activity of cAMP [47,48]. In turn, non-small-cell lung cancer shows a lower expression level of PRKACB than normal tissues, and an increase in PRKACB transcriptional activity reduces the number of proliferative, colony-forming and invasive cells and increases the incidence of apoptosis [49].
In a study devoted to the development of predictors of GBM treatment, high expression of hsa-let-7a-5p was interpreted as a favourable prognostic factor associated with survival [50]. In breast cancer, hsa-let-7a-5p mediates sensitivity to the proteasome inhibitor bortezomib, regardless of the molecular subtype [51]. The potential significance of let-7a-5p was determined by Zhou X and colleagues in a study of non-small-cell lung cancer [52].
Using layer-by-layer data filtering, the hsa-let-7a-5p and CYP46A1 correlation pair was identified. The CYP46A1 genetic locus encodes cholesterol 24hydroxylase, an enzyme responsible for eliminating cholesterol in 24S-hydroxycholesterol (24S-OHC, oxysterol) that is specific to the brain [53]. A decrease in CYP46A1 expression was demonstrated in a study by Han M et al., in which in silico analysis of DEGs and subsequent in vitro verification were performed for GBM samples. Single-cell transcriptome sequencing data demonstrated prevailing CYP46A1 expression in neurons, astrocytes and oligodendrocyte progenitor cells compared to tumour cells [54]. Overexpressed CYP46A1 modulated a decrease in colony formation, proliferation and generation of tumour spheroids in vitro due to reduced cholesterol accumulation [54].
The level of hsa-miR-466 was found to be decreased in osteosarcoma and correlated with IRS1 expression. Increased expression of has-miR-466 inhibited the in vitro proliferation and invasion of osteosarcoma cells [55]. Similarly, a decrease in hsa-miR-466 expression has been demonstrated in hepatocellular carcinoma cells, and an increase in hsa-miR-466 expression inhibits proliferation, induces apoptosis and reduces the metastatic potential of the tumour by targeting MTDH [56]. In prostate cancer, a decrease in hsa-miR-466 expression was demonstrated due to the activity of the long noncoding RNA TUC338. The authors also revealed the antiproliferative potential of hsa-miR-466 [57]. However, an increase in hsa-miR-466 transcriptional activity was observed in malignant tumours of the cervix uteri [58].
In our study, increased expression of hsa-miR-466 negatively correlated with the mRNA level of the CKS2 oncogene. CKS2 encodes a regulatory subunit of 2 cyclin-dependent kinases. Increased expression of CKS2 is associated with tumour progression in bladder cancer, hepatocellular carcinoma, breast cancer and GBM [59][60][61][62].
We did not find microRNA-mRNA correlation pairs in the GBM group. However, thirteen expressed miR-NAs were specific for GBM. These miRNAs are also very promising for screening panels. It is worth noting that in the TCGA GBM project, there is an extreme lack of miRNA sequencing data. Future studies using next-generation sequencing technology will help uncover new links and mechanisms in GBM oncogenesis.
Previously, studies have been carried out aimed at finding genetic markers for the molecular classification of glial tumours. In the study of Shai R et al. [63], cluster analysis of 170 DEGs made it possible to differentiate the subtypes of astrocytomas, oligodendrogliomas and glioblastomas. However, the number of samples included in the study was relatively small-35 samples of gliomas of various grades of malignancy. In addition, the percentage of errors in predicting survival was relatively high (22%) [63]. Later, a transcriptomic study of 225 samples of glial tumours was carried out. The authors were able to find a relationship between a favourable prognosis and a mutation in IDH1/2, which was most common in G1 tumours. Nonetheless, no subtypespecific genetic signature was found [64]. There are also a number of other research papers devoted to the study of the molecular characteristics of one or two types of glial tumours [65][66][67]. Our study is characterized by a large number of samples and tumour subtypes (939 samples of 5 different subtypes of gliomas were analysed) and an integrative approach, which consists in the analysis of signalling pathways, survival and the search for subtype-specific markers among DEGs and DE miRNAs.
The results obtained open up broad prospects in the early diagnosis of certain subtypes of gliomas and in predicting the overall survival of patients using miRNA and mRNA panels. We hope that further experimental studies of expressions subtype-specific miRNAs and mRNAs will confirm the effectiveness of appropriate genetic tests, which will allow clinicians to timely stratify patients depending on the tumour subtype and choose the optimal treatment regimen.

Conclusions
As a result of our study, the molecular genetic characteristics of five glioma subtypes were described: astrocytoma G2, astrocytoma G3, oligodendroglioma G2, oligodendroglioma G3 and GBM G4. We have shown that each subtype has a unique expression profile of genes and miRNAs and has differences in key signalling pathways and regulatory networks of miRNA-RNA interactions. Based on the analysis of differential expression and overall survival, we formed a panel of genes and miRNAs that allows the identification of different glioma subtypes and the evaluation of prognosis in terms of overall survival. The identified regulatory miRNA-mRNA pairs unique to the five glioma subtypes can stimulate the development of new therapeutic approaches based on subtype-specific regulatory mechanisms.