Transcriptome data of 72 normal renal tissue specimens and 539 ccRCC specimens were downloaded from The Cancer Genome Atlas database (TCGA¹). We then used the edgeR package² to preprocess the raw data and identify the differentially expressed RBPs based on | log₂ fold change (FC)| > 1.0 and false discovery rate (FDR) < 0.05. We also downloaded the E-MTAB-1980 dataset from the ArrayExpress database³ and downloaded the transcriptome data of 436 ccRCC patients containing DFS information from the cBioportal database⁴.

We used the WEB-based Gene Set Analysis Toolkit (WebGestalt⁵) online analysis tool to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these differentially expressed RBPs (Liao et al., 2019). The GO terms including biological process, cellular component, and molecular function. All analysis results were screened according to the criteria of P < 0.05 and gene number > 5.

To identify RBPs with important prognostic significance, we first performed univariate Cox regression analysis of all these differentially expressed RBPs. The least absolute shrinkage and selection operator (LASSO) regression analysis was then used for further screening. Finally, multivariate Cox regression analysis was used to further screen out RBPs with important prognostic value. A P < 0.05 was considered significant.

We constructed a multivariate Cox proportional hazards regression model to predict the prognosis of ccRCC patients based on these prognostic related RBPs. The risk score for each patient in the model was calculated using the following formula:

In this formula, Exp represents the expression value of each gene, and β represents the corresponding regression coefficient. We then divided ccRCC patients from the TCGA cohort into low-risk and high-risk subgroups based on the median risk score, and compared OS between the two groups to initially assess the predictive power of the model. In addition, we used the Survival ROC R package to establish the ROC curve to assess the prognostic efficacy of the model and used the rms R package to draw the nomogram to predict OS. Finally, we divided the 539 samples in the TCGA cohort into the training group and the validation group as internal validation and the E-MTAB-1980 cohort with 101 sample information as external validation to evaluate the stability and predictive efficacy of the model.

To explore the clinical significance of the prognostic model in different clinical parameters, we stratified the patients according to the different clinical parameters and performed survival analysis. We also explored the relationship between these eight prognostic RBPs and clinical parameters. A P < 0.05 was considered significant.

We divided the patients into low-risk and high-risk groups based on the median risk score of the prognostic model, and then performed gene set enrichment analysis (GSEA) by using GSEA_4.0.3 software⁶. A P < 0.05 and FDR < 0.25 were considered to be significant differences.

We used The Human Protein Atlas (HPA⁷) online database to verify the protein expression levels of these prognostic related RBPs. And the Kaplan–Meier plotter⁸ online tool was used to assess the prognostic significance of these prognostic related RBPs in ccRCC patients.

Since DFS is also important for the prognosis of tumor patients, we constructed a prognostic model for DFS. We downloaded transcriptomic data from the cBioportal database for 436 ccRCC patients with DFS information. Then the prognostic RBPs were screened by Cox regression analysis and LASSO regression analysis and a prognostic model for DFS was constructed.

R software (Version 4.0.0) was used for statistical analysis. The differentially expressed genes in tumor tissues and normal tissues were analyzed by “edgeR” package. Cox regression analysis was used to screen for genes associated with prognosis. The OS and DFS of patients were analyzed by Kaplan–Meier method and log-rank test. The “survival ROC” package was used to analyze the ROC curve. The “rms” package was used to draw the nomogram. The Student’s t-test or non-parametric Mann–Whitney rank sum test was used to compare the correlation between risk score, prognostic genes, and clinicopathological variables. P < 0.05 was considered statistically significant.

The analysis process of this study was shown in Figure 1. Transcriptome data of ccRCC patients were downloaded from the TCGA database, including 72 normal renal tissue samples and 539 tumor tissue samples (Supplementary Table S1). The edger R package was used to process the data and identify the differentially expressed RBPs. Of the 1542 RBPs (Gerstberger et al., 2014), 40 met our criteria (| log2 FC| > 1.0, FDR < 0.05), including 10 down-regulated RBPs and 30 up-regulated RBPs. Figure 2 showed the expression and distribution of these differentially expressed RBPs.

We performed GO and KEGG enrichment analyses for these differentially expressed RBPs using the WebGestalt online analysis tool to investigate the biological functions and molecular mechanisms of these genes. The analysis results were shown in Table 1. The biological processes analysis showed that these RBPs were significantly enriched in RNA catabolic process, posttranscriptional regulation of gene expression, translational initiation, regulation of cellular amide metabolic process, protein localization to endoplasmic reticulum, meiotic cell cycle, gene silencing, cellular process involved in reproduction in multicellular organism, transposition, and regulation of mRNA metabolic process. The cellular component showed that these RBPs were significantly enriched in polysome, ribosome, ribonucleoprotein granule, cytosolic part, and rough endoplasmic reticulum. In terms of molecular function, these RBPs were significantly enriched in mRNA binding, catalytic activity, acting on RNA, structural constituent of ribosome, helicase activity, nuclease activity, translation regulator activity, snRNA binding, double-stranded RNA binding, ATPase activity, and nucleotidyltransferase activity. Moreover, KEGG analysis showed that these RBPs were mainly enriched in ribosome, RNA transport, influenza A, mRNA surveillance pathway, herpes simplex infection, and ribosome biogenesis in eukaryotes.

We performed a univariate Cox regression analysis on all these differentially expressed RBPs and obtained 25 prognostic related RBPs (Figure 3). We further performed LASSO regression analysis on these 25 genes to screen the RBPs with prognostic significance, and obtained 9 RBPs including APOBEC3G, AUH, DAZL, DDX47, EIF4A1, IGF2BP3, NR0B1, RPL36A, and TRMT1 (Supplementary Figure S2). And multivariate Cox regression analysis showed that 8 of the 9 RBPs, namely, APOBEC3G, AUH, DAZL, EIF4A1, IGF2BP3, NR0B1, RPL36A, and TRMT1 independently predicted prognosis of ccRCC patients.

We used these eight genes screened from multivariate Cox regression analysis to establish a prognostic model for OS (Table 2). Each ccRCC patient’s risk score was calculated according to the following formula:

Risk score = (0.0951 × Exp APOBEC3G) + (−0.1621 × Exp AUH) + (0.0945 × Exp DAZL) + (0.1571 × Exp EIF4A1) + (0.1190 × Exp IGF2BP3) + (0.0998 × Exp NR0B1) + (0.1722 × Exp RPL36A) + (0.2380 × Exp TRMT1) Based on the median risk score, 539 ccRCC patients in the TCGA cohort were divided into low-risk and high-risk subgroups for survival analysis to assess the predictive power of the model. Survival analysis showed that patients in the high-risk group had lower OS than those in the low-risk group (P = 9.556e-13, Figure 4A). We then performed the time-dependent receiver operating characteristic (ROC) analysis to further evaluate the predictive performance of the eight RBPs signature, and the area under the ROC curve (AUC) of the model was 0.729 at 1 year, 0.702 at 3 years, and 0.726 at 5 years (Figure 4B). Figure 4C showed the survival status of each patient in the TCGA cohort assessed by risk score. Subsequently, to evaluate the applicability and stability of the prognostic model for OS, these 539 ccRCC patients in the TCGA cohort were randomly divided into a training data set and a validation data set. We then used the same formula to calculate the risk score of each patient to assess the predictive performance of the model. The results showed that patients in the high-risk group in the training data set had worse OS than those in the low-risk group (P = 1.908e-05, Figure 5A). We found that the AUC was 0.750 at 1 year, 0.697 at 3 years, and 0.759 at 5 years (Figure 5B). And patients in the validation data set had similar results (Figures 5C,D). In addition, to assess whether the model has similar predictive power in other ccRCC patient cohorts, the same risk score formula was used for the E-MTAB-1980 dataset. Survival analysis also showed that patients in the high-risk group had lower OS than those in the low-risk group (P = 0.00033, Figure 6A), and the AUC of the model was 0.788 at 1 year, 0.783 at 3 years, and 0.795 at 5 years (Figure 6B). And Figure 6C showed the survival status of each patient in the E-MTAB-1980 cohort assessed by risk score. These results showed that the signature of these eight RBPs has good predictive performance and stability.

Next, we used the risk score in the prognostic model as the label to explore the functional differences between the two subgroups by conducting GSEA. The results showed that ribosome was mainly enriched in the high-risk group (Figure 7A), indicating that the regulation of RBPs was mainly involved in high-risk ccRCC patients. In addition, we performed a univariate Cox regression analysis for different clinical characteristics of ccRCC patients to evaluate their respective predictive significance. The results showed that age, tumor grade, tumor stage, primary tumor location, regional lymph node invasion, distant metastasis, and risk score were all associated with the OS of ccRCC patients (Figure 7B). However, multiple regression analysis showed that only age (P < 0.001), tumor grade (P = 0.020), tumor stage (P < 0.001), and risk score (P < 0.001) were independent prognostic factors related to OS of ccRCC patients (Figure 7C). Moreover, to establish a quantitative prognostic approach for ccRCC patients, we drew a nomogram based on the risk score and other clinical variables (Figure 7D). By drawing a vertical line between each prognosis axis and the total point axis, we can predict the survival probability of ccRCC patients at 1, 3, and 5 years. We also constructed calibration curves to evaluate the predictive performance of the nomogram, and the results showed that there was high consistency between the predicted results and the actual results (Figures 7E–G). And we used the TCGA and E-MTAB-1980 cohorts to verify the accuracy and stability of nomogram to expand its clinical application and availability. Survival analysis showed that nomogram could better distinguish ccRCC patients with low survival rates in TCGA and E-MTAB-1980 cohorts (P < 0.001 and P = 2.32e-05, Figures 7H,J). Based on the nomogram, the AUC in the TCGA cohort was 0.867 at 1 year, 0.806 at 3 years and 0.778 at 5 years (Figure 7I), and the AUC in the E-MTAB-1980 cohort was 0.910 at 1 year, 0.917 at 3 years, and 0.892 at 5 years (Figure 7K).

To explore the clinical significance of the signature based on these eight RBPs in the ccRCC patients stratified by different clinical parameters, we stratified ccRCC patients from TCGA database according to age, gender, grade, stage, T stage, M stage, and N stage. Kaplan–Meier survival curve analysis showed that the OS was significantly shorter for the ccRCC patients in the high-risk group compared to the low-risk group ccRCC patients (Figure 8). These results indicate that the signature of these eight RBPs can predict the prognosis of ccRCC patients without considering clinical parameters.

We analyzed the correlation between the prognostic model based on these eight RBPs and clinical parameters to explore whether the prognostic model might influence the progression of ccRCC. The results showed no significant correlation between age and prognostic model (Figure 9A). However, the risk score of females was significantly lower than that of male (Figure 9B), the risk score of G1-2 was significantly lower than that of G3-4 (Figure 9C), the risk score of stage I-II was significantly lower than that of stage III-IV (Figure 9D), the risk score of T1-2 was significantly lower than that of T3-4 (Figure 9E), the risk score of M0 was significantly lower than that of M1 (Figure 9F) (The N1 in the N stage is very small and cannot be analyzed). These results showed that prognostic model for OS was significantly correlated with ccRCC tumor progression.

Based on the above results, we also analyzed the relationship between prognostic RBPs for OS and clinical parameters to further investigate the role of prognostic RBPs in ccRCC. The results showed that AUH, EIF4A1, IGF2BP3, and RPL36A were significantly correlated with gender; APOBEC3G, AUH, IGF2BP3, RPL36A, and TRMT1 were significantly correlated with grade; APOBEC3G, AUH, DAZL, IGF2BP3, RPL36A, and TRMT1 were significantly correlated with stage; APOBEC3G, AUH, DAZL, IGF2BP3, NR0B1, RPL36A, and TRMT1 were significantly correlated with T stage; APOBEC3G, AUH, IGF2BP3, RPL36A, and TRMT1 were significantly correlated with M stage. However, there was no significant correlation between NR0B1 and these clinical parameters (Table 3).

To assess the prognostic significance of these prognostic related RBPs in ccRCC patients, we used the Kaplan–Meier plotter online tool to confirm the relationship between these genes and OS. The results showed that all the eight RBPs were related to the OS in ccRCC patients (Figure 10). Subsequently, we used the HPA online database to verify the protein expression levels of these prognostic related RBPs, the results showed that APOBEC3G, EIF4A1, and TRMT1 were significantly increased in ccRCC tissue compared with normal renal tissue (Figures 11A,D,G). And AUH, DAZL, IGF2BP3, and RPL36A were significantly reduced in ccRCC tissue compared with normal renal tissue (Figures 11B,C,E,F). However, the protein expression level of NR0B1 was not available on the HPA online database.

In view of the important influence of DFS on the prognosis of ccRCC, we also constructed a prognostic model for DFS. The expression data of 436 ccRCC patients and the corresponding DFS information were download from the cBioportal database. We then identified 9 prognostic RBPs including APOBEC3G, AUH, DDX47, IGF2BP3, MOV10L1, NANOS1, PIH1D3, TDRD9, and TRMT1 by univariate Cox regression analysis, LASSO regression analysis and multivariate Cox regression analysis. We then constructed a prognostic model for DFS based on these nine prognostic RBPs and calculated each patient’s risk score based on the following formula: Risk score = (0.0852 × Exp APOBEC3G) + (−0.3683 × Exp AUH) + (0.4195 × Exp DDX47) + (0.1445 × Exp IGF2BP3) + (−0.2077 × Exp MOV10L1) + (0.4206 × Exp NANOS1) + (0.7675 × Exp PIH1D3) + (−0.1011 × Exp TDRD9) + (0.2895 × Exp TRMT1). Based on the median risk score, these 436 ccRCC patients were divided into high-risk and low-risk groups for survival analysis to assess the predictive performance of the prognostic model. The results showed that patients in the high-risk group had worse DFS than those in the low-risk group (P = 1.110e-16, Figure 12A). We found that the AUC for DFS was 0.729 at 1 year, 0.764 at 3 years, and 0.782 at 5 years (Figure 12D). These results showed that the RBPs associated prognostic model for DFS has good predictive performance.

In addition, we randomly divided the whole dataset into a training data set (n = 218) and a validation data set (n = 218) to assess the applicability and stability of the prognostic model for DFS. We used the same formula to calculate each patient’s risk score. Survival analysis showed that patients in the high-risk group in the training data set had worse DFS than those in the low-risk group (P = 1.127e-10, Figure 12B). The AUC for DFS was 0.718 at 1 year, 0.763 at 3 years, and 0.813 at 5 years (Figure 12E). Patients in the validation data set had similar results (Figures 12C,F).

Moreover, the prognostic value of the prognostic model for DFS and different clinical parameters were evaluated by Cox regression analysis. The results indicated that the tumor grade, tumor stage, primary tumor location, distant metastasis, and risk score of ccRCC patients were significantly correlated with DFS (P < 0.001, Figure 12G). However, multiple regression analysis revealed that tumor grade, tumor stage, and risk score were independent prognostic factors associated with DFS (P < 0.001, Figure 12H). These results suggested that the RBPs associated prognostic model for DFS was also a good predictor of ccRCC patient outcomes.