Abstract RNA binding proteins (RBPs) play a key role in post-transcriptional gene regulation. They have been shown to be dysfunctional in a variety of cancers and are closely related to the occurrence and progression of cancers. However, the biological function and clinical significance of RBPs in clear cell renal carcinoma (ccRCC) are unclear. In our current study, we downloaded the transcriptome data of ccRCC patients from The Cancer Genome Atlas (TCGA) database and identified differential expression of RBPs between tumor tissue and normal kidney tissue. Then the biological function and clinical value of these RBPs were explored by using a variety of bioinformatics techniques. We identified a total of 40 differentially expressed RBPs, including 10 down-regulated RBPs and 30 up-regulated RBPs. Eight RBPs (APOBEC3G, AUH, DAZL, EIF4A1, IGF2BP3, NR0B1, RPL36A, and TRMT1) and nine RBPs (APOBEC3G, AUH, DDX47, IGF2BP3, MOV10L1, NANOS1, PIH1D3, TDRD9, and TRMT1) were identified as prognostic related to overall survival (OS) and disease-free survival (DFS), respectively, and prognostic models for OS and DFS were constructed based on these RBPs. Further analysis showed that OS and DFS were worse in high-risk group than in the low-risk group. The area under the receiver operator characteristic curve of the model for OS was 0.702 at 3 years and 0.726 at 5 years in TCGA cohort and 0.783 at 3 years and 0.795 at 5 years in E-MTAB-1980 cohort, showing good predictive performance. Both models have been shown to independently predict the prognosis of ccRCC patients. We also established a nomogram based on these prognostic RBPs for OS and performed internal validation in the TCGA cohort, showing an accurate prediction of ccRCC prognosis. Stratified analysis showed a significant correlation between the prognostic model for OS and ccRCC progression. Keywords: clear cell renal cell carcinoma, RNA binding proteins, prognostic model, survival analysis, bioinformatics Introduction Renal cell carcinoma (RCC) accounts for 2.4% of all malignancies, with an estimated 400,000 new cases and 175,000 deaths worldwide each year ([45]Bray et al., 2018; [46]Siegel et al., 2018). the clear cell renal cell carcinoma (ccRCC) is the most common subtype of RCC, accounting for approximately 70–80% and presents a high risk of heterogeneity and metastasis ([47]Rini et al., 2009; [48]Ljungberg et al., 2019). Although surgical resection can effectively resolve the early stage of ccRCC, 30% of patients still have recurrence or metastasis after surgery, and the late stage of ccRCC has a high mortality rate due to insensitivity to traditional radiotherapy or chemotherapy ([49]Battaglia and Lucarelli, 2015; [50]Moch et al., 2016; [51]Tamma et al., 2019). Therefore, further understanding of the molecular mechanisms of ccRCC and the discovery of more effective molecular biomarkers are essential for early screening, diagnosis, monitoring for metastasis, recurrence, and quality of life in patients. Post-transcriptional regulation of RNA is an important aspect of gene expression regulation. RNA binding proteins (RBPs) are a class of proteins widely expressed in cells, which form ribonucleoprotein (RNP) complexes through binding at different sites or random interaction with target RNA, thus strictly regulating RNA metabolism ([52]Iadevaia and Gerber, 2015; [53]Hentze et al., 2018). Currently, there are 1,542 RBP coding genes, accounting for 7.5% of all human protein-coding genes, which have been verified by experiments ([54]Gerstberger et al., 2014). These RBPs regulate a variety of biological processes including RNA processing, splicing, mRNA stability, output, localization, and translation, thus maintaining the physiological balance of the cell ([55]Masuda and Kuwano, 2019). Given this, it comes as no surprise that RBPs dysfunction has been linked to a variety of human diseases. Ribosomal diseases caused by ribosomal protein and rRNA biogenic factor defects, such as Diamond–Blackfan anemia and Shwachman–Diamond syndrome, affect the same tissues and exhibit similar pathology precisely because RBPs bind to the same type of RNA ([56]Narla and Ebert, 2010). Mutations in mRBPs or their targets in neurons lead to abnormal aggregation of proteins or RNA, resulting in a variety of neurodegenerative and neuromuscular diseases ([57]Scheper et al., 2007). However, the role of RBPs in tumor genesis and development is rare. Some studies have shown that RBPs are abnormally expressed in tumor tissues compared with normal tissues and are associated with patient prognosis ([58]Patry et al., 2003; [59]Busà et al., 2007; [60]Ortiz-Zapater et al., 2011). In lung cancer, QKI inhibits tumor cell proliferation by competing with the splicing factor SF1 ([61]Zong et al., 2014). In melanoma, CPEB4 promotes tumor cell proliferation by regulating polyadenylation and promoting the translation of melanoma drivers ([62]Pérez-Guijarro et al., 2016). Knockdown SAM68 in breast cancer cells inhibited tumor cell proliferation by upregulation of cell cycle inhibitors P21 and CDKN1B/P27 ([63]Song et al., 2010). However, in the field of ccRCC, existing studies only described the effect of RBPs on the overall survival (OS) of ccRCC patients ([64]Hua et al., 2020; [65]Zhu et al., 2020), and few RBPs models can be used to predict the prognosis of ccRCC patients. The development of new RBPs models has gradually become an effective method to explore new therapeutic targets. Therefore, in our current study, we systematically and deeply analyzed the molecular biological function and clinical significance of RBPs in ccRCC to promote our understanding of ccRCC progress, and established risk score models for OS and disease-free survival (DFS), which may provide new biomarkers for disease diagnosis and treatment prognosis. Materials and Methods Preprocessing Data and Identifying Differential Expression RBPs Transcriptome data of 72 normal renal tissue specimens and 539 ccRCC specimens were downloaded from The Cancer Genome Atlas database (TCGA^[66]1). We then used the edgeR package^[67]2 to preprocess the raw data and identify the differentially expressed RBPs based on | log[2] fold change (FC)| > 1.0 and false discovery rate (FDR) < 0.05. We also downloaded the E-MTAB-1980 dataset from the ArrayExpress database^[68]3 and downloaded the transcriptome data of 436 ccRCC patients containing DFS information from the cBioportal database^[69]4. Function and Pathway Enrichment Analysis We used the WEB-based Gene Set Analysis Toolkit (WebGestalt^[70]5) online analysis tool to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these differentially expressed RBPs ([71]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. Selection of Prognostic Related RBPs 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. Construction and Evaluation of Prognostic Model for OS 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: [MATH: Riskscor< /mi>e=∑i=1n Expiβi, :MATH] 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. Correlation Between Prognostic Model for OS, Prognostic RBPs and Clinical Parameters 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. Gene Set Enrichment Analysis 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^[72]6. A P < 0.05 and FDR < 0.25 were considered to be significant differences. Express Level and Prognostic Significance Verification of Prognostic Related RBPs We used The Human Protein Atlas (HPA^[73]7) online database to verify the protein expression levels of these prognostic related RBPs. And the Kaplan–Meier plotter^[74]8 online tool was used to assess the prognostic significance of these prognostic related RBPs in ccRCC patients. Construction and Evaluation of Prognostic Model for DFS 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. Statistical Analysis 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. Results Screening Differentially Expressed RBPs in ccRCC The analysis process of this study was shown in [75]Figure 1. Transcriptome data of ccRCC patients were downloaded from the TCGA database, including 72 normal renal tissue samples and 539 tumor tissue samples ([76]Supplementary Table S1). The edger R package was used to process the data and identify the differentially expressed RBPs. Of the 1542 RBPs ([77]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. [78]Figure 2 showed the expression and distribution of these differentially expressed RBPs. FIGURE 1. [79]FIGURE 1 [80]Open in a new tab The research diagram for analyzing RBPs in ccRCC. FIGURE 2. [81]FIGURE 2 [82]Open in a new tab The expression and distribution of differentially expressed RBPs in ccRCC. (A) Heatmap of the 40 differentially expressed RBPs; (B) Volcano plot of 455 RBPs (C) visualization of the expression levels of the 40 differentially expressed RBPs. Function and Pathway Enrichment Analysis 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 [83]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. TABLE 1. KEGG pathway and GO enrichment analysis of differentially expressed RNA binding proteins. GO term P-value Biological processes RNA catabolic process 2.81e-12 Posttranscriptional regulation of gene expression 2.34e-10 Translational initiation 3.72e-8 Regulation of cellular amide metabolic process 5.31e-8 Protein localization to endoplasmic reticulum 0.000002 Meiotic cell cycle 0.000004 Gene silencing 0.000016 Cellular process involved in reproduction in multicellular organism 0.000034 Transposition 0.000057 Regulation of mRNA metabolic process 0.000073 Cellular component Polysome 3.38e-10 Ribosome 4.88e-8 Ribonucleoprotein granule 6.13e-7 Cytosolic part 0.000001 Rough endoplasmic reticulum 0.019853 Molecular function mRNA binding 4.25e-10 Catalytic activity, acting on RNA 2.15e-8 Structural constituent of ribosome 1.84e-7 Helicase activity 0.000052 Nuclease activity 0.000228 Translation regulator activity 0.000561 snRNA binding 0.005137 Double-stranded RNA binding 0.018541 ATPase activity 0.031468 Nucleotidyltransferase activity 0.049070 KEGG pathway Ribosome 1.03e-8 RNA transport 0.000494 Influenza A 0.000553 mRNA surveillance pathway 0.001136 Herpes simplex infection 0.008472 Ribosome biogenesis in eukaryotes 0.014958 [84]Open in a new tab GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Prognostic Related RBPs Selection We performed a univariate Cox regression analysis on all these differentially expressed RBPs and obtained 25 prognostic related RBPs ([85]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 ([86]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. FIGURE 3. [87]FIGURE 3 [88]Open in a new tab Univariate Cox regression analysis of differentially expressed RBPs. Prognostic Related Risk Score Model for OS Construction and Evaluation We used these eight genes screened from multivariate Cox regression analysis to establish a prognostic model for OS ([89]Table 2). Each ccRCC patient’s risk score was calculated according to the following formula: TABLE 2. Multivariate Cox regression analysis to identify prognosis-related RNA binding proteins. Gene Coef Exp(coef) se(coef) z Pr (> | z|) APOBEC3G 0.0951 1.0998 0.0844 1.1264 0.2600 AUH −0.1621 0.8504 0.1318 −1.2299 0.2187 DAZL 0.0945 1.0991 0.0695 1.3597 0.1739 EIF4A1 0.1571 1.1701 0.0710 2.2138 0.0268 IGF2BP3 0.1190 1.1264 0.0346 3.4376 0.0006 NR0B1 0.0998 1.1050 0.0366 2.7241 0.0064 RPL36A 0.1722 1.1879 0.1493 1.1532 0.2488 TRMT1 0.2380 1.2687 0.1632 1.4583 0.1448 [90]Open in a new tab Coef, coefficient. 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, [91]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 ([92]Figure 4B). [93]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, [94]Figure 5A). We found that the AUC was 0.750 at 1 year, 0.697 at 3 years, and 0.759 at 5 years ([95]Figure 5B). And patients in the validation data set had similar results ([96]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, [97]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 ([98]Figure 6B). And [99]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. FIGURE 4. [100]FIGURE 4 [101]Open in a new tab Risk score analysis of the eight RBPs prognostic model for OS in the TCGA cohort. (A) Kaplan–Meier OS curve analysis; (B) time dependent ROC curve analysis; (C) survival status of each patient. FIGURE 5. [102]FIGURE 5 [103]Open in a new tab Risk score analysis of the eight RBPs prognostic model for OS in the training and validation data set. (A) Kaplan–Meier OS curve analysis in the training data set; (B) time dependent ROC curve analysis in the training data set; (C) Kaplan–Meier OS curve analysis in the validation data set; (D) time dependent ROC curve analysis in the validation data set. FIGURE 6. [104]FIGURE 6 [105]Open in a new tab Risk score analysis of the eight RBPs prognostic model for OS in the E-MTAB-1980 cohort. (A) Kaplan–Meier OS curve analysis; (B) time dependent ROC curve analysis; (C) survival status of each patient. 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 ([106]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 ([107]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 ([108]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 ([109]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 ([110]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, [111]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 ([112]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 ([113]Figure 7K). FIGURE 7. [114]FIGURE 7 [115]Open in a new tab Construction of a nomogram and assessment the prognostic significance of different clinical characteristics in ccRCC patients. (A) Gene set enrichment analysis comparing the high-risk and low-risk groups based on the TCGA cohort; (B) univariate Cox regression analysis of correlations between risk score for OS and clinical parameters; (C) multivariate Cox regression analysis of correlations between risk score for OS and clinical parameters; (D) nomogram for predicting the 1- year, 3-year, and 5-year OS of ccRCC patients; (E–G) calibration curves of the nomogram to predict OS at 1, 3, and 5 years; (H) Kaplan–Meier OS curve analysis in the TCGA cohort based the nomogram; (I) time dependent ROC curve analysis in the TCGA cohort based the nomogram; (J) Kaplan–Meier OS curve analysis in the E-MTAB-1980 cohort based the nomogram; (K) Time dependent ROC curve analysis in the E-MTAB-1980 cohort based the nomogram. Prognostic Value of the Prognostic Model for OS Stratified by Clinical Parameters 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 ([116]Figure 8). These results indicate that the signature of these eight RBPs can predict the prognosis of ccRCC patients without considering clinical parameters. FIGURE 8. [117]FIGURE 8 [118]Open in a new tab Kaplan–Meier survival curves analysis stratified by different clinical parameters. (A) Age ≤ 65; (B) Age > 65; (C) Female; (D) Male; (E) Grade1-2; (F) Grade3-4; (G) Stage I-II; (H) Stage III-IV; (I) T stage1-2; (J) T stage3-4; (K) M stage0; (L) M stage1; (M) N stage0; (N) N stage1-X. Relationship Between Prognostic Model for OS and 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 ([119]Figure 9A). However, the risk score of females was significantly lower than that of male ([120]Figure 9B), the risk score of G1-2 was significantly lower than that of G3-4 ([121]Figure 9C), the risk score of stage I-II was significantly lower than that of stage III-IV ([122]Figure 9D), the risk score of T1-2 was significantly lower than that of T3-4 ([123]Figure 9E), the risk score of M0 was significantly lower than that of M1 ([124]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. FIGURE 9. [125]FIGURE 9 [126]Open in a new tab Relationship between prognostic model for OS and clinical parameters. (A) Age; (B) Gender; (C) Grade; (D) Stage; (E) T stage; (F) M stage. Correlation Between Prognostic RBPs and Clinical Parameters 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 ([127]Table 3). TABLE 3. The relationship between prognostic related RNA binding proteins and clinicopathologic parameters. Gene Gender __________________________________________________________________ Grade __________________________________________________________________ Stage __________________________________________________________________ T stage __________________________________________________________________ M stage __________________________________________________________________ Female Male G1-2 G3-4 I-II III-IV T1-T2 T3-T4 M0 M1 N 186 353 249 282 331 205 349 190 428 78 APOBEC3G t-value 1.432 5.900 5.688 5.095 4.057 P-value 0.153 <0.001 <0.001 <0.001 <0.001 AUH t-value 2.799 NA* 6.545 5.595 4.589 P-value 0.005 <0.001 <0.001 <0.001 <0.001 DAZL t-value 0.181 NA* NA* NA* 0.817 P-value 0.857 0.256 0.047 0.049 0.415 EIF4A1 t-value 2.652 0.947 1.545 1.783 0.845 P-value 0.008 0.344 0.123 0.075 0.398 IGF2BP3 t-value 2.566 6.141 NA* NA* NA* P-value 0.011 <0.001 <0.001 <0.001 <0.001 NR0B1 t-value 0.951 0.355 NA* NA* NA* P-value 0.342 0.723 0.087 0.030 0.569 RPL36A t-value NA* 3.359 4.685 3.871 2.464 P-value <0.001 0.001 <0.001 <0.001 0.014 TRMT1 t-value 0.526 3.356 2.443 2.059 2.225 P-value 0.599 <0.001 0.015 0.040 0.027 [128]Open in a new tab NA, not available. *Non-parametric Mann–Whitney rank sum test. Express Level and Prognostic Significance Verification of Prognostic Related RBPs 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 ([129]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 ([130]Figures 11A,D,G). And AUH, DAZL, IGF2BP3, and RPL36A were significantly reduced in ccRCC tissue compared with normal renal tissue ([131]Figures 11B,C,E,F). However, the protein expression level of NR0B1 was not available on the HPA online database. FIGURE 10. [132]FIGURE 10 [133]Open in a new tab Validation the prognostic value of the prognostic RBPs for OS in ccRCC by Kaplan–Meier plotter. (A) APOBEC3G; (B) AUH; (C) DAZL; (D) EIF4A1; (E) IGF2BP3; (F) NR0B1; (G) RPL36A; (H) TRMT1. FIGURE 11. [134]FIGURE 11 [135]Open in a new tab The expression status of the prognostic RBPs proteins in ccRCC and normal renal tissues in the HPA database. (A) APOBEC3G; (B) AUH; (C) DAZL; (D) EIF4A1; (E) IGF2BP3; (F) RPL36A; (G) TRMT1. Construction of a Prognostic Model for DFS 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, [136]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 ([137]Figure 12D). These results showed that the RBPs associated prognostic model for DFS has good predictive performance. FIGURE 12. [138]FIGURE 12 [139]Open in a new tab The RBPs-related prognostic model for DFS. (A) Kaplan–Meier DFS curve analysis in the cBioportal cohort; (B) Kaplan–Meier DFS curve analysis in the training data set; (C) Kaplan–Meier DFS curve analysis in the validation data set (D) time dependent ROC curve analysis in the cBioportal cohort; (E) time dependent ROC curve analysis in the training data set; (F) time dependent ROC curve analysis in the validation data set; (G) univariate Cox regression analysis of correlations between risk score for DFS and clinical parameters; (H) multivariate Cox regression analysis of correlations between risk score for DFS and clinical parameters. 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, [140]Figure 12B). The AUC for DFS was 0.718 at 1 year, 0.763 at 3 years, and 0.813 at 5 years ([141]Figure 12E). Patients in the validation data set had similar results ([142]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, [143]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, [144]Figure 12H). These results suggested that the RBPs associated prognostic model for DFS was also a good predictor of ccRCC patient outcomes. Discussion Malignant tumor is a kind of complex heterogeneous diseases, apart from the classic view that affect cancer or tumor suppressor gene signal channel change decision, It has also been found to be associated with post-transcriptional hijacking by tumor cells, enabling them to rapidly and stably regulate protein expression levels in response to intracellular and extracellular signaling changes to adapt to local microenvironments ([145]Pereira et al., 2017). RBPs are a key player in post-transcriptional events, participating in almost all post-transcriptional regulation, controlling intracellular transcript metabolism and function, and maintaining homeostasis. Multiple studies have reported that RBPs are dysregulated in cancers and regulate cancers progression through a variety of mechanisms, including transcriptional and posttranscriptional regulation, genomic change, and posttranslational modification ([146]Patry et al., 2003; [147]Busà et al., 2007; [148]Ortiz-Zapater et al., 2011; [149]Song et al., 2010; [150]Zong et al., 2014; [151]Pérez-Guijarro et al., 2016). However, the expression pattern and role of RBPs in ccRCC are rarely reported. In this study, we systematically analyzed the transcriptome data of ccRCC patients from TCGA database, and identified differential expression RBPs between tumor tissue and normal kidney tissue. We then performed functional enrichment analysis to evaluate their biological function, and performed univariate Cox regression analysis, LASSO regression analysis and multivariate Cox regression analysis to screen prognostic related RBPs and constructed a prognostic risk score model for OS based on these RBPs. In addition, we also built a prognostic model for DFS to predict ccRCC prognosis based on prognostic related RBPs. The biological function and pathway enrichment analysis of differentially expressed RBPs showed that these genes were significantly enriched in posttranscriptional regulation of gene expression, translational initiation, transposition, protein localization to endoplasmic reticulum, RNA catabolic process, regulation of cellular amide metabolic process, regulation of mRNA metabolic process, gene silencing, ribonucleoprotein granule, mRNA binding, ribosome, polysome, catalytic activity, acting on RNA, translation regulator activity, nuclease activity, double-stranded RNA binding, nucleotidyltransferase activity, RNA transport, and mRNA surveillance pathway, which involved RNA processing, splicing, localization, RNA metabolism and subsequent translation regulation. Previous studies have shown that multiple RBPs regulatory mechanisms have been identified in cancers, including transcriptional and posttranscriptional regulation, genomic change, and posttranslational modification ([152]Patry et al., 2003; [153]Trabucchi et al., 2009; [154]Ishii et al., 2014; [155]Preca et al., 2015). In lung adenocarcinoma, splice regulator RBM10 inhibits tumor cell proliferation and Notch signaling activity ([156]Bechara et al., 2013). Cancer transcription factor MYC up-regulates the mRNA expression of hnRNPA1 and hnRNPA2 in gliomas, which promotes the synthesis of pyruvate kinase M subtype 2 (PKM2) and participates in glycolytic transformation ([157]Clower et al., 2010). IMP1 has been reported to be elevated in multiple tumors, and reduced IMP1 expression can impair the normal transmission and local translation of adhesive and motif-related target mRNAs ([158]Gu et al., 2009). EIF4E is a key factor in mRNA cycling and translation, and it has been found that EIF4E is overexpressed in a variety of tumors and is associated with poor prognosis ([159]Ruggero et al., 2004). These results suggest that RBPs may influence the occurrence and progression of tumors by regulating multiple biological processes including RNA processing, RNA metabolism, RNA transport, translation regulation and mRNA surveillance pathway. In addition, we performed univariate Cox regression analysis, LASSO regression analysis and multivariate Cox regression analysis on these differentially expressed RBPs, and 8 prognostic related RBPs including APOBEC3G, AUH, DAZL, EIF4A1, IGF2BP3, NR0B1, RPL36A, and TRMT1 were selected. APOBEC3G, a member of the Apolipoprotein B mRNA editing enzyme-catalyzed polypeptide (APOBEC) family, was found to be overexpressed in renal carcinoma tissues and cell lines ([160]Komohara et al., 2007), consistent with our results. [161]Olson et al. (2018) found that this family is the source of somatic mutations in tumor cells that drive tumor evolution and may be associated with tumor cell recurrence, metastasis, and treatment resistance. AUH was found to be under-expressed in RCC and significantly associated with poorer survival in patients ([162]Zhang et al., 2019), which is similar to our results. The DAZL mutation was found to be associated with testicular cancer ([163]Ruark et al., 2013). The main function of EIF4A1 is to release mRNA structure in combination with other translation factors ([164]Qi et al., 2013). EIF4A1 has been reported to be associated with malignant phenotypes of tumor cells, tumor-specific survival, and susceptibility to therapeutic drugs ([165]Nagel et al., 2010; [166]Liang et al., 2014). [167]Wei et al. (2019) found that miR-1284 inhibited the progression of gastric cancer by targeting EIF4A1. IGF2BP3 has been found to be overexpressed in a variety of tumors including lung ([168]Wang et al., 2003), colon ([169]Li et al., 2009), and liver cancers ([170]Jeng et al., 2008). Accumulating studies have shown that IGF2BP3 is a promising prognostic factor for a variety of cancers including gastric cancer and RCC ([171]Kim et al., 2014; [172]Tschirdewahn et al., 2019). NR0B1 is a member of the orphan receptor family and is normally expressed mainly in the adrenal cortex, ovaries and support cells ([173]Ikeda et al., 1996). Studies have found that NR0B1 is abnormally expressed in endometrial cancer, prostate cancer, lung cancer and other cancers, and plays an important role ([174]Saito et al., 2005; [175]Seo et al., 2007; [176]Nakamura et al., 2009). [177]Oda et al. (2009) found that NR0B1 mainly affects tumor cell invasion, colony formation and tumorigenic activity, and is related to the malignant potential of lung adenocarcinoma. RPL36A mainly encodes ribosomal protein L36a. [178]Kim et al. (2004) found that overexpression of RPL36A in hepatocellular carcinoma was associated with enhanced cell proliferation, and RPL36A may be a potential target for anticancer therapy for hepatocellular carcinoma. [179]Alshabi et al. (2019) also found that high expression of RPL36A was associated with the tumorigenesis of glioblastoma multiform. [180]Nagel et al. (2010) found that TRMT1 was involved in the activation of LYL1 in leukemia cells and thus affected the differentiation of lymphocytes. GSEA analysis results showed that the regulation of RBPs was mainly concentrated in patients in the high-risk group, indicating that RBPs mainly regulates and affects patients in the high-risk group. However, the exact molecular mechanisms are unknown, and further exploration of possible mechanisms may be valuable. Subsequently, we constructed a prognostic model for OS based on these 8 RBPs to predict the prognosis of ccRCC patients. Survival analysis and ROC curve analysis showed that the model has good predictive performance. We then plotted a nomogram to establish a quantitative assessment method to predict the survival probability of ccRCC patients. According to our prognostic model for OS, patients with poor prognosis can be screened out, which may be conducive to timely adjustment of treatment regimens and individualized treatment. Further analysis showed that the prognostic model for OS could independently predict the prognosis of ccRCC patients and was associated with the progression of ccRCC tumors. And the results of Kaplan-Meier Plotter online tool analysis showed that all 8 prognostic RBPs were related to OS in ccRCC patients. Moreover, we constructed an RBPS-related prognostic model for DFS, showing that this prognostic model can also independently predict the prognosis of ccRCC patients. Overall, our study provides new insights into the occurrence and progress of ccRCC. In addition, the prognostic models for OS and DFS based on prognostic RBPs have good predictive performance, which are helpful to improve the clinical treatment decision and monitor the prognosis of patients. However, there are limitations in our study. First, our study is mainly based on a single bioomics information, and different characteristics of different platforms may lead to patient heterogeneity. Second, the model construction and validation of this study were designed by retrospective analysis, and the model still needs to be validated through a prospective clinical cohort. Moreover, the lack of clinical prognostic information in the study analysis may reduce the reliability of statistics. Finally, the prognostic models for OS and DFS based on prognostic RBPs showed good predictive performance. However, the exact molecular mechanisms of these prognostic RBPs involved in the occurrence, progression, and prognosis of renal cancer are still unclear, and the possible molecular mechanism and biological function need to be further explored. Conclusion We systematically analyzed the biological function and prognostic value of RBPs in ccRCC by using a variety of bioinformatics techniques. These RBPs may be involved in the pathogenesis, progression and metastasis of tumors. For the first time, we established prognostic risk score models for OS and DFS based on prognostic RBPs, and revealed they are independent prognostic factors related to OS and DFS in ccRCC patients. Our results are helpful to understand the molecular mechanism of ccRCC from a new perspective and to develop new prognostic markers or therapeutic targets. Data Availability Statement Publicly available datasets were analyzed in this study. This data can be found here: [181]https://portal.gdc.cancer.gov/. Author Contributions YW designed the study and performed the data analysis. XW, HF, BH, and BL performed the data analysis. YL, YR, XL, ZL, SW, and JL performed the data analysis and revised the manuscript. TW designed the study and revised the manuscript. All the authors read and approved the final manuscript. Conflict of Interest The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Abbreviations RCC renal cell carcinoma ccRCC clear cell renal cell carcinoma TCGA The Cancer Genome Atlas RBPs RNA binding proteins GO Gene Ontology KEGG Kyoto Encyclopedia of Genes and Genomes FC fold change OS overall survival DFS disease-free survival LASSO least absolute shrinkage and selection operator ROC receiver operating characteristic AUC area under the receiver operating characteristic curve FDR false discovery rate. Funding. This study was supported by the Medical Youth Top Talent Program of Hubei Provincial. ^1 [182]https://portal.gdc.cancer.gov/ ^2 [183]http://www.bioconductor.org/packages/release/bioc/html/edgeR.html ^3 [184]https://www.ebi.ac.uk/arrayexpress/ ^4 [185]https://www.cbioportal.org/datasets ^5 [186]http://www.webgestalt.org/ ^6 [187]https://www.gsea-msigdb.org/gsea/downloads.jsp ^7 [188]https://www.proteinatlas.org/ ^8 [189]https://kmplot.com/analysis/ Supplementary Material The Supplementary Material for this article can be found online at: [190]https://www.frontiersin.org/articles/10.3389/fgene.2020.617872/ful l#supplementary-material [191]Click here for additional data file.^ (299.3KB, TIF) [192]Click here for additional data file.^ (14.3KB, DOCX) References