Skip to main content

Construction and evaluation of a prognostic model based on the expression of the metabolism-related signatures in patients with osteosarcoma

Abstract

Background

The aim of this study was to screen three major substance metabolism-related genes and establish a prognostic model for osteosarcoma.

Methods

RNA-seq expression data for osteosarcoma were downloaded from The Cancer Genome Atlas (TCGA) and GEO databases. Differentially expressed (DE) RNAs were selected, followed by the selection of metabolic-related DE mRNAs. Using Cox regression analysis, prognostic DE RNAs were identified to construct a prognostic model. Subsequently, independent prognostic clinical factors were screened, and the functions of the long non-coding RNAs (lncRNAs) were analyzed. Finally, the expression of signature genes was further tested in osteosarcoma cells using quantitative reverse transcription quantitative real-time polymerase chain reaction (qRT-PCR) and western blotting.

Results

A total of 432 DE RNAs, comprising 79 DE lncRNAs and 353 DE mRNAs were obtained, and then 107 metabolic-related DE mRNAs. Afterwards signature genes (LINC00545, LINC01537, FOXC2-AS1, CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB) served as optimal combinations, and a prognostic score model was successfully proposed. Three verification datasets (GSE16091, GSE21257, and GSE39055) showed that the model had high specificity and sensitivity. In addition, two independent prognostic clinical factors (age and tumor metastasis) were identified. Finally, the concordance rate between the in silico analysis, qRT-PCR, and western blotting analysis was 88.89% (8/9), suggesting the robustness of our analysis.

Conclusions

The prognostic model based on the nine signature genes accurately predicted the prognosis of patients with osteosarcoma; CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB may serve as metabolism-related biomarkers in osteosarcoma.

Peer Review reports

Background

Osteosarcoma is the most common cancer of bones and joints, accounting for approximately 20% of all benign and malignant bone neoplasia [1]. Approximately 800 new cases of osteosarcoma are diagnosed each year, half of which involve adolescents and young adults [2]. The 10-year overall survival rates are 60%–70% for all patients with osteosarcoma and only 20% for patients with metastatic disease. In particular, the survival of metastatic patients has not improved significantly in recent years, mainly because of limited therapeutic options for metastatic disease [3]. Therefore, there is an urgent need to develop predictive methods to improve the survival of patients [4]. Prognostic biomarkers can predict the possible outcomes of cancers associated with disease progression [5], which could provide considerable assistance in clinical practice for patient stratification, treatment management, and disease status monitoring.

Many factors, such as miRNA-191 [6], long non-coding RNA (lncRNA) HOTTIP [7], and survival [8], have demonstrated prognostic significance in osteosarcoma; however, their detection is inconvenient and costly. Therefore, identification of easily assessed prognostic factors is required. Metabolites are an objective measurement of the phenotype of a cell, as well as of the interaction between the genotype and environment of a cell [9]. Metabolomic signatures have been reported in the liver, colon, breast, lung, ovarian, and pancreatic cancers [9, 10]. Importantly, metabolic alterations have been used to classify osteosarcoma [11, 12]. Amino acids, derivatives, carbohydrates, and lipids are the three metabolites and their metabolism is associated with osteosarcoma progression. Dysregulation of serine and glycine metabolism has been identified as a metabolic regulator that supports tumor cell growth. A previous study showed the anticancer effects of high concentrations of glycine and D-serine in osteosarcoma cells; however, when used at low physiological concentrations, amino acids can induce the proliferation and migration of osteosarcoma cells [13]. Cancer cells can undergo glycolysis, metabolizing glucose into lactic acid, rather than oxidative phosphorylation. The metabolic transition from oxidative phosphorylation to glycolysis is usually considered a hallmark of osteosarcoma [14]. Huang et al. [15] built a risk model based on seven glycolytic genes to effectively evaluate the prognosis of osteosarcoma. Reprogramming of lipid metabolism has been viewed as a new hallmark of tumor malignancy, and disorders of lipid metabolism play key roles in tumorigenesis, tumor progression, and treatment [16]. A recent study showed that the expression of lipid metabolism-related genes is closely related to the immune microenvironment of patients with osteosarcoma and can be used to accurately predict their prognosis [17]. In addition, amino acids and their derivatives, lipids, and carbohydrates regulate gene expression through molecules that sense these macronutrients and act as transcription factors [18]. Alterations in metabolism-related genes offer not only new biomarkers for diagnosis and prognosis but also potential new targets for cancer therapy [19]. However, the expression patterns of the metabolism-related genes in osteosarcoma, and their prognostic features remain unclear.

In this study, we downloaded the RNA-seq expression data of osteosarcoma from The Cancer Genome Atlas (TCGA) database, which served as the training set. Additionally, GSE21257, GSE39055, and GSE16091 were used as validation datasets. Based on these datasets, we selected three major metabolism-related genes and established a prognostic model for osteosarcoma.

Methods

Data sources and preprocessing

RNA-seq expression data of osteosarcoma were downloaded from the TCGA database, and RNA sequence testing was done using the Illumina HiSeq 2000 platform. This dataset contained 265 samples that corresponded to the clinical information of the osteosarcoma samples, and 176 samples with metastasis, clinical survival, and prognosis information were included in this analysis. These samples and their corresponding RNA-Seq expression data were used as training datasets.

Additionally, with “Osteosarcoma and Homo sapiens” as keywords, all public profiles in the NCBI GEO database were searched. Datasets that met the following criteria were included in our study: 1) datasets that were gene expression profile data; 2) datasets that were about the expression profile of solid osteosarcoma tissues; 3) the samples in the dataset had survival and prognosis information; 4) the total sample size was no less than 30; and 5) the samples contained clear classification information for osteosarcoma with and without metastasis. Each dataset was required to at least meet requirements 1 to 4. Finally, three datasets were analyzed as follows.

  1. A.

    GSE21257 [20]: contained 53 osteosarcoma samples, including 34 metastatic and 19 non-metastatic osteosarcoma samples, and all samples had clinical prognostic information. The detection platform was an Illumina human-6 v2.0 expression beadchip. This dataset simultaneously met requirements 1 to 5.

  2. B.

    GSE39055 [21]: contained 37 osteosarcoma samples with clinical survival and prognosis information and was detected on the Illumina HumanhT-12 WG-DASL V4.0 R2 expression beadchip. This dataset simultaneously met requirements 1 to 4.

  3. C.

    GSE16091 [22]: a total of 34 osteosarcoma samples with clinical survival and prognostic information were included. The detection platform used was the Affymetrix Human Genome U133A Array. This dataset simultaneously met criteria 1 to 4.

Three datasets were used as validation datasets for the analysis.

Expression profile reannotation and differentially expressed RNAs screening

According to the RefSeq ID information, the lncRNAs and mRNA in the downloaded datasets were annotated based on 4641 lncRNAs and 19,198 protein-coding genes in the HUGO Gene Nomenclature Committee database [23].

Osteosarcoma samples in the TCGA dataset were divided into metastatic and non-metastatic groups. The limma package (version 3.34.7) [24] in R3.6.1 was used to screen for significantly differentially expressed (DE) RNAs (mRNAs and lncRNAs) between metastatic and non-metastatic groups with thresholds of false discovery rate (FDR) < 0.05, and |log2 fold change (FC)|> 0.5. According to the expression values of DE RNAs in the training dataset samples, bidirectional hierarchical clustering based on the centered Pearson correlation algorithm [25] was performed using pheatmap (version 1.0.8) [26] in R3.6.1.

Screening of DE mRNAs associated with metabolism of three substances

From the Gene Set Enrichment Analysis (GSEA) database [27], the genes associated with metabolism of amino acids and derivatives, carbohydrates, and lipids, were downloaded. These genes were compared with DE mRNAs, and the intersections served as metabolism-related DE mRNAs.

Construction of a prognostic model

Screening for differentially expressed prognostic RNAs

In osteosarcoma samples from the TCGA training set, the DE mRNAs and lncRNAs associated with survival prognosis were selected based on their relationship to metabolic processes, using univariate Cox regression analysis in R3.6.1 survival pack version 2.41–1 [28]. Further multivariate Cox regression analysis was performed on these DE RNAs to screen for those associated with overall survival prognosis and independent prognosis. The log-rank test (p < 0.05) was used as the threshold.

Screening for optimal combinations of DE RNAs

Based on the identified DE RNAs associated with overall survival prognosis and independent prognosis, the optimal combination of DE RNAs was screened using the least absolute shrinkage and selection operator (LASSO) Cox regression model [29] in the R3.6.1 penalized package version 0.9.50 [30]. According to the prognostic factors of each element in the optimized DE RNAs combination and the expression levels of DE RNAs in the training samples, the prognostic score (PS) model was constructed as follows:

$$\mathrm{Prognostic}\;\mathrm{score}\;(\mathrm{PS})={\textstyle\sum_{}}{\mathrm\beta}_{\mathrm{DE}\;\mathrm{RNAs}'\;}\times{\mathrm{Exp}}_{\mathrm{DE}\;\mathrm{RNAs}'}$$

Where βDERNAs represents the prognostic factor of signature DE RNAs, and Exp DERNAs represents the expression level of target DE RNAs in the training data set.

Efficacy evaluation and comparison of prognostic risk prediction models

The PS values of each sample in the TCGA dataset were calculated. Then, according to the PS median value, the samples in the training set were divided into high- (PS ≥ PS median) and low- (PS < PS median) risk groups, and the correlation between the high- and low- groups and actual survival prognosis information was evaluated using the Kaplan–Meier (KM) curve method using the survival pack version 2.41–1 [28]. Meanwhile, in the validation set, the expression values of the target DE mRNAs were extracted, and the PS value of each sample was calculated. The validation samples were also divided into high- and low-risk groups to evaluate the correlation between high- and low-risk groups and actual survival prognosis information.

Independent prognostic clinical factors screening and stratification analysis

In the osteosarcoma samples from the TCGA training set, univariate and multivariate Cox regression analyses in the R3.6.1 survival package (version 2.41–1) [28] were used to screen the independent prognostic factors, and the log-rank p value < 0.05 was selected as the threshold value of significant correlation.

To further study the relationship between independent prognostic clinical factors and risk groups, we performed a stratification analysis for these independent prognostic clinical factors. Briefly, the samples were divided into different groups according to clinical factors, and a correlation analysis of the risk prognosis model was performed in these groups.

Nomogram of the 3- and 5-year survival models

To further investigate the correlation between independent prognostic clinical factors, the PS model, and survival prognosis, we constructed a nomogram of the 3- and 5-year survival rate prediction models based on the obtained independent prognostic factors combined with the risk information identified by the prognostic prediction model using the RMS package version 5.1–2 [31, 32] of R3.6.1.

Comparative analysis of multiple models

Based on the independent prognostic clinical factors obtained above, a prognostic model of clinical factors was constructed and compared with the PS model using the following two parameters: the area under the receiver operating characteristic (ROC) curve and the C-index. The ROC curve is a classification model and the area under the ROC curve is a quantitative index of the ROC curve. The value distribution of the area under the ROC curve is between 0.5 and 1, such that the closer to 1, the better is the performance of the classifier. This evaluation index was approved by R3.6.1 pROC version 1.14.0 [33]. The C-index is used to evaluate the predictive power of the model and estimate the probability that the predicted results are consistent with the actual observed results [34]. A C-index above 0.70 indicates a good model, where as a score of approximately 0.50 indicates a random background. The evaluation index was calculated using R3.6.1 survcomp version 1.34.0 [35].

Functional analysis of signature lncRNAs

With respect to the lncRNAs used for model construction, combined with the metabolism-related DE mRNAs, the Pearson correlation coefficients between the expression (based on the expression level in the TCGA training set) of signature lncRNAs and target DE mRNAs were calculated using the cor function in R3.6.1. A co-expression network between significant DE signature lncRNAs and metabolism-related genes was constructed based on expression correlation and visualized using Cytoscape 3.6.1 [36]. The gene nodes in the network were subjected to biological process of gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using DAVID version 6.8 [37, 38], and a p value < 0.05 was selected as the threshold of enrichment significance.

Cell culture

Human osteosarcoma cells (U2OS and Saos2) and osteoblasts hFOB1.19 were purchased from iCell Bioscience Inc. (Shanghai, China), and cultured in McCoy’s 5A medium containing 10% fetal bovine serum and 1% penicillin/streptomycin. Cells were cultured and passaged at 37 °C in a 5% CO2 humidified atmosphere.

Reverse transcription quantitative real-time PCR (qRT-PCR)

Total RNA was extracted from the cells using RNAiso Plus (Takara, Tokyo, Japan) and reverse transcribed into cDNA using PrimeScript RT Master Mix (Takara) following the manufacturer’s instructions. Power SYBR Green PCR Master Mix (Thermo Fisher, Waltham, USA) was added for gene amplification using a fluorescence ratio PCR instrument (LongGene, Hangzhou, China), with GAPDH as an internal reference. The qRT-PCR reaction was initiated at 50 °C for 3 min, followed by 95 °C for 30 s, 40 cycles at 95 °C for 5 s, and 60 °C for 30 s. The sequences of all primers are listed in Table 1. The relative expression levels of lncRNAs FOXC2-AS1, LINC00545, and LINC01537, and mRNAs CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB were calculated using the 2Ct method.

Table 1 The sequences of all primers used in this study

Western blotting

Total protein was isolated from the cells using RIPA lysis buffer and quantified using a BCA protein assay kit (Beyotime, Shanghai, China) according to the manufacturer’s protocols. Protein samples (20 μg) were then separated by SDS-PAGE and transferred onto polyvinylidene fluoride membranes. After blocked with 5% skim milk at 37 ℃ for 2 h, the membranes were respectively incubated with the anti-PFKFB4 antibody (1:1000, Proteintech, Wuhan, China), anti-PHYKPL antibody (1:1000, Proteintech), anti-PXMP2 antibody (1:1000, Abcam, Cambridge, UK), anti-XYLB antibody (1: 1000, Thermo Fisher Scientific), and anti-GAPDH antibody (1:1000, Proteintech) at 4 ℃ overnight. Thereafter, the membranes were incubated with the HRP-conjugated secondary antibodies (1:5000; Jackson ImmunoResearch, Lancaster, PA, USA) at 37 ℃ for 1 h. After washing, protein bands were visualized using an enhanced chemiluminescence assay kit (Thermo Fisher Scientific).

Statistical analysis

Each experiment was repeated three times, and data were expressed as the mean ± standard deviation (SD). Statistical analysis was performed using GraphPad Prism 5 (San Diego, CA, USA), and Student’s t-test was used to compare the differences between normal osteoblasts and osteosarcoma cells. Statistical significance was set at p < 0.05. level.

Results

Identification of DE RNAs between the metastatic and non-metastatic samples

A total of 18,497 mRNAs, and 2528 lncRNAs were identified based on the RefSeq ID information provided by the downloaded data. The 176 osteosarcoma samples in the TCGA training set were divided into metastatic (56 samples) and non-metastatic (120 samples) groups. Based on the thresholds of FDR < 0.05 and |log2 FC|> 0.5, 432 DE RNAs, including 79 DE lncRNAs and 353 DE mRNAs were obtained (Fig. 1A). The bidirectional hierarchical clustering of the identified DE RNAs showed that these DE RNAs could significantly differentiate metastatic from non-metastatic samples (Fig. 1B).

Fig. 1
figure 1

Identification of differentially expressed (DE) RNAs and metabolic-related DE mRNAs. A. Effect size (log2FC)–log10 (false discovery rate) volcano plots. Red and blue dots represent DE RNAs, horizontal dashed lines represent false discovery rates < 0.05, two vertical dashed lines represent |log2 fold change (FC)|> 0.5, and the size of each point represents the absolute value of logFC. B. Bidirectional hierarchical clustering heatmap based on upregulated and downregulated DE-RNAs (top50). C. Venn diagram of the DE mRNAs and metabolic genes

Screening of metabolic-related DE mRNAs

Genes related to the metabolism of the amino acids, carbohydrates, and lipids were downloaded from the GSEA database, and 372 genes associated with the metabolism of amino acids and derivatives, 293 associated with the metabolism of carbohydrates, and 738 related to the metabolism of lipids were identified. After comparison with the DE RNAs, 107 overlapping genes were selected as metabolism-related DE mRNAs (Fig. 1C).

Construction and assessment of the proposed prognostic model

Based on univariate Cox regression analysis, 44 DE RNAs significantly correlated with disease prognosis were screened, including 20 DE lncRNAs and 24 DE mRNAs. Further, multivariate Cox regression analysis identified 13 independent prognostic DE RNAs (six lncRNAs and seven mRNAs). The Cox proportional hazards model based on the L1-penalized regularized regression algorithm in the penalized package was used to screen the optimized DE RNAs combinations, and nine DE RNAs (LINC00545, LINC01537, FOXC2-AS1, CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB) were the optimal combinations (Table 2).

Table 2 The optimized signature differentially expressed RNAs combination

Based on the LASSO prognostic coefficients of the nine optimized DE RNAs combinations and the expression levels of DE RNAs in the TCGA training set, a PS model was constructed. The PSs of the samples in the training dataset were calculated. All samples in the TCGA training dataset (Fig. 2A) and the three verification datasets (GSE16091, GSE21257, and GSE39055) (Fig. 2B–D) were divided into high- and low-risk score groups. There was a significant correlation between the high and low grouping and the actual disease prognosis information evaluated using the survival package (version 2.41–1 KM curve in R3.6.1), as shown in Fig. 2.

Fig. 2
figure 2

Construction and assessment of the prognostic model. Kaplan–Meier (KM) curves and receiver operating characteristic (ROC) curves for A TCGA, B GSE16091, C GSE21257, and D GSE39055. Left: KM curves based on the prognosis score prediction model. Green and red curves represent low- and high-risk samples, respectively. Right panel: ROC curves based on the prognostic model. The figures in parentheses indicate the specificity and sensitivity of the ROC curve

Independent prognostic clinical factors associated with survival and prognosis

Age, tumor metastasis, and PS model status were independent prognostic factors (Table 3). The prognosis-related KM curves for age and tumor metastasis are shown in Fig. 4A and B. The prognosis of young patients and those without metastatic osteosarcoma was better, which is consistent with the actual situation.

Table 3 Independent prognostic clinical factors screening

Then according to the age information, the sample was divided into ≤ 65 years old and > 65 years old groups for correlation analysis of risk prognosis model. The results showed that in the two age groups, there was a significant correlation between the different risk groups based on the PS model prediction and the actual prognosis. Besides, the correlation was more pronounced in the age group ≤ 65 years old, which indicated, to some extent, that the PS model was more beneficial to the accuracy of prediction of survival and prognosis of young patients (Fig. 3A).

Fig. 3
figure 3

Screening of independent prognostic clinical factors. A. The KM curves of age. Left: age in the TCGA prognostic KM curve. Green and red represent the KM curves of the sample group below or above 65 years of age, respectively. Middle and right: KM curves related to prognosis based on the prognosis score prediction model for the sample groups aged under and over 65 years old. Green and red curves represent low- and high-risk samples, respectively. B. The KM curves of metastasis. Left: tumor metastatic factors in TCGA sample prognostic KM curves. Green and red represent the KM curves of the sample group below or above 65 years of age, respectively. Middle and right: KM curves related to prognosis based on the prognosis score prediction model for the sample groups with and without tumor metastasis. Green and red curves represent low- and high-risk samples, respectively

Based on tumor metastasis information, the samples were divided into metastatic and non-metastatic osteosarcoma groups. Correlation analysis of the risk prognosis model in the two sample groups showed a significant correlation between the different risk groups based on the PS model prediction and actual prognosis (Fig. 3B).

Establishment and assessment of the nomogram survival models

To further analyze the correlation between age, tumor metastasis, PS model status factors, and survival prognosis, a nomogram survival rate model of the TCGA samples was constructed and analyzed (Fig. 4A). The “total points” axis of the nomogram combined various clinical indicators to predict the survival of the samples. The consistency between the 3- and 5-year survival rates predicted by the nomogram survival rate model and the actual 3- and 5-year survival rates was analyzed and verified, as shown in Fig. 4B. The 3- and 5-year C-index indices were 0.772 and 0.823, respectively.

Fig. 4
figure 4

Establishment and evaluation of nomogram models. A Nomogram of the survival prediction model based on the nomogram of independent prognostic factors. B Consistent plot of 3- and 5-year survival rates predicted using actual survival rates. The horizontal axis represents the predicted osteosarcoma rate, the vertical axis represents the actual osteosarcoma rate, and black and red represent the 3- and 5-year predicted line charts, respectively. C ROC curve of the prognostic model based on various factors in the TCGA training set (left) and validation set GSE21257 (right)

Moreover, based on the two independent prognostic factors obtained, a prognostic model of clinical factors was constructed separately, which was then compared with the previously constructed PS prognostic model based on the TCGA training set and GSE21257 dataset (there were two factors of age and tumor metastatic in this dataset, which could be used to verify multiple model results). The ROC curves for each type of model are shown in Fig. 4C. The parameters of the ROC curve are listed in Table 4.

Table 4 ROC curve parameter of each model

Co-expression of lncRNAs-mRNAs, and functional analysis

For the three DE lncRNAs used for model construction, combined with the 107 DE mRNAs related to the metabolism of the three types of substrates, the Pearson correlation coefficients between the three signature lncRNAs and the 107 DE mRNAs were calculated based on their expression levels in the TCGA training set. Connection pairs with significantly positive correlation coefficients were retained and 169 connection pairs were obtained. A co-expression network of the significant DE signature lncRNAs and metabolic genes was constructed by expressing this correlation, as shown in Fig. 5A. The network contained 91 nodes, including 3 lncRNAs and 88 metabolism-related DE mRNAs. Functional and pathway enrichment analyses were performed for DE mRNAs with a positive correlation with the signature lncRNAs contained in the co-expression network. Nineteen biological processes (leukotriene biosynthetic process, cyclooxygenase pathway, oxidation–reduction process, etc.) and seven KEGG signaling pathways (metabolic pathways, arachidonic acid metabolism, and glycerophospholipid metabolism) were screened (Fig. 5B).

Fig. 5
figure 5

Co-expression of lncRNAs and metabolism-related genes and functional analysis. A Co-expression network of signature lncRNAs and metabolism-related genes. Squares and rounds represent lncRNA and mRNA, respectively; green to red represent the change in the degree of differential expression from negative to positive. Oversized circular nodes are the signature DE genes used to build the model. B Biological processes of GO terms and KEGG signaling pathways enriched by genes in the co-expression network. The horizontal axis represents the number of genes and the vertical axis represents the item name. The size of the point represents the number of genes and the color represents significance

Validation of the signature genes

The expression levels of the nine signatures, including three lncRNAs (LINC00545, FOXC2-AS1, and LINC01537) and six mRNAs (CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB) in osteosarcoma cells were determined by qRT-PCR and western blotting. We found that the levels of LINC00545 and LINC01537 were significantly decreased, while the FOXC2-AS1 level was evidently increased in the osteosarcoma cells (U2OS and Saos2) compared with the osteoblasts hFOB1.19 (p < 0.05, Fig. 6A), which was consistent with the bioinformatics analysis. In addition, compared to osteoblast hFOB1.19, the expression of PFKFB4, PXMP2, and XYLB was significantly upregulated, whereas the expression of CYP27B1 and PHYKPL was downregulated in osteosarcoma cells (U2OS and Saos2) (p < 0.05, Fig. 6B). However, there was no significant difference in PHKG1 expression between osteoblasts hFOB1.19 and osteosarcoma cells (U2OS and Saos2) (p > 0.05, Fig. 6B). Finally, the protein expression levels of PFKFB4, PHYKPL, PXMP2, and XYLB in osteosarcoma cells were determined by western blotting. The expression of PFKFB4, PHYKPL, PXMP2, and XYLB proteins in hFOB1.19, U2OS, and Saos2 cells was in accordance with their mRNA expression determined by qRT-PCR (Fig. 6C). These results indicate that the concordance rate between the in silico analysis and qRT-PCR/western blotting was 88.89% (8/9), suggesting the robustness of our analysis.

Fig. 6
figure 6

Validation of signature genes in osteosarcoma cells. A The expression levels of three lncRNAs (LINC00545, FOXC2-AS1, and LINC01537) were determined by qRT-PCR. B The expression levels of six mRNAs (CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB) were measured by qRT-PCR. C The protein expression levels of PFKFB4, PHYKPL, PXMP2, and XYLB in osteosarcoma cells were detected by western blotting. * p < 0.05, ** p < 0.01, compared with the osteoblasts hFOB1.19

Discussion

A growing body of evidence shows that the occurrence, progression, and prognosis of osteosarcoma are correlated with metabolites of amino acids and their derivatives, carbohydrates, and lipids [39, 40]. In this study, we first identified 107 metabolism-related DE mRNAs and successfully proposed a PS model based on nine genes (LINC00545, LINC01537, FOXC2-AS1, CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB). A performance evaluation of this PS model based on the TCGA training dataset and three verification datasets (GSE16091, GSE21257, and GSE39055) showed that the PS model had high specificity and sensitivity. In addition, two independent prognostic clinical factors (age and tumor metastasis) were identified. Subsequently, a co-expression network, including three lncRNAs and 88 metabolism-related DE mRNAs, was constructed, and functional analysis showed that the co-expressed genes LINC00545 and LINC01537 were enriched in the cyclooxygenase pathway.

lncRNAs, a class of non-coding RNAs longer than 200 nt, are involved in multiple physiological and pathological processes [41]; they have been reported to play critical roles in tumorigenesis, progression, metastasis and prognosis [42]. In this study, three signature lncRNAs (LINC00545, LINC01537, and FOXC2-AS1) were selected to construct a PS model, suggesting that they might play critical roles in osteosarcoma prognosis. Recent studies have demonstrated that some natural antisense lncRNAs play key roles in regulating cancer biology [43, 44]. FOXC2-AS1 is a single antisense oligonucleotide RNA transcribed from the negative strand of FOXC2 (forkhead box C2), which facilitates the proliferation of many human cancers [45, 46]. Zhang et al. [47] demonstrated that the lncRNA FOXC2-AS1 was positively upregulated in adriamycin-resistant osteosarcoma cell lines and tissues, was associated with poor prognosis of osteosarcoma, and promoted adriamycin resistance in osteosarcoma cells in vitro and in vivo. Furthermore, the co-expression genes of LINC00545 and LINC01537, such as PTGIS, PTGES and PTGS1, were associated with cyclooxygenase pathway function. Cyclooxygenase, an enzyme that catalyzes the conversion of arachidonic acid to prostaglandin, is not only a key enzyme in arachidonic acid metabolism but is also involved in the regulation of cell proliferation, apoptosis, and angiogenesis [48]. Importantly, it has been suggested that cyclooxygenase-2 expression is correlated with the progression of some human malignancies, including osteosarcoma [49, 50]. Wang et al. [50] recently reported that cyclooxygenase-2 expression may correlate with metastasis and poor prognosis in osteosarcoma. Taken together, we speculate that LINC00545, LINC01537, and FOXC2-AS1 may play important roles in osteosarcoma occurrence and development, and that LINC00545 and LINC01537 may be associated with osteosarcoma prognosis via the cyclooxygenase pathway.

Since metabolic differences between tumor cells and adjacent normal cells were first reported in 1956, metabolic reprogramming has become a popular topic in tumor biology research [51]. The present study focused on three major metabolism-related genes and identified six signature metabolism-related genes (CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB) as osteosarcoma biomarkers. qRT-PCR and western blotting showed the downregulation of CYP27B1 and PHYKPL and the upregulation of PFKFB4, PXMP2, and XYLB in osteosarcoma cells. CYP27B1 (cytochrome P450 family 27, subfamily B, polypeptide 1) is a member of the cytochrome P450 superfamily that plays critical roles in xenobiotic metabolism and calcium homeostasis [52]. CYP27B1 level is upregulated in breast tumors compared to that in normal tissues [53]. PFKFB4 (6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase-4) encodes a bifunctional metabolic enzyme that synthesizes important sugar phosphate metabolites stimulating glycolysis [54]. PFKFB4 is a molecular fulcrum that activates transcriptional activation coupled with glucose metabolism by stimulating steroid receptor coactivator-3, which is critical for promoting metastatic tumors [55]. PHKG1, (phosphorylase kinase subunit G1) is the catalytic gamma subunit of glycogen phosphorylase kinase. As the only enzyme known to catalyze glycogen phosphorylase activation, glycogen phosphorylase kinase catalyzes the rate-limiting step in glycogen breakdown. It can regulate glycogenolysis by activating glycogen phosphorylase, thus leading to ATP production via glycolysis [56]. This extra burst of energy allows the tumor cells to perform high-energy tasks. PHYKPL (5-phosphohydroxy-L-lysine phospho-lyase), XYLB (xylulokinase) and CYP27B1 are involved in various metabolic pathways. Metabolic reprogramming is a hallmark of cancer [57]. Oncogene activation and loss of tumor suppressors promote the metabolic reprogramming of cancer cells, thereby enhancing nutrient uptake to provide energy and to supply biosynthetic pathways [58]. PXMP2 (peroxisomal membrane protein 2) is associated with inflammatory responses and bone remodeling [59]. These reports, together with our findings, suggest that CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB may be biomarkers of osteosarcoma progression and prognosis, and that the PS model constructed from the nine signatures (three lncRNAs and six mRNAs) can be effectively applied for risk classification and prognosis prediction of osteosarcoma.

Furthermore, in addition to genetic biomarkers, two independent prognostic clinical factors (age and tumor metastasis) were identified in our study. In many malignancies, age at diagnosis is a well-known prognostic factor. A previous study demonstrated that younger patients with osteosarcoma have significantly better outcomes than older patients [60]. Regarding the other factor (tumor metastasis), it has been reported that the prognosis for patients with metastatic osteosarcoma is poor compared to that for non-metastatic tumors [61]. Therefore, we concluded that age and tumor metastasis are closely associated with the development and prognosis of osteosarcoma.

However, this study has some limitations. First, the effects of the identified signature lncRNAs on metabolism require further exploration, and the specific roles and mechanisms of the screened signatures in osteosarcoma should be further investigated using in vitro and in vivo assays. In addition, the performance of the proposed PS model based on nine signatures should be evaluated and validated clinically.

Conclusions

In conclusion, we successfully constructed a PS model based on the expression levels of LINC00545, LINC01537, FOXC2-AS1, CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB. This PS model could accurately predict the prognosis of patients with osteosarcoma. In addition, CYP27B1, PFKFB4, PHKG1, PHYKPL, PXMP2, and XYLB may serve as metabolism-related biomarkers of osteosarcoma. These findings help expand our understanding of metabolism-related genes associated with osteosarcoma and provide a theoretical basis for the proposed PS model based on nine signatures to aid in the diagnosis, therapy, and prognostic assessment of osteosarcoma patients.

Data availability

The datasets during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Abbreviations

DE:

Differentially expressed

FC:

Fold change

FDR:

False discovery rate

GO:

Gene ontology

GSEA :

Gene Set Enrichment Analysis

KEGG:

Kyoto Encyclopedia Genes and Genomes

KM:

Kaplan–Meier

LASSO:

Least absolute shrinkage and selection operator

lncRNA:

Long non-coding RNA

PS:

Prognostic score

ROC:

Receiver-operating characteristic

TCGA:

The Cancer Genome Atlas

References

  1. Rickel K, Fang F, Tao J. Molecular genetics of osteosarcoma. Bone. 2017;102:69–79.

    Article  CAS  PubMed  Google Scholar 

  2. Mirabello L, Troisi RJ, Savage SA. Osteosarcoma incidence and survival rates from 1973 to 2004: data from the Surveillance, Epidemiology, and End Results Program. Cancer. 2009;115(7):1531–43.

    Article  PubMed  Google Scholar 

  3. Flores RJ, Kelly AJ, Li Y, Nakka M, Barkauskas DA, Krailo M, et al. A novel prognostic model for osteosarcoma using circulating CXCL 10 and FLT 3 LG. Cancer. 2017;123(1):144–54.

    Article  CAS  PubMed  Google Scholar 

  4. Moore DD, Luu HH. Osteosarcoma. Orthopaedic oncology. Springer; 2014. p. 65–92.

  5. Ballman KV. Biomarker: predictive or prognostic? J Clin Oncol. 2015;33(33):3968–71.

    Article  CAS  PubMed  Google Scholar 

  6. Wang T, Ji F, Dai Z, Xie Y, Yuan D. Increased expression of microRNA-191 as a potential serum biomarker for diagnosis and prognosis in human osteosarcoma. Cancer Biomark. 2015;15(5):543–50.

    Article  CAS  PubMed  Google Scholar 

  7. Li F, Cao L, Hang D, Wang F, Wang Q. Long non-coding RNA HOTTIP is up-regulated and associated with poor prognosis in patients with osteosarcoma. Int J Clin Exp Pathol. 2015;8(9):11414.

    PubMed  PubMed Central  Google Scholar 

  8. Liu Y, Teng Z, Wang Y, Gao P, Chen J. Prognostic significance of survivin expression in osteosarcoma patients: a meta-analysis. Med Sci Monit. 2015;21:2877.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Min L, Choy E, Tu C, Hornicek F, Duan Z. Application of metabolomics in sarcoma: From biomarkers to therapeutic targets. Crit Rev Oncol Hematol. 2017;116:1.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Fritsche-Guenther R, Zasada C, Mastrobuoni G, Royla N, Rainer R, Roßner F, et al. Alterations of mTOR signaling impact metabolic stress resistance in colorectal carcinomas with BRAF and KRAS mutations. Sci Rep. 2018;8(1):1–17.

    Article  CAS  Google Scholar 

  11. Ren L, Hong ES, Mendoza A, Issaq S, Hoang CT, Lizardo M, et al. Metabolomics uncovers a link between inositol metabolism and osteosarcoma metastasis. Oncotarget. 2017;8(24):38541.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Li W, Zhang H, Assaraf YG, Zhao K, Xu X, Xie J, et al. Overcoming ABC transporter-mediated multidrug resistance: Molecular mechanisms and novel therapeutic drug strategies. Drug Resist Updates. 2016;27:14–29.

    Article  CAS  Google Scholar 

  13. Gorska-Ponikowska M, Perricone U, Kuban-Jankowska A, Lo Bosco G, Barone G. 2-methoxyestradiol impacts on amino acids-mediated metabolic reprogramming in osteosarcoma cells by its interaction with NMDA receptor. J Cell Physiol. 2017;232(11):3030–49. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/jcp.25888.

    Article  CAS  PubMed  Google Scholar 

  14. Zhu R, Li X, Ma Y. miR-23b-3p suppressing PGC1α promotes proliferation through reprogramming metabolism in osteosarcoma. Cell Death Dis. 2019;10(6):381. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41419-019-1614-1.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Huang W, Xiao Y, Wang H, Chen G, Li K. Identification of risk model based on glycolysis-related genes in the metastasis of osteosarcoma. Front Endocrinol (Lausanne). 2022;13:1047433. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fendo.2022.1047433.

    Article  PubMed  Google Scholar 

  16. Cao Y. Adipocyte and lipid metabolism in cancer drug resistance. J Clin Investig. 2019;129(8):3006–17. https://doiorg.publicaciones.saludcastillayleon.es/10.1172/jci127201.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Qian H, Lei T, Hu Y, Lei P. Expression of Lipid-Metabolism Genes Is Correlated With Immune Microenvironment and Predicts Prognosis in Osteosarcoma. Front Cell Dev Biol. 2021;9: 673827. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fcell.2021.673827.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Bravo-Ruiz I, Medina M, Martínez-Poveda B. From Food to Genes: Transcriptional Regulation of Metabolism by Lipids and Carbohydrates. Nutrients. 2021;13(5). https://doiorg.publicaciones.saludcastillayleon.es/10.3390/nu13051513.

  19. Oermann EK, Wu J, Guan K-L, Xiong Y. Alterations of metabolic genes and metabolites in cancer. Seminars in cell & developmental biology: Elsevier; 2012. p. 370–80.

  20. Buddingh EP, Kuijjer ML, Duim RA, Bürger H, Agelopoulos K, Myklebost O, et al. Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents. Clin Cancer Res. 2011;17(8):2110–9.

    Article  CAS  PubMed  Google Scholar 

  21. Kelly AD, Haibe-Kains B, Janeway KA, Hill KE, Howe E, Goldsmith J, et al. MicroRNA paraffin-based studies in osteosarcoma reveal reproducible independent prognostic profiles at 14q32. Genome Med. 2013;5(1).

  22. Paoloni M, Davis S, Lana S, Withrow S, Sangiorgi L, Picci P, et al. Canine tumor cross-species genomics uncovers targets linked to osteosarcoma progression. BMC Genomics. 2009;10(625):1471–2164.

    Google Scholar 

  23. Wright MW. A short guide to long non-coding RNA gene nomenclature. Hum Genomics. 2014;8(1):7. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/1479-7364-8-7.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):20.

    Article  Google Scholar 

  25. Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A. 1998;95(25):14863–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Wang L, Cao C, Ma Q, Zeng Q, Wang H, Cheng Z, et al. RNA-seq analyses of multiple meristems of soybean: novel and alternative transcripts, evolutionary and functional implications. BMC Plant Biol. 2014;14(169):1471–2229.

    Google Scholar 

  27. 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 U S A. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Wang P, Wang Y, Hang B, Zou X, Mao JH. A novel gene expression-based prognostic scoring system to predict survival in gastric cancer. Oncotarget. 2016;7(34):55343–51.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997;16(4):385–95.

    Article  CAS  PubMed  Google Scholar 

  30. Goeman JJ. L1 penalized estimation in the Cox proportional hazards model. Biom J. 2010;52(1):70–84. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/bimj.200900028.

    Article  PubMed  Google Scholar 

  31. Anderson WI, Schlafer DH, Vesely KR. Thyroid follicular carcinoma with pulmonary metastases in a beaver (Castor canadensis). J Wildl Dis. 1989;25(4):599–600.

    Article  CAS  PubMed  Google Scholar 

  32. Eng KH, Schiller E, Morrell K. On representing the prognostic value of continuous gene expression biomarkers with the restricted mean survival curve. Oncotarget. 2015;6(34):36308–18.

    Article  PubMed  PubMed Central  Google Scholar 

  33. 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):1471–2105.

    Google Scholar 

  34. Mayr A, Schmid M. Boosting the concordance index for survival data--a unified framework to derive and evaluate biomarker combinations. Plos One. 2014;9(1).

  35. Shan S, Chen W, Jia JD. Transcriptome analysis revealed a highly connected gene module associated with cirrhosis to hepatocellular carcinoma development. Front Genet. 2019;10:305.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. 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(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  PubMed  Google Scholar 

  38. da Huang W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.

    Article  PubMed  Google Scholar 

  39. Li Z, Jia M, Wu X, Cui J, Pan A, Li L. Overexpression of Trps1 contributes to tumor angiogenesis and poor prognosis of human osteosarcoma. Diagn Pathol. 2015;10(1):167.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Ren L, Khanna C. Role of ezrin in osteosarcoma metastasis. Current Advances in Osteosarcoma. Springer; 2014. p. 181–201.

  41. Majidinia M, Yousefi B. Long non-coding RNAs in cancer drug resistance development. DNA Repair. 2016;45:25–33.

    Article  CAS  PubMed  Google Scholar 

  42. Zhang J, Zhang B, Wang T, Wang H. LncRNA MALAT1 overexpression is an unfavorable prognostic factor in human cancer: evidence from a meta-analysis. Int J Clin Exp Med. 2015;8(4):5499.

    PubMed  PubMed Central  Google Scholar 

  43. Jadaliha M, Gholamalamdari O, Tang W, Zhang Y, Petracovici A, Hao Q, et al. A natural antisense lncRNA controls breast cancer progression by promoting tumor suppressor gene mRNA stability. PLoS Genet. 2018;14(11):e1007802.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Deng S-J, Chen H-Y, Zeng Z, Deng S, Zhu S, Ye Z, et al. Nutrient stress–dysregulated antisense lncRNA GLS-AS impairs GLS-mediated metabolism and represses pancreatic cancer progression. Can Res. 2019;79(7):1398–412.

    Article  CAS  Google Scholar 

  45. Chen Y, Gu M, Liu C, Wan X, Shi Q, Chen Q, et al. Long noncoding RNA FOXC2-AS1 facilitates the proliferation and progression of prostate cancer via targeting miR-1253/EZH2. Gene. 2019;686:37–42.

    Article  CAS  PubMed  Google Scholar 

  46. Yang H, Chen T, Xu S, Zhang S, Zhang M. Long noncoding RNA FOXC2-AS1 predicts poor survival in breast cancer patients and promotes cell proliferation. Oncol Res. 2019;27(2):219–26.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Zhang CL, Zhu KP, Ma XL. Antisense lncRNA FOXC2-AS1 promotes doxorubicin resistance in osteosarcoma by increasing the expression of FOXC2. Cancer Lett. 2017;396:66–75. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.canlet.2017.03.018.

    Article  CAS  PubMed  Google Scholar 

  48. Lee EJ, Choi EM, Kim SR, Park JH, Kim H, Ha KS, et al. Cyclooxygenase-2 promotes cell proliferation, migration and invasion in U2OS human osteosarcoma cells. Exp Mol Med. 2007;39(4):469–76.

    Article  CAS  PubMed  Google Scholar 

  49. Sicking I, Rommens K, Battista MJ, Böhm D, Gebhard S, Lebrecht A, et al. Prognostic influence of cyclooxygenase-2 protein and mRNA expression in node-negative breast cancer patients. BMC Cancer. 2014;14(1):952.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Wang S, Gao H, Zuo J, Gao Z. Cyclooxygenase-2 expression correlates with development, progression, metastasis, and prognosis of osteosarcoma: a meta-analysis and trial sequential analysis. FEBS Open Bio. 2019;9(2):226–40.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Warburg O. On the origin of cancer cells. Science. 1956;123(3191):309–14.

    Article  CAS  PubMed  Google Scholar 

  52. Dunbar DR, Khaled H, Evans LC, Al-Dujaili EA, Mullins LJ, Mullins JJ, et al. Transcriptional and physiological responses to chronic ACTH treatment by the mouse kidney. Physiol Genomics. 2010;40(3):158–66.

    Article  CAS  PubMed  Google Scholar 

  53. Lopes N, Sousa B, Martins D, Gomes M, Vieira D, Veronese LA, et al. Alterations in Vitamin D signalling and metabolic pathways in breast cancer progression: a study of VDR, CYP27B1 and CYP24A1 expression in benign and malignant breast lesions Vitamin D pathways unbalanced in breast lesions. BMC Cancer. 2010;10(1):483.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Pilkis SJ, Claus TH, Kurland IJ, Lange AJ. 6-Phosphofructo-2-kinase/fructose-2, 6-bisphosphatase: a metabolic signaling enzyme. Annu Rev Biochem. 1995;64(1):799–835.

    Article  CAS  PubMed  Google Scholar 

  55. Dasgupta S, Rajapakshe K, Zhu B, Nikolai BC, Yi P, Putluri N, et al. Metabolic enzyme PFKFB4 activates transcriptional coactivator SRC-3 to drive breast cancer. Nature. 2018;556(7700):249–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Gu J, Qi Y, Lu Y, Tao Q, Yu D, Jiang C, et al. Lung adenocarcinoma-derived vWF promotes tumor metastasis by regulating PHKG1-mediated glycogen metabolism. Cancer Sci. 2022;113(4):1362–76. https://doiorg.publicaciones.saludcastillayleon.es/10.1111/cas.15298.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74.

    Article  CAS  PubMed  Google Scholar 

  58. Boroughs LK, DeBerardinis RJ. Metabolic pathways promoting cancer cell survival and growth. Nat Cell Biol. 2015;17(4):351–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Montalvany-Antonucci CC, Zicker MC, Ferreira AVM, Macari S, Ramos-Junior ES, Gomez RS, et al. High-fat diet disrupts bone remodeling by inducing local and systemic alterations. J Nutr Biochem. 2018;59:93–103. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.jnutbio.2018.06.006.

    Article  CAS  PubMed  Google Scholar 

  60. Hagleitner MM, Hoogerbrugge PM, van der Graaf WT, Flucke U, Schreuder HB, te Loo DMW. Age as prognostic factor in patients with osteosarcoma. Bone. 2011;49(6):1173–7.

    Article  PubMed  Google Scholar 

  61. Mialou V, Philip T, Kalifa C, Perol D, Gentet JC, Marec-Berard P, et al. Metastatic osteosarcoma at diagnosis: prognostic factors and long-term outcome the-French pediatric experience. Cancer. 2005;104(5):1100–9.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

None.

Funding

None.

Author information

Authors and Affiliations

Authors

Contributions

TLW carried out the conception and design of the research, participated in the acquisition of data. XYW carried out the analysis and interpretation of data. TLW performed the statistical analysis. TLW drafted the manuscript. XYW revised the manuscript for important intellectual content. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Xingyi Wu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wu, T., Wu, X. Construction and evaluation of a prognostic model based on the expression of the metabolism-related signatures in patients with osteosarcoma. BMC Musculoskelet Disord 26, 303 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12891-025-08439-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12891-025-08439-9

Keywords