Skip to main content

Integrated investigation of the clinical implications and targeted landscape for RNA methylation modifications in hepatocellular carcinoma



RNA methylation (RM) is a crucial post-translational modification (PTM) that directs epigenetic regulation. It mostly consists of N1-methyladenosine (m1A), 5-methylcytosine (m5C), N3-methylcytidine (m3C), N6-methyladenosine (m6A), and 2′-O-methylation (Nm). The “writers” mainly act as intermediaries between these modifications and associated biological processes. However, little is known about the interactions and potential functions of these RM writers in hepatocellular carcinoma (HCC).


The expression properties and genetic alterations of 38 RM writers were assessed in HCC samples from five bioinformatic datasets. Two patterns associated with RM writers were identified using consensus clustering. Then, utilizing differentially expressed genes (DEGs) from different RM subtypes, we built a risk model called RM_Score. Additionally, we investigated the correlation of RM_Score with clinical characteristics, tumor microenvironment (TME) infiltration, molecular subtypes, therapeutic response, immunotherapy effectiveness, and competing endogenous RNA (ceRNA) network.


RM writers were correlated with TME cell infiltration and prognosis. Cluster_1/2 and gene.cluster_A/B were shown to be capable of distinguishing the HCC patients with poor prognosis after consensus and unsupervised clustering of RNA methylation writers. Additionally, we constructed RNA modification pattern-specific risk model and subdivided the cases into RM_Score high and RM_Score low subgroups. In individual cohorts or merged datasets, the high RM_Score was related to a worse overall survival of HCC patients. RM_Score also exhibited correlations with immune and proliferation related pathways. In response to anti-cancer treatments, the RM_Score had a negative correlation (drug sensitive) with drugs that focused on the MAPK/ERK and metabolism signaling, and a positive correlation (drug resistant) with compounds targeting RKT and PI3K/mTOR signaling pathway. Notably, the RM_Score was connected to the therapeutic effectiveness of PD-L1 blockage, implying that RM writers may be the target of immunotherapy to optimize clinical outcomes. Additionally, a ceRNA network was generated including 2 lncRNAs, 4 miRNAs, and 7 mRNAs that was connected to RM writers.


We thoroughly investigated the potential functions of RNA methylation writers and established an RM_patterns-based risk model for HCC patients. This study emphasized the critical functions of RM modification in TME infiltration, targeted therapy, and immunotherapy, providing potential targets for HCC.


Hepatocellular carcinoma (HCC), accounting for approximately 90% of liver cancer, is the major histologic type of liver cancer. Annually, more than a million individuals are newly diagnosed with liver cancer, raising concerns on the world’s health [1]. As an inflammation-associated malignancy, HCC is frequently caused by a number of risk factors, including drinking, nonalcoholic steatohepatitis (NASH), and hepatitis B virus (HBV)/hepatitis C virus (HCV) infection [2]. In contrast to other solid cancers, HCC is commonly diagnosed at an advanced stage, until which it is too late for transplantation, surgical resection, or local ablation. By decades, long-term improvements have been made in the management of advanced HCC. Sorafenib was first approved as an effective systemic option for advanced HCC. More recently, other targeted drugs such as lenvatinib, regorafenib, cabozantinib, nivolumab, and ramucirumab have also displayed clinical efficacy [3]. However, heterogeneity revealed by multiple omics poses a significant challenge for precision medicine and increases survival of HCC patients. Understanding the probable mechanisms of hepatocarcinogenesis, as well as searching for robust biomarkers and targets, will thus tremendously benefit HCC treatment.

Benefiting from high-throughput sequencing and large-scale profiling, post-translational modification (PTM) has been uncovered in a variety of RNAs, such as messenger RNA (mRNA), transfer RNA (tRNA), ribosomal RNA (rRNA) long noncoding RNA (lncRNA), and enhancer RNAs (eRNAs). RNA methylation (RM), including includes N1-methyladenosine (m1A), 5-methylcytosine (m5C), N3-methylcytidine (m3C), N6-methyladenosine (m6A), and 2′-O-methylation (Nm), is described as a crucial PTM in governing RNA maturation, splicing, stability, and translation. Among all these RNA methylation modifications, N6-methylation of adenosines is the most common and abundant RNA methylation modification that occurs at stop codons and within 3′ UTRs. The m1A mutation is likewise found at the first position of the adenine base. Recently, m1A is also reported enriched at translation start sites (5′ UTRs). Nm, which is found on the 2′ hydroxyl ribose moiety of ribonucleosides, has been found in all major eukaryotic RNA [4]. In contrast, m5C and m3C, two isomeric cytidine methylation modifications, present a methyl group on the nucleobase of cytidines. The newly identified RNA methylation type known as m3C differs from other well-known modification patterns. It was first discovered in rRNA and tRNA, in which it took part in modeling RNA structure and protein–RNA interactions [5, 6]. However, the m3C modification in mRNA of mice and human was also recently discovered [7].

Specific “writers” (methyltransferase), “readers” (binding proteins), and “erasers” (demethylases), in that order, can deposit, recognize, and remove labeled RNA methylation modifications, respectively. These enzymes have been connected to immune responses, cell differentiation, DNA replication, DNA repair, and other common biological processes. However, a growing number of studies have shown the critical functions of dysregulated modification modulators, particularly “writers”, in modulating the malignant phenotypes of various cancer types. For HCC, the oncogenic roles of m6A writers have drawn more attentions. Through the HuR-ETS1-p21/p27 axis, WTAP-induced m6A alteration enhanced the progression of HCC [8]. METTL3-mediated m6A modification contributed to sorafenib resistance of HCC cells in the hypoxic tumor microenvironment [9]. METTL14, on the other hand, impeded the invasion of HCC cells by suppressing EGFR/PI3K/AKT signaling [10]. Given the biological features of m6A writers, these inhibitors are promising drug candidates for tumor therapy, including FB23-2, R-2HG, BTYNB, and CS1 [11]. By recruiting H19 lncRNA to methylate, the typical m5C methylase NSUN2 accelerated carcinogenesis, cell migration, and invasion of HCC cells. 2′-O-methylation writer Fibrillarin (FBL) was correlated with advanced TNM stage and poor survival of HCC cases [12]. According to the most recent study, TRMT6/TRMT61A-mediated m1A methylation enhanced liver tumorigenesis and CSC development via activating Hedgehog and PPAR signaling [13]. More intriguingly, despite the fact that m3C has been specifically identified in the HBV-infected HCC cell line Huh7, the relationship between m3C writers and the evolution of HCC has not been fully understood [14]. It has been postulated that this methylation modification may play an important role in liver cancer brought on by viral infection.

Given the aforementioned, RNA methylation modification authors may be essential for the development of HCC. Though recent studies highlighted the significance of RNA methylation modification in tumorigenesis, the comprehensive analyses of the RNA methylation writers for HCC remain limited [15, 16]. In the current study, we investigated thoroughly the writers’ expression and mutation levels of 5 RNA methylation types (m1A, m5C, m3C, m6A, and Nm) in 5 HCC cohorts. Additionally, utilizing various clustering algorithms, we constructed an RNA methylation model to forecast tumor microenvironment (TME) alterations, immunotherapy response, and survival status. As a result, an integrated investigation of the patterns of RNA methylation may uncover the interaction regulatory mechanisms and provides potential HCC targets.


Data collection and processing

The study’s flowchart was depicted in Fig. 1 for reference. Gene expression levels and clinical information were downloaded from Gene-Expression Omnibus (GEO), the Cancer Genome Atlas (TCGA) database and ICGC database. mRNA expression, miRNA expression, lncRNA expression, somatic mutation, SCNAs were retrieved from TCGA database ( 2 GEO liver cancer cohorts (GSE76427 and GSE54236) [17, 18], 2 ICGC cohorts (LIRI-JP and LICA-FR), and TCGA-LIHC cohort were included for combined analysis (Additional file 2: Table S1). The log transformed FPKM values from TCGA, LIRI-JP, and LICA-FR were integrated with the gene array data (GSE76427 and GSE54236). Batch effects were corrected by using “ComBat” algorithm of svaPackage. IMvigor210 cohort, including 348 metastatic urinary tract transitional cell carcinoma cases administrated with atezolizumab, was included to evaluate the predictive value of RM_Score in predicting immunotherapy [19, 20]. The genomic, transcriptomic, and clinical data were obtained from

Fig. 1
figure 1

Flowchart of this study. The expression and mutation features of RNA methylated modification writers were investigated in 5 HCC datasets. Following clustering algorithms, a RNA methylation-based model was constructed to forecast tumor microenvironment alterations, immunotherapy response, targeted drug sensitivity, and survival status of HCC individuals

Mutation status and copy number variations of the RM writers

A total of 38 RM writers were involved in this analysis (Additional file 2: Table S2), including m6A writers (METTL3, METTL14, METTL16, ZCCHC4, WTAP, RBM15, RBM15B, ZC3H13, KIAA1429), m5C writers (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, TRDMT1), Nm writers (CMTR1, CMTR2, FBL, TRMT13, TRMT44, TARBP1, HENMT1, FTSJ1, FTSJ3, MRM1, FTSJ2, RNMTL1), m1A writers (TRMT6, TRMT61A, TRMT61B, TRMT10C, RRP8), and m3C writers (METTL2A, METTL2B, METTL6, METTL8). The “maftools” package was utilized to demonstrate the mutation levels of the RNA methylation writers in pan-cancers of TCGA dataset. The patients were split into groups with and without mutations, which were further compared for their overall survival. Copy number variations was calculated to evaluate the expression levels and copy number.

Clustering expression pattern of RNA methylation “writers”

Following excluding genes missed in any of the datasets (NSNU7 missed in LIRI-JP, RBM15B and MRM1 missed in LICA-FR), a total of 35 writers covering the 5 types of RNA methylated modification were included for further analyses, including 11 Nm writers (CMTR1, CMTR2, FBL, TRMT13, TRMT44, TARBP1, HENMT1, FTSJ1, FTSJ3, FTSJ2, RNMTL1), 8 m6A modification enzymes (METTL3, METTL14, METTL16, ZCCHC4, WTAP, RBM15, ZC3H13, KIAA1429), 7 m5C writers (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, TRDMT1), 4 m3C writers (METTL2A, METTL2B, METTL6, METTL8), and 5 m1A writers (TRMT6, TRMT61A, TRMT61B, TRMT10C, RRP8). Consensus clustering algorithm was performed to cluster analysis of 35 RM writers in combined liver cancer cohorts by performing “ConsensusClusterPlus” package with 100 repetitions.

Gene set variation analysis and TME analysis

By using the "GSVA" R package, gene set variation analysis (GSVA) was performed to investigate the impact of RNA modification clusters on biological processes and pathways [21]. “c2.cp.kegg.v7.1.symbols.gmt” and “h.all.v7.4.entrez.gmt.txt” for GSVA analysis was downloaded from the MSigDB database ( 35 RM genes were functionally annotated using the clusterProfiler R Package. he levels of infiltration of 28 different types of immune cells were measured using single sample gene set enrichment analysis (ssGSEA) in the "GSVA" R package. By analyzing 547 immune-cell signature genes, CIBERSORT ( was used to determine the relative abundance of 22 different immune cell types.

Construction of the RM_Score model

Using the Limma R package, DEGs associated to RNA methylation phenotype were identified between Cluster 1 and Cluster 2. Then we performed unsupervised clustering to obtain gene.cluster_A and gene.cluster_B. The prognostic value of DEGs was further determined using a univariate Cox regression model. Then, to establish a scoring system, DEGs related to survival were retrieved. RM_Score of each patient was calculated as follows: RMScore = \(\sum {PC}_{1}+\sum {PC}_{2}\). RMScore is equal to PC 1 plus PC 2. The optimal cut-off value for the RM_Score, which separated the cases into high RM_Score group and low RM_Score group, was defined using the surv_cutpoint function in the survminer R package.

Association analysis of RM_Score and drug sensitivity

The transcription profiles for about 1000 cancer cell lines, drug response measurements as AUC for antitumor drugs in cancer cell lines, and targets/pathways of drugs are downloaded from Genomics of Drug Sensitivity in Cancer (GDSC, [64]. We performed Spearman correlation analysis to calculate the correlation between drug sensitivity and RM_Score. |Rs|> 0.3 and FDR < 0.05 (Benjamini and Hochberg adjusted) were considered as significant correlation.

Construction of competing endogenous RNA (ceRNA) network

miRNA–mRNA interactions were retrieved from the miRTarBase(, TargetScan ( and miRDB ( miRNA–mRNA interactions overlapped in the three databases were utilized for further analysis. lncRNA–miRNA interactions were predicted by TargetScan. Cytoscape were performed to generate mRNA–miRNA–lncRNA network.


The data were analyzed using R (version 3.6.2) and R Bioconductor packages. Receiver operating characteristic (ROC) curve, Kaplan–Meier method and univariate/multivariate Cox regression model were used to verify the validity of the model. RCircos package was used to present the distribution of the RM writers in the chromosome. p value less than 0.05 was considered as statistically significant.


Genetic alterations of five types of RM writers in HCC

38 “writers” of the five most prevalent RNA methylation modifications were analyzed in this study, including 5 m1A writers, 8 m5C writers, 9 m6A writers, 4 m3C writers, and 12 Nm writers. 69 of the 377 liver cancer samples tested positive for the mutations of writers mentioned above. Among them, KIAA1429 displayed relatively higher mutation frequency (Fig. 2A). Despite the fact that there were no appreciable differences in survival between HCC patients with or without mutations (Additional file 1: Fig. S1A, B), the poor prognosis of HCC patients was substantially associated with mutations of the 20 RNA methylated modification writers (FTSJ2, FTSJ3, FTSJD1, HENMT1, KIAA1429, METTL2A, METTL2B, METTL8, METTL14, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, RBM15, TRDMT1, TRMT6, TRMT44, and TRMT61B; Additional file 1: Fig. S1C, P < 0.05). It was speculated that the writers' genetic alterations may have contributed to the development of HCC.

Fig. 2
figure 2

Comprehensive landscape of RNA methylated modification writers in liver cancer. A The mutation frequency of RNA modification writers in liver cancer patients from the TCGA cohort. B GSEA elucidating the samples with or without writer mutation. C The mutation landscape of the RNA modification writers among pan-cancers in TCGA cohort. D Different CNV types of RNA modification writers in the TCGA-LIHC cohort. E The expression of RNA modification writers in paired normal and cancer tissues in TCGA-LIHC cohort

To assess the differences between the individuals with and without mutation, Gene Set Variation Analysis (GSVA) was conducted using hallmark gene sets. As shown in Fig. 2B, the mutant group exhibited increased expression of cancer-related signals such as MYC targets, E2F targets, DNA repair, oxidative phosphorylation, and G2M checkpoint pathways. In contrast, TNFα signaling, inflammatory response, and interferon-γ response were downregulated. Next, we examined somatic copy number alterations of the RM writers. Of them, TARBP1, KIAA1429, MTEEL3, RBM15, RBM15B, and TARBP1 had extensive somatic mutation among pan-cancers (Fig. 2C). For HCC patients, the highest frequency of copy number variation (CNV) gain was observed in TARBP1 and KIAA1429, while the highest frequency of CNV loss was identified in RNMTL1 and METTL16 (Fig. 2D). Apart from the genetic variations, the majority of the writers had considerably higher levels of mRNA expression in HCC samples compared to that in normal tissues. In contrast, m5C writer NSUN6 displayed downregulated expression levels in HCC tissues (Fig. 2E). In order to comprehend more about the discrepancy between CNV and mRNA expression, we examined the subgroups of CNV levels, such as normal tissues, CNV gain, CNV loss, and/or non-significant CNV alterations (Additional file 1: Fig. S2). Indeed, patients with CNV gain showed frequently higher mRNA expression than those with CNV loss. However, ZC3H13, NSUN6, NSUN7, and RNMTL1 showed significant downregulation or non-significant alteration in CNV loss group compared to normal tissues. Thus, it suggested that CNV could partially modulate mRNA expression of the writers in HCC, indicating a high heterogeneity of genetic landscape and RM writer expression in HCC.

Distinct patterns of RM writers associated with cancer hallmarks and immune infiltration

As shown in Fig. 3A, the “writers” of the RNA methylation modifications exhibited distinct gene distributions. To comprehensively understand the expression pattern of the “writers” engaged in carcinogenesis, we selected five HCC datasets (TCGA, LIRI-JP, LICA-FR, GSE76427, and GSE54236) containing 955 HCC samples and 179 normal samples. Following excluding genes missed in any of the datasets, 5 writers covering 5 RM categories were recruited for further analysis. According to PCA analysis, these writers were able to effectively discriminate tumor samples and normal samples (Fig. 3B). Univariate Cox regression showed that 27 of 35 RM writers were substantially correlated with prognosis of HCC patients (Additional file 1: Fig. S3A). Additionally, consensus clustering showed that these writers had potentially negative or positive correlations (Fig. 3C).

Fig. 3
figure 3

Cancer hallmarks and immune infiltration of the RNA methylation writers. A The genomic distribution of the RNA methylation writers. B PCA of the RNA methylation writers for distinguishing tumor and normal patients in TCGA-LIHC cohort. C The regulatory network of the RNA methylation writers based on the consensus analysis. D Spearman correlation analysis showed positive or negative correlation among RNA methylation writers in HCC. E Consensus clustering of RNA methylation writers in 5 cohorts. F Kaplan–Meier curves demonstrating the overall survival of Cluster_1 and Cluster_2. G GSVA enrichment analysis shows the correlation of the two clusters with biological pathways and cancer hallmarks. H The difference of immune cell infiltration in TME between Cluster_1 and Cluster_2. I The expression difference M0 macrophage and CD4 T cell between Cluster_1 and Cluster_2

The spearman correlation analysis found that the majority of these writers had favorable correlations. However, NSUN6 and ZC3H13 had weak or negative associations with other writers (Fig. 3D). It was hypothesized that the writers of the various RNA methylation modification patterns might interact with one another. After consensus clustering, 955 patients from the pooled datasets were stratified into Cluster_1 and Cluster_2 (Fig. 3E; Additional file 2: Table S3). Remarkably, patients at Cluster 1 had better prognosis than that of Cluster 2 modification pattern (Fig. 3F, p < 0.0001). Subsequently, we conducted GSVA enrichment analysis to explore the biological implications of these different RM patterns (Additional file 2: Table S4, Table S5). According to the top 60 enriched pathways from KEGG analysis, DEGs of the two Clusters were mainly enriched in metabolism and pathways associated with proliferation, spliceosome and RNA modification. In the Hallmark analysis, DEGs were substantially associated with a number of immunological and oncogenic pathways, including Kras signaling, PI3K/AKT/mTOR signaling, IL2/IL6 signaling, interferon α/γ response, and TNFA signaling (Fig. 3G), suggesting the potential role of “writers” in immune and TME. Actually, there may be a specific association between TME cell infiltration and RM “writers” (Additional file 1: Fig. S3B). For example, NSUN6, METTL14, and ZC3H13 were significantly negatively correlated with T Cells CD4 memory activated, while FTSJ1, METTL6, NSUN4 and WTAP were positively correlated with such cells. Then, differences in TME cell infiltration were calculated between the two RM clusters (Additional file 2: Table S6). As observed in Fig. 3H, the infiltration of macrophages M0 (p = 0.027) and T cells CD4 memory resting (p = 9.03 × 10− 3) was higher in Cluster_1. In line with this, M0 macrophage and T cell CD4 memory resting marker gene expression was considerably elevated in Cluster_1 (Fig. 3I). It was hypothesized that RM patterns impacted the degree of infiltration by various immune cell types.

Construction of RM writer signature

To further characterize the functional role of the two RM patterns above, we identified 62 RM phenotype-related DEGs and performed enrichment analysis. We discovered that these genes were enriched in metabolic activities, including drug catabolic process, alpha-amino acid metabolic process, and steroid metabolism (Additional file 1: Fig. S4A), which was also implicated in pentose phosphate pathway, glycolysis, PPAR signaling pathway (Additional file 1: Fig. S4B). The patients were then categorized into two genomic subtypes known as gene.cluster_A and gene.cluster_B using unsupervised clustering analyses based on the DEGs (Additional file 1: Fig. 4A; Additional file 2: Table S7). As shown in the survival analysis in Fig. 4B, patients at the subgroup of gene.cluster_ B had a worse prognosis than those in gene.cluster_A (p < 0.0001). We then developed an RM_Score algorithm to characterize the RNA modification profile of individual HCC patients based on these RM-related DEGs. We found that RM_Score of Cluster_2 and gene.cluster_B were significantly higher than Cluster_1 and gene.cluster_A (Fig. 4C, D). To further assess the clinical relevance of the RM_Score, we divided patients into RM_Score-low and -high groups with the cut-off value of 0.05859662 determined by the “surv_cutpoint” algorithm of “survminer” package (Additional file 2: Table S8). We performed overlap analysis of these three different classifiers based on the Wayne diagram and the histogram of frequency distribution (Additional file 2: Table S9). As shown in Fig. 4E-H, these three computational techniques of categorization had a high degree of coincidence.

Fig. 4
figure 4

Construction of RM writers-related signature. A Unsupervised clustering of the RNA methylation-related genes. B Kaplan–Meier curves demonstrating the overall survival of gene.cluster_A and gene.cluster_B. C Comparison of RM_Score between Cluster_1 and Cluster_2. D Comparison of RM_Score between gene.cluster_A and gene.cluster_B. E–H The overlap and frequency of classifiers of RM_Score model with Cluster_1/2 and gene.cluster_A/B. ns, p > 0.05; *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Clinical characteristics associated with RM_Score in HCC

HCC can be divided into 3 molecular subtypes, named iCluster1–3 with distinct molecular features. Alluvial diagrams were plotted to display the association between different classifiers and the subtypes (Fig. 5A, Additional file 2: Table S10). To examine the association between the RM_Score and subtypes, we then calculated the RM_Scores of different subtypes. As shown in Fig. 5B, RM_Scores dramatically varied among subtypes, with iCluster 1 achieving the highest score. Additionally, more advanced stages and grades of HCC cases were associated with higher RM_Score (Fig. 5C, D), implying potential relevance for RM patterns with tumor progression. As elucidated in Fig. 5E-H, patients with higher RM_Score experienced worse survival in the pooled cohorts or single cohort of the TCGA, ICGC, and GSE54236. At 3, 6, and 12 months of overall survival, the AUCs of ROC curves based on RM_Score were 0.63, 0.72, and 0.71, respectively (Fig. 5I). Then RM_Score was examined using multivariate Cox regression analysis as an independent prognostic factor. In the combined cohorts, RM_Score (Fig. 5J; HR = 2.78, 95% CI 2.04–3.8, p < 0.001), together with stage, were independent prognostic biomarkers. Consistently, RM_Score (Fig. 5K; HR = 2.8, 95% CI 1.87–4.1, p < 0.001) and stage were also defined as independent markers in TCGA cohort.

Fig. 5
figure 5

The clinical implications and prognostic value of the RM_Score system. A The correlation of molecular subtypes with Cluster_1/2, gene.cluster_A/B, and RM_Score high/low. B-D The comparison of RM_Score in patients at different iClusters, stages, and grades. E Kaplan–Meier curves demonstrating the overall survival of RM_Score high and RM_Score low subgroups in combined cohorts. F–H Kaplan–Meier curves demonstrating the overall survival of RM_Score high and RM_Score low subgroups in TCGA, LIRI, and GSE54236 cohorts, respectively. I The ROC curves of the RM_Score in predicting the overall survival of the liver cancer patients. J, K the Multivariate Cox regression model analysis to evaluated RM_Score as an independent prognostic factor in combined cohorts and TCGA-LIHC cohort. ns, p > 0.05; *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Pathways and immune activity with RM_Score

Then, we initially compared the two RM_Score subgroups-related pathways in the combined cohorts. DNA replication, cell cycle, immune checkpoints, DNA repair, the Wnt pathway, and DNA replication were enriched in the RM_Score-high group, whereas angiogenesis and antigen processing were involved in the RM_Score-low group (Fig. 6A). In TCGA cohort, the DEGs between RM_Score high and low groups were also predominantly enriched in proliferation- and immune-related pathways (Fig. 6B). Furthermore, we also explored into whether the two groups differed in terms of mutation. Although the two groups had identical altered genes, as illustrated in Fig. 6C, the mutation frequency and type significantly varied, particularly for TP53, TTN, CTNNB1, and MUC4/16. Additionally, we identified a connection between RM_Score and immune cells in TME. Neutrophil and activated CD4 T cell proportions showed the most notable differences in cell proportion (Additional file 1: Fig. S5). The immune score and estimate score were both positively correlated with RM_Score, whereas tumor purity was considerably adversely correlated with this score (Fig. 6D; Additional file 2: Table S11). Tracking Tumor Immunophenotype (TIP) discovered an obvious distinction between RM_Score high and low groups (Fig. 6E). The RM_Score-high group was triggered in Step1.Releasing of cancer cell antigens, Step4.Th22 cell.recruiting, Step4.Neutrophil.recruiting, Step4.Th17 cell.recruiting, and Step4.MDSC.recruiting. Contrarily, Step6.Recognition of cancer cells by T cells, Step7. Killing of cancer cells were implicated in RM_Score-low group.

Fig. 6
figure 6

Potential characteristics underlying the RM_Score. A, B The differences of enriched signaling pathways between RM_Score-high and -low subgroups in combined datasets and the TCGA-LIHC dataset. C The mutational genes with top frequency in RM_Score high and low subgroups. D The correlation of RM_Score with STROMAL score, Immune score, Estimate score, and tumor purity. E the Tracking Tumor Immunophenotype analysis of RM_Score high and low subgroups

The RM_Score model predicts response to immunotherapy

Developing biomarkers to forecast immunotherapy response has been a top priority. The expression of the checkpoint genes was then discovered in the two groups. CTLA4, PDCD1 and SLAMF7 were significantly elevated in high RM_Score group (Fig. 7A). Considering the potential association between RM_Score and immunological milieu, we examined the predictive capacity of the RM_Score towards ICB therapy. In the IMvigor210 cohort, we discovered that patients with low RM_Score (cut-off value = 0.002636961) showed substantially therapeutic benefits and had a noticeably extended overall survival (Fig. 7B, p = 0.00022). Additionally, the 348 individuals in the IMvigor210 cohort responded to anti-PD-L1 blockers with varying degrees (Additional file 2: Table S12). As shown in Fig. 7C, patients with high RM_Scores exhibited larger ratios of stable disease (SD) and progressive disease (PD) than patients with low RM_Scores. Additionally, the patients with complete response (CR) and partial response (PR) showed reduced RM_Score than SD and PD patients (Fig. 7D; Additional file 2: Table S13). Submap analysis indicated that the low RM_Score group was more inclined to respond to anti-PD-1 treatment (Fig. 7E). We further calculated the RM_Score of the three immune subtypes of IMvigor210, including “immune inflamed”, “immune excluded”, and “immune desert”. The results revealed that the “immune inflamed” type had lower RM_Score than two other groups (Fig. 7F). Likewise, the TMB and neoantigen burden were considerably higher in the group with a low RM_Score than in the group with a high RM_Score (Fig. 7G, H), indicating a potential correlation between the RM_Score model and immunotherapy.

Fig. 7
figure 7

The capacity of RM_Score model in predicting immunotherapy response. A The expression of immune-related genes in RM_Score-high and -low subgroups. B Kaplan–Meier curves show overall survival in the RM_Score-high and -low subgroups after the PD-L1 immunotherapy in the IMvigor210 cohort. C The proportion of patients in the IMvigor210 cohort with different responses to PD-L1 blockade immunotherapy. D The difference of clinical outcomes with anti-PD-L1 treatment in the RM_Score high and low subgroups in the IMvigor210 cohort. E The similarity of gene expression profiles between RM patterns and melanoma patients treated with ICB (n = 47). F The RM_Score in three immune subtypes of IMvigor210, including “immune inflamed”, “immune excluded”, and “immune desert”. G, H The TMB and neoantigen burden in the RM_Score-low and RM_Score subgroups. ns, p > 0.05; *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Potential therapeutic value of the RM_Score and ceRNA network development

Given the availability of targeted therapy for HCC, we assessed sorafenib sensitivity in both groups. To ensure the stability of the prediction results, prediction models were trained using the Genomics of Drug Sensitivity in Cancer (GDSC) cell line dataset by ridge regression and validated using tenfold crossover. For each sample, subsequently, IC50 values were calculated and the discrepancies were compared. We identified that the high score group in the discovery cohort was more likely to be responsive to sorafenib (Fig. 8A, p = 3.2 × 10−7). In the GDSC database, we found 28 strongly associated pairs between RM_Score and drug sensitivity using the Spearman correlation analysis (Fig. 8B). Among them, 9 pairs showed that drug resistance correlated with the RM_Score, e.g., FLT3 inhibitor (quizartinib, Rs = 0.46, p = 2.28 × 10−20), BCL2 inhibitor (venetoclax, Rs = 0.37, p = 3.41 × 10−13), mTOR inhibitor (OSI-207, Rs = 0.41, p = 3.94 × 10−16; PI-103, Rs = 0.31, p = 8.0 × 10−10; rapamycin, Rs = 0.35, p = 4.86 × 10−12) and ALK inhibitor (LDN-193189. Rs = 0.39, p = 1.24 × 10−14). 19 pairs showed drug sensitivity, e.g., MEK inhibitor (trametinib, Rs =  − 0.59, p = 2.57 × 10− 35; refametinib, Rs =  − 0.44, p = 5.18 × 10− 19; PD0325901, Rs =  − 0.47, p = 5.70 × 10− 22; CI-1040, Rs =  − 0.47, p = 1.37 × 10− 21) and ERK inhibitor (SCH772984, Rs =  − 0.54, p = 7.44 × 10− 30; ulixertinib, Rs =  − 0.48, p = 1.26 × 10− 22). Further, we analyzed the signaling pathways of the genes targeted by these drugs. We discovered that compounds that were sensitive to high RM_Score primarily targeted the signaling pathways for MAPK/ERK and metabolism. In contrast, low RM_Score associated drugs targeted RKT and PI3K/mTOR signaling pathway (Fig. 8C). Thus, it indicated that RNA methylation patterns are correlated with drug sensitivity, which provides potential treatment strategies. In order to investigate the RM_Score-related potential key axis, we further constructed the ceRNA network based on the identification and validation of the clinical implications. By mincing the TCGA dataset, 735 differentially expressed mRNA (|logFC|> 1, p < 0.05), 165 differentially expressed miRNA (|logFC|> 1, p < 0.05) and 302 differentially expressed lncRNA (|logFC|> 1, p < 0.05) were identified between high RM_Score and low RM_Score samples (Fig. 8D). By utilizing Cytoscape, 2 DElncRNAs, 4 DEmiRNAs and 7 DEmRNAs were used to construct the ceRNA network (Fig. 8E, Additional file 2: Table S14).

Fig. 8
figure 8

RM_Score-related Therapeutic sensitivity and ceRNA network. A IC50 values of sorafenib in RM_Score high and low groups in Genomics of Drug Sensitivity in Cancer (GDSC) cell line dataset. B The Spearman correlation analysis discovered the association between the RM_Score and drug sensitivity in the GDSC database. C The signaling pathways targeted by these RM_Score sensitive drugs. D The volcano map and heatmaps elucidating the differentially expressed mRNA, miRNA, lncRNA between RM_Score high cases and RM_Score low cases E The ceRNA network consisting of 2 DElncRNAs, 4 DEmiRNAs and 7 DEmRNAs


Methylation modification is a common procedure that occurs in practically all biological processes. Methylation modifications on DNA and proteins have been thoroughly investigated during the past few decades [22, 23]. The characteristics and mechanisms of various RNA methylation modifications on eukaryotic RNA, however, are still largely unknown. The involvement of RNA methylation modifications dysregulation in several diseases, particularly malignancies, are being highlighted by a growing body of studies [24]. There have been several “writers” of RNA methylation modifications that have been associated with the onset or progression of liver cancer [25, 26]. A thorough study of the characteristics of RNA methylation modifications might enable to improve survival prediction and uncover potential molecular targets for HCC. In the current study, for the first time, we provided a comprehensive overview of five RNA methylation alteration types in HCC, with an emphasis on the writers. According to the mutation data of TCGA, 20 mutated writers were correlated with poor survival of HCC patients, of which KIAA1429 and FTSJD2 had highest mutation frequency. DEGs between HCC patients with and without writer mutation were mainly enriched in proliferation and inflammatory pathways. While two m5C writers, NSUN6 and TRDMT1, were downregulated in HCC tissues, the majority of writers were overexpressed in HCC samples. Subsequent CNV analysis revealed significant genetic heterogeneity and RNA modification "writer" expression discrepancies between the normal and HCC samples.

Following the exclusion of missing genes, 35 writers were recruited to investigate the expression pattern in five HCC datasets. Positive interactions were discovered among the authors, indicating a cross-talk of the various RNA methylation modification types. HCC individuals were categorized into two clusters based on consensus clustering, with distinct overall survival. Additionally, Cluster_1 was enriched in RNA processing pathways, whereas Cluster_2 was implicated in metabolic pathways. Next, we examined the probable functions of “writers” in TME. Diverse immune cell types were infiltrated into two clusters, indicating that RNA modification was involved in the infiltration of particular immune cell types.

Then, using unsupervised clustering, we established two genomic subgroups using RNA methylation modification-related DEGs. A poorer prognosis was attributed to gene.cluster B. Additionally, we generated an RNA methylation modification DEGs-based model to score HCC individuals in the pooled cohort, which was then divided into two groups based on the RM_Score. This stratification using the RM_Score was particularly related to the other two computational methods. The relationship between the RM_Score and clinical parameters was then assessed. Molecular subtype iCluster 1 had a higher RM_Score than iCluster 2/3. According to a previous study, iCluster 1 had a higher tumor grade, more macrovascular invasions, and a considerably worse prognosis than the other two clusters [27]. RM_Score is consistently elevated in HCC cases compared with higher grades or stages. Moreover, ROC and multivariate Cox regression analysis showed that the RM_Score model had remarkable prediction performance for the survival of HCC patients. In addition, patients with high RM_Score have significantly shorter overall in pooled or either dataset. It was hypothesized that the high RM_Score might be related to the clinical prognosis and HCC features. P53, widely acknowledged as a tumor suppressor, is inactivated with mutation, which promotes tumor progression [28]. It is interesting to note that the prevalence of TP53 mutations in the RM_Score-high group was more than twice as high as that in the RM_Score-low group, pointing to a possible link between aggressive traits related to RM_Score and TP53 mutation.

On the basis of the enrichment analysis, we further investigated the relationship between RM_Score and tumor immunity. Although there was no discernible difference in the proportion of immune cells between the groups, RM_Score was favorably associated with immune score, estimate score, and tumor immunophenotypes including Th22/Th17 cell recruitment, release of cancer cell antigens, and neutrophil/MDSC recruitment. It also implied a connection between RM_Score and tumor immune. Nowadays, immunotherapy is gaining popularity as a promising treatment option for advanced HCC with immune escape [29, 30]. In light of this, we assessed the RM_Score model's accuracy in predicting immunotherapy response. Patients with low RM_Scores may benefit from ICB treatment and have significantly longer overall survival. Patients with high RM_Score, on the other hand, showed less response after therapy, indicating a correlation between RM_Score and ICB resistance. In addition, RM_Score was varied in different levels of TMB and neoantigen burden, which were frequently correlated with immunotherapy efficacy. Thus, the RM_Score is a candidate model for predicting ICB response.

Tyrosine kinase inhibitor (TKI) sorafenib was approved by the USA FDA in 2007 for the treatment of advanced HCC. Increasing oral multi-targeted TKIs were subsequently awarded permission for use in treating HCC worldwide [31, 32]. Furthermore, we evaluated the relationship between RM_Score and Sorafenib sensitivity. Interestingly, the group with a high RM_Score demonstrated enhanced sorafenib sensitivity. The administration of sorafenib might be advantageous for those who experienced ICB resistance. For other targeted drugs, according to the analyses in GDSC database, high RM_Score drug indicated resistance to BCL2 inhibitor, JAK inhibitor, mTOR inhibitor and ALK inhibitor, while low RM_Score showed drug sensitivity of MEK inhibitor and ERK inhibitor. Thus, it indicated that RNA methylation patterns are correlated with sensitivity of targeted drug, providing potential treatment strategies for HCC patients. Additionally, a ceRNA network associated with RM_Score was developed to further explore clinical implications or malignant behaviors. It should be noted that LINC00998 and DHRS4-AS1 were found to be hub lncRNAs that were associated with miRNAs and mRNAs. Multiple tumor types have been found to be suppressed by LINC00998 [33, 34]. However, LINC00998-encoded peptide accelerated HCC tumorigenesis and facilitated aggressive behaviors [35]. DHRS4-AS1 was also defined as an anti-oncogenic ncRNA, which could inhibit proliferation of HCC cells [36].

As previously stated, the findings presented here demonstrate that these RM networks enable connections between clinical and RNA modification by stratifying patients' prognosis and therapy responses. However, this study had some limitations as well. Though the promising conclusions were based on the integrative bioinformatics in multiple levels, future functional and mechanistic research on these RNA writers will reveal clinical phenotypes that are caused by RNA methylation writers. To fully illustrate the network of RNA methylation modifications, future studies should also include “readers” and “erasers” and other RNA methylation like m7G and A-I modification. In addition, larger and multi-center clinical cohorts should be employed to assess the accuracy of RNA methylation-based classifier and the RM_Score prognostic model.


In conclusion, we have demonstrated that RNA methylation modification patterns may serve crucial role in the progression of liver malignancy and immune dysfunction. Additionally, the subgroup classification based on RNA methylation writers may distinguish individuals with poor survival and forecast the efficacy of immunotherapy or targeted therapy. Moreover, RNA methylation modification regulators are also highlighted as robust biomarkers or potential targets for HCC.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.



Hepatocellular carcinoma


Hepatitis B virus


Hepatitis C virus


The Cancer Genome Atlas Liver Hepatocellular Carcinoma


National Center for Biotechnology Information


Gene Expression Omnibus


RNA methylation


Post-translational modification












Differentially expressed genes


Competing endogenous RNA


Gene set variation analysis


Copy number variation


Genomics of Drug Sensitivity in Cancer


  1. Vogel A, Meyer T, Sapisochin G, Salem R, Saborowski A. Hepatocellular carcinoma. Lancet. 2022;400(10360):1345–62.

    Article  CAS  Google Scholar 

  2. Lee YT, Fujiwara N, Yang JD, Hoshida Y. Risk stratification and early detection biomarkers for precision HCC screening. Hepatology. 2022.

    Article  Google Scholar 

  3. Yang C, Zhang H, Zhang L, Zhu AX, Bernards R, Qin W, et al. Evolving therapeutic landscape of advanced hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2022.

    Article  Google Scholar 

  4. Cantara WA, Crain PF, Rozenski J, McCloskey JA, Harris KA, Zhang X, et al. The RNA Modification Database, RNAMDB: 2011 update. Nucleic Acids Res. 2011;39(Suppl_1):D195-201.

    Article  CAS  Google Scholar 

  5. Liu F, He C. A new modification for mammalian messenger RNA. J Biol Chem. 2017;292(35):14704–5.

    Article  CAS  Google Scholar 

  6. Laptev I, Shvetsova E, Levitskii S, Serebryakova M, Rubtsova M, Zgoda V, et al. METTL15 interacts with the assembly intermediate of murine mitochondrial small ribosomal subunit to form m4C840 12S rRNA residue. Nucleic Acids Res. 2020;48(14):8022–34.

    Article  CAS  Google Scholar 

  7. Xu L, Liu X, Sheng N, Oo KS, Liang J, Chionh YH, et al. Three distinct 3-methylcytidine (m(3)C) methyltransferases modify tRNA and mRNA in mice and humans. J Biol Chem. 2017;292(35):14695–703.

    Article  CAS  Google Scholar 

  8. Chen Y, Peng C, Chen J, Chen D, Yang B, He B, et al. WTAP facilitates progression of hepatocellular carcinoma via m6A-HuR-dependent epigenetic silencing of ETS1. Mol Cancer. 2019;18(1):127.

    Article  Google Scholar 

  9. Lin Z, Niu Y, Wan A, Chen D, Liang H, Chen X, et al. RNA m(6) A methylation regulates sorafenib resistance in liver cancer through FOXO3-mediated autophagy. EMBO J. 2020;39(12): e103181.

    Article  CAS  Google Scholar 

  10. Shi Y, Zhuang Y, Zhang J, Chen M, Wu S. METTL14 inhibits hepatocellular carcinoma metastasis through regulating EGFR/PI3K/AKT signaling pathway in an m6A-dependent manner. Cancer Manag Res. 2020;12:13173–84.

    Article  CAS  Google Scholar 

  11. Li W, Hao Y, Zhang X, Xu S, Pang D. Targeting RNA N(6)-methyladenosine modification: a precise weapon in overcoming tumor immune escape. Mol Cancer. 2022;21(1):176.

    Article  CAS  Google Scholar 

  12. Zhang J, Yang G, Li Q, Xie F. Increased fibrillarin expression is associated with tumor progression and an unfavorable prognosis in hepatocellular carcinoma. Oncol Lett. 2021;21(2):92.

    Article  CAS  Google Scholar 

  13. Wang Y, Wang J, Li X, Xiong X, Wang J, Zhou Z, et al. N(1)-methyladenosine methylation in tRNA drives liver tumourigenesis by regulating cholesterol metabolism. Nat Commun. 2021;12(1):6314.

    Article  CAS  Google Scholar 

  14. McIntyre W, Netzband R, Bonenfant G, Biegel JM, Miller C, Fuchs G, et al. Positive-sense RNA viruses reveal the complexity and dynamics of the cellular and viral epitranscriptomes during infection. Nucleic Acids Res. 2018;46(11):5776–91.

    Article  CAS  Google Scholar 

  15. Zhao K, Li W, Yang Y, Hu X, Dai Y, Huang M, et al. Comprehensive analysis of m(6)A/m(5)C/m(1)A-related gene expression, immune infiltration, and sensitivity of antineoplastic drugs in glioma. Front Immunol. 2022;13: 955848.

    Article  CAS  Google Scholar 

  16. Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, et al. Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer. 2021;20(1):29.

    Article  CAS  Google Scholar 

  17. Grinchuk OV, Yenamandra SP, Iyer R, Singh M, Lee HK, Lim KH, et al. Tumor-adjacent tissue co-expression profile analysis reveals pro-oncogenic ribosomal gene signature for prognosis of resectable hepatocellular carcinoma. Mol Oncol. 2018;12(1):89–113.

    Article  CAS  Google Scholar 

  18. Villa E, Critelli R, Lei B, Marzocchi G, Camma C, Giannelli G, et al. Neoangiogenesis-related genes are hallmarks of fast-growing hepatocellular carcinomas and worst survival. Results from a prospective study. Gut. 2016;65(5):861–9.

    Article  CAS  Google Scholar 

  19. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.

    Article  CAS  Google Scholar 

  20. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, et al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet. 2017;389(10064):67–76.

    Article  CAS  Google Scholar 

  21. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  Google Scholar 

  22. Murakami S, Jaffrey SR. Hidden codes in mRNA: Control of gene expression by m(6)A. Mol Cell. 2022;82(12):2236–51.

    Article  CAS  Google Scholar 

  23. Boulias K, Greer EL. Means, mechanisms and consequences of adenine methylation in DNA. Nat Rev Genet. 2022;23(7):411–28.

    Article  CAS  Google Scholar 

  24. An Y, Duan H. The role of m6A RNA methylation in cancer metabolism. Mol Cancer. 2022;21(1):14.

    Article  CAS  Google Scholar 

  25. Liu L, Gu M, Ma J, Wang Y, Li M, Wang H, et al. CircGPR137B/miR-4739/FTO feedback loop suppresses tumorigenesis and metastasis of hepatocellular carcinoma. Mol Cancer. 2022;21(1):149.

    Article  CAS  Google Scholar 

  26. Dai YZ, Liu YD, Li J, Chen MT, Huang M, Wang F, et al. METTL16 promotes hepatocellular carcinoma progression through downregulating RAB11B-AS1 in an m(6)A-dependent manner. Cell Mol Biol Lett. 2022;27(1):41.

    Article  CAS  Google Scholar 

  27. Cancer Genome Atlas Research Network. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017;169(7):1327-1341 e1323.

    Article  Google Scholar 

  28. Zhu M, Wu M, Bian S, Song Q, Xiao M, Huang H, et al. DNA primase subunit 1 deteriorated progression of hepatocellular carcinoma by activating AKT/mTOR signaling and UBE2C-mediated P53 ubiquitination. Cell Biosci. 2021;11(1):42.

    Article  CAS  Google Scholar 

  29. Cheng K, Cai N, Zhu J, Yang X, Liang H, Zhang W. Tumor-associated macrophages in liver cancer: from mechanisms to therapy. Cancer Commun (Lond). 2022.

    Article  Google Scholar 

  30. Donne R, Lujambio A. The liver cancer immune microenvironment: therapeutic implications for hepatocellular carcinoma. Hepatology. 2022.

    Article  Google Scholar 

  31. Huang A, Yang XR, Chung WY, Dennison AR, Zhou J. Targeted therapy for hepatocellular carcinoma. Signal Transduct Target Ther. 2020;5(1):146.

    Article  Google Scholar 

  32. Sun Y, Zhang W, Bi X, Yang Z, Tang Y, Jiang L, et al. Systemic therapy for hepatocellular carcinoma: Chinese consensus-based interdisciplinary expert statements. Liver Cancer. 2022;11(3):192–208.

    Article  CAS  Google Scholar 

  33. Cai H, Yu Y, Ni X, Li C, Hu Y, Wang J, et al. LncRNA LINC00998 inhibits the malignant glioma phenotype via the CBX3-mediated c-Met/Akt/mTOR axis. Cell Death Dis. 2020;11(12):1032.

    Article  CAS  Google Scholar 

  34. Fang X, Pan X, Mai H, Yuan X, Liu S, Wen F. LINC00998 functions as a novel tumor suppressor in acute myeloid leukemia via regulating the ZFP36 ring finger protein/mammalian target of rapamycin complex 2 axis. Bioengineered. 2021.

    Article  Google Scholar 

  35. Pang Y, Liu Z, Han H, Wang B, Li W, Mao C, et al. Peptide SMIM30 promotes HCC development by inducing SRC/YES1 membrane anchoring and MAPK pathway activation. J Hepatol. 2020;73(5):1155–69.

    Article  CAS  Google Scholar 

  36. Zhou Y, Li K, Zou X, Hua Z, Wang H, Bian W, et al. LncRNA DHRS4-AS1 ameliorates hepatocellular carcinoma by suppressing proliferation and promoting apoptosis via miR-522–3p/SOCS5 axis. Bioengineered. 2021.

    Article  Google Scholar 

Download references


Not applicable.


This study was supported by grants from National Natural Science Foundation (82070622, 82272839, 81702419), the Key Research and Development Plan of Jiangsu Province (BE2020668), Science and Technology Project of Nantong City (JCZ21041).

Author information

Authors and Affiliations



WZ, JZ and CC designed and supervised this study. JZ, MH, and JG conducted the analyses and drafted of the manuscript. SX performed the statistics. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Chun Cheng, Wenjie Zheng or Jie Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1

. The survival status in patients with or without RM writer mutation. A, B The OS and DFS in HCC patients at mutation type and non-mutation type in TCGA-LIHC cohort. C1-C38, The OS in HCC patients with or without RM writers mutations in TCGA-LIHC cohort. Figure S2. The distribution of correlation coefficient between writer expression and CNV in HCC. The mRNA expression of the 38 RM writers in Normal, CNV_loss, None_CNV, and CNV_gain groups. Figure S3. The prognostic analysis of RM writers and correlation with TME cells. A The Univariate cox analysis to evaluate the correlation of RM writers with overall survival of HCC patients. B Heatmap showed the positive (red) and the negative (blue) correlation between TME infiltration and RM writers in HCC. Figure S4. The biological functions and pathways underlying the RM phenotype-related DEGs. A GO enrichment of the 62 RM phenotype-related DEGs. B KEGG enrichment analysis of the 62 RM phenotype-related DEGs. Figure S5. The correlation of RM_Score with TME infiltration. Heatmap shows the differences in TME infiltration between RM_Score-high and -low groups in the combined cohorts.

Additional file 2: Table S1.

Information of HCC cohorts. Table S2. RNA Modification Writers. Table S3. Samples clustering in 5 HCC cohorts. Table S4. GSVA_KEGG analysis of the DEGs between the Cluster 1/2. Table S5. GSVA_HALLMARK analysis of the DEGs between the Cluster1/2. Table S6.  Difference of TME infiltration characteristics between Cluster_1 /2. Table S7. Gene clustering in the pooled HCC cohorts. Table S8. RM_ScoreGroup stratification of the pooled cohort. Table S9. Sample information of the three computational methods of classification in 5 HCC cohorts. Table S10. Association of subtype with three computational methods of classification in TCGA LIHC cohorts. Table S11. The immune-related score of the samples in 5 HCC cohort. Table S12. Clinical information of IMvigor210 cohort. Table S13. The RM_Score of the samples in IMvigor210 cohort. Table S14. The RM_Score related ceRNA network.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, J., Gao, J., Hu, M. et al. Integrated investigation of the clinical implications and targeted landscape for RNA methylation modifications in hepatocellular carcinoma. Eur J Med Res 28, 46 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: