Exploring the potential biomarkers for prognosis of glioblastoma via weighted gene co-expression network analysis

Data information and construction of WGCNA

The Series Matrix file data File of GSE145128 was downloaded from the NCBI GEO public database, which were contained 15 GBM patients and sets of transcriptional data, including untreated group (n = 7) and recurrent group (n = 8), for the construction of WGCNA co-expression network.

We constructed a weighted gene co-expression network to find co-expressed gene modules, and clarified the relationship between the gene network and phenotype and hub genes. The WGCNA-R package was used to construct a co-expression network of genes in the GSE145128 dataset, where the soft-thresholding power was set to 16. The weighted adjacency matrix is converted into a topological overlap matrix (TOM) to estimate the network connectivity, and a hierarchical clustering method is used to construct a clustering tree structure of the TOM matrix. Different gene modules are represented by different cluster tree branches and colors. All genes are divided into multiple modules through gene expression patterns, and genes with similar expression patterns are divided into one module based on weighted correlation coefficients and expression patterns.

Enrichment analysis of gene module function

In order to obtain the biological functions and signaling pathways involved in the interest module of WGCNA, the Metascape database (http://www.metascape.org) was used for annotation and visualization, and Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on the genes in the specific module. Min overlap ≥3 & p ≤ 0.01 was considered statistically significant.

Identifying key genes

To determine key genes, the most important thing is whether they have an impact on tumor prognosis. According to WGCNA theory, key genes have the highest connectivity in the module, which determines the biological significance of the module (Chen et al., 2019). So we think that key genes must exist in the interested module of WGCNA. Combined with the above two points, we analyzed the Kaplan–Meier survival rate of all genes in the interest module of WGCNA. We believe that the genes that can affect the prognosis of GBM patients are the key genes. And then,the next step is to explore and verify the specific molecular mechanism of key genes.

Download and pre-processing data from TCGA

TCGA database as the biggest cancer gene information database, including gene expression data, the miRNA expression data and copy number variation, DNA methylation, SNPS and other data. We downloaded the processed original mRNA expression data of GBM. A total of 159 specimens were collected (Table S1).

Co-expression analysis

The co-expression of the key genes were analyzed. The correlation coefficient filter condition was 0.3 and the p-value was 0.001. After screening the genes with the most significant expression of key genes, the correlation analysis circles of key genes were plotted using “corrplot” and “circlize” packages.

GSCALite and GDSC

GSCALite is a set analysis platform for cancer genes. It integrates cancer genomics data of 33 cancer types from TCGA, drug response data from GDSC and CTRP, and normal tissue data from GTEx, and conducts gene set analysis in the data analysis process. Our study through this analysis was carried out on the key genes.Then base on the largest publicly available pharmacogenomics database (GDSC, the Genomics of Drug Sensitivity in Cancer, https://www.cancerrxgene.org/), we used the R packet “pRophetic” to predict the chemosensitivity of each tumor sample, and the estimated IC_ (50) of each specific chemotherapeutic drug treatment was obtained by ridge regression. The prediction accuracy was measured by 10-fold cross-validation with the GDSC training set. Select default values for all parameters, including “combats” to remove batch effects, “allSoldTumours” for tissue types, and average values for summarizing repetitive gene expressions (Liu et al., 2018).

TIMER database analysis

TIMER is a website for systematically testing the molecular characteristics of tumor-immune interactions (Li et al., 2017b). This website has incorporated 10,897 samples ranging from 32 different kinds of cancer types from the TCGA dataset (Shi et al., 2020). In this study, TIMER was used to explore the relationship between key genes and the contents of immune cells and to compare the infiltration levels between tumors with different somatic copy number changes of key genes.

Gene functional analysis

GSVA uses a non-parametric and unsupervised method, and bypasses the traditional method of explicitly modeling phenotypes in affluent scoring algorithms (Hanzelmann, Castelo & Guinney, 2013). By comprehensively scoring the gene set of interest, GSVA converts the gene level change into the pathway level change, and then judges the biological function of the sample. In this study, the gene sets were downloaded from the Molecular signatures database (v7.0 version),and used the algorithm of GSVA comprehensive score of each gene set, evaluating the potential biological function change different samples.

Statistical analysis

All statistical analyses were performed in R language (version 3.6). All the statistical tests were bilateral, and p < 0.05 was statistically significant.

Identification of gene co-expression modules

The Series Matrix File data File of GSE145128 was downloaded from the NCBI GEO public database. A total of 15 transcriptional data sets, including untreated group (n = 7) and recurrent group (n = 8), were used to construct the WGCNA co-expression network. In order to determine whether the 15 samples in GSE145128 were suitable for network analysis, a sample dendrogram and similar clinical features were studied. We confirmed that all samples were included in the group (Fig. 1A). The soft-thresholding power was set as 16 for the subsequent construction of co-expression network (Figs. 1B, 1C). The clustering tree structure of TOM matrix was constructed by hierarchical clustering method. The different branches and colors represent different gene modules (Fig. 1D). The network heatmap was used to analyze the interaction of 41 modules (Fig. 1E). The results showed that each module was independent of each other, indicating that each module was highly personalized and the gene expression of each module was relatively independent.

Correlation of modules and clinical traits

In order to study the relationship between these modules and the prognosis of GBM, we investigated the correlation between each module and the prognosis of GBM.We found that the darkturquoise module had the highest correlation with disease relapse (r =  − 0.62, p = 0.01); (Figs. 2A, 2B). We used Metascape to analyze the function and pathway of the darkturquoise module. Metascape can identify the enrichment process in the gene list and the association between enrichment processes (Tripathi et al., 2015a; Shannon et al., 2003) by querying many databases,such as GO functional, Hallmark Gene Sets, and KEGG pathways (Tripathi et al., 2015; Zheng et al., 2018). Based on GO enrichment analysis, it was found that the co-expressed genes within the modules of interest s mainly related to the steady state of cellular transition metal ion homeostasis, protein hydroxylation, intrinsic apoptotic signaling pathway, ncRNA metabolic process (Fig. 2C). The KEGG pathway analysis revealed that the co-expressed genes within the modules of interest was mostly enriched in the ‘Ferroptosis’ (Fig. 2C). In addition, the enrichment processes were highly connected and could be clustered into a complete network (Fig. 2C). These results indicated that these functions were related in the occurrence and development of GBM.

Identification of key genes in darkturquoise module

According to the WGCNA theory, key genes have the highest connectivity in the module, which determine the biological significance of the module and therefore influence the survial of patients intensively (Chen et al., 2019). Therefore, we searched for hub genes in the darkturquoise module. We analyzed the Kaplan–Meier survival curves of 139 genes in the darkturquoise module. Fourteen genes with significant survival analysis results (p < 0.05) were selected for sequencing (Table 1). Further determination of the key gene requires a combination of its expression, typicality in previous studies, and previous research experience in our lab. Among the 14 candidate genes we selected, most genes have been confirmed to be related to the prognosis of GBM, such as: SPAG4, FKBP1B, DLEU1, PRKAR2B, NRL, CD24, GAS6; some genes have not been clearly studied to be related to the occurrence of any known tumors, such as CAMSAP2; the remaining genes have not been significantly expressed in GBM, such as CORO6 (Lo et al., 2009; Sun et al., 2021; Wang et al., 2019; Zhao et al., 2019). Finally, we found that four genes (DARS/GDI2/P4HA2/TRUB1) had an impact on the survival rate of patients with GBM, and were also confirmed to be related to tumorigenesis. At the same time, they were not confirmed to participate in the occurrence and development of GBM, which met the conditions for our further research (Fig. 3). Survival analysis showed that the patients with low expression of DARS/GDI2/TRUB1 and high expression of P4HA2 had poor prognosis.

Analysis of the co-expression of key genes

It was clear that the key genes can affect the process of disease progression by regulating related genes. It can be assumed that DARS/GDI2/P4HA2/TRUB1 was associated with the most abundant pathways and genes and could regulate more biological processes. In order to assess the gene correlation of DARS/GDI2/P4HA2/TRUB1, we analyzed the co-expression of DARS/GDI2/P4HA2/TRUB1 through Pearson correlation analysis(cor>0.3, p < 0.001). We screened the 10 genes with the strongest correlation with the expression of DARS/GDI2/P4HA2/TRUB1, drew the correlation analysis map and heat map of DARS/GDI2/P4HA2/TRUB1 (Figs. 4A–4H), and found that the correlation between DARS and CDC73 was the highest, and the correlation between GDI2 and CDC123 was the highest. P4HA2 and B4GALT1 have the highest correlation, TRUB1 and CUL2 have the highest correlation (Figs. 4I–4L). Then,we verified the modules of these four genes in WGCNA. We found that CDC73, CDC123 and CUL2 all exist in the darkturquoise module, which is consistent with DARS/GDI2/P4HA2/TRUB1 (Table S2). Although B4GALT1 does not exist in the darkturquoise module, studies have confirmed that B4GALT1 can affect the development of GBM by regulating the apoptosis and autophagy (Wang, Li & Xie, 2020a). Among them, CDC73 and CDC123 are cyclins of cell division (Sun et al., 2017), B4GALT1 is one of seven β - 1, 4-galactosyltransferases (B4GALT). CUL2 contributes to form E3 ubiquitin ligase that can recognize numerous substrates and is involved in a variety of cellular processes (Liu et al., 2019). These four genes have been shown to have a close relationship with many kinds of cancers (Cao et al., 2020a; Dou et al., 2020; Li et al., 2019a), such as thyroid carcinoma (Sarquis et al., 2019), breast cancer, etc.

Cancer-related pathways and drug sensitivity analysis of key genes

First, we investigated the role of key genes in all well-known cancer-related pathways, as the following: TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, Hormone ER, Hormone AR, EMT, DNA Damage Response, Cell Cycle, Apoptosis pathways. The results found that DARS participated in the TSC/mTOR pathway activation; GDI2 was involved in Apoptosis, TSC/mTOR pathway activation; P4HA2 was involved in DNA Damage Response, EMT, Hormone AR, Hormone ER, RAS/MAPK and TSC/mTOR pathway; TRUB1 was involved in Apoptosis, DNA Damage Response, EMT, Hormone AR, PI3K/AKT and TSC/mTOR pathway (Fig. 5A). To investigate whether the expression of DARS/GDI2/P4HA2/TRUB1 in GBM had an impact on treatment (e.g., chemotherapies), we constructed a predictive model on six commonly used chemo drugs (i.e., AKT.inhibitor, Cisplatin, Dasatinib, Erlotinib, Gefitinib, and Gemcitabine) and confirmed that high expression of DARS was less sensitive to Cisplatin (p = 0.00026) and Gemcitabine (p = 0.0024), high expression of GDI2 was less sensitive to AKT.inhibitor (p = 5.2e −05), Cisplatin (p = 0.00067), Dasatinib (p = 0.012) and Gemcitabine (p = 3.5e −06), low expression of P4HA2 was less sensitive to Cisplatin (p = 0.0032), and Gemcitabine (p = 0.00017), and high expression of TRUB1 was less sensitive to AKT. inhibitor (p = 0.00043)(Fig. 5B).

Immune cells infiltration analysis

In view of the obvious prognostic value of DARS/GDI2/P4HA2/TRUB1 gene, we used the TIMER database to determine whether there was an association between tumor-infiltrating and immune cells and DARS/GDI2/P4HA2/TRUB1 expression.Results showed that DARS expression was positively correlated with tumor purity, P4HA2 and B cells (partial cor =−0.239, p = 7.89e−07), P4HA2 and CD8+ T cells (partial cor =−0.158, p = 1.19e−03), TRUB1 and CD8+ T cells (partial cor =−0.206, p = 1.87e−02), DARS and neutrophils (partial cor = 0.245, p = 3.73e−07), GDI2 and neutrophils (partial cor =0.184, p = 1.62e−04), GDI2 and Dendritic cells (partial cor =0.167, p = 6.21e−04) P4HA2 and B cells (partial cor =−0.239, p = 7.89e−07), P4HA2 and CD8+ T cells (partial cor =−0.158, p = 1.19e−03), TRUB1 and CD8+ T cells (partial cor =−0.206, p = 1.87e−02), DARS and neutrophils (partial cor = 0.245, p = 3.73e−07), GDI2 and neutrophils (partial cor =0.184, p = 1.62e−04), GDI2 and Dendritic cells (partial cor =0.167, p = 6.21e−04) (Fig. 6A). We also explored the correlation between tumor immune cell infiltration and somatic copy number alterations (SCNAs). The samples were divided into four types according to the copy number of genes.The distribution of infiltrating immune cells among the four types of samples was compared, as shown in Fig. 6B. We found that the various forms of mutations carried by the DARS/GDI2/P4HA2 / TRUB1 gene can usually suppress immune infiltration, including CD8+T cells, neutrophils, dendritic cells, macrophages, CD4+T cells, and B cells. Also, we found that these four pivotal genes had a greater effect on immune infiltration than alterations in the genes.

Genomic alterations of DARS/GDI2/P4HA2/TRUB1 in GBM

We then used the cBioPortal tool to determine the types and frequency of DARS/GDI2/P4HA2/TRUB1 alterations based on DNA sequencing data from GBM patients. The genetic variation rates of DARS/GDI2/P4HA2/TRUB1 ranged from 0% to 4% (DARS was 4%, GDI2 was 1.4%, P4HA2 was 4%, TRUB1 was 0.0%). These alterations include Missense Mutation, mRNA High, mRNA Low, Amplification (AMP), and Deep Deletion (Fig. 7). In view of this, DARS and P4HA2 show potentially stronger cancer-driving properties at a higher mutation frequency. In contrast, TRUB1 is genetically stable and could potentially act as a stable biomarker.

Gene functional analysis

We downloaded the DARS/GDI2/P4HA2/TRUB1 gene sets from the Molecular signatures database (v7.0 version) and comprehensively evaluated the gene sets through GSVA. Our analysis showed that in the DARS gene set, 17 gene sets were up-regulated (t>, 1) and 14 gene sets were down-regulated (t<1). In GDI2, 13 gene sets were up-regulated and 21 gene sets were down-regulated. In P4HA2, 11 gene sets were up-regulated and 30 gene sets were down-regulated. In TRUB1, 21 gene sets were down-regulated and 14 gene sets were down-regulated (Figs. 8A–8D).