To retrieve the transcriptomic expression levels of SOX9 in GBM, we utilized the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) database. Additionally, to validate protein-level expression, we performed western blotting using GBM tumor tissues and adjacent normal brain tissues collected from clinical samples. The pan-cancer RNA-seq data were obtained from TCGA (https://portal.gdc.cancer.gov/) and GTEx (https://gtexportal.org/). HTSeq-FPKM and HTSeq-Count data of the GBM samples were acquired from the TCGA (https://portal.gdc.cancer.gov/repository) for further analysis. The study fully complies with the guidelines provided by TCGA and GTEx.

LinkedOmics (http://www.linkedomics.org/) (17) was used to assess and draw a heatmap of the top 35 positively/negatively related genes with SOX9. Related genes with the adjusted P-value <0.05 were applied for functional enrichment analysis. We used Metascape (https://metascape.org) (18) to visualize the enriched biological process (BP), cellular composition (CC), molecular function (MF) and KEGG pathway terms of the SOX family and its co-expressed genes.

The DESeq2 R package was adopted to compare expression data of low- and high-expression of COMMD7 (cut-off value of 50%) in AML samples (HTseq-Count) to identify DEGs (19). R package ggplot2 (3.3.6) was used to visualize volcano plots of the result.

The PPI network of DEGs was predicted using the Search Tool for the Retrieval of Interacting Genes (STRING) database (20). The interaction score threshold of 0.4 was set as the cut-off criterion. The PPI network was mapped using Cytoscape (version 3.7.1) (21), and the most significant modules in the PPI network were identified using MCODE (version 1.6.1) (22). Selection criteria were as follows: MCODE scores >5, degree cut-off = 2, node score cut- off = 0.2, Max depth = 100, and k-score = 2.

DEGs with the threshold for | log fold change (logFC)| >2 and adjusted P-value (adj P-value) <0.05 were applied for functional enrichment analysis. Gene Ontology (GO) functional analysis comprising cellular component (CC), molecular function (MF), and biological process (BP), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were implemented using the ClusteProfiler package in R (23).

R package ClusteProfiler (3.14.3) was used for GSEA to elucidate the functional and pathway differences between the high- and low-expression groups of SOX9 (23). The gene set was permutated 1,000 times for each analysis. Adjusted P-value < 0.05 and FDR q-value < 0.25 were considered to be statistically significant.

SOX9 expression was compared between tumors and normal tissues by receiver operating characteristic (ROC) analysis to test the predictive value of SOX9 for GBM diagnosis. Clinicopathological characteristics were compared for high- and low- SOX9 expression groups using the Wilcoxon rank sum test (continuous variables) or Spearman chi-square test (rank variables). The correlation between GBM expression and clinicopathological characteristics was evaluated by logistic analysis. Kaplan–Meier (K-M) analysis, univariate, and multivariate Cox regression analysis were employed for prognosis analysis.

In order to individualize the prediction of overall survival (OS) in THCA patients, a nomogram was generated using the RMS R package (version 5.1-3), which included genes screened by LASSO coefficient filtering, prominent clinical characteristics and calibration plots. The calibration curves were evaluated graphically by mapping the nomogram-predicted probabilities against the observed rates, and the 45°line represented the best predictive values. Concordance index (C-index) was used to determine the discrimination of the nomogram, and the bootstrap approach was used to calculate 1000 resamples. In addition, C-index and receiver operating characteristic (ROC) were used to analyze and compare the predictive accuracy of the nomogram and separate prognostic factors. All statistical tests were double-tailed with 0.05 as the statistical significance level.

The ssGSEA package and ESTIMATE package in the GSVA package [version 1.34.0] were used for immuno-infiltration correlation analysis of SOX9. The statistical significance of the difference was evaluated by Spearman’s test. In addition, Wilcoxon rank sum test was used to analyze the correlation between SOX9 expression and immune checkpoints expression in GBM. The correlation between SOX9 expression and immune checkpoints expression was evaluated by Spearman chi-square test.

We used the GSCALite database (http://bioinfo.life.hust.edu.cn/web/GSCALite/) (24) and TISDIB database (http://cis.hku.hk/TISIDB/index.php) to analyze the correlation between SOX9 expression and sensitivity to current chemotherapeutic or targeted drugs for GBM.

Tumor and brain tissues were lysed by 1% SDS lysis buffer and boiled at 95°C for 10 min, then centrifuged at 12,000 rpm for 10 min at room temperature. The protein concentration was determined by a BCA protein assay reagent kit (Beyotime Institute of Biotechnology, Jiangsu, China). The protein (35 µg) was separated by SDS-PAGE gels and transferred onto PVDF membranes (Merck, Darmstadt, Germany). The membranes were blocked for 2h by 5% fat free milk in TBST at room temperature, followed by incubation overnight at 4°C with the primary antibodies for anti-SOX9 (ABclonal, Wu Han, China; Cat# A19710) and anti-GAPDH (Cell Signaling Technology, Danvers, MA, USA; Cat# 2118). The membranes were washed with TBST and incubated for 1h with the corresponding HRP-conjugated second antibodies. Bands were visualized with ECL Reagents (Beyotime Institute of Biotechnology, Jiangsu, China; P0018S). Each Western blot experiment was independently repeated four times (N = 4) using patient-derived GBM and adjacent normal brain tissues.

In this study, R (4.2.1) and corresponding R packages were utilized for statistical analysis. In all tests, P value < 0.05 was considered statistically significant.

SOX9 is highly expressed in many tissues in the human body, including the central nervous system ( Figure 1A ). RNA-seq data was downloaded in TCGA and GTEx formats processed uniformly through the toil process. By comparing the expression of SOX9 normal samples in TCGA/GTEx with the corresponding tumor samples in TCGA, SOX9 was found significantly high expressed in 14 types of cancer ( Figure 1B ), including GBM ( Figure 1C ). The Western blot experiment has confirmed that the expression of SOX9 is significantly higher in GBM compared to adjacent normal tissues ( Figure 1D ). Quantitative analysis confirmed that SOX9 protein levels were significantly elevated in GBM tissues compared to adjacent normal tissues ( Figure 1E ).

Co-expressed genes of the SOX9 in GBM were obtained from the LinkedOmics database, and heatmaps of the top 35 positively/negatively related genes were drawn ( Figures 2A, B ). Then, we analyzed the GO and KEGG pathway terms of these genes using the Metascape database. The results of positively related gens showed significant enrichment of the biological processes (BP) terms “cellular process” and some processes related to “growth and development” ( Figure 2C ). The results showed significant enrichment of top 20 clusters terms “chromatin organization”, “mRNA metabolic process”, “DNA damage response”, “mitotic cell cycle”, “regulation of cell cycle process” and “regulation of cellular response to stress” ( Figure 2D ). The results of negatively related gens were shown in Figures 2E, F .

The high- and low-expression groups’ gene expression profiles were analyzed for differences in the median mRNA expression. A total of 731 DEGs, including 38 upregulated and 693 down-regulated, were identified statistically significant between SOX9 high- and low-expressed groups (|log FC| >2 and adj P <0.05) ( Figure 3A ). The network of TREM1-related DEGs was constructed by STRING, with a threshold of 0.4. The PPI network was displayed by Cytoscape-MCODE ( Figure 3B ). The results of up-regulated gens showed significant enrichment of terms “cell fate specification” ( Figure 3C ). The results of down-regulated gens showed significant enrichment of terms “complement activation”, “cytokine−cytokine receptor interaction” and “cytokine receptor” et al. ( Figure 3D ). To further understand the biologic pathways involved in GBM with different SOX9 expression levels, GSEA was performed between low- and high-SOX9 expression datasets to identify critical signaling pathways. Significant differences (FDR <0.05, adj P <0.05) were observed in the enrichment of MSigDB Collection (C2.all.v7.0.symbols.gmt) of these pathways. G2/M Checkpoint was enriched in SOX9 high-expression phenotype ( Figure 3E ). In the low expression of SOX9 phenotypes, epithelial mesenchymal transition, inflammatory response, coagulation, TNF-α signaling via NF-κB, KRAS signaling up, and IL-6-JAK STAT3 signaling presented significantly enriched ( Figure 3F ).

The main clinical characteristics of GBM in TCGA was shown in Table 1 . In this study, a total of 168 cases (59 females and 109 males) were analyzed. The median SOX9 expression [log2(TPM+1)] was regarded as the cut-off value. Correlation analysis ( Table 2 ) suggested that SOX9 expression was significantly correlated with IDH status (P <0.05). Logistic analysis was applied to further verify the relationship between GBM the factor of IDH status and the SOX9 high-low dichotomy. As a result, high expression of SOX9 showed a significant positive correlation with IDH status (odds ratio [OR], 0.168; P <0.05). Furthermore, the potential value of SOX9 in GBM patients was examined by ROC curve analysis, with the AUC of 0.916, revealing that SOX9 was a potential biomarker ( Figure 4A ). Besides, the Wilcoxon Rank SUM test was used to compare the expression of SOX9 in patients with different IDH status features. SOX9 expression was significantly higher in IDH-wildtype GBM compared to IDH-mutant samples (P < 0.05, Wilcoxon rank-sum test), as shown in Figure 4B . The relationship between SOX9 expression and prognosis was analyzed in GBM patients by using Kaplan-Meier. Patients with high expression of SOX9 had a worse prognosis than those with low SOX9 expression (over survival [OS], hazard ratio [HR]=1.38 (0.97 – 1.95), P = 0.070; disease specific survival [DSS], HR=1.24(0.86 -1.79), P=0.252, progress free interval [PFI], HR= 1.26 (0.89 − 1.78), P = 0.202) ( Figures 4C–E ). In order to further explore the accuracy of SOX9 in the evaluation of the prognosis of GBM patients, we conducted ROC curve analysis in both the training cohort and the validation cohort. As shown in Figure 4F , in the TCGA cohort, the area under the curve (AUC) values at 1, 3 and 5 years were 0.514, 0.440, and 0.857, respectively. In addition, the Kaplan-Meier analysis was used to generate survival curves for SOX9 expression in the IDH-mutant subtype. Patients with high expression of SOX9 had a worse prognosis than those with low SOX9 expression: IDH status: WT (OS, HR = 1.00 (0.70 − 1.43), P = 0.989), IDH status: Mut (OS, HR = 1.60 (0.35 − 7.37),P = 0.543) ( Figures 4G, H ).

Screening was performed by LASSO coefficients to select non-zero variables that satisfied the coefficients of lambda. min, and 4 genes related with SOX9 were screened out. The risk score is calculated as follows: risk score = (0.09*EREG expression level - 0.027*TPTEP1 expression level - 0.309*OR4K2 expression level - 0.329*TSPY2 expression level) ( Figures 5A–B ). Multifactorial Cox regression analysis of these genes was performed using the survival R package and presented as forest plots ( Figure 5C ). Kaplan-Meier curves showed that OR4K2 high expression is a positive prognostic indicator ( Figure 5D ). The predictive power of the OS risk score assessed by the ROC curve over time was 0.391 AUC at 1 year, 0.457 at 3 years and 0.365 at 5 years with OR4K2 high expression ( Figure 5E ). A nomogram was constructed based on the Cox regression analysis results using the RMS R package ( Figure 5F ). IDH status, SOX9 expression and OR4K2 expression were included in the model. The points of each variable were accumulated and recorded as the total points. The probability of GBM patient survival at 1-, 3-, and 5-year was determined by drawing a line from the total point axis straight down to the outcome axis. The 1-year survival probability was determined by drawing a vertical line downward on the total point axis along the 140-direction ending axis, suggesting the probability of 1-year survival < 50%, both of the probability of 3- and 5-year < 5%. The prediction results of the nomogram calibration curve of OS were consistent with all patients’ observation results ( Figure 5G ).

Spearman correlation analysis showed that the expression level of TREM1 in the THCA microenvironment was correlated with the immune cell infiltration level quantified by SSGSEA. Specifically, SOX9 was positively associated with NK cells, Tcm, NK CD56bright and Tgd, but T cells inversely ( Figures 6A–F ).

PD1/PD-L1 and other immunization checkpoints are key immunization checkpoints that are responsible for tumor immune escape. Considering the potential oncogenic role of SOX9 in GBM, the relationship of SOX9 with PDCD1 and other immunization checkpoints was assessed. This correlation is visualized in the overall expression heatmaps and scatter plots ( Figures 7A, B ). The results showed that MSMO1 was negatively correlated with the expression of CD274 (P<0.05) ( Figure 7C ). The results showed that MSMO1 was negatively correlated with the expression of HAVCR2, PDCD1 and TIGIT (P <0.05) ( Figures 7D–F ). These results demonstrate that tumor immune escape might be involved in SOX9 mediated carcinogenesis of GBM. Finally, we used GSCALite online tool and TISDB database to analyze the relationship between the expression of SOX9 and sensitivity to current therapies. The result indicated that the expression levels of SOX9 were positively correlated with sensitivity to the most current cancer-targeted drugs or chemotherapy drugs ( Figure 7G ). Thus, SOX9 could represent a new target for predicting drug sensitivity and for developing multitarget combined therapy.