Identification of robust diagnostic and prognostic gene signatures in different grades of gliomas: a retrospective study

Data resource and preprocessing

In order to obtain the genes related to the prognosis of GBM patients, the phenotype information and gene expression data of tumor samples and those of normal samples were successively collected from TCGA database and The Genotype-Tissue Expression (GTEx) project, according to overall survival (OS) of glioblastoma (GBM) patients. The standardized form of TCGA data were collected from UCSC Xena (http://xena.ucsc.edu), where gene expression was recorded as the log2 transformed FPKM (Reads Per Kilobase of exon model per Million mapped reads). To validate whether the gene biomarkers were valid in Chinese people, two datasets were downloaded from CGGA (http://www.cgga.org.cn/). The following three kinds of samples were excluded: (1) patients with incomplete information of either phenotype or gene expression, (2) patients with OS time less than 30 days, (3) patients with recurrent tumor. In the TCGA database, GBM sample included 147 samples and 196 normal tissue that selected from GTEx database. In the CGGA database, three datasets were included in this database. Of these datasets, CGGA_639 and CGGA_325 were tested by RNA-Seq. Another one was performed by Affymatrix chip. So, we selected two datasets that used the same platform.

Differential expression analysis

To identify the genes related to overall survival of GBM patients, the differentially expressed mRNAs between GBM patients and normal individuals were selected at first. The Limma package of R 3.6.1 was employed to analyze the data (https://www.r-project.org/). The threshold was set as follows: the adjusted p value was less than 0.001, abstract of log2-fold-change was larger than 1, B>5 and AveExpr>5. The p value was calculated via Student-t test and adjusted by Benjamini–Hochberg method (Puoliväli, Palva & Palva, 2020). The differentially expressed genes (DEGs) could be explained as genes expressed differently between GBM and normal tissue.

The predictive prognostic genes

Univariate Cox (Uni-Cox) proportional hazard regression and least absolute shrinkage and selection operator (LASSO) were employed to screen the genes that are related to survival among TCGA GBM patients. Uni-Cox was applied to identify the independent effect related to overall survival of each DEG. The hazard ratio (HR) of each mRNA was calculated according to following equation: (1)${\displaystyle \mathrm{HR}=e^{\beta }}$ where β represents the coefficient from Uni-Cox. Here Survival package of R was applied.

Then to simplify the predictive genes, least absolute shrinkage and selection operator (LASSO) method was adopted, which further filtered the statistically significant mRNAs further (Zhang et al., 2019). Package glmnet of R was applied.

Prognostic index model construction

Based on the identified prognostic genes, a risk value named prognosis index (PI) was calculated for every patient as an integrated signature. PI was calculated as follows: (2)${\displaystyle \mathrm{PI}=\sum _{i=1}^m{\beta }_i\times E_i}$ where β_i represents the coefficient of the involved gene i, and E_i represents the corresponding gene’s expression level. Then normalized PI was calculated as follows: (3)${\displaystyle \text{normalized}P{I}_i=\frac{P{I}_i-mean\left(PI\right)}{sd\left(PI\right)}}$ where mean(PI) represents the mean and sd(PI) represents the standard deviation of all PIs, respectively. The PIs below all referred to normalized PI. In order to classify the patients of high risks and low risks, the median of PIs was set as a cut-off. If the PI of a patient is larger than the cut-off, the patient is going to be assigned to high risk group, and be predicted with a bad OS, otherwise the patient will be predicted to have low risks of death.

Analysis of clinical confounding factors with PI

Several clinical characteristics might correlate with the OS of GBM patients. In TCGA database, age, KPS score, cancer status, gender and race were considered as the main clinical confounding factors affecting prognosis (Xiong et al., 2014). Therefore, we explored the role of age, KPS score, cancer status, gender, race and PI in TCGA datasets via Uni-Cox. Following this, multivariable Cox (Multi-Cox) proportional hazard regression was also conducted to explore the joint effect of these clinical factors and the previously calculated PI, which allowed comparison of the prognostic value of PI to that of each clinical factor. In the CGGA database, clinical variables contain more information, which is helpful to understand the influence of various confounding factors on Cox regression. In the CGGA database, it mainly includes grade, gender, age, radiotherapy status, chemotherapy status, IDH mutation status and 1p19q codeletion status. The inclusion and exclusion criteria of the cases to enter into the multivariable model were listed as follows: (1) patients with clinical data were selected; (2) untreated, primary (de novo) GBM patients were selected; (3) patients with history of neoadjuvant treatment were not admissible; (4) patients more than 30-day survival were selected. Variables with p-values < 0.05 were selected as candidates entering into the multivariable model.

Survival analysis and model test

To evaluate the predictive effect of PI, Kaplan–Meier survival curves were created. According to the log-rank test, the PI is considered as a good predictor if the p value is less than 0.05. Package Survminer of R was used.Then the time-dependent receiver operating characteristic (ROC) curve analysis was introduced to assess the model. Package survival ROC of R was applied and the area under the curve (AUC) was calculated. If AUC equals to 0.5, it indicates that the predictive effect of the model is tantamount to random allocation of patients. But if AUC is more than 0.5, it implies that its predictive effect is superior to random allocation (Kottas, Kuss & Zapf, 2014).

Geno ontology (GO) functional enrichment

To analyze the basic biological function of genes identified in our model, package ClusterProfiler (Yu et al., 2012) of R was adopted to conduct the GO analysis of functional process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The p-value cut-off is less than 0.05 as significant enrichment threshold.

Validation in CGGA dataset

Because most of the patients in the sample set is white, and only 8 samples are Asians (accounting for 5.5% of all samples), whether the model is valid for Chinese patients is unknown. Therefore, the validation of the model was conducted in two datasets collected from CGGA, one of which named CGGA_325 (a dataset with 325 Chinese GBM patients) and the other is CGGA_693 (a dataset with 693 Chinese GBM patients), respectively. The two datasets contain different sample numbers, among which CGGA_693 contains 693 samples, CGGA_325 contains 325 samples. The two datasets also contain different types of samples. CGGA_693 database contains primary LGG, recurrent LGG, primary GBM, recurrent GBM. CGGA_325 database contains primary LGG, recurrent LGG, primary GBM, recurrent GBM and secondary GBM. The PI of each patient was calculated according to the mRNAs identified in the model, and then K-M curves and ROC curves analysis were performed. Meanwhile, several clinical factors recorded in CGGA datasets were added into multivariable COX model to validate the predictive effects of PI when adjusted for confounders.

Human tissue samples

A total of 30 glioma tissue specimens and 10 cases of peritumor brain tissue (used control) resected by corticostomy were collected from Lanzhou University Second Hospital between October 2019 and May 2020. All these patients did not receive preoperative radiotherapy, chemotherapy or other immunotherapy, and were confirmed by surgical pathology. According to the WHO classification, there were 15 cases of WHO Grade II. Among them, 9 cases are of males, 6 cases are of females and the age ranges from 21 to 59 (40.5 ± 9.6). There were 9 cases with diffuse astrocytoma (A), 4 cases with oligodendroglioma (O), 2 cases with oligoastrocytoma (OA). There were 15 cases of WHO Grade III/IV, among which 6 cases are of males, 9 cases are of females, and the age ranges from 29 to 69 (47.9 ± 12.9). There were 4 cases with Anaplastic astrocytoma (AA), 1 case with anaplastic oligodendroglioma (AO), 4 cases with anaplastic oligoastrocytoma (AOA), 6 cases with glioblastoma (GBM). Among control groups, 5 cases are of males, 5 cases are of females, and the age ranges from 31 to 69 (49.7 ± 11.4). This study was approved by the Ethics Committee of Lanzhou University Second Hospital (2020A-147), and all patients signed informed consent before surgery.

RNA extraction and qRT-PCR

Total RNA from all cells and tumor tissues were isolated using Trizol reagent (Takara, Dalian, China) and the first-strand cDNA was converted using the PrimeScript RT reagent Kit with genomic DNA Eraser (Takara). Then, TB Green Premix ExTaq (Takara) was used to perform qRT-PCR (Bio-Rad CFX96). Based on the results of differential analysis, GAPDH was selected as reference gene in glioma samples. Relative gene expression was evaluated by a comparative CT method (2-ΔΔCtmethod). Statistical significance of qRT-PCR data was analyzed using IBM SPSS Statistics 22.0 software(Armonk, NY, USA) and determined by the Student’s t-test. p < 0.05 was considered to be statistically significant. Primers used were as following: HSPA5 forward: 5′-GACATCAAGTTCTTGCCGTTCA-3′, HSPA5 reverse: 5′-CCAGCAATAGTTCCAGCG TCTT-3′; MTPN forward: 5′-CGGAGACTTGGATGAGGTGAA-3′, MTPN reverse: 5′-AGAGCTTTGATTGCCTGGTTG-3′; GAPDH forward: 5′-GGAAGCTTGTCATCAAT GGAAATC-3′, GAPDH reverse: 5′-TGATGACCCTTTTGGCTCCC-3′.

Western blot assay

Tissue samples were lysed in RIPA lysis buffer, and lysates were harvested by centrifugation (12,000 rpm) at 4 °C for 30 min. Western blot was performed in accordance with the protocols as described above, using β-actin as the internal control. Briefly, the protein sample was separated by SDS-PAGE electrophoresis and then transferred to a PVDF membrane. After blocking nonspecific binding sites for 60 min with 5% non-fat milk, the membrane was incubated with the primary antibody at 4 °C overnight. Membranes were washed three times with tris buffered saline with 1‰tween-20 and incubated with horseradish peroxidase-conjugated secondary antibody at 37 °C for 1 h. After 3 washes, the bands were detected by an enhanced chemiluminescence system (WBKLS0500, Merck KGaA). Band density was measured using ImageJ software (National Institutes of Health, Bethesda, MD) and standardized to that of β-actin. Antibodies against proteins and dilution multiples are as follows: HSPA5(ab21685, abcam, USA) (1:1000), MTPN(bs-11891R,bioss,China) (1:1000), β-actin (service, China)(1:3000).

Immunohistochemical assay

Two-Step method was adopted to detect the expression of HSPA5 and MTPN. Tissue specimens were fixed promptly with 100g/L formaldehyde solution, embedded in paraffin and cut into 3 µm sections. All sections were dehydrated with graded alcohol, repaired with antigen, quenched with peroxidase solution, put into tap water, and washed with PBS. The slides were incubated in a moist chamber with HSPA5 or MTPN rabbit polyclonal antibody (1:100) at 37 °C for 30min. Then they were incubated in a moist chamber with the goat polyclonal antibody against rabbit at 37 °C for 30min. After being washed in PBS completely, the slides were developed in 0.05% freshly prepared diaminobenzedine solution (DAB) for 10 min, and then counterstained with hematoxylin. Finally, the slides were dehydrated in ascending concentrations of ethanol, airdried, and mounted. The expression of HSPA5 and MTPN were scored according to the degree of staining and the number of stained cells. Degree of staining: 0 for non-staining, 1 for light yellow, 2 for brownish, and 3 for tan. Stained cell counts are defined as follows: Under high power microscope, 5 fields were randomly selected from each section to count the percentage of stained cells; 0 point is given to match the stained cells which represent less than 5%; 1 point matches 5% to 25%; 2 points matches 26%∼50%; 3 points matches 51%∼75%; 4 points matches more than 75%. The product of the two scores is scored:0 is negative, 1-3 is weakly positive, 4-5 is moderately positive, and ≥6 is strongly positive. Negative and weak positive were classified as negative expression, while moderate positive and strong positive were classified as positive expression.

Statistical analysis

Statistical analysis of human tissue data was expressed as mean ± SD (standard deviation) from three independent experiments. Differences between groups were estimated using the Student’s t-test. Discrete data values were expressed as rate and analyzed by chi-squared test. All these analyses were performed using IBM SPSS Statistics 22.0 software (Armonk, NY, USA) and a two-tailed value of P < 0.05 was considered statistically significant.

DEGs Between GBM samples and Non-GBM samples

By preprocessing, 147 GBM patients with 19,199 genes in TCGA and 196 non-GBM individuals with 19199 genes in GTEx were collected to do the following analysis.

Through the analysis of differential gene expressions, 581 DEGs were left in the model according to the threshold. Among them, 138 mRNAs were highly expressed and 443 mRNAs were lowly expressed. Figure 1 shows the volcano plots of the results of DEGs. For further biological experiment, we filtered average expression more than 5 as threshold. The gene significantly expression showed in Fig. 1A. Filtering genes were showed in Fig. 1B.

The survival related genes of GBM

Through Uni-cox regression, 50 mRNAs were identified as survival related genes according to p value less than 0.05 in TCGA database. Then LASSO method was employed and 20 mRNAs were finally selected. The information of these 20 genes were listed in Table 1. Among them, 16 mRNAs were high-risk genes with HRs from 1.230 to 2.039, which meant that the higher their expression, the worse OS the patient might have. In contrast, 4 genes were protective factors with HRs from 0.789 to 0.552.

PI construction

According to Eqs. (2) and (3), the PI of each patient was calculated and the cut-off was set to 0.207, which implied patients with PIs larger than 0.207 were considered with high risks of overall survival. Ordered PIs were depicted in Fig. 2. The HR of PI, which was displayed in the last line of Table 1, was 2.653 with 95% CI from 1.976 to 4.122. The p-value of Uni-Cox of PI was less than 0.001.

Clinical factors and their joint effect with PI

Table 2 summarizes the results of clinical factors analysis. Five clinical factors were included in the analysis according to the characteristics of the original training set. They are age, karnofsky performance score (KPS), cancer status, gender and race of patients. Based on Uni-cox and survival analysis, it showed that age, KPS and cancer status were independent significant factors to OS of GBM patients. Then, through multivariate analysis, only KPS, cancer status and PI were statistically significant. HR of PI is 2.653 (95% CI [1.780–3.955], p < 0.001), which demonstrated PI is significant risk factor.

Survival analysis and model test

Figure 3 shows the results of K-M curves of PI. The red line was the survival curve of the low-risk patients and the green line was that of high-risk patients. For 500 days, 31 patients (41%, 31/74) were alive in the low-risk group and 11 patients (15%, 11/73) in the high-risk group. For 1,000 days, all patients in the high-risk groups were dead, while 10 patients (13%, 10/74) in the low-risk group were still alive. And the p value of log-rank test was less than 0.001. The right term of the figure was the ROC curves of our model. The AUCs of 3 years and 5 years were 0.824 and 0.820, respectively, which represented that our model had a good prognostic effect for GBM.

Results of GO enrichment analysis

According to GO functional enrichment analysis, 8 biological processes were enriched with adjusted p value less than 0.05 (Fig. 4). They were cell adhesion molecule binding, chaperone binding, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding, cadherin binding, translation initiation factor activity, translation initiation factor activity, translation factor activity-RNA binding and unfolded protein binding. DEGs mainly involved in cell adhesion molecule binding, chaperone binding, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding and cadherin binding. These genes are engaged in the biological processes of various brain functions such as cerebellum morphogenesis, hindbrain morphogenesis, cerebellar cortex development and etc. In molecular function, these genes play a part in binding function.

Validation in CGGA dataset

In validation section, we applied above results from TCGA to test in CGGA datasets. The Uni-Cox results of 20 mRNAs were listed in Table 3. It showed us that among 20 identified genes through TCGA, 4 genes (MTPN, RTN4, MAP1LC3A and DKK3) were significantly associated with OS. Although MTPN is not remarkably associated with OS, the results of calculation showed that HR>1 of MTPN and p-value showed a significant trend.

Seven of these genes only were statistically significant via Uni-Cox in either CGGA_693 or CGGA_325 (LY6E, FLII, PLD3, EIF4A2, RNF10, RPS19 and NDUFB2). The other 9 genes (TTYH3, TAGLN2, FN1, HSPA5, EI3L, SCG5, SERPINE2, FDPS and GRN) were still prognostic genes in both CGGA_693 and CGGA_325. Among these nine genes, SCG5 and SERPINE2 showed the opposite effects in CGGA and TCGA, which might because of racial difference. And the seven genes mentioned above showed the same effects (protective factor or risk factor) between TCGA and CGGA, which might be the real prognostic genes of GBM. The HRs of PI in CGGA_693 and CGGA_325 were 5.200 (95%CI [3.872–6.983]) and 6.524 (95%CI [4.496–9.466]), respectively. Thus our 20-gene model still had a good predictive effect among Chinese GBM patients. As is shown in Table 4, WHO clinical grade and 1p19q codelection status were also significant high-risk clinical factors of GBM patients for overall survival. Accounting for clinical confounders, PI was still a critical prognostic predictor in both CGGA_693 and CGGA_325 datasets.

Verification gene signature performance of CGGA dataset

Gene signature performance was analyzed in two CGGA datasets. Survival curve and ROC curve were employed to test prediction performance of the gene signature. Figure 5 showed that the p-values of log-rank test of K-M curves were less than 0.0001 both in CGGA_693 and CGGA_325. In CGGA_693, 185 patients (74%, 185/248) were alive in the low-risk group and 44 patients (27%, 44/161) in the high-risk group for 1,000 days. And all patients in the high-risk group were dead, while 11 patients (4%, 11/248) in the low-risk group were still alive for 3,000 days. In CGGA_325, 104 patients (79%, 104/131) were alive in the low-risk group and 18 patients (19%, 18/93) in the high-risk group for 1,000 days. And all patients in the high-risk group were dead, while 15 patients (11%, 15/131) in the low-risk group were still alive for 4,000 days. This implied that the biomarkers could efficiently classify the Chinese GBM patients into good and poor prognosis groups. ROC analysis represented that the AUCs of CGGA_693 were 0.831 and 0.808 for 3 and 5 years, respectively, and AUCs of CGGA_325 were 0.907 and 0.912 for 3 and 5 years, respectively. Therefore, our forecasting model for GBM had a good classificaion ability in CGGA datasets. To test the predictive effects of our 20-gene model among Chinese GBM patients further, K-M curves and ROC analysis were conducted.

The expression of mRNA HSPA5 and MTPN by PCR

By real-time quantitative PCR detection, compared with normal brain tissue (0.96 ± 0.28), the expression of mRNA HSPA5 in WHO II gliomas and WHO III/IV glioma tissue both increased significantly (1.25 ± 0.32,1.85 ± 0.70), and the differences were statistically significant (P < 0.05, Fig. 6).

Compared with normal brain tissue (0.60 ± 0.28), the expression of mRNA MTPN in WHO II gliomas and WHO III/IV glioma tissues also increased greatly (0.95 ± 0.19,1.01 ± 0.39), and the difference was statistically significant (P < 0.05, Fig. 6).

We filtered all gene signature in three datasets. All genes were significantly changed in three datasets. These genes included six genes (TAGLN2, HSPA5, FN1, TTYH3, GRN and MTPN). Of these genes, TAGLN2 (Beyer et al., 2018; Han et al., 2017), FN1 (Yu et al., 2020; Gu, Gu & Shou, 2014; Guo, Heller & Thorslund, 2016; Liao et al., 2018), TTYH3 (Weinberg et al., 2020), GRN, (Ness, Riemenschneider & Baches., 2009; Trigos et al., 2019) have been reported in GBM.

Western blot analysis

Western blotting showed HSPA5 protein level in WHO II glioma tissue and WHO III/IV glioma tissue were both elevated obviously compared with normal brain tissue. MTPN protein expression in WHO II glioma tissue and WHO III/IV glioma tissue were also significantly higher than normal brain tissue (Fig. 7).

Immunohistochemical analysis

HSPA5 protein was mainly expressed in the cytoplasm while MTPN protein was expressed in the nucleus and cytoplasm.

Immunohistochemical results indicated that compared with normal brain tissue, the positive expression of HSPA5 protein in gliomas rose obviously. For the most part, the expression of protein was strongly positive in grade III/IV glioma tissue, moderately positive in grade II glioma tissue, and mostly negative or weakly positive in normal brain tissue (Figs. 8A–8C, 8G). It is equally true of the expression of MTPN protein in gliomas (Figs. 8D–8F, 8H).