Prognostic Signatures of Metabolic Genes and Metabolism-Related Long Non-coding RNAs Accurately Predict Overall Survival for Osteosarcoma Patients

In this study, we identified eight survival-related metabolic genes in differentially expressed metabolic genes by univariate Cox regression analysis based on the therapeutically applicable research to generate effective treatments (n = 84) data set and genotype tissue expression data set (n = 396). We also constructed a six metabolic gene signature to predict the overall survival of osteosarcoma (OS) patients using least absolute shrinkage and selection operator (Lasso) Cox regression analysis. Our results show that the six metabolic gene signature showed good performance in predicting survival of OS patients and was also an independent prognostic factor. Stratified correlation analysis showed that the metabolic gene signature accurately predicted survival outcomes in high-risk and low-risk OS patients. The six metabolic gene signature was also verified to perform well in predicting survival of OS patients in an independent cohort (GSE21257). Then, using univariate Cox regression and Lasso Cox regression analyses, we identified an eight metabolism-related long noncoding RNA (lncRNA) signature that accurately predicts overall survival of OS patients. Gene set variation analysis showed that the apical surface and bile acid metabolism, epithelial mesenchymal transition, and P53 pathway were activated in the high-risk group based on the eight metabolism-related lncRNA signature. Furthermore, we constructed a competing endogenous RNA (ceRNA) network and conducted immunization score analysis based on the eight metabolism-related lncRNA signature. These results showed that the six metabolic gene signature and eight metabolism-related lncRNA signature have good performance in predicting the survival outcomes of OS patients.

Introduction

Osteosarcoma (OS) is a primary malignant bone cancer and commonly occurs in adolescents and children. The overall annual incidence of OS is 3.4 million worldwide (Mirabello et al., 2009; Pingping et al., 2019; Czarnecka et al., 2020; Mirabello et al., 2020). As a malignant tumor, OS typically occurs in the metaphysis of the long bones, such as the distal femur (43%), proximal tibia (23%), or humerus (10%) (Isakoff et al., 2015). Previous reports suggest that the 5-year survival rate of patients with nonmetastatic OS is 70–75%, but the long-term survival rate of metastatic OS patients is only 30% (Anwar et al., 2020). In addition, multidrug resistance is a major challenge in OS treatment (Strauss et al., 2010). Hence, there is an urgent need to identify novel targets and biomarkers for the diagnosis and prognosis of OS to improve the survival rate of OS patients.

Cancer metabolism is one of the oldest areas of research in cancer biology, and targeting metabolism has been an effective cancer treatment modality for decades (Cairns et al., 2011; Luengo et al., 2017). An increasing amount of evidence shows that changes in cell metabolism contribute to cancer development and progression (Vander Heiden and DeBerardinis, 2017; Kreuzaler et al., 2020). Additionally, several studies show that the tumor suppressor p53, MYC oncogene, pyruvate kinase isozymes M2 (PKM2), and hypoxia-inducible factor 1(HIF1) regulate cancer metabolism and are involved in the prognosis of cancers (Denko, 2008; Dayton et al., 2016; Gomes et al., 2018; Khan et al., 2020). Long noncoding RNAs (lncRNAs) are noncoding RNAs more than 200 nucleotides in length and play important roles in transcriptional regulation, epigenetic gene regulation, and disease (Mercer et al., 2009; Kumar and Goyal, 2017). The main difference between lncRNAs and mRNAs is that lncRNAs lack reading frames encoding proteins (Dinger et al., 2011). However, investigations using advanced molecular techniques suggest that some lncRNAs contain short open reading frames (sORFs) and can interact with ribosomes to encode proteins (Chen et al., 2020; Vergara et al., 2020). Recent studies show that lncRNAs play important roles in glucose, protein, lipid, and nucleic acid metabolism by directly or indirectly targeting enzymes and oncogenic signaling pathways (Fan et al., 2017; Zeng et al., 2018; Lu et al., 2019). Tang et al. (2019) found that lncRNA GLCC1 promotes carcinogenesis and glucose metabolism by stabilizing c-Myc, resulting in poor prognosis in colorectal cancer. Another study showed that lncRNA ANRIL is involved in regulating AML development by modulating the glucose metabolism pathway of AdipoR1/AMPK/SIRT1 (Sun et al., 2018). Although many metabolic biomarkers have been identified for the diagnosis and prognosis of human cancers, research on metabolic biomarkers for the prognosis of OS is limited. Therefore, in this study, we sought to identify novel metabolic signatures that are related to the diagnosis and prognosis of OS patients.

In this study, we identified a six metabolic gene signature using bioinformatics that shows good performance in predicting survival of OS patients. Stratified correlation analysis shows that the metabolic gene signature accurately predicted survival outcomes in both high- and low-risk OS patients. The six metabolic gene signature was also verified to perform well in predicting survival of OS patients in an independent data set (GSE21257). Furthermore, we also identified an eight metabolism-related lncRNA signature that is related to overall survival in OS and shows good performance in predicting overall survival of OS. Gene set variation analysis (GSVA) revealed that multiple metabolic processes and signaling pathways were significantly enriched in the high-risk groups. Finally, we constructed a competing endogenous RNA (ceRNA) network based on the eight metabolism-related lncRNA signature and analyzed the immunization scores in high- and low-risk OS patients. Our results suggest that the six metabolic gene signature and eight metabolism-related lncRNA signature identified show robust performance in predicting the survival outcomes of OS patients.

Materials and Methods

Data Downloading and Processing

A total of 944 metabolic genes were extracted from metabolism-related Kyoto encyclopedia of genes and genomes (KEGG) pathways. RNA sequence data and clinical information on the OS patients were downloaded from the therapeutically applicable research to generate effective treatments (TARGET¹). Transcriptome data of normal human tissues were downloaded from the genotype tissue expression (GTEx²) database. Expression profiling and clinical information on a test independent data set (GSE21257) was downloaded from the Gene Expression Omnibus (GEO³). The clinical characteristics of all OS patients in the two data sets are listed in Table 1. The R software⁴ sva (Li et al., 2018) package was used to merge the raw data of the two sets (TARGET OS n = 84 and GTEx OS n = 396) and eliminate batch-to-batch differences.

TABLE 1

Table 1. Clinical data for all patients.

Function Annotation

A total of 2,282 differentially expressed genes were selected based on a cutoff value of | log₂FC | > 1 and a P value less than 0.05 in OS tissue and normal muscle tissue using the Limma (Ritchie et al., 2015) package, and 64 differentially expressed genes were screened from the 2,282 differentially expressed genes. Next, we conducted gene ontology (GO) and KEGG enrichment analysis of the 64 differentially expressed metabolic genes using the clusterProfiler (Kanehisa et al., 2017) package, established a protein–protein interaction (PPI) network, and identified 10 hub metabolic genes using the Search Tool for the Retrieval of Interacting Genes online tool (STRING⁵) and Cytoscape (Shannon et al., 2003).

Identification and Construction of Prognostic Signature

Univariate Cox regression analysis was used to identify differentially expressed metabolic genes whose expression levels were significantly associated (P < 0.05) with overall survival of OS patients in the training data set (TARGET OS). Next, we conducted least absolute shrinkage and selection operator (Lasso) Cox regression analysis to identify metabolic genes related to OS prognosis using the glmnet (Friedman et al., 2010) package, and the OS patients were divided into high- (n = 42) and low-risk (n = 42) groups based on the median risk score. Finally, we constructed a six metabolic gene signature. To identify the metabolism-related lncRNAs, we performed Pearson correlation analysis between the lncRNAs and metabolic genes related to OS prognosis using | R| ≥ 0.4 and P < 0.05 as the selection criteria. Thereafter, we used the same method to screen a metabolism-related lncRNA signature from 147 metabolism-related lncRNAs.

Evaluation and Verification of the Prognostic Signature

The OS patients were classified into high-risk and low-risk groups based on their prognostic risk score using the median risk score in the metabolic gene signature and metabolism-related lncRNA signature. First, we compared the overall survival of the high-risk and low-risk groups of patients in the two kinds of signatures using the Kaplan–Meier survival curve. According to receiver operating characteristic (ROC) curves, we evaluated the diagnostic efficacy of each clinicopathological characteristic and the prognostic signature for OS patients. Furthermore, univariate and multivariate Cox regression analyses were performed to evaluate whether the risk score was independent of other clinical variables, such as age, gender, and metastasis, in the prognostic signature. We constructed a nomogram by integrating the traditional clinical variables, such as age, gender, and metastasis as well as the risk score derived from the prognostic signature of metabolic genes to analyze the probable 1, 3, and 5-year overall survival of the OS patients. To test whether the prognostic signature of the metabolic genes has a robust ability to predict patient survival in an independent data set, we verified it with an independent data set (GSE21257).

Gene Set Enrichment Analysis (GSEA) and GSVA

To assess the important functional phenotypes between the high- and low-risk groups based on the eight metabolism-related lncRNA signature, we performed GSEA and GSVA. We used “c5.all.v7.0.symbols.gmt” as the reference gene sets and performed GSEA enrichment analysis using GSEA software (version: 4.0.3). We used “h.all.v7.0.symbols.gmt” as the reference gene set for the GSVA analysis, and the adjusted p-value <0.05 was considered statistically significant.

ceRNA Network and Immunization Scores Analysis

Using startbase databases, three of eight metabolism-related lncRNAs were extracted to construct the ceRNA network. We then used three databases (miRTarBase, miRDB, and TargetScan) to search target mRNAs based on the 39 miRNAs. Next, we extracted 88 target mRNAs from the differentially expressed genes. Then, we obtained the immunization scores of immune cells and immune-related functions in the high- and low-risk groups using ssGSEA based on the eight metabolism-related lncRNA signature.

Results

Functional Annotation of Differentially Expressed Metabolic Genes

First, we used TARGET (n = 84) and GTEx (n = 396) data sets to compare OS and normal muscle tissues and identified 2,282 differentially expressed genes. We then identified 64 differentially expressed metabolic genes in the 2,282 differentially expressed genes, of which 12 metabolic genes were upregulated in normal tissues, and 52 metabolic genes were upregulated in OS tissues (Figures 1A,B). GO annotation showed that the 64 differentially expressed metabolic genes were involved in various biological functions, including the carboxylic acid biosynthetic process, organic acid biosynthetic process, DNA polymerase complex, ficolin-1-rich granule, coenzyme binding, and transferase activity. KEGG annotation revealed that the 64 differentially expressed metabolic genes were related to signaling pathways, including biosynthesis of amino acids, carbon metabolism, pentose phosphate pathway, and glycolysis/gluconeogenesis (Figure 1C). Last, we constructed a PPI network (Figure 1D) and obtained 10 differentially expressed metabolic genes using Cytoscape: PHGDH, GLUD2, PYGM, ALDOA, ALDOC, PFKM, FBP2, FBP1, GPD1, and GLUL (Figure 1E).

FIGURE 1

Figure 1. Functional annotation of differentially expressed metabolic genes. (A) Heat map of 64 differentially expressed metabolic genes. (B) Volcano plot of 944 metabolic genes. (C) GO and KEGG analysis of 64 differentially expressed metabolic genes. (D) PPI of 64 differentially expressed metabolic genes. (E) Ten hub differentially expressed metabolic genes.

Identification of Prognostically Significant Metabolic Genes in OS Patients

First, based on the 64 differentially expressed metabolic genes that were screened, we obtained clinical information for the 84 OS samples from the TARGET data set. Then, using univariate Cox regression analysis, we found that eight of the 64 differentially expressed metabolic genes were significantly correlated with overall survival of OS patients (P < 0.05; Figure 2A). These genes are the following: PYGM, CKMT2, NAT1, AADAT, FADS2, GPX7, PHOSPHO1, and CHST13.

FIGURE 2

Figure 2. Construction of metabolic gene signature. (A) Univariate Cox regression analysis showed that eight out of the 64 differentially expressed metabolic genes significantly correlated with the overall survival of OS patients (P < 0.05). (B,C) Lasso Cox regression analysis showing that six out of the eight metabolic genes were good candidates for constructing the prognostic signature. (D) Coefficients of six out of the eight metabolic genes. (E) Six metabolic genes signature constructed by iterative Lasso cox regression analysis. (F) ROC curve to evaluate 1, 3, and 5-year prediction efficiency of the six metabolic genes signature.

Construction and Evaluation of the Metabolic Genes Signature

To identify a metabolic gene signature for predicting overall survival of OS patients, we conducted Lasso Cox regression for the eight metabolic genes that correlated with the overall survival of OS patients. Lasso Cox regression analysis showed that six of the eight metabolic genes were good candidates for constructing a prognostic signature (Figures 2B–D). Next, we used iterative Lasso Cox regression analysis to construct an optimal prognostic signature of six metabolic genes composed of PYGM, CKMT2, NAT1, AADAT, GPX7, and CHST13 (Figure 2E). ROC analysis indicated that prediction efficiency of the 6 metabolic genes signature was pretty good (1-year AUC = 0.759, 3-year AUC = 0.775, and 5-year AUC = 0.755) and was robust in predicting OS prognosis (Figure 2F).

Based on the risk score of each OS patient in the TARGET data set, the patients were divided into high-risk (n = 42) and low-risk (n = 42) groups. Kaplan–Meier survival curve analysis showed that the overall survival of OS patients with high-risk scores was significantly shorter than those of patients with low-risk scores (Figure 3A). Then, based on the metabolic gene signature, we obtained the risk score distribution, survival status, and a heat map of the six metabolic genes (Figure 3B). Principal component analysis (PCA) based on the six metabolic genes showed two markedly different distribution patterns between high-risk and low-risk groups (Figure 3C). Next, using univariate Cox analyses, we found that the metastatic and metabolic gene prognostic risk scores were significantly associated with overall survival (P < 0.001, Figure 3D). Multivariate Cox analyses showed that the prognostic risk score for metabolic genes was significantly associated with overall survival (P < 0.001, Figure 3E). As shown in Figure 3F, the ROC curve analysis demonstrated that the AUC values of age, gender, and metastatis were less than the prognostic signature of the metabolic genes (AUC = 0.755). We further performed a stratification analysis to investigate the prognostic value of the metabolic genes (Figure 4A). Finally, we constructed a nomogram and calibration curve analysis to accurately estimate the 1, 3, and 5-year survival probabilities using the risk scores calculated from the metabolic gene prognostic signature and other clinicopathological factors, including age, gender, and metastatis (Figure 4B). These results suggest that the prognosis signature could accurately determine the prognosis of patients with OS.

FIGURE 3

Figure 3. Evaluation of the metabolic gene signature. (A) Kaplan–Meier curves of patients in the low-risk and high-risk groups based on the six metabolic genes signature. (B) Risk score distribution and survival status of OS patients based on the six metabolic genes signature and heat map of the 6 metabolic genes signature expression pattern. (C) PCA based on the confirmed six metabolic genes signature. (D) Univariate Cox regression analysis showing that the metastatic and metabolic genes signature risk scores are significantly associated with overall survival. (E) Multivariate Cox regression analysis showing that the metabolic genes signature risk score is an independent prognostic indicator for overall survival of OS patients. (F) ROC curve analysis with the prognostic accuracy of age, gender, and metastatis.

FIGURE 4

Figure 4. Analysis of prognosis signature with other clinical characteristics. (A) Stratification analysis to investigate the prognostic value of the six metabolic genes signature (age, metastatis, gender, histologic response). (B) Nomogram for predicting 1, 3, and 5-year survival rates of OS patients based on the six metabolic genes signature risk score and other clinical characteristics, such as age, gender, and metastatis. Calibration curve of the nomogram for 5-year survival rates.

Verification of Metabolic Genes Signature in an Independent Cohort

To further examine the prognostic value of the six metabolic genes, we verified the six metabolic genes signature in an independent cohort (GSE21257). First, based on the risk score for each OS patient in the GSE21257 data set, OS patients were divided into high-risk (n = 26) and low-risk (n = 27) groups. Kaplan–Meier survival curve analysis showed that the overall survival of OS patients with high-risk scores was significantly shorter than that of patients with low-risk scores (p = 1.469e-03, Figure 5A). Next, we obtained the risk score distribution, survival status, and a heat map of the six metabolic genes (Figure 5B) based on the metabolic gene prognosis signature. ROC revealed that the AUCs for 1-, 3-, and 5-year survival were 0.745, 0.781, and 0.819, respectively, in the independent cohort (GSE21257; Figure 5C). These results demonstrate that the six metabolic genes signature can also predict the survival of OS patients in other independent cohorts.

FIGURE 5

Figure 5. Verification of metabolic gene signature in an independent cohort. (A) Risk score distribution and survival status of OS patients based on the six metabolic genes signature and a heat map of the six metabolic genes signature expression pattern in one independent cohort (GSE21257). (B) Kaplan–Meier curves of patients in the low- and high-risk groups based on the 6 metabolic genes signature in an independent cohort (GSE21257). (C) ROC curve to evaluate 1, 3, and 5-year prediction efficiency of the six metabolic genes signature in an independent cohort (GSE21257).

Identification of Metabolism-Related lncRNAs

To explore metabolism-related lncRNAs, we first identified 148 metabolism-related lncRNAs by performing Pearson correlation analysis between the lncRNAs and the metabolism-related genes using | R| ≥ 0.4 and P < 0.05 as the selection criteria. Univariate Cox regression analysis showed that expression of 17 metabolism-related lncRNAs were significantly correlated with the overall survival of OS patients (P < 0.05; Table 2).

TABLE 2

Table 2. The 17 metabolism-related prognostic lncRNAs.

Construction and Evaluation of the Metabolism-Related lncRNA Signature

First, we performed Lasso Cox regression analysis for 17 metabolism-related lncRNAs correlated with the overall survival of OS patients. The analysis showed that eight of the 17 metabolism-related lncRNAs were good candidates for constructing the prognostic signature, including JMJD1C-AS1, FLJ45513, FAM155A-IT1, LINC00837, MANCR, LINC00506, CACNA1C-AS1, and IL10RB-DT (Figures 6A–D). ROC analysis suggested that the metabolism-related lncRNA signature has good performance in predicting overall survival of OS patients, and the AUCs for 1, 3, and 5-year survival were 0.813, 0.814, and 0.802, respectively (Figure 6E).

FIGURE 6

Figure 6. Construction of metabolism-related lncRNA signature. (A) Heat map of the correlations between metabolic genes and the 17 prognostic metabolism-related lncRNAs (**p < 0.01 and ***p < 0.001). (B,C) Lasso Cox regression analysis showed that eight of the 17 prognostic metabolism-related lncRNAs were good candidates for constructing the prognostic signature. (D) Coefficients of the eight metabolism-related lncRNA signature. (E) Eight metabolism-related lncRNA signature constructed using iterative Lasso Cox regression analysis. (F) ROC curve to evaluate 1, 3, and 5-year prediction efficiency of the eight metabolism-related lncRNA signature.

Based on the risk score of each OS patient in the TARGET data set, OS patients were divided into high-risk (n = 42) and low-risk (n = 42) groups. Kaplan–Meier survival curve analysis showed that overall survival of OS patients with high-risk scores was significantly shorter than those with low-risk scores (p = 1.113e-04, Figure 6F). Based on the eight metabolism-related lncRNAs prognosis signature, we obtained the risk score distribution, survival status, and a heat map of the eight metabolism-related lncRNAs (Figure 7A). Using univariate Cox analyses revealed that the metastatic and metabolism-related lncRNA prognostic risk scores was significantly associated with overall survival (P < 0.001; Figure 7B). Multivariate Cox analyses showed that the prognostic risk score of the metabolism-related lncRNAs was significantly associated with overall survival (P < 0.001; Figure 7C). Finally, Kaplan–Meier curves showed that patients with different expression levels of the eight metabolism-related lncRNAs had different overall survival (Figure 7D).

FIGURE 7

Figure 7. Evaluation of metabolism-related lncRNA signature. (A) Risk score distribution and survival status of OS patients based on the eight metabolism-related lncRNAs and a heat map of the eight metabolism-related lncRNAs expression pattern. (B) Univariate Cox regression analysis showing the metastatic and the eight metabolism-related lncRNA signature risk score were significantly associated with overall survival. (C) Multivariate Cox regression analysis showing that the eight metabolism-related lncRNA signature risk score was an independent prognostic indicator for overall survival of OS patients. (D) Kaplan–Meier curves showing that the expression of the eight metabolism-related lncRNA had different overall survival for OS patients.

GSEA and GSVA

To examine the potential biological processes involved, we carried out GSEA based on the eight metabolism-related lncRNA signature. The top five biological processes in the high- and low-risk groups are shown in Figures 8A,B. GSVA results show that apical surface and bile acid metabolism, epithelial mesenchymal transition, and P53 pathway were activated in the high-risk OS patients (Figure 8C).

FIGURE 8

Figure 8. Pathway enrichment, ceRNA network, and immunization scores based on metabolism-related lncRNA signature. (A,B) Top five biological processes in the high- and low-risk groups based on GSEA results. (C) GSVA results showing that apical surface and bile acid metabolism, epithelial mesenchymal transition, and P53 pathway were activated in high-risk OS patients. (D) ceRNA network. (E,F) Immunization scores of immune cells and immune-related functions based on the eight metabolism-related lncRNA signature. *p < 0.05, **p < 0.01, ^nsp > 0.05.

Construction of the ceRNA Network and Immunization Score Analysis

To explore the functions of the eight metabolism-related lncRNAs, we constructed a ceRNA network based on the eight metabolism-related lncRNAs (Figure 8D). First, one of the eight metabolism-related lncRNAs were extracted to construct the ceRNA network. We then used three databases (miRTarBase, miRDB, and TargetScan) to search for target mRNAs based on the 11 miRNAs and extracted 44 target mRNAs from differentially expressed genes. Next, we performed immunization scores of immune cells and immune-related functions in the high-risk and low-risk groups using ssGSEA. Results of immunization score analysis show that the immunization scores of B cells, CD8+ T cells, Mast cells, and Th1 cells were significant in the high-risk and low-risk groups (Figure 8E). Furthermore, results of immunization score analysis show that the immunization scores of type II IFN response were significant in the high-risk and low-risk groups (Figure 8F).

Discussion

In this study, eight metabolic genes were found to be significantly correlated with OS based on univariate Cox regression analysis. Then, 6 metabolic genes (PYGM, CKMT2, NAT1, AADAT, GPX7, and CHST13) were selected to construct a prognostic signature based on their performance using the Lasso Cox regression analysis. OS patients with high-risk scores showed shorter survival times compared with those with low-risk scores based on the six metabolic genes signature. Univariate and multivariate Cox analyses suggest that the risk score based on the six metabolic gene prognostic signature was an independent prognostic factor. Stratified correlation analysis shows that the metabolic gene signature accurately predicted survival outcomes in high- and low-risk OS patients. In addition, we constructed a nomogram that was good at predicting the 1, 3, and 5-year survival probabilities using the risk score and other clinicopathological factors, including age, gender, and metastatis. Taken together, these results confirm that the six metabolic genes prognostic signature has good performance in predicting the survival outcomes of OS patients in our study.

In recent years, a mounting body of evidence suggests that reprogramming of metabolism in cancer cells has an important effect on cancer development and progression (Hanahan and Weinberg, 2011). Moreover, an increasing number of studies suggests that protein, lipid, and nucleic acid govern cell growth and are activated in cancer cells via tumorigenic mutations, resulting in cancer development and progression (DeBerardinis and Chandel, 2016; Merino Salvador et al., 2017). However, the association between metabolism and OS progression remains unclear. In our study, we identified the six metabolic genes (PYGM, CKMT2, NAT1, AADAT, GPX7, and CHST13) prognostic signature that showed good performance in predicting survival outcomes of OS patients. However, there are few reports on the role of these genes in OS. PYGM is significantly downregulated in head and neck squamous cell carcinoma (HNSCC) and correlates with worse prognosis of HNSCC (Jin and Yang, 2019). In hepatocellular carcinoma, downregulation of CHST13 regulates the metastasis and chemosensitivity of human hepatocellular carcinoma cells via the mitogen-activated protein kinase (MAPK) pathway (Zhou et al., 2016). Epigenetic inactivation of GPX7 may be an important mechanism of esophageal cancer (Peng et al., 2009; Peppelenbosch et al., 2014). In Barrett’s esophagus, GPX7 suppresses bile salt-induced expression of pro-inflammatory cytokines to inhibit Barrett’s carcinogenesis and is also related to gastroesophageal reflux disease–associated Barrett’s carcinogenesis (Peng et al., 2014a, b). From these data, we know that the six metabolic genes identified in our study play different roles in human cancers. However, research regarding the role of these genes is limited in OS. It is, therefore, worth exploring the functions of these six metabolic genes in OS.

Additionally, 17 metabolism-related lncRNAs were found to be significantly correlated with the overall survival of OS patients in this study. Lasso Cox regression analysis showed that eight metabolism-related lncRNAs (JMJD1C-AS1, FLJ45513, FAM155A-IT1, LINC00837, MANCR, LINC00506, CACNA1C-AS1, and IL10RB-DT) were found to be good candidates for the construction of a prognostic signature. Based on the eight metabolism-related lncRNA signature, the clinical outcome of OS patients with high-risk scores were significantly worse than that for patients with low-risk scores. Like the six metabolic genes signature, the eight metabolism-related lncRNAs signature was also independent of other clinical variables, such as age, gender and metastatis. GSVA results showed that bile acid metabolism, epithelial mesenchymal transition, and P53 pathway were activated in high-risk OS patients. Immunization scores analysis suggests that there was a lower score in the high-risk group than in the low-risk group, and high immunity correlated with good prognosis.

Dysregulation of lncRNAs is known to be involved in tumor growth, metabolism, and metastasis (Lin, 2020). Increasing numbers of studies show that dysregulated lncRNAs have an important effect on glucose, lipid, and cholesterol metabolism by regulating mitochondrial function and oxidative stress (van Solingen et al., 2018; Zeng et al., 2018; Lu et al., 2020). Research on the prognostic potential of metabolism-related lncRNAs is limited in OS. Our results suggest that the eight metabolism-related lncRNAs identified in this study can accurately predict overall survival of OS patients. In the eight metabolism-related lncRNAs, MANCR, LINC00837, LINC00506, and IL10RB-DT are relatively well characterized. Tahmouresi et al. (2020) suggests that MANCR is a potential diagnostic biomarker for breast carcinoma (BC) and is associated with aggressive clinical parameters of BC. The study also shows that MANCR was functionally associated with cell proliferation, viability, and genomic stability and represented a potential therapeutic target for BC (Tracy et al., 2018). It is well known that the P53 pathway is related to the metabolism of cancer cells in multiple cancers (Chen et al., 2017; Goyal et al., 2019), and GSVA results also suggest that the P53 pathway was activated in high-risk OS patients. Taken together, these results also suggest that the eight metabolism-related lncRNA signature has strong ability to predict the prognosis of OS patients.

Despite the identification of six metabolic gene and eight metabolism-related lncRNA prognostic signatures, few reports regarding these signatures have been reported in OS previously. In the future, it will be necessary to explore the molecular biological functions of the six metabolic genes and eight metabolism-related lncRNA signatures in OS tumorigenesis and progression, such as in cell proliferation, cell viability, cell metabolism, cell motility, tumor angiogenesis, and drug resistance. In addition, to further validate these metabolic genes and metabolism-related lncRNA signatures in accurate OS diagnosis and prognosis, more clinical evidence, including prospective large-scale cohorts related to these signatures, is crucial. Last, for the potential application of these signatures in the personalized treatment of OS, we believe that it is important to identify the best biomarker and target from these signatures to improve cellular metabolism and immunotherapy of OS. There are several limitations to our study. First, owing to the limited availability of OS sample size (n = 84) and clinical data, the subgroup analysis based on other clinical characteristics was limited. Second, the OS sample sizes and clinical data of the independent cohort are also hampered. Last, biological functions of the metabolic genes and metabolism-related lncRNAs need to be verified in the future.

In conclusion, we identify six metabolic gene and eight metabolism-related lncRNA prognostic signatures that show good performance in predicting the survival outcomes of OS patients and are independent of other clinical risk factors. Overall, our study suggests that the six metabolic genes and eight metabolism-related lncRNAs are promising prognostic and diagnostic biomarkers for OS therapy and diagnosis.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: TARGET; https://ocg.cancer.gov/programs/target and Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/), GSE21257.

Author Contributions

NW and ZH-H worked on design and conception of this study. GC-Y, TR, SY-Q, and LT-C collected the data, performed the data analysis, and drafted the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by Natural Science Foundation of China (31960175), Natural Science Foundation of Gansu Province (18JR3RA331), and Cuiying Project of Lanzhou University Second Hospital (CY2017-MS06).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

1. ^ https://ocg.cancer.gov/programs/target
2. ^ https://www.gtexportal.org/home/
3. ^ https://www.ncbi.nlm.nih.gov/geo/
4. ^ https://www.r-project.org/
5. ^ https://www.string-db.org/

References

Anwar, M. A., El-Baba, C., Elnaggar, M. H., Elkholy, Y. O., Mottawea, M., Johar, D., et al. (2020). Novel therapeutic strategies for spinal osteosarcomas. Semin. Cancer Biol. 64, 83–92. doi: 10.1016/j.semcancer.2019.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Cairns, R. A., Harris, I. S., and Mak, T. W. (2011). Regulation of cancer cell metabolism. Nat. Rev. Cancer 11, 85–95.

Google Scholar

Chen, J., Brunner, A. D., Cogan, J. Z., Nuñez, J. K., Fields, A. P., Adamson, B., et al. (2020). Pervasive functional translation of noncanonical human open reading frames. Science 367, 1140–1146. doi: 10.1126/science.aay0262

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R., Liu, Y., Zhuang, H., Yang, B., Hei, K., Xiao, M., et al. (2017). Quantitative proteomics reveals that long non-coding RNA MALAT1 interacts with DBC1 to regulate p53 acetylation. Nucleic Acids Res. 45, 9947–9959. doi: 10.1093/nar/gkx600

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

Dayton, T. L., Jacks, T., and Vander Heiden, M. G. (2016). PKM2, cancer metabolism, and the road ahead. Embo Rep. 17, 1721–1730.

Google Scholar

DeBerardinis, R. J., and Chandel, N. S. (2016). Fundamentals of cancer metabolism. Sci. Adv. 2:e1600200. doi: 10.1126/sciadv.1600200

PubMed Abstract | CrossRef Full Text | Google Scholar

Denko, N. C. (2008). Hypoxia, HIF1 and glucose metabolism in the solid tumour. Nat. Rev. Cancer 8, 705–713. doi: 10.1038/nrc2468

PubMed Abstract | CrossRef Full Text | Google Scholar

Dinger, M. E., Gascoigne, D. K., and Mattick, J. S. (2011). The evolution of RNAs with multiple functions. Biochimie 93, 2013–2018. doi: 10.1016/j.biochi.2011.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, C., Tang, Y., Wang, J., Xiong, F., Guo, C., Wang, Y., et al. (2017). Role of long non-coding RNAs in glucose metabolism in cancer. Mol. Cancer 16:130.

Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Goyal, N., Tiwary, S., Kesharwani, D., and Datta, M. (2019). Long non-coding RNA H19 inhibition promotes hyperglycemia in mice by upregulating hepatic FoxO1 levels and promoting gluconeogenesis. J. Mol. Med. 97, 115–126. doi: 10.1007/s00109-018-1718-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of cancer: the next generation. Cell 144, 646–674. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Isakoff, M. S., Bielack, S. S., Meltzer, P., and Gorlick, R. (2015). Osteosarcoma: current treatment and a collaborative pathway to success. J. Clin. Oncol. Offic. J. Am. Soc. Clin. Oncol. 33, 3029–3035. doi: 10.1200/jco.2014.59.4895

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., and Yang, Y. (2019). Bioinformatics-based discovery of PYGM and TNNC2 as potential biomarkers of head and neck squamous cell carcinoma. Biosci. Rep. 39:BSR20191612.

Google Scholar

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361.

Google Scholar

Khan, A., Valli, E., Lam, H., Scott, D. A., Murray, J., Hanssen, K. M., et al. (2020). Targeting metabolic activity in high-risk neuroblastoma through Monocarboxylate Transporter 1 (MCT1) inhibition. Oncogene 39, 3555–3570. doi: 10.1038/s41388-020-1235-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreuzaler, P., Panina, Y., Segal, J., and Yuneva, M. (2020). Adapt and conquer: metabolic flexibility in cancer growth, invasion and evasion. Mol. Metab. 33, 83–101. doi: 10.1016/j.molmet.2019.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, M. M., and Goyal, R. (2017). LncRNA as a therapeutic target for angiogenesis. Curr. Top. Med. Chem. 17, 1750–1757. doi: 10.2174/1568026617666161116144744

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, N., Zhan, X., and Zhan, X. (2018). The lncRNA SNHG3 regulates energy metabolism of ovarian cancer by an analysis of mitochondrial proteomes. Gynecol. Oncol. 150, 343–354. doi: 10.1016/j.ygyno.2018.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J.-H. (2020). Crosstalk of lncRNA and cellular metabolism and their regulatory mechanism in cancer. Int. J. Mol. Sci. 21:2947. doi: 10.3390/ijms21082947

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Q., Guo, P., Liu, A., Ares, I., Martínez-Larrañaga, M.-R., Wang, X., et al. (2020). The role of long noncoding RNA in lipid, cholesterol, and glucose metabolism and treatment of obesity syndrome. Med. Res. Rev. 1–24. doi: 10.1002/med.21775

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, W., Cao, F., Wang, S., Sheng, X., and Ma, J. (2019). LncRNAs: the regulator of glucose and lipid metabolism in tumor cells. Front. Oncol. 9:1099. doi: 10.3389/fonc.2019.01099

PubMed Abstract | CrossRef Full Text | Google Scholar

Luengo, A., GuiM, D. Y., and Vander Heiden, G. (2017). Targeting metabolism for cancer therapy. Cell Chem. Biol. 24, 1161–1180.

Google Scholar

Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159.

Google Scholar

Merino Salvador, M., Gómez, de Cedrón, M., Moreno Rubio, J., Falagán Martínez, S., Sánchez Martínez, R., et al. (2017). Lipid metabolism and lung cancer. Crit. Rev. Oncol. Hematol. 112, 31–40.

Google Scholar

Mirabello, L., TroisiS, R. J., and Savage, S. A. (2009). Savage, International osteosarcoma incidence patterns in children and adolescents, middle ages and elderly persons. Int. J. Cancer 125, 229–234. doi: 10.1002/ijc.24320

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirabello, L., Zhu, B., Koster, R., Karlins, E., Dean, M., Yeager, M., et al. (2020). Frequency of pathogenic germline variants in cancer-susceptibility genes in patients with osteosarcoma. JAMA Oncol. 6, 724–734.

Google Scholar

Peng, D. F., Hu, T. L., Soutto, M., Belkhiri, A., and El-Rifai, W. (2014a). Glutathione Peroxidase 7 suppresses bile salt-induced expression of pro-inflammatory cytokines in Barrett’s Carcinogenesis. J. Cancer 5, 510–517. doi: 10.7150/jca.9215

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, D. F., Hu, T. L., Soutto, M., Belkhiri, A., and El-Rifai, W. (2014b). Loss of glutathione peroxidase 7 promotes TNF-α-induced NF-κB activation in Barrett’s carcinogenesis. Carcinogenesis 35, 1620–1628. doi: 10.1093/carcin/bgu083

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, D. F., Razvi, M., Chen, H., Washington, K., Roessner, A., Schneider-Stock, R., et al. (2009). DNA hypermethylation regulates the expression of members of the Mu-class glutathione S-transferases and glutathione peroxidases in Barrett’s adenocarcinoma. Gut 58, 5–15. doi: 10.1136/gut.2007.146290

PubMed Abstract | CrossRef Full Text | Google Scholar

Peppelenbosch, M. P., Spaander, M. C., and Bruno, M. J. (2014). Glutathione peroxidase 7 prevents cancer in the oesophagus. Gut 63, 537–538. doi: 10.1136/gutjnl-2013-304906

PubMed Abstract | CrossRef Full Text | Google Scholar

Pingping, B., Yuhong, Z., Weiqi, L., Chunxiao, W., Chunfang, W., Yuanjue, S., et al. (2019). Incidence and mortality of sarcomas in Shanghai, China, During 2002-2014. Front. Oncol. 9:662. doi: 10.3389/fonc.2019.00662

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Strauss, S. J., Ng, T., Mendoza-Naranjo, A., Whelan, J., and Sorensen, P. H. (2010). Understanding micrometastatic disease and Anoikis resistance in ewing family of tumors and osteosarcoma. oncologist 15, 627–635. doi: 10.1634/theoncologist.2010-0093

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L. Y., Li, X. J., Sun, Y. M., Huang, W., Fang, K., Han, C., et al. (2018). LncRNA ANRIL regulates AML development through modulating the glucose metabolism pathway of AdipoR1/AMPK/SIRT1. Mol. Cancer 17:127.

Google Scholar

Tahmouresi, F., Razmara, E., Pakravan, K., Mossahebi-Mohammadi, M., Rouhollah, F., Montazeri, M., et al. (2020). Upregulation of the long noncoding RNAs DSCAM-AS1 and MANCR is a potential diagnostic marker for breast carcinoma. Biotechnol. Appl. Biochem. [Epub ahead of print]. doi: 10.1002/bab.2048.

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, J., Yan, T., Bao, Y., Shen, C., Yu, C., Zhu, X., et al. (2019). LncRNA GLCC1 promotes colorectal carcinogenesis and glucose metabolism by stabilizing c-Myc. Nat. Commun. 10:3499.

Google Scholar

Tracy, K. M., Tye, C. E., Ghule, P. N., Malaby, H. L. H., Stumpff, J., Stein, J. L., et al. (2018). Mitotically-associated lncRNA (MANCR) affects genomic stability and cell division in aggressive breast cancer. Mol. Cancer Res. 16, 587–598. doi: 10.1158/1541-7786.mcr-17-0548

PubMed Abstract | CrossRef Full Text | Google Scholar

van Solingen, C., Scacalossi, K. R., and Moore, K. J. (2018). Long noncoding RNAs in lipid metabolism. Curr. Opin. Lipidol. 29, 224–232.

Google Scholar

Vander Heiden, M. G., and DeBerardinis, R. J. (2017). Understanding the intersections between metabolism and cancer biology. Cell 168, 657–669. doi: 10.1016/j.cell.2016.12.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Vergara, D., Verri, T., Damato, M., Trerotola, M., Simeone, P., Franck, J., et al. (2020). A hidden human proteome signature characterizes the epithelial mesenchymal transition program. Curr. Pharmaceut. Design 26, 372–375. doi: 10.2174/1381612826666200129091610

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, Y., Ren, K., Zhu, X., Zheng, Z., and Yi, G. (2018). Long noncoding RNAs: advances in lipid metabolism. Adv. Clin. Chem. 87, 1–36. doi: 10.4137/atv.s29816

CrossRef Full Text | Google Scholar

Zhou, H., Li, Y., Song, X., Zhao, Y., Cheng, L., Zhao, L., et al. (2016). CHST11/13 regulate the metastasis and chemosensitivity of human hepatocellular carcinoma cells via mitogen-activated protein kinase pathway. Digest. Dis. Sci. 61, 1972–1985. doi: 10.1007/s10620-016-4114-5

PubMed Abstract | CrossRef Full Text | Google Scholar