Translate this page into:
Prognostic model for osteosarcoma: RNA diagnostics, tumor microenvironment and immunotherapy response
, Zhong Liua,
⁎Corresponding author: Hui-Xiang Tian. tianhuixiang99@163.com
⁎⁎Corresponding author: Zhong Liu. 104665296@qq.com
-
Received: ,
Accepted: ,
This article was originally published by Reed Elsevier India Pvt. Ltd. and was migrated to Scientific Scholar after the change of Publisher.
Abstract
Abstract
Osteosarcoma, an aggressive malignancy with poor prognosis, requires reliable prognostic models to predict treatment efficacy and toxicity. SnoRNAs, due to their stability and ubiquity, are promising biomarkers for tumor prognosis.
Gene expression and clinical data were downloaded from TCGA/GTEx/GEO. Then prognosis, survival analysis, gene differential expression, functional enrichment, immune cells infiltration, ESTIMATE/immune/stromal scores, and drug sensitivity analysis were performed. Finally, the drug sensitivity and gene expression were verified by in vitro experiments.
In this study, we established the 5 snoRNAs prognostic model (SNORA2B, SNORA12, SNORD99, SNORD123 and SNORD11B). It showed superior prognostic accuracy, but its prediction for the immune microenvironment was limited. Furthermore, the snoRNA-lncRNA-mRNA network was also used to construct a prognostic model comprising 4 RNAs (PARD6G-AS1, DLX2, TPD52 and GRAMD1B), which demonstrated significant correlation with immune-infiltrating tumor microenvironment. The drug sensitivity prediction of patients showed high consistency in the two models, and in vitro experiments proved that osteosarcoma cells were sensitive to simvastatin. Additionally, the differential expression of GRAMD1B, DLX2 and SNORD99 in the prognostic model was also validated between tumor and adjacent tissues.
Our findings revealed that independent snoRNAs have potential biological advantages as a prognostic biomarker, and the prognostic model combined with other RNAs has offered more information for immune microenvironment and immune function. The drug sensitivity prediction of patients showed high consistency in the two models. Overall, our study provides reliable prognostic and immune-related biomarkers for osteosarcoma patients, and valuable insights for therapeutic strategy and therapy toxicity.
Keywords
Osteosarcoma
SnoRNA
Prognostic model
Immune microenvironment
Drug sensitivity
1 Introduction
Osteosarcoma (OS) is a malignancy that originates from mesenchymal tissue,1 and accounts for about 35 % of all primary malignant bone tumors. It primarily affects children and adolescents, and poses a significant clinical challenge.2 Although extensive treatments such as chemotherapy, surgery, immunotherapy and targeted therapies still have some side effects, the overall survival rate of OS patients has increased to 70 %, and the prognosis has also improved.3 Despite these achievements, metastasis remains a primary cause of poor prognosis, contributing to the persistently low five-year survival rate of OS patients, which currently stands at only 20 %.4 Therefore, the establishment of new prognostic models and drug prediction models are still a research focus for reducing cancer therapy toxicity and improving the survival rate of OS patients.
RNA was first discovered in the 1940s and played an important role in protein translation.5 However, it was not until 1950s when the non-coding RNAs (ncRNAs) such as ribosomal RNAs (rRNAs) and transfer RNAs (tRNAs) were first discovered, establishing that the RNAs were not just an intermediary in the process of protein production.6,7 In 1960, researchers discovered small nucleolar RNAs (snoRNAs) in mammalian cells,8 a class of small non-coding RNA widely present in the nucleolus of eukaryotic cells, with a length of 60-300 nt, which can combine with nucleolar ribonucleoproteins to form small nucleolar ribonucleoprotein complexes (snoRNPs).9 SnoRNAs are mainly involved in the processing of rRNAs, regulating the messenger RNAs (mRNAs) splicing and translation, and responding to oxidative stress.10 Increasing reports have shown that snoRNAs are abnormally regulated in multiple tumors and are also involved in the process of genetic diseases, hematopoiesis and metabolism.11,12 For example, SNORD88C promoted the proliferation and metastasis of non-small cell lung cancer by enhancing SCD1 translation and inhibiting autophagy.13 Moreover, SNORA13 and SNORA28 could reduce the cytotoxicity of doxorubicin in osteosarcoma.14 Therefore, based on the diversity of snoRNAs in cellular functions, some researchers believe that tumor-specific snoRNAs may serve as cancer-related biomarkers.15 For instance, SNORA42 was proved to be associated with poor prognosis in prostate cancer.16 And in the non-small cell lung cancer, SNORD33, SNORD66, and SNORD76 are significantly upregulated in both tumors and plasma, offering the potential to differentiate individuals with non-small cell lung cancer based on these snoRNAs.17 Taken together, these findings suggest that snoRNAs hold promise as biomarkers for both diagnosis and prognosis of cancer.
SnoRNAs not only perform the function of splicing mRNAs, but also have a close association with other types of ncRNAs with a wide range of functional similarities. In 2012, Chen's research group at the Shanghai Institute of Biochemistry first reported a box H/ACA snoRNA-ended long noncoding RNA (lncRNA) that enhances pre-rRNA transcription. These RNAs sequences contain complete snoRNAs sequences and are essential for the stability and subcellular localization of sno-lncRNAs.18 Additionally, studies have revealed unexpected interactions between snoRNAs and mRNAs.19 The other study found that the small nucleolar RNA host gene 1 (SNHG1) interacted with miR-154-5p, affecting colorectal cancer cell growth.20 These studies indicated that snoRNAs may play an important role in regulating the interactions between different types of RNAs. Similarly, there are studies suggesting a potential role for snoRNAs in the regulation of the immune microenvironment. For example, SNORD46 inhibitors improved viability of obese NK and anti-tumor immunity of CAR-NK cell therapy.21 Cai et al. identified tumour immune infiltration-associated snoRNAs and predicted their response to immunotherapy.22
In our study, we identified the prognostic snoRNAs, lncRNAs and mRNAs in OS using the RNA sequencing data obtained from the Cancer Genome Atlas (TCGA) database. Specifically, we developed two prognostic models derived from snoRNAs and integrated sonRNA/lncRNA/mRNA correlation network, respectively, and compared the function of prognostic signatures. Our study confirmed the advantages of snoRNAs in prognostic accuracy and found the value of comprehensive snoRNAs, lncRNAs and mRNAs in predicting immune microenvironment and immune function. In addition, we also predicted the drug sensitivity of patients, which has important indicative value for the cancer therapy toxicity and effectiveness. Finally, the drug sensitivity and the expression of prognostic genes were verified through in vitro experiments. Taken together, our results offer new perspectives into the molecular mechanisms underlying OS, which may inform new avenues for clinical research, therapeutic interventions and reduction the incidence of drug toxicity. Additionally, our study provides further considerations for the establishment of disease prognosis models.
2 Materials and methods
2.1 Data sources
RNA sequencing data was retrieved from TCGA (https://portal.gdc.cancer.gov/)23 for 88 OS samples from the TCGA-OS cohort. The quantification of gene expression profiles were achieved via the measurement of fragments per kilobase of transcript per million mapped reads (FPKM). Matched clinicopathological data, including age, gender, primary tumor site, metastasis, and survival information were collected to facilitate comprehensive analysis of the datasets. Samples with complete survival data and clinicopathological characteristics were used for subsequent analysis. Then, the transcriptome data of snoRNAs, mRNAs, miRNAs and lncRNAs were extracted and low abundance genes with an average expression level <0.05 were filtered. Additionally, transcriptomic data from normal muscle tissues samples in the GTEx V8 cohort (https://www.gtexportal.org/home/) were downloaded and integrated with osteosarcoma data from TCGA for differential expression analysis. The GSE253548 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE253548) from the GEO database, comprising 50 osteosarcoma samples and 40 healthy tissues samples, was also incorporated. Differential expression analysis was performed using the limma package.
2.2 Construction of snoRNA-related competing endogenous RNA (ceRNA) network
Univariate Cox analysis was employed to detect prognostic significance of snoRNAs. Further, mRNAs and lncRNAs with p < 0.05 and cor >0.5 in relation to the snoRNAs. Prognostic-related mRNAs and lncRNAs were picked out to construct the prognostic-related snoRNA-lncRNA-mRNA network. To explore potential mechanisms underlying this network, ceRNA network was constructed based on the fact that lncRNAs can interact with micro RNAs (miRNAs) to modulate mRNAs activity. The ENCORI database (starBase v3.0, http://starbase.sysu.edu.cn/)24 stores millions of RNA-RNA interactions. Thus, the target miRNAs related to lncRNAs were directly predicted in the ENCORI database. The overlapping miRNA-mRNA relationship pairs in the miRmap, microT, miRanda, TargetScan databases from ENCORI were selected to obtain mRNA-related miRNAs. Cytoscape software visualized the prognosis-related snoRNA-lncRNA-mRNA correlation network and ceRNA network. The snoRNA-related ceRNA network, which was obtained by the intersection of the above two networks, was used for subsequent analysis. The visualization of this network was accomplished using the “ggalluvial” R package (version: 0.9.1).
2.3 Model establishment for prognosis prediction
To assess the prognostic significance of snoRNAs and snoRNA-associated ceRNA networks, the “glmnet” R package was employed to perform LASSO regression analysis and construct a prognostic model. An associated regression coefficient for each gene was used in the model to determine patients' risk score, which was calculated using the formula: ∑ (βi × Expi), where βi represents the coefficient corresponding to the gene and Expi represents the gene expression level. Patients were then classified into low-risk or high-risk group, based on the median risk score as the cut-off point. To assess the difference in overall survival between the low-risk and high-risk groups, we utilized the Kaplan-Meier (K-M) method and performed log-rank tests. This analysis helps determine whether the risk stratification based on the prognostic model significantly impacts overall survival. To evaluate the performance of the prognostic model in predicting overall survival, we employed the “timeROC” package to plot the receiver operating characteristic (ROC) curve and calculate the area under the curve (AUC). Additionally, we utilized the “scatterplot3d” package to generate principal component analysis (PCA) plots, visually demonstrating the effectiveness of our risk stratification approach.
2.4 Evaluation of prognostic models
Univariate and multivariate Cox regression analyses to determine if the risk score was an independent prognostic factor for overall survival. Furthermore, we assessed the relationship between the risk score and clinicopathological characteristics. The “rms” R package was used to construct nomograms for predicting 1-, 3- and 5-year overall survival, incorporating all independent prognostic factors identified by multivariate Cox regression analysis. Calibration curves were plotted to assess the nomograms' predictive performance.
2.5 Enrichment analysis
To investigate the biological functions and pathways associated with low-risk and high-risk groups, Wilcoxon test was employed to identify differential expressed genes (DEGs) between these two groups. Subsequently, the “clusterProfiler” R package performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.25
2.6 Assessment of immune infiltration and tumor microenvironment (TME)
To explore the potential association between the risk score and TME, ssGSEA was used to evaluate the level of immune cell infiltration and the activity of immune-related pathways in both the high-risk and low-risk groups.26 Moreover, we employed CIBERSORT to determine the abundance of 22 different immune cell subsets in heterogeneous samples, and investigated the correlation between these immune cells and the prognostic signatures.27 We evaluated the disparities in Stromal score, Immune score, and ESTIMATE score between the high-risk and low-risk groups by utilizing the powerful “estimate” R package.28
2.7 Drug sensitivity and immunotherapy response prediction
To screen out potentially effective drugs for the treatment of OS and evaluate the effectiveness of these two prognostic models for clinical application, we used the “oncopredict” R package to estimate half-maximal inhibitory concentrations (IC50) of drugs, and performed statistical analyzes of differences between the two risk groups. Furthermore, Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu), an online site that predicts immunotherapy efficacy and immune escape, was used to calculate the TIDE score of OS.
2.8 Detection of cell proliferation
Human osteosarcoma cell lines SaOS-2 and U2OS, obtained from Cell Bank, Chinese Academy of Sciences, were cultured in high-glucose DMEM media supplemented with 10 % foetal bovine serum, and 1 % penicillin and streptomycin. hBMSC cell line was purchased from Cyagen Biosciences Co., LTD. It was cultured in MEMα supplemented with 10 % FBS. Cells were planted in 96-well plates and then treated with different concentrations of simvastatin (5–80 μmol/L) for 24h, 48h or 72h. Cell viability was assayed by adding 10 % CCK8 solution and cell viability was calculated based on the absorbance.
2.9 Real-time quantitative polymerase chain reaction (RT-qPCR)
Osteosarcoma cells and hBMSC were collected. RNA was extracted by a Total RNA Rapid Extraction Kit. RNA was reverse transcribed and the expression of each gene was determined using SYBR Green and specific primers. U6 was used as reference gene for snoRNA expression, and GAPDH was used as reference gene for lncRNA and mRNA expression. The primers were listed in Supplementary Table S1.
2.10 Statistical analysis
All statistical analysis and graphs were performed using R 4.2.1. The Wilcoxon test was used to compare the data of the high-risk and low-risk groups respectively. The correlational assessment was performed using the “cor_test” and “stat_compare_means” modules with default parameters. The univariate and multivariate Cox regression model estimated hazard ratios (HR) and 95 % confidence intervals (CI). If not specified above, p < 0.05 was considered statistically significant.
3 Results
3.1 Discovery of prognostic snoRNAs
The study design workflow is depicted in Fig. 1. We obtained transcriptome data and paired clinical information of 85 OS patients for subsequent analysis. Table S2 summarized the clinical parameters of the TCGA-OS patients. We analyzed the transcriptome data and extracted the data of snoRNAs, lncRNAs, mRNAs and miRNAs to construct the heatmap (Fig. S1). It was found that the prognosis of OS patients was greatly related to 18 snoRNAs by univariate Cox regression analysis (P < 0.05, as demonstrated in Fig. 2A). It should be emphasized that all 18 snoRNAs were identified as predictors for OS patients (HR > 1). This finding strongly implied that they had the potential to serve as valuable prognostic markers.
3.2 Construction of snoRNA prognostic model
Based on the 18 prognostic snoRNAs, we construct the OS prognostic model. Employing the LASSO regression with the best fit regression coefficient and 10-fold cross-validation, the prognostic signatures of 5 snoRNAs (SNORA2B, SNORA12, SNORD99, SNORD123 and SNORD11B) were identified (Fig. 2B and C). The risk score of prognostic signatures was calculated using the following formula: risk score = (0.054 × SNORA2B) + (0.039 × SNORA12) + (0.031 × SNORD99) + (0.108 × SNORD123) + (0.070 × SNORD11B). We employed the survival curve to evaluate the prognostic signature's predictive ability and the results showed that patients in the high-risk group had significantly worse overall survival than those in the low-risk group (P < 0.001) (Fig. 2D). At 1-, 3- and 5-year time points, AUC scores of the ROC graph were determined to be 0.739, 0.802, and 0.856, respectively (Fig. 2E). PCA further confirmed the excellent discriminatory performance of the prognostic signatures within the TCGA-OS cohort (Fig. 2F). These findings collectively indicate that the snoRNA-based prognostic model exhibits a robust ability to accurately predict the prognosis of OS patients.
3.3 Immune infiltration analysis of snoRNAs prognostic model
It has been demonstrated in previous studies that TME, particularly in relation to the immune cell infiltration, plays a crucial role in tumor progression. Hence, we conducted further inquiry into the correlation between the prognostic signatures and TME. Initially, we utilized the ssGSEA algorithm to discern dissimilarities in immune cells and immune functionality between cohorts categorized as low-risk and high-risk. Unfortunately, our model failed to accurately reflect the differences in immunocytes and immunoregulation between the high-risk and low-risk groups (Fig. 2G and H). We further used the CIBERSORT algorithm for validation analysis and observed that the risk score was related to naive T cells CD4 memory and resting dendritic cells (Fig. 2I). Finally, we validated our results using the ESTIMATE algorithm. The Stromal score, Immune score, and ESTIMATE score had no significant difference in the high-risk and low-risk groups (Fig. 2J). These results demonstrate that while the snoRNA-based prognostic model alone can effectively predict the patients' prognosis, it still has limitations in accurately distinguishing the immune landscape within the TME.
3.4 Establishment of snoRNA-lncRNA-mRNA network related to prognosis
Although the potential significance of snoRNAs in osteosarcoma (OS) has been recognized, their underlying mechanisms have not been fully elucidated. To further explore the mechanism of snoRNAs in OS, we included lncRNAs, mRNAs and snoRNAs for OS prognostic analysis. Firstly, we conducted a correlation analysis between the 18 prognosis-related snoRNAs and lncRNAs and mRNAs. Our findings revealed that among the 18 snoRNAs, 17 of them exhibited 256 co-expression relationships with 169 lncRNAs and 470 co-expression relationships with 365 mRNAs (Fig. S2). Notably, all of these co-expression relationships were positively correlated. Subsequently, we performed univariate Cox analysis on these 169 lncRNAs and 365 mRNAs, and found that 74 lncRNAs and 115 mRNAs were significantly related with OS prognosis (Table S3). Moreover, 12 of the prognostic snoRNAs were associated with both the prognosis-related lncRNAs and mRNAs. Therefore, we established a relevance network graph consisting of 12 snoRNAs - 73 lncRNAs - 108 mRNAs which were the OS prognosis-related RNAs (Fig. S3).
3.5 Construction of snoRNA-related ceRNA network
As it is well established, miRNAs can effectuate gene silencing through binding to mRNAs, whereas lncRNAs modulate gene expression by competing with miRNAs. Such ceRNA crosstalk is widely observed in different biological processes and diseases.29,30 Using this information, we developed a ceRNA network consisting of lncRNAs and mRNAs that are associated with prognosis. A total of 71,952 lncRNA-miRNA relationship chains were downloaded from the ENCORI database, and 685 relationship chains between 21 lncRNAs and 400 miRNAs were found among the 73 lncRNAs (Fig. S4). In the meantime, we downloaded 38,580 mRNA-miRNA relationship chains that concomitantly exist in miRmap/microT/miRanda/TargetScan databases, and obtained 307 relationship chains between 46 mRNAs and 130 miRNAs among the 108 mRNAs (Fig. S5). Of these, 96 miRNAs were predicted simultaneously. In order to accurately demonstrate the functional mechanism of snoRNAs, a ceRNA network containing 113 axes was constructed by 12 lncRNAs, 62 miRNAs, and 20 mRNAs. And in the network, the lncRNAs and mRNAs were significantly correlated with at least one same snoRNA (Fig. S6).
3.6 Construction of snoRNA-lncRNA-mRNA network prognostic model
Considering that mRNAs, lncRNAs, and snoRNAs all have unique molecular mechanisms, we integrated the prognostic-related snoRNAs, lncRNAs and mRNAs to construct the OS prognostic model. The prognostic signatures of 4 RNAs (PARD6G-AS1, DLX2, TPD52 and GRAMD1B) were identified (Fig. 3A and B) and the risk score of prognostic signatures was calculated using the following formula: risk score = (0.015 × PARD6G-AS1) + (0.024 × DLX2) + (0.018 × TPD52) + (0.078 × GRAMD1B). To assess the predictive power of our prognostic signatures, we utilized survival analysis. The results indicated that patients in the high-risk group had significantly worse overall survival than those in the low-risk group (P < 0.001) (Fig. 3C). Scatter plots and risk curves exhibited that patients in the high-risk group had elevated risk scores and higher mortality rates (Fig. 3D). Furthermore, compared with the snoRNAs prognostic model, the snoRNA-lncRNA-mRNA prognostic model exhibited a weaker predictive accuracy for the AUC value of ROC curve. We calculated the AUC scores for the ROC graph using the TCGA-OS cohort. Our results showed that the AUC values were 0.698, 0.724, and 0.716 for the 1-, 3- and 5-year time points, respectively (Fig. 3E). PCA confirmed the excellent discriminatory performance of the prognostic signatures within the TCGA-OS cohort (Fig. 3F). Finally, the risk heatmap was employed to visualize the expression of genes in the model for the low-risk and high-risk groups (Fig. 3G). Notably, we observed that the 4 genes included in the model were highly expressed in the high-risk group (P < 0.05) (Fig. 3H). Taken together, these findings showed that the snoRNA-lncRNA-mRNA prognostic signatures also demonstrated the ability to accurately predict the prognosis of OS patients.
3.7 Independent prognostic analysis and clinical relevant analysis
Both univariate and multivariate Cox regression analysis were performed to ascertain independent risk factor associated with OS prognosis. The univariate Cox regression analysis exposed that metastasis (P < 0.001, HR = 4.740, 95 % CI: 2.271–9.895) and risk score (P < 0.001, HR = 5.854, 95 % CI: 1.523–2.536) were both significant prognostic factors for OS (as shown in Fig. S7A). Furthermore, the multivariate Cox regression analysis demonstrated that metastasis (P < 0.001, HR = 3.697, 95 % CI: 1.526–7.892) and risk score (P < 0.001, HR = 4.117, 95 % CI: 1.281–2.196) remained as independent prognostic factors for OS (Fig. S7B). These findings indicated that the prediction of OS prognosis was not affected by other clinical factors in the model. Furthermore, we performed a correlation analysis of the risk score and four clinical parameters, namely age, gender, primary tumor site, and metastasis, using the OS cohort. It showed that the risk score was not significantly correlated with any of the aforementioned clinical factors, which was consistent with the findings of the independent prognosis (Table S4). This suggested that the risk score had the potential to predict OS prognosis without being confounded by the clinicopathological characteristics of patients.
3.8 Nomogram and calibration curves
To further improve individual prognosis prediction, researchers often incorporate conventional clinical parameters into the model and establish a nomogram. In our study, we identified two independent prognostic factors and utilized them to develop a nomogram that predicted the likelihood of 1-, 3- and 5-year overall survival in patients (Fig. S8A). Each prognostic parameter was assigned a score, and the sum of the two scores corresponded to the predicted 1-, 3- and 5-year overall survival. An enhanced total score indicated a worse prognosis. Additionally, the calibration curve for the 1-, 3- and 5-year overall survival demonstrated comparable performance to the ideal model, suggesting that the nomogram had strong predictive discrimination and accuracy (Fig. S8B). The results indicated that the nomogram we developed could be a useful tool in the clinical management of OS.
3.9 Enrichment analysis
We analyzed the differences between low-risk and high-risk groups and identified 602 DEGs, as shown in Fig. S9A. Subsequently, we performed GO and KEGG enrichment analysis on the DEGs, with P < 0.05 considered statistically significant. KEGG enrichment analysis revealed significant enrichment of neuroactive ligand-receptor interaction and calcium signaling pathway, as depicted in Fig. S9B. In GO enrichment analysis, the DEGs were enriched in biological processes (BP) related to extracellular matrix organization and extracellular structure organization. Regarding cellular component (CC), these genes exhibited significant enrichment in collagen-containing extracellular matrix and endoplasmic reticulum lumen. The molecular function (MF) enrichment analysis further identified enrichment in receptor ligand activity and signaling receptor activator activity (Fig. S9C).
3.10 Correlation of TME and prognostic signatures
Due to the limitations of snoRNAs prognostic model in immune prediction, we further explored the ability of snoRNA-lncRNA-mRNA to predict immune infiltration. At first, we also utilized the ssGSEA algorithm to discern dissimilarities in immune cells and immune functionality. Our findings indicated a marked reduction in the concentration of CD8+ T cells, Tfh cells, and Th2 cells within the high-risk group, in contrast to the low-risk group (Fig. 4A). Regarding immune functions, our analysis demonstrated a noteworthy decline in APC_co_inhibition, CCR, Inflammation-promoting, Parainflammation, T_cell_co-inhibition, T_cell_co-stimulation and Type_II_IFN_ response within the high-risk group (Fig. 4B). This observation suggested that there might be a relationship between high-risk patients and inhibitory immune microenvironment. In addition, employing the CIBERSORT algorithm, we identified a robust association between the risk score and immune infiltration cells (Fig. 4C), including naive B cells, CD8+ T cells, gamma delta T cells, resting Dendritic cells and Neutrophils. We conducted a correlation analysis between immune-related cells and risk scores and discovered a positive correlation between the Immune score and naive B cells, CD4 naive T cells as well as resting dendritic cells. However, we found that Tregs, T follicular helper cells and CD8+ T cells were negatively associated with the Immune score (Fig. 4D). Importantly, we validated our results using the ESTIMATE algorithm and found that Stromal score, and ESTIMATE score were significantly lower in the high-risk group compared to the low-risk group (Fig. 4E). Nonetheless, we found no significant difference of the Immune score, which might be attributed to the immune evasion and immunosuppression characteristics of OS. We then explored how the prognostic signatures were associated with different types of immune cells. Our primary results showed that the expression levels of DLX2, GRAMD1B, and PARD6G-AS1 exhibited a positive correlation with resting B cells and resting dendritic cells, while a negative correlation with CD8+ T cells. (Fig. 4F). These findings were consistent with our previous results and suggested that these genes might modulate the tumor-infiltrating microenvironment, ultimately affecting tumor growth and progression. Despite that, we did not find any significant relationship between TPD52 and any of the immune cells that were analyzed.
3.11 Prediction of drug sensitivity and immunotherapy response
Furthermore, we also predicted the drug sensitivity of the patients. We scored and compared 545 drugs in the high-risk and low-risk groups using both the snoRNAs prognostic model and the snoRNA-lncRNA-mRNA network prognostic model. Our analysis revealed that 77 drugs exhibited differences in sensitivity between the high-risk and low-risk groups in the snoRNAs prognostic model, while there were 111 drugs in the snoRNA-lncRNA-mRNA network prognostic model. Interestingly, 25 drugs showed consistent results in both prognostic models (Fig. S10), and some representative drugs are shown in Fig. 5A and B. Among them, patients in the high-risk group showed lower drug sensitivity to oxaliplatin, belinostat etc., indicating that the therapy toxicity in high-risk patients may be higher. Several statins (simvastatin and lovastatin) were also found to show potential benefit in patients in the high-risk group. In addition, we verified the sensitivity of simvastatin to osteosarcoma cells. The results showed that simvastatin had a certain tumor cell inhibitory effect (Fig. 5C). Finally, the TIDE scores were calculated. The results showed that patients in the high-risk group had a higher response to immunotherapy in the snoRNA-lncRNA-mRNA network prognostic model, which was consistent with the prediction of the immune microenvironment. However, the snoRNAs prognostic model does not distinguish the immunotherapy response of patients in the high- and low-risk groups (Fig. 5D).
3.12 Prognostic RNAs verification
Finally, to validate the expression of these prognostic RNAs between tumor and normal tissues, we conducted differential expression analysis using RNA expression data from TCGA, GTEx and GEO. Firstly, we extracted and filtered RNA expression data from 803 normal muscle tissues in the GTEx database and combined them with RNA data from 85 osteosarcoma patients’ tumor tissues used in the previous analysis for normalization and comparative analysis (Fig. S11). The results found that SNORD99 in the snoRNAs prognostic model was significantly upregulated in tumor tissues, while DLX2 and GRAMD1B from the snoRNA-lncRNA-mRNA prognostic model also exhibited higher expression in tumor tissues (Fig. 6A). Additionally, the GSE253548 dataset in the GEO database was identified, which includes sequencing data from 50 osteosarcoma tumors and 40 healthy tissues. After initial normalization and correction of the data, the results revealed 2313 DEGs between tumor and adjacent tissues. None of the five snoRNAs in the snoRNAs prognostic model showed significant expression differences, while GRAMD1B from the snoRNA-lncRNA-mRNA prognostic model were found to be highly expressed in tumor tissues (Fig. 6B). Notably, GRAMD1B showed consistent results across both datasets. To further validate these findings, we included hBMSC as normal control cells and analyzed the expression of RNAs from the two prognostic models in normal versus osteosarcoma cells. The results showed that SNORD99, DLX2 and GRAMD1B had higher expressions in osteosarcoma cells (Fig. 6C).
4 Discussion
Osteosarcoma (OS) is the most prevalent primary bone malignancy known for its invasive and metastatic capabilities.4,31 Despite notable progress in the treatment of OS over recent years, the overall survival has hit a plateau, hindering further progress.32 More and more studies have used high-throughput sequencing technology and bioinformatics methods to find prognostic markers33 and predict the efficacy34 and therapy toxicity35 of OS, which opens up a new avenue for the therapy and investigation of OS. However, there is still a lack of reliable biomarkers and prognostic models to help improve clinical outcomes.36
Since its discovery, RNAs have been widely recognized to be involved in the whole stage of disease development, and have important potential research value for drug response and patient prognosis. With the advancement of detection techniques, some small RNAs have been discovered, and snoRNA is one of them.37 Compared with other RNAs, the metabolism of snoRNAs is very stable, and they can be widely detected in tissues, cells, serum, sputum and urine.10,38 Especially for the body fluid examination in clinical, this low invasive, simple detection method, low cost and convenient census make snoRNAs of great value in clinical detection. In recent years, snoRNAs have emerged as the promising predictive biomarkers, a diagnostic and prospective therapeutic target for various cancers.22,39,40 However, their potential in OS is yet to be fully explored and requires further investigation. Therefore, in our study, we first focused on the research value of snoRNAs in the OS prognosis. To this end, we screened 18 snoRNAs related to prognosis through TCGA-OS data, of which SNORD1B was identified as a tumor-promoting RNA for non-small cell lung cancer.41 SNORA14B was abnormally expressed in pancreatic cancer patients' serum.42 Moreover, Kangyu Wang et al. regarded SNORD83A in plasma as a potential biomarker for early diagnosis of non-small cell lung cancer.43 SNORD3B-1 was identified as an immune-related gene in acute myeloid leukemia.44 It is reported that rpL13a snoRNAs are induced in response to LPS-mediated liver injury and the knockdown of the snoRNAs protects against H2O2 cytotoxicity.45 Although there are many reports on the role of snoRNAs in disease diagnosis and cancer therapy toxicity, and the mechanism of action involves a wide range of processes such as immune regulation,46 metabolism47 and epigenetic modification,48 the specific mechanism has not been discussed in detail. To further explore, we must know the origin, aetiology and carcinogenic mechanism of the diseaseand.49 This will provide an important reference for early disease screening, tumour prognosis and drug sensitivity prediction.
In addition, we also established sonRNAs prognostic model. It contained only 5 RNAs (SNORA2B, SNORA12, SNORD99, SNORD123 and SNORD11B) and had strong prognostic accuracy. Among them, SNORA12 has been selected as one of the diagnostic markers for the identification of cervical cancer.50 Some studies have shown that some SNORNAs, including SNORD99, may be related to the pathogenesis of Alzheimer's disease.51 SNORD123 was found to be associated with age and immune stimulation.52 Moreover, SNORD11B was proved to promote cell proliferation and invasion and inhibits apoptosis in colorectal carcinogenesis.53 For the prediction of immune landscape, the model was unable to differentiate between differences in the immune microenvironment of high-risk and low-risk groups. It may be due to the limited annotation of snoRNAs function at present.54 It is also possible that as the ncRNAs, snoRNAs don't have direct function, and they may need to cooperate with other RNAs to exert micro-regulation in the body.55
Accordingly, we constructed the snoRNA-lncRNA-mRNA network related to prognosis. Given the correlation among the snoRNAs, lncRNAs and mRNAs, we believe that it is possible to find prognostic signals with clinical indicative value from the snoRNA-lncRNA-mRNA network. Finally, a prognostic model including 4 RNAs (PARD6G-AS1, DLX2, TPD52, and GRAMD1B) was established to predict the prognosis of OS patients in this study. It is reported that DLX2 promotes epithelial-mesenchymal transition and doxyorubicin resistance in osteosarcoma by affecting HOXC8 together with CDH2.56 TPD52 could be a possible new disease biomarker and therapeutic target because TPD52 knockdown inhibited the migration and invasive of prostate cancer cells and the growth of tumor.57 GRAMD1B, considered as a novel biomarker in breast cancer, inhibited cell migration by JAK/STAT and Akt signaling.58 In our study, the model's prediction accuracy of 1-, 3- and 5-year overall survival was lower than the snoRNAs prognostic model, which suggested the advantage of snoRNAs as prognostic markers compared with complex RNA systems. Additional model evaluation indicated that both tumor metastasis and risk score were recognized as independent prognostic factors for OS survival. The functional enrichment analysis of prognostic characteristics revealed that they were closely related to immune-related functional items. We further used ssGSEA, Cibersort and ESTIMATE algorithms, determining that prognostic characteristics were significantly correlated with immune-infiltrating TME.
Meanwhile, we also predicted and verified the drug sensitivity of this model. Current advances in drug therapy for OS, including immunotherapy, chemotherapy, and targeted drugs, have significantly improved the survival and prognosis of patients.32 The value of snoRNAs in drug sensitivity prediction has been demonstrated.14 This suggests a direct link between snoRNA and drug response. On the other hand, drug treatments may lead to a series of toxicities,59 such as ototoxicity, nephrotoxic and myelosuppression, which are also associated with snoRNAs. Therefore, it is important to find therapeutic drugs that are more suitable for patients and reduce the toxicity caused by drug treatment. In this study, among 545 drugs, we predicted that cabozantinib, simvastatin and lovastatin may be candidate drugs for high-risk group patients, and oxaliplatin, belinostat and lomeguatrib may not be ideal for high-risk group patients, which may be helpful for clinical effective drug use and improving toxic events in tumor therapy. In addition, the TIDE score showed that the high-risk group of snoRNA-lncRNA-mRNA network prognostic model was more sensitivity to immune checkpoint inhibitors. Although the use of immune drugs in OS patients is limited at present, this undoubtedly gives some hints that immune checkpoint inhibitors may be promising for the high-risk population.
In light of these findings, we believe that snoRNAs have a unique advantage in indicating the OS prognosis, cancer therapy efficacy and toxicity. It may be due to their unique stability and tissue specificity. But their prediction of TME still lack sufficient evidence. The other functional RNAs such as mRNAs and lncRNAs, are related to immunity and also show an intimate correlation with snoRNAs. Prognostic analysis of multiple prognosis-related RNAs in vivo may have better predictive sensitivity for immune-related regulation and response in patients.
Finally, this study incorporated the GSE253548 dataset from the GEO database and RNA data from normal muscle tissues in the GTEx database. The differences in expression patterns of these nine RNAs between normal and tumor tissues were investigated, and the results were validated by qPCR. Ultimately, GRAMD1B, DLX2, and SNORD99 were confirmed not only to hold significant prognostic value for patients but also to exhibit potential in distinguishing osteosarcoma. This suggests their potential involvement in the pathological mechanisms of osteosarcoma, warranting further exploration.
5 Conclusions
In summary, we identified two prognostic models for OS based on snoRNAs and snoRNA-lncRNA-mRNA network, and compared their performance. It is confirmed that the snoRNAs prognostic model had better prognostic accuracy, and the model of comprehensive snoRNA-lncRNA-mRNA had more advantages in predicting the TME of patients. Moreover, the two models showed high consistency in drug sensitivity prediction, including the prediction results of statins. Among them, GRAMD1B, DLX2, and SNORD99 were validated to demonstrate differential expression between tumor and adjacent tissues. This study provides prospective insights for the construction of prognostic models, drug sensitivity and therapy toxicity prediction with RNA disease.
Guardian/patient's consent
Not applicable.
Availability of data and materials
The datasets analyzed for this study can be found in the TCGA-OS project (https://portal.gdc.cancer.gov/), GTEx V8 cohort (https://www.gtexportal.org/home/) and GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE253548).
Consent for publication
All authors agreed to publish the study.
Ethical statement
Not applicable.
Credit author statement
Ding-Chao Rong: Concept and design, Data acquisition, Data analysis. Xu Rong: Data acquisition, Statistical analysis. Yuan-Shen Chen: Data analysis. Yu-Ligh Liou: Statistical analysis. Jie Mei: Concept and design. Hui-Xiang Tian: Literature search, Manuscript preparation, Manuscript editing. Zhong Liu: Funding supporting and Manuscript review. All authors approved the final manuscript.
Funding statement
This work was supported by Hunan Provincial Natural Science Foundation of China [Grant No. 2022JJ50192] and Science and Technology Plan Project of Shaoyang City [Grant No. 2024PT4054].
References
- DNA damage response and repair in osteosarcoma: defects, regulation and therapeutic implications. DNA Repair (Amst). 2021;102
- [Google Scholar]
- Immunotherapy for osteosarcoma: fundamental mechanism, rationale, and recent breakthroughs. Cancer Lett. 2021;500:1-10.
- [Google Scholar]
- History, discovery, and classification of lncRNAs. Adv Exp Med Biol. 2017;1008:1-46.
- [Google Scholar]
- snoRNAs: functions and mechanisms in biological processes, and roles in tumor pathophysiology. Cell Death Discov. 2022;8(1):259.
- [Google Scholar]
- Emerging functions for snoRNAs and snoRNA-Derived fragments. Int J Mol Sci. 2021;22(19)
- [Google Scholar]
- Beyond microRNA--novel RNAs derived from small non-coding RNA and their implication in cancer. Cancer Lett. 2013;340(2):201-211.
- [Google Scholar]
- SNORD88C guided 2'-O-methylation of 28S rRNA regulates SCD1 translation to inhibit autophagy and promote growth and metastasis in non-small cell lung cancer. Cell Death Differ. 2023;30(2):341-355.
- [Google Scholar]
- Small nucleolar RNAs determine resistance to doxorubicin in human osteosarcoma. Int J Mol Sci. 2020;21(12)
- [Google Scholar]
- Targeting snoRNAs as an emerging method of therapeutic development for cancer. Am J Cancer Res. 2019;9(8):1504-1516.
- [Google Scholar]
- SNORA42 enhances prostate cancer cell viability, migration and EMT and is correlated with prostate cancer poor prognosis. Int J Biochem Cell Biol. 2018;102:138-150.
- [Google Scholar]
- Small nucleolar RNA signatures as biomarkers for non-small-cell lung cancer. Mol Cancer. 2010;9:198.
- [Google Scholar]
- SLERT regulates DDX21 rings associated with pol I transcription. Cell. 2017;169(4):664-678.
- [Google Scholar]
- The long noncoding RNA SNHG1 regulates colorectal cancer cell growth through interactions with EZH2 and miR-154-5p. Mol Cancer. 2018;17(1):141.
- [Google Scholar]
- Molecular mechanisms of snoRNA-IL-15 crosstalk in adipocyte lipolysis and NK cell rejuvenation. Cell Metab. 2023;35(8):1457-1473.
- [Google Scholar]
- Identification of tumour immune infiltration-associated snoRNAs (TIIsno) for predicting prognosis and immune landscape in patients with colon cancer via a TIIsno score model. EBioMedicine. 2022;76
- [Google Scholar]
- The cancer genome atlas pan-cancer analysis project. Nat Genet. 2013;45(10):1113-1120.
- [Google Scholar]
- starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-seq data. Nucleic Acids Res. 2014;42:D92-D97.
- [Google Scholar]
- clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-287.
- [Google Scholar]
- Identification and validation of immune-related lncRNA prognostic signature for breast cancer. Genomics. 2020;112(3):2640-2646.
- [Google Scholar]
- Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453-457.
- [Google Scholar]
- Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
- [Google Scholar]
- Systematic analysis of lncRNA-miRNA-mRNA competing endogenous RNA network identifies four-lncRNA signature as a prognostic biomarker for breast cancer. J Transl Med. 2018;16(1):264.
- [Google Scholar]
- ceRNA cross-talk in cancer: when ce-bling rivalries go awry. Cancer Discov. 2013;3(10):1113-1121.
- [Google Scholar]
- Immune checkpoints in osteosarcoma: recent advances and therapeutic potential. Cancer Lett. 2022;547
- [Google Scholar]
- Circular RNA circ_001422 promotes the progression and metastasis of osteosarcoma via the miR-195-5p/FGF2/PI3K/Akt axis. J Exp Clin Cancer Res. 2021;40(1):235.
- [Google Scholar]
- Emerging next-generation sequencing-based discoveries for targeted osteosarcoma therapy. Cancer Lett. 2020;474:158-167.
- [Google Scholar]
- miR-34a predicts survival of Ewing's sarcoma patients and directly influences cell chemo-sensitivity and malignancy. J Pathol. 2012;226(5):796-805.
- [Google Scholar]
- Disorders and roles of tsRNA, snoRNA, snRNA and piRNA in cancer. J Med Genet. 2022;59(7):623-631.
- [Google Scholar]
- Regulatory role of small nucleolar RNAs in human diseases. BioMed Res Int. 2015;2015
- [Google Scholar]
- Targeting SNORA38B attenuates tumorigenesis and sensitizes immune checkpoint blockade in non-small cell lung cancer by remodeling the tumor microenvironment via regulation of GAB2/AKT/mTOR signaling pathway. J Immunother Cancer. 2022;10(5)
- [Google Scholar]
- Small nucleolar RNA and its potential role in breast cancer - a comprehensive review. Biochim Biophys Acta Rev Cancer. 2021;1875(1)
- [Google Scholar]
- The function of non-coding RNAs in lung cancer tumorigenesis. Cancers (Basel). 2019;11(5)
- [Google Scholar]
- Circulating pancreatic cancer exosomal RNAs for detection of pancreatic cancer. Mol Oncol. 2019;13(2):212-227.
- [Google Scholar]
- Plasma SNORD83A as a potential biomarker for early diagnosis of non-small-cell lung cancer. Future Oncol. 2022;18(7):821-832.
- [Google Scholar]
- Increased NFATC4 correlates with poor prognosis of AML through recruiting regulatory T cells. Front Genet. 2020;11
- [Google Scholar]
- Small nucleolar RNAs U32a, U33, and U35a are critical mediators of metabolic stress. Cell Metab. 2011;14(1):33-44.
- [Google Scholar]
- Dicer-independent snRNA/snoRNA-derived nuclear RNA 3 regulates tumor-associated macrophage function by epigenetically repressing inducible nitric oxide synthase transcription. Cancer Commun (Lond). 2021;41(2):140-153.
- [Google Scholar]
- Loss of SNORA73 reprograms cellular metabolism and protects against steatohepatitis. Nat Commun. 2021;12(1):5214.
- [Google Scholar]
- Long noncoding RNA ZFAS1 promoting small nucleolar RNA-mediated 2'-O-methylation via NOP58 recruitment in colorectal cancer. Mol Cancer. 2020;19(1):95.
- [Google Scholar]
- Cellular models of development of ovarian high-grade serous carcinoma: a review of cell of origin and mechanisms of carcinogenesis. Cell Prolif. 2021;54(5)
- [Google Scholar]
- Screening plasma exosomal RNAs as diagnostic markers for cervical cancer: an analysis of patients who underwent primary chemoradiotherapy. Biomolecules. 2021;11(11)
- [Google Scholar]
- Small RNA sequencing in the Tg4-42 mouse model suggests the involvement of snoRNAs in the etiology of alzheimer's disease. J Alzheimers Dis. 2022;87(4):1671-1681.
- [Google Scholar]
- Unique transcriptome changes in peripheral B cells revealed by comparing age groups from naive or vaccinated mice, including snoRNA and Cdkn2a. J Gerontol A Biol Sci Med Sci. 2020;75(12):2326-2332.
- [Google Scholar]
- SNORD11B-mediated 2'-O-methylation of primary let-7a in colorectal carcinogenesis. Oncogene. 2023;42(41):3035-3046.
- [Google Scholar]
- Interaction and cross-talk between non-coding RNAs. Cell Mol Life Sci. 2018;75(3):467-484.
- [Google Scholar]
- DLX2 promotes osteosarcoma epithelial-mesenchymal transition and doxorubicin resistance by enhancing HOXC8-CDH2 axis. iScience. 2023;26(11)
- [Google Scholar]
- Tissue microarray profiling and integrative proteomics indicate the modulatory potential of Maytenus royleanus in inhibition of overexpressed TPD52 in prostate cancers. Sci Rep. 2021;11(1)
- [Google Scholar]
- GRAMD1B regulates cell migration in breast cancer cells through JAK/STAT and akt signaling. Sci Rep. 2018;8(1):9511.
- [Google Scholar]
- Targeting molecular mechanisms underlying treatment efficacy and resistance in osteosarcoma: a review of current and future strategies. Int J Mol Sci. 2020;21(18)
- [Google Scholar]