Identification of the Prognostic Signatures of Glioma With Different PTEN Status.

The high-grade glioma is characterized by cell heterogeneity, gene mutations, and poor prognosis. The deletions and mutations of the tumor suppressor gene PTEN (5%-40%) in glioma patients are associated with worse survival and therapeutic resistance. Characterization of unique prognosis molecular signatures by PTEN status in glioma is still unclear. This study established a novel risk model, screened optimal prognostic signatures, and calculated the risk score for the individual glioma patients with different PTEN status. Screening results revealed fourteen independent prognostic gene signatures in PTEN-wt and three in the -50PTEN-mut subgroup. Moreover, we verified risk score as an independent prognostic factor significantly correlated with tumor malignancy. Due to the higher malignancy of the PTEN-mut gliomas, we explored the independent prognostic signatures (CLCF1, AEBP1, and OS9) for a potential therapeutic target in PTEN-mut glioma. We further separated IDH wild-type glioma patients into GBM and LGG to verify the therapeutic target along with PTEN status, notably, the above screened therapeutic targets are also significant prognostic genes in both IDH-wt/PTEN-mut GBM and LGG patients. We further identified the small molecule compound (+)-JQ1 binds to all three targets, indicating a potential therapy for PTEN-mut glioma. In sum, gene signatures and risk scores in the novel risk model facilitate glioma diagnosis, prognosis prediction, and treatment.

Introduction

The most common primary brain tumor, glioma, starts with glial support cells around nerve cells (1, 2). Clinically, glioma is associated with high mortality and recurrence rates and poor prognosis (3, 4). Although surgery, radiotherapy, and alkylated chemotherapies are used in treatment, patient survival rate did not significantly improve in the past decade (5). Genomics, transcriptomics, and epigenetic analyses gave rise to new concepts for the molecular classification and treatment of gliomas (6). Thus, the molecular signature is important for the tumor’s diagnosis, patient stratification, and personalized treatment. In 2016, the World Health Organization (WHO) incorporated morphology and genetic variations to update guidelines for classifying brain tumors, particularly gliomas (7). Many studies have shown that glioma patients with IDH (Isocitrate dehydrogenase) mutations have a better prognosis (8–14). Thus, the WHO included IDH mutation as an essential diagnostic and classification criterion for glioma (7). Additionally, MGMT (O6-methylguanine-DNA methyltransferase) promoter methylation is another example of a glioma signature. MGMT was initially identified as a prognostic and predictive signature for glioma diagnosis in patients treated with temozolomide (15). Therefore, signatures in glioma have been highly valued as prognostic, predictive, and clinical application targets.

WHO divided glioma into grades I, II, III, and IV. The higher the glioma grade, the worst the prognosis (16–19). Even though grade IV glioma (glioblastoma multiforme, GBM) is considered different from low-grade glioma (LGG) in general, many studies still combine LGG and GBM dataset to screen tumorigenesis marker or therapeutic target by considering GBM is the most malignant form of glioma (6, 20). Notably, 50% of GBM harbor somatic alterations in the phosphatidylinositol-3-kinase (PI3K) pathway (21, 22). The tumor suppressor PTEN is the most critical negative regulator of the PI3K pathway (23). PTEN is significantly altered in GBMs (30–40%), and mutational loss of PTEN function is an important malignancy event in glioma. The meta-analysis by Han et al. showed that PTEN mutations are strongly associated with shorter survival in glioma patients, suggesting that PTEN status strongly correlates with the prognosis of patients (24). Additionally, the loss of PTEN increases drug resistance. For example, (1) GBM patients with PTEN mutation have no significant response to anti-PD-1 (Programmed cell death 1) immunotherapy due to the mutation-induced changes in the immune microenvironment (25); (2) patients with PTEN-negative GBM have a shorter survival time after initiation of bevacizumab than those with PTEN-positive GBM (26); (3) loss of PTEN leads to clinical resistance to PI3Kα inhibitors (27); (4) fibroblast growth factor receptor 2-mediated phosphorylation of PTEN at tyrosine 240 contributes to the radioresistance of glioma (28).

Detecting the PTEN status of glioma patients and treating them separately according to their unique signatures may reduce drug resistance and improve the survival of patients. In this study, we divided glioma patients into PTEN-wt and PTEN-mut subgroups. We observed that survival of the PTEN-mut is significantly lower than the PTEN-wt subgroup in both TCGA and the CGGA glioma datasets, indicating the correlation between PTEN status and the prognosis of gliomas patients. We established prognostic risk models by fitting the L1-penalized (LASSO) Cox-PH regression model, obtained the optimal prognostic signatures, and calculated risk scores in these two glioma subgroups. Due to the higher malignancy of the PTEN-mut glioma, we look for potential therapies targeting the prognostic signatures (CLCF1, AEBP1, and OS9) in this subgroup. We then focused on IDH-wt/PTEN-mut GBM and LGG patients and found that the aforementioned screened genes also significantly prognostic in both subgroups. Finding tissue/cancer-specific miRNAs, TFs (transcription factors), and the small molecular compound is critical in the tumor treatment (29–33). Herein, we provide some miRNAs and TFs binding to the optimal prognostic genes, and small molecular compounds binding to the respective proteins in the PTEN-mut subgroup. Importantly, we found that (+)-JQ1 compound interacting with these three optimal prognostic genes as a potential drug.

Materials and Methods

Data Source

The Cancer Genome Atlas Program (TCGA) glioma project is the training dataset. The RNA-seq data was downloaded from https://tcga.xenahubs.net/download/TCGA.GBMLGG.sampleMap/HiSeqV2.gz, including 697 glioma samples and 5 Normal samples. This data, which combines TCGA brain lower grade glioma and glioblastoma multiforme datasets, shows the gene-level transcription estimates, as in log2 (x+1) transformed RSEM Normalized count. The clinical information of 1148 patients was downloaded from https://tcga.xenahubs.net/download/TCGA.GBMLGG.sampleMap/GBMLGG_clinicalMatrix. Finally, the mutation data was download from https://gdc.xenahubs.net/download/TCGA-LGG.varscan2_snv.tsv.gz and https://gdc.xenahubs.net/download/TCGA-GBM.varscan2_snv.tsv.gz, including 506 low-grade gliomas samples and 390 glioblastomas samples, respectively. Two types of somatic mutations in the dataset are single nucleotide polymorphisms (SNPs) and small insertion/deletions (INDELs). Each patient must have clinical, RNA-seq, and mutation information; thus, the 653 overlapping patients in three data is our training and research dataset ( ).

The Chinese Glioma Genome Atlas (CGGA) glioma project is the validation dataset. The RNA-seq dataset of CGGA was download from http://cgga.org.cn/download?file=download/20200506/CGGA.mRNAseq_693.RSEM-genes.20200506.txt.zip&type=mRNAseq_693&time=20200506 (expression data from STAR+RSEM), and corresponding clinical information was downloaded from http://cgga.org.cn/download?file=download/20200506/CGGA.mRNAseq_693_clinical.20200506.txt.zip&type=mRNAseq_693_clinical&time=20200506, which contains 693 glioma samples. The mutation dataset was download from http://cgga.org.cn/download?file=download/20200506/CGGA.WEseq_286.20200506.txt.zip&type=WEseq_286&time=20200506, containing 286 glioma samples. A total of 144 overlapping patient samples from the three datasets were used as a validation dataset ( ).

Differentially Expressed Genes

Patients were divided into two glioma subgroups according to the status of PTEN ( ), including PTEN-wt and PTEN-mut patients, whose clinical features are respectively presented in . Differentially expressed genes (DEGs) between gliomas of our dataset and Normal samples (DEGs-all from TCGA glioma vs. Normal, 653 vs. 5; DEGs-wt from TCGA PTEN-wt vs. Normal, 575 vs. 5; DEGs-mut from TCGA PTEN-mut vs. Normal, 78 vs. 5) were screened ( ) using the limma package (version 3.40.6) of R 3.6.1 (34). Fold discovery rate (FDR) < 0.05 and |log2FC (fold change) | > 1 was as the selection criteria. Based on the expression values of the DEGs, the clustering analysis was used by heatmap (version 1.0.12) (35).

Figure 1

Data analysis flow (Materials and Methods).

Table 1

Clinical information of the PTEN-wt and PTEN-mut subgroups (*** significant difference).

Group	PTEN_wt	PTEN_mut	p-value
Patient(number)	575	78
Age(median age)	44	59	<0.0001***
OS.time(median days)	576	461	<0.0001***
Male	57% (330)	59% (44)
Female	43% (245)	44% (34)
Chemotherapy	19% (108)	14% (4)
Radiation therapy	9% (50)	3% (1)
Astrocytoma	30% (171)	30% (24)
Oligoastrocytoma	22% (124)	14% (11)
Oligodendroglioma	31% (180)	8% (6)
Glioblastoma	17% (100)	63% (49)
IDH1_mutation	20% (116)	12% (9)

Prognostic Risk Models

We established prognostic risk models using DEGs in different PTEN statuses and machine learning. In brief, we obtained prognostic DEGs (PR-DEGs) that were significantly linked to survival time by calculating the log-rank test and univariate Cox regression ( ). A gene with log-rank p < 0.05 and Wald test p < 0.05 was considered as a significant PR-DEG. These identified PR-DEGs were then used to fit L1-penalized (LASSO) Cox-PH regression model (35) for the selection of the optimal prognostic DEGs (OPR-DEGs) with glmnet package (version 3.0.2) of R ( ). The optimal parameter of ‘lambda’ in the selection model was calculated via the cross-validation likelihood (CVL) method 1000 times.

A prognostic risk model is established from the DEGs’ expression levels and Cox regression coefficients ( ). The risk score of every patient was established by the following formula (36):

βRNAn and ExpRNAn represent the Cox-PH coefficient and the RNA expression level of ORP-DEGs.

Each subgroup set was classified into high-risk patients and low-risk patients by median risk score. After that, the Kaplan-Meier method was used to draw the survival curves between subgroups. Additionally, the accuracy for survival prediction using these models was evaluated by the area under the receiver operating characteristic curve (AUC). The CGGA dataset was used to validate the prognostic risk model ( ).

Multivariate Cox Regression Analysis and Nomogram Survival Rate Model

Risk scores and clinical features were used for Multivariate Cox regression analysis (survival package, 3.1-8 in R) and Nomogram survival rate model (rms package, 5.1-4 in R) (37), to construct a nomogram 1-year, 3-year, 5-year, 10-year, and 15-year survival rate model ( ). Bootstrap resamples (1000) in the TCGA dataset were used for internal validation, and the CGGA dataset was used for external validation to evaluate the predictive effect. The value of the C-index indicates the accuracy of the predictive ability.

Generation of ESTIMATEScore and Tumor Purity Score

The ESTIMATE algorithm loaded in the estimate package (1.0.13) in R is used to calculate the ratio of the immune/stromal component in the tumor microenvironment, exhibited by four kinds of scores: Immune Score, Stromal Score, ESTIMATEScore, and TumorPurityScore ( ).

Function Analysis, Protein-Protein Interaction, Gene-miRNA and TFs Interaction and Protein-Chemical Interaction

In function analysis using the clusterProfiler package (version 3.12.0) in R, we analyzed the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for identified DEGs (p < 0.05). The protein-protein interactions (PPIs) of OPR-DEGs were analyzed based on the STRING database (Version: 11.0) (38). The gene-miRNA, gene-TFs, and protein–chemical interactions were constructed by NetworkAnalyst (39) and visualized using Cytoscape (Version 3.6.0) (40).

Dock Between Proteins With Small Molecular Compound

Autodock (Version 4.2) was used to dock proteins and small molecular compounds, and VDM (Version 1.9.3) was used to visualize docking sites.

Cell Culture

Colony formation experiment was used to analyze the effect of (+)-JQ1 on GBM cell proliferation. Cells (1000/well) of U251 (GBM cell line, PTEN-deficiency) and U343 (GBM cell line, PTEN-wt) were seeded into 12-well plates. After overnight culture, fresh medium with different concentrations of (+)-JQ1 was replaced. After 2 weeks of treatment, the medium was pumped, the cells were flushed with PBS 3 times, 2ml 4% paraformaldehyde was added to fix the cells (10 minutes), and the cells were washed with PBS for 3 times. After dyeing with 0.5% crystal violet for 10 minutes, colony formed.

Results

The PTEN-mut With Worse Survival Than the PTEN-wt Subgroup

The PTEN-mut group presented the worst survival. The TCGA dataset showed an 18% PTEN mutation rate in glioma ( ). Of the top five mutation genes in glioma, including TP53, IDH1, TTN, ATRX, and PTEN, the PTEN mutation is the most significant in reducing survival ( and ). As expected, the CGGA dataset showed worse survival in PTEN-mut gliomas compared to the PTEN-wt subgroup ( ). To identify survival differences among mutated PTEN glioma patients, we extracted the PTEN mutation sites list ( ) of TCGA glioma patients. The survival rates of the glioma patients with loss of function mutation (57%) and patients with non-reported mutation sites (43%) have no difference ( ). Specifically, there was no significant difference (p = 0.78) in survival time between the truncated PTEN protein at the arginine 130 (R130, Complete loss of function; median survival time: 345; n = 4) and point mutation glutamine (R130Q, Partial loss of function; median survival time: 360; n = 3).

Figure 2

The characteristic of PTEN-wt and PTEN-mut subgroup: (A) The top 5 gene mutation rate in TCGA glioma; (B) Survival analysis (Kaplan Meier Curve) of PTEN-wt (blue, n = 575) and PTEN-wt subgroup (yellow, n = 78) in TCGA dataset; (C) Survival analysis (Kaplan Meier Curve) of PTEN-wt (blue, n = 134) and PTEN-wt subgroup (yellow, n = 9) in CGGA dataset; (D) The principal component analysis of TCGA PTEN-wt (blue, n = 575) and TCGA PTEN-mut subgroup (yellow, n = 78 ); (E) Pathway analysis of DEGs-all (DEGs from TCGA glioma vs. normal tissue, 653 vs 5); (F) Pathway analysis of DEGs-wt (DEGs from TCGA PTEN-wt glioma vs. normal tissue, 575 vs 5); (G) Pathway analysis of DEGs-mut (DEGs from TCGA PTEN-mut glioma vs. normal tissue, 78 vs 5); (H) The overlap of DEGs-all, DEGs-wt, and DEGs-mut; (I) The GESA enrichment plot of cellular component for DEGs only in PTEN-wt (n = 147); (J) The GESA enrichment plot of cellular component for DEGs only in PTEN-mut (n =1938). DEGs-all: DEGs from TCGA glioma vs. normal tissue ; DEGs-wt: DEGs from TCGA PTEN-wt glioma vs. normal tissue; DEGs-mut: DEGs from TCGA PTEN-mut glioma vs. normal tissue. ***p < 0.001.

The clinical information of the PTEN-wt and the PTEN-mut subgroups was statistically analyzed ( ). The median age of the PTEN-wt subgroup was 44, and the PTEN-mut subgroup was 59, showing a significant statistical difference between the two subgroups ( ). In addition, PTEN-wt contains 17% GBM, while the poorer survival PTEN-mut subgroup includes 63% GBM, indicating that malignant GBM is more prone to PTEN mutation ( ). Thus, all data show that PTEN is a critical indicator of glioma prognosis.

Patients in different subgroups showed significant differences in gene expression as indicated by the principal component analysis ( ). To identify the subgroup-specific DEGs of the PTEN-wt and PTEN-mut subgroups, the DEGs-all, DEGs-wt, and DEGs-mut (see Materials and methods) were calculated ( ). Notably, the pathway analysis of DEGs ( ) showed the enriched mitogen-activated protein kinase (MAPK) signaling pathway in the PTEN-mut subgroup ( ). MAPK signaling pathway is a set of the evolution of conservative serine/threonine-protein kinase, which controls many physiological activities, such as inflammation, apoptosis, canceration, invasion, and metastasis of tumor cells (41). The intersection of DEGs-all, DEGs-wt, and DEGs-mut in the distribution diagram showed that there were 3537 overlapping DEGs in both subgroups. However, 147 PTEN-wt group-specific DEGs and 1938 PTEN-mut subgroup-specific DEGs were also found ( , and ), suggesting gliomas with PTEN mutation are non-homeostatic. Gene set enrichment analysis (GSEA) showed 147 DEGs were mainly down-regulated only in DEGs-wt ( , ) without significant enrichment in the Cellular Component term. Notably, most of the 1938 DEGs were upregulated only in DEGs-mut ( , ), having significant enrichment in the extracellular region of Cellular Component term, which is closely related to tumor progression ( , FDR < 0.05). All of these results indicate that the PTEN-mut glioma presents a more malignant state.

The Signature and Risk Score in the PTEN-wt Subgroup

To establish the prognostic risk model of the PTEN-wt subgroup, log-rank test and univariate Cox analysis were performed using subgroup-specific DEGs (n = 147) of PTEN-wt to obtain PR-DEGs (n = 127, log-rank p < 0.05 and Wald test p < 0.05). Then, the PR-DEGs were used to fit L1-penalized (LASSO) Cox-pH regression model to obtain the OPR-DEGs (n = 44, ). The AUC value of the L1-penalized (LASSO) Cox-pH regression model was equal to 0.8807 ( ), indicating the reliability of this model. We calculated the patients’ risk scores based on 44 OPR-DEGs according to the calculation formula cited in the materials and methods. To calculate the hazard of risk score, the risk score and clinical features of patients were analyzed by multivariate Cox analysis. The Concordance Index obtained by multivariate Cox analysis was 0.89, demonstrating the accuracy of our constructed risk score of polygenes and the prediction model in the PTEN-wt subgroup ( ). Notably, the multivariate Cox analysis showed the risk scores of patients were independent prognostic factors (p < 0.001), and the hazard ratio of patients with a high-risk score is 2.45-fold ( ). Finally, we calculated the nomogram by clinical features and the risk scores to predict clinical survival ( ). As shown in calibration plots by internal validation of Bootstrap Resamples (1000) of the TCGA-PTEN-wt dataset and external validation of the CGGA-PTEN-wt dataset ( ), consistent results were observed between predicted and actual survival probability at 1-year, 3-year, and 5-year times (C-index of internal validation = 0.894, C-index of external validation = 0.808).

Figure 3

Prognostic risk model in TCGA PTEN-wt subgroup: (A) Multivariate Cox forest plot of risk score and clinical features in TCGA PTEN-wt subgroup (n = 575 , AIC = 1588.11, C-index = 0.89); (B) A nomogram incorporating risk score based on 44 OPR-DEGs ( , risk score in ) and clinical features in TCGA PTEN-wt subgroup; (C) Survival analysis (Kaplan Meier Curve) of high-risk score patients (PTEN_wt_RS_H, blue, n = 287) and low-risk score patients (PTEN_wt_RS_L, yellow, n = 287) in TCGA PTEN-wt subgroup (high-risk score and low-risk score separated from the median risk score of TCGA PTEN-wt subgroup); (D) Survival analysis of high-risk score patients (blue, n = 67) and low-risk score patients (yellow, n = 67) in CGGA PTEN-wt subgroup (high-risk score and low-risk score separated from the median risk score of CGGA PTEN-wt subgroup); (E) The heatmap of expression of DEGs (n = 1468, ) between PTEN_wt_RS_H and PTEN_wt_RS_L, and the risk score is ranked from low to high; (F) Pathway analysis of DEGs (n = 1468, ) between PTEN_wt_RS_H and PTEN_wt_RS_L; (G) Expression of PTEN in PTEN_wt_RS_H and PTEN_wt_RS_L; (H–L) Expression of TDH, SSTR5, SEL1L3, LGR6, HMGN5 in PTEN_wt _RS_H and PTEN_wt _RS_L in TCGA. PTEN_wt_RS_H: high-risk score patients in TCGA PTEN-wt; PTEN_wt _RS_L: low-risk score patients in TCGA PTEN-wt; ***p < 0.001.

PTEN-wt subgroup in TCGA and CGGA dataset was separated into high-risk and low-risk patients by the median value of risk score. The Kaplan-Meier survival curve shows the prognosis of low-risk is significantly better than that of high-risk patients ( ; p < 0.001). We then screened DEGs between the high-risk patients and the low-risk patients and obtained 658 down-regulated, and 810 upregulated DEGs in the TCGA PTEN-wt subgroup ( , selection criteria: FDR < 0.05 and |log2FC| > 1, ). The patients’ risk score is ranked from low to high, and the DEGs between the high-risk patients and the low-risk patients are presented ( , ). These upregulated DEGs in the high-risk patients were significantly enriched in 10 functionally KEGG signaling pathways, including neuroactive ligand-receptor interaction, focal adhesion, extracellular matrix (ECM)-receptor interaction, and cell cycle ( , ). The clinical features of a high-risk and low-risk patient are shown in . We found the expression level of PTEN in high-risk patients was lower than low-risk patients in the PTEN-wt subgroup (p < 0.001), which showed that the risk score was also correlated with the level of PTEN expression ( ).

Table 2

Clinical information about risk high and risk low patients in PTEN-wt subgroups (*** significant difference).

Group	PTEN_wt Risk score high	PTEN_wt Risk score low	p-value
Patient(number)	287	288
Age(median age)	53	39	<0.0001***
OS.time(median days)	480	765	<0.0001***
Male	60%(172)	55%(158)
Female	40%(115)	45%(130)
Chemotherapy	24%(70)	14%(4)
Radiation therapy	10%(28)	3%(1)
Astrocytoma	31%(90)	28%(81)
Oligoastrocytoma	15%(43)	28%(81)
Oligodendroglioma	20%(56)	43%(124)
Glioblastoma	34%(98)	0.01%(2)

OPR-DEGs, which best fit the patient’s survival time, are potential targets for driving glioma progression. The 44 OPR-DEGs in the TCGA PTEN-wt subgroup have 14 OPR-DEGs (Multivariate Cox: p < 0.05, ) that are independent prognostic genes and have a strong relationship with patient survival. Thus, 14/44 OPR-DEGs are optimal survival and risk prediction indicators in the TCGA PTEN-wt subgroup, for example, SSTR5 and TDH ( ). Notably, 5/14 OPR-DEGs, including TDH, SSTR5, HMGN5, LGR6, and SEL1L3, also have different expressions between high-risk patients and low-risk patients in the TCGA PTEN-wt subgroup ( , p < 0.001), suggesting that they are associated with the prognosis of PTEN-wt glioma.

The Signature and Risk Score in the PTEN-mut Subgroup

We used the same prognostic risk model in the PTEN-mut subgroup and obtained 11 OPR-DEGs ( ). The AUC value was equal to 0.9414, indicating the model’s fine applicability ( ). Notably, the multivariate Cox analysis showed that patients’ risk score in the PTEN-mut subgroup is an independent prognostic factor. The hazard ratio of patients with a high-risk score is 3.08 ( ). The 11 OPR-DEGs in the PTEN-mut subgroup have three independent OPR-DEGs (Multivariate Cox: p < 0.05), including CLCF1, AEBP1, and OS9 ( ). Meanwhile, the nomogram for the PTEN-mut subgroup was also constructed to predict individual patient survival rates ( ). The calibration plots by internal validation of Bootstrap Resamples (1000) in the TCGA-PTEN-mut subgroup demonstrate the high availability of predictive models ( , C-index of internal validation = 0.822). As the CGGA-PTEN-mut subgroup has only nine patients, external validation could not be performed.

Figure 4

Prognostic risk model in TCGA PTEN-mut subgroup: (A) Multivariate Cox forest plot of risk score and clinical feature in TCGA PTEN-mut subgroup (n = 78 , AIC = 322.17, C-index = 0.83); (B) A nomogram incorporating risk score based on 11 OPR-DEGs ( , risk score in ) and clinical features in TCGA PTEN-mut subgroup; (C) Survival analysis (Kaplan Meier Curve) of high-risk score patients (PTEN_mut_RS_H, blue, n = 39) and low-risk score patients (PTEN_mut_RS_L, yellow, n = 39) in TCGA PTEN-mut subgroup (high-risk score and low-risk score separated from the median risk score of TCGA PTEN-mut subgroup); (D) Survival analysis of high-risk score patients (blue, n = 5) and low-risk score patients (yellow, n = 4) in CGGA PTEN-mut subgroup (high-risk score and low-risk score separated from the median risk score of CGGA PTEN-mut subgroup); (E) The heatmap of expression of DEGs (n = 476, ) between PTEN_mut_RS_H and PTEN_mut_RS_L, and the risk score is ranked from low to high; (F) Pathway analysis of DEGs (n = 476, ) between PTEN_mut_RS_H and PTEN_mut_RS_L; (G–I) Expression of CLCF1 (G), AEBP1 (H) and OS9 (I) in PTEN_mut _RS_H and PTEN_mut _RS_L. PTEN_mut_RS_H: high-risk score patients in TCGA PTEN-mut; PTEN_mut _RS_L: low-risk score patients in TCGA PTEN-mut; ***p < 0.001.

As we predicted, the prognosis of low-risk patients was significantly better than that of high-risk patients ( ; p < 0.0001). The clinical features of the high-risk and low-risk patients are shown in . There are 273 upregulated and 203 down-regulated DEGs between the high-risk patients and the low-risk patients in the TCGA-PTEN-mut subgroup ( , selection criteria: FDR < 0.05 and |log2FC| > 1, ). The heatmap presents DEGs’ expression between the high-risk patients and the low-risk patients ( , ). These upregulated DEGs in high-risk patients were significantly enriched in 5 KEGG signaling pathways, consisting of neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, complement, and coagulation cascades, staphylococcus aureus infection, and nicotine addiction ( , ). The expression of CLCF1, AEBP1, and OS9 is significantly associated with the prognosis of PTEN-wt glioma, indeed which are significantly different between high-risk score patients and low-risk patients in the PTEN-mut subgroup ( , p < 0.001).

Table 3

Clinical information about risk high and risk low patients in PTEN-mut subgroups (*** significant difference).

Group	PTEN_mut Risk score high	PTEN_mut Risk score low	p-value
Patient(number)	39	39
Age(median age)	62	57	0.09
OS.time(median days)	360	535	<0.0001***
Male	67%(26)	46%(18)
Female	33%(13)	54%(21)
Chemotherapy	33%%(13)	14%(4)
Radiation therapy	10%(4)	3%(1)
Astrocytoma	15%(6)	33%(13)
Oligoastrocytoma	0	10%%(4)
Oligodendroglioma	3%(1)	13%(5)
Glioblastoma	82%(32)	43%(17)

The Correlation Between Risk Score and Tumor Malignancy

To verify the correlation between risk score and glioma malignancy, we calculated ESTIMATEScore and TumorPurityScore for gliomas. Tumor microenvironment (TME) is closely related to the progression and prognosis of glioma (42–44). Studies have shown that TME significantly affects treatment response and clinical outcomes in cancer patients (45, 46). The components of TME are mainly resident stromal cells and immune cells, which are involved in the development of tumors (47–49). ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) is used for predicting tumor purity and the ratio of infiltrating stromal/immune cells (49). Higher ImmuneScore and StromalScore represent more immune/stromal components in TME. ESTIMATEScore, the sum of ImmuneScore and StromalScore, indicates the combined ratio of the immune and stromal components in TME (50). Ke-Wei et al. found that ESTIMATEScore was positively correlated with the survival rate at the tumor stage, suggesting that immune and stromal components were related to the invasion and metastasis of Lung adenocarcinoma (51). Glioma purity is highly correlated with major clinical and molecular characteristics, and low-purity glioma is more likely to be diagnosed as malignant. It is independently associated with reduced survival time (52). Thus, ESTIMATEScore and TumorPurityScore can bet the entry point to verify the malignant degree of glioma.

The ESTIMATE algorithm was applied for 653 glioma patients in TCGA. There was a high correlation between the risk score of patients and tumor purity score, ESTIMATEScore, StromalScore, and ImmuneScore ( ). These patients with the high-risk score in the PTEN-mut subgroup have the highest ESTIMATEScore and lowest tumor purity score, which indicates the accuracy of the risk score based on OPR-DEGs in two subgroups. Patients in each subgroup were visually represented linearly ( ). The worst survival patients with the high-risk score in the PTEN-mut subgroup were the most aggressive GBM, mostly older than 60, entirely had high ESTIMATEScore, and altogether had died at the last follow-up. On the contrary, the best survival patients with low-risk in the PTEN-wt subgroup belong to the low-grade glioma; their ages were between 20 and 60; most patients had low ESTIMATEScore, and were still alive until the last follow-up. These results demonstrate consistency between risk score and clinical survival rate.

Validation of Risk Score and Prognostic Genes in Clinical Significance

To verify the clinical value of our newly established risk model and prognostic genes, we firstly calculated and found that risk scores were higher in the PTEN-mut than in the PTEN-wt subgroup both in GBM and LGG patients ( ). We used PCA (Principal Component Analysis) to convert the dimensionality reduction of all genes into three-dimensional coordinates. Surprisingly, a similar pattern was obtained from PCA analysis by using all genes or the prognostic genes from indicated glioma groups, suggesting prognostic genes identified by our risk model reflect the variation among different types of patients ( ,  and ).

Secondly, WHO identified IDH1 mutation as an essential glioma classification criteria in 2016 (7). We divided the TCGA-glioma dataset into IDH wild-type (IDH-wt) and mutant types (IDH-mut). Only 2% of patients in the IDH-mut group had PTEN mutations (n = 7), compared with 28% in IDH-wt (n = 71, ). We then focused on the patients with IDH-wt and showed that PTEN-mut had higher risk scores than PTEN-wt groups ( ). PCA analysis showed prognostic genes identified by our risk model represent the variation of IDH -wt patients with different PTEN status, strongly indicating the accuracy of risk model ( , and ). We then concentrate on glioblastoma only or lower grades only to analyze the effect of PTEN states in IDH-wt, including 65% IDH-wt/GBM/PTEN-wt, 35% IDH-wt/GBM/PTEN-mut, 81% IDH-wt/LGG/PTEN-wt, and 19% IDH-wt/LGG/PTEN-mut ( ). Consistently, the risk score of PTEN-mut was higher than that of PTEN-wt in both IDH_wt/GBM and IDH_wt/LGG patients ( ). PCA results also supported the idea that the risk model plays a role in clinical significance evidenced by a similar pattern shown in more defined types of glioma ( , and ). To verify that previous identified prognostic genes of TCGA-glioma PTEN-mut are also the prognostic genes in IDH1-wt patients with PTEN-mut, we observed that AEBP1, CLCF1, and OS9 were significantly prognostic genes in IDH1-wt, IDH1-wt/GBM, or IDH-wt/LGG gliomas with PTEN-mut ( ). Altogether, these data demonstrate the reliability of the risk model and the potential therapeutic targets of the three genes in PTEN-mut glioma for clinical use of IDH-wt patients.

Figure 5

(A) The number of patients, survival rates, and the risk score established by our model in different types of IDH-wt glioma dataset. There are 258 IDH-wt glioma patients, 142 (55%) IDH-wt/GBM, 116 (45%) IDH-wt/LGG, 93 (65%) IDH-wt/GBM/PTEN-wt, 49 (35%) IDH-wt/GBM/PTEN-mut, 94 (81%) IDH-wt/LGG/PTEN-wt and 22 (19%) IDH-wt/LGG/PTEN-mut; (B) PCA plots of different type of IDH-wt glioma (IDH-wt/GBM/PTEN-wt, IDH-wt/GBM/PTEN-mut, IDH-wt/LGG/PTEN-wt and IDH-wt/LGG/PTEN-mut) by the expression of all genes and optimal prognostic DEGs (OPR-DEGs of PTEN-wt and PTEN-mut in and , n = 55); (C–F) Survival analysis of the expression level of AEBP1, CLCF1 and OS9 in different types of glioma with PTEN-mut, high expression and low expression of these three genes in TCGA PTEN-mut and TCGA IDH-wt/PTEN-mut were grouped by median values of expression ( C: AEBP1_high: n = 39, AEBP1_low: n = 39, CLCF1_high: n = 39, CLCF1_low: n = 39, OS9_high: n = 39; OS9_low: n = 39 in TCGA PTEN-mut; D: AEBP1_high: n = 35, AEBP1_low: n = 35, CLCF1_high: n = 35, CLCF1_low: n = 35, OS9_high: n = 35; OS9_low: n = 35 in IDH-wt/PTEN-mut), in IDH-wt/GBM/PTEN-mut and IDH-wt/LGG/PTEN-mut grouped by the best dividing points form tinyarray package of R ( (E): AEBP1_high: n = 39, AEBP1_low: n = 10, CLCF1_high: n = 8, CLCF1_low: n = 41, OS9_high: n = 14; OS9_low: n = 35 in IDH-wt/GBM/PTEN-mut; F: AEBP1_high: n = 11, AEBP1_low: n = 11, CLCF1_high: n = 10, CLCF1_low: n = 12, OS9_high: n = 14; OS9_low: n = 8 in IDH-wt/LGG/PTEN-mut).

Targeting CLCF1, AEBP1, and OS9 Genes in PTEN-mut Glioma to Explore the Mechanism and Potential Treatment

Due to the higher malignancy of the PTEN mutant glioma, we focused on the mechanism and treatment of the PTEN-mut subgroup. Subgroup-specific CLCF1, AEBP1, and OS9 are the independent OPR-DEGs and the differential expression genes between high-risk patients and low-risk patients in the PTEN-mut subgroup ( and ). Thus, they are associated with tumor progression and may be critical prognostic targets for gliomas with mutant PTEN, and more precisely in patients with IDH1-wt/PTEN-mut ( ). To explore the underlying mechanism of PTEN-mut glioma progression, the subgroup-specific and independent indicators, CLCF1, AEBP1, and OS9, are analyzed for PPI (Protein-Protein Interactions, ). In the PPI network of CLCF1 ( ), the signal transducer and activator of transcription 3 (STAT3) is upregulated in PTEN-mut but not in the PTEN-wt subgroup (PTEN-mut vs. Normal, PTEN-wt vs. Normal), showing STAT3 is correlated with PTEN-mut glioma. Similarly, in the PPI network of AEBP1 ( ), ATP/GTP binding protein 1 (AGTPBP1), tumor protein P53 target 5 (TP53TG5), leucine-rich repeat transmembrane neuronal 4 (LRRTM4), transcription factor 3 (TCF3) and aspartoacylase (ASPA) are DEGs between PTEN-mut subgroup and Normal patients and may be associated with poor prognosis of PTEN-mut subgroup. In the PPI network of OS9 ( ), OS9 and derlin 2 (DERL2) involve in protein processing in the endoplasmic reticulum and complicate CFTR (Cystic fibrosis transmembrane conductance regulator) causes cystic fibrosis pathways, which may be associated with the prognosis of glioma with mutant PTEN.

For the targeted treatment of PTEN-mut glioma, we constructed the miRNA-gene ( ), TFs-genes ( ), and protein-chemical network ( ). CLCF1, AEBP1, and OS9 have different binding miRNA and TFs, respectively. Notably, the small compound (+)-JQ1 can bind to all these three proteins (CLCF1, AEBP1, and OS9). To verify that (+)-JQ1 binds with CLCF1, AEBP1, and OS9, we predicted the structure of CLCF1, AEBP1, and OS9, and docked the binding sites of (+)-JQ1 with these three proteins. Docking results showed that (+)-JQ1 binds to all three proteins stably ( ). Furthermore, we treated U251(PTEN-deficient) and U343 (PTEN-wt) cells with different concentrations of (+)-JQ1 ( ). Results showed that U251 cells were more sensitive to (+)-JQ1. This experimental data indicated that (+)-JQ1 indeed had a potential medicinal effect on the PTEN-mut subgroup. Therefore, (+)-JQ1 potentially competes for and binds to the CLCF1, AEBP1, and OS9 proteins, resulting in the interruption of the pro-oncogenic pathway in PTNE-mut gliomas. These results highlight the therapeutic potential of (+)-JQ1 in treating PTEN-mut gliomas and demand extensive in vitro and in vivo studies for experimental verifications.

Figure 6

Targeted treatment for PTEN-mut glioma: (A) The miRNA-gene network of CLCF1, AEBP1 and OS9 by NetworkAnalyst (Materials and methods 2.6); (B) The TFs-gene network of CLCF1, AEBP1 and OS9 by NetworkAnalyst (Materials and methods 2.6); (C) The protein-chemical network of CLCF1, AEBP1, and OS9 by NetworkAnalyst (Materials and methods 2.6); (D–F) The docking site between CLCF1, AEBP1, OS9 with (+)-JQ1 by Autodock (Materials and methods 2.7); (G) Colony formation assay was performed in U251(PTEN- deficiency) cells treated by (+)-JQ1 with different concentration (Materials and methods 2.8); (H) Colony formation assay was performed in U343 cells treated by (+)-JQ1 with different concentration (Materials and Methods 2.8).

Discussion

Gliomas are the most lethal primary brain tumors lacking effective treatment, and the PTEN mutation event significantly reduces the survival rates of glioma patients. Depending on the PTEN status, glioma can be divided into two apparent subgroups: the PTEN-wt and PTEN-mut. In this study, we established prognostic risk models for the two subgroups, respectively. In summary, we first analyzed the expression of subgroup-specific tumor DEGs and the enriched pathways. Based on the subgroup-specific DEGs, we successfully screened subgroup-specific optimal prognostic signature (OPR-DEGs) and calculated risk scores for individual patients based on their OPR-DEGs. CLCF1, AEPB1, and OS9 in the PTEN-mut subgroup are the optimal and independent prognostic signatures and can be used as potential targets for the diagnosis, prediction, and treatment of PTEN-mut glioma. Finally, we found miRNAs and TFs that interacted with CLCF1, AEPB1, and OS9 genes and chemicals that interacted with CLCF1, AEPB1, and OS9 proteins. We showed the therapeutic potential of the gene or protein interacting agents, especially (+)-JQ1.

Mutational loss of PTEN is an established malignancy event and PTEN significantly influences therapeutic efficacy in glioma (23, 25, 27, 28). PTEN status of glioma patients according to their unique signatures improve the survival rate. Chen et al. reported symbiotic macrophage-glioma cell interactions to cause synthetic lethality in PTEN-Null glioma (53). In 2006, Parsa et al. had shown tumor-specific T cells more effectively target human glioma expressing wild-type PTEN than those expressing mutant PTEN (54). Conversely, we demonstrated that the oncogenic role of Smurf1 promotes GBM growth by mediating PTEN ubiquitylation and degradation (55). Additionally, we also realized that PTEN status is a prognostic marker in all grades of gliomas and LGG only, but not in the GBM group only. These may due to PTEN defect is an initial cause of malignancy in glioma patients. The progression of GBM leads to acquired secondary mutations (6). In this rationale, we combined the RNA-seq data of GBM and LGG to build a risk model. To verify the clinical significance of our risk model, we focused on IDH-wt glioma first, then concentrate on GBM only or LGG only. We identified that PTEN status is also a critical prognostic factor in terms of clinical classification of patients, such as in IDH1-wt, IDH1-wt/GBM, and IDH1-wt/LGG ( ).

In the PTEN-wt subgroup, 44 OPR-DEGs were fitted with the best survival time. Among the 44 OPR-DEG genes, 14 genes were independent prognostic signatures (Multivariable Cox, p < 0.05, ). SSTR5 has the cytostatic effects of somatostatin in C6 glioma cells by activating PTPeta (protein tyrosine phosphatase eta) and inhibiting extracellular signal-regulated kinase (ERK)1/2 activity (56). HMGN5, as a potential oncogene, is highly expressed in breast cancer and hormone-induced mouse uterine adenocarcinoma. Qu et al. found suppression of HMGN5 caused a delay in the growth of human glioma cells (57). MEGF10, a critical IDH mutation predictor, plays an important role in cell migration, apoptosis, and proliferation (58). SSTR5, HMGN5, and MEGF10 have been associated with glioma, but they are not associated with PTEN status. HSPC159 induces epithelial-mesenchymal transition, activates the PI3K/protein kinase B (Akt) pathway, and promotes proliferation and metastasis in breast cancer (59). Studies reported that upregulated SLC6A6 induces tumorigenesis and reduces clinical outcomes in gastric cancer (60). Pei et al. found that loss of LAMC2 changes epithelial-mesenchymal transition and inhibits angiogenesis in cholangiocarcinoma via inactivation of the epidermal growth factor receptor (EGFR) signaling pathway (61). LGR6, as a tumor suppressor, belongs to the rhodopsin-like seven-transmembrane domain receptor superfamily and is a high-affinity receptor of R-spondins with potential functions (62). These genes have been reported as markers for other tumors but have not been studied in gliomas and may be potential targets for gliomas, especially the PTEN-wt subgroup.

In the PTEN-mut subgroup, 11 OPR-DEGs were fitted with the best survival time. Out of 11, three OPR-DEGs (CLCF1, AEBP1, and OS9) are independent prognostic signatures and are upregulated in PTEN-mut high-risk patients ( and ). CLCF1 is a neurotrophic and B cell-stimulating factor belonging to the interleukin-6 (IL-6) family. Studies have found that CLCF1 induced phosphorylation of signal transducer and activator of transcription 3 (STAT3) ( ) (63). STAT3 is upregulated in PTEN-mut but not in the PTEN-wt subgroup (PTEN-mut vs. Normal, PTEN-wt vs. Normal), showing STAT3 is correlated with PTEN-mut glioma ( ). Notably, STAT3 participates in three signaling pathways, including IL-6 family signaling, Janus kinase (JAK)/STAT signaling pathway, and EGFR/RTK pathway. Continuous activation of STAT3 in GBM with PTEN loss can induce cell proliferation, anti-apoptosis, maintenance of glioma stem cells, tumor invasion, angiogenesis, and immune evasion (64, 65). Inhibition of CLCF1 to reduce phosphorylation of STAT3 may be an effective strategy for treating glioma with mutant PTEN. AEBP1 is involved in the regulation of adipogenesis, mammary gland development, inflammation, macrophage cholesterol homeostasis, and atherogenesis in various human tumors (66–70). Swati Sinha et al. reported that the downregulation of AEBP1 in PTEN-deficient cells activated cell death through a caspase-independent pathway that is different from PTEN-wt glioma ( ) (66, 67). Our study further confirms the importance and potential therapeutic intervention implications of CLCF1 and AEBP1 in PTEN mutant glioma. We also predicted that OS9 is the best independent prognostic signature of the PTEN-mut subgroup. The protein encoded by OS9 is highly expressed in osteosarcoma and binds to hypoxia-inducible factor 1 (HIF-1), a key regulator of hypoxia response and angiogenesis (71–73). Studies on OS9 in glioma have not been reported, but our model demonstrated the strong correlation between OS9 and survival time, indicating that OS9 is also a potential target for glioma. We analyzed the gene-miRNA, gene-TF, and protein-chemical interaction networks to guide the targeted therapy ( ). (+)-JQ1 compound is a triazole-diazepine compound family member, which functions as a pan-BET (bromodomain and extra-terminal motif) family inhibitor (74). (+)-JQ1 is known to suppress cell proliferation and can be used as a therapeutic drug for many cancers, including multiple myeloma and acute myeloid leukemia (75). Following a single intraperitoneal dose of (+)-JQ1 (50 mg/kg) in male mice, Matzuk et al. measured the (+)-JQ1 concentration in serum, testis, and brain. They observed nearly uniform blood-brain barrier permeability after single-dose pharmacokinetic studies (AUCbrain/AUCplasma = 98%) (76). Korb et al. in 2015 examined whether JQ1 affects brain function in mice. They injected wild-type adult male mice with (+)-JQ1 (50 mg/kg) daily for 1 week or 3 weeks before performing behavioral tests, and the result proved that (+)-JQ1 has excellent blood-brain barrier permeability (77). (+)-JQ1 compound may combine and inhibit the CLCF1, AEBP1, and OS9 interacting proteins, and hence, it can be a potential treatment option for PTEN-mut gliomas.

Conclusion

In conclusion, we established a prognostic risk model and risk score in glioma with different PTEN status and obtained 14 independent prognostic signatures in PTEN-wt glioma and 3 independent prognostic signatures in PTEN-mut glioma. Thus, in the PTEN-mut glioma, CLCF1, AEBP1, and OS9, which are significantly associated with survival time, may induce glioma progression and are critical targets for diagnosis, prognosis prediction, and treatment, giving therapeutic recommendations to glioma with mutant PTEN.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding authors.

Author Contributions

PZ contributed to the conception, design, and drafting of the manuscript. QX provided useful comments and suggestions. LD and QX revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the Beijing Natural Science Foundation (Z190018), the National Natural Science Foundation of China (81870123), the China Postdoctoral Science Foundation Grant (2018M641206), the National Science Foundation for Young Scientists of China (81902545), Beijing Institute of Technology Research Fund Program for Young Scholars (XSQD-202110002).

Conflict of Interest

Authors LD, PZ and QX applied for a patent in China, the patent application number is 202010977490.6.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.