The expression profile and corresponding clinical data of OS patients were downloaded from the TARGET database (22).¹ The dataset includes 85 OS patients with survival information, which were then randomly separated into training (n = 43) and testing (n = 42) cohorts at cut-off 5:5. Moreover, the GSE39055 dataset which contains 36 OS patients with survival information were downloaded from the GEO database² as a validation cohort. Furthermore, we extracted 8 glycometabolic gene sets from the molecular signatures database (MSigDB) (23).³ As shown in Supplementary Table 1, the entire glycometabolic gene set contained 291 genes after removing overlapping genes. There were 282 genes left after removing genes whose expression is 0 in 50% of the samples. The workflow chart of this study is shown in Figure 1.

Univariate Cox analysis was performed to assess the prognostic value of glycometabolic genes. Genes with a P-value < 0.05 were considered as potential prognostic genes. Next, the LASSO regression algorithm was used to construct an optimal glycometabolic gene signature. The risk score for each patient was calculated as follows: risk score = ∑i=1nβi×Expgene(i), where n is the number of genes in this prognosis model, β represents the regression coefficient, and Expgene is the expression level of each gene. We divided the patients from the training cohort (n = 43) into high- and low-risk groups based on the median value of the risk score, that individuals in the high-risk groups suffered from a lower survival probability and higher risk of death compared to that in low-risk group. Differences of survival probability between these two risk groups were assessed using Kaplan–Meier (KM) curves. Receiver operating characteristic (ROC) curves were used to evaluate the prognostic capacity of the gene signature. The gene signature was further tested and validated in the testing (n = 42) and validation cohort (n = 36), respectively.

The univariate and multivariate Cox regression analyses were used to identify independent prognostic factors in the TARGET dataset. The correlation between risk score and metastasis was evaluated using chi-squared test. The results of multivariate regression analysis were used to establish a nomogram to predict the 1-, 3-, and 5-year overall survival. Moreover, the discrimination and accuracy of the nomogram were assessed by the concordance index (C-index), calibration curves, and decision curves, respectively.

The “limma” package was used to analyze differentially expressed genes (DEGs) between the high- and low-risk groups in the TARGET dataset. Gene set enrichment analysis (GSEA) was performed to determine the biological functions of DEGs. Gene Ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses were performed using the ClusterProfiler package (24, 25). Moreover, we also explored the single-sample GSEA (ssGSEA) scores between the high- and low-risk groups based on gene expression profiles involved in the top 10 GO and KEGG terms.

To analyze the differences in the proportion of 28 immunocytes between the high- and low-risk groups, ssGSEA was performed using TARGET dataset (26). We also investigated the correlation between the proportion of the 28 immunocytes and prognostic genes. Tumor immune infiltration scores including immune score, stromal score, ESTIMATE score, and tumorpurity were assessed using the “estimate” and “limma” packages.

The R package “pRRophetic” was used to calculate the half-maximal inhibitory concentration (IC50) of chemotherapy drugs, and the Wilcoxon signed-rank test was employed to compare the differences of IC50 between the high- and low-risk groups.

To explore the interactions among these prognostic genes, a PPI network was constructed using the online tool search tool for the retrieval of interacting genes (STRING) database (27). Furthermore, a plug-in of Cytoscape software (version 1.6.20), molecular complex detection (MCODE), was used to screen the significant modules in the PPI network (28).

MiRNAs that can regulate prognostic genes were predicted based on miRanda software. Then, lncRNAs that can regulate the predicted miRNAs were predicted through miRanda software. To improve the accuracy of the competing endogenous RNA (ceRNA) network, we further screened the results using the following criteria: combined score > 200 was set as the screening criteria of lncRNA-miRNA interactions, combined score > 200 and minimum free energy (MFE) score < −200 were set as the screening criteria of miRNA-gene interactions. Finally, a ceRNA network was constructed using Cytoscape.

This study was approved by the Ethics Committee of Affiliated Hospital of Guizhou Medical University (No. GZYD003-201753035). The tumor samples and corresponding adjacent normal tissues were obtained from the patients who were diagnosed with OS. All the clinical information of patients were shown in Supplementary Table 2. Human samples were fixed with 4% paraformaldehyde, embedded in paraffin, and sliced into 5-μm sections. Sections were deparaffinized in xylene, hydrated with a graded ethanol series at room temperature. Next, sections were treated with 3% H₂O₂ to block endogenous peroxidase activity, and then blocked with 5% bull serum albumin (BSA) for 30 min at room temperature. The specimens were incubated with the appropriate primary antibodies at 4 °C overnight and incubated with goat anti-rabbit secondary antibodies (1:200, Proteintech, Wuhan, China, SA00001-2) at 37 °C for 2 h at room temperature. Primary antibodies used in this experiment included antibodies against PRKACB (1:200, Proteintech, 12232-1-AP), SEPHS2 (1:200, Proteintech, 14109-1-AP), GPX7 (1:200, Proteintech, 13501-1-AP), and PFKFB3 (1:200, Proteintech, 13763-1-AP). The DAB substrate system (Solarbio, China) was used for color development, and hematoxylin staining was used to reveal the cell nuclei. Images were obtained under a light microscope.

All the above analyses were completed by the R software (Version 4.0.5). A time-dependent ROC analysis was performed by the “pROC” package (29). Chi-squared test was used to explore the correlation between risk score and metastasis of OS. And a nomogram was constructed by the “rms” packages. Results with a P-value < 0.05 were considered as statistical significance.

In order to establish a glycometabolic gene signature, we identified 19 glycometabolic genes which were significantly associated with the overall survival of the patients with OS in the training cohort using univariate Cox regression analysis (Figure 2A). Subsequently, LASSO Cox regression analysis was performed to screen candidate genes for constructing the gene signature. The results indicated that 4 glycometabolic genes, including protein kinase CAMP-activated catalytic subunit beta (PRKACB), selenophosphate synthetase 2 (SEPHS2), glutathione peroxidase 7 (GPX7), and 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 (PFKFB3), were identified according to the optimum penalty parameter (λ) value (Figures 2B, C). The coefficient of each candidate gene is shown in Figure 2D. Combining the gene expression with corresponding regression coefficient, a risk score model was then established using a formula as follows: risk score = expression level of PRKACB × (−0.10894429) + expression level of SEPHS2 × (−0.01563243) + expression level of GPX7 × 0.32559346 + expression level of PFKFB3 × 0.55009921. According to the median value of the risk scores, the patients in the training cohort were divided into high- (n = 21) and low-risk group (n = 22). The expression of GPX7 and PFKFB3 had a positive correlation with risk scores, while the expression of SEPHS2 and PRKACB was negatively associated with risk scores (Figure 2E). However, the expression of these four genes was not associated with age and gender. As shown in Figure 2F, the patients with high-risk scores seemed to have high mortality rates and shorter survival time than those with low-risk scores. Similarly, the KM curves showed that patients in the low-risk group had a better clinical prognosis than those in the high-risk group (Figure 2G). Furthermore, the respective area under the ROC curve (AUC) of 1-, 2-, 3-, 4- and 5-year survival was 0.904, 0.939, 0.900, 0.900, and 0.919, suggesting that this glycometabolic gene signature showed a favorable prognostic value of overall survival (Figure 2H).

According to the median value of the risk scores, the 42 patients in the testing cohort were equally separated into high- and low-risk groups. The heatmap visualized the expression of GPX7, PFKFB3, SEPHS2, and PRKACB in high- and low-risk groups. The results showed that the expression patterns of these 4 prognostic genes in testing cohort were similar to that in training cohort (Figure 3A). The risk scatter plots and KM survival analysis also indicated that the patients with high-risk scores had worse prognoses than those with low-risk scores (Figures 3B, C). The AUC was 0.687 at 1-year, 0.730 at 2-year, 0.729 at 3-year, 0.729 at 4-year, and 0.718 at 5-year (Figure 3D). In order to verify the robustness of this prognostic signature, GSE39055 which contains 36 OS patients was used as an external validation cohort. The risk score of each patient in validation cohort was calculated using the same formula determined in training cohort. The 36 patients in the validation cohort were equally separated into high- and low-risk groups. The expression patterns of these 4 prognostic genes were depicted on the heatmap which also showed similarity to that in training cohort (Figure 3E). As shown in Figures 3F, G, poor prognosis was found in patients with high-risk scores, and the percentages of dead patients were 50% (9/18) and 5.6% (1/18) in the high- and low-risk groups, respectively. The respective AUC was 0.907, 0.825, 0.743, 0.743 and 0.743 for 1, 2-, 3-, 4-, and 5-years survival (Figure 3H). All these results indicated that this glycometabolic gene signature had a robust performance in predicting the prognosis of patients with OS.

We firstly analyzed the correlation between risk score and metastasis. The results revealed that the proportion of individuals with metastasis in the high-risk group was more than that in low-risk group (Supplementary Figure 1). To further assess the independent prognostic value of the glycometabolic gene signature, we conducted univariate and multivariate Cox regression analyses based on the risk score and clinical characteristics in TARGET dataset. Univariate Cox regression analysis identified that the risk score could serve as a survival related variable (Figure 4A). Subsequently, multivariate Cox regression analysis was used to determine the independent prognostic value of the risk score, and the results showed that the risk score was an independent prognostic factor (Figure 4B). Thus, these results demonstrated that the risk score was significantly associated with metastasis, and could serve as an independent prognostic factor for patients with OS.

Besides, to better provide an applicable quantitative tool for clinic practice, we constructed a new nomogram based on the independent prognostic factor to predict the prognosis of patients with OS in TARGET dataset at different years after diagnosis. Then, total points of each patient could be calculated according to the survival rate at 1, 3-, and 5-years (Figure 5A). The results showed that the overall survival of patients at 1, 3-, and 5-years reduced along with the increase of total scores. The C-index was 0.76 and the calibration curve was similar to the ideal curve (Figure 5B). Furthermore, decision curves were drawn to assess the clinical utility of the risk score model. The results showed that the risk score model could yield more net benefit for predicting the 5-year survival rates than both treat-all and treat-none strategies (Figure 5C). All these findings indicated that the nomogram showed favorable predictive ability and application value.

In order to understand the prognostic values of these four glycometabolic genes, we firstly investigated the expression of each gene between high- and low-risk groups in TARGET dataset. The expression of PFKFB3 and GPX7 was upregulated while the expression of PRKACB and SEPHS2 was downregulated in high-risk group when compared to that in low-risk group (Figure 6A). The expression trends of these prognostic genes in validation cohort were similar to that in TARGET dataset except for SEPHS2 (Figure 6B). Then the prognostic values of glycometabolic genes were evaluated in TARGET dataset. As a result, the elevated expression of PRKACB was associated with favorable clinical outcomes (Figure 6C), whereas the patients with high expression of GPX7 had poorer prognosis (Figure 6E). Nevertheless, there was no distinct difference between the expression level and overall survival in terms of SEPHS2 and PFKFB3 (Figures 6D, F). These findings showed that glycometabolic genes including PRKACB and GPX7 were significantly correlated with the prognosis of patients with OS.

In order to understand the potential biological functions, the DEGs based on risk score were explored in the TARGET dataset. As shown in Supplementary Table 3, a total of 940 DEGs including 336 upregulated genes and 604 downregulated genes were identified between risk groups. Subsequently, we conducted GSEA to explore the discrepancies of biological processes and pathways between two risk groups. The results of GO enrichment analysis indicated that DEGs were enriched in multiple biological processes including activation of innate immune response, adaptive immune response, antigen processing and presentation, and antigen receptor-mediated signaling pathway (Supplementary Table 4). The top 10 enriched biological processes are shown in Figure 7A. Results from KEGG analysis showed that the DEGs were mainly involved in allograft rejection, antigen processing and presentation, autoimmune thyroid disease, chemokine signaling pathway, complement and coagulation cascades, and cytokine-cytokine receptor interaction (Supplementary Table 5). The top 10 enriched pathways are shown in Figure 7B. The above results suggested that these DEGs were involved in various immune related processes and pathways, indicating widespread correlation between glycometabolism and immune status. In addition, we further calculated the scores of each biological process and pathway in each patient using ssGSEA. The heatmaps depicted the top 10 biological processes and pathways with significant differences of scores between the two risk groups, respectively (Figures 7C, D). The results revealed that several immune related processes and pathways were downregulated in high-risk group.

In order to further evaluate the correlation between immune and the glycometabolic gene signature, we next carried out ssGSEA to analyze the immune infiltration in TARGET dataset. As shown in Figure 8A, the heatmap described the abundance of 28 immunocytes between the two risk groups. A total of 26 immunocytes were significantly decreased in the patients with high-risk scores in comparison to the patients with low-risk scores (Figure 8B). Higher level of Type 17 T helper cell was found in low-risk group while there was no significant difference between the two risk groups. Moreover, the immune infiltration scores and tumorpurity were also assessed between two risk groups. As shown in Figures 8C–F, the immune score, stromal score, and ESTIMATE score were higher while the tumorpurity was lower in low-risk group when compared to those in high-risk group. Besides, the correlation between these 4 prognostic genes and the immunocytes was further analyzed. The results showed that the SEPHS2 was positively correlated with 14 immune cells, PRKACB was positively correlated with 17 immune cells, GPX7 was negatively correlated with 11 immune cells, and PFKFB3 was negatively correlated with 2 immune cells, respectively (Figure 8G).

The IC50 values of chemotherapy drugs were compared in order to provide treatment options for patients with OS. The IC50 values of cisplatin, methotrexate, and paclitaxel showed no significance between two risk groups (Figures 9A–C). Significantly, the IC50 value of doxorubicin was lower in high-risk group, suggesting that the patients in high-risk group may be more sensitive to doxorubicin (Figure 9D).

Both GO and KEGG analyses were performed to further identify the potential biological functions of these 4 prognostic genes. Several glycometabolism related biological processes such as fructose metabolic process, and carbohydrate phosphorylation were enriched (Figure 10A). KEGG analyses showed that these 4 genes were related to the pathway of thyroid hormone synthesis, fructose and mannose metabolism, and hedgehog signaling pathway (Figure 10B). In the PPI network, we found that these 4 prognostic genes could directly or indirectly interact with other 50 genes (Figure 10C). Furthermore, a lncRNA-miRNA-mRNA regulatory networks was conducted to identify the lncRNAs and miRNAs regulating the expression of prognostic genes. Finally, a ceRNA network composed by 4 glycometabolic genes, 148 miRNAs, and 91 lncRNAs, were identified (Figure 10D and Supplementary Table 6).

The expression of these 4 prognostic genes were verified using immunohistochemistry (Figure 11). It was shown that GPX7 and PFKFB3 were up-regulated, and SEPHS2 was down-regulated in OS tissues in comparison to those in adjacent normal tissues. However, the expression of PRKACB appeared to be no significant difference between OS tissue and adjacent normal tissue. These findings indicated that GPX7, PFKFB3, and SEPHS2 may play critical roles in the development of OS.