The gene sequencing data and clinical information of 84 osteosarcoma patients were screened from the TARGET database (https://ocg.cancer.gov/programs/target). Gene expression data of musculoskeletal samples from 396 healthy humans were collected from the GTEx database (https://gtexportal.org/). To eliminate the platform data difference between TCGA and GTEx databases, the gene transcriptional expression data of each sample were transformed into log2 (FPKM value +1). Subsequently, the combat function from the “sva” R package was used to integrate the GTEx and TARGET datasets into one dataset.

The microarray dataset, named the GSE21257, was obtained from the high-throughput microarray expression profile database (Gene Expression Omnibus database, GEO, https://www.ncbi.nlm.nih.gov/geo/). GSE21257 contained the gene expression data and related clinical information of 53 osteosarcoma patients, which were used as the verification cohort for follow-up analysis and model verification.

The RNA-seq data was converted from FPKM to transcripts per million (TPM) using the algorithm described in the previous study for subsequent analysis (Li et al., 2010; Wagner et al., 2012). The 1,542 RBPs were obtained from previously published research (Gerstberger et al., 2014). The research process is shown in Figure 1A.

The “limma” package was used to screen the RBPs that were differentially expressed between osteosarcoma and normal tissues. The false discovery rate (FDR) < 0.05 and log2 fold change >1 were set as critical values. Then Gene Ontology (GO) function enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway analyses were performed to explore the potential molecular mechanisms of the DERBPs. Analyzed the interaction between the DERBPs by the STRING online tool (http://www.string-db.org/) (Udhaya Kumar et al., 2020). Cytoscape software (version 3.7.2) was used to construct and visualize the PPI network. The expression levels of ZC3HAV1 in several common cancers were visualized by using the TIMER online tool (https://cistrome.shinyapps.io/timer/) and GEPIA online tool (Gene Expression Profiling Interactive Analysis, http://gepia.cancer-pku.cn/) (Kumar et al., 2019).

Taking the TARGET dataset as the training cohort, the “survival” package in R was used to perform univariate Cox regression analysis to screen genes related to prognosis, which were then used as candidate genes for constructing the model. Based on the above-mentioned prognostic-related genes, we performed LASSO Cox regression analysis through the “glmnet” package to determine the best penalty value, and the corresponding gene was selected as the modeling gene. The formula for calculating the risk value of each patient was: R i s k   s c o r e = ∑ i = 1 n ( C o e f i ∗ x i )

Here, C o e f i   represents the RBPs coefficient, and x i represents the RBPs expression.

The above formula was used to calculate the risk score of each patient and divide them into high-risk groups and low-risk groups accordingly. Using the “survminer” package in R, we drawn the Kaplan–Meier (KM) curve. The “ggplot2” package was used to draw the distribution map of the risk score and survival status of the two patient groups and the modeled gene expression heatmap. The time-dependent ROC curve was produced through the “survivalROC” package in R. In parallel, combined with clinical information, such as age, gender, and tumor metastasis of osteosarcoma patients in the training cohort, the univariate and multivariate Cox regression model was used to analyze whether the risk score was an independent factor for judging the poor prognosis of osteosarcoma. The “rms” package in R was used to build and verify the nomogram model. The “calibrate” function of the rms software package was used to draw the calibration curve.

To verify the functional differences between the low-risk and high-risk groups, we used GSEA to compare the enrichment of tumor characteristic gene sets (hallmark gene sets) between the two groups. The tumor characteristic gene set comes from the Molecular Signatures database (MSigDB) database and contains 50 sets of genes involved in inflammation and hypoxia.

The osteosarcoma cell lines (U2OS, 143B) and the osteoblast hFOB 1.19 cell line were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The cells were grown in DMEM medium containing 10% fetal bovine serum (Gibco, United States) and 1% penicillin/streptomycin (solarbio, China). Osteosarcoma cells were grown at 37°C and 5% CO₂. Osteoblasts were cultured at 34°C with a volume fraction of 5% CO₂.

We collected six osteosarcoma tissues and six matched adjacent normal tissues. The samples came from patients who underwent surgery at The Second Affiliated Hospital of Nanchang University and were pathologically diagnosed with osteosarcoma. After the surgical resection, the sample is immediately stored in liquid nitrogen until the RNA or protein is extracted. All patients signed an informed consent form, and the study was approved by the Research Ethics Committee of the Second Affiliated Hospital of Nanchang University.

Using TRIzol™ Reagent (Thermo Fisher Scientific, United States) to extract total RNA from cells and tissues, and using PrimeScript™ RT reagent Kit (Perfect Real Time) reverse transcription kit (Takara, Japan) to reverse transcription into cDNA, GAPDH as internal control, and real-time PCR was performed. According to the 2^−ΔΔCt method for relative quantitative analysis. The primer sequence was shown in Supplementary Table S1.

RIPA lysis buffer (Beyotime, China) was used to extract proteins from cells and osteosarcoma tissues. The equal amount of protein was separated by 10% SDS-PAGE gel electrophoresis, and then transferred to PVDF membrane by electroblotting (BioRad, United States). Then it was blocked with 5% skim milk, and incubated overnight with anti-ZC3HAV1 (1:1,000, Proteintech, China) and anti-GAPDH (1:2000, Proteintech, China) primary antibodies at 4 °C. After washing with TBST, incubated with secondary antibody (1:2000, Proteintech, china) for 1 h at room temperature. The blot was observed using ECL (enhanced chemiluminescence) and analyzed using ImageJ software.

All statistical analyses were conducted using R. 4.0.4 (https://www.r-project.org/) and SPSS Statistics 25 (https://www.ibm.com/products/software). Two groups were compared with Student’s t-tests. p-value < 0.05 is considered to be statistically significant.

We extracted the RBPs expression matrix from the three datasets, and the Venn diagram displayed 1,093 intersected RBPs in all datasets (Figure 1A). Based on the dataset merged by TARGET and GTEx, we identified 142 DERBPs using the “limma” package and visualized by heatmap (Figure 1B). Next, functional analysis was performed on these DERBPs. GO analysis showed that, in terms of biological process (BP), DERBPs were mainly enriched in RNA splicing, regulation of translation, and mRNA metabolic process. In terms of cellular component (CC), cytoplasmic ribonucleoprotein granule, ribosome, and nuclear speck were enriched. In terms of molecular function (MF), DERBPs were mainly related to translation regulator activity, acting on RNA catalytic activity, and mRNA 3′-UTR binding (Figure 1C). Enrichment analysis of the KEGG pathway showed that the DERBPs mainly act on RNA transport (Figure 1D). We also constructed a PPI network. Figure 1E shows the interaction between the DERBPs. The darker the node, the more interacting proteins.

First, in the TARGET cohort, 142 DERBPs were analyzed using univariate regression analysis to obtain 11 RBPs related to prognosis (Figure 2A). LASSO regression analysis was performed on these 11 RBP genes (Figures 2B,C). When the log lambda reached the minimum, the parameters corresponding to the best modeling parameters and 10 model genes (TDRD6, TLR8, NXT2, EIF4E3, RPS27L, CPEB3, RBM34, TERT, RPS29, and ZC3HAV1), and their coefficients were obtained (Figure 2D).

We evaluated and verified the signature in the TARGET cohort and the GSE21257 cohort. Survival analysis showed a significant difference in survival between the high- and low-risk groups of the training and validation datasets (Figures 3A,B). For predictions of 1–5 years survival, the AUC values ​​in the training and test cohorts were 0.84, 0.87, and 0.88 and 0.72, 0.75, and 0.81, respectively (Figures 3C,D). As the risk score increased, the number of deaths gradually increased, and the survival time was significantly reduced (Figures 3E,F). PCA analysis showed that the distribution patterns of patients in different risk groups were significantly different (Figures 3G,H). These results show that the RBP signature has excellent forecasting capabilities for patient prognosis and has been verified in independent data.

In the TARGET cohort, univariate and multivariate Cox regression analysis showed that gender and age do not predict the prognosis of osteosarcoma patients, although metastasis status and risk score can be used as an independent factor in prognosis (Figures 4A,B). The difference in clinical characteristics and RBP gene expression patterns between the two groups is revealed by the heatmap (Figures 4C,D). The risk score was significantly related to metastasis, and high-risk patients were more likely to have tumor metastasis, as shown by the box plot in Figures 4E–J.

The expression levels of the 10 prognostic RBPs were displayed as boxplots, grouped according to the status of metastasis and risk (Figures 5A–D). Most RBPs showed obvious differential expression between the two groups, and the trends were consistent in the training and validation cohorts. We also established a clinical nomogram to clinically predict the survival of patients (Figure 5E). The C-index was used to evaluate the nomograms of the two datasets, and both sets of results showed robust predictive power (the C index of the TARGET dataset was 0.86, and the GSE21257 dataset was 0.74). A calibration plot showed that the predicted 3–5 years overall survival was in good agreement with the overall survival observed in the TARGET cohort (Figures 5F,G) and GSE21257 cohort (Supplementary Figure S1A,B). These data indicate that the nomogram is stable in predicting the survival of patients with osteosarcoma.

GSEA results showed that in the TARGET dataset, the signal pathways or biological processes related to the tumor and immunity were enriched in the high-risk group. Such as allograft rejection, complement, inflammatory response, IL6-JAK-STAT3 signaling, and the interferon-gamma response had higher negative enrichment scores (NES) in the high-risk groups (Figure 6).

We verified the expression levels of these prognostic RBPs genes, by using qRT-PCR and western blot. Our results showed that compared with osteoblasts, CPEB3, EIF4E3, RBM34, RPS27L, RPS29 and TDRD6 were down-regulated in U2OS and 143B, while NXT2, TERT, TLR8 and ZC3HAV1 were up-regulated in osteosarcoma cells (Figure 7A). The RT-qPCR results of the 10 prognostic RBPs genes expression levels are consistent with the RNA-sequence data (Figure 7B). Among these 10 RBPs, ZC3HAV1 is an antiviral protein, recent studies have shown that it is related to the occurrence of a variety of cancers (Lin et al., 2014; Todorova et al., 2015). It is up-regulated in breast cancer, cervical squamous cell carcinoma, cholangiocarcinoma, esophageal cancer, glioma, head and neck squamous cell carcinoma and other malignant tumors (Supplementary Figure S1C, D). However, the role of ZC3HAV1 in osteosarcoma has not been reported yet. Therefore, we further quantified the expression of ZC3HAV1 in normal osteoblasts (hFOB1.19 cells) and 143B and U2OS osteosarcoma cell lines. Western blotting showed that compared with normal osteoblasts, ZC3HAV1 protein levels in osteosarcoma cells were significantly up-regulated (Figure 7C). Consistently, ZC3HAV1 mRNA and protein levels in osteosarcoma tissues were significantly higher than adjacent normal tissues (Figures 7D,E).