The gene expression profiles and clinical annotations were downloaded from the Cancer Genome Atlas (TCGA) portal (http://cancergenome.nih.gov/). Data cohorts with missing information were removed. A total of 356 cases of TCGA-LIHC were used for further analysis. The R Bioconductor package and R (version 3.6.2) were employed for data analysis.

Drug sensitivity data were collected from The Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerRxgene.org) (Yang et al., 2013). Spearman’s correlation analysis was utilized to evaluate the association between the scoring model and drug reaction, where |Rs| > 0.2, and FDR <0.05 was considered significant correlation.

The immunotherapy dataset IMvigor210 cohort was used to explore the immunotherapy response and prognosis of HCC patients with different WM-Score values. The standardized RNA-sequencing data of 1111 HCC patients with detailed clinicopathological data were downloaded from http://research-pub.gene.com. The data were analyzed using the IMvigor210CoreBiologies R package.

In order to explore the robust clustering of HCC cases, we employed the unsupervised clustering approach to analyze the gene profiles of RNA modification writers. A total of 26 RNA modification regulators, including seven m⁶A modification enzymes (KIAA1429, METTL14, ZC3H13, METTL3, WTAP, RBM15B, and RBM15), 12 APA modification enzymes (CPSF1, CPSF2, CPSF3, CPSF4, CSTF1, CSTF2, CSTF3, CF1, PCF11, CLP1, NUDT21, and PABPN1), four m¹A modification enzymes (TRMT10C, TRMT6, TRMT61A, and TRMT61B), and three A-I modification enzymes (ADARB1, ADARB2 and ADAR) were analyzed. An NMF-based consistent clustering algorithm was used to determine RNA modification patterns based on the mRNA expression of analyzed regulators. Unsupervised cluster analysis was performed by The Consensus Cluster Plus package as previously described (Wilkerson and Hayes, 2010).

GSVA is a gene set enrichment method that provides increased power to estimate changes of subtle pathway activity over a sample population in an unsupervised manner (Hänzelmann et al., 2013). We conducted GSVA analysis to explore the association between RNA modifications and biological processes. The gene set “h.all.v7.2” and “c2.cp.kegg.v7.1”were derived from the MSigDB database (Zhu et al., 2020). The functional annotation of 26 “writer” genes was conducted by the clusterProfiler R package, with a cutoff value of FDR <0.05. An adjusted P with value < 0.05 was considered as indicative of statistical significance.

CIBERSORT is a method that can accurately estimate the fraction of diverse cell subsets in gene expression profiles from complex tissues (http://cibersort.stanford.edu) (Newman et al., 2015). To predict the immune subset composition of HCC samples from gene expression profiles, CIBERSORT was used to estimate the relative abundance of 22 types of immune cells (model = absolute, permutation = 1,000, disable quantile normalization for RNA-Seq data as recommended).

Firstly, the RNA modification-related differentially expressed genes (DEGs) among distinct RNA modification clusters were collected using “limma” package of R software. Next, we performed univariate cox regression model to analyze the correlation of each gene with overall survival, and the significant prognosis DEGs were used for further analysis. Subsequently, distinct genomic subtypes were determined by unsupervised clustering analyses. In addition, the prognostic analysis was performed for each genomic subtype and extract principal component 1 and 2 as the signature scores. Finally, the RNA modification score was defined using a method similar to that used in analyzing gene-gene interactions (GGIs): WM-Score = Σ (PC1i + PC2i), which is defined as the expression of final RNA modification phenotype-related genes (Sotiriou et al., 2006; Zeng et al., 2019).

A Wilcoxon rank-sum test was utilized to compare differences between two groups, and Kruskal-Wallis test was used for comparisons of multiple groups. The discrimination accuracy of the WM-Score model was described by receiver operating characteristic (ROC) analysis. Kaplan-Meier method estimate curves were generated for prognostic analysis, and the differences between groups were evaluated by a log-rank test. Univariate and multivariate analyses were further carried out to assess independent risk factors. All data were analyzed by the R 4.0.1 software. A two-tailed p < 0.05 was considered as statistically significant.

A total of 26 RNA modification “writers” were selected in this study, which included seven m⁶A modification “writers”, three A-I modification “writers”, 12 APA modification “writers”, and four m¹A modification “writers” (Supplementary Table S1) (Li et al., 2016; Tang S. J. et al., 2020; Shen et al., 2021). To explore the genetic alterations in RNA modification writers, we examined the incidence of somatic mutations and somatic copy number variation (CNV) for all “writers” based on the TCGA database. Among 356 samples from TCGA-LIHC, 42 (11.8%) exhibited genetic changes of these writers, and the details was shown in the Figure 1A. The highest mutation frequency was presented in CPSF1, followed by ADARB2 and KIAA1429 (Figure 1A), while METTL3, METTL14, TRMT61A, TRMT61B, CSTF3, and NUDT2 did not show any mutations in tumor samples. Next, we used the hallmark gene set to perform gene set variation analysis (GAVA) to compare the mutation groups and those without mutation in “writers”. The GSVA indicated significantly enriched carcinogenic activation pathways in the mutation group, such as those of E2F targets, G2M checkpoint, MYC, and MTORC1 signaling pathway (Figure 1B). Furthermore, the investigation of CNV alteration in 26 regulators showed that ADAR, CPSF1, CPSF4, TRMT10C and KIAA1429 had a widespread frequency of CNV gain, while ZC3H13, CF1, METTL14, NUDT21, and WTAP had a significant CNV loss (Figure 1C). To explore whether the above CNV alterations affected the expression of the 26 RNA modification regulators, we compared the expression level of these regulators between tumor samples and paired normal samples. The results showed increased mRNA levels of most “writers” in tumor samples in comparison to normal samples (Figure 1D), suggesting that CNV might be the major factor leading to the aberrant expression of medication regulators. Notably, the mRNA levels of some “writers” were increased, while the frequencies of CNV loss for those were high. Therefore, further investigations were performed. According to the CNV value, patients were divided into 3 groups, including CNV amplification group, CNV deletion group, and normal group, and the mRNA expression of “writers” were compared between these groups (Figure 1E). The results showed mostly elevated expression for the group of patients with CNV amplification compared with the other groups with CNV deletion or normal CNV in these “writers”. Taken together, we mapped the genetic alterations of the 26 RNA modification “writers” between control tissues and tumor tissues, suggesting that these changes might play vital functions in HCC tumorigenesis and progression.

In order to further understand the role of RNA modification “writers” in HCC, we performed univariate analysis of the 26 regulators based on the TCGA-LIHC cohort. We found that 16 of 26 “writers” were markedly correlated with the OS of HCC patients (Figure 2A). Next, we explored the relationship among “writers” and found that most were positively or negatively correlated with each other (Figure 2B). Thus, it is suspected that the crosstalk between different “writers” may have a vital function in the different modification patterns of HCC.

We performed consensus clustering to classify patients into distinct RNA modification patterns based on the mRNA expression of “writers” (Supplementary Table S2). Eventually, two RNA modification patterns with 204 cases were determined in pattern 1 (cluster 1), and 113 cases in pattern 2 (cluster 2) (Figure 2C). Subsequently, “GSVA” enrichment analysis was employed to further understand the biological behaviors between the distinct two clusters. Our results indicated that cluster 1 was significantly enriched in metabolism and drug metabolism pathways, such as sulfur metabolism, primary bile acid biosynthesis, tyrosine metabolism, tryptophan metabolism, drug metabolism cytochrome P450, drug metabolism other enzymes, renin angiotensin system, while cluster 2 enrichment pathways were mainly linked to proliferation and signal transduction, including cell cycle, calcium conduction, etc. (Figure 2D).

Emerging evidence suggests that RNA modifications interact with the tumorigenic environment, thus affecting tumor occurrence, development, and prognosis (Jiang et al., 2020; Chen H. et al., 2021; Chong et al., 2021). Therefore, the function of the RNA methylations in the TME were further explored. The association analysis using the CIBERSORT method revealed that the identified RNA modification regulators might have close links with immune cell infiltration from the TME (Figure 3A). For instance, METTL14, ZC3H13, CSTF3, and ADAR were markedly negatively associated with Mo macrophage differentiation, while their positive association was observed with METTL3, RBM15B, KIAA1429, TRMT61A, TRMT6, CPSF1, and NUDT21. Moreover, we analyzed the difference in immune cell infiltration from TME between cluster 1 and cluster 2. The results revealed that the infiltration of M2, T cells, mast cells, and monocytes was higher in cluster 1. Notably, though, the infiltration of M1, regulatory T cells and follicular helper T cells was higher in cluster 2 (Figures 3B,C). Overall, cluster 2 was usually enriched in immunosuppressive cells, indicating a poor prognosis, whereas cluster 1 was characterized by immune cell activity, indicating a beneficial prognosis. These findings suggested that RNA modification “writers” play crucial roles in immune cell infiltration and TME formation.

Our results above demonstrated the important role of RNA modification in TME formation and patient prognosis, while these findings were based on RNA modification patterns and could not accurately evaluate the capacity of the RNA modification as a prognostic predictor in individual HCC patients. The underlying genetic alterations in these two RNA modification patterns were still unclear. Based on these queries, we examined the transcriptional expression change between the two patterns. A total of 273 DEGs related to RNA modification patterns were identified, and the further enrichment analysis showed that these DEGs were enriched in many essential biological processes, including DNA-binding transcription activator activity, signaling receptor activator activity, and multicellular organismal response to stress (Figures 4A–C). Subsequently, according to unsupervised clustering analysis based on the 273 DEGs, patients were classified into two stable transcriptomic subtypes: cluster A and cluster B (Figure 4D), with 242 and 75 of the 317 HCC patients, respectively. The prognosis of patients in cluster B was poorer than those in gene cluster A (Figure 4E; p < 0.0001, log-rank test).

Furthermore, we developed a score model based on the DEGs between gene clusters. As described in the Methods section, a scoring model named writers of RNA modification-score (WM-Score) was constructed. We discovered that cluster 1 had a higher WM-Score value than cluster 2 (Figure 4F). Consistently with this, cluster A also showed a higher score value than cluster B (Figure 4G). To evaluate the association of WM-Score value with TME, we further calculated the abundance of immune cell infiltration for the low and the high WM-Score value groups. We found that the infiltration rate of M0 macrophages, monocytes, and TfCD8 was higher in the high WM-Score value group, and that of activated NK-activated cells and M1 macrophages was higher in the low WM-Score value group (Figure 4H).

After confirming the efficacy of the WM-Score model in predicting patient prognosis, we investigated whether this model could be applied to determine the tumorigenesis, progression, invasion and metastasis of HCC. The prognostic efficiency of the scoring model was explored through classifying patients into low and high score groups the using “survminer” package. As expected, patients with high score demonstrated a poorer prognosis than those with low score in the TCGA-LIHC cohort (Figure 5A). We used ROC curve analysis to determine the discrimination accuracy of the scoring model in predicting patient prognosis. The area under the ROC curves (AUCs) of WM-Score values were 0.84, 0.76 and 0.79 at 1, 3 and 5 years overall survival, respectively (Figure 5B). Multivariate analysis for the TCGA-LIHC cohort also demonstrated that the WM-Score could serve as an independent prognostic predictor in HCC (Figure 5C). All of these results indicated that the WM-Score model has accurate prognostic value for HCC patients. The analysis of difference in WM-Scores between different TNM grades and clinical grades in the TCGA database indicated that samples with higher clinical grades and TNM stages usually have higher WM-Score values (Figures 5D,E). In addition, considering the EMT-related pathways, the samples with different WM-Score value had different pathway characteristics. For the TCGA database, samples with high WM-Score value were significantly related to cell cycle, DNA damage repair, and DNA replication, while samples with low WM-Score value were related to EMT, WNT target, and cell cycle regulators (Figure 5F).

For several years, sorafenib has been approved a treatment option for advanced HCC patients, while efficacy of sorafenib is limited by drug resistance (Gnoni et al., 2019). Aiming to further investigate whether the WM-Score value affected drug sensitivity, we evaluated the correlation between the scoring model and the drug response of tumor cell lines. Using Spearman’s correlation analysis, 15 significant correlation pairs were identified in the Cancer Drug Sensitivity Genomics (GDSC) database between scoring model and drug reaction (Yang et al., 2013). Among them, eight pairs of drug sensitivity were related to WM-Score value, and seven pairs showed resistance related to WM-Score value (Figure 6A). In addition, we also analyzed the signaling pathways of these drugs to determine target genes. We found that drugs associated with high WM-Score value mainly target KIT, CLAP, and cell cycle signaling pathways. In contrast, drugs related to low WM-Score value mostly target apoptosis regulation and cell cycle signaling pathways (Figure 6B). Taken together, these findings indicate that the WM-Score values are related to drug reaction, and thus might offer a framework to guide the treatment strategy of HCC.

In recent years, immune checkpoint inhibitors (ICIs) have made breakthroughs in the treatment of advanced HCC, while biomarkers that could effectively predict the efficacy of immunotherapy are still lacking. Herein, we explored whether the WM-Score model could predict therapeutic response to ICI therapy in HCC patients. For the IMvigor210 cohort, the therapeutic efficacy was significantly better in patients with low WM-Score value compared to those with high WM-Score value (Figure 6C). The frequency of response to anti-PD-1/L1 treatment in the low WM-Score value group was higher than that in the high WM-Score value group (Figure 6D). We also analyzed the survival difference of all samples of IMvigor210 and those under different stages. The results showed that total samples (Figure 6E), or Stage I + II samples (Figure 6F), or Stage III + IV samples (Figure 6G) all exhibited a marked difference in survival between samples with high and low WM-Score value. Especially in the prediction of high-stage clinical samples, the WM-Score value demonstrated extremely high power. Collectively, our results proved that the WM-Score model might serve as a potential predictor of response to anti-PD-1/L1 immunotherapy.