BioMed Research International
Volume 2019 (2019), Article ID 6581576, 13 pages
https://doi.org/10.1155/2019/6581576

Research Article

Identification of Potential Biomarkers in Glioblastoma through Bioinformatic Analysis and Evaluating Their Prognostic Value

Yangmei Zhou, Li Yang, Xiaoxi Zhang, Rui Chen, Xiuqiong Chen, Wenhua Tang, and Mengxian Zhang

Department of Oncology, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China

Correspondence should be addressed to Mengxian Zhang; zhangmx73@163.com

Received 8 January 2019; Revised 3 March 2019; Accepted 11 March 2019; Published 15 April 2019

Academic Editor: Steven De Vleeschouwer

Copyright © 2019 Yangmei Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Glioblastoma is a common malignant tumor in the central nervous system with an extremely poor outcome; understanding the mechanisms of glioblastoma at the molecular level is essential for clinical treatment. In the present study, we used bioinformatics analysis to identify potential biomarkers associated with prognosis in glioblastoma and elucidate the underlying mechanisms. The result revealed that 552 common genes were differentially expressed between glioblastoma and normal tissues based on TCGA, GSE4290, and GSE 50161 datasets. Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and protein-protein interaction (PPI) network were carried out to gain insight into the actions of differentially expressed genes (DEGs). As a result, 20 genes (CALB1, CDC20, CDCA8, CDK1, CEP55, DLGAP5, KIF20A, KIF4A, NDC80, PBK, RRM2, SYN1, SYP, SYT1, TPX2, TTK, VEGFA, BDNF, GNG3, and TOP2A) were found as hub genes via CytoHubba in Cytoscape and functioned mainly by participating in cell cycle and p53 signaling pathway; among them, RRM2 and CEP55 were considered to have relationship with the prognosis of glioblastoma, especially RRM2. High expression of RRM2 was consistent with shorter overall survival time. In conclusion, our study displayed the bioinformatic analysis methods in screening potential oncogenes in glioblastoma and underlying mechanisms. What is more is that we successfully identified RRM2 as a novel biomarker linked with prognosis, which might be expected to be a promising target for the therapy of glioblastoma.

1. Introduction

Glioblastoma (GBM) remains one of the most common aggressive tumors in the central nervous system [1]. Though the incidence of GBM is about 3.19/100,000 population [2], the prognosis is extremely poor because of the complex biological characteristic and limited treatment options. After standard therapy, maximal safe surgical resection is followed by combined radiochemotherapy; still approximately 70% of patients die within two years, while over 90% die within five years [3–6]. Thus, it is critical to comprehend the mechanisms of GBM and develop some more effective therapeutic strategies to improve the outcome of patients.

Over the past years, tremendous studies have been conducted to explore the potential molecular mechanisms, genetics, and pathways responsible for the development and progression of GBM. However, the precise mechanism of GBM remains unclear. Recently, biological information analysis has caused extensive attention and made consistent breakthrough in searching the oncogenes; various biomarkers for diagnosis and prognosis of cancer have been identified [7–9], so it could also be expanded in GBM for better understanding of the underlying molecular mechanism in GBM and finding a clue for new therapeutic targets.

In this study, we performed biological information analysis among 3 profiles (TCGA, GSE4290, and GSE50161), which were downloaded from the TCGA database and the GEO database and identified the differentially expressed genes (DEGs) between the tumor and normal samples. Besides, functional enrichments and the protein-protein interaction (PPI) networks were applied to annotate gene function and screen hub genes. Prognostic value was evaluated by survival analysis. Findings of our study may hint at a potential prognostic biomarker and therapeutic target.

2. Materials and Methods

2.1. Microarray Data

The mRNA expression profiles and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/), which contains 169 GBM samples and 5 normal samples; the other two profiles, GSE4290 and GSE50161, were downloaded from the Gene Expression Omnibus (GEO) database (http://www.earthobservations.org/index2.php); we extracted 81 GBM samples and 23 normal samples from GSE4290, and 34 GBM samples and 13 normal samples from GSE50161. Both of the GEO data were generated using the GPL570 (HG-U133_PLUS_2) Affymetrix Human Genome U133 Plus 2.0 microarray platform data.

2.2. Data Preprocessing and DEGs Screening

We utilized the R statistical software (version 3.4.4; https://www.r-project.org/) and Bioconductor analysis tools (http://www.bioconductor.org/) to process the raw data. The edgeR package in R was used to screen DEGs between tumor and normal tissues with the threshold of false discovery rate (FDR)<0.01, >2 in TCGA datasets; the corresponding survival information of 159 glioblastoma samples was obtained after deleting the incomplete statistics. For the GEO datasets, we converted the probe-level data into gene expression profiles based on their platform annotation files; the gene expression values were averaged when analogous to multiple probe; DEGs only with the >2 and FDR<0.05 were extracted with the limma package after background correction and data normalization. The differences between the two groups were calculated by Student’s t test; the P-value was adjusted by Benjamini-Hochberg (BH) method.

Venn analysis was utilized to choose the overlap of the DEGs among the three datasets above. As a result, 552 overlapped DEGs in GBM samples were set as a cohort for further exploring.

2.3. Function and Pathway Enrichment Analysis

The GO analysis was undertaken from the following three categories: molecular function (MF), biological process (BP), and cellular component (CC). KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was also performed to search the key biological pathways the candidate genes were involved in. The functional analysis was conducted through The Database for Annotation, Visualization and Integrated Discovery (DAVID: https://david-d.ncifcrf.gov/) online tools [10] and separately organized by upregulated and downregulated groups. The p-value was set at <0.05, which means an enrichment score and represented the significance of the GO or pathways terms.

2.4. PPI Network Construction and Hub Genes Identification

The Search Tool for the Retrieval of Interacting Genes [11] (STRING: https://string-db.org/) is an online tool to show the known and predicted interaction relationship between the candidate genes; we used it to develop the overlapped DEGs’ s interaction network. Cytoscape software (version 3.6.1, http://www.cytoscape.org/) was applied to visualize the network of these gene interaction pairs with coefficients of and P<0.05, and the hub genes in the network were identified by utilizing the cytoHubba application in it.

2.5. Survival Analysis

We performed survival analysis through R language survival package; Kaplan-Meier survival curves were plotted after adjusting relative hazards from the Cox proportion hazards model; P<0.05 was considered to have statistical significance.

2.6. Correlation Analysis

The correlation between the expression of RRM2 and other hub genes was evaluated through GEPIA based on TCGA database. Pearson’s correlation coefficient analysis was used to define correlations.

3. Results

3.1. Identification of DEGs

We filtered a total of 2932 DEGs including 1341 upregulated and 1591 downregulated genes between the tumors and adjacent normal samples from TCGA database, and a total of 2039 DEGs (816 upregulated and 1223 downregulated) and 823 DEGs (212 upregulated and 611 downregulated) from the expression profile datasets GSE50161 and GSE4290 (Figure 1). Finally, 552 common DEGs from three datasets were determined via Venn analysis (Figure 2), including 128 upregulated genes and 424 downregulated genes, which were subsequently used for further investigation.

(a)

(b)

(c)

Figure 1: Volcano plot of microarray. Horizontal axis represents log2-fold change, and vertical axis represents adjusted p value. (a) GSE4290. (b) GSE50161. (c) TCGA.

Figure 2: Venn analysis of DEGs. A comparison of 2039, 2932, and 823 DEGs from three datasets revealed 552 common DEGs between GBM and normal tissues.

3.2. GO and KEGG Pathway Analysis for Overlapped DEGs

Those common DEGs were well enriched in abundant of functional groups and pathways; the top 5 terms of BP, CC, MF, and KEGG pathway were shown in Tables 1 and 2. In upregulated group, the most enriched terms in each category were GO:0030198 extracellular matrix organization (BP, P-value 3.88E-12), GO:0031012~extracellular matrix (CC, P-value 6.61E-13), GO:0005201~extracellular matrix structural constituent (MF, P-value 2.58E-07), and hsa04512:ECM-receptor interaction (KEGG, P-value 2.22E-09), in downregulated group were GO:0007268~synaptic transmission (BP P-value 1.62E-31), GO:0045202~synapse (CC P-value 1.62E-35), GO:0022836~gated channel activity (MF P-value 6.73E-16), and hsa04080:neuroactive ligand-receptor interaction (KEGG P-value 3.52E-10).

Table 1: The GO and KEGG terms enriched by upregulated DEGs.

Category
Term
Count
P-Value

GOTERM_BP_FAT
GO:0030198~extracellular matrix organization
14
3.88E-12

GOTERM_BP_FAT
GO:0043062~extracellular structure organization
16
7.24E-12

GOTERM_BP_FAT
GO:0007049~cell cycle
28
2.27E-10

GOTERM_BP_FAT
GO:0000279~M phase
19
3.12E-10

GOTERM_BP_FAT
GO:0007067~mitosis
16
5.21E-10

GOTERM_CC_FAT
GO:0031012~extracellular matrix
21
6.61E-13

GOTERM_CC_FAT
GO:0005578~proteinaceous extracellular matrix
19
1.80E-11

GOTERM_CC_FAT
GO:0005576~extracellular region
42
1.61E-10

GOTERM_CC_FAT
GO:0044421~extracellular region part
28
8.32E-10

GOTERM_CC_FAT
GO:0044420~extracellular matrix part
12
9.61E-10

GOTERM_MF_FAT
GO:0005201~extracellular matrix structural constituent
9
2.58E-07

GOTERM_MF_FAT
GO:0048407~platelet-derived growth factor binding
4
6.59E-05

GOTERM_MF_FAT
GO:0004857~enzyme inhibitor activity
9
9.87E-04

GOTERM_MF_FAT
GO:0050840~extracellular matrix binding
4
0.001069556

GOTERM_MF_FAT
GO:0001871~pattern binding
7
0.001084076

KEGG_PATHWAY
hsa04512: ECM-receptor interaction
10
2.22E-09

KEGG_PATHWAY
hsa04510: Focal adhesion
10
4.37E-06

KEGG_PATHWAY
hsa04610: Complement and coagulation cascades
6
8.15E-05

KEGG_PATHWAY
hsa04110: Cell cycle
5
0.009122958

KEGG_PATHWAY
hsa04115: p53 signaling pathway
4
0.010179414

BP: biological process; CC: cellular component; MF: molecular function.

Table 2: The GO and KEGG terms enriched by downregulated DEGs.

Category
Term
Count
P-Value

GOTERM_BP_FAT
GO:0007268~synaptic transmission
52
1.62E-31

GOTERM_BP_FAT
GO:0019226~transmission of nerve impulse
54
5.18E-30

GOTERM_BP_FAT
GO:0007267~cell-cell signaling
57
4.96E-21

GOTERM_BP_FAT
GO:0006836~neurotransmitter transport
23
2.20E-18

GOTERM_BP_FAT
GO:0050877~neurological system process
72
3.01E-15

GOTERM_CC_FAT
GO:0045202~synapse
63
1.62E-35

GOTERM_CC_FAT
GO:0044456~synapse part
51
7.00E-32

GOTERM_CC_FAT
GO:0043005~neuron projection
50
4.79E-24

GOTERM_CC_FAT
GO:0008021~synaptic vesicle
25
1.28E-20

GOTERM_CC_FAT
GO:0030054~cell junction
49
1.50E-15

GOTERM_MF_FAT
GO:0022836~gated channel activity
36
6.73E-16

GOTERM_MF_FAT
GO:0005216~ion channel activity
39
3.23E-15

GOTERM_MF_FAT
GO:0015267~channel activity
40
4.98E-15

GOTERM_MF_FAT
GO:0022803~passive transmembrane transporter activity
40
5.43E-15

GOTERM_MF_FAT
GO:0022838~substrate specific channel activity
39
8.74E-15

KEGG_PATHWAY
hsa04080: Neuroactive ligand-receptor interaction
25
3.52E-10

KEGG_PATHWAY
hsa04020: Calcium signaling pathway
21
4.34E-10

KEGG_PATHWAY
hsa04720: Long-term potentiation
9
8.19E-05

KEGG_PATHWAY
hsa04012: ErbB signaling pathway
8
0.00232

KEGG_PATHWAY
hsa04010: MAPK signaling pathway
13
0.010508

BP: biological process; CC: cellular component; MF: molecular function.

3.3. PPI Network and Hub Genes

We mapped the overlapped DEGs to the STRING website, setting P<0.05, coefficients of as a strict threshold, constructed a co-expression network consisting of 401 nodes and 1890 edges (those disconnected genes in the network were deleted); Figure 3 displayed the network visualized by Cytoscape. The distribution of node degrees complied with exponential distribution; R square and correlation coefficient are 0.796 and 0.963 (Figure 4). Then we applied the cytoHubba application to screen out the hub genes in the network; the partial result calculated by multiform algorithm was shown in Table 3. Since the degree is the number of the edges of a gene in the network and represents the interaction pairs with others, it is more likely that the genes locate in a core position and act as a significant function; we adopted the top 20 genes ranked by degree as hub genes; they are CALB1, CDC20, CDCA8, CDK1, CEP55, DLGAP5, KIF20A, KIF4A, NDC80, PBK, RRM2, SYN1, SYP, SYT1, TPX2, TTK, VEGFA, BDNF, GNG3, and TOP2A. Figure 5 shows the interactions among the 20 hub genes.

Table 3: The statistical results of the connectivity of the network.

Gene
MCC
DMNC
MNC
Degree

TOP2A
9.22E+13
0.87325
40
56

CDK1
9.22E+13
0.68829
47
47

VEGFA
5763
0.18902
40
43

SYT1
4.37E+07
0.37425
40
42

CDC20
9.22E+13
1.01161
34
37

KIF4A
9.22E+13
0.98669
34
36

CALB1
4586
0.27262
33
35

BDNF
7597
0.33146
32
35

SYP
3760934
0.39583
33
35

DLGAP5
9.22E+13
1.1534
33
35

KIF20A
9.22E+13
1.05547
35
35

TTK
9.22E+13
1.10629
34
35

NDC80
9.22E+13
1.13243
33
35

GNG3
7267294
0.41708
32
34

CEP55
9.22E+13
1.07639
34
34

SYN1
3799485
0.41397
31
34

PBK
9.22E+13
1.09881
34
34

CDCA8
9.22E+13
1.1038
34
34

TPX2
9.22E+13
1.11408
33
34

RRM2
9.22E+13
1.03281
33
34

Figure 3: Protein-protein interaction work of overlapped DEGs. Nodes with higher degree are diplayed in larger size and bright orange colour and nodes with lower degree are shown in smaller size and dark blue colour. The edge size is consistent with the coexpression intensity.

Figure 4: Scatter plot of node degree distribution for overlapped DEGs, R square=0.796, and correlation=0.963.

Figure 5: Protein-protein interaction work for 20 hub genes. The shape of triangle represents upregulated genes, the shape of “V” represents downregulated genes.

3.4. Function and Survival Analysis of the Hub Genes

By searching the DAVID website, we found the hub genes mainly enriched in the following GO terms: M phase; nuclear division; mitosis; M phase of mitotic cell cycle; organelle fission; microtubule cytoskeleton; mitotic cell cycle; spindle; cell cycle phase; and cell cycle process. CDK1, TTK, CDC20, and RRM2 were found to participate in the cell cycle and p53 signaling pathway (Table 4). Then we utilized the univariate Cox regression analysis and log rank test to draw the survival curves to explore the prognostic related value of the hub genes based on TCGA datasets. Samples with survival information were divided into high expression group and low expression according to the median expression of the candidate genes. CEP55 and RRM2 were discovered to have survival differences (P<0.05) (Figures 6(a) and 6(b)); highly expressing CEP55 and RRM2 tends to have poor survival outcomes. To test and verify the result, we did survival analysis on GSE74187, another GEO datasets which contain 60 GBM samples with overall survival (OS) and progress-free survival (FPS) time, with the same calculation strategy; it indicated that the expression level of RRM2 might be related to both OS and PFS (P<0.05) (Figures 6(c) and 6(d)), which is consistent with earlier investigation. However, the CEP55 effect cannot be seen in GSE74187, maybe because of the fact that interference comes from limited samples or other confounding factors, so we choose RRM2 for the next step of exploration.

Table 4: The GO and KEGG terms enriched by 20 hub genes.

Category
Term
count
P-Value
Genes

GOTERM_BP_FAT
GO:0000279~M phase
9
6.73E-09
CDK1, CDCA8, DLGAP5, TPX2, TTK, CDC20, NDC80, CEP55, PBK

GOTERM_BP_FAT
GO:0000280~nuclear division
8
1.17E-08
CDK1, CDCA8, DLGAP5, TPX2, CDC20, NDC80, CEP55, PBK

GOTERM_BP_FAT
GO:0007067~mitosis
8
1.17E-08
CDK1, CDCA8, DLGAP5, TPX2, CDC20, NDC80, CEP55, PBK

GOTERM_BP_FAT
GO:0000087~M phase of mitotic cell cycle
8
1.32E-08
CDK1, CDCA8, DLGAP5, TPX2, CDC20, NDC80, CEP55, PBK

GOTERM_BP_FAT
GO:0048285~organelle fission
8
1.54E-08
CDK1, CDCA8, DLGAP5, TPX2, CDC20, NDC80, CEP55, PBK

GOTERM_CC_FAT
GO:0015630~microtubule cytoskeleton
10
1.60E-08
CDK1, CDCA8, KIF4A, DLGAP5, TPX2, TTK, CDC20, CEP55, TOP2A, KIF20A

GOTERM_BP_FAT
GO:0000278~mitotic cell cycle
9
1.69E-08
CDK1, CDCA8, DLGAP5, TPX2, TTK, CDC20, NDC80, CEP55, PBK

GOTERM_CC_FAT
GO:0005819~spindle
7
3.46E-08
CDK1, CDCA8, KIF4A, DLGAP5, TPX2, TTK, CDC20

GOTERM_BP_FAT
GO:0022403~cell cycle phase
9
4.04E-08
CDK1, CDCA8, DLGAP5, TPX2, TTK, CDC20, NDC80, CEP55, PBK

GOTERM_BP_FAT
GO:0022402~cell cycle process
9
4.43E-07
CDK1, CDCA8, DLGAP5, TPX2, TTK, CDC20, NDC80, CEP55, PBK

KEGG_PATHWAY
hsa04110: Cell cycle
3
0.008429
CDK1, TTK, CDC20

KEGG_PATHWAY
hsa04115: p53 signaling pathway
2
0.077638
CDK1, RRM2

BP: biological process; CC: cellular component; MF: molecular function.

(a)

(b)

(c)

(d)

Figure 6: Kaplan-Meier survival curves of GBM patients grouped by median expression level of relevant gene. (a) CEP55 in TCGA. (b) RRM2 in TCGA. (c. d) RRM2 in GSE74187.

3.5. Correlation between RRM2 and Other Hub Genes

Correlation analysis implied that RRM2 expression was positively correlated with DLGAP5 (R=0.79, P<0.05), KIF4A (R=0.78, P<0.05), CDCA8 (R=0.75, P<0.05), TPX2 (R=0.75, P<0.05), KIF20A (R=0.72,P<0.05), PBK (R=0.72,P<0.05), CDC20 (R=0.68, P<0.05), CDK1 (R=0.68, P<0.05), NDC80 (R=0.68, P<0.05), TTK (R=0.68, P<0.05), CEP55 (R=0.62, P<0.05), and TOP2A (R=0.58, P<0.05), while there was insignificant correlation with SYP, VEGFA, BDNF,CALB1, SYT1, GNG3, and SYN1 (supplementary material 1).

3.6. Exploration of RRM2

According to the GEPIA website (http://gepia.cancer-pku.cn/), we observed that RRM2 was differently expressed in many types of tumor tissues including GBM compared to normal tissues (Figure 7(a)), which further indicated it may contribute to the occurrence and development of cancer. Figure 7(b) shows the different expression level of RRM2 between the GBM samples from TCGA and normal brain samples matched from GTEx datasets (https://gtexportal.org/home/), which contain more abundant normal samples compared to TCGA database. Apparently GBM patients expressed higher level of RRM2. Besides, RRM2 were highly expressed in GBM cancer cell lines (Figure 7(c)) as demonstrated from CCLE analysis (https://portals.broadinstitute.org/ccle/), while the prognostic was worse according to previous research. To explore the probable regulated mechanisms of RRM2 in GBM, we explored the cBioPortal (http://www.cbioportal.org/) and found that RRM2 was altered in 4% of GBM cases from TCGA cohort; the main type of genetic alteration is mRNA upregulation. Correspondingly, those with RRM2 altered seem to have poor outcome (Figure 8). We annotated the biological processes of RRM2 by employing GeneMANIA (http://genemania.org/), a freely available and effective web interface for functional prediction of genes; as is shown in Figure 9, RRM2 interacted with 20 proteins and was mainly involved in the following functions: nucleobase-containing small molecule interconversion; regulation of transcription involved in G1/S transition of mitotic cell cycle; G1/S transition of mitotic cell cycle; and deoxyribonucleotide metabolic process. Considering that, we infer that RRM2 may affect the processes of GBM via regulating mitotic cell cycle and nucleotide metabolism.

(a)

(b)

(c)

Figure 7: The expression level of RRM2 in different types of cancers and cell lines. (a) In pan-cancer tissues and normal tissues. (b) In GBM tissues and normal tissues. (c) In multitumor cell lines.

(a)

(b)

Figure 8: Kaplan-Meier survival curves of GBM patients grouped by RRM2 alteration according to cBioPortal. (a). OS curves. (b). PFS curves.

Figure 9: PPI network of RRM2 generated by using GeneMANIA. The network edges refer to the interaction types between gene pairs and the colour filled in the nodes refers to a pathway the gene node participated in.

4. Discussion

Glioblastoma (GBM) is the most malignant glial tumor in brain with devastating prognosis because of its complex biological behaviors and limited strategies for therapy; people with GBM suffer from tumor progression or recurrence. Recently, advances in genomic analysis help us form a comprehensive insight into the mechanisms of GBM and identify potential biomarkers for GBM prognosis and therapy. In the present study, we performed gene expression profile analysis among three datasets: TCGA, GSE4290, and GSE50161, totally included 284 GBM and 41 normal samples. At last, the overlapped 552 DEGs among three datasets were identified, of which 128 were upregulated and 424 were downregulated. GO annotation indicated that the upregulated genes were mainly manifested in cell cycle, mitosis, and extracellular structure organization; the downregulated genes were primarily related to neurological system process, cell-cell signaling, and synaptic transmission. KEGG pathway shows that the upregulated genes were associated with ECM-receptor interaction, focal adhesion, complement and coagulation cascades, cell cycle, and p53 signaling pathway; the downregulated genes were associated with neuroactive ligand-receptor interaction, calcium signaling pathway, long-term potentiation, ErbB signaling pathway, and MAPK signaling pathway. These results indicated that the pathogenesis of tumor is a complex biological process driven by specific genes and epigenetic changes. Abnormal regulation of multiple genes will promote the occurrence and development of GBM through different mechanisms.

To further screen core genes, we constructed a protein-protein interaction network consisting of 401 nodes and 1890 edges through STRING; then we focus our attention on the 20 core genes in the network picked out by cytoHubba according to their degree. GO annotation indicated the 20 hubgenes mainly took part in cell mitosis and effect the cell cycle; among them, CDK1, TTK, and CDC20 were enriched in cell cycle pathway; CDK1 and RRM2 were enriched in p53 signaling pathway. As is known to all, those pathways play a vital role in carcinoma. p53 signaling pathway is composed of a set of genes and their products that response to a variety of intrinsic and extrinsic stress signals in apoptosis, cellular senescence, or cell cycle arrest style [12]; dysregulation of the genes in the p53 network might disrupt the fidelity of DNA replication and cell division and present a greater risk of carcinogenesis [13]. Evidently, rapid proliferation is indispensable for tumor cells to maintain growth and invasion; disorders in the regulation of cell cycle contribute to the initiation and progression of tumor [14]. The hallmarks of rapid proliferation and high aggressiveness of GBM make identifying cell cycle pathway and p53 signaling pathway as a predominant signalling pathway in GBM reasonable. Further unravelling of the precise molecular mechanism about those pathways and relevant genes may be necessary for future efforts.

Present survival analysis revealed that RRM2, one of the hub genes in the coexpression network, was associated with the prognosis of GBM patients; patients with high expression of RRM2 showed shorter overall survival time than those with low expression level. Besides, 4% of TCGA GBM samples detected out RRM2 alteration according to cBioPortal; however, this alteration seems to relate to poor outcome. Hence, we hypothesized that RRM2 could be a potential prognostic factor of patients with GBM. Ribonucleotide reductase (RR), a rate-limiting enzyme which catalyzes the formation of deoxyribonucleotides from ribonucleotides, consists of ribonucleotide reductase subunit M1 (RRM1) and ribonucleotide reductase subunit M2 (RRM2). Over the years, studies continually found that RRM2 was overexpressed in various tumors, such as pancreatic cancer [15], neuroblastoma [16], thyroid malignant tumor [17], breast cancer [18], melanoma [19], lung cancer [20], prostate carcinoma [21], and hepatocellular carcinoma [22]; elevated RRM2 activity played a vital role in tumorigenesis, progression, and invasion and can serve as prognostic markers in some instances. Inhibition of RR can break the balance of ribonucleotide and deoxyribonucleotide level, effecting the DNA synthesis and repair and inducting cell cycle arrest and apoptosis [23–25], making RR an attracting cancer therapeutic target and encouraged incessant exploring. Gemcitabine was one of the famous RRM1 inhibitors and was approved clinically as first-line drug for anticancer therapy in various cancers. There are other efficient RRM2 inhibitors such as hydroxyurea (HU), already used in the treatment of hematological malignances [26, 27], and Triapine (3AP), evaluated in a number of clinical trials and showed encouraging results in anticancer treatment [28, 29]. However, RRM2 is underexplored and recognized as a viable therapeutic option. With the recent advances in genetic strategies for RRM2 inhibition, some new treatments were reported and highlight the potential of gene therapy. Knockdown of RRM2 through specific small interfering RNA (siRNA) displayed effective antitumor activity in various solid tumors, like head and neck cancer [30], including oral squamous cell carcinoma [31], ovarian cancer [32], gastric adenocarcinoma [33], hepatocellular carcinoma [34], and colorectal cancers [35]. What is more is that therapies targeting or decreasing RRM2 expression through antisense cDNA [36] and RRM2-specific siRNA [37] have displayed a reversal of gemcitabine resistance. Even though the extensive exploration of RRM2 in GBM was rare, Li C. et al. illustrated that RRM2 was overexpressed in human GBM cells and can promote proliferation, migration, and invasion but inhibit apoptosis of GBM cells at experiment level [38]; however, the study only focused on RRM2 simply and the signaling pathway about how RRM2 was involved in was not mentioned. Only Rasmussen et al. demonstrated the RRM2-BRCA1 interaction mechanism and found that BRCA1 protects GBM cells from endogenous replication stress and promotes tumorigenicity [39]. Consistent with previous study, we found RRM2 was strongly high expressed in GBM samples compared to normal brain tissues via comprehensive bioinformation analysis; meanwhile RRM2 was a negative prognostic factor informing worse clinical survival. Importantly, we noted RRM2 may play a role in GBM mainly by effecting cell cycle and nucleotide metabolism; knockdown RRM2 could suppress the proliferation of GBM cells, induce cell arrest in the G1 phase, and promote cell apoptosis [38], supporting present hypothesis. Notably, we identified 20 core genes which were considered to play a central role in GBMs, among which, CDK1, CDC20, CDCA8, CEP55, DLGAP5, KIF20A, KIF4A, NDC80, PBK, TPX2, TTK, VEGFA, and TOP2A had a direct interaction with RRM2; subsequent correlation analysis showed that RRM2 was highly positively correlated with DLGAP5, KIF4A, CDCA8, and TPX2 (R≥0.75, P<0.05). DLGAP5 (Discs Large Homolog Associated Protein 5), also known as HURP or DLG7, is a mitotic spindle protein that promotes the formation of tubulin polymers [40]; it has been documented that DLGAP5 may lead to cancer by allowing cells to progress through the regulation of M phase progression by modulating the function of the spindle apparatus and its organization [41]. KIF4A (Kinesin Family Member 4A) is a microtubule-based motor protein involved in maintaining chromosome integrity during cell mitosis [42, 43]; it was implicated as an potential biomarker in types of cancers [44, 45]. CDCA8 (Cell Division Cycle Associated 8) is one of the components of chromosomal passenger complex (CPC) in mitosis and cell division; in fact, CDCA8 was considered as a putative oncogene for it was overexpressed in tumor tissues but had low or undetectable expression in normal tissues [46, 47]. TPX2 (Microtubule Nucleation Factor) is well known for its critical role in mitotic spindle assembly [48]; a number of studies displayed the efficiency of the strategy of decreasing TPX2 level in anticancer treatment [49]. Moreover, recent works on TPX2 in DNA damage response opened an extended therapeutic window for TPX2-targeted therapies in cancer [50]. Apparently, these genes take part in cell mitosis and cell cycle and are cancer related, even though precise mechanisms of how this PPI network worked are confused; this information adds value to the notion that the RRM2-related network is relevant in carcinogenesis and deserves further exploration in GBM; it may provide a novel avenue worth pursuing when developing effective drugs by combining RRM2 inhibitors (such as 3AP) with other agents. Collecting all, we confirmed that RRM2 is a prognostic factor and a promising therapeutic target for GBM treatment.

5. Conclusion

In conclusion, the present study presented a comprehensive bioinformatics analysis of DEGs between glioblastoma and normal tissues and successfully screened several crucial genes and certain associated pathways such as cell cycle and p53 signaling pathway. Our data suggested that pathogenesis and development of glioblastoma are regulated by complex gene network through different biological pathways; RRM2 is a core gene in this network and is associated with the prognosis of glioblastoma; it could be a promising target for efficient suppression of glioblastoma; however, further experiments and more efforts are needed to confirm it and illustrate the specific molecular biological mechanism.

Data Availability

The data used for analysis in this study are available from the Cancer Genome Atlas and the Gene Expression Omnibus database freely.

Conflicts of Interest

All authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported by the National Nature Science Foundation of China (81772680).

Supplementary Materials

Supplementary Figure 1: correlation between the expression of RRM2 and other hub genes. Supplementary Table 1: correlation between the expression of RRM2 and other hub genes. (Supplementary Materials)

References

D. N. Louis, A. Perry, G. Reifenberger et al., “The 2016 world health organization classification of tumors of the central nervous system: a summary,” Acta Neuropathologica, vol. 131, no. 6, pp. 803–820, 2016.

Q. T. Ostrom, H. Gittleman, P. Farah et al., “CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2006–2010,” Neuro-Oncology, vol. 15, supplement 2, pp. ii1–ii56, 2013.

R. Stupp, M. E. Hegi, W. P. Mason et al., “Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial,” Lancet Oncology, vol. 10, no. 5, pp. 459–466, 2009.

M. Weller, M. van den Bent, K. Hopkins et al., “EANO guideline for the diagnosis and treatment of anaplastic gliomas and glioblastoma,” The Lancet Oncology, vol. 15, no. 9, pp. e395–e403, 2014.

E. G. Van Meir, C. G. Hadjipanayis, A. D. Norden, H.-K. Shu, P. Y. Wen, and J. J. Olson, “Exciting new advances in neuro-oncology: the avenue to a cure for malignant glioma,” CA: A Cancer Journal for Clinicians, vol. 60, no. 3, pp. 166–193, 2010.

T. C. Ryken, S. N. Kalkanis, J. M. Buatti, and J. J. Olson, “The role of cytoreductive surgery in the management of progressive glioblastoma,” Journal of Neuro-Oncology, vol. 118, no. 3, pp. 479–488, 2014.

V. Kulasingam and E. P. Diamandis, “Strategies for discovering novel cancer biomarkers through utilization of emerging technologies,” Nature Clinical Practice Oncology, vol. 5, no. 10, pp. 588–599, 2008.

P. Dong, B. Yu, L. Pan, X. Tian, and F. Liu, “Identification of key genes and pathways in triple-negative breast cancer by integrated bioinformatics analysis,” BioMed Research International, vol. 2018, Article ID 2760918, pp. 1–10, 2018.

L. Zhuang, Z. Yang, and Z. Meng, “Upregulation of BUB1B, CCNB1, CDC7, CDC20, and MCM3 in tumor tissues predicted worse overall survival and disease-free survival in hepatocellular carcinoma patients,” BioMed Research International, vol. 2018, Article ID 7897346, pp. 1–8, 2018.

D. W. Huang, B. T. Sherman, and R. A. Lempicki, “Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources,” Nature Protocols, vol. 4, no. 1, pp. 44–57, 2009.

D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, pp. D447–D452, 2015.

S. Jin and A. J. Levine, “The p53 functional circuit,” Journal of Cell Science, vol. 114, no. 23, pp. 4139-4140, 2001.

B. Vogelstein, D. Lane, and A. J. Levine, “Surfing the p53 network,” Nature, vol. 408, no. 6810, pp. 307–310, 2000.

M. Malumbres and M. Barbacid, “Cell cycle, CDKs and cancer: a changing paradigm,” Nature Reviews Cancer, vol. 9, no. 3, pp. 153–166, 2009.

H. Sun, L. Zhao, K. Pan, Z. Zhang, M. Zhou, and G. Cao, “Integrated analysis of mRNA and miRNA expression profiles in pancreatic ductal adenocarcinoma,” Oncology Reports, vol. 37, no. 5, pp. 2779–2786, 2017.

J. Li, J. Pang, Y. Liu et al., “Suppression of RRM2 inhibits cell proliferation, causes cell cycle arrest and promotes the apoptosis of human neuroblastoma cells and in human neuroblastoma RRM2 is suppressed following chemotherapy,” Oncology Reports, vol. 40, no. 1, pp. 355–360, 2018.

E. Castelblanco, C. Zafon, J. Maravall et al., “APLP2, RRM2, and PRC1: New Putative Markers for the Differential Diagnosis of Thyroid Follicular Lesions,” Thyroid, vol. 27, no. 1, pp. 59–66, 2017.

R. Bell, R. Barraclough, and O. Vasieva, “Gene expression meta-analysis of potential metastatic breast cancer markers,” Current Molecular Medicine, vol. 17, no. 3, pp. 200–210, 2017.

N. Fatkhutdinov, K. Sproesser, C. Krepler et al., “Targeting RRM2 and mutant BRAF is a novel combinatorial strategy for melanoma,” Molecular Cancer Research, vol. 14, no. 9, pp. 767–775, 2016.

F. Grossi, M. G. Dal Bello, S. Salvi et al., “Expression of ribonucleotide reductase subunit-2 and thymidylate synthase correlates with poor prognosis in patients with resected stages I–III non-small cell lung cancer,” Disease Markers, vol. 2015, Article ID 302649, 18 pages, 2015.

Z. He, F. Tang, Z. Lu et al., “Analysis of differentially expressed genes, clinical value and biological pathways in prostate cancer,” American Journal of Translational Research, vol. 10, no. 5, pp. 1444–1456, 2018.

B. Lee, S. Y. Ha, D. H. Song, H. W. Lee, S. Y. Cho, and C. Park, “High expression of ribonucleotide reductase subunit m2 correlates with poor prognosis of hepatocellular carcinoma,” Gut and Liver, vol. 8, no. 6, pp. 662–668, 2014.

D. A. Albert and L. J. Gudas, “Ribonucleotide reductase activity and deoxyribonucleoside triphosphate metabolism during the cell cycle of S49 wild-type and mutant mouse T-lymphoma cells,” The Journal of Biological Chemistry, vol. 260, no. 1, pp. 679–684, 1985.

J. Guo, Q. Chen, C. W. Lam et al., “Profiling ribonucleotide and deoxyribonucleotide pools perturbed by gemcitabine in human non-small cell lung cancer cells,” Scientific Reports, vol. 6, no. 1, Article ID 37250, 2016.

H. Fan, C. Villegas, A. Huang, and J. A. Wright, “The mammalian ribonucleotide reductase R2 component cooperates with a variety of oncogenes in mechanisms of cellular transformation,” Cancer Research, vol. 58, no. 8, pp. 1650–1653, 1998.

R. Hehlmann, H. Heimpel, J. Hasford et al., “Randomized comparison of busulfan and hydroxyurea in chronic myelogenous leukemia: prolongation of survival by hydroxyurea. The German CML Study Group,” in Blood, vol. 82, pp. 398–407, 2 edition, 1993.

Y. Sterkers, C. Preudhomme, J.-L. Laï et al., “Acute myeloid leukemia and myelodysplastic syndromes following essential thrombocythemia treated with hydroxyurea: high proportion of cases with 17p deletion,” Blood, vol. 91, no. 2, pp. 616–622, 1998.

C. M. Nutting, C. M. L. Van Herpen, A. B. Miah et al., “Phase II study of 3-AP Triapine in patients with recurrent or metastatic head and neck squamous cell carcinoma,” Annals of Oncology, vol. 20, no. 7, pp. 1275–1279, 2009.

J. F. Zeidner, J. E. Karp, A. L. Blackford et al., “A phase II trial of sequential ribonucleotide reductase inhibition in aggressive myeloproliferative neoplasms,” Haematologica, vol. 99, no. 4, pp. 672–678, 2014.

M. A. Rahman, A. R. M. R. Amin, X. Wang et al., “Systemic delivery of siRNA nanoparticles targeting RRM2 suppresses head and neck tumor growth,” Journal of Controlled Release, vol. 159, no. 3, pp. 384–392, 2012.

K. Iwamoto, K.-I. Nakashiro, H. Tanaka, N. Tokuzen, and H. Hamakawa, “Ribonucleotide reductase M2 is a promising molecular target for the treatment of oral squamous cell carcinoma,” International Journal of Oncology, vol. 46, no. 5, pp. 1971–1977, 2015.

M. Zhang, J. Wang, R. Yao, and L. Wang, “Small interfering RNA (siRNA)-mediated silencing of the M2 subunit of ribonucleotide reductase: a novel therapeutic strategy in ovarian cancer,” International Journal of Gynecological Cancer, vol. 23, no. 4, pp. 659–666, 2013.

W. Kang, J. H. M. Tong, A. W. H. Chan et al., “Targeting ribonucleotide reductase M2 subunit by small interfering RNA exerts anti-oncogenic effects in gastric adenocarcinoma,” Oncology Reports, vol. 31, no. 6, pp. 2579–2586, 2014.

J. Gao, H. Chen, Y. Yu et al., “Inhibition of hepatocellular carcinoma growth using immunoliposomes for co-delivery of adriamycin and ribonucleotide reductase M2 siRNA,” Biomaterials, vol. 34, no. 38, pp. 10084–10098, 2013.

X. Liu, H. Zhang, L. Lai et al., “Ribonucleotide reductase small subunit M2 serves as a prognostic biomarker and predicts poor survival of colorectal cancers,” Clinical Science, vol. 124, no. 9, pp. 567–578, 2013.

Y. Lee, A. Vassilakos, N. Feng et al., “GTI-2040, an antisense agent targeting the small subunit component (R2) of human ribonucleotide reductase, shows potent antitumor activity against a variety of tumors,” Cancer Research, vol. 63, no. 11, pp. 2802–2811, 2003.

M. S. Duxbury, H. Ito, M. J. Zinner, S. W. Ashley, and E. E. Whang, “RNA interference targeting the M2 subunit of ribonucleotide reductase enhances pancreatic adenocarcinoma chemosensitivity to gemcitabine,” Oncogene, vol. 23, no. 8, pp. 1539–1548, 2004.

C. Li, J. Zheng, S. Chen et al., “RRM2 promotes the progression of human glioblastoma,” Journal of Cellular Physiology, vol. 233, no. 10, pp. 6759–6767, 2018.

R. D. Rasmussen, M. K. Gajjar, L. Tuckova et al., “BRCA1-regulated RRM2 expression protects glioblastoma cells from endogenous replication stress and promotes tumorigenicity,” Nature Communications, vol. 7, p. 13398, 2016.

R. A. Santarella, M. D. Koffa, P. Tittmann, H. Gross, and A. Hoenger, “HURP wraps microtubule ends with an additional tubulin sheet that has a novel conformation of tubulin,” Journal of Molecular Biology, vol. 365, no. 5, pp. 1587–1595, 2007.

A.-P. Tsou, C.-W. Yang, C.-Y. F. Huang et al., “Identification of a novel cell cycle regulated gene, HURP, overexpressed in human hepatocellular carcinoma,” Oncogene, vol. 22, no. 2, pp. 298–307, 2003.

Y. M. Lee, S. Lee, E. Lee et al., “Human kinesin superfamily member 4 is dominantly localized in the nuclear matrix and is associated with chromosomes during mitosis,” Biochemical Journal, vol. 360, no. 3, pp. 549–556, 2001.

C. Zhu and W. Jiang, “Cell cycle-dependent translocation of PRC1 on the spindle by Kif4 is essential for midzone formation and cytokinesis,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 102, no. 2, pp. 343–348, 2005.

M. Taniwaki, A. Takano, N. Ishikawa et al., “Activation of KIF4A as a prognostic biomarker and therapeutic target for lung cancer,” Clinical Cancer Research, vol. 13, no. 22, pp. 6624–6631, 2007.

G. Hou, C. Dong, Z. Dong et al., “Upregulate KIF4A enhances proliferation, invasion of hepatocellular carcinoma and indicates poor prognosis across human cancer types,” Scientific Reports, vol. 7, article 4148, 2017.

J.-L. Chang, T.-H. Chen, C.-F. Wang et al., “Borealin/Dasra B is a cell cycle-regulated chromosomal passenger protein and its nuclear accumulation is linked to poor prognosis for human gastric cancer,” Experimental Cell Research, vol. 312, no. 7, pp. 962–973, 2006.

S. Hayama, Y. Daigo, T. Yamabuki et al., “Phosphorylation and activation of cell division cycle associated 8 by aurora kinase B plays a significant role in human lung carcinogenesis,” Cancer Research, vol. 67, no. 9, pp. 4113–4122, 2007.

P. Wadsworth, “TPX2,” Current Biology, vol. 25, no. 24, pp. R1156–R1158, 2015.

R. Satow, M. Shitashige, Y. Kanai et al., “Combined functional genome survey of therapeutic targets for hepatocellular carcinoma,” Clinical Cancer Research, vol. 16, no. 9, pp. 2518–2528, 2010.

G. Neumayer, C. Belzil, O. J. Gruss, and M. D. Nguyen, “TPX2: of spindle assembly, DNA damage response, and cancer,” Cellular and Molecular Life Sciences, vol. 71, no. 16, pp. 3027–3047, 2014.