A four-lncRNA risk signature for prognostic prediction of osteosarcoma

Abstract

Aim: Osteosarcoma is the most common primary malignant tumor of bone. However, our understanding of the prognostic indicators and the genetic mechanisms of the disease progression are still incomplete. The aim of this study was to identify a long noncoding RNA (lncRNA) risk signature for osteosarcoma survival prediction.

Methods: RNA sequencing data and relevant clinical information of osteosarcoma patients were downloaded from the database of Therapeutically Applicable Research to Generate Effective Treatments (TARGET). We analyzed the differentially expressed lncRNAs between deceased and living patients by univariate and multivariate Cox regression analysis to identify a risk signature. We calculated a prognostic risk score for each sample according to this prognosis signature, and divided patients into high-risk and low-risk groups according to the median value of the risk score (0.975). Kaplan–Meier analysis and receiver operating characteristic (ROC) curve statistics were used to evaluate the performance of the signature. Next, we analyzed the signature’s potential function through Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene-set enrichment analysis (GSEA). Lastly, qRT-PCR was used to validate the expression levels of the four lncRNAs in clinical samples.

Results: Twenty-six differentially expressed lncRNAs were identified between deceased and living patients. Four of these lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) were identified as independent prognostic factors, and a risk signature of these four lncRNAs for osteosarcoma survival prediction was constructed. Kaplan–Meier analysis showed that the five-year survival time in high-risk and low-risk groups was 33.1% and 82.5%, and the area under the curve (AUC) of the ROC was 0.784, which demonstrated that the prognostic signature was reliable and had the potential to predict the survival of patients with osteosarcoma. The expression level of the four lncRNAs in osteosarcoma tissues and cells was determined by qRT-PCR. Functional enrichment analysis suggested that the signature might be related to osteosarcoma through regulation of the MAPK signaling pathway, the PI3K-Akt signaling pathway, and the extracellular matrix and also provided new insights into the study of osteosarcoma, including the role of papillomavirus infection, olfactory receptor activity, and olfactory transduction in osteosarcoma.

Conclusion: We constructed a novel lncRNA risk signature that served as an independent biomarker for predicting the prognosis of osteosarcoma patients.

Introduction

As the most common malignant tumor of bone, osteosarcoma has a poor prognosis and a high incidence rate in children and adolescents (Kansara et al., 2014). In most cases, osteosarcomas can be found in the long bones of limbs, near the metaphyseal growth plate. Femur, tibia, and humerus are common sites, while skull, jaw, and pelvis are much less common (Ritter and Bielack, 2010; Simpson and Brown, 2018). Current therapies incorporate surgical resection with further combinational chemotherapy including doxorubicin, methotrexate, cisplatin, or ifosfamide, which not only shrinks the tumor but also eradicates micro-metastases and cures about 70% of patients (Lamoureux et al., 2007). However, for the past 30 years, the five-year survival rate of 20% for patients with metastatic or relapsed osteosarcoma has remained virtually unchanged (Collins et al., 2013). Despite the fact that it was known that tumor size, localized sites, present stage, and reaction to chemotherapy were of great significance in predicting overall survival, our perception and understanding of the prognostic indicators as well as the genetic mechanisms of disease progression were still incomplete and needed further exploration. It has been reported that genomic complexity and instability played a crucial role in the prognosis of osteosarcoma. Therefore, it was necessary to optimize the early detection, early diagnosis, and prognosis prediction of osteosarcoma from the perspective of molecular genetics.

Over the last decade, non-coding RNA has become a hot topic in research on diagnostic biomarkers and mechanism of many kinds of tumors. Long non-coding RNA (lncRNA) is a type of non-coding RNA with a length >200 nucleotides (Quinn and Chang, 2016), which has attracted worldwide attention because of its role in numerous diseases. Recent next-generation and high-throughput sequencing techniques have led to major breakthroughs in lncRNA identification and characterization. Accumulating evidence revealed that various lncRNAs were dysregulated in cancers, and were related to cancer recurrence, metastasis, and poor prognosis through epigenetic, transcriptional, and posttranscriptional mechanisms in various types of biological processes, including tumor initiation, growth, and metastasis (Kopp and Mendell, 2018; Sanchez Calle et al., 2018). Some studies have proved that lncRNAs were related to osteosarcoma progression (Wang et al., 2018; Chi et al., 2019; Zhu et al., 2019; Shen and Cheng, 2020). Nevertheless, whether lncRNAs could be used as prognostic markers of osteosarcoma was still unknown. Thus, further study of the clinical prognostic value of lncRNAs was definitely needed to offer osteosarcoma patients better therapeutic options. Existing and current studies were marshalled to categorize lncRNAs as robust predictors for osteosarcoma patients’ prognosis and to establish a prognostic risk signature based on comprehensive analysis of RNA sequencing data.

Our goal was to establish a prognostic model for osteosarcoma survival prediction, making use of RNA expression data as the foundation of a broad spectrum of genetic and clinical factors. A large cohort of osteosarcoma patients from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) project was used to identify potential lncRNAs related to the overall survival of patients with osteosarcoma. The lncRNA expression profiles and clinical characteristics were analyzed with univariable Cox regression and multivariable Cox regression, and four lncRNAs related to prognosis of osteosarcoma were identified and selected for the study. In addition, we constructed a four-lncRNA signature that could act as an independent prognostic factor for use in predicting overall survival in patients with osteosarcoma. Altogether, the current study presents a novel and effective method for the prognosis of patients with osteosarcoma.

Methods

Data collection and preprocessing

RNA sequencing (RNA-seq) data on mRNAs and lncRNAs and relevant clinical information of osteosarcoma patients from the TARGET project were downloaded from the website (https://ocg.cancer.gov/programs/target) in June 2021. RNA-seq data files were merged into a matrix file using the merge script of the Perl language (Supplementary Material S1), and the merged RNA matrix is shown in Supplementary Material S2. The gene names were converted from the Ensembl ID to the matrix of the gene symbol through the Ensembl database (http://asia.ensembl.org/index.html). A total of 101 osteosarcoma TARGET RNA-seq data samples were downloaded from the database. The RNA-Seq data were merged with clinical information about the patients, such as vital status, survival time, and gender. There remained 93 samples after excluding those with duplicate data and samples without vital status and survival time. In order to identify survival-related lncRNAs, the 93 samples were divided into two groups according to the vital status. The R package, edgeR (Robinson et al., 2010), was used to identify genes that were differentially expressed between 55 samples from living patients and 38 samples from deceased patients, with a false discovery rate (FDR) < 0.05 and log₂FC (FC) >1 as the threshold.

Identification of a prognostic lncRNA signature

Univariate Cox proportional hazard regression analysis was performed using the R survival package (Lee, 2019) to screen out potential lncRNAs related to the prognosis of osteosarcoma patients. Only those lncRNAs with p < 0.05 were considered candidate prognostic factors of osteosarcoma. Multivariate Cox proportional hazards regression was performed using the R package, “survival” (Lee, 2019), to identify significant prognostic lncRNAs that affect the survival of osteosarcoma patients. In the multivariate Cox model, sex, age, and tumor metastasis status could be covariates but were not considered in our study. We mainly analyzed the 11 lncRNAs identified in the univariate Cox. Lastly, we chose four lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) according to p<0.05 and constructed a prognosis signature of osteosarcoma patients, based on the expression level of the four lncRNAs as a variable. The risk score of the prognosis signature was calculated as follows:Here, lncRNAi is the identifier of the selected lncRNAs. A high-risk score indicates poor survival for osteosarcoma patients.

Survival analysis and ROC evaluation

One essential characteristic of a reasonable prognostic signature is that it should have no correlation with the currently used clinicopathological prognostic factors. In order to evaluate and measure the independence as well as applicability of such a signature, we calculated the prognostic risk score for each sample based on the signature and categorized patients into high-risk and low-risk groups, taking the median value of the risk score into consideration. Kaplan–Meier analysis was also performed to create overall survival curves. Meanwhile, log-rank tests were utilized to evaluate surviving discrepancy of both high-risk and low-risk patient groups. By calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curve (Liu et al., 2020; Yu et al., 2020), the performance of lncRNA prognostic signature was evaluated. In this section, Kaplan–Meier analysis was performed using the R package, “survival” (Lorent et al., 2014), and the ROC curve calculation was performed using the R package, “survivalROC”(Lorent et al., 2014).

Functional enrichment analysis

On the basis of univariate and multivariate analyses, there were four prognostic lncRNAs that could be correlated with the prognosis of osteosarcoma patients. LncRNA-correlated protein-coding genes (PCGs) underlie the foundation of the analysis of functional enrichment. The mRNA sequence data of osteosarcoma samples were obtained from the TARGET data matrix. Pearson correlation coefficients among the expression profiles of all four prognostic lncRNAs and PCGs were calculated to distinguish the co-expression relationships between lncRNAs and PCGs. PCGs with the Pearson correlation coefficient >0.40 and p<0.001 were regarded as lncRNA-related PCGs. To determine the biological processes and pathways through which the lncRNAs may play their role, we used functional analysis of Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2021) via the database for both annotation and visualization, as well as the integrated discovery (http://metascape.org/) internet tool. Significant functional categories were defined as those in which p<0.05. GO and KEGG pathway analysis was performed using the R package, “clusterProfiler” (Yu et al., 2012).

Protein–protein interaction network construction

In order to explore the interactions between proteins and the most likely target proteins, protein–protein interaction (PPI) analysis of lncRNA-correlated PCGs was performed using the STRING protein database 11.0 (http://string-db.org/). We set the cut-off criterion for the interaction score as >0.7.

Gene set enrichment analysis

To explore the important functional phenotypes between the high- and low-risk groups based on the four-lncRNAs signature, we performed gene set enrichment analysis (GSEA). Using c2.all.v7.0.symbols.gmt and c5.all.v7.0.symbols.gmt as reference gene sets, enrichment analysis was performed using GSEA software (ver 4.0.3). A nominal p<0.05 and a false discovery rate <0.05 were considered statistically significant.

Cell culture

The human skeletal muscle cell line (HSMC) and osteosarcoma cell line (OS) were obtained from the American Type Culture Collection (ATCC, VA, USA). Cells were incubated at 37°C with 5 % CO₂ in Dulbecco’s modified Eagle’s medium (DMEM) (HyClone, UT, USA) supplemented with 10% fetal bovine serum (Gibco, CA, USA) and 1 % penicillin–streptomycin (HyClone, UT, USA).

Tissue specimens

Tumor samples were obtained from 15 osteosarcoma patients who underwent tumorectomy in Shandong Provincial Hospital from January 2020 to July 2022. Inclusion criteria were as follows: all patients had a definite pathological diagnosis of osteosarcoma and underwent tumor biopsy or resection. Exclusion criteria included history of other types of cancer and history of evidence of severe acute or chronic infection or infectious illness. Patients treated previously for osteosarcoma and patients not willing to give consent were also excluded from the study. Skeletal tissues from 10 trauma patients without any tumor were collected as controls. All tissue samples were stored in liquid nitrogen and briefly stored in a −80°C refrigerator when being processed. The clinical characteristics of the 15 osteosarcoma patients are shown in Supplementary File S3.

RNA extraction and qRT-PCR analysis

RNA extraction and qRT-PCR analysis were carried out as described in our previous study (Zhang et al., 2021a). Total RNA from cultured cells was extracted using SparkZol reagent (SparkJade, Shandong, China), and cDNA was synthesized with the Evo M-MLV RT Premix for qPCR (AG11706, Accurate Biotechnology, Hunan, China). The qRT-PCR analysis was carried out using the SYBR Green Premix Pro Taq HS qPCR Kit (AG11701, Accurate Biotechnology, Hunan, China) as directed by the manufacturer. Primers were synthesized by BioSune Co., Ltd. (Shanghai, China) and the sequences are shown in Table 1. The qRT-PCR was performed to evaluate the expression levels of the four lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) using a LightCycler480 system (Roche Diagnostics, Switzerland). The lncRNAs expression was normalized to ACTB, and the relative expression level was determined by the 2^-∆∆CT method.

The study design and workflow of the present research is shown in Figure 1A. This study was approved by the Ethics Committee of Shandong Provincial Hospital (the ethics approval number is SZZJJ:NO.2021-419) and was carried out in accordance with relevant guidelines and regulations.

FIGURE 1

Results

Sample characteristics

The RNA-Seq data were obtained from the TARGET-OS database and merged with patient clinical information, such as vital status, survival time, and gender. There were a total of 93 samples including 55 from living patients and 38 from deceased patients, after excluding duplicate data and samples without vital status and survival time. The clinical characteristics are presented in Table 2. The specific clinical characteristics of each sample are shown in Supplementary File S4.

Identification of survival-related, differentially expressed lncRNAs

Differential gene expression between deceased and living patients was compared with the edgeR package using R studio. In total, 26 differentially expressed lncRNAs were identified, of which 18 were upregulated and 8 were downregulated in the deceased group, with thresholds for FDR <0.05 and log₂FC (FC) >1 (Supplementary File S5). Volcano plot analysis was used to show variations in lncRNA expression between samples form living and deceased patients (Figure 1B).

Construction of an lncRNA-based multivariable Cox model

To identify prognostic lncRNAs, expression data of 26 differentially expressed lncRNAs and patients’ survival data were subjected to univariate Cox proportional hazards regression analysis. Eleven lncRNAs were found to be significantly associated with overall survival (p<0.05) (Supplementary File S6). Multivariate Cox proportional hazards regression analysis was next used to select the optimal lncRNAs as significant independent prognostic factors. Lastly, we identified four lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) according to p<0.05 and constructed a prognosis signature of osteosarcoma, based on the expression levels of these four lncRNAs. Details about the four lncRNAs are shown in Table 3. The expression of CTB-4E7.1 and RP11-24N18.1 was upregulated and the expression of RP11-553A10.1 and PVRL3-AS1 was down-regulated in deceased compared with living patients. The results of the four prognostic lncRNAs in the univariable and multivariable Cox regression models are shown in Table 4. A four-lncRNA survival-related risk signature for osteosarcoma survival prediction was constructed. This signature was developed as a linear combination of the expression levels of the four lncRNAs weighted by their relative regression coefficients by multivariate Cox regression as follows:

Based on our results, the formula was:

Osteosarcoma survival prediction using the four-lncRNA risk signature

Kaplan–Meier analysis of the four lncRNAs was performed and the expression of CTB-4E7.1 showed no influence on patient survival (p>0.05) (Figure 2A). The expression of RP11-553A10.1 and PVRL3-AS1 were positively correlated with survival of osteosarcoma patients, while RP11-24N18.1 showed the opposite effect (Figures 2B–D). Next, we calculated the prognostic risk score for each sample according to this prognosis signature, and the median value of the risk score was 0.975 (Supplementary File S7). Then we divided patients into high-risk and low-risk groups according to the median value of the risk score (Figure 3A). The survival status and survival time of the osteosarcoma patients is shown in the dot plot (Figure 3B), and indicates that as the risk score increased, the risk of death also increased. The specific survival times associated with samples from the deceased and the length of follow-up for the living patient samples are shown in Supplementary File S4. Kaplan–Meier analysis was performed to compare the survival rates between the high-risk and low-risk groups. The overall survival of the two groups differed significantly, suggesting that the survival rate of high-risk patients was lower than that in low-risk patients at any point (p = 3.566e-07). The overall five-year survival in the high-risk and low-risk groups was 33.1% and 82.5%, respectively (Figure 3C). In addition, we constructed a ROC curve to assess the predictive performance of the four-lncRNA prognostic risk model (Figure 3D). The AUC value was 0.784, proving that the prognostic signature was convincing and had the potential to evaluate the survival of osteosarcoma patients.

FIGURE 2

FIGURE 3

Functional enrichment analysis

A total of 746 PCGs related to the four lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) with Pearson correlation coefficients >0.40 and p<0.001 were identified. To explore the functional implications of the four-lncRNA signature, we performed GO and KEGG functional enrichment analysis to reveal the specific functional categories of lncRNA-correlated PCGs. Figure 4A highlights the most significantly enriched GO terms of the PCGs. The top three were extracellular matrix structural constituents, collagen binding, and integrin binding. After KEGG analysis, we found that the MAPK signaling pathway, the PI3K-Akt signaling pathway, and human papillomavirus infection were the most three enriched pathways (Figure 4B).

FIGURE 4

Protein–protein interaction network construction

To deduce the connections among the lncRNA-correlated PCGs, a PPI network was constructed using the STRING protein database, including 353 nodes and 739 edges (Supplementary File S8). We calculated the number of genes connected to each node (gene) in the network and the average was 4.19. The top 10 genes with the most connections are shown in Figure 5B and include EGFR, GNG12, BMP4, and GNG4, which might play important roles in the PPI network. EGFR had the most connections to other genes, and Figure 5A shows the EGFR network.

FIGURE 5

Gene set enrichment analysis

To further understand the potential biological processes involved and the underlying signaling pathways, GSEA was carried out to determine the associated signaling pathways between the two risk groups. Through KEGG pathway analysis, the top five pathways enriched in the high-risk group were olfactory receptor activity, sensory perception of smell, sensory perception of chemical stimulus, detection of chemical stimuli, and detection of stimuli involved in sensory perception. In contrast, in the low-risk group five other pathways, including DNA templated transcription termination, nucleosome binding, bone cell development, termination of RNA polymerase II transcription, and basement membrane organization were significantly enriched (Table 5 and Figure 6A). GO gene sets analysis showed that the high-risk groups were positively associated with olfactory transduction, ribosomes, retinol metabolism, drug metabolism, cytochrome P450 drug metabolism, and metabolism of xenobiotics by cytochrome P450, and negatively related to lysosomes, focal adhesion, Fc gamma R-mediated phagocytosis, ubiquitin-mediated proteolysis, and glycosaminoglycan biosynthesis heparan sulfate (Table 6 and Figure 6B).

FIGURE 6

Validation of lncRNA expression

The expression of the four selected lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) was validated by RT-qPCR in osteosarcoma tumors. Skeletal muscle from trauma patients served as the control tissue. The results showed that CTB-4E7.1 and RP11-24N18.1 were upregulated in osteosarcoma tissues, while RP11-553A10.1 and PVRL3-AS1 were down-regulated in osteosarcoma compared with controls. (Figure 7). In addition, we detected the expression of the four lncRNAs in osteosarcoma cells compared with human skeletal muscle cell (HSMC) by RT-qPCR. The results were consistent with the clinical samples (Figure 8).

FIGURE 7

FIGURE 8

Discussion

Osteosarcoma is the most common primary malignant tumor of bone, with a poor prognosis and a high incidence rate in children and adolescents. Research has shown that a series of clinical characteristics, such as tumor size, localized sites, stage at presentation, and response to chemotherapy are significant for predicting overall survival of osteosarcoma patients; however, the lack of specific and sensitive biomarkers of prognosis is still an urgent problem in the management of osteosarcoma patients. In recent years, with the development of microarray and high-throughput sequencing technology, the expression patterns of mRNAs and lncRNAs have been widely studied. Many mRNAs have been identified as prognostic markers for osteosarcoma patients, and several prognostic signatures based on mRNAs have been constructed (Zhang et al., 2020a; Rothzerg et al., 2021; Wu et al., 2021; Yin et al., 2021). Although several previous studies have identified a series of lncRNAs that play an important role in osteosarcoma, proof of the efficacy of lncRNAs as prognostic markers of osteosarcoma is still limited. SiYuan Yu et al. identified a five-lncRNA, metastasis-related risk signature for osteosarcoma survival prediction, and the AUC for predicting 5-year survival was 0.745 (Yu et al., 2021). Gao et al. (2020) constructed a multivariable Cox regression model based on seven lncRNAs for prognosis prediction, among which AC011442.1 was assumed to act as an oncogenic driver in osteosarcoma. Our study identified a four-lncRNAs risk signature for osteosarcoma survival prediction.

In this study, we first found 26 survival-related differentially expressed lncRNAs between deceased and living patients based on RNA sequencing data and clinical information obtained from the TARGET database. These 26 differentially expressed lncRNAs and the patient survival data were next subjected to univariate and multivariate Cox proportional hazards regression, and four lncRNAs (CTB-4E7.1, RP11-553A10.1, RP11-24N18.1, and PVRL3-AS1) were identified. A four-lncRNA risk signature for osteosarcoma survival prediction was constructed and the subsequent log-rank tests and Kaplan–Meier analyses further confirmed the predictive value of the risk score. A high-risk score was related to poor prognosis. Based on the risk scores, we divided the patients into a high-risk group (>0.975) and a low-risk group (<0.975), and then calculated the five-year survival rate by Kaplan–Meier analysis. In comparison with other signatures, our AUC value of 0.784 was larger than the 0.745 AUC value of SiYuan Yu’s signature (Yu et al., 2021). This confirmed that our prognostic signature was convincing and had the potential to evaluate the survival of patients with osteosarcoma.

To explore the functional implication of the four-lncRNAs signature, we performed GO and KEGG functional enrichment analysis to reveal particular functional categories of lncRNA-correlated PCGs. The most highly enriched GO terms were extracellular matrix structural constituents, collagen binding, and integrin binding. Collagen and integrin are important constituents of the extracellular matrix (Zeltz and Gullberg, 2016). The extracellular matrix in the tumor microenvironment is a complex network of structural and instructional molecules tightly connected to the tumor cells, and the disorganization of extracellular matrix structure promotes tumor cellular transformation and metastasis (Jain et al., 2013; Schaefer and Reinhardt, 2016). The extracellular matrix microenvironment plays an important role in regulating tumor cell behavior, and many studies have proved that the extracellular matrix undergoes numerous changes in composition and organization in cancer development and progression (Crotti et al., 2017). Based on this and our GO enrichment results, we speculated that the four-lncRNAs signature might reflect the risk of osteosarcoma progression by regulating the extracellular matrix and tumor microenvironment.

After KEGG analysis, the MAPK signaling pathway, PI3K-Akt signaling pathway, and human papillomavirus infection were the most highly enriched pathways. As two well-studied pathways, the MAPK and PI3K-Akt signaling pathways have been proved to play important roles in osteosarcoma, by regulating cell proliferation, apoptosis, migration, and metastasis (Cheng et al., 2017; Han et al., 2019; Zhang et al., 2020b; Zhang et al., 2020c). Thus, our four-lncRNAs signature might influence osteosarcoma prognosis through effects on the MAPK and PI3K-Akt signaling pathways, which supports the credibility of the signature from its functional linkage. Interestingly, the third most enriched pathway had to do with human papillomavirus infection. Papillomavirus infection is a critical factor in the development of cervical cancer in women. However, a relationship between papillomavirus infection and osteosarcoma is still not clear, and this finding points the way to a new direction in the study of osteosarcoma etiology.

We created a PPI network to identify the potential biological interactions among co-expressed lncRNA-associated PCGs (Quo et al., 2012; Li et al., 2019; Chen et al., 2020). The results showed that EGFR, GNG12, BMP4, and GNG4 were key players in the PPI network. Several studies have shown a link between EGFR and BMP4 in osteosarcoma (Li et al., 2020a; Nakajima et al., 2021). However, a role for GNG12 and GNG4 in the mechanism of osteosarcoma progression has not been reported, which also provides a new direction for osteosarcoma research.

In further investigating the lncRNA signature, we carried out GSEA to determine the associated signaling pathways between high-risk and low-risk groups. The most enriched pathways in the high-risk group were in olfactory receptor activity and the most enriched GO term was olfactory transduction in the high-risk group. What is the relationship between olfactory function and osteosarcoma? This result is intriguing. Olfactory receptor activity and olfactory transduction have been found to be involved in the regulation of cell-cell recognition, migration, proliferation, the apoptotic cycle, exocytosis, and pathfinding processes. Additionally, accumulating evidence showed that olfactory receptors were highly expressed in different cancer tissues and they have the potential to serve as diagnostic and therapeutic targets (Massberg and Hatt, 2018). However, the role of olfactory receptor activity and olfactory transduction in osteosarcoma has not been studied and needs further investigation.

Lastly, we verified the expression of the four lncRNAs in clinical samples and osteosarcoma cell lines, since these four lncRNAs have never been reported, and they are only listed in the database. We validated their authenticity by qRT-PCR on clinical samples and cell lines.

The potential biological mechanisms of the four identified biomarkers for prognosis prediction include a number of possibilities. Based on the GO enrichment results, the four-lncRNAs signature might influence osteosarcoma by regulating the extracellular matrix and tumor microenvironment (Cohen et al., 2016; Zhou et al., 2016; Yan et al., 2018; Zheng et al., 2019), while the KEGG data suggested that MAPK and PI3K-Akt signaling may be important through their effects on cell proliferation, apoptosis, migration, and metastasis (Xiang et al., 2014; Li et al., 2015; Zi et al., 2015; Li et al., 2020b; Yang et al., 2020; Naz et al., 2021). Olfactory receptor activity and the olfactory transduction pathway may participate in the progression of osteosarcoma by RNA and DNA modifications (Wang et al., 2015; Jin et al., 2020; Wang et al., 2020). In the future, we will explore the mechanism of the four lncRNAs through cell and animal experiments, and the causal effects of these our lncRNAs on the prognosis will be further investigated through examination of additional multi-omics data sets (Hou et al., 2020; Wu et al., 2020; Wang et al., 2021; Zhang et al., 2021b).

The signature presented in this study has several advantages over previous prognostic models. First, the four lncRNAs were selected from survival-related lncRNAs, which were differently expressed in living vs. deceased patients. The other studies mainly selected metastasis-related or immune-related genes. The survival-related lncRNAs were better able to predict osteosarcoma prognosis. Compared with other signatures, our AUC value (0.784) was larger. These facts demonstrated that our prognostic signature was most convincing and had the potential to evaluate the survival of patients with osteosarcoma. In addition, the functional analysis of the signature indicated that the lncRNAs might influence osteosarcoma by regulating the extracellular matrix, the tumor microenvironment, or the MAPK and PI3K-Akt signaling pathways.

Some limitations of this study should be pointed out, however. First, the dysregulated survival-related lncRNAs in osteosarcoma were analyzed using previously released data sets and validated by RT-qPCR in osteosarcoma tissues and cells, but functional verification was still lacking. Future experimental validation is needed to verify our findings. Second, all the data in our research were collected from the TARGET database, and an external validation cohort for the prognostic model is needed. The TARGET database is dedicated to comprehensive genetic information on childhood tumors, like osteosarcoma. The data in the TARGET database include mRNA, miRNA, lncRNA, and detailed clinical information, and in our study, we used lncRNA-seq data, mRNA-seq data, and clinical information. For this reason, we only used TARGET data and did not include other public data. In the future, we will explore the public data related to osteosarcoma from the GEO and TCGA databases. However, there may be duplications in these databases, so careful screening will be necessary before using data from them. We hope that, when more new datasets become available in the near future, we can further confirm the accuracy and therapeutic efficacy of our signature. In addition, with the rapid development of artificial intelligence (AI), it is being used more and more in clinical medicine and bioinformatics. AI provides new ideas and methods for diagnosis and prognosis prediction of some diseases. In planned studies, we will adapt and apply AI to construct and test new clinical prediction models (Chen et al., 2022; Gao et al., 2022).

Conclusion

We established a novel risk signature and demonstrated its credibility and rationale for osteosarcoma survival prediction. The signature not only serves as a new biomarker for predicting osteosarcoma prognosis but also provides new insights into the molecular mechanisms of osteosarcoma progression. We will further confirm the signature’s performance using external validation cohorts in the near future. It is hoped that these studies will be applied to the diagnostic and prognostic assessment of osteosarcoma patients, provide new insights into the pathogenesis of osteosarcoma, and promote the development of more effective biomarkers.

References

• 1

  ChenJ.ZhaoX.CuiL.HeG.WangX.WangF.et al (2020). Genetic regulatory subnetworks and key regulating genes in rat hippocampus perturbed by prenatal malnutrition: Implications for major brain disorders. Aging (Albany NY)12 (9), 8434–8458. 10.18632/aging.103150

  • CrossRef
  • Google Scholar

• 2

  ChenT.ZhengW.HuH.LuoC.ChenJ.YuanC.et al (2022). A corresponding region fusion framework for multi-modal cervical lesion detection. IEEE/ACM Trans. Comput. Biol. Bioinform.2022, 1. 10.1109/TCBB.2022.3178725

  • CrossRef
  • Google Scholar

• 3

  ChengD. D.ZhuB.LiS. J.YuanT.YangQ. C.FanC. Y. (2017). Down-regulation of RPS9 inhibits osteosarcoma cell growth through inactivation of MAPK signaling pathway. J. Cancer8 (14), 2720–2728. 10.7150/jca.19130

  • CrossRef
  • Google Scholar

• 4

  ChiY.WangD.WangJ.YuW.YangJ. (2019). Long non-coding RNA in the pathogenesis of cancers. Cells8 (9), 1015. 10.3390/cells8091015

  • CrossRef
  • Google Scholar

• 5

  CohenO. S.WeickertT. W.HessJ. L.PaishL. M.McCoyS. Y.RothmondD. A.et al (2016). A splicing-regulatory polymorphism in DRD2 disrupts ZRANB2 binding, impairs cognitive functioning and increases risk for schizophrenia in six Han Chinese samples. Mol. Psychiatry21 (7), 975–982. 10.1038/mp.2015.137

  • CrossRef
  • Google Scholar

• 6

  CollinsM.WilhelmM.ConyersR.HerschtalA.WhelanJ.BielackS.et al (2013). Benefits and adverse events in younger versus older patients receiving neoadjuvant chemotherapy for osteosarcoma: Findings from a meta-analysis. J. Clin. Oncol.31 (18), 2303–2312. 10.1200/JCO.2012.43.8598

  • CrossRef
  • Google Scholar

• 7

  CrottiS.PiccoliM.RizzolioF.GiordanoA.NittiD.AgostiniM. (2017). Extracellular matrix and colorectal cancer: How surrounding microenvironment affects cancer cell behavior?J. Cell. Physiol.232 (5), 967–975. 10.1002/jcp.25658

  • CrossRef
  • Google Scholar

• 8

  GaoH.GuoY.ZhangM.YiZ. (2020). Comprehensive characterization of prognostic long noncoding RNAs in osteosarcoma. Biomed. Res. Int.2020, 6725753. 10.1155/2020/6725753

  • CrossRef
  • Google Scholar

• 9

  GaoH.XiaoJ.YinY.LiuT.ShiJ. (2022). A mutually supervised graph attention network for few-shot segmentation: The perspective of fully utilizing limited samples. IEEE Trans. Neural Netw. Learn Syst.2022, 1. 10.1109/TNNLS.2022.3155486

  • CrossRef
  • Google Scholar

• 10

  HanD.WangM.YuZ.YinL.LiuC.WangJ.et al (2019). FGF5 promotes osteosarcoma cells proliferation via activating MAPK signaling pathway. Cancer Manag. Res.11, 6457–6466. 10.2147/CMAR.S200234

  • CrossRef
  • Google Scholar

• 11

  HouL.XuM.YuY.SunX.LiuX.LiuL.et al (2020). Exploring the causal pathway from ischemic stroke to atrial fibrillation: A network mendelian randomization study. Mol. Med.26 (1), 7. 10.1186/s10020-019-0133-y

  • CrossRef
  • Google Scholar

• 12

  JainP.BaranwalS.DongS.StruckhoffA. P.WorthylakeR. A.AlahariS. K. (2013). Integrin-binding protein nischarin interacts with tumor suppressor liver kinase B1 (LKB1) to regulate cell migration of breast epithelial cells. J. Biol. Chem.288 (22), 15495–15509. 10.1074/jbc.M112.418103

  • CrossRef
  • Google Scholar

• 13

  JinG.XuM.ZouM.DuanS. (2020). The processing, gene regulation, biological functions, and clinical relevance of N4-acetylcytidine on RNA: A systematic review. Mol. Ther. Nucleic Acids20, 13–24. 10.1016/j.omtn.2020.01.037

  • CrossRef
  • Google Scholar

• 14

  KanehisaM.FurumichiM.SatoY.Ishiguro-WatanabeM.TanabeM. (2021). Kegg: Integrating viruses and cellular organisms. Nucleic Acids Res.49 (D1), D545–D551. 10.1093/nar/gkaa970

  • CrossRef
  • Google Scholar

• 15

  KansaraM.TengM. W.SmythM. J.ThomasD. M. (2014). Translational biology of osteosarcoma. Nat. Rev. Cancer14 (11), 722–735. 10.1038/nrc3838

  • CrossRef
  • Google Scholar

• 16

  KoppF.MendellJ. T. (2018). Functional classification and experimental dissection of long noncoding RNAs. Cell172 (3), 393–407. 10.1016/j.cell.2018.01.011

  • CrossRef
  • Google Scholar

• 17

  LamoureuxF.TrichetV.ChipoyC.BlanchardF.GouinF.RediniF. (2007). Recent advances in the management of osteosarcoma and forthcoming therapeutic strategies. Expert Rev. Anticancer Ther.7 (2), 169–181. 10.1586/14737140.7.2.169

  • CrossRef
  • Google Scholar

• 18

  LeeD. K. (2019). Survival analysis: Part II - applied clinical data analysis. Korean J. Anesthesiol.72 (5), 441–457. 10.4097/kja.19183

  • CrossRef
  • Google Scholar

• 19

  LiH.LanM.LiaoX.TangZ.YangC. (2020). Circular RNA cir-ITCH promotes osteosarcoma migration and invasion through cir-ITCH/miR-7/EGFR pathway. Technol. Cancer Res. Treat.19, 1533033819898728. 10.1177/1533033819898728

  • CrossRef
  • Google Scholar

• 20

  LiH.WangX.LuX.ZhuH.LiS.DuanS.et al (2019). Co- expression network analysis identified hub genes critical to triglyceride and free fatty acid metabolism as key regulators of age-related vascular dysfunction in mice. Aging (Albany NY)11 (18), 7620–7638. 10.18632/aging.102275

  • CrossRef
  • Google Scholar

• 21

  LiQ.LinJ.ZhangY.LiuX.ChenX. Q.XuM. Q.et al (2015). Differential behavioral responses of zebrafish larvae to yohimbine treatment. Psychopharmacol. Berl.232 (1), 197–208. 10.1007/s00213-014-3656-5

  • CrossRef
  • Google Scholar

• 22

  LiT.LiF.LinJ.ZhangY.ZhangQ.SunY.et al (2020). Deletion of c16orf45 in zebrafish results in a low fertilization rate and increased thigmotaxis. Dev. Psychobiol.62 (8), 1003–1010. 10.1002/dev.21984

  • CrossRef
  • Google Scholar

• 23

  LiuM.FanL.YanH.WangK.MaY.ShenL.et al (2020). A multi-model deep convolutional neural network for automatic hippocampus segmentation and classification in Alzheimer's disease. NeuroImage208, 116459. 10.1016/j.neuroimage.2019.116459

  • CrossRef
  • Google Scholar

• 24

  LorentM.GiralM.FoucherY. (2014). Net time-dependent ROC curves: A solution for evaluating the accuracy of a marker to predict disease-related mortality. Stat. Med.33 (14), 2379–2389. 10.1002/sim.6079

  • CrossRef
  • Google Scholar

• 25

  MassbergD.HattH. (2018). Human olfactory receptors: Novel cellular functions outside of the nose. Physiol. Rev.98 (3), 1739–1763. 10.1152/physrev.00013.2017

  • CrossRef
  • Google Scholar

• 26

  NakajimaK.KidaniT.MiuraH. (2021). Molecular profiling of bone remodeling occurring in musculoskeletal tumors. J. Orthop. Res.39 (7), 1402–1410. 10.1002/jor.24879

  • CrossRef
  • Google Scholar

• 27

  NazA.ZhangS.AnL.SongZ.ZiZ.WuJ.et al (2021). Muscle-specific programmed cell death 5 deletion attenuates cardiac aging. Int. J. Cardiol.345, 98–104. 10.1016/j.ijcard.2021.10.142

  • CrossRef
  • Google Scholar

• 28

  QuinnJ. J.ChangH. Y. (2016). Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet.17 (1), 47–62. 10.1038/nrg.2015.10

  • CrossRef
  • Google Scholar

• 29

  QuoC. F.KaddiC.PhanJ. H.ZollanvariA.XuM.WangM. D.et al (2012). Reverse engineering biomolecular systems using -omic data: Challenges, progress and opportunities. Brief. Bioinform.13 (4), 430–445. 10.1093/bib/bbs026

  • CrossRef
  • Google Scholar

• 30

  RitterJ.BielackS. S. (2010). Osteosarcoma. Ann. Oncol.21, vii320–5. 10.1093/annonc/mdq276

  • CrossRef
  • Google Scholar

• 31

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics26 (1), 139–140. 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 32

  RothzergE.XuJ.WoodD.KoksS. (2021). 12 Survival-related differentially expressed genes based on the TARGET-osteosarcoma database. Exp. Biol. Med.246, 2072–2081. 10.1177/15353702211007410

  • CrossRef
  • Google Scholar

• 33

  Sanchez CalleA.KawamuraY.YamamotoY.TakeshitaF.OchiyaT. (2018). Emerging roles of long non-coding RNA in cancer. Cancer Sci.109 (7), 2093–2100. 10.1111/cas.13642

  • CrossRef
  • Google Scholar

• 34

  SchaeferL.ReinhardtD. P. (2016). Special issue: Extracellular matrix: Therapeutic tools and targets in cancer treatment. Adv. Drug Deliv. Rev.97, 1–3. 10.1016/j.addr.2016.01.001

  • CrossRef
  • Google Scholar

• 35

  ShenP.ChengY. (2020). Long noncoding RNA lncARSR confers resistance to Adriamycin and promotes osteosarcoma progression. Cell Death Dis.11 (5), 362. 10.1038/s41419-020-2573-2

  • CrossRef
  • Google Scholar

• 36

  SimpsonE.BrownH. L. (2018). Understanding osteosarcomas. JAAPA31 (8), 15–19. 10.1097/01.JAA.0000541477.24116.8d

  • CrossRef
  • Google Scholar

• 37

  WangC. M.YeH. D.LiY. R.HongQ. X.TangL. L.ZhouA. N.et al (2015). Lack of an association between matrix metalloproteinase polymorphisms and coronary heart disease in a Han Chinese population. Genet. Mol. Res.14 (4), 12254–12261. 10.4238/2015.October.9.14

  • CrossRef
  • Google Scholar

• 38

  WangX.FangX.ZhengW.ZhouJ.SongZ.XuM.et al (2021). Genetic support of A causal relationship between iron status and type 2 diabetes: A mendelian randomization study. J. Clin. Endocrinol. Metab.106 (11), e4641–e4651. 10.1210/clinem/dgab454

  • CrossRef
  • Google Scholar

• 39

  WangX.JiaoX.XuM.WangB.LiJ.YangF.et al (2020). Effects of circulating vitamin D concentrations on emotion, behavior and attention: A cross-sectional study in preschool children with follow-up behavior experiments in juvenile mice. J. Affect. Disord.275, 290–298. 10.1016/j.jad.2020.06.043

  • CrossRef
  • Google Scholar

• 40

  WangY.ZengX.WangN.ZhaoW.ZhangX.TengS.et al (2018). Long noncoding RNA DANCR, working as a competitive endogenous RNA, promotes ROCK1-mediated proliferation and metastasis via decoying of miR-335-5p and miR-1972 in osteosarcoma. Mol. Cancer17 (1), 89. 10.1186/s12943-018-0837-6

  • CrossRef
  • Google Scholar

• 41

  WuB.WangZ.LinN.YanX.LvZ.YingZ.et al (2021). A panel of eight mRNA signatures improves prognosis prediction of osteosarcoma patients. Med. Baltim.100 (14), e24118. 10.1097/MD.0000000000024118

  • CrossRef
  • Google Scholar

• 42

  WuY.CaoH.BaranovaA.HuangH.LiS.CaiL.et al (2020). Multi-trait analysis for genome-wide association study of five psychiatric disorders. Transl. Psychiatry10 (1), 209. 10.1038/s41398-020-00902-6

  • CrossRef
  • Google Scholar

• 43

  XiangY.ZhangJ.LiQ.ZhouX.WangT.XuM.et al (2014). DNA methylome profiling of maternal peripheral blood and placentas reveal potential fetal DNA markers for non-invasive prenatal testing. Mol. Hum. Reprod.20 (9), 875–884. 10.1093/molehr/gau048

  • CrossRef
  • Google Scholar

• 44

  YanX.ZhaoX.LiJ.LinH.XuM. (2018). Effects of early-life malnutrition on neurodevelopment and neuropsychiatric disorders and the potential mechanisms. Prog. Neuropsychopharmacol. Biol. Psychiatry83, 64–75. 10.1016/j.pnpbp.2017.12.016

  • CrossRef
  • Google Scholar

• 45

  YangQ.WuF.MiY.WangF.CaiK.YangX.et al (2020). Aberrant expression of miR-29b-3p influences heart development and cardiomyocyte proliferation by targeting NOTCH2. Cell Prolif.53 (3), e12764. 10.1111/cpr.12764

  • CrossRef
  • Google Scholar

• 46

  YinC. D.HouY. L.LiuX. R.HeY. S.WangX. P.LiC. J.et al (2021). Development of an immune-related prognostic index associated with osteosarcoma. Bioengineered12 (1), 172–182. 10.1080/21655979.2020.1864096

  • CrossRef
  • Google Scholar

• 47

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 48

  YuH.PanR.QiY.ZhengZ.LiJ.LiH.et al (2020). LEPR hypomethylation is significantly associated with gastric cancer in males. Exp. Mol. Pathol.116, 104493. 10.1016/j.yexmp.2020.104493

  • CrossRef
  • Google Scholar

• 49

  YuS.ShaoF.LiuH.LiuQ. (2021). A five metastasis-related long noncoding RNA risk signature for osteosarcoma survival prediction. BMC Med. Genomics14 (1), 124. 10.1186/s12920-021-00972-5

  • CrossRef
  • Google Scholar

• 50

  ZeltzC.GullbergD. (2016). The integrin-collagen connection - a glue for tissue repair?J. Cell Sci.129 (6), 1284. 10.1242/jcs.188672

  • CrossRef
  • Google Scholar

• 51

  ZhangF.BaranovaA.ZhouC.CaoH.ChenJ.ZhangX.et al (2021b). Causal influences of neuroticism on mental health and cardiovascular disease. Hum. Genet.140 (9), 1267–1281. 10.1007/s00439-021-02288-x

  • CrossRef
  • Google Scholar

• 52

  ZhangL.HanB.LiuH.WangJ.FengX.SunW.et al (2021a). Circular RNA circACSL1 aggravated myocardial inflammation and myocardial injury by sponging miR-8055 and regulating MAPK14 expression. Cell Death Dis.12 (5), 487. 10.1038/s41419-021-03777-7

  • CrossRef
  • Google Scholar

• 53

  ZhangT.NieY.XiaH.ZhangY.CaiK.ChenX.et al (2020). Identification of immune-related prognostic genes and LncRNAs biomarkers associated with osteosarcoma microenvironment. Front. Oncol.10, 1109. 10.3389/fonc.2020.01109

  • CrossRef
  • Google Scholar

• 54

  ZhangY.WengQ.ChenJ.HanJ. (2020). Morusin inhibits human osteosarcoma via the PI3K-akt signaling pathway. Curr. Pharm. Biotechnol.21 (13), 1402–1409. 10.2174/1389201021666200416093457

  • CrossRef
  • Google Scholar

• 55

  ZhangY.WengQ.HanJ.ChenJ. (2020). Alantolactone suppresses human osteosarcoma through the PI3K/AKT signaling pathway. Mol. Med. Rep.21 (2), 675–684. 10.3892/mmr.2019.10882

  • CrossRef
  • Google Scholar

• 56

  ZhengS.ZhaoT.YuanS.YangL.DingJ.CuiL.et al (2019). Immunodeficiency promotes adaptive alterations of host gut microbiome: An observational metagenomic study in mice. Front. Microbiol.10, 2415. 10.3389/fmicb.2019.02415

  • CrossRef
  • Google Scholar

• 57

  ZhouX.LiQ.XuJ.ZhangX.ZhangH.XiangY.et al (2016). The aberrantly expressed miR-193b-3p contributes to preeclampsia through regulating transforming growth factor-β signaling. Sci. Rep.6, 19910. 10.1038/srep19910

  • CrossRef
  • Google Scholar

• 58

  ZhuK. P.ZhangC. L.MaX. L.HuJ. P.CaiT.ZhangL. (2019). Analyzing the interactions of mRNAs and ncRNAs to predict competing endogenous RNA networks in osteosarcoma chemo-resistance. Mol. Ther.27 (3), 518–530. 10.1016/j.ymthe.2019.01.001

  • CrossRef
  • Google Scholar

• 59

  ZiZ.SongZ.ZhangS.YeY.LiC.XuM.et al (2015). Rubicon deficiency enhances cardiac autophagy and protects mice from lipopolysaccharide-induced lethality and reduction in stroke volume. J. Cardiovasc. Pharmacol.65 (3), 252–261. 10.1097/FJC.0000000000000188

  • CrossRef
  • Google Scholar