Dissecting Prognosis Modules and Biomarkers in Glioblastoma Based on Weighted Gene Co-Expression Network Analysis

Introduction

Gliomas contain heterogeneous brain tumors which are traditionally classified based on distinct histological appearance and malignancy degrees.¹ In 2016 World Health Organization (WHO) classification, glioma have been assigned WHO grades I–IV where grade IV is regarded as glioblastoma.² Glioblastoma multiforme (GBM) is known as the most aggressive primary intracranial tumor, with short median survival time, excessive proliferation, and invasiveness. Hence, developing efficient prognosis signatures is crucial for glioblastoma patients. Past molecular analyses of glioblastoma reported many well-known molecular markers, such as p53, PI3K, RAF, isocitrate dehydrogenase 1 and 2 (IDH1/2), epidermal growth factor receptor (EGFR), X-linked alpha thalassemia mental retardation syndrome gene (ATRX), platelet-derived growth factor receptor alpha (PDGFRA), and insulin-like growth factor 1 (IGF-1).^3–7 Another prognosis marker, codeletion of chromosome arms 1p and 19q (1p/19q loss), displays a better prognosis.⁸ Additionally, O(6)-methylguanine-DNA methyltransferase (MGMT) is associated with the prognosis of patients and sensitivity of temozolomide therapy.⁹ Lack of NF1 could cause increased tumor-associated microglia infiltration.¹⁰

Benefited from progresses of high-throughput sequencing and bioinformatics, huge volumes of data make it easier to detect GBM-related hub genes, networks, and pathways. Nonetheless, identification of appropriate biological features remains challenging. Weighted gene co-expression network analysis (WGCNA) serves as a system biology and statistical tool to describe the correlation patterns among genes, as well as to recognize highly correlated modules.¹¹ In this study, we performed WGCNA for RNA-Seq data derived from the Chinese Glioma Genome Atlas (CGGA) and reconstructed gene co-expression networks. Then, we identified glioblastoma-associated gene modules and delineated several prognostic biomarkers related to the survival time of glioblastoma samples both in CGGA and The Cancer Genome Atlas (TCGA) datasets.

The independent prognostic model were identified by Log rank test along with Kaplan–Meier survival analysis. The results found that GRPR, CXCL5, and CXCL11 has been put forward as biomarkers for further guiding malignant glioma diagnosis and prognosis. Overall, this finding provides molecular considerations and points directions toward the molecular mechanism concerning gliomagenesis.

Materials and Methods

Data Acquisition

The gene expression data of 693 glioma samples was downloaded from the CGGA database (http://cgga.org.cn/). mRNA-sequencing procession methods were available for public in CGGA website (http://www.cgga.org.cn/about.jsp#mRNA_seq). In WGCNA, those genes without significant changes in expression between samples would be closely correlated. Next, the mRNA expression levels of 23,987 genes were calculated by variance analysis, and top 25% most variable genes (5,997) were subjected in the following steps. The clinical metadata were also derived and filtered for useful information in Table S1 (including CGGA_ID, progression status, histology, grade, gender, age, overall survival (OS), censor, and IDH mutation status). Among them, WHO grade was closely associated with OS.¹² The work flow diagram of this study was depicted in Figure S1.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was recruited via the R package ‘WGCNA’.¹¹ Firstly, RNA-Seq data were filtered to reduce outliers. In brief, soft threshold power of β (β=3) was established using scale-free topology criterion and the developed weighted adjacency matrix. Then, we chose the bottom-up algorithm, dynamic tree cut method, to recognize co-expression gene modules. To quantify module similarity and further identify biological interesting modules closely associated with tissues type, modules eigengenes (MEs) were estimated by the first principal component.¹³ Statistical significance was determined by a correlation test p-value <0.05. The module P-reservation function in the ‘WGCNA’ package was applied to calculate the Z-summary to estimate the conserved modules.

Comprehensive Functional Enrichment Analysis

Functional enrichment analyses including GO (Gene ontology), KEGG (Kyoto Encyclopedia of Genes and Genomes), ClinVar, Wiki pathway, and REACTOME of candidate genes in modules were carried out using the ClueGO plug-in of Cytoscape and p-value<0.05 (received by two-sided tests based on hypergeometric distribution) was the cut-off criterion.¹⁴ Enriched ontological terms and pathways were visualized by Cytoscape 3.8.0.

Integration of the PPI Network and Hub Genes Selection

The protein–protein interaction (PPI) network in gene module was developed using the STRING 11.0 database (http://string-db.org).¹⁵ P-value was received by a hypergeometric test. A combined score >0.4 was set in the PPI networks. A protein interaction relationship network table was visualized by Cytoscape 3.8.0 software (http://www.cytoscape.org/index.html).

The cytoHubba plug-in was employed to explore important nodes using several topological algorithms, such as Degree, Maximum Neighborhood Component (MNC) and centralities based on shortest paths, including bottleneck (BN), betweenness, and radiality.¹⁶ Proteins with high degree and betweenness might be important candidate genes which have key physiological regulatory functions.

Validation of Key Candidate Genes Expression Levels via RNA-Seq

The expression levels of core putative genes were next validated by CCGA online tools and the Gene Expression Profiling Interactive Analysis (GEPIA 2.0, http://gepia2.cancer-pku.cn/#index) which could analyze the RNA-Seq expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the Genotype-Tissue Expression (GTEx) projects.^17,18 Hub genes with |log₂Fold Change (FC)| > 1 and p-value < 0.05 (received by one-way ANOVA) were considered as statistically significant in GEPIA. p-Value was tested by t-test in CGGA tools.

Patients and Glioma Tissue Collection and qRT-PCR

The study protocol was approved by the Clinical Research Ethics Committee from the affiliated hospital of Zunyi medical college, which was conducted in accordance with the Declaration of Helsinki. (Zunyi, Guizhou, China). Written informed consent was obtained from all participants. GBM (n=25), LGG (n=25), and control brain tissue (n=25, sampled during surgical procedures for epilepsy) were obtained from the affiliated hospital of Zunyi medical college, immediately frozen in liquid nitrogen, and stored at −80°C until use.¹⁹ All tissues were subjected to quantitative reverse transcriptase PCR (qRT-PCR) to measure key candidate genes expression. The qRT-PCR reaction (10 µL) was formulated using the 2X SYBR Green qPCR Master Mix (US Everbright^® Inc., Suzhou, China). All qRT-PCRs were carried out on a CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The primer sequences designed by Primer 5.0 were shown in Table S2. The relative mRNA expression level was quantified using the 2^-ΔΔCt method.

Protein Expression Level and Immune Infiltrates of Candidate Genes

The protein expression levels of hub genes in normal and glioma tissues were investigated. The images were obtained from the Human Protein Atlas (HPA, http://www.proteinatlas.org) database. The purity-corrected Spearman’s rho of CXCL11, CXCL5, and GRPR expressions with the abundance of immune infiltrates, including dendritic cells, CD4+ T cells, neutrophils, macrophages, and NK cell were detected in GBM patients on TIMER (Tumor Immune Estimation Resource) online tools (https://cistrome.shinyapps.io/timer).

Survival Analysis of Key Genes

We assessed the association between hub genes and glioblastoma overall survival (OS) and disease free survival (DFS) in TCGA cohort by GEPIA 2.0.¹⁷ Survival data of the GBM subgroup were extracted based on Kaplan–Meier analysis of OS and DFS. Set the median risk score as cutoff, all patients in the TCGA set were classified into a high-risk group and a low-risk group using Log rank test. The prognostic value of screened hub genes in GBM was also analyzed using PrognoScan database (http://dna00.bio.kyutech.ac.jp/PrognoScan/index.html).²⁰

Statistical Analysis

The two-sided t-tests were used to perform differential analysis. All statistical tests and graphing were implemented via the R language, version 3.6.1. In the figures, statistical significance was indicated as follows: *, **, and *** represent p-value < 0.05, 0.01, and 0.001, respectively.

Results

Co-Expression Network Construction and Module Mining

WGCNA was utilized to frame a gene co-expression network to mine biologically meaningful gene modules. We selected top 25% (5,997/23,987) most variable genes which derived from 693 samples in CGGA datasets to construct a weighted gene co-expression network. Validation of co-expression modules might boost identification of hub genes that maintain important functions. WGCNA defined gene modules as a set of genes with topological overlaps. Based on scale-free topology criterion, the soft threshold power of β was set as 3 while scale-free topology model-fit R² was equal to 0.87 (Figure S2). After establishing a hierarchical clustering dendrogram, co-expression modules were represented by branches of different colors in the dendrogram (Figure 1A and B).

Figure 1 Weighted gene co-expression network in glioma. (A) Hierarchical clustering dendrogram of genes based on co-expression network analysis. Each piece of the leaves (short vertical lines) corresponds to individual gene. The branches correspond to modules of highly interconnected genes. Different colors below the dendrograms represent different gene modules. (B) Dendrogram of consensus module eigengenes. (C) The correlation coefficient heatmap of the co-expression module genes. Each bright spot corresponds to the correlation between each gene and other genes. The deeper the colors, the stronger is the connectivity between the two genes in the corresponding row and column. (D) Correlation between the gene module and clinical traits. The clinical traits include different WHO grades of glioma. The correlation coefficient in each cell represented the correlation between the gene module and the clinical traits, which decreased in size from red to blue. The corresponding P value is also annotated.

To determine module preservation, the Z-summary score was next calculated. Modules with a Z-summary score > 10 were considered as preserved.^21,22 Ultimately, 10 gene modules were distinguished with module size ranging from 54 genes in magenta module to 2,288 genes in turquoise module (Table S3). The genes which were not categorized into any gene modules finally were assigned to grey module (Figure 1C and Figure S3). With ranked in a descending order of Z-score, top five modules were turquoise (Z-score=39.57), brown (Z-score=37.68), blue (Z-score=37.10), yellow (Z-score=24.04), and grey (Z-score=20.28) modules. This observation indicated that these modules might provide better knowledge underlying the pathological mechanisms of GBM.

To exploit the relationship between gene modules and clinical indexes, the module eigengene (ME) was developed as the overall gene expression level of corresponding modules. We calculated correlations with ME and clinical phenotypes, such as primary and recurrent glioma in WHO II, WHO III, and WHO IV grades (Figure 1D). Genes in blue and yellow modules possessed similar expression patterns which were significantly accumulated in grade IV (GBM) and declined in other grades. Moreover, hierarchical clustering analysis of modules found that the blue and yellow modules were on the same branch (Figure 1B). The findings illustrated that the blue and yellow modules exhibited similar gene expression modes. Thus, we focused on analysis of the blue and yellow modules in the follow-up analyses.

Functional Enrichment Analysis of Gene in Modules

For a more in-depth understanding of genes in two candidate modules, GO functional enrichment, KEGG/wiki/REACTOME pathway enrichment, REACTOME reaction enrichment, and CLINVAR human diseases enrichment analyses were conducted using ClueGO plug-in of Cytoscape software (Table S4).^21,22 As shown in Figure 2A, genes in blue module were particularly enriched in neuroactive ligand-receptor interaction, intrinsic component of presynaptic membrane, taste transduction, and Golgi cis cisterna. Meanwhile, genes in yellow module were mainly involved in detection of chemical stimulus related to sensory perception of bitter taste, and hematopoietic stem cell differentiation (Figure 2B).

Figure 2 Functional Enrichment analysis of genes in blue (A) and yellow (B) modules. Different color displayed distinct GO terms. Statistical significance was derived by two-sided tests.

Hub Genes and Pathway Determination with PPI and Modular Analyses

Protein–protein interaction relationships of proteins in two modules were established using STRING. Finally, 164, 73 nodes and 79, 32 edges in blue (PPI enrichment p-value=6.62e-14) and yellow (PPI enrichment p-value=6.01e-08) modules were obtained, separately (Figure 3A and B). The degree of connectivity of genes was calculated using the cytoHubba plug-in of Cytoscape software. Genes with top 10 degree of connectivity were displayed in Figure 3C and D. Considering the degree and betweenness, HTR1B, P2RY4, GRM5, GRPR, and NMBR in blue module network and CXCL11, CXCL5, and IL5 in yellow module network were selected as hub genes (Table S5). Those candidate genes possibly have an important impact in glioma progression and can be regarded as diagnostic biomarkers.

Figure 3 Module visualisation and Identification of hub genes. Visualisation of the gene co-expression network of the blue (A) and yellow modules (B) were generated using Cytoscape. Top 10 hub genes in blue (C) and yellow (D) modules were calculated by cytoHubba ranked by degree. The darkness of red in nodes represents the depth of degree.

Transcriptional Expression Patterns of Key Candidate Genes in TCGA and GTEx

To investigate if the eight hub genes were potential prognostic biomarkers of glioma, we first compared these mRNA expression levels in GTEx dataset containing 207 normal tissue and in TCGA cohort which comprised 163 GBM and 518 LGG specimens. Consequently, we found that GRPR, CXCL5, and CXCL11 mRNA were more abundant in glioma than in control specimens both in GBM and LGG (Figure S4A-C). Whereas, mRNA expression of other five genes did not show significant difference even enriched highly in control compared to glioma samples. After checking the expression patterns about different progression status of glioma of these three central genes in CGGA database, the findings implicated that CXCL11 expression was relatively higher in recurrent than primary type significantly, especially in glioma and WHO II grade (Figure 4A and B and S4G).

Figure 4 Expressions patterns about different grades and progression status of CXCL5 (A) and CXCL11 (B) in CGGA. Expression levels of CXCL5 (A) and CXCL11 (B) among normal, GBM and LGG samples. as well as in different molecular subtypes of GBM. P-value was tested by t-test. Immunohistochemical analyses of CXCL5 (C) and CXCL11 (D) protein expression levels in normal brain tissues (cerebral cortex) and glioma cancer specimens. Images were obtained from the HPA resource.

As GRPR, CXCL5, and CXCL11 were mainly up-regulated in GBM, we explored these expression across distinct subtypes according to Verhaak’s classification by GEPIA 2.0.²³ We compared the mRNA levels between four molecular stratifications, namely classical (n=40), mesenchymal (n=55), neural (n=28), and proneural (n=37). GRPR and CXCL11 mRNA were enriched in classical glioblastoma in comparison with other subtypes (Figure S4D and S4F). Meanwhile, mesenchymal glioblastoma highly expressed CXCL5 than other subtypes (Figure S4E). Immunohistochemical analysis was also investigated and the results revealed that the protein expression level of CXCL11 was accumulated in glioma tumor tissues compared to non-tumor tissues. Remarkably, we noted that protein expression of CXCL5 showed the opposite trend which need more investigation. Immunohistochemistry Images were captured from the Human Protein Atlas (http://www.proteinatlas.org) online database (Figure 4C and D).

Survival Prediction of GBM Patients

To get an indication of the relationship between higher GRPR, CXCL5, and CXCL11 mRNA abundance and survival situation of glioma patients, we analyzed Kaplan–Meier survival curves of the all glioma cohorts in TCGA. The results suggested that a high level of GRPR mRNA was associated with clearly shorter overall survival (OS) as well as disease free survival (DFS) (Figure 5A). Similarly, high expressions of CXCL5 and CXCL11 were linked with a worse prognosis of glioma patients (Figure 5B and C).

Figure 5 Survival curves between groups with low and high GRPR (A), CXCL5 (B), and CXCL11 (C) expressions in TCGA. The red and green lines represent cases with high and low expressions, separately. The X axis indicates survival time (months) and the Y axis indicates the present survival (%).

To further confirm the prognostic roles of GRPR, CXCL5, and CXCL11, two GEO (Gene Expression Omnibus) datasets (GSE4412 and GSE7696) of GBM and glioma were evaluated by PrognoScan database.²⁰ As shown in Figure S5, the higher expression of these hub genes might figure a worse OS in GBM and glioma, especially for CXCL5 and CXCL11 whose showed statistical significance. These findings offered evidence for the prognostic power of the candidate implications GRPR, CXCL5, and CXCL11.

Association of Immune Infiltration Level and Immune Checkpoints

Tumor infiltrating lymphocytes (TILs) have been found to be an independent predictor of survival in cancers.²⁴ Thus, we researched whether GRPR, CXCL5, and CXCL11 expression was correlated with immune infiltration levels in GBM by TIMER tools. The findings implied that CXCL5 and CXCL11 had marginal positively associations with infiltrating levels of neutrophil, B and dendritic cell, while negatively correlated with CD8⁺ T cell and CD4⁺ T cell in GBM. GRPR expression had relatively negative relationships with these immune cells except B cell and CD8⁺ T cell (Figure 6A–C). Then, Pearson correlation analysis was performed between three candidates and several immune checkpoints, which have been developed as therapeutic targets in clinical characterization.²⁵ As shown in Figure 6D, the outcomes revealed that CXCL11, CXCL5, and GRPR was positively associated with 20 checkpoint members, suggesting possible synergistic effects of the hub genes with immune checkpoints, especially with CD28, CD80, IDO1, and PDCD1LG2. Notably, co-expression of PD-1 with CXCL5 was consistent with the finding in past study.²⁶

Figure 6 Correlation of CXCL5 (A), CXCL11 (B), and GRPR (C) expressions with immune infiltration level in GBM. (D) Corrplot map of three DEGs and immune checkpoint markers in TCGA-GBM database.

Verification of Functional Roles Using qRT-PCR

To gain better insight on mRNA expression levels of three core DEGs in vivo, qRT-PCR validation was performed based on our study cohort patients. The Clinical data of 75 patients were tabulated in Table S6. The qRT-PCR results revealed that the relative expression of GRPR, CXCL5, and CXCL11 were distinctly increased in the glioma group as compared to the control treatment (Figure 7A–C). Moreover, the same trends were observed between LGG and control specimens (Figure 7D–F). To sum up, the qRT-PCR results were consistent with our integrative meta-analysis, suggesting the critical role of these core DGEs might play necessary roles in the mechanism of glioma, especially in GBM.

Figure 7 GRPR (A), CXCL5 (B), and CXCL11 (C) expression levels correlated in glioma tissue. qRT-PCR analysis of GRPR (A and D), CXCL5 (B and E), and CXCL11 (C and F) in GBM (n=25), LGG (n=25), and control brain tissues (n=25). P-value was tested by t-test.

Discussion

Glioblastoma multiforme (GBM) is a common aggressive form of primary brain tumor, characterized by rapid cell proliferation, resistance to apoptosis, and extensive invasion.²⁷ Due to resistance to currently available therapeutic regimens, the prognosis for patients with GBM is remarkably poor. As a key factor of glioma classification, glioma grade is involved in metastasis and tumor relapse while usually estimated by tumor imaging and biopsy.²⁸ To be a supplement to existing techniques, molecular biomarker analysis has a wide application perspective. Although there are several published studies using WGCNA on glioma, most of them has focused on the relationship between module and traits like age, gender, overall survival, IDH mutation, and so on.^21,29–31 In this case, our study paid more attention to identify genes which are clustered with primary and recurrent types of four WHO grade I to IV. Similarly, Mukherjee’s work has analyzed relationship between module and WHO grade II to IV.³² The enriched GO terms in above previous research contain negative regulation of response to cytokine stimulus, fall in line with our findings about chemokines.

Chemokines, a multifunctional type of small-conserved proteins, are reported to be essential for many biological processes, especially tumor cell proliferation, migration and invasion.^33,34 Chemokines and its receptors have been found to dictate cancer development and progression through epithelial-mesenchymal transition.³⁵ Chemokines act as principal components which express at tumor site and interact with immune cells.³⁶ For sample, CXCL11 is elevated in diverse solid tumors mediating tumor progression, metastasis, and lymphocyte infiltration.³⁷ Inhibition of CXCL11 by neutralizing antibodies might prominently decrease the degree of morphine tolerance and attenuate the activation of astrocytes.³⁸ As for another important chemokine CXCL5, overexpressed CXCL5 is found in various cancers and also facilitates tumor progression and metastasis.^39,40 CXCL5 promotes the metastatic potential of cells by upregulation of Snail in breast cancer.⁴¹ In nasopharyngeal carcinoma, CXCL5 induces epithelial-mesenchymal transition (EMT) and activates neutrophils to promotes nasopharyngeal carcinoma, lung, and gastric cancer metastases.^42,43 Furthermore, CXCL5 from cancer-related fibroblasts enhances PD-L1 expression in colorectal carcinoma and melanoma cells.²⁶ Of note, recent studies demonstrate that CXCL chemokines play important roles in glioblastoma biology, involved in cell proliferation, death and migration.⁴⁴ CXCL5 stimulates glioma cells proliferation and migration in autocrine- and paracrine-dependent ways.⁴⁵

Besides, gastrin-releasing peptide receptor (GRPR) acts as proliferative and inflammatory factors in living systems.⁴⁶ The activation of GRPR promotes cell proliferation, whereas GRPR antagonists induce antiproliferative effects in both in vitro and in vivo experimental GBM models. In addition, according to the knockdown of GRPR induces senescence in GBM cells, GRPR have the potential to serve as a therapeutic target in brain tumors.⁴⁷

These studies implicate that GRPR, CXCL11, and CXCL5 play diverse roles in oncoma development and progression. As mentioned, GRPR, CXCL11, and CXCL5 are the centermost upregulated gene in our network analysis. Expression of GRPR, CXCL11, and CXCL5 were higher in glioma patients than health control. Thus, our findings that GRPR, CXCL11, and CXCL5 are tightly linked to glioma is consistent with above earlier reports, suggesting our analysis method is reliable and practical. Moreover, expression of CXCL11 involves in mitogen-activated protein kinase (MAPK) pathway,⁴⁸ known as hyperactivated by genetic alteration in GBM. We consequently investigated the Pearson correlation between the expression of three candidate genes and that of Top 10 altered genes in GBM (Figure S6). The significant associations were observed between CXCL11 with EGFR, CXCL5 with ATRX and PTEN, as well as GRPR with FLG, NF1, TP53, and TTN. The observations implied high expression of these DEGs might mirror the existence of the oncogenic mutations frequently existed in GBM. We next tested whether these DEGs have an intrinsic relationship in glioma by analyzing the STRING database. We found that there is one significant connection between CXCL11 and CXCL5. However, there are some limitations need to be addressed in this study. Since the glioma patient samples are usually hard to get, we just checked a small cohort of 25 GBM patients in this study. Functional verification of the selected genes in large-scale patient samples would support our results better. Moreover, the knockout mouse model of each candidate DEG will facilitate a clear understanding of its role in glioma pathological mechanisms.