Workflow schedule of this study

The technology schedule of this study was presented in Fig. 1. Firstly, we explored the differentially expressed RBPs between paracancerous and CESC, and the corresponding functional enrichment analysis was also conducted. Secondly, the TCGA-CESC cohort was separated into a training cohort (208 patients) and a testing cohort (88 patients) according to the ratio of 7:3. Univariate Cox regression, lasso Cox regression and stepwise multivariate Cox regression were employed to select candidate hub RBPs and a 10-RBPs signature was finally constructed in the training cohort. The risk signature was evaluated through survival analysis, multivariate Cox analysis and tAUC in the training cohort and then validated in the testing cohort. Thirdly, the training cohort and testing cohort were integrated into a pooled cohort and the patients with full scale of information were retained for further analysis. The decision tree was established to optime the risk stratification of CESC patients and the nomogram was constructed to quantify the risk value of CESC patients. Then the nomogram was assessed through calibration curves, tAUC and DCA. Fourth, the study investigated the mutational landscape, effect on immunotherapy and chemotherapy response of the RBPs’ signature through multiple bioinformatics methods and online webtools. In addition, the expression of PRPF40B was detected in CESC cell lines and clinical tissues, respectively. Finally, the biological functions of PRPF40B were explored by CCK-8, wound healing and transwell assay in vitro experiments.

Fig. 1

Schematic diagram of this research

Construction of 10-RBPs signature in CESC

Firstly, a total of 398 DEGs were identified between the paracancerous and tumor samples. The heatmap and volcano plot were showed in Fig. 2A, B. Functional enrichment analysis indicated that the DEGs are mainly associated with mRNA processing and RNA splicing biological process, in addition, KEGG results presented these DEGs mainly participate in mRNA surveillance and spliceosome pathway (Fig. 2C, D). Secondly, we performed a series of bioinformatics methods to select candidate genes and establish the risk signature. TCGA-CESC cohort was randomly separated into a training cohort (208 patients) and a testing cohort (88 patients), the detailed information was listed in Additional file 1: Table S3. Univariate Cox regression showed that 26 RBPs were significantly associated with prognosis among CESC (Fig. 2E, Additional file 1: Table S4). Lasso regression indicated that the partial likelihood deviance was minimal when log(lambda) equal to − 3.846, and 16 genes were obtained for further analysis (Fig. 2F, Additional file 1: Figure S1). The 10-RBPs signature was finally constructed by multivariate stepwise Cox regression (Fig. 2G). The detailed information of Cox result was presented in Table 1 and the risk score was calculated as mentioned above.

Fig. 2

Construction the risk signature in training cohort. A Heatmap was drawn to show the DEGs between normal cervical tissues and CESC with the following criterion: |logFC|> 1 and adjust P value < 0.05; B A total of 398 DEGs were presented in volcano plot, up-regulated DEGs were displayed as red points and down-regulated DEGs were displayed as green points; C Top 10 enriched biological process of these 398 DEGs; D KEGG enrichment results of these 398 DEGs; E 26 DEGs showed statistical significance (P < 0.05) through univariate Cox regression, DEGs with HR < 1 were presented as blue bars and DEGs with HR > 1 were presented as red bars, the bars’ length represent the 95% confidence interval; F 16 DEGs were selected through lasso Cox regression by tenfold cross validation, red represent the risk DEGs and blue represent the protective DEGs; G 10 DEGs were finally gathered through step-wise multivariate Cox regression to construct the risk signature

Table 1 Multivariate Cox results of the 10 RBPs in the risk signature10-RBPs signature performed satisfactory prognostic ability both in training cohort and testing cohort

Next, we evaluated the 10-RBPs signature in the training cohort firstly. As shown in Fig. 3A, the transcript levels of these 10 RBPs were presented in a heatmap, the five protective genes were higher expressed in low-risk group, and the remaining five risk genes were higher expressed in high-risk group. The risk score was significantly higher in dead patients compared with alive patients (P < 0.001). Meanwhile, the proportion of dead patients in high-risk group was much higher. Survival plot indicated that the patients in high-risk group exhibit significantly shorter overall survival time (P < 0.001). Multivariate Cox regression suggested that the 10-RBPs signature served as an independent prognostic factor among CESC patients (P < 0.001). tAUC plot showed the 10-RBPs signature was much more accurate than other clinical parameters including age, stage and grade. The bar plot was used to show the quantitative results of tAUC. Furthermore, the 10-RBPs signature was validated in testing cohort and similarity results were presented in Fig. 3B.

Fig. 3

Evaluation the 10-RBPs signature in training cohort (A) and testing cohort (B). The transcript difference of 10 RBPs between high-risk group and low-risk group were presented in heatmap; Boxplot to show the difference of risk value between alive patients and dead patients; Stacked column chart was employed to present the proportion of alive and dead patients between high-risk group and low-risk group; Survival plot was conducted to exhibit the prognostic difference between high-risk group and low-risk group; Forest plot presented the multivariate Cox regression of 10-RBPs signature and clinical parameters; Prediction ability of 10-RBPs signature and other clinical factors were displayed in tAUC curves; Barplot to show the difference of AUC towards various indicators

The nomogram could effectively predict the prognosis among TCGA-CESC patients

To explore the clinical application of the 10-RBPs signature, we combined the training cohort and testing cohort and 263 CESC patients with full-scale clinical information were extracted. Multivariate Cox regression indicated the pathological stage (P = 0.020) and risk score (P < 0.001) were both remarkably associated with prognosis among TCGA-CESC patients (Fig. 4A). Hence, the stage and risk scores were submitted to establish the decision tree. As shown in Fig. 4B, three different subgroups were identified based on two major components in which the risk score showed the most powerful effect. Patients with lower risk scores were labeled as “low-risk”, while higher risk scores & stage I–III and higher risk & stage IV were labeled as “intermediate-risk” and “high-risk”, respectively. Survival analysis presented significant differences between the subgroups (Fig. 4C). To quantify the individual risk assessment with CESC, a nomogram was developed to achieve this goal (Fig. 4D). The red arrow showed an illustration of the calculation process. Calibration curves indicated that the nomogram performed good performance in predicting prognosis (Fig. 4E). tAUC illustrates the nomogram acts with much higher accuracy in prognosis prediction than else clinicopathological features (Fig. 4F). DCA curves indicated the nomogram brings about more net benefit for CESC patients than other clinical parameters in different cut-offs of survival time (Fig. 4G).

Fig. 4

Assessment the clinical application of the 10-RBPs signature in TCGA-CESC cohort. A Forest plot to show the multivariate Cox regression of age, grade, stage and risk signature; B A decision tree was established by combining the risk signature and pathological stage for the purpose of optime the risk stratification; C Significant survival difference was observed between the three subgroups; D Details of the nomogram, red arrow represent the calculate process of a random selected patient; E Calibration curves of the nomogram; F Comparison of the tAUC between nomogram, age, grade and stage; G DCA curves to show the net-benefit of nomogram at 3-year, 5-year and 10-year time point

10-RBPs signature may participate in various cancer-related pathways and play a vital role in tumor progression

To comprehensively explore the biological functions of the risk signature, we performed enrichment analysis in TCGA-CESC based on Hallmark gene sets. Firstly, the 50 Hallmark gene sets were quantified by “ssgsea” and the results showed most gene sets were significantly higher in high-risk group compared with low-risk group, including hypoxia, glycolysis and EMT, et al. (Fig. 5A). GSEA results indicated most cancer-related pathways were enriched in high-risk group, including angiogenesis, EMT, hypoxia, et al., and the immune-related pathways were enriched in low-risk group, including IFN-alpha response and IFN-gamma response (Fig. 5B). Survival analysis presented that these cancer-related pathways were positively associated with poor prognosis, which is consistent with the above results (Fig. 5C). In addition, the landscape of somatic mutation was also explored between different risk groups. The top 20 mutated genes in different groups were presented and compared (Additional file 1: Figure S2A, B). PIK3CA showed the highest mutation frequency in high-risk group, while TTN showed the highest mutation frequency in low-risk group. Fisher’s exact test revealed that PDE3A, STK11, APC, BAHCC1, DSCAML1, TCOF1, SPEN and DNAH3 were higher mutated in the high-risk group compared with low-risk group (P < 0.05, Additional file 1: Figure S2C). The lollipop plot showed the different mutational sites of PDE3A (Additional file 1: Figure S2D). Besides, co-occurrence and mutually exclusive mutations were also investigated. From the results, we know that most co-occurrence patterns happened in high-risk group, which indicated a potential common effect induced by their mutations (Additional file 1: Figure S2E, F).

Fig. 5

Functional enrichment analysis between high-risk group and low-risk group. A Hallmark gene sets were quantified by “GSVA” and compared between high-risk group and low-risk group; B GSEA to explore the difference of pathway between high-risk group and low-risk group; C Cancer-related pathways presented significant survival difference in CESC

10-RBPs signature showed a significant relationship with immune cell infiltration and immunotherapy response

A total of 3608 DEGs were identified between high-risk group and low-risk group, the results were presented in volcano plot (Fig. 6A). Enrichment analysis indicated the up-regulated genes were mainly participate in cancer-related pathways, including glycolysis/gluconeogenesis, HIF-1 signaling pathway, focal adhesion and Hippo signaling pathway, et al. (Fig. 6B). The down-regulated genes were enriched in immune-related pathways, for example, the T cell receptor signaling pathway (Fig. 6C). These results indicated the 10-RBPs signature may influence the tumor-immune interaction. Hence, we quantified and compared the infiltration of 28 immune cells between different risk groups. As shown in Fig. 6D, most immune cells presented higher infiltration in low-risk group. Similarly, the immune checkpoint genes were significantly higher expressed in low-risk group (Fig. 6E). Furthermore, the IPS, IPS-CTLA4 blocker, IPS-PD1/PD-L1/PD-L2 blocker and IPS-CTLA4- PD1/PD-L1/PD-L2 blocker were remarkably elevated in low-risk group, indicated the patients in low-risk group may present more sensitivity to immunotherapy response (Fig. 6F). Finally, the IMvigor210 cohort was used to evaluate the prognostic ability of 10-RBPs signature. As shown in Fig. 6G, boxplot showed the risk score was significantly higher in “SD/PD” group compared with “CR/PR” group, from which we know the risk score was negatively associated with immunotherapy response. Survival plot declared the 10-RBPs signature also act good performance in IMvigor210 cohort (P < 0.001). As expected, the high-risk group possesses a higher proportion of patients who showed “SD/PD” response to immunotherapy.

Fig. 6

The effect of 10-RBPs on immunotherapy. A A total of 3608 DEGs were showed in volcano plot (P < 0.05); B KEGG enrichment of the up-regulated DEGs; C KEGG enrichment of the down-regulated DEGs; D Comparison of the 28 immune cells between high-risk group and low-risk group; E Comparison the transcript difference of 10 immune checkpoints between high-risk group and low-risk group; F Violin plot presented the differences of IPS scores between high-risk group and low-risk group; G Exploration of the immunotherapy response of 10-RBPs signature in IMvigor210 cohort

10-RBPs signature served as a potential biomarker in chemotherapy resistance

GSEA results indicated that two chemotherapy-related pathways were enriched in high-risk group, including gefitinib resistance pathway and endocrine therapy resistance pathway (Fig. 7A). GSCALike web tool provided the Spearman correlation between gene expression and IC50 drug data, the positive correlation means the gene’s high expression is resistant to the drug, and vice versa. The results presented that DDX26B, RBM38 and ANGEL2 were negatively correlated with IC50 data, while PRPF40B was positively correlated with IC50 data, providing certain evidence for the effect of 10-RBPs on chemotherapy resistance (Fig. 7B). In addition, we performed chemotherapy prediction in TCGA-CESC and the results indicated the high-risk patients showed higher IC50 in several chemotherapy drugs, including Roscovitine, BMS.536924, PF.02341066, Rapamycin, Sunitinib, VX.680, Bortezomib, LFM.A13, Metformin, NVP.TAE684, MS.275 and Methotrexate (Fig. 7C). Besides, a total of 222 TCGA-CESC patients who received adjuvant therapy were selected and significant survival difference was observed between different risk groups among these patients. We further divided these 222 patients into three subgroups labeled “CR”, “PR/SD” and “PD”, respectively. Results indicated the “PD” subgroup possessed highest risk score compared with other subgroups, while the “CR” subgroup presented the lowest risk score. Meanwhile, the proportions of “PD” and “CR” patients in high-risk group and low-risk group were significantly different (Fig. 7D).

Fig. 7

The effect of 10-RBPs on chemotherapy. A GSEA showed the high-risk group was positively associated with chemotherapy resistance pathways; B Correlation between transcript expression and IC50 data retrieved from GDSCLike database; C Comparison the estimated IC50 data between high-risk group and low-risk group in TCGA-CESC; D Evaluate the chemotherapy response of 10-RBPs in CESC

PRPF40B was up-regulated in CESC tissues and cell lines

From the above results, we observed PRPF40B may act as a prognostic biomarker in CESC. Hence, the study verified the expression level of PRPF40B in clinical tissues and CESC cell lines. Firstly, survival analysis indicated that higher expression of PRPF40B showed a positive association with worse OS in TCGA database (Fig. 8A). Secondly, results of RT-qPCR indicated that the relative mRNA expression of PRPF40B was significantly higher in CESC tissues compared with paracancerous tissues (P < 0.05, Fig. 8B). The translation level of PRPF40B was verified through western blot (P < 0.05, Fig. 8C, D). In addition, the protein level of PRPF40B was also increased in CESC cell lines Hela and Siha (P < 0.05, Fig. 8E, F).

Fig. 8

Verify the expression level of PRPF40B in CESC tissues and cell lines. A Survival plot of PRPF40B in TCGA database; B the transcript level of PRPF40B in CESC tissues and paired normal tissues; C, D the translation level of PRPF40B in CESC tissues and paired normal tissues, and the quantitative results; E, F the translation level of PRPF40B in four cervical cell lines and the quantitative results. Values are expressed as the mean ± SD (n = 3); *P < 0.05 versus vehicle-treated control group; A.U. arbitrary units

PRPF40B accelerates the proliferation, migration and invasion of CESC in vitro

To further evaluate the biological functions of PRPF40B, we conduct experimental validation for this hypothesis. Western-blot presented that the PRPF40B was successfully knockdown in Hela and Siha cell lines (Fig. 9A, G). CCK-8 assay suggested that the proliferation of Hela and Siha cells was significantly suppressed in the siPRPF40B group (P < 0.05, Fig. 9B, H). The wound healing assay indicated that siPRPF40B significantly inhibits the growth rate of Hela and Siha cells (Fig. 9C, D, I, J). Consistent with the above results, transwell assay suggested the migration and invasion ability of these two cells were obviously suppressed after siPRPF40B was treated (Fig. 9E, F, K, L). These results declared that PRPF40B acts as an oncogene in CESC, which is in accordance with the prognostic effect, even though the underlying mechanism is unknown.

Fig. 9

Investigate the biological function of PRPF40B in CESC cell lines. A, G Western-blot displayed the siPRPF40B model in Hela and Siha cell lines; Proliferation assay of vector and siPRPF40B in Hela cell (B) and Siha cell (H); Wound-healing assay of vector and siPRPF40B in Hela cell (C, D) and Siha cell (I, J); Transwell assay of vector and siPRPF40B in Hela cell (E, F) and Siha cell (K, L). Values are expressed as the mean ± SD (n = 5); *P < 0.05 versus vehicle-treated control group