RNA-seq transcriptome data and complete clinical information of 374 HCC and 50 control specimens were retrieved from the Cancer Genome Atlas (TCGA) database¹. Among them, 108 HBV-related HCC were extracted from TCGA dataset, which were employed as the discovery set. Expression profiling and clinical features of 224 HBV-related HCC samples were obtained from the GSE14520 dataset in the Gene Expression Omnibus (GEO;²) repository. This GSE14520 dataset was based on the GPL571 and GPL3921 platforms. This dataset was used as the validation set. Table 1 listed the clinical characteristics of HBV-related HCC patients. Each probe was transformed to the corresponding gene symbol. If multiple probes matched the same gene symbol, the average value was determined as the expression value of this gene. Based on published articles, a list of 1,542 RBPs was collected (Supplementary Table 1).

Differences in expression levels of RBPs were analyzed between 108 HBV-related HCC and 50 control samples via limma package based on RNA-seq transcriptome data (Ritchie et al., 2015). Under the threshold of | fold-change| > 2 and adjusted p < 0.001, up- and down-regulated RBPs were screened for HBV-related HCC.

Gene ontology (GO) enrichment analysis primarily contains three categories: biological process (BP), cellular component (CC) and molecular function (MF). Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) can provide signaling pathways involved in RBPs. Here, GO, and KEGG enrichment analyses of differentially expressed RBPs were carried out via clusterProfiler package (Ashburner et al., 2000). Terms with false discovery rate (FDR) < 0.05 were significantly enriched by above RBPs. Functional association between HBV-related RBPs was analyzed through the STRING database³ (Szklarczyk et al., 2017). A protein-protein interaction (PPI) network was conducted utilizing Cytoscape software (Doncheva et al., 2019).

Univariate Cox regression analysis was applied for analyzing the associations between overall survival (OS) of HBV-related HCC patients and differentially expressed RBPs utilizing survival package. Prognosis-related RBPs were determined with the threshold of p < 0.05. Then, key prognosis-related RBPs were screened through LASSO regression analysis with glmnet package (Engebretsen and Bohlin, 2019). Normalized regression coefficients of key prognosis-related RBPs were calculated based on multivariate regression analysis.

A risk score was established on the basis of the key prognosis-related RBPs for HBV-related HCC patients in the TCGA dataset. This study calculated the risk score of each patient, according to the formula: risk score = expression of RBP1 × regression coefficient of RBP1 + expression of RBP2 × regression coefficient of RBP2 + … + expression of RBPn × regression coefficient of RBPn. According to the median risk score, HBV-related HCC patients were assigned into high- and low-risk groups. Hierarchical clustering analysis was applied for showing the associations between expression patterns of above RBPs and clinical characteristics (stage, grade, gender, and age). Kaplan-Meier curve analyses were presented for investigating the differences in OS, disease-free survival (DFS) and progression-free interval (PFI) between two groups using Survival package. P values were determined with log-rank tests. Time-dependent receiver operating characteristic (ROC) curves under one-, three- and five-year OS, DFS, and PFI were generated by SurvivalROC package (Heagerty et al., 2000). Then, the area under the curve (AUC) was calculated to evaluate the predictive usefulness of the risk score.

The associations between risk score, age, gender, grade, stage, and prognosis were evaluated utilizing univariate cox regression analysis. Prognostic factors with p < 0.05 were incorporated for multivariate cox regression analysis. Independent prognostic factors were then identified for HBV-related HCC. Time-independent ROCs of risk score, age, gender, grade, and stage were separately constructed and the predictive power was compared. In published literature, Fang and Chen (2020) proposed a two-m⁶A RNA methylation regulator prognostic model (HNRNPA2B1 and RBM15) for HBV-related HCC prognosis. By ROCs, the predictive efficacy was compared with our prognostic model. For evaluating the clinical generalizability of the model, this model was verified in the GSE14520 dataset. With the same formula, the risk scores of patients were calculated in this dataset. Based on the median risk score, OS differences between high- and low-risk groups were assessed by Kaplan-Meier curves and log-rank tests. The predictive performance was verified by ROCs.

HBV-related HCC subjects in TCGA dataset were stratified into subgroups according to clinical characteristics, as follows: age ≥ 65 and < 65 subgroups, female and male subgroups, grade 1–2 and 3–4 subgroups and stage I-II and stage III-IV subgroups. In each subgroup, Kaplan-Meier curves of OS were conducted between high- and low-risk patients. Survival differences were estimated with log-rank tests.

A nomogram by incorporating gender, age, grade, stage, and risk score was created for predicting one-, three-, and five-year survival probabilities of HBV-related HCC patients in the TCGA dataset. The accuracy in predicting prognosis was evaluated by ROC curves, calibration curves and decision curve analysis (DCA) (Vickers and Elkin, 2006). ROC curves were conducted for evaluating the predictive performance of this nomogram for one-, three- and five-year OS. By calibration curves, discrepancy between nomogram-estimated and actual one-, three-, and five-year survival duration was analyzed. DCA was applied for quantifying the clinical practical use with survival outcomes of a decision considered.

HBV-related HCC specimens were stratified into high- and low-risk groups. By Gene Set Enrichment Analyses (GSEA)⁴, potential mechanisms of this prognosis-related model were elucidated (Subramanian et al., 2005). GSEA was carried out for finding enriched KEGG pathways, with KEGG gene set “c2.cp.kegg.v7.0.symbols.gmt” as a reference. Pathways with nominal p < 0.05 were distinctly enriched.

By applying CIBERSORT algorithm⁵, the proportions of 22 immune cells in HBV-related HCC and control specimens were quantified based on their expression profiling (Newman et al., 2015). The proportions of all immune cells in each specimen were equal to 1. The LM22 signature was employed as a reference set. The permutations were set as 1,000. Samples with p < 0.05 were screened for further analyses. The correlations between risk score and infiltrations of immune cells and HLA family expression were assessed with Spearson correlation analysis.

This study collected 1,542 RBPs and analyzed their expression in 108 HBV-related HCC and 50 control tissues. | Fold-change| > 2 and adjusted p < 0.001 were set as the screening thresholds. As a result, 340 RBPs were up-regulated in HBV-related HCC than control specimens (Figures 1A,B). The first 20 up-regulated RBPs were shown in Table 2. Meanwhile, there were six down-regulated RBPs (GSPT2, PAIP2B, ANG, AZGP1, CPEB3, and ZFP36) in HBV-related HCC compared to controls (Figures 1A,B and Table 2). These RBPs were primarily enriched in mRNA metabolic biological processes such as nuclear-transcribed mRNA catabolic process, RNA catabolic process, nuclear-transcribed mRNA catabolic process, mRNA catabolic process and ncRNA processing (Figure 1C). Our KEGG analysis demonstrated that above RBPs were distinctly related to key pathways including ribosome, spliceosome, RNA transport, mRNA surveillance pathway, ribosome biogenesis in eukaryotes and RNA degradation, confirming their roles on post-transcriptional gene regulation (Figure 1D). A PPI network was conducted for revealing the interactions between these HCC-related RBPs (Figure 1E).

By univariate analyses, we investigated which dysregulated RBPs were in relation to HBV-related HCC patients’ survival outcomes. With the cutoff of p < 0.05, 54 RBPs were distinctly correlated to prognosis (Table 3). These prognosis-related RBPs were utilized for LASSO analysis. As a result, 10 key RBPs were screened for constructing a prognostic model (Figures 2A,B). The regression coefficients of above RBPs were as follows: POLR2L = 0.453206443450835; MRP S12 = 0.00334358058970182; DYNLL1 = 0.0556454917799665; ZFP36 = 0.215686692231002; PPIH = 0.300574508736105; RARS = 1.0223374037773; SRP14 = −1.08473362012663; DDX 41 = 0.00122492284215762; EIF2B4 = 0.348820298750318; NOL12 = 0.155516694418799. The risk score of each subject from TCGA dataset was determined according to expressions and regression coefficients of RBPs. Heatmap depicted the correlations between expression patterns of POLR2L, MRPS12, DYNLL1, ZFP36, PPIH, RARS, SRP14, DDX41, EIF2B4, and NOL12 and clinical characteristics of HBV-related HCC patients (Figure 2C). Patients from the TCGA dataset were assigned into high- and low-risk groups. As a result, high-risk subjects were indicative of unfavorable OS (p = 3.734e-06; Figure 2D), DFS (p = 6.495e-02; Figure 2E) and PFI (p = 2.135e-02; Figure 2F). The ROCs were plotted for evaluation of the predictive performance. The AUCs under one-, three-, and five-year OS were separately 0.900, 0.945 and 0.886, demonstrating that this model could be precisely predictive of OS probabilities (Figure 2G). Meanwhile, the AUCs under one-, three-, and five-year DFS were 0.686, 0.729 and 0.668 (Figure 2H) as well as the AUCs under one-, three-, and five-year PFI (Figure 2I) were 0.673, 0.784, and 0.667, which were indicative that this model possessed the well performance on predicting HBV-related HCC recurrence and progression.

By univariate analyses, we investigated the associations between survival outcomes and risk score and clinical parameters in the TCGA dataset. In Figure 3A, risk score (p < 0.001) and stage (p = 0.009) were both risk factors of HBV-related HCC. Following multivariate analyses, the risk score was independently predictive of survival outcomes (p < 0.00l; Figure 3B). In comparison to other clinical parameters, this risk score exhibited the highest AUC of OS (0.940). This demonstrated that the risk score possessed more excellent predictive performance than other parameters (Figure 3C). Compared to the published prognostic model (HNRNPA2B1 and RBM15), higher AUC value was investigated in our model (Figure 3D). We further evaluated the clinical generalizability of this model. In the GSE14520 dataset, high-risk scores were also indicative of unfavorable survival outcomes (p = 2.999e-03; Figure 3E) and the AUC value was 0.656 (Figure 3F). Collectively, this model displayed the independent and well power in predicting prognosis.

For investigating the predictive sensitivity of this model in HBV-related HCC prognosis, survival analyses were carried out in different subgroups. Our data demonstrated that high-risk scores were distinctly predictive of poorer survival outcomes than low-risk score in age ≥ 65 (p = 0.067; Figure 4A) and < 65 (p < 0.001; Figure 4B) subgroups, female (p = 0.042; Figure 4C) and male (p < 0.001; Figure 4D) subgroups, grade 1–2 (p = 0.027; Figure 4E) and grade 3–4 (p < 0.001; Figure 4F) subgroups and stage I-II (p < 0.001; Figure 4G) and stage III-IV (p = 0.044; Figure 4H) subgroups.

To facilitate personalized treatment, we constructed a nomogram by incorporating gender, age, grade, stage, and risk score in TCGA dataset. This nomogram was utilized for estimating one-, three-, and five-year survival probabilities (Figure 5A). ROC curves demonstrated the well performance on predicting one-, three- and five-year (Figure 5B). As shown in our calibration plots, nomogram-estimated one-, three- and five-year survival probabilities were close to actual survival consequences (Figures 5C–E). Meanwhile, DCA showed that this nomogram exhibited the best net benefit for one-, three- and five-year survival duration (Figures 5F–H). Hence, the nomogram model could assist clinical management and decision.

For observing potential pathways involved in unfavorable survival outcomes, GSEA was carried out. We found that endocytosis (Figure 6A), RNA degradation (Figure 6B), spliceosome (Figure 6C) and ubiquitin-mediated proteolysis (Figure 6D) were distinctly activated in high-risk HBV-related HCC specimens.

CIBERSORT algorithm was applied for inferring the proportions of 22 immune cells in HBV-related HCC tissues, including B cells naïve, B cells memory, plasma cells, T cells CD8, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, monocytes, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils and neutrophils (Figure 7A). There was the heterogeneity in the immune microenvironment among subjects. The close crosstalk between immune cells was found, as shown in Figure 7B. Also, the risk score was associated with mast cells resting, NK cells resting, neutrophils, mast cells activated, macrophages M0, Tregs and eosinophils. Moreover, the risk score displayed the significant correlations to HLA expression in HBV-related HCC specimens (Figure 7C). The data indicated that the risk score was related to the immune microenvironment of HBV-related HCC.