The gene mutation, expression and clinical information of HCC patients was downloaded from the TCGA and ICGC database.

There were huge number of RNA regulators identified in eukaryotes. However, we only choose RNA regulators which were proved to play critical roles in HCC and other cancers (Delaunay and Frye, 2019; Barbieri and Kouzarides, 2020; Han et al., 2021). As a result, a total of one hundred RNA regulators which belong to eight types of RNA modifications (m⁶A, m¹A, m⁵C, Nm, m⁷G, Ψ, A-to-I, and mcm⁵s²U) were included in our study. They were listed in Supplementary Table S1.

Unpaired t-test were used to compare the expression of 100 RNA regulators between tumor and adjacent normal tissues. The “Maftools” R package was used to analyze the mutation and interaction of 100 RNA modification genes. The protein-protein interaction between 100 RNA modification regulators were analyzed by the STRING database and visualized by using Cytoscape software (Szklarczyk et al., 2019).

To figure out the relationship between the expression of 100 RNA modification regulator and HCC subtypes, consensus clustering analysis was conducted in TCGA-LIHC patients by using the “ConsensusClusterPlus” R package.

To reveal the biological function enriched in different group of patients, the Gene Set Variation Analysis (GSVA) was performed by utilizing the R package “GSVA”. The annotated gene set “h.all.v7.4.symbols.gmt” was used for GSVA analysis. The enrichment scores of immune cells in sample from different groups were determined with single-sample gene set enrichment analysis (ssGSEA).

To construct the RMScore, we first identified 1,427 differentially expressed genes (DEGs) in cluster 1 and cluster 2 groups. We divided TCGA-LIHC into training set and validation set. In the training set, univariate cox regression analysis revealed that 266 DEGs were identified to be significantly related to the prognosis of HCC patients (p < 0.001). Then LASSO regression analysis was conducted to select survival-associated DEGs. Multivariate cox regression analysis was conducted to construct the final prognostic model. The RMScore was calculated by the following equation: RMScore = Exp1* Coe 1 + Exp2* Coe 2 + … Exp N* Coe N.

The validation dataset of TCGA-LIHC, whole TCGA-LIHC dataset, and ICGC dataset were utilized as external validation dataset. The RMScore of each sample was calculated by the formula presented above. Patients were divided into low RMScore and high RMScore groups by the median value of the RMScore. The “Survival” R package was conducted to analysis the survival of two groups. The “SurvivalROC” R package was used to plot the ROC curves.

To predict the sensitivity to chemotherapeutic drugs, we use the R package “pRRophetic” to estimate the half-maximal inhibitory concentration (IC50) in individual patients. The R package “ggplot2” was utilized to generate the plots.

The stage and RMScore, both of which were identified to be an independent prognostic value by multivariate Cox proportional hazards analysis, were included to construct the nomogram. The calibration curve was plotted to observe the prediction probabilities of the nomogram against the observed rates. The prognostic ability of the nomogram was evaluated by receiver operating (ROC) curves and decision curve analysis (DCA) curves.

We choose a total of one hundred RNA modification regulators which were reported to play roles in at least one kind of cancer to perform following analysis (Supplementary Table S1).

In order to characterize the expression profile of 100 RNA regulators, unpaired t-test was conducted to compare the gene expression between tumor and normal tissues. Notably, nearly 90% RNA regulators exhibited significantly higher expression in tumors than normal tissues (Figure 1A). Among the eight types of RNA modification regulators, all regulators in mcm⁵s²U, Ψ and A-to-I showed elevated expression in tumor tissues. The dysregulation of these RNA regulators may play key roles in the development of HCC.

The mutation profiles of 100 hundred regulators in HCC were depicted together. As a results, 94 of 371 samples (about 25.34%) showed RNA modification regulator mutations. M⁵C eraser TET1 and A-to-I writer ADARB2 had a mutation rate of 2%. About 40% RNA regulators had a mutation rate of 1% (Figure 1B). Genetic interaction analysis showed that mutation of m⁵C eraser TET3 was positively correlated with mutations of TRMT10B, ADAT2, PRRC2A, YTHDC1 and TET2. Mutation of m⁶A reader HNRNPA2B1 was positively correlated with mutations of WTAP, ELAVL1, MRM2, HNRNPC, HNRNPA1, YTHDC1, PRRC2A. Mutation of m⁶A writer ZC3H13 was positively correlated with mutations of TRMT10B, TET3, TRMT13, PRRC2A, YTHDC1, ADAR, YTHDC2 and ZCCHC4, while negatively correlated with mutations of METTL1 and CTU2 (Figure 1C).

To reveal the interactive relationships between proteins encoded by 100 RNA regulators, we constructed a PPI network using STRING and Cytoscape software. Correlations with adjusted p-value <0.05 and |r| ≥ 0.6 were shown in Figure 1D. We can see that these RNA regulators, especially regulators of m⁶A and m⁵C, showed close correlations with each other.

Results above collectively demonstrated that there were huge interaction and crosstalk between these 100 RNA regulators. They may form a complicated network to synergistically mediate the development and progression of HCC.

In order to excavate the RNA modification patterns of HCC, consensus clustering analysis were conducted in 374 patients in TCGA-LIHC dataset. Patients in TCGA-LIHC dataset were separated into diverse clusters (K = 2–9) according to the expression of 100 RNA modification genes via an unsupervised consensus clustering analysis. The cumulative distribution function (CDF) curve and the area under the CDF curve was plotted, which indicated that the cohort was well distributed when K = 2 (Figures 2A–C). Kaplan-Meier survival analysis showed that patients in cluster 1 had a significantly better prognosis than patients in cluster 2 (p < 00.001) (Figure 2D). GSVA and ssGSEA enrichment analysis were performed to analyze biological process in two clusters. As shown in Figure 2E, most signaling pathways including DNA repair, MYC targets, G2M checkpoint, E2F targets were significantly enriched in Cluster 2. These may lead to the high degree of malignancy and poor prognosis of tumors in Cluster 2. Signaling pathways such as oxidative phosphorylation, reactive oxygen species pathway, cholesterol homeostasis and fatty acid metabolism were enriched in Cluster 1. There were more NK cells, eosinophils, neutrophils, effector memory CD8 T cells in Cluster 1, while more activated CD4 T cells, T follicular helper cells, type 2 T helper cell and effector memory CD4 T cells were accumulated in Cluster 2. This result implied two clusters had different immune infiltrating patterns.

Results above collectively suggested that 100 RNA modification regulators defined two RNA modification patterns in HCC, which exhibited different malignancy, immune microenvironment, and prognosis.

To further investigate the potential genetic alterations between two RNA modification patterns, differential expression analyses were carried out between two clusters (Figure 3A). A total of 1,472 differentially expressed genes (DEGs) were identified (Supplementary Table S2). Next, we divided the TCGA-LIHC queue into a training cohort and a verification cohort. Among the 1,427 DEGs, 266 genes were associated with the prognosis of HCC patients in training cohort by univariate COX regression analysis (Supplementary Table S3). Then, the LASSO regression analysis was conducted to further narrow the survival-related DEGs (Supplementary Table S4). Finally, six genes were identified by stepwise multivariate regression analysis and subsequently used to construct the RNA modification score (RMScore). The six genes identified were CYP2C9, CBX2, NDRG1, EPS8L3, FAM83D and MYCN. The RMScore = 0.1034*Exp of CBX2 + 0.0087*Exp of NDRG1 + 0.0392*Exp of EPS8L3 + 0.0517*Exp of FAM83D + 0.0596*Exp of MYCN−0.0031*Exp of CYP2C9.

By Kaplan-Meier analysis, we found that patients in low-risk group had a clear survival advantage over the high-risk group in training set (p < 0.001) (Figure 3B). The same survival advantage was also observed in the validation cohort (p = 0.033) and the ICGC cohort (p < 0.001) (Figures 3D, F). The time-dependent ROC curves were utilized in assessing the performance of RMScore, which showed the RMScore had good sensitivity and specificity in prognosis prediction (Figures 3C, E, G). These results suggested that the RMScore had good prediction ability for the survival of HCC patients.

The RMScore together with clinicopathological factors, such as age, gender, grade, and stage were included for univariate and multivariate Cox regression analysis. In the univariate analysis, the RMScore was significantly related to the prognosis of HCC patients (hazard ratio [HR], 1.053; 95% confidence interval [CI], 1.037–1.069; p < 0.001) (Figure 4A). Multivariate analysis also demonstrated that the RMScore was an independent risk of the prognosis of HCC patients (hazard ratio [HR], 1.045; 95% confidence interval [CI], 1.028–1.062; p < 0.001) (Figure 4B). In the heatmap incorporating gene expression and clinicopathological traits, the majority of patients in Cluster 1 were gathered in low RMScore group, which demonstrated that the RMScore well represented the traits of patients in Cluster 1. Among the six gene which constructing the RMScore, only CYP2C9 was upregulated in the Cluster 1 as well as low RMScore group. The other five genes (CBX2, NDRG1, EPSBL3, FAM83D and MYCN) were upregulated in Cluster 2 as well as high RMScore group (Figure 4C). We explored the pathway enrichment associated with the RMScore to identify the internal mechanisms involved. Similar to the pathway enrichment results between Cluster 1 and Cluster 2, the high RMScore group was mainly enriched in signaling pathways related to DNA repair, MYC targets, G2M checkpoint, and E2F targets, while the low RMScore group was enriched in oxidative phosphorylation, reactive oxygen species pathway, cholesterol homeostasis and fatty acid metabolism (Figure 4D).

We investigated the association between the RMScore and the sensitivity to chemotherapeutics in HCC patients. As a result, significant differences in the estimated IC50 between the two risk groups were observed for common chemotherapy or targeted drugs including paclitaxel (Figure 5A), sorafenib (Figure 5B), 5-flurouracil (Figure 5C), doxorubicin (Figure 5D), gemcitabine (Figure 5E) and erlotinib (Figure 5F). Lower RMScore was related to higher IC50 among paclitaxel, sorafenib, 5-flurouracil, gemcitabine, suggesting that patients with low RMScore were less sensitive to these treatments. However, patients with low RMScore had a lower IC50 to erlotinib, suggesting that they may be more sensitive to erlotinib than high RMScore patients.

To facilitate the clinical applicability and availability of the RMScore, a predictive nomogram for 1-, 3-, 5-year OS combined with tumor stage and RMScore was developed (Figure 6A). The calibration curves showed that actual overall survival of patients was almost consistent with the predicted overall survival of patients (Figure 6B). The DCA plot suggested that the nomogram had a superior efficiency than RMScore in predicting patient survival outcomes (Figure 6C). The AUCs of the nomogram in the 1-, 3-, and 5-year ROC curves were 0.750, 0.743 and 0.735, respectively (Figure 6D).