Identification of biomarkers for the transition from low-grade glioma to secondary glioblastoma by an integrated bioinformatic analysis

Abstract

Secondary glioblastoma (sGBM) is a type of glioblastoma multiforme that evolves from low-grade glioma (LGG). However, the mechanism of this transition still remains poorly understood. In this study, we used weighted gene co-expression network analysis (WGCNA) on the gene expression profiles of glioma samples from the Chinese Glioma Genome Atlas (CGGA) database to identify key genetic module related to distinguish histological characteristics. Here, the brown module was highly correlated with histological characteristics and was selected as the hub module. By applying functional annotation analysis, we found that biological processes related to the cell-cycle and DNA-replication were enriched in the genes of the brown module. After constructing a protein-protein interaction (PPI) network, validation of differential gene expression, and survival analyses, we ultimately identified five hub genes: CCNB2 (Cyclin B2), KIF2C (Kinesin Family Member 2C), CDC20 (Cell Division Cycle 20), TPX2 (TPX2 Microtubule Nucleation Factor), and PLK1 (Polo Like Kinase 1). In addition, a computational risk model was developed for predicting the clinical outcomes of sGBM patients by combining gene expression levels. This gene signature was demonstrated to be an independent predictor of survival by univariate and multivariable Cox regression analysis. Finally, we used the Genomics of Drug Sensitivity in Cancer (GDSC) database to predict the responses of sGBM patients to routine chemotherapeutic drugs. Patients from the high-risk group were more sensitive to common chemotherapies during clinical treatment. Our findings based on comprehensive analyses might advance the understanding of sGBM transition and aid the development of novel biomarkers for diagnosing and predicting the survival of sGBM patients.

Keywords: Secondary glioblastoma; bioinformatics analysis; molecular mechanism; transition; weighted gene coexpression network analysis (WGCNA).

Figure 1

Clustering dendrogram of samples and selection of soft-threshold power. A. The clustering was based on the expression data from the CGGA RNA-sequencing dataset. The top 5000 genes with the highest MAD values were used for WGCNA analysis. The displayed colors correspond to the histological characteristics of samples. B. Analysis of the scale-free fit index for different soft-thresholding powers. C. Analysis of the mean connectivity for different soft-threshold powers. D. Histogram of connectivity distribution when β = 11. E. Linear model fitting of the R² index showing a good quality of fit (R² = -0.94).

Figure 2

Identification of modules associated with the histopathological features of the samples. A. Clustering tree (dendrogram) of genes based on co-expression network analysis. Genes were clustered based on dissimilarity measure (1-TOM). Bars below correspond to modules of genes with high interconnectivity. B. Heatmap of correlations between the modules’ eigengenes and histological characteristics of the samples. Each row corresponds to a specific module color. The upper number in each cell is the correlation coefficient of each module with histology, and the lower number is the p-value. Color is coded according to the correlation coefficient. C. Distribution of average gene significance and errors in the modules associated with the histological characteristics of the samples.

Figure 3

GO functional and KEGG pathway enrichment analysis of genes in the brown module. A. Top 10 significantly enriched biological process annotations. B. Top 10 significantly enriched cellular component annotations. C. Top 10 significantly enriched molecular function annotations. D. Top 10 significantly enriched KEGG pathways. The x-axis represents the number of genes in the corresponding gene term and the y-axis shows the gene terms. The adjusted p-value of each term is colored according to the legend.

Figure 4

Detection and validation of candidate hub genes. A. PPI network of 22 genes with high connectivity and top-20 MCC genes in the brown module. Nodes colored in yellow are candidate genes identified both in the co-expression network and the PPI network. B. Boxplots of the expression levels of candidate genes in LGG and sGBM samples in the CGGA RNA-sequencing and microarray datasets. **; P < 0.01, ***; P < 0.001. Two-tailed Student’s t-test was used to evaluate the statistical significance of differences. C. ROC curves measuring the predictive value of each candidate gene in the CGGA database. The X-axis shows the false-positive rate, shown as “1-Specifcity”. The Y-axis indicates the true positive rate, shown as “Sensitivity”.

Figure 5

Overall survival analyses based on seven candidate hub genes using the log-rank test. Candidate genes in the brown module and significant results of survival analysis (P < 0.05 was defined as statistically significant). These were CCNB2 (A), CDC20 (B), PLK1 (C), KIF2C (D) and TPX2 (E), respectively. HR: hazard ratio, CI: confidence interval.

Figure 6

Protein expression levels of five hub genes in LGG and sGBM tissues. A. Representative photographs of IHC staining for CCNB2, CDC20, PLK1, KIF2C, and TPX2 in clinical human samples of LGG and sGBM. Magnification, × 200. Scale bar = 50 μm. B. Quantification of IHC staining intensity for each protein in the specimens. The integrated optical density (IOD) and the area were quantified using Image-Pro Plus software. Significance tested by Student’s t-test, *; P < 0.05, **; P < 0.01, ***; P < 0.001.

Figure 7

Construction of the risk signature for sGBM based on five hub genes. The patients were classified into high- and low-risk groups based on the median value of the risk scores. A. Risk score distribution, survival status of sGBM patients and expression heatmap of the five hub genes. Red indicates a high expression level of a given gene, whereas blue indicates a low expression level. B. Kaplan-Meier survival curves for overall survival in the high- and low-risk group. C. Time-dependent ROC curves for predicting one-year survival of sGBM patients in the RNA-sequencing dataset based on the signature and IDH1 status. The AUC of the gene signature was 0.85, and the AUC of the IDH1 mutation was 0.55.

Figure 8

Differentially regulated pathways and predicted responses to chemotherapy in the high- and low-risk group. All transcripts were ranked according to the log2 (fold change) value derived from differential gene expression analysis between the two groups. Various pathways enriched in the high- (A) and low-risk group (B) were plotted. The NES and FDR value of each term were shown. (C) Boxplots of the estimated IC50 values of TMZ and cisplatin for tumor cells from two groups. Wilcoxon test (Mann-Whitney test) was used for comparison across groups.