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)(9)(10)(11)(12)(13)(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)(17)(18)(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 PI3Ka 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 L1penalized (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)(30)(31)(32)(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.

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)

Differentially Expressed Genes
Patients were divided into two glioma subgroups according to the status of PTEN (Figure 1), including PTEN-wt and PTENmut patients, whose clinical features are respectively presented in Table 1. 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 ( Figure 1) 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).
prognostic DEGs (PR-DEGs) that were significantly linked to survival time by calculating the log-rank test and univariate Cox regression ( Figure 1). 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 ( Figure 1). 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 ( Figure 1). The risk score of every patient was established by the following formula (36):  Expression Risk Score = o bRNAn Â ExpRNAn b RNAn and Exp RNAn 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 (Figure 1).

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 ( Figure 1). 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 ( Figure 1).

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 proteinchemical 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 ( Figure 2A). 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 (Supplementary Figures 1A-D and Figure 2B). As expected, the CGGA dataset showed worse survival in PTEN-mut gliomas compared to the PTEN-wt subgroup ( Figure 2C). To identify survival differences among mutated PTEN glioma patients, we extracted the PTEN mutation sites list (Supplementary Table 1) 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 (Supplementary Figure 1E). 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).
The clinical information of the PTEN-wt and the PTEN-mut subgroups was statistically analyzed ( Table 1). 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 ( Table 1). 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 (Table 1). 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 ( Figure 2D). To identify the subgroup-specific DEGs of the PTEN-wt and PTEN-mut subgroups, the DEGs-all, DEGswt, and DEGs-mut (see Materials and methods) were calculated ( Figure 1). Notably, the pathway analysis of DEGs ( Figures 2E-G) showed the enriched mitogen-activated protein kinase (MAPK) signaling pathway in the PTEN-mut subgroup ( Figure 2G). 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 ( Figure 2H Table 3), having significant enrichment in the extracellular region of Cellular Component term, which is closely related to tumor progression ( Figure 2J, 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, Supplementary Table 5). The AUC value of the L1-penalized (LASSO) Cox-pH regression model was equal to 0.8807 (Supplementary Figure 2A), 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 ( Figure 3A). 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 ( Figure 3A). Finally, we calculated the nomogram by clinical features and the risk scores to predict clinical survival ( Figure 3B). As shown in calibration plots by internal validation of Bootstrap Resamples (1000) Table 6). 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 ( Figure 3E, Supplementary Table 6). 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 ( Figure 3F, Supplementary Table 6). The clinical features of a high-risk and low-risk patient are shown in Table 2. We found the expression level of PTEN in high-risk patients was lower than lowrisk patients in the PTEN-wt subgroup (p < 0.001), which showed that the risk score was also correlated with the level of PTEN expression ( Figure 3G).
OPR-DEGs, which best fit the patient's survival time, are potential targets for driving glioma progression.  Figure 2E). 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 (Figures 3H-L, 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 (Supplementary Table 7). The AUC value was equal to 0.9414, indicating the model's fine applicability (Supplementary Figure 3A). 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 ( Figure 4A). The 11 OPR-DEGs in the PTEN-mut subgroup have three independent OPR-DEGs (Multivariate Cox: p < 0.05), including CLCF1, AEBP1, and OS9 (Supplementary Table 7). Meanwhile, the nomogram for the PTEN-mut subgroup was also constructed to predict individual patient survival rates ( Figure 4B). The calibration plots by internal validation of Bootstrap Resamples (1000) in the TCGA-PTEN-mut subgroup demonstrate the high availability of predictive models (Supplementary Figure 3B, C-index of internal validation = 0.822). As the CGGA-PTEN-mut subgroup has only nine patients, external validation could not be performed.
As we predicted, the prognosis of low-risk patients was significantly better than that of high-risk patients (Figures 4C, D; p < 0.0001). The clinical features of the high-risk and low-risk patients are shown in Table 3. There are 273 upregulated and 203 down-regulated DEGs between the high-risk patients and the lowrisk patients in the TCGA-PTEN-mut subgroup (Supplementary Figure 3C, selection criteria: FDR < 0.05 and |log2FC| > 1, Supplementary Table 8). The heatmap presents DEGs' expression between the high-risk patients and the low-risk patients ( Figure 4E, Supplementary Table 6). 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 ( Figure 4F, Supplementary Table 8). 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 lowrisk patients in the PTEN-mut subgroup (Figures 4G-I, p < 0.001).

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)(43)(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)(48)(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 ( Supplementary Figures 4A-D). 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 (Supplementary Figure 4E). 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 (Supplementary Figure 5A). 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 (Supplementary Tables 5, 7 and Supplementary Figure 5B).
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, Supplementary Figure 5C). We then focused on the patients with IDH-wt and showed that PTEN-mut had higher risk scores than PTEN-wt groups (Supplementary Figure 5C). 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 (Supplementary Tables 5, 7 and Supplementary Figure 5D). 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/PTENwt, 35% IDH-wt/GBM/PTEN-mut, 81% IDH-wt/LGG/PTEN-wt, and 19% IDH-wt/LGG/PTEN-mut ( Figure 5A). Consistently, the risk score of PTEN-mut was higher than that of PTEN-wt in both IDH_wt/GBM and IDH_wt/LGG patients ( Figure 5A). 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 (Supplementary Tables 5, 7 and Figure 5B). To verify that previous identified prognostic genes of TCGA-glioma PTEN-mut are also the prognostic genes in IDH1wt 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 ( Figures 5C-F). Altogether, these data demonstrate the reliability of the risk model and the potential therapeutic targets of the three genes in PTENmut glioma for clinical use of IDH-wt patients.
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 PTENmut subgroup ( Figures 4G-I and Supplementary Table 7). 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 ( Figures 5C-F).
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   Table 7, risk score in Supplementary  Table 1) 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, Supplementary Table 8) 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, Supplementary Table 8 Figure 6C), 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.

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 ( Figure 5).
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, Supplementary Table 5). 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 seventransmembrane 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 ( Figures 4G-I and Supplementary Table 7). CLCF1 is a neurotrophic and B cellstimulating 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) (Supplementary Figure 6D) (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 (Supplementary Figure 6A). 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, antiapoptosis, 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)(67)(68)(69)(70). Swati Sinha et al. reported that the downregulation of AEBP1 in PTENdeficient cells activated cell death through a caspase-independent pathway that is different from PTEN-wt glioma (Supplementary Figure 6E) (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)(72)(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 ( Figures  6A-C). (+)-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/Supplementary Material. 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.