A glycometabolic gene signature associating with immune infiltration and chemosensitivity and predicting the prognosis of patients with osteosarcoma

Background: Accumulating evidence has suggested that glycometabolism plays an important role in the pathogenesis of tumorigenesis. However, few studies have investigated the prognostic values of glycometabolic genes in patients with osteosarcoma (OS). This study aimed to recognize and establish a glycometabolic gene signature to forecast the prognosis, and provide therapeutic options for patients with OS.

Methods: Univariate and multivariate Cox regression, LASSO Cox regression, overall survival analysis, receiver operating characteristic curve, and nomogram were adopted to develop the glycometabolic gene signature, and further evaluate the prognostic values of this signature. Functional analyses including Gene Ontology (GO), kyoto encyclopedia of genes and genomes analyses (KEGG), gene set enrichment analysis, single-sample gene set enrichment analysis (ssGSEA), and competing endogenous RNA (ceRNA) network, were used to explore the molecular mechanisms of OS and the correlation between immune infiltration and gene signature. Moreover, these prognostic genes were further validated by immunohistochemical staining.

Results: A total of four genes including PRKACB, SEPHS2, GPX7, and PFKFB3 were identified for constructing a glycometabolic gene signature which had a favorable performance in predicting the prognosis of patients with OS. Univariate and multivariate Cox regression analyses revealed that the risk score was an independent prognostic factor. Functional analyses indicated that multiple immune associated biological processes and pathways were enriched in the low-risk group, while 26 immunocytes were down-regulated in the high-risk group. The patients in high-risk group showed elevated sensitivity to doxorubicin. Furthermore, these prognostic genes could directly or indirectly interact with other 50 genes. A ceRNA regulatory network based on these prognostic genes was also constructed. The results of immunohistochemical staining showed that SEPHS2, GPX7, and PFKFB3 were differentially expressed between OS tissues and adjacent normal tissues.

Conclusion: The preset study constructed and validated a novel glycometabolic gene signature which could predict the prognosis of patients with OS, identify the degree of immune infiltration in tumor microenvironment, and provide guidance for the selection of chemotherapeutic drugs. These findings may shed new light on the investigation of molecular mechanisms and comprehensive treatments for OS.

1. Introduction

Osteosarcoma (OS) is derived from primitive osteogenic mesenchymal cells, which is recognized as the most common type of malignant bone tumor in childhood and adolescence (1). OS prefers to occur in the metaphysis of long bones, especially around the knees, and is characterized by high rates of metastasis and progression (2). With the improvement of neoadjuvant therapy and surgical resection, more than two-thirds of patients with localized lesions are likely to achieve long-term survival (3, 4). However, approximately 30% of the patients with non-metastasis at diagnosis suffer from lung metastasis after comprehensive treatments (5). What’s worse, the 5-year survival rate of the patients with distant metastasis at diagnosis is still unfavorable, and patients with metastatic OS or chemoresistance are lacking of the effective therapeutic interventions (6). It is believed that early detection and intervention is of great importance to improve the overall survival of patients with OS. Thus, novel prognostic biomarkers are urgently needed for the diagnosis and treatment of OS.

Glycometabolism is a universal pathway involved in cell growth and survival. Normal and non-proliferating cells obtain energy though oxidative phosphorylation under aerobic conditions (7). Whereases cancer cells mainly depend on aerobic glycolysis for adenosine triphosphate (ATP) production for the requirement of rapid proliferation and invasion, and display enhanced glucose uptake for the compensation against the low energy yield of aerobic glycolysis (8, 9). Increasing studies have shown that elevated glucose uptake and aerobic glycolysis in cancer cells is associated with distant metastasis and unfavorable prognosis (10). Recent studies also have found that dysregulation of glycometabolism is involved in tumorigenesis and treatment of OS (11–13). Increased aerobic glycolysis facilitates cell growth, metastasis, and chemoresistance in OS (14, 15). Targeting aerobic glycolysis may be an attractive therapeutic option for the treatment of OS (16, 17). However, the mechanisms of glycometabolic genes in OS remain largely unknown. Recently, several signatures based on glycometabolic genes and lncRNAs have been established in multiple tumors, and these signatures can contribute to elucidate the association between glycometabolism and prognosis (18–21). Nevertheless, the study on glycometabolic gene signature in OS still remains limited.

In this study, we obtained the expression profiles and corresponding clinical information of OS patients from the therapeutically applicable research to generate effective treatments (TARGET) database and the gene expression omnibus (GEO) database and then constructed a novel glycometabolic gene signature to predict the clinical outcomes of patients with OS. Meanwhile, the immune status and chemotherapy drug sensitivity between high- and low-risk groups were also evaluated. Finally, the potential mechanisms of these prognostic genes in tumorigenesis of OS and their up-regulatory network were further investigated. These findings may provide novel prognostic biomarkers and molecular mechanisms for the diagnosis and treatment of OS.

2. Materials and methods

2.1. Selection of datasets and data acquisition

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.

FIGURE 1

Figure 1. The flowchart of this study.

2.2. Construction and validation of a glycometabolic gene signature

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 = ${\displaystyle {\sum }_{i=1}^n\beta i\times 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.

2.3. Independent prognostic analysis and construction of a nomogram

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.

2.4. Functional enrichment analysis

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.

2.5. Analysis of immune infiltration

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.

2.6. Drug sensitivity analysis

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.

2.7. Protein-protein interaction (PPI) network

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).

2.8. Construction of a ceRNA network

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.

2.9. OS samples and immunohistochemical staining

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.

2.10. Statistical analysis

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.

3. Results

3.1. Construction of a glycometabolic gene signature in OS

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).

FIGURE 2

Figure 2. Construction of a glycometabolic gene signature in the training cohort. (A) Univariate Cox regression analysis of glycometabolic genes. (B) LASSO Cox regression analysis of glycometabolic genes. (C) Selection of the optimal penalty parameter for LASSO Cox regression. (D) Coefficients of the four glycometabolic genes though LASSO Cox regression. (E) Heat map visualizing the expression of the four glycometabolic genes between risk groups. (F) The scatter plot of patients’ risk score and survival time. (G) Kaplan–Meier curves of patients in low- and high-risk groups. (H) Time-dependent ROC curves of 1-, 2-, 3-, 4-, and 5-years survival.

3.2. Validation of the glycometabolic gene signature

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.

FIGURE 3

Figure 3. Validation of the glycometabolic gene signature in testing and validation cohorts. (A) Heat map visualizing the expression of the four glycometabolic genes between risk groups in testing cohort. (B) The scatter plot of patients’ risk score and survival time in testing cohort. (C) Kaplan–Meier curves of patients in testing cohort. (D) Time-dependent ROC curves of 1-, 2-, 3-, 4-, and 5-years survival in testing cohort. (E) Heat map visualizing the expression of the four glycometabolic genes between risk groups in validation cohort. (F) The scatter plot of patients’ risk score and survival time in validation cohort. (G) Kaplan–Meier curves of patients in validation cohort. (H) Time-dependent ROC curves of 1-, 2-, 3-, 4-, and 5-years survival in validation cohort.

3.3. Independent prognostic value of the glycometabolic gene signature and construction of a nomogram

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.

FIGURE 4

Figure 4. Independent prognostic value of glycometabolic gene signature. (A) Univariate Cox regression analysis and (B) multivariate Cox regression analysis of the correlations between prognosis and age, gender, and risk score.

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.

FIGURE 5

Figure 5. Construction of a nomogram based on risk score. (A) A nomogram for predicting 1/3/5-year survival rates of OS patients. (B) 1/3/5-year calibration curves of the nomogram. (C) Decision curve of the risk score.

3.4. The expression of prognostic genes and their prognostic values

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.

FIGURE 6

Figure 6. The expression of the prognostic genes and their prognostic values. (A) The expression of the prognostic genes in TARGET dataset. (B) The expression of the prognostic genes in validation cohort. (C–F) The Kaplan–Meier curves for overall survival of osteosarcoma patients in the high and low PRKACB (C), SEPHS2 (D), GPX7 (E), and PFKFB3 (F) expression groups from the TARGET dataset.

3.5. Functional enrichment analyses of DEGs based on risk score

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.

FIGURE 7

Figure 7. Functional enrichment analyses of DEGs between risk groups. (A) The top 10 GO terms enriched in the low-risk group using GSEA analysis. (B) The top 10 KEGG terms enriched in the low-risk group using GSEA analysis. (C) Heatmap of the top 10 GO terms between risk groups using ssGSEA analysis. (D) Heatmap of the top 10 KEGG terms between risk groups using ssGSEA analysis.

3.6. Evaluation of immune microenvironment characteristics between risk groups

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).

FIGURE 8

Figure 8. Evaluation the differences of immune microenvironment between risk groups. (A) Heatmap of immune cell populations between risk groups. (B) Proportion of immune cell subpopulations between risk groups. (C–F) Immune scores (C), stromal score (D), ESTIMATE score (E), and tumorpurity score (F) between risk groups. (G) Heatmap visualizing the correlation between the four prognostic genes and immune cell subpopulations.

3.7. Analysis of chemotherapy drugs sensitivity between risk groups

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).

FIGURE 9

Figure 9. Comparation the differences of chemotherapy drug sensitivity between risk groups. (A) Cisplatin. (B) Methotrexate. (C) Paclitaxel. (D) Doxorubicin. The differences were assessed by Wilcoxon tests.

3.8. Functions of prognostic genes and their regulatory mechanisms

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).

FIGURE 10

Figure 10. Functional analyses of prognostic genes. (A) Biological processes of GO functional analysis. (B) KEGG pathway analysis. (C) Protein-protein interaction network. (D) ceRNA regulatory network.

3.9. Validation the expression of prognostic genes in OS tissues

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.

FIGURE 11

Figure 11. Validation the expression of prognostic genes in OS tissues. The immunohistochemical staining was performed to detect the expression of PRKACB, SEPHS2, GPX7, and PFKFB3 between OS tissues and the corresponding adjacent tissues (magnification 20 and 40×).

4. Discussion

Glycometabolism in tumor cells is characterized by the enhanced glucose uptake and aerobic glycolysis (30). Aerobic glycolysis allows the conversion of glucose into pyruvate eventually contributing to the production of lactate. This energy metabolic reprogramming promotes energy generation and thus facilitates the proliferation, invasion, and chemoresistance of OS cells (31). Several factors, including glycolytic enzymes (e.g., GLUT1), oncogenes (e.g., KRT17), transcription factors (e.g., HIF1α), tumor suppressors (e.g., p53), and related signaling pathways, have been reported to regulate the glycometabolism of OS cells (32–34). Deep insight into these factors may help to formulate therapeutic strategies. In this study, a total of 4 glycometabolism related prognostic genes, including PRKACB, SEPHS2, GPX7, and PFKFB3, were identified using the univariate Cox regression analysis followed by a LASSO regression analysis. Then a glycometabolic gene signature was established based on these genes. Further analysis showed that the glycometabolic gene signature had a good performance and application value in predicting the prognosis of patients with OS. Meanwhile, this prognostic signature was significantly associated with immune microenvironment, and could provide options for the application of chemotherapy drugs. Functional enrichment showed that these 4 genes were enriched in multiple biological processes and pathways which have been demonstrated involving in the pathogenesis of OS. PPI and ceRNA network revealed that these 4 prognostic genes could interact with several genes, and their expression might be regulated by a variety of miRNAs and lncRNAs. These findings indicated that the novel glycometabolic gene signature may serve as a new diagnostic biomarker, and have the potential to predict the prognosis of patients with OS, which also may shed new light on the underlying mechanisms of OS from the perspective of glycometabolism.

Osteosarcoma (OS) is one of the highly aggressive tumors with high rates of metastasis and recurrence (1). The development of biomarkers for early diagnosis and predicting prognosis will facilitate the comprehensive treatment of patients with OS (35). Recently, several gene signatures have been established from different perspectives which could serve as diagnostic biomarkers and predict the prognosis of OS. These prognostic signatures can be divided into several categories, including immune related gene signatures (36, 37), hypoxic gene signatures (38, 39), cell death related gene signature (40–43), noncoding RNA related prognostic signatures (44, 45), epigenetics related gene signatures (46, 47), and metastasis associated gene signatures (48, 49). In addition, some metabolism related prognostic signatures have also been developed (50, 51). Glycometabolism, as one of the most important metabolic pathways in tumor, is of great importance to the initiation and progression of OS (10). Previous studies have showed that glycolysis related gene signatures are able to predict the prognosis of patients undergoing OS (52, 53). It is well known that the acquisition of energy in tumor cells also depends on oxidative phosphorylation (10). Hence, comprehensive analysis the roles of glycometabolism related genes will be more helpful to elucidate the correlation between glycometabolism and OS. In this study, we identified a novel prognostic gene signature based on 4 glycometabolic genes including PRKACB, SEPHS2, GPX7, and PFKFB3. Comprehensive analyses such as KM survival analysis, ROC curve, and nomogram indicated that this gene signature was robust and had a good performance in predicting the prognosis of patients with OS. These findings may enrich the study on the association between glycometabolic genes and the prognosis of OS.

In current study, bioinformatics analysis showed these 4 glycometabolic genes were differentially expressed between high- and low- risk groups. Significantly, the results of IHC staining also indicated the distinct difference of prognostic genes including SEPHS2, GPX7, and PFKFB3 between tumor and adjacent normal tissues. Hence, we speculate these glycometabolic genes may play critical roles in the occurrence, and development of OS. However, due to the small sample size of OS tumor tissues included in this study, we failed to observe the significant difference of PRKACB between tumor and adjacent normal tissues. What is more, increasing studies have investigated their roles in the tumorigenesis. For instance, Lower PRKACB expression was found in tumor tissues and significantly associated with unfavorable overall survival in patients with colorectal carcinoma (54). It was shown that SEPHS2 was elevated in breast tumor samples, and its overexpression was correlated with advanced tumor grade, suggesting that SEPHS2 may serve as a prognostic marker and therapeutic target for patients with breast cancer (55). GPX7 is a member of the glutathione peroxidase (GPx) family with weak GPx activity (56). GPX7 has been confirmed to inhibit tumorigenesis, and may function as a tumor suppressor in multiple tumors (57, 58). PFKFB3, as a key regulator of glycolysis, has been implicated in tumorigenesis, angiogenesis chemoresistance, and tumor microenvironment (59). Furthermore, it was shown that upregulated PFKFB3 could accelerate cell growth and metastasis, which may be a potential biomarker for OS (60, 61). Taking together, it is believed that the prognostic value of this signature may be attributed to the potential impacts of these glycometabolic genes to the tumorigenesis of OS.

The tumor immune microenvironment has long been shown to be strongly correlated with tumor development, recurrence and metastasis (62). Previous studies have suggested that Warburg effect participates in immunomodulation in the tumor microenvironment, and promotes immune evasion by against macrophage immunosurveillance (63, 64). However, little is known about the association between glycometabolism and immune microenvironment in OS. Herein, GSEA was performed and found that DEGs were enriched several immune related biological functions and pathways, such as in activation of innate immune response, adaptive immune response, antigen processing and presentation, natural killer (NK) cell mediated cytotoxicity, T cell receptor signaling pathway, and B cell receptor signaling pathway. ssGSEA analysis was conducted to further determine the correlations between glycometabolism and immune infiltration. The ssGSEA results indicated that a total of 26 immunocytes, including activated B cell, activated CD8 T cell, NK cell, and activated dendritic cell, were lower infiltration in high-risk OS samples than that in low-risk OS samples. Further analysis showed that significant links were found between these 4 glycometabolic genes and multiple immunocytes. Among these cells, NK cells, as a type of lymphocytes of the innate immune system, are able to recognize OS cells, release cytokines, and induce the cell lysis via multiple mechanisms (65). Enhancing NK cell-versus-OS effect has become a promising treatment for OS (66). Dendritic cells (DCs) are known as the central regulators of the adaptive immune response, and are essential for T cell mediated cancer immunity. It has been demonstrated that OS immunotherapy using a DCs-fused tumor vaccine could significantly facilitate T cells proliferation, and promote the tumor-cytotoxic activity of cytotoxic T cells (67). The combination of DCs and anti-glucocorticoid-induced tumor necrosis factor receptor antibodies had the ability to enhance systemic immune responses, and inhibit primary OS growth (68). In addition, we found that the decrease in immune score, stromal score, and ESTIMATE score while the increase in tumorpurity was detected in high-risk group. Stromal cells and immune cells are considered as the main components of tumor microenvironment, which is essential for tumorigenesis, invasion and immune infiltration (69). The increase of both stromal cells and immune cells in tumor environment facilitates immune resistance and immune escape, and the OS patients with higher stromal and immune scores may be appropriate for immune checkpoint inhibitor treatment (70, 71). Hence, these findings indicated that the patients with low- risk scores may obtain more benefit from immune checkpoint inhibitor treatment. Taking together, this glycometabolic gene signature may contribute to estimate the immune microenvironment of OS samples, and thus help formulate individualized anti-tumor immunotherapies.

It is known that 4 agents, including cisplatin, methotrexate, paclitaxel, and doxorubicin, are regarded as the first-line chemotherapy drug regarding OS treatment. Nevertheless, numerous studies have found that the resistance to these drugs can lead to unfavorable outcomes of patients with OS (72). Currently, the prediction of the responses to chemotherapies is still limited by lack of effective biomarkers. In present study, we investigated the drug sensitivity between two risk groups, and the results showed that the patients in high-risk group were more sensitive to doxorubicin. As mentioned above, an increase of tumorpurity score was found in high-risk group, indicating that there were much more OS tumor cells in high-risk samples when compared to that in low-risk samples. Meanwhile, the correlation between risk scores and metastasis suggested that higher proportion of patients with metastasis was found in high-risk group in comparation to that in low-risk group. Hence, we speculate that the patients with high tumorpurity score or metastatic OS may be more effective to the chemotherapy with doxorubicin. These findings indicated that the risk score model may help guide the selection of chemotherapy drugs. Although multiple studies have demonstrated the pharmacological effect on OS treatment, it is necessary to conduct large sample size of cohort study to explore the clinical effectiveness of doxorubicin for the patients with high-risk scores.

What is more, we also tried to identify more details about these 4 prognostic genes in the pathogenesis of OS. GO functional and KEGG analyses showed that these 4 prognostic genes were enriched in the biological processes of fructose metabolic process, carbohydrate phosphorylation, and the pathways of thyroid hormone synthesis, fructose and mannose metabolism, and hedgehog signaling pathway. These biological processes and pathways may participate in the regulation of glycometabolism in OS. The PPI network suggested that these 4 prognostic genes were able to directly or indirectly interact with other 50 genes. Some of these genes have been shown involving the initiation and progression of OS. For instance, Wang and Sun (73) demonstrated that inhibition of FOXO1 was able to suppress both proliferation and metastases in OS cells, indicating that FOXO1 may be a potential target for the treatment of OS. In addition, CREB1 could lead to OS progression and metastasis through promoting epithelial-mesenchymal transition (74). Hence, deeper insight into these genes in the network may help uncover the underlying mechanisms of these 4 prognostic genes in the development of OS. CeRNA network is regarded to be a prevalent form of post-transcriptional regulation of gene expression in mammals (75). Message RNA (mRNA), long non-coding RNA (lncRNA), pseudogene, and circular RNA (circRNA) can affect the stability or translation activity of target RNAs by competitive combination with microRNA (miRNA). Increasing studies have confirmed that ceRNA regulatory network is widely involved in the occurrence and development of OS (76). To further understand the upstream mechanisms of prognostic genes, a ceRNA network was performed to identify the target miRNAs and the corresponding lncRNAs. The results showed that a total of 148 miRNAs and 91 lncRNAs were linked with these four prognostic genes. These findings may facilitate a better understanding of these 4 prognostic genes in the pathogenesis of OS.

Some limitations of this study should be interpreted. First, both training cohort and testing cohorts used in current study contain relatively small sample size. Meanwhile, some clinical features including tumor size and grade of OS patients are missing in TARGET dataset. Hence, we did not perform stratification analysis to investigate the correlations between risk score and tumor size and grade. Second, due to limited clinical tumor samples used in current study, we failed to find distinct difference of PRKACB between tumor and adjacent normal tissues. Given that low expression of PRKACB has been demonstrated leading to unfavorable clinical outcomes in patients undergoing breast cancer and colorectal carcinoma (54, 77), we speculate that the low expression of PRKACB may contribute to the occurrence and development of OS, thereby bringing about inferior clinical outcomes. However, no studies have reported the functions of these prognostic genes including PRKACB, SEPHS2, and GPX7 in the pathogenesis of OS. Therefore, more OS tumor samples and multiple OS cell lines should be applied to determine the expression of these prognostic genes and their roles in the tumorigenesis of OS.

In summary, our research identified a novel signature based on four glycometabolic genes, and constructed a risk score model that can predict the survival, immune infiltration, and chemosensitivity of patients with OS. These findings may provide new insights into the role of glycometabolic genes in the molecular mechanisms of OS, and help develop novel diagnostic and therapeutic strategies.

Data availability statement

The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors (78).

Ethics statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Affiliated Hospital of Guizhou Medical University (No. GZYD003-201753035). The patients/participants provided their written informed consent to participate in this study.

Author contributions

HS and KL conceived the study. FW and KY performed the bioinformatic and statistical analysis. RP participated in collecting the data from online database. ZX and PL participated in figures and charts drawing. YX supervised the implement of the current study. HS, FW, and KL provided fund and prepared the original draft. HS reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the Science and Technology Fund of Guizhou Provincial Health Commission (gzwjkj2020-1-120), the Youth Fund Cultivation Program of National Natural Science Foundation of Affiliated Hospital of Guizhou Medical University (gyfynsfc-2021-12), and Graduate Scientific Research Fund Project of Guizhou [YJSKYJJ (2021)157].

Conflict of interest

The 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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2023.1115759/full#supplementary-material

Supplementary Figure 1 | Non-metastasis and metastasis cases between low- and high-risk groups in TARGET dataset.

Footnotes

1. ^ https://ocg.cancer.gov/programs/target
2. ^ https://www.ncbi.nlm.nih.gov/
3. ^ http://www.broad.mit.edu/gsea/msigdb

References

1. Czarnecka AM, Synoradzki K, Firlej W, Bartnik E, Sobczuk P, Fiedorowicz M, et al. Molecular biology of osteosarcoma. Cancers. (2020) 12:2130.

Google Scholar

2. Misaghi A, Goldin A, Awad M, Kulidjian AA. Osteosarcoma: a comprehensive review. SICOT J. (2018) 4:12.

Google Scholar

3. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: current treatment and a collaborative pathway to success. J Clin Oncol. (2015) 33:3029–35. doi: 10.1200/JCO.2014.59.4895

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Whelan JS, Davis LE. Osteosarcoma, chondrosarcoma, and chordoma. J Clin Oncol. (2018) 36:188–93.

Google Scholar

5. Harrison DJ, Geller DS, Gill JD, Lewis VO, Gorlick R. Current and future therapeutic approaches for osteosarcoma. Expert Rev Anticancer Ther. (2018) 18:39–50.

Google Scholar

6. He H, Ni J, Huang J. Molecular mechanisms of chemoresistance in osteosarcoma (Review). Oncol Lett. (2014) 7:1352–62.

Google Scholar

7. Wilson DF. Oxidative phosphorylation: regulation and role in cellular and tissue metabolism. J Physiol. (2017) 595:7023–38.

Google Scholar

8. Ancey PB, Contat C, Meylan E. Glucose transporters in cancer - from tumor cells to the tumor microenvironment. FEBS J. (2018) 285:2926–43. doi: 10.1111/febs.14577

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Lunt SY, Vander Heiden MG. Aerobic glycolysis: meeting the metabolic requirements of cell proliferation. Annu Rev Cell Dev Biol. (2011) 27:441–64.

Google Scholar

10. Ghanavat M, Shahrouzian M, Deris Zayeri Z, Banihashemi S, Kazemi SM, Saki N. Digging deeper through glucose metabolism and its regulators in cancer and metastasis. Life Sci. (2021) 264:118603. doi: 10.1016/j.lfs.2020.118603

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Deng B, Deng J, Yi X, Zou Y, Li C. ROCK2 promotes osteosarcoma growth and glycolysis by up-regulating HKII via Phospho-PI3K/AKT signalling. Cancer Manag Res. (2021) 13:449–62. doi: 10.2147/CMAR.S279496

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shen Y, Zhao S, Wang S, Pan X, Zhang Y, Xu J, et al. S1P/S1PR3 axis promotes aerobic glycolysis by YAP/c-MYC/PGAM1 axis in osteosarcoma. EBioMedicine. (2019) 40:210–23. doi: 10.1016/j.ebiom.2018.12.038

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wan J, Liu Y, Long F, Tian J, Zhang C. circPVT1 promotes osteosarcoma glycolysis and metastasis by sponging miR-423-5p to activate Wnt5a/Ror2 signaling. Cancer Sci. (2021) 112:1707–22. doi: 10.1111/cas.14787

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Shen Y, Xu J, Pan X, Zhang Y, Weng Y, Zhou D, et al. LncRNA KCNQ1OT1 sponges miR-34c-5p to promote osteosarcoma growth via ALDOA enhanced aerobic glycolysis. Cell Death Dis. (2020) 11:278. doi: 10.1038/s41419-020-2485-1

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang Q, Wu J, Zhang X, Cao L, Wu Y, Miao X. Transcription factor ELK1 accelerates aerobic glycolysis to enhance osteosarcoma chemoresistance through miR-134/PTBP1 signaling cascade. Aging. (2021) 13:6804–19. doi: 10.18632/aging.202538

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hu R, Chen S, Yan J. Blocking circ-CNST suppresses malignant behaviors of osteosarcoma cells and inhibits glycolysis through circ-CNST-miR-578-LDHA/PDK1 ceRNA networks. J Orthop Surg Res. (2021) 16:300. doi: 10.1186/s13018-021-02427-0

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wang Q, Liu MJ, Bu J, Deng JL, Jiang BY, Jiang LD, et al. miR-485-3p regulated by MALAT1 inhibits osteosarcoma glycolysis and metastasis by directly suppressing c-MET and AKT3/mTOR signalling. Life Sci. (2021) 268:118925. doi: 10.1016/j.lfs.2020.118925

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bai Y, Lin H, Chen J, Wu Y, Yu S. Identification of prognostic glycolysis-related lncRNA signature in tumor immune microenvironment of hepatocellular carcinoma. Front Mol Biosci. (2021) 8:645084. doi: 10.3389/fmolb.2021.645084

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bi J, Bi F, Pan X, Yang Q. Establishment of a novel glycolysis-related prognostic gene signature for ovarian cancer and its relationships with immune infiltration of the tumor microenvironment. J Transl Med. (2021) 19:382. doi: 10.1186/s12967-021-03057-0

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li W, Xu M, Li Y, Huang Z, Zhou J, Zhao Q, et al. Comprehensive analysis of the association between tumor glycolysis and immune/inflammation function in breast cancer. J Transl Med. (2020) 18:92. doi: 10.1186/s12967-020-02267-2

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yu S, Hu C, Cai L, Du X, Lin F, Yu Q, et al. Seven-gene signature based on glycolysis is closely related to the prognosis and tumor immune infiltration of patients with gastric cancer. Front Oncol. (2020) 10:1778. doi: 10.3389/fonc.2020.01778

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chen X, Ji ZL, Chen YZ. TTD: therapeutic target database. Nucleic Acids Res. (2002) 30:412–5.

Google Scholar

23. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. (2011) 27:1739–40. doi: 10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. (2005) 102:15545–50.

Google Scholar

25. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lin R, Fogarty CE, Ma B, Li H, Ni G, Liu X, et al. Identification of ferroptosis genes in immune infiltration and prognosis in thyroid papillary carcinoma using network analysis. BMC Genomics. (2021) 22:576. doi: 10.1186/s12864-021-07895-6

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. (2015) 43:D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. (2011) 12:77. doi: 10.1186/1471-2105-12-77

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Spencer NY, Stanton RC. The Warburg effect, lactate, and nearly a century of trying to cure cancer. Semin Nephrol. (2019) 39:380–93. doi: 10.1016/j.semnephrol.2019.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Li G, Li Y, Wang DY. Overexpression of miR-329-3p sensitizes osteosarcoma cells to cisplatin through suppression of glucose metabolism by targeting LDHA. Cell Biol Int. (2021) 45:766–74. doi: 10.1002/cbin.11476

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Gomes AS, Ramos H, Soares J, Saraiva L. p53 and glucose metabolism: an orchestra to be directed in cancer therapy. Pharmacol Res. (2018) 131:75–86. doi: 10.1016/j.phrs.2018.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yan X, Yang C, Hu W, Chen T, Wang Q, Pan F, et al. Knockdown of KRT17 decreases osteosarcoma cell proliferation and the Warburg effect via the AKT/mTOR/HIF1alpha pathway. Oncol Rep. (2020) 44:103–14. doi: 10.3892/or.2020.7611

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Yi W, Tu MJ, Liu Z, Zhang C, Batra N, Yu AX, et al. Bioengineered miR-328-3p modulates GLUT1-mediated glucose uptake and metabolism to exert synergistic antiproliferative effects with chemotherapeutics. Acta Pharm Sin B. (2020) 10:159–70. doi: 10.1016/j.apsb.2019.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Smrke A, Anderson PM, Gulia A, Gennatas S, Huang PH, Jones RL. Future directions in the treatment of osteosarcoma. Cells. (2021) 10:172.

Google Scholar

36. Fan L, Ru J, Liu T, Ma C. Identification of a novel prognostic gene signature from the immune cell infiltration landscape of osteosarcoma. Front Cell Dev Biol. (2021) 9:718624. doi: 10.3389/fcell.2021.718624

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Liu W, Xie X, Qi Y, Wu J. Exploration of immune-related gene expression in osteosarcoma and association with outcomes. JAMA Netw Open. (2021) 4:e2119132.

Google Scholar

38. Fu Y, Bao Q, Liu Z, He G, Wen J, Liu Q, et al. Development and validation of a hypoxia-associated prognostic signature related to osteosarcoma metastasis and immune infiltration. Front Cell Dev Biol. (2021) 9:633607. doi: 10.3389/fcell.2021.633607

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wu F, Xu J, Jin M, Jiang X, Li J, Li X, et al. Development and verification of a hypoxic gene signature for predicting prognosis, immune microenvironment, and chemosensitivity for osteosarcoma. Front Mol Biosci. (2022) 8:705148. doi: 10.3389/fmolb.2021.705148

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Hong J, Li Q, Wang X, Li J, Ding W, Hu H, et al. Development and validation of apoptosis-related signature and molecular subtype to improve prognosis prediction in osteosarcoma patients. J Clin Lab Anal. (2022) 36:e24501. doi: 10.1002/jcla.24501

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lei T, Qian H, Lei P, Hu Y. Ferroptosis-related gene signature associates with immunity and predicts prognosis accurately in patients with osteosarcoma. Cancer Sci. (2021) 112:4785–98.

Google Scholar

42. Qi W, Yan Q, Lv M, Song D, Wang X, Tian K. Prognostic signature of osteosarcoma based on 14 autophagy-related genes. Pathol Oncol Res. (2021) 27:1609782. doi: 10.3389/pore.2021.1609782

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhang Y, He R, Lei X, Mao L, Jiang P, Ni C, et al. A novel pyroptosis-related signature for predicting prognosis and indicating immune microenvironment features in osteosarcoma. Front Genet. (2021) 12:780780. doi: 10.3389/fgene.2021.780780

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Luo X, Tang J, Xuan H, Liu J, Li X. Identification and validation of a potent multi-miRNA signature for prediction of prognosis of osteosarcoma patients. Med Sci Monit. (2020) 26:e919272. doi: 10.12659/MSM.919272

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ying T, Dong JL, Yuan C, Li P, Guo Q. The lncRNAs RP1-261G23.7, RP11-69E11.4 and SATB2-AS1 are a novel clinical signature for predicting recurrent osteosarcoma. Biosci Rep. (2020) 40:BSR20191251. doi: 10.1042/BSR20191251

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Tian W, Li Y, Zhang J, Li J, Gao J. Combined analysis of DNA methylation and gene expression profiles of osteosarcoma identified several prognosis signatures. Gene. (2018) 650:7–14. doi: 10.1016/j.gene.2018.01.093

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Yang K, Wang F, Li K, Peng G, Yang H, Xu H, et al. N6-methyladenosine modification-related long non-coding RNAs are potential biomarkers for predicting the prognosis of patients with osteosarcoma. Technol Cancer Res Treat. (2022) 21:15330338221085354. doi: 10.1177/15330338221085354

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Cao MD, Song YC, Yang ZM, Wang DW, Lin YM, Lu HD. Identification of osteosarcoma metastasis-associated gene biomarkers and potentially targeted drugs based on bioinformatic and experimental analysis. Onco Targets Ther. (2020) 13:8095–107. doi: 10.2147/OTT.S256617

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Shi Y, He R, Zhuang Z, Ren J, Wang Z, Liu Y, et al. A risk signature-based on metastasis-associated genes to predict survival of patients with osteosarcoma. J Cell Biochem. (2020) 121:3479–90. doi: 10.1002/jcb.29622

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Chen X, Ye Z, Lou P, Liu W, Liu Y. Comprehensive analysis of metabolism-related lncRNAs related to the progression and prognosis in osteosarcoma from TCGA. J Orthop Surg Res. (2021) 16:523. doi: 10.1186/s13018-021-02647-4

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Li LQ, Zhang LH, Yuan YB, Lu XC, Zhang Y, Liu YK, et al. Signature based on metabolic-related gene pairs can predict overall survival of osteosarcoma patients. Cancer Med. (2021) 10:4493–509. doi: 10.1002/cam4.3984

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Tian K, Qi W, Yan Q, Lv M, Song D. Signature constructed by glycolysis-immune-related genes can predict the prognosis of osteosarcoma patients. Invest New Drugs. (2022) 40:818–30. doi: 10.1007/s10637-022-01228-4

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Yang M, Ma X, Wang Z, Zhang T, Hua Y, Cai Z. Identification of a novel glycolysis-related gene signature for predicting the prognosis of osteosarcoma patients. Aging. (2021) 13:12896–918. doi: 10.18632/aging.202958

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Yao X, Hu W, Zhang J, Huang C, Zhao H, Yao X. Application of cAMP-dependent catalytic subunit beta (PRKACB) low expression in predicting worse overall survival: a potential therapeutic target for colorectal carcinoma. J Cancer. (2020) 11:4841–50. doi: 10.7150/jca.46156

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Nunziata C, Polo A, Sorice A, Capone F, Accardo M, Guerriero E, et al. Structural analysis of human SEPHS2 protein, a selenocysteine machinery component, over-expressed in triple negative breast cancer. Sci Rep. (2019) 9:16131. doi: 10.1038/s41598-019-52718-0

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Chen YI, Wei PC, Hsu JL, Su FY, Lee WH. NPGPx (GPx7): a novel oxidative stress sensor/transmitter with multiple roles in redox homeostasis. Am J Transl Res. (2016) 8:1626–40.

PubMed Abstract | Google Scholar

57. Chen Z, Hu T, Zhu S, Mukaisho K, El-Rifai W, Peng DF. Glutathione peroxidase 7 suppresses cancer cell growth and is hypermethylated in gastric cancer. Oncotarget. (2017) 8:54345–56. doi: 10.18632/oncotarget.17527

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Peng D, Hu T, Soutto M, Belkhiri A, Zaika A, El-Rifai W. Glutathione peroxidase 7 has potential tumour suppressor functions that are silenced by location-specific methylation in oesophageal adenocarcinoma. Gut. (2014) 63:540–51. doi: 10.1136/gutjnl-2013-304612

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Shi L, Pan H, Liu Z, Xie J, Han W. Roles of PFKFB3 in cancer. Signal Transduct Target Ther. (2017) 2:17044.

Google Scholar

60. Deng X, Deng J, Yi X, Zou Y, Liu H, Li C, et al. Ubiquitin-like protein FAT10 promotes osteosarcoma glycolysis and growth by upregulating PFKFB3 via stabilization of EGFR. Am J Cancer Res. (2020) 10:2066–82.

PubMed Abstract | Google Scholar

61. Deng X, Yi X, Deng J, Zou Y, Wang S, Shan W, et al. ROCK2 promotes osteosarcoma growth and metastasis by modifying PFKFB3 ubiquitination and degradation. Exp Cell Res. (2019) 385:111689. doi: 10.1016/j.yexcr.2019.111689

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. (2013) 14:1014–22.

Google Scholar

63. Chen J, Cao X, Li B, Zhao Z, Chen S, Lai SWT, et al. Warburg effect is a cancer immune evasion mechanism against macrophage immunosurveillance. Front Immunol. (2020) 11:621757. doi: 10.3389/fimmu.2020.621757

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Urata K, Kajihara I, Miyauchi H, Mijiddorj T, Otsuka-Maeda S, Sakamoto R, et al. The Warburg effect and tumour immune microenvironment in extramammary Paget’s disease: overexpression of lactate dehydrogenase A correlates with immune resistance. J Eur Acad Dermatol Venereol. (2020) 34:1715–21. doi: 10.1111/jdv.16145

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Tullius BP, Setty BA, Lee DA. Natural killer cell immunotherapy for osteosarcoma. Adv Exp Med Biol. (2020) 1257:141–54.

Google Scholar

66. Omer N, Nicholls W, Ruegg B, Souza-Fonseca-Guimaraes F, Rossi GR. Enhancing natural killer cell targeting of pediatric sarcoma. Front Immunol. (2021) 12:791206. doi: 10.3389/fimmu.2021.791206

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Fang X, Jiang C, Xia Q. Effectiveness evaluation of dendritic cell immunotherapy for osteosarcoma on survival rate and in vitro immune response. Genet Mol Res. (2015) 14:11763–70. doi: 10.4238/2015.October.2.10

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kawano M, Tanaka K, Itonaga I, Iwasaki T, Miyazaki M, Ikeda S, et al. Dendritic cells combined with anti-GITR antibody produce antitumor effects in osteosarcoma. Oncol Rep. (2015) 34:1995–2001. doi: 10.3892/or.2015.4161

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Bejarano L, Jordao MJC, Joyce JA. Therapeutic targeting of the tumor microenvironment. Cancer Discov. (2021) 11:933–59.

Google Scholar

70. Gill J, Gorlick R. Advancing therapy for osteosarcoma. Nat Rev Clin Oncol. (2021) 18:609–24.

Google Scholar

71. Klemm F, Joyce JA. Microenvironmental regulation of therapeutic response in cancer. Trends Cell Biol. (2015) 25:198–213.

Google Scholar

72. Hattinger CM, Patrizio MP, Fantoni L, Casotti C, Riganti C, Serra M. Drug resistance in osteosarcoma: emerging biomarkers, therapeutic targets and treatment strategies. Cancers. (2021) 13:2878.

Google Scholar

73. Wang J, Sun G. FOXO1-MALAT1-miR-26a-5p feedback loop mediates proliferation and migration in osteosarcoma cells. Oncol Res. (2017) 25:1517–27. doi: 10.3727/096504017X14859934460780

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Ma H, Su R, Feng H, Guo Y, Su G. Long noncoding RNA UCA1 promotes osteosarcoma metastasis through CREB1-mediated epithelial-mesenchymal transition and activating PI3K/AKT/mTOR pathway. J Bone Oncol. (2019) 16:100228. doi: 10.1016/j.jbo.2019.100228

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Grull MP, Masse E. Mimicry, deception and competition: the life of competing endogenous RNAs. Wiley Interdiscip Rev RNA. (2019) 10:e1525. doi: 10.1002/wrna.1525

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Wang JY, Yang Y, Ma Y, Wang F, Xue A, Zhu J, et al. Potential regulatory role of lncRNA-miRNA-mRNA axis in osteosarcoma. Biomed Pharmacother. (2020) 121:109627. doi: 10.1016/j.biopha.2019.109627

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Tang J, Luo Y, Wu G. A glycolysis-related gene expression signature in predicting recurrence of breast cancer. Aging. (2020) 12:24983–94. doi: 10.18632/aging.103806

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Wang F, Yang K, Li K, Peng G, Yang H, Xiang Y, et al. The signature of glycometabolism-related genes in predicting the prognosis of patients with osteosarcoma. Res Sq. [Preprint]. (2022). doi: 10.21203/rs.3.rs-1270483/v1

CrossRef Full Text | Google Scholar