Identification of key modules and hub genes in glioblastoma multiforme based on co-expression network analysis

Abstract

Glioblastoma multiforme (GBM) is the most malignant primary tumour in the central nervous system, but the molecular mechanisms underlying its pathogenesis remain unclear. In this study, data set GSE50161 was used to construct a co-expression network for weighted gene co-expression network analysis. Two modules (dubbed brown and turquoise) were found to have the strongest correlation with GBM. Functional enrichment analysis indicated that the brown module was involved in the cell cycle, DNA replication, and pyrimidine metabolism. The turquoise module was primarily related to circadian rhythm entrainment, glutamatergic synapses, and axonal guidance. Hub genes were screened by survival analysis using The Cancer Genome Atlas and Human Protein Atlas databases and further tested using the GSE4290 and Gene Expression Profiling Interactive Analysis databases. The eight hub genes (NUSAP1, SHCBP1, KNL1, SULT4A1, SLC12A5, NUF2, NAPB, and GARNL3) were verified at both the transcriptional and translational levels, and these gene expression levels were significant based on the World Health Organization classification system. These hub genes may be potential biomarkers and therapeutic targets for the accurate diagnosis and management of GBM.

Keywords: TCGA; WGCNA; biomarkers; glioblastoma multiforme; survival.

Fig. 1

Clustering of samples and determination of soft‐thresholding power. (A) Clustering of samples. (B) The scale‐free topology fitting index (R ², y‐axis) as a function of the soft‐thresholding power (x‐axis). The red line indicates R ² = 0.9. (C) The mean connectivity (degree, y‐axis) is displayed as a function of the soft‐thresholding power (x‐axis). Red Arabic numbers in the panels denote different soft‐thresholds. β = 8, there was a trade‐off between maximising the scale‐free topology model fitting index (R ²) and maintaining a high mean number of connections.

Fig. 2

Cluster of module eigengenes and cluster dendrogram. (A) Cluster dendrogram of module eigengenes. (B) The cluster dendrogram of genes in GSE50161. Each branch of the dendrogram represents a gene. Each designated colour represents a co‐expression module.

Fig. 3

Identification of key modules. (A) Network heatmap plot of co‐expression genes. Different colours of the horizontal axis and the vertical axis represent different modules. The progressively saturated yellow and red colours indicate the high co‐expressed interrelation in the heatmap. (B) Heatmap of the correlation between module eigengenes and GBM. Each row represents a module eigengene and each column represents trait. The table is coloured by correlation according to the colour legend. The turquoise module is the most negatively correlated with GBM, and the brown module is the most positively correlated with GBM. (C) The hierarchical clustering dendrogram and the adjacency heatmap of each module. The top is the hierarchical clustering of the module crucial genes (labelled by colours). The bottom is the adjacency heatmap, where each column and row correspond to one module crucial gene (labelled by colour) or trait labelled by y. Red represents high adjacency (positive correlation), while blue colour represents low adjacency (negative correlation). Squares of red colour along the diagonal are the meta‐module.

Fig. 4

Module eigengenes in key modules. (A) Scatter plot of module eigengenes in the brown module. (B) Scatter plot of module eigengenes in the turquoise module.

Fig. 5

GO enrichment analysis of key modules. (A) GO enrichment analysis of genes in the brown module (top 10 in BP, CC, and MF are listed). The y‐axis depicts names of terms in BP, CC, and MF, respectively, and the x‐axis depicts gene ratio in the module. The circle size represents the count and colours represent the P‐value; (B) GO enrichment analysis of turquoise module genes (top 10 in BP, CC and MF are listed). CC, cellular component; MF, molecular function.

Fig. 6

KEGG pathway enrichment and top 30 genes of key modules. (A) The KEGG pathways of the brown module. The y‐axis depicts the names of the terms of the pathways, and the x‐axis represents the count. Colours represent P‐values. (B) The KEGG pathways of the turquoise module. (C) In the brown module, genes with node degrees in the top 30 are displayed. (D) In the turquoise module, genes with node degrees in the top 30 are displayed. The higher the rank of the genes, the deeper the colour of the genes.

Fig. 7

Survival analysis of hub genes. (A–D) In the brown module, genes with top 30 node degrees and significant results of survival analysis in HPA are SHCBP1, NUSAP1, NUF2 and KNL1 (P < 0.05 was regarded as significant). (E–H) In the turquoise module, genes with top 30 node degrees and significant results of survival analysis in TCGA are SULT4A1, SLC12A5, NAPB and GARNL3 (P < 0.05 was regarded as significant).

Fig. 8

Screening and validation of hub genes at the transcriptional level. (A) Screening hub genes in GSE50161 (including 34 paediatric GBM samples and 13 normal samples). The expression status of four hub genes (NUSAP1, NUF2, SHCBP1 and KNL1) was positively correlated with disease status. This was also consistent with the results for the brown module of WGCNA. The expression status of four hub genes (SULT4A1, SLC12A5, NAPB and GARNL3) was negatively correlated with disease status. This was consistent with the results in the turquoise module of WGCNA. (B) Validation of eight hub genes in GSE4290 (including 81 paediatric GBM samples and 23 normal samples), and the results were the same as earlier. (C) The differential expression levels of hub genes were demonstrated in TCGA by GEPIA. These results were consistent with the above results. These results fully demonstrated the reliability of our findings. (t‐test; *P < 0.05; **P < 0.01; ***P < 0.001).

Fig. 9

Validation of hub genes at the translational level by The HPA database. The translational expression levels of NUSAP1, SHCBP1 and KNL1 in GBM were higher than those in normal tissue. The translational expression levels of SULT4A1, SLC12A5, NAPB and GARNL3 in GBM were lower than those in normal tissue.

Fig. 10

Hub gene expression distribution in the World Health Organization classification (188 GBM samples in WHO grade II, 255 samples in grade III and 249 samples in grade IV). The expression levels of NUSAP1, SHCBP1, KNL1 and NUF2 genes of the brown module and SULT4A1, SLC12A5, NAPB and GARNL3 genes of the turquoise module were significantly different in the different World Health Organization classification. (t‐test; *P < 0.05; **P < 0.01; ***P < 0.001).

Fig. 11

Knockout NUSAP1, SHCBP1, NUF2 and KNL1 gene can decrease the proliferation and cloning of U87 cells. (A) The level of NUSAP1, SHCBP1, NUF2 and KNL1 mRNA was determined by RT‐PCR after siRNA transfection. (B) The proliferation ability of U87 cells was evaluated by CCK8 method. (C) Clone formation analysis was used to evaluate the cloning ability of U87 cells (unpaired t‐test, N = 3, mean ± SD, *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001).