- Research
- Open access
- Published:
Construction and evaluation of a prognostic model based on the expression of the metabolism-related signatures in patients with osteosarcoma
BMC Musculoskeletal Disorders volume 26, Article number: 303 (2025)
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.
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.
-
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.
-
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.
-
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:
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 2–△△Ct method.
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).
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).
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.
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.
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).
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.
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.
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).
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.
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
Rickel K, Fang F, Tao J. Molecular genetics of osteosarcoma. Bone. 2017;102:69–79.
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.
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.
Moore DD, Luu HH. Osteosarcoma. Orthopaedic oncology. Springer; 2014. p. 65–92.
Ballman KV. Biomarker: predictive or prognostic? J Clin Oncol. 2015;33(33):3968–71.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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).
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.
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.
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.
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.
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.
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.
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.
Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997;16(4):385–95.
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.
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.
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.
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.
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).
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.
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.
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.
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.
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.
Ren L, Khanna C. Role of ezrin in osteosarcoma metastasis. Current Advances in Osteosarcoma. Springer; 2014. p. 181–201.
Majidinia M, Yousefi B. Long non-coding RNAs in cancer drug resistance development. DNA Repair. 2016;45:25–33.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Warburg O. On the origin of cancer cells. Science. 1956;123(3191):309–14.
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.
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.
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.
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.
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.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74.
Boroughs LK, DeBerardinis RJ. Metabolic pathways promoting cancer cell survival and growth. Nat Cell Biol. 2015;17(4):351–9.
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.
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.
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.
Acknowledgements
None.
Funding
None.
Author information
Authors and Affiliations
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
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/.
About this article
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
Received:
Accepted:
Published:
DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12891-025-08439-9