This article is mentioned in:

Introduction

As a biologically aggressive subtype of malignant glioma, glioblastoma multiforme (GBM) is the most common and lethal brain tumour in adults (1). GBM diffusely infiltrates the brain at an early time point, and its complexity and distinct pathophysiology are facilitated and dictated by the unique brain milieu and cellular interactions (2,3). In total, >80% of GBM cases are primary tumours, which typically respond poorly to current therapeutic approaches (4). Although previous studies detected a number of events associated with initiation or progression (5,6), identification of the key molecular targets underlying the regulatory mechanisms of GBM remains to be elucidated.

At present, the wealth of molecular information generated by massively parallel sequencing technologies provides a great opportunity for discovering novel biomarkers in different types of cancer (7,8). Transcriptome analysis involving protein-coding genes and noncoding microRNAs (miRNAs/miRs) has proven to be a valuable first step for studying the genetic characteristics of various complex human diseases, especially in malignant gliomas (9,10). For example, Lai et al (11) found that miRNA-210 could serve as a potentially non-invasive biomarker for the diagnosis and prognosis of GBM. An increased expression level of the protein-coding gene PLAUR was determined to be a valuable predictor of the mesenchymal GBM subtype (12). However, it is unlikely that a single dysfunctional gene would lead to abnormal phenotypes. Multiple aberrant gene-gene interactions may have co-operative effects in the development and progression of complex diseases such as cancer (13,14).

Understanding the molecular events associated with activation of the risk pathways in patients with GBM may facilitate the development and clinical testing of potential therapeutic targets.

In the present study, differentially expressed genes and miRNAs were identified by analysing a large protein-coding gene and miRNA expression dataset comprised of normal and tumour tissues from GBM samples obtained from The Cancer Genome Atlas (TCGA). A GBM-specific regulatory network involving transcription factors (TFs) and miRNAs was built. The GBM-specific regulatory network showed a scale-free network with a small-world property. The key molecules in the regulatory network were significantly associated with immune-related functions. Linear risk pathways of GBM were identified through systems-level analyses of the regulatory network. The identified molecular risk pathways may be a potential resource for understanding the pathogenesis and aetiology underlying GBM.

Materials and methods

Acquiring known GBM-related genes

The GeneCards database provides a comprehensive map of manually curated human disease genes (15). In total, 32 GBM-associated genes were obtained by searching the GeneCards database according to all descriptions and aliases of GBM, including ‘GBM’, ‘glioblastoma multiforme’, ‘malignant brain tumour’, ‘gliomatosis’ and ‘high-grade glioma’. In addition, 41 experimentally validated GBM-related miRNAs from the miR2Disease database were retrieved (16), a manually curated database of miRNA deregulation in various human diseases.

Identification of differentially expressed genes and miRNAs

The mRNA and miRNA expression profiles for 378 patients with GBM were downloaded from TCGA database (https://portal.gdc.cancer.gov/). The dataset contained 378 tumour samples and 10 patient-matched normal samples. The genes were mapped to the Entrez Gene IDs for mRNA expression data using biomaRt software version 2.40.0 (17). In total, 11,273 genes and 470 miRNAs were included in the subsequent analysis. A Student's t-test was used to identify significantly differentially expressed genes and miRNAs from the GBM expression profiles. All P-values were adjusted using Bonferroni correction to account for multiple comparisons. Only genes and miRNAs with an adjusted P-value <0.01 and a fold change >1.2 were considered statistically significant.

Extraction of the GBM-specific regulatory network

A TF and miRNA regulatory network was constructed by integrating the molecular interaction associations between genes and miRNAs from four databases [miRTarBase (18); The Transcription Factor Database (19); miRecords (20); and TarBase (21)]. To ensure the reliability of the relationships between the curated regulations, each individual interaction presented is supported by the literature.

Although differential expression analysis is often used to screen disease candidate genes, previous studies have shown that a number of disease-associated genes do not always exhibit differential expression patterns in microarray experiments (14,22). Therefore, a candidate risk regulatory network was constructed by connecting all of the risk-associated differentially expressed nodes with their immediate neighbouring non-differentially expressed nodes.

Mining for GBM-specific risk regulatory linear pathways

Linear risk pathways potentially associated with GBM were obtained according to the following criteria: i) The beginning of a risk pathway was defined as a gene or miRNA with a 0-indegree (it could not be regulated by other genes or miRNAs); and ii) the end of a risk pathway was defined as a gene or miRNA with a 0-outdegree (it could not regulate other genes or miRNAs). In the GBM-specific regulatory network, all linear regulatory pathways were defined as pathways from the beginning node to the end node and were detected using R (23).

Functional enrichment analysis

The R package ‘clusterProfiler’ (24) was used to examine the significantly enriched gene ontology (GO) biological functions (25) of the differentially expressed genes. The adjusted P-values were calculated using the Benjamini and Hochberg method.

Results

Differential expression analysis of mRNAs and miRNAs

To detect potential key molecules underlying the pathology of GBM, the mRNA and miRNA gene profiles of patients with GBM from TCGA database were downloaded and differential expression analysis was performed. In total, 1,827 and 30 significantly differentially expressed genes and miRNAs were detected, respectively. The number of upregulated genes, (1,085) was greater than the number of downregulated genes (742). Conversely, the number of upregulated miRNAs (11) was smaller than the number of downregulated miRNAs, 19 (Fig. 1). These results suggested that a number of candidate molecules may contribute to the carcinogenesis of GBM, and that miRNAs may serve a negative regulatory role.

Figure 1. Volcano plot of the differentially expressed genes and miRNAs associated with glioblastoma multiforme. A total of 1,827 and 30 significantly differentially expressed genes and miRNAs, respectively were detected. Red dots represent genes and miRNAs with ≥1.2-fold change increase and blue dots represent genes and miRNAs with ≤1.2-fold decrease. miRNA, microRNA.

Systematic characterization of a regulatory network specific to GBM

To examine the systems-level characteristics of the ‘true’ biological processes, including interactions and other combined activities associated with GBM, a comprehensive regulatory network among TFs, miRNAs and genes was built by combining four manually curated databases that focused on both transcriptional and post-transcriptional regulation. The background regulatory network was composed of 6,036 regulatory relationships. Of the background regulatory relationships, 1,827 differentially expressed genes and 30 differentially expressed miRNAs were considered as the risk seeds and mapped into the regulatory network to obtain a GBM-specific sub-network (Fig. 2). Given that networks provide knowledge that can be used to infer the information flow between regulators and target genes, a GBM-specific sub-network consisting of differentially expressed genes/miRNAs and their immediate neighbouring nodes from the TF-miRNA-gene regulatory network was further extracted. The GBM-specific regulatory network contained 160 interactions involving 23 TFs, 17 miRNAs and 105 target genes.

Figure 2. Glioblastoma multiforme-specific regulatory network. Green square nodes represent transcription factors, orange pointed nodes represent microRNAs and blue round nodes represent target genes.

The network showed dense local interconnectivity and complied with a power-law distribution suggesting that the GBM-specific sub-network was scale-free, similar to most biological networks (26) (Fig. 3). The biological functions of the GBM-specific network were investigated. Functional enrichment analysis was performed, and the top 30 significantly enriched Gene Ontology functions are presented in Fig. 3. Notably, the results demonstrated a number of immune response-associated functions, including ‘neutrophil-mediated immunity’, ‘neutrophil activation involved in immune response’ and ‘leukocyte differentiation’. Numerous studies have reported profound and generalised immunosuppression of GBM tumours, particularly within the context of cell-mediated immunity (27,28). In addition, a number of other types of functions were found, such as ‘axon development’, ‘haemostasis’ and ‘response to peptide’, which might be unique or important to the brain. Together, these results support the possibility that the GBM-specific sub-network may be a potential resource for GBM-related research.

Figure 3. Characterization of the GBM-specific regulatory network. Top 30 significantly enriched biological processes in the GBM-specific regulatory network. The colour bar on the right represents the adjusted P-value and the circle size represents the number of genes.

GBM-specific risk pathways

A biological interaction network contains a large amount of information, and its complex structure may hinder the deduction of new insights for interested risk pathways. The risk pathways obtained from the regulatory network can shed light on the molecular mechanisms of GBM. Detecting disease-associated pathways in a clear and simple manner is an important goal. In the present study, 615 linear risk pathways from the GBM-specific regulatory network were identified. Notably, a TF named EGFR was not only a known GBM disease gene, but also a differentially expressed gene (29). A total of 15-risk pathways were connected by EGFR and its up-stream regulator hsa-miR-7 (Fig. 4), which is a miRNA that has been reported as a prospective risk factor for GBM (30,31).

Figure 4. Risk regulatory pathways connected by EGFR and hsa-miR-7. EGFR, epidermal growth factor; miR, microRNA.

Discussion

Glioblastoma is an astrocyte-derived tumour with a propensity for malignancy and is caused by an abnormal disorder of multiple interacting genes rather than a single gene acting in isolation (32). Activation of molecular signalling cascades have been shown in GBM tumourigenesis (33). Identification of a cascade regulatory pathway, which usually has a clear and simple linear structure (34), may help reveal molecular mechanisms of the disease. A network-based approach to analyse risk regulatory pathways in GBM was used in the present study. Noncoding RNAs are emerging as key regulators of diverse biological processes and are important in cancer pathogenesis (35,36). Therefore, noncoding RNAs involved in pathway regulation were also considered.

Through differential expression analysis, a total of 1,827 genes and 30 miRNAs were identified in normal and tumour tissues from GBM samples. When these candidate risk factors were taken as seed nodes, a GBM-specific regulatory network involving TFs, miRNAs and target genes was identified. Functional analysis showed that the GBM-specific network was significantly associated with several important functions of the immune response, axon development, haemostasis and response to peptides. From these, 15 reliable risk pathways specific to GBM were identified, which may be a useful resource for increasing our fundamental understanding of the molecular role and regulation of GBM. These results showed that the 15 GBM-associated risk pathways were connected by a TF named EGFR that was both a known GBM gene and a differentially expressed gene. Numerous studies have highlighted the vital role played by EGFR during the emergence and development of GBM. Shao et al (37) showed that ~80% of patients with GBM presented with EGFR genomic amplification or overexpression. Mischel et al (38) reported that EGFR-overexpressing GBMs promoted GBM cell proliferation, survival and angiogenesis. Wheeler et al (39) showed that EGFR inhibitors may be a novel therapeutic strategy for treating patients with GBM. In the present study, EGFR indirectly regulates 15 GBM-related genes (Fig. 4) and forms a cascade regulatory pathway with hsa-miR-7 and hsa-miR-21, participating in the occurrence and development of GBM. Regarding the downstream target genes, FAS-induced expression of chemokines in human glioma cells (40) and TGFβ1 gene expression are potential signatures for the mesenchymal high-grade glioma subtype (41). Therefore, the risk regulatory pathways miR-7, EGFR, miR-21 and FAS, and miR-7, EGFR, miR-21 and TGFB1 may cause disease by affecting the abnormal function of glioma cells.

In conclusion, accumulating evidence indicates that a single mutated or aberrantly expressed gene is not sufficient to result in an abnormal phenotype and that complex diseases are caused by the accumulative and co-operative effects of multiple dysregulated genes and gene interactions (42,43). Risk regulatory pathways that consist of interacting genes with a linear chain structure can help reveal the molecular mechanisms of diseases. Through integration of transcriptional and post-transcriptional regulation, risk regulatory pathways in GBM that might serve as potential targets for clinical treatment were identified.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the author on reasonable request.

Authors' contributions

JL, YJX, YW, and YY collected and analyzed all data. CZ, YX and JXW performed statistical analysis. ZYL, JL and JXW contributed to the concept and design of the study and drafted the manuscript. JL, YW, and YY revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References