Prognostic Markers Identification in Glioma by Gene Expression Profile Analysis

Abstract

This study aimed to explore more gene markers associated with glioma or its prognosis. The glioma-related RNAseq data from the Gene Expression Omnibus database and The Cancer Genome Atlas dataset in UCSC Xena database were downloaded. There was a total of 971 tumor samples and 102 normal samples in the 2 datasets. The differentially expressed genes (DEGs) data between tumor and normal samples were analyzed, on which were then performed function and pathway enrichment analyses. Pearson correlation coefficient between DEGs was calculated to construct the coexpression network. Finally, prognostic genes were screened. A total of 634 upregulated and 769 downregulated DEGs were identified between tumor and control groups. These DEGs were significantly involved in 15 upregulated pathways, such as p53 signaling pathway, and 16 downregulated pathways, such as neuroactive ligand–receptor interaction, and cell adhesion molecules. In the coexpression network, pseudouridine synthase 7 (PUS7), EFR3 homolog B (EFR3B), and neuronal cell adhesion molecule (NRCAM) had the top three highest degrees. Additionally, 17 prognostic genes were selected, such as thrombospondin-1 (THBS1), caspase-8 (CASP8), glutamate ionotropic receptor AMPA type subunit 2 (GRIA2), GRIA4, and ADCYAP receptor type I (ADCYAP1R1). Pathways of p53 signaling pathway and neuroactive ligand–receptor interaction may play important roles in glioma progression. PUS7, EFR3B, and NRCAM may be potential biomarkers of glioma. THBS1, CASP8, GRIA2, GRIA4, and ADCYAP1R1 may serve as prognostic markers in glioma.

1. Introduction

Gliomas are primary brain cancers that arise from glial cells, which account for the majority of malignant brain tumors in adults (Goodenberger and Jenkins, 2012). Gliomas are traditionally classified as astrocytic, ependymal, or oligodendroglial tumors based on their histological appearance (Ferluga and Debinski, 2014; Weller et al., 2015). The astrocytic tumors are the most common and include the most malignant type of glioblastoma that had a very poor survival time of ∼15 months (CBTRUS, 2011). It is necessary to advance our mechanistic understanding of glioma tumorigenesis.

Like cancer in general, gliomas develop due to genetic alterations that accumulate with tumor progression (Riemenschneider et al., 2010). Genome-wide molecular-profiling studies have revealed comprehensive mutational landscapes for the major types of gliomas (Cancer Genome Atlas Research Network et al., 2015; Korshunov et al., 2015; Ceccarelli et al., 2016). For instance, mutation of the isocitrate dehydrogenase 1 (IDH1) gene, or the related IDH2 gene, has been identified in the majority of world health organization grades II and III astrocytic, oligodendroglial, and oligoastrocytic gliomas (Balss et al., 2008; Yan et al., 2009). Mutation of the tumor suppressor gene, TP53, has been found in more than half of the grade II diffuse astrocytomas (Riemenschneider et al., 2010). Moreover, primary glioblastoma shows frequent EGFR amplification and PTEN mutation, and secondary glioblastoma is characterized by frequent mutations in the TP53 and IDH1 genes (Kleihues and Ohgaki, 1999; Ohgaki and Kleihues, 2007). These molecular characteristics can be used to improve the prediction of outcomes of patients, and to guide individualized treatment (Reifenberger et al., 2016).

In this study, we downloaded the glioma-related RNAseq datasets from the Gene Expression Omnibus (GEO) database and The Cancer Genome Atlas (TCGA) database to explore more gene markers associated with glioma. Additionally, we performed survival analysis to screen prognostic genes of glioma.

2. Methods

2.1. Data collection

We downloaded the glioma-related RNAseq data from the GEO database (Barrett et al., 2013; www.ncbi.nlm.nih.gov/geo) and TCGA dataset in UCSC Xena database (Haeussler et al., 2019; https://xenabrowser.net/datapages). The GEO dataset (GSE57872) contained 441 tumor samples and 102 normal samples. On this dataset was performed the gene expression quantification using RSEM method. The RNAseq data (IlluminaHiSeq pancan normalized) downloaded from UCSC Xena database included 530 tumor samples.

2.2. Differential expression analysis

There was a total of 971 tumor samples and 102 normal samples in the 2 datasets. After normalization of the data matrixes using betaqn method in R, the differentially expressed genes (DEGs) data between tumor and normal samples were analyzed using limma package (version 3.10.3; Ritchie et al., 2015). Paired t test was used for significance test to obtain the p value. Benjamini–Hochberg (BH) method was used for multiple test adjustment to obtain the adjusted p value. The thresholds for DEGs were set as adjusted p value <0.05 and |log fold change (FC)| > 1.

2.3. Function and pathway enrichment analyses of DEGs

The functional annotations, including biological process (BP), molecular function (MF), and cellular component (CC) enriched by the DEGs, were analyzed using Database for Annotation, Visualization, and Integrated Discovery (Huang da et al., 2009; version 6.7, https://david-d.ncifcrf.gov). Results were visualized using GOplot (Walter et al., 2015). Based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database, gene set enrichment analysis (GSEA, version 3.0; Damian and Gorfine, 2004) was used to conduct KEGG pathway analysis. Adjusted p value <0.05 was used as the threshold.

2.4. Coexpression analysis of genes

Pearson correlation coefficient between DEGs was calculated, and the correlation test was performed using the corr.test method (parameter ci = F, adjust = “BH”) of psych (Jason, 2013) in R. BH method was used for multiple test adjustment. Gene pairs with r > 0.8 and p value <0.05 were screened and Cytoscape software (Shannon et al., 2003; version 3.7.0) was used to construct the coexpression network.

2.5. Prognostic gene screening

The prognostic information of patients, including overall survival (OS) and OS status were collated according to the downloaded clinical data. Based on the obtained DEGs as the candidate genes, the patients were divided into high-expression and low-expression group according to the median of gene expression value. Then combining the prognostic information, the Kaplan–Meier (K-M) survival curve was drawn with K-M survival analysis method. The significance p value was calculated using the logrank test. The gene with p value <0.05 was preliminarily considered to be associated with prognosis.

3. Results

3.1. DEGs identification

The number of annotated genes in the TCGA dataset and GEO dataset were 20,530 and 5948, respectively. There were a total of 5678 common genes in the two datasets, among which there were 1403 DEGs (634 upregulated and 769 and downregulated; Fig. 1A). The top 10 upregulated and downregulated genes are shown in Figure 1B.

FIG. 1.

(A) Volcano map of DEGs. The thresholds were |log₂(FC)| > 1 and adjusted p value <0.05. Black dots indicate upregulated genes and grey dots indicate downregulated genes. (B) Heatmap of top 10 upregulated and downregulated genes (according to an ascending order of adjusted p value). FC, fold change; DEGs, differentially expressed genes.

3.2. Gene ontology enrichment analysis

Gene ontology (GO) enrichment analysis was performed on these DEGs, and the results showed that a total of 24 BP, 23 CC, and 5 MF were significantly enriched. As shown in Figure 2, DEGs were significantly enriched in functions related to cell projection organization, morphogenesis, vesicle, membrane-bound vesicle, cytoskeletal protein, and actin. These functions may be related to aberration of cancer cells, acquisition of foreign substances, intercellular communication, and metastasis of cancer cells.

FIG. 2.

Function enrichment analysis results. (A) The most significant five terms were selected for BP, MF, and CC, respectively. (B) For the top 8 GO terms selected, the phylogenetic tree was constructed by clustering the genes according to their occurrence in the term. BP, biological process; CC, cellular component; FDR, false discovery rate; GO, gene ontology; MF, molecular function.

3.3. KEGG pathway enrichment analysis

According to the enrichment analysis by GSEA, a total of 31 significantly enriched pathways we identified, including 15 upregulated pathways [normalized enrichment score (NES) >0] and 16 downregulated ones (NES <0; Fig. 3A, B). p53 signaling pathway had the largest NES enrichment score among the upregulated pathways (Fig. 3C). Additionally, cell cycle pathway, and small-cell lung cancer were also significantly upregulated. Among the downregulated pathways, neuroactive ligand–receptor interaction had the highest NES (Fig. 3D). Many signal transduction pathways, such as phospholipase D signaling pathway, peroxisome proliferators-activated receptor (PPAR) signaling pathway, and cyclic adenosine monophosphate (CAMP) signaling pathway were also downregulated.

FIG. 3.

Enrichment results of GSEA pathway. (A) The dotplot of partially enriched pathways. (B) The joyplot of the upregulated and downregulated pathways. (C) The upregulated KEGG_P53_SIGNALING_PATHWAY. (D) The downregulated KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.4. Gene coexpression network

For the 1403 DEGs, we calculated the correlation coefficient and significance degree between two genes based on their expression amount in different samples, and screened the gene pairs with r > 0.8 and adjusted p value <0.05. Cystoscape was used for visualized analysis on the gene pairs, as shown in Figure 4. There were 119 nodes (genes) in this coexpression network, among which 54 genes were upregulated and 65 genes were downregulated. As can be seen from Figure 4, there were 5 connected networks in the coexpression network, and the corresponding gene numbers were 2, 3, 3, 6, and 105, respectively. In this coexpression network, pseudouridine synthase 7 (PUS7), EFR3 homolog B (EFR3B), and neuronal cell adhesion molecule (NRCAM) had high connectivity degree (hub genes).

FIG. 4.

Gene coexpression network. Dark grey indicates upregulated genes and light grey indicates downregulated genes.

3.4. Prognosis analysis

The downloaded TCGA clinical data contained the information of 515 patients. The information of survival status and corresponding survival time was extracted, and the mismarked cases were deleted, and finally the clinical information of 514 patients was used. In GO and GSEA enrichment analyses, respectively, 779 and 884 DEGs were analyzed, and 286 genes were common genes. According to logFC values of these genes, the top 10 upregulated and downregulated genes were selected. The K-M survival curves were drawn by combining their expression values and clinical information (Fig. 5). Based on the significance threshold of p value <0.05, the expressions of 17 genes (except for ATP1A2, CNTN2, and MYH11) were significantly correlated with the survival time of the patients.

FIG. 5.

Gene expression related survival curve of 20 genes. The black curve and the grey curve represent the low expression group and the high expression group, respectively. HR, hazard ratio.

4. Discussion

In this study, a total of 634 upregulated and 769 and downregulated DEGs were identified between tumor and control groups. These DEGs were significantly involved in 15 upregulated pathways, such as p53 signaling pathway, and 16 downregulated pathways, such as neuroactive ligand–receptor interaction, and cell adhesion molecules (CAMs). In the coexpression network, PUS7, EFR3B, and NRCAM had the top three highest degrees. Additionally, 17 prognostic genes were selected, such as thrombospondin-1 (THBS1), caspase-8 (CASP8), glutamate ionotropic receptor AMPA type subunit 2 (GRIA2), GRIA4, and ADCYAP receptor type I (ADCYAP1R1).

TP53 is one of the most commonly mutated tumor suppressor genes in human cancers (Slattery et al., 2009; Yan et al., 2016). p53 is normally expressed at low levels to maintain the homoeostasis of cell cycle and cell death. The p53 protein is activated in response to stress factors, such as oxidative free radicals, hypoxia, DNA damage, or oncogenes. Activated p53 can either repress or transactivate the downstream target genes that in turn regulate apoptosis, cell cycle arrest, angiogenesis, and metastasis (Harris and Levine, 2005). In this study, the p53 signaling pathway was significantly upregulated, suggesting that it may promote the progression of glioma.

Moreover, many of the dysregulated genes in the p53 signaling pathway were found to be correlated with the survival time of the patients, such as THBS1 and CASP8. THBS1 is an extracellular matrix molecule involved in migration, invasion, and suppression of angiogenesis in glioma (Sargiannidou et al., 2001). Previous study has reported that THBS1 and transforming growth factor-beta 2 are correlated with glioma grade, showing higher expression levels in malignant high-grade gliomas, than in low-grade gliomas (Kawataki et al., 2010). CASP8 belongs to the cysteine aspartyl-specific protease family, an important defense system against tumor progression and malignant transformation (Bethke et al., 2008). It has been reported that advanced glioblastoma exhibits dysregulated CASP8 expression (Martinez et al., 2007). Taken together, THBS1 and CASP8 may be implicated in glioma development through p53 signaling pathway.

The pathway of neuroactive ligand–receptor interaction was the most significantly downregulated pathway. A recent study reveals that defective neuroactive ligand–receptor interaction pathway is a poor prognosticator in glioblastoma patients (Pal et al., 2017). GRIA2, GRIA4, and ADCYAP1R1 that associated with prognosis were enriched in this pathway. GRIA2 and GRIA4 encode AMPA-selective ionotropic glutamate receptor subunit. Glutamate is a major neurotransmitter in the central nervous system. There is evidence that glutamate and its receptors play roles in regulating tumor growth (Stepulak et al., 2005). Blocking expression of glutamate receptor subunits inhibits proliferation of cancer cells (Watanabe et al., 2008). ADCYAP1R1 (also named PAC1) is present in many neurologic cancers, such as neuroendocrine, brain, and pheochromocytomas, and play an important role in the modulation of growth and differentiation of these cancers (Reubi et al., 2000; Martinez and Vicente, 2013). Therefore, we speculated that pathway of neuroactive ligand–receptor interaction as well as its enriched genes, GRIA2, GRIA4, and ADCYAP1R1, may act as potential prognostic markers in glioma.

In the coexpression network, PUS7, EFR3B, and NRCAM had the top three highest degrees, besides NRCAM was involved in pathway of CAMs. PUS7 gene is reported to be a potential target in the immunomodulation of cancer (Qadir et al., 2018). The role of EFR3B has not been fully investigated to our knowledge. NRCAM was initially discovered in central and peripheral nervous system (Grumet et al., 1991). It was recently reported to be aberrantly expressed in human cancers (Chan et al., 2011; Tsourlakis et al., 2013). Overexpression of NRCAM promotes cell proliferation and motility. Previous study has reported that CAMs play an important role in the nervous system development (Sakurai, 2016). Taken together, considering the critical roles of the three genes in the coexpression network, we speculated that they may involve in glioma progression.

5. Conclusions

In conclusion, pathways of p53 signaling pathway and neuroactive ligand–receptor interaction may play important roles in glioma progression. PUS7, EFR3B, and NRCAM may be potential biomarkers of glioma. THBS1, CASP8, GRIA2, GRIA4, and ADCYAP1R1 may serve as prognostic markers in glioma.

Footnotes

Acknowledgment

This work was supported by the National Natural Science Foundation for Youth (Project No. 81400404).

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

References

Balss

, Meyer

, Mueller

, et al. 2008. Analysis of the IDH1 codon 132 mutation in brain tumors. Acta. Neuropathol. 116, 597–602.

Barrett

, Wilhite

S.E.

, Ledoux

, et al. 2013. NCBI GEO: Archive for functional genomics data sets—update. Nucleic Acids Res. 41, D991–D995.

Bethke

, Sullivan

, Webb

, et al. 2008. The common D302H variant of CASP8 is associated with risk of glioma. Cancer Epidemiol. Biomarkers Prevent. 17, 987–989.

Cancer Genome Atlas Research Network, Brat

D.J.

, Verhaak

R.G.

2015. Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N. Engl. J. Med. 372, 2481–2498.

CBTRUS. 2011. CBTRUS Statistical Report: Primary Brain and Central Nervous System Tumors Diagnosed in the United States in 2004–2007. Central Brain Tumor Registry of the United States, Hinsdale, IL.

Ceccarelli

, Barthel

F.P.

, Malta

T.M.

, et al. 2016. Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma. Cell. 164, 550–563.

Chan

J.Y.

, Ong

C.W.

, and Manuel

S.T.

2011. Overexpression of neurone glial-related cell adhesion molecule is an independent predictor of poor prognosis in advanced colorectal cancer. Cancer Sci. 102, 1855–1861.

Damian

, and Gorfine

2004. Statistical concerns about the GSEA procedure. Nat. Genet. 36, 663; author reply 663.

Ferluga

, and Debinski

2014. Ephs and Ephrins in malignant gliomas. Growth Factors. 32, 190–201.

10.

Goodenberger

M.L.

, and Jenkins

R.B.

2012. Genetics of adult glioma. Cancer Genet. 205, 613–621.

11.

Grumet

, Mauro

, Burgoon

M.P.

, et al. 1991. Structure of a new nervous system glycoprotein, Nr-CAM, and its relationship to subgroups of neural cell adhesion molecules. J. Cell Biol. 113, 1399–1412.

12.

Haeussler

, Zweig

A.S.

, Tyner

, et al. 2019. The UCSC Genome Browser database: 2019 update. Nucleic Acids Res. 47, D853–D858.

13.

Harris

S.L.

, and Levine

A.J.

2005. The p53 pathway: Positive and negative feedback loops. Oncogene. 24, 2899.

14.

Huang da

, Sherman

B.T.

, and Lempicki

R.A.

2009. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57.

15.

Jason

2013. Psych issues. JEMS. 38, 14.

16.

Kastin

2013. Handbook of Biologically Active Peptides, 2nd ed. Academic Press, Massachusetts, USA.

17.

Kawataki

, Naganuma

, Sasaki

, et al. 2010. Correlation of thrombospondin-1 and transforming growth factor-beta expression with malignancy of glioma. Neuropathology. 20, 161–169.

18.

Kleihues

, and Ohgaki

1999. Primary and secondary glioblastomas: From concept to clinical diagnosis. Neurooncology. 1, 44.

19.

Korshunov

, Ryzhova

, Hovestadt

, et al. 2015. Integrated analysis of pediatric glioblastoma reveals a subset of biologically favorable tumors with associated molecular prognostic markers. Acta Neuropathol. 129, 669–678.

20.

Martinez

, Setien

, Voelter

, et al. 2007. CpG island promoter hypermethylation of the pro-apoptotic gene caspase-8 is a common hallmark of relapsed glioblastoma multiforme. Carcinogenesis. 28, 1264–1268.

21.

Ohgaki

, and Kleihues

2007. Genetic pathways to primary and secondary glioblastoma. Am. J. Pathol. 170, 1445–1453.

22.

Pal

, Patil

, Kumar

, et al. 2017. Genetic landscape of glioma reveals defective neuroactive ligand receptor interaction pathway as a poor prognosticator in glioblastoma patients. AACR. 77(13 Suppl.), Abstract nr 2454.

23.

Qadir

, Aziz

M.A.

, Sari

C.P.

, et al. 2018. Transcriptome reprogramming by cancer exosomes: Identification of novel molecular targets in matrix and immune modulation. Mol. Cancer. 17, 97.

24.

Reifenberger

, Wirsching

H.G.

, Knobbethomsen

C.B.

, et al. 2016. Advances in the molecular genetics of gliomas—implications for classification and therapy. Nat. Rev. Clin. Oncol. 14, 434–452.

25.

Reubi

J.C.

, Läderach

, Waser

, et al. 2000. Vasoactive intestinal peptide/pituitary adenylate cyclase-activating peptide receptor subtypes in human tumors and their tissues of origin. Cancer Res. 60, 3105–3112.

26.

Riemenschneider

M.J.

, Jeuken

J.W.

, Wesseling

, et al. 2010. Molecular diagnostics of gliomas: State of the art. Acta Neuropathol. 120, 567–584.

27.

Ritchie

M.E.

, Phipson

, Wu

, et al. 2015. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47.

28.

Sakurai

2016. The role of cell adhesion molecules in brain wiring and neuropsychiatric disorders. Mol. Cell. Neurosci. 81, 4–11.

29.

Sargiannidou

, Zhou

, and Tuszynski

G.P.

2001. The role of thrombospondin-1 in tumor progression. Exp. Biol. Med. 226, 726–733.

30.

Shannon

, Markiel

, Ozier

, et al. 2003. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504.

31.

Slattery

M.L.

, Curtin

, Wolff

R.K.

, et al. 2009. A comparison of colon and rectal somatic DNA alterations. Dis. Colon Rectum. 52, 1304.

32.

Stepulak

, Sifringer

, Rzeski

, et al. 2005. NMDA antagonist inhibits the extracellular signal-regulated kinase pathway and suppresses cancer growth. Proc. Natl. Acad. Sci. U.S.A. 102, 15605–15610.

33.

Tsourlakis

M.C.

, Walter

, Quaas

, et al. 2013. High Nr-CAM expression is associated with favorable phenotype and late PSA recurrence in prostate cancer treated by prostatectomy. Prostate Cancer Prostatic Dis. 16, 159–164.

34.

Walter

, Sanchez-Cabo

, and Ricote

2015. GOplot: An R package for visually combining expression data with functional analysis. Bioinformatics. 31, 2912–2914.

35.

Watanabe

, Kanno

, Oshima

, et al. 2008. The NMDA receptor NR2A subunit regulates proliferation of MKN45 human gastric cancer cells. Biochem. Biophys. Res. Commun. 367, 487–490.

36.

Weller

, Wick

, Aldape

, et al. 2015. Glioma. Nat. Rev. Dis. Primers. 1, 15017.

37.

Yan

, Parsons

D.W.

, Jin

, et al. 2009. IDH1 and IDH2 mutations in gliomas. N. Engl. J. Med. 360, 765–773.

38.

Yan

, Wang

, Zhao

, et al. 2016. MicroRNA-34a affects chondrocyte apoptosis and proliferation by targeting the SIRT1/p53 signaling pathway during the pathogenesis of osteoarthritis. Int. J. Mol. Med. 38, 201–209.