Skip to main content
Download PDF

Development and validation of a novel prognostic signature based on m6A/m5C/m1A-related genes in hepatocellular carcinoma

BMC Medical Genomics volume 16, Article number: 177 (2023) Cite this article

Abstract

Background

RNA methylation modification plays an important role in cancers. This study sought to examine the association between m6A/m5C/m1A-related genes and hepatocellular carcinoma (HCC).

Methods

Gene expression and clinical data of HCC patients were obtained from the TCGA database. Unsupervised consensus clustering was performed according to the expression of m6A/m5C/m1A-related genes in HCC. The relationships among prognosis, clinicopathological features and molecular subtypes were analyzed. Least absolute shrinkage and selection operator (LASSO) regression analysis was used to establish the m6A/m5C/m1A-related gene prognostic signature. Furthermore, the prognostic signature was validated based on the ICGC dataset. RT‒qPCR was used to detect the expression of the model genes in HCC. Clinicopathological features, functional enrichment, gene mutations, immune cell infiltration, and immunotherapy response in different risk groups were analyzed. A nomogram based on risk score and stage was constructed to predict HCC patient prognosis.

Results

Two m6A/m5C/m1A-related molecular subtypes were identified in HCC, and the prognosis of cluster C1 was worse than that of cluster C2 (p < 0.001). Highly expressed genes in cluster C1 are significantly correlated with G3-4, T3-4, stage III-IV (p < 0.05). An m6A/m5C/m1A-related prognostic signature was established and validated. The RT‒qPCR results showed that the risk signature genes were significantly upregulated in liver cancer tissue (p < 0.05). The prognosis of HCC patients in the high-risk group was worse than that of those in the low-risk group (p < 0.05). Multivariate Cox analysis indicated that the risk score was an independent factor predicting prognosis in HCC patients. ssGSEA revealed that the risk score correlated with the tumor immune microenvironment in HCC. Gene mutation analysis showed that the tumor mutation burden of patients in the high-risk group was much higher (p < 0.05), and the prognosis of HCC patients with high risk scores and high mutation burden was the worst (p = 0.007). A nomogram combining risk scores with clinicopathological features showed performed well in predicting HCC prognosis.

Conclusions

The m6A/m5C/m1A-related genes could predict the prognosis and tumor microenvironment features of HCC and can be important biomarkers relevant to the immunotherapy response.

Peer Review reports

Background

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors in the world. In 2020, there were approximately 906,000 newly diagnosed cases of HCC worldwide and 830,000 deaths from HCC, making it the sixth most common cancer and the third leading cause of cancer-related death [1, 2]. HCC has insidious, rapid onset and an extremely high degree of malignancy [3]. The 5-year survival rate is less than 18% [1, 4], which seriously affects public health. Considering the limitations of HCC treatment, new therapeutic targets are needed to improve the prognosis of HCC patients. Therefore, it is urgent to find a prognosis-related diagnostic model that could provide new directions for developing feasible targeted therapy approaches and improving the survival and prognosis of patients.

Epigenetic modifications include chemical modifications of DNA, RNA and proteins characterized by altered gene expression and function without any changes in the gene sequence. In addition to well-established DNA and protein epigenetic modifications, reversible RNA methylation has led to the third wave of studies in the epigenetic field. The main forms of RNA methylation are N6-methyladenosine (m6A), 5-methylcytosine (m5C) and N1-methyladenosine (m1A). m6A is the most abundant internal RNA modification in eukaryotic cells.

A large number of scholars have previously conducted research regarding the mechanism of m6A/m5C/m1A-related genes in HCC. LY6/PLAUR domain 1 (LYPD1) can promote tumorigenesis, ALKB Homolog 5 (ALKBH5) mediated m6A demethylation leads to posttranscriptional repression of LYPD1, and dysregulation of the ALKBH5/LYPD1 axis leads to HCC progression [5]. Wang et al. [6] found that high expression level of circ-KIAA1429 in hepatoma cells and KIAA1429 acts as an oncogene to promote HCC invasion and migration by altering the methylation of m6A in Inhibitor of DNA binding 2(ID2) and GATA-binding protein(GATA3) mRNAs [7]. YTH N6-methyladenosine RNA binding protein F3 (YTHDF3) can increase the stability of zinc finger e-box binding homeobox 1 (Zeb1) mRNA and participate in the occurrence and development of liver cancer [6]. Methyltransferase-like 3 (METTL3) and Methyltransferase-like 14 (METTL14) are the two core molecules. However, METTL3 and METTL14 exert opposing regulatory roles in HCC [8]. Hepatitis B virus X-interacting protein (HBXIP) interference suppressed the malignant behavior of HCC and suppressed the Warburg effect in HCC cells. METTL3 is upregulated in HCC tissues and positively regulated by HBXIP. HBXIP-mediated METTL3 promotes metabolic reprogramming and proliferation, invasion and metastasis of HCC cells [9]. The expression of ubiquitin-specific peptidase 48 (USP48) is significantly reduced in HCC, and methyltransferase-like 14 (Mettl14)-induced m6A modification is involved in the regulation of USP48 in HCC by maintaining USP48 mRNA stability, thereby inhibiting the development of HCC [10]. YTH N6-methyladenosine RNA binding protein 1 (YTHDF1) can promote the occurrence and development of liver cancer through various pathways. Significant overexpression of YTHDF1 in HCC tissues is associated with a poor prognosis, and YTHDF1 deficiency inhibits HCC autophagy, growth and metastasis [11]. YTHDF1 can accelerate translational export of FZD5 mRNA in an m6A-dependent manner and function as an oncogene through the WNT/β-catenin pathway [12]. The m6A demethylase FTO promotes hepatocellular carcinoma tumorigenesis by mediating pyruvate kinase M2 (PKM2) demethylation [13]. NOP2/Sun RNA methyltransferase 2 (NSUN2) is an RNA methyltransferase responsible for m5C modification of multiple RNAs. The H19 lncRNA is a specific target of NSUN2 modifiers. m5C-modified H19 lncRNA may promote tumorigenesis and development by recruiting the G3BP1 protein [14]. NOP2/Sun RNA methyltransferase 4 (NSUN4) is significantly upregulated in tissues and cells of HCC patients and is closely related to the occurrence and development of HCC [15]. Aly/REF export Factor (ALYREF) is significantly upregulated in liver cancer tissues and liver cancer cell lines and is significantly associated with poor prognosis in HCC [16]. TRNA methyltransferase 6 non-catalytic subunit (TRMT6) and TRNA methyltransferase 61A (TRMT61A) form the m1A methyltransferase complex, which is highly expressed in advanced HCC tumors and negatively correlates with HCC survival [17]. Emerging reports confirm that dysregulation of RNA methylation contributes to a variety of human diseases, particularly hepatocellular carcinoma. However, no study has comprehensively analyzed the relationship between these several forms of RNA methylation and HCC.

Therefore, based on the existing research on RNA methylation in HCC, we collected m6A/m5C/m1A-related genes that are closely correlated with the occurrence and development of HCC. We aimed to explore the relationship between the expression levels of m6A/m5C/m1A-related genes and prognosis of HCC patients. Furthermore, we tried to construct and validate a m6A/m5C/m1A-related gene prognostic signature and explore its relationship with the clinicopathological features, immune microenvironment and immunotherapy of HCC, in order to provide individualized strategies for clinical treatment of HCC.

Methods

Data sources

We extracted 10 m6A/m5C/ m1A-related genes that are closely related to the occurrence and development of HCC, and the information on m6A/m5C/m1A-related genes were listed in Table 1. HCC RNA-sequencing (RNA-seq) data and clinical data (Table S1) were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Moreover, we downloaded HCC RNA-sequencing (RNA-seq) data and clinical data (Table S2) from the ICGC database (https://dcc.icgc.org/) as external validation data.

Table 1 Information of m6A/m5C/m1A-related genes in HCC
Full size table

Copy number variation and differential expression analysis of m6A/m5C/m1A-related genes in HCC

The copy number variation frequency information of TCGA HCC samples were downloaded from the UCSC Xena website (https://xena.ucsc.edu/). We used the “limma” package to distinguish differentially expressed genes (DEGs) with p values < 0.05. Then, we constructed a PPI network associated with m6A/m5C/m1A genes using the STRING website (https://string-db.org/) for analysis of interacting genes. Furthermore, the correlation analysis based on each DEG was performed using the R software "reshape2" package and "igraph" packages.

The m6A/m5C/m1A-related molecular subtypes identification

We used the R software "ConsensusClusterPlus" and "limma" packages to classify HCC patients according to the expression of m6A/m5C/m1A-related genes. Using the R software "survival" package, we analyzed the prognosis of patients with different molecular subtypes. Furthermore, the "limma" package was used to analyze the relationship between molecular subtypes and clinicopathological characteristics.

Establishment and validation of m6A/m5C/m1A-related gene prognostic signature

The R software "survival" R package was used to perform univariate Cox analysis on HCC patients in TCGA to study the effect of m6A/m5C/m1A-related genes on prognosis. Then, using the R software package "glmnet", the prognostic genes screened by univariate Cox analysis were used to construct a prognostic signature by LASSO Cox regression analysis. An individualized risk score was obtained based on the expression level of the prognostic genes and the estimated regression coefficients in LASSO Cox regression analysis. The risk score for each HCC patient was calculated by the following formula:

$$\mathrm{Risk}\;\mathrm{Score}=\sum\nolimits_{\mathrm i=1}^{\mathrm n}\left(\mathrm{Expi}\ast\mathrm{Coei}\right)$$

HCC patients were divided into high- and low-risk groups according to the median risk score, and Kaplan‒Meier analysis was performed to compare 1-, 2-, and 3-year overall survival. The R packages "survival", "survminer" and "time ROC" were used to conduct a receiver operating characteristic (ROC) curve analysis across time periods of 1-, 2-, and 3- years. The distribution and survival of patients are displayed by the "bioRiskPlot" function according to the risk score. Principal component analysis (PCA) was used to assess the samples in each risk group based on their similarities and differences. Patients in the ICGC cohort were separated into low- and high-risk groups based on the median risk score from the TCGA cohort. Then, the ICGC cohort was used to validate the prognostic signature.

Prognostic signature genes mRNA expression analysis via RT‒qPCR

Tissue microarray of human liver cancer and paired adjacent normal tissues (MecDNA-HLivH060PG02) was purchased from Shanghai Outdo Biotech Company (Shanghai, China). The study was approved by the Ethics Committee of Shanghai Outdo Biotech Company. The tissue microarray contains 30 paired cancer and non-cancerous liver tissues. To validate gene expression, quantitative real-time PCR (qRT-PCR) analyses were performed. β-Actin served as an internal reference. Primer sequences are listed in Table S3.

Identification of independent prognostic factors

The R software “survival” package was used for univariate and multivariate analyses. In this analysis, variable factors included age, gender, grade, stage and risk score. Furthermore, we used the “limma” and “ggpubr” packages to analyze the relationship between risk scores and clinicopathological features.

Functional enrichment and immune microenvironment analysis

HCC patients were divided into two groups according to the median risk score. DEGs were screened between the high- and low-risk groups (FDR < 0.05, |log2FC|≥ 1). On this basis, the "clusterProfiler" software package was used for GO and KEGG analysis [18,19,20]. The "GSVA" package was used for single-sample gene set enrichment analysis (ssGSEA) to calculate the scores of infiltrating immune cells and the activity of immune-related pathways in high- and low-risk groups of HCC patients. We further used the R software "limma" package to analyze the differential expression of immune checkpoints between high and low risk groups.

Gene mutation analysis

We further downloaded the HCC simple nucleotide variation data from the TCGA GDC database (https://portal.gdc.cancer.gov/) and then used the R software "maftools" to analyze the gene mutations of the high- and low-risk groups and draw a waterfall diagram. Furthermore, we analyzed the relationship between risk group and tumor mutation burden.

Development and validation of a nomogram

We used the "rms" package to construct a prediction nomogram based on the clinical features and prognostic signature risk score. In the nomogram scoring approach, each variable was given a score, and the total score was computed by adding the scores for all variables in each sample. The calibration plots of the nomogram were used to highlight the predictive value by comparing the predictions vs. actual outcomes in terms of 1-, 2-, and 3-year survival. Area under the curve of ROC was used to evaluate the ability of different prognostic factors to predict the 1-, 2-, and 3-year survival of HCC patients.

Statistical analyses

R version 4.1.0 was used for all statistical analyses. The level of statistical significance was set at p < 0.05.

Results

The landscape of m6A/m5C/m1A-related genes in HCC

Figure 1 illustrates our study process. We investigated copy number variation (CNV) in m6A/m5C/m1A-related genes and discovered that CNV was prevalent in all genes. ALKBH5, METTL3, YTHDF3, YTHDF1, FTO, NSUN2 and ALYREF exhibited noticeable acquired CNV, but FTO, NSUN4, and METTL14 had reduced CNV (Fig. 2a). Figure 2b depicts the position of the m6A/m5C/m1A-related genes on their respective chromosomes. Based on TCGA data including 376 HCC and 50 normal liver tissues, we compared the expression levels of 10 m6A/m5C/m1A-related genes and finally found 9 DEGs (p < 0.05). ALKBH5, METTL3, YTHDF3, YTHDF1, FTO, NSUN2, NSUN4, ALYREF and TRMT6 were significantly upregulated in HCC tissues vs. normal tissues (Fig. 2c). The PPI network revealed the interaction of m6A/m5C/m1A-related genes (Fig. 2d). Figure 2e shown the correlation network of m6A/m5C/m1A-related genes.

Fig. 1
figure 1

Flow chart of this study

Full size image
Fig. 2
figure 2

Expression and interaction of m6A/m5C/m1A-related genes. a Frequencies of CNV gain, loss, and no CNV among m6A/m5C/m1A-related genes. b Locations of CNV alterations in m6A/m5C/m1A-related genes on 23 chromosomes. c Heatmap of the genes concerning m6A/m5C/m1A between HCC tissues and normal tissues. d The PPI network indicated the interaction of genes concerning m6A/m5C/m1A (interaction score = 0.4). e Relevant network of the genes concerning m6A/m5C/m1A. CNV, copy number variation; HCC, hepatocellular carcinoma. PPI, protein–protein interaction. p values: *p < 0.05, ** p < 0.01, and *** P < 0.001

Full size image

Cluster analysis based on the expression of m6A/m5C/m1A-related genes in HCC

To better understand the expression characteristics of m6A/m5C/m1A-related genes in HCC, we used a consensus clustering method to classify HCC patients based on the expression levels of 10 m6A/m5C/m1A-related genes. Our results indicated that k = 2 appeared to be the best choice to classify all HCC patients into the cluster C1 (n = 157) and cluster C2 (n = 213) (Fig. 3a-c). The Kaplan‒Meier curve revealed that the overall survival/prognosis of the cluster C2 was significantly better than that of the cluster C1 (p < 0.001) (Fig. 3d). DEGs of the two molecular subtypes were identified (Fig. 3e). According to the DEGs expression of subtypes, we analyzed the correlation between molecular subtypes and clinicopathological features. Highly expressed genes in cluster C1 are significantly correlated with G3-4, T3-4, stage III-IV (p < 0.05) (Fig. 3f).

Fig. 3
figure 3

Identification of m6A/m5C/m1A-related subtypes. a A total of 370 HCC patients were classified into two subtypes based on consensus cluster analysis. b CDF for k = 2 to 9. c Relative changes in the areas under the CDF curve for k = 2 to 9. d OS of different m6A/m5C/m1A-related subtypes. e Identification of DEGs between the two subtypes. f The relationship between m6A/m5C/m1A-related subtypes and clinicopathological features in HCC. CDF, cumulative distribution function; OS, overall survival; DEGs, differentially expressed genes; HCC, hepatocellular carcinoma

Full size image

Construction of a prognostic signature based on the TCGA cohort

We used univariate Cox analysis to screen out prognostic m6A/m5C/m1A-related genes for subsequent analysis. A total of 6 genes were identified (Fig. 4a). Furthermore, we used the TCGA dataset as the training set (n = 370). LASSO regression analysis was performed on 6 prognostic m6A/m5C/m1A-related genes to further select the best prognostic indicators (Fig. 4b-c). Finally, 4 prognostic genes (METTL3, YTHDF1, NSUN4, and TRMT6) were retained with the least partial likelihood deviation. We developed a risk score model based on the following formula: risk score = (0.02* expression of METTL3) + (0.43* expression of YTHDF1) + (0.46* expression of NSUN4) + (0.40* expression of TRMT6). A total of 370 HCC patients in the training set were divided into high- and low-risk groups based on the median risk score, and the difference in OS between the high- and low-risk groups was statistically significant (p = 0.002, Fig. 4d). The areas under the ROC curve of the risk signature for the 1-, 2-, and 3-year periods were 0.739, 0.649, and 0.664, respectively (Fig. 4e). Compared with patients in the low-risk group, patients in the high-risk group showed shorter survival time and higher risk of death (Fig. 4f-g). PCA and t-SNE (Fig. 4h-i) showed that the patients were well divided into two risk groups.

Fig. 4
figure 4

Establishing a prognostic signature based on m6A/m5C/m1A-related genes. a Univariate Cox regression analysis of m6A/m5C/m1A-related genes. b-c LASSO regression analysis and partial likelihood deviance of the prognostic genes. d OS analysis of different risk groups. e ROC curves about the prognostic signature predictive power. f Distribution of survival status in high- and low- risk groups. g Risk score curve about the high- and low-risk groups. h PCA and (i) tSNE results of the prognostic signature. LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; PCA, principal component analysis; tSNE, t-distributed stochastic neighbor embedding

Full size image

Validation of the prognostic signature in an external cohort

A total of 232 HCC patients in the ICGC database were used as the external validation cohort. According to the median risk score in the training set, 192 patients in the ICGC cohort were assigned to the high-risk group, and 40 patients were assigned to the low-risk group. There was a significant difference in OS between the high- and low-risk groups (p = 0.008, Fig. 5a). Time-dependent ROC analysis was used to assess the sensitivity and specificity of the prognostic model, and areas under curves for 1-, 2- and 3-year survival were 0.630, 0.671, and 0.689, respectively (Fig. 5b). Compared with patients in the low-risk group, patients in the high-risk group had a higher risk of death and shorter survival time (Fig. 5c-d). PCA and t-SNE (Fig. 5e-f) showed that the patients were well divided into two risk groups.

Fig. 5
figure 5

Validation of the prognostic signature based on the ICGC cohort. a OS analysis of risk groups. b ROC curves about the signature predictive power. c Risk score curve about the high- and low-risk groups. d Distribution of patients’ survival status. e PCA and (f) tSNE results of risk groups. ICGC, international cancer genome consortium; OS, overall survival; ROC, receiver operating characteristic; PCA, principal component analysis; tSNE, t-distributed stochastic neighbor embedding

Full size image

Independent prognostic value of the prognostic signature

We used univariate and multivariate Cox regression analyses to assess whether the risk score of the m6A/m5A/m1A-related gene prognostic signature could serve as an independent prognostic factor. Univariate Cox regression analysis showed the risk score was an independent factor for predicting poorer survival (HR = 2.995, 95% CI: 1.903–4.713, Fig. 6a). Multivariate analysis showed that the risk score was a prognostic factor after adjustment for other confounding factors (HR = 2.577, 95% CI: 1.613–4.117, Fig. 6b). Additionally, we drew a heatmap to analyze the relationship between m6A/m5C/m1A-related genes and clinical characteristics in the TCGA cohort (Fig. 6c) and found that the risk score was correlated with grade and stage in HCC patients (p < 0.05).

Fig. 6
figure 6

Correlation analysis between risk score and clinicopathological features. a Univariate and (b) multivariate Cox regression analyses. c Relationship between risk score and clinicopathological features

Full size image

Validation of the expression of the 4 prognostic genes

The expression levels of 4 prognostic genes were determined using RT‒qPCR in 30 HCC and 30 neighboring normal tissues. As shown in Fig. 7a-d, TRMT6, NSUN4, METTL3 and YTHDF1 were significantly overexpressed in HCC tissues (p < 0.05).

Fig. 7
figure 7

Validation of the expression of the prognostic signature’s genes. The expression of (a) TRMT6, (b) NSUN4, (c) METTL3 and (d) YTHDF in liver cancer and normal tissues. *p < 0.05, **p < 0.01, and ***p < 0.001

Full size image

Functional enrichment analysis based on the different risk groups

To further explore differences in gene function and pathways between risk groups, we employed the "limma" R package to extract DEGs with FDR < 0.05 and |log2FC|≥ 1. In the TCGA cohort, a total of 187 DEGs were identified between the low- and high-risk groups. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed based on these DEGs. The results showed that DEGs were associated with small molecule catabolic process, steroid metabolic process, carboxylic acid catabolic process, organic acid catabolic process, olefinic compound metabolic process, retinol metabolism (Fig. 8a-b).

Fig. 8
figure 8

Functional enrichment analysis. a GO and (b) KEGG analysis about DEGs in the low- and high-risk groups. DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genomes

Full size image

Comparison of the immune microenvironment between high- and low-risk groups

We used the ssGSEA method to compare the enrichment scores of 8 immune cell activities and 10 immune-related pathways between the high- and low-risk groups. The levels of B cells, neutrophils, DCs, Th2 cells and TIL expression levels were significantly increased in the low-risk group (Fig. 9a). The high-risk subgroup generally had lower activity for all immune pathways, except the MHC I pathway (Fig. 9b). High-risk group patients showed significant overexpression of HAVCR2, PDCD1, CTLA4, CD274, and TIGIT (p < 0.05) (Fig. 9c).

Fig. 9
figure 9

Correlation analysis between risk score and immune microenvironment. Boxplots show the scores of immune cells (a) and immune-related functions (b) in different risk groups. c Expression of immune checkpoints in differential risk groups. *, P < 0.05; **, P < 0.01; ***, P < 0.001

Full size image

Mutation analysis

Gene mutation analysis showed that the frequency of gene mutations was higher in the high-risk group than in the low-risk group (Fig. 10a-b). Patients in the high-risk group had increased TMB compared with those in the low-risk group (p < 0.05) (Fig. 10c). Patients in the high-TMB group had worse OS than those in the low-TMB group (p < 0.05) (Fig. 10d). The survival/prognosis of patients with high TMB in the high-risk group was significantly worse than that of patients with low TMB (Fig. 10e).

Fig. 10
figure 10

Gene mutation analysis based on the prognostic signature. Analysis of gene mutation frequencies in the low-risk (a) and high-risk (b) groups. c Analysis of TMB differences between the high-risk and low-risk groups. d Survival analysis based on tumor mutation burden. e Survival analysis based on TMB combined with the risk score. TMB, tumor mutation burden

Full size image

Development of a nomogram to predict survival

We established a nomogram containing the risk score and clinicopathological characteristics to predict overall survival (Fig. 11a). Calibration curves of the nomogram for predicting 1-, 2- and 3-year OS suggested that the performance of the proposed nomogram was ideal (Fig. 11b). In predicting the survival/prognosis of patients, the 1-, 2- and 3-year ROC AUCs of the nomogram were better than those of clinicopathological features (Fig. 11c-e).

Fig. 11
figure 11

Establishment and validation of the nomogram based on clinicopathological characteristics and risk score. a Nomogram for risk score and clinical characteristics in patients with HCC. b Calibration analysis of the risk score containing nomogram for 1-, 2-, 3- year OS. c 1-year, (d) 2-year and (e) 3-year ROC analysis for the nomogram. *p < 0.05, **p < 0.01, and ***p < 0.001. HCC, hepatocellular carcinoma; OS, overall survival; ROC, receiver operating characteristic curve

Full size image

Discussion

Nonmutational epigenetic reprogramming is considered to be an important cancer hallmark, and epigenetic modification is closely related to the occurrence and development of tumors [21]. Recent studies have demonstrated that RNA modifications may affect RNA metabolism, splicing, stability and translation, thereby affecting gene expression [22, 23]. RNA methylation is an important epigenetic modification that does not alter gene sequence but may play an essential role in multiple biological processes, such as gene expression, genome editing, and cell differentiation. With advances in RNA detection, various forms of RNA methylation, including m6A, m5C, m1A,, can be assessed.

A large number of studies have shown that RNA methylation modification was associated with the development of various malignant tumors. Cui et al. [24] found that m6A is involved in the regulation of glioma stem cell (GSC) self-renewal and differentiation. Knocking down the expression of the methyltransferase METTL3 and METTL14 proteins in cells can induce the mRNA expression of the proto-oncogenes ADAM19, EPHA3 and KLF4 in GSCs and promote the growth, self-renewal and differentiation of GSCs. Liu et al. [25] found that m6A methylation of EphA2 and VEGFA can promote angiogenesis in colorectal cancer. Wu et al. [26] found that protein arginine methyltransferase 5 (PRMT5) contributes to doxorubicin resistance in breast cancer by enhancing nuclear translocation of the RNA demethylase ALKBH5. In addition, a number of studies have confirmed that m6A modification is closely related to tumor immunotherapy [27,28,29], chemo-radiotherapy [30], and targeted therapy [31, 32] resistance. Hu et al. [33] found that m5C RNA methyltransferase NSUN2 involved in the proliferation, invasion and migration of gastric cancer cells. Wang et al. [34] found that overexpression of ALYREF can promote bladder cancer cell proliferation through PKM2-mediated glycolysis. The known demethylases for m1A modification include ALKBH3 and ALKBH1; methyltransferase enzymes include the TRMT family; and RNA-binding proteins include YTH family proteins [35, 36]. A large number of studies have shown that m1A also closely related to tumor proliferation, invasion, metabolism, and immune microenvironment. In most studies, m1A writers TRMT6 and eraser ALKBH3 are associated with poorer prognosis in multiple cancers [37,38,39,40] In this study, we comprehensively analyzed the relationship between RNA methylation types and HCC. In this study, two different molecular subtypes were identified based on the expression of m6A/m5C/m1A-related genes. The overall survival and prognosis of HCC patients in the two subtypes were significantly different. Also, the clinicopathologic features of different molecular subtypes were significantly different. The risk prognosis model based on m6A/m5C/m1A related genes showed that patients in the high-risk group had a significantly worse prognosis. Our results further confirmed that RNA modification is widely involved in the occurrence and development of liver cancer and affects the prognosis of patients.

Previous studies have shown that RNA methylation closely related to the tumor microenvironment. Researchers found that ALKBH5 promotes HCC growth, metastasis and macrophage recruitment through the ALKBH5/MAP3K8 axis [41]. Liu et al. [42] found that METTL3 mediated m6A methylation modification of circIGF2BP3 to suppress CD8 + T cell responses and promote immune escape in non-small cell lung cancer. Gao et al. [43]found that m1A modification related to the immune microenvironment of colon cancer through bioinformatics research. An important aspect of our study is that we explored the correlation between m6A/m5C/m1A-related genes and tumor immune microenvironment in HCC patients. Interestingly, we found that tumor-infiltrating lymphocyte (TIL) and B-cell infiltration levels were significantly increased in the low-risk group. Previous studies have confirmed that the presence of TILs and B cells is associated with improved survival in HCC patients [44, 45]. Furthermore, except for the MHC response pathway, the activation level of immune pathways in the high-risk group was lower than that in the lower-risk group. Based on these findings, the poor survival outcomes of HCC patients in the high-risk group may be due to reduced levels of antitumor immunity.

With the further study of tumor immunology and molecular biology, immunotherapy provides a new direction for the treatment of tumors. This immunotherapy includes immune checkpoint inhibitors (ICIs), therapeutic antibodies, and cell therapy. Research on ICIs for CTLA-4, PD-1, and PD-L1 is booming, and clinical studies have demonstrated safety and efficacy [46, 47]. We evaluated immune checkpoint-associated gene expression in the low- and high-risk groups. The expression of immune checkpoints was significantly increased in the high-risk group. These findings indicated that the high-risk group of HCC patients with suppressed immune function may benefit from immune checkpoint inhibitor therapy. Accumulating evidence suggests that tumor mutational burden (TMB) can serve as an important prognostic marker in relation to immunotherapy [48, 49]. We also compared differences in TMB and observed that subgroups with low risk scores tended to have lower TMB, with a significant survival benefit in patients with low TMB compared with those with high TMB. These results suggest that the m6A/m5C/m1A-related prognostic signature can be used for individualized treatment of HCC patients.

However, our study still has some limitations. First, prognostic markers were established and validated using retrospective data, and their clinical applicability needs to be validated with prospective data. Second, the potential biological functions and specific molecular mechanisms of these three m6A/m5C/m1A-related genes need to be further investigated. Third, the correlation between the risk score and tumor immunity was not experimentally demonstrated.

Conclusions

In conclusion, our study further confirmed that m6A, m5C and m1A are closely related to HCC. Furthermore, in the TCGA and ICGC cohorts, the risk score generated by the risk model based on the 4 m6A/m5C/m1A-related genes was an independent risk factor for predicting the overall survival/prognosis of HCC patients. DEGs between the low-risk and high-risk groups were associated with tumor immunity. The combined analysis of RNA methylation-related genes in this study provides a new prognostic signature for HCC patients and provides an important basis for further research on the relationship between RNA methylation and immune function in HCC patients.

Availability of data and materials

The datasets analyzed during the current study are available in The Cancer Genome Atlas (TCGA) repository (https://portal.gdc.cancer.gov/), the International Cancer Genome Consortium (ICGC) repository (https://dcc.icgc.org/) and the STRING website (https://string-db.org/).

Abbreviations

m6A:

N6-methyladenosine

m5C:

5-Methylcytosine

m1A:

N1-methyladenosine

HCC:

Hepatocellular cancer

TCGA:

The Cancer Genome Atlas database

ICGC:

International Cancer Genome Consortium

CNV:

Copy number variation

PPI:

Protein–protein interaction

CDF:

Cumulative distribution function

OS:

Overall survival

LASSO:

Least absolute shrinkage and selection operator

ROC:

Receiver operating characteristic

DEGs:

Differentially expressed genes

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genomes

PCA:

Principal component analysis

tSNE:

T-distributed stochastic neighbor embedding

TMB:

Tumor mutation burden

References

  1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021;71:7–33.

    Article  PubMed  Google Scholar 

  2. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49.

    Article  PubMed  Google Scholar 

  3. Yang JD, Hainaut P, Gores GJ, Amadou A, Plymoth A, Roberts LR. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol. 2019;16:589–604.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Henley SJ, Ward EM, Scott S, Ma J, Anderson RN, Firth AU, et al. Annual report to the nation on the status of cancer, part I: National cancer statistics. Cancer. 2020;126:2225–49.

    Article  PubMed  Google Scholar 

  5. Chen Y, Zhao Y, Chen J, Peng C, Zhang Y, Tong R, et al. ALKBH5 suppresses malignancy of hepatocellular carcinoma via m6A-guided epigenetic inhibition of LYPD1. Mol Cancer. 2020;19:123.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Wang M, Yang Y, Yang J, Yang J, Han S. circ_KIAA1429 accelerates hepatocellular carcinoma advancement through the mechanism of m6A-YTHDF3-Zeb1. Life Sci. 2020;257:118082.

    Article  CAS  PubMed  Google Scholar 

  7. Cheng X, Li M, Rao X, Zhang W, Li X, Wang L, et al. KIAA1429 regulates the migration and invasion of hepatocellular carcinoma by altering m6A modification of ID2 mRNA. Onco Targets Ther. 2019;12:3421–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Liu X, Qin J, Gao T, Li C, Chen X, Zeng K, et al. Analysis of METTL3 and METTL14 in hepatocellular carcinoma. Aging (Albany NY). 2020;12:21638–59.

    Article  CAS  PubMed  Google Scholar 

  9. Yang N, Wang T, Li Q, Han F, Wang Z, Zhu R, et al. HBXIP drives metabolic reprogramming in hepatocellular carcinoma cells via METTL3-mediated m6A modification of HIF-1α. J Cell Physiol. 2021;236:3863–80.

    Article  CAS  PubMed  Google Scholar 

  10. Du L, Li Y, Kang M, Feng M, Ren Y, Dai H, et al. USP48 is upregulated by Mettl14 to attenuate hepatocellular carcinoma via regulating SIRT6 stabilization. Cancer Res. 2021;81:3822–34.

    Article  CAS  PubMed  Google Scholar 

  11. Li Q, Ni Y, Zhang L, Jiang R, Xu J, Yang H, et al. HIF-1α-induced expression of m6A reader YTHDF1 drives hypoxia-induced autophagy and malignancy of hepatocellular carcinoma by promoting ATG2A and ATG14 translation. Signal Transduct Target Ther. 2021;6:76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Liu X, Qin J, Gao T, Li C, He B, Pan B, et al. YTHDF1 facilitates the progression of hepatocellular carcinoma by promoting FZD5 mRNA translation in an m6A-dependent manner. Mol Ther Nucleic Acids. 2020;22:750–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Li J, Zhu L, Shi Y, Liu J, Lin L, Chen X. m6A demethylase FTO promotes hepatocellular carcinoma tumorigenesis via mediating PKM2 demethylation. Am J Transl Res. 2019;11:6084–92.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Sun Z, Xue S, Zhang M, Xu H, Hu X, Chen S, et al. Aberrant NSUN2-mediated m5C modification of H19 lncRNA is associated with poor differentiation of hepatocellular carcinoma. Oncogene. 2020;39:6906–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Cui M, Qu F, Wang L, Liu X, Yu J, Tang Z, et al. m5C RNA methyltransferase-related gene NSUN4 stimulates malignant progression of hepatocellular carcinoma and can be a prognostic marker. Cancer Biomark. 2022;33:389–400.

    Article  CAS  PubMed  Google Scholar 

  16. Xue C, Zhao Y, Li G, Li L. Multi-omic analyses of the m5C regulator ALYREF reveal its essential roles in hepatocellular carcinoma. Front Oncol. 2021;11: 633415.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51:D587–92.

    Article  CAS  PubMed  Google Scholar 

  20. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Hanahan D. Hallmarks of cancer: new dimensions. Cancer Discov. 2022;12:31–46.

    Article  CAS  PubMed  Google Scholar 

  22. Xue C, Zhao Y, Li L. Advances in RNA cytosine-5 methylation: detection, regulatory mechanisms, biological functions and links to cancer. Biomark Res. 2020;8:43.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Nombela P, Miguel-López B, Blanco S. The role of m6A, m5C and Ψ RNA modifications in cancer: novel therapeutic opportunities. Mol Cancer. 2021;20:18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, et al. m6A RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells. Cell Rep. 2017;18:2622–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Liu X, He H, Zhang F, Hu X, Bi F, Li K, et al. m6A methylated EphA2 and VEGFA through IGF2BP2/3 regulation promotes vasculogenic mimicry in colorectal cancer via PI3K/AKT and ERK1/2 signaling. Cell Death Dis. 2022;13:483.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Wu Y, Wang Z, Han L, Guo Z, Yan B, Guo L, et al. PRMT5 regulates RNA m6A demethylation for doxorubicin sensitivity in breast cancer. Mol Ther. 2022;30:2603–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yin H, Zhang X, Yang P, Zhang X, Peng Y, Li D, et al. RNA m6A methylation orchestrates cancer growth and metastasis via macrophage reprogramming. Nat Commun. 2021;12:1394.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Li N, Kang Y, Wang L, Huff S, Tang R, Hui H, et al. ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment. Proc Natl Acad Sci U S A. 2020;117:20159–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Qiu X, Yang S, Wang S, Wu J, Zheng B, Wang K, et al. M(6)A demethylase ALKBH5 regulates PD-L1 expression and tumor immunoenvironment in intrahepatic cholangiocarcinoma. Cancer Res. 2021;81:4778–93.

    Article  CAS  PubMed  Google Scholar 

  30. Zhou S, Bai ZL, Xia D, Zhao ZJ, Zhao R, Wang YY, et al. FTO regulates the chemo-radiotherapy resistance of cervical squamous cell carcinoma (CSCC) by targeting β-catenin through mRNA demethylation. Mol Carcinog. 2018;57:590–7.

    Article  CAS  PubMed  Google Scholar 

  31. Xiao P, Liu YK, Han W, Hu Y, Zhang BY, Liu WL. Exosomal delivery of FTO confers Gefitinib resistance to recipient cells through ABCC10 regulation in an m6A-dependent manner. Mol Cancer Res. 2021;19:726–38.

    Article  CAS  PubMed  Google Scholar 

  32. Yan F, Al-Kali A, Zhang Z, Liu J, Pang J, Zhao N, et al. A dynamic N(6)-methyladenosine methylome regulates intrinsic and acquired resistance to tyrosine kinase inhibitors. Cell Res. 2018;28:1062–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Hu Y, Chen C, Tong X, Chen S, Hu X, Pan B, et al. NSUN2 modified by SUMO-2/3 promotes gastric cancer progression and regulates mRNA m5C methylation. Cell Death Dis. 2021;12:842.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Wang JZ, Zhu W, Han J, Yang X, Zhou R, Lu HC, et al. The role of the HIF-1α/ALYREF/PKM2 axis in glycolysis and tumorigenesis of bladder cancer. Cancer Commun (Lond). 2021;41:560–75.

    Article  PubMed  Google Scholar 

  35. Li X, Xiong X, Zhang M, Wang K, Chen Y, Zhou J, et al. Base-resolution mapping reveals distinct m1A methylome in nuclear- and mitochondrial-encoded transcripts. Mol Cell. 2017;68:993-1005.e9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Dai X, Wang T, Gonzalez G, Wang Y. Identification of YTH domain-containing proteins as the readers for N1-methyladenosine in RNA. Anal Chem. 2018;90:6380–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Woo HH, Chambers SK. Human ALKBH3-induced m(1)A demethylation increases the CSF-1 mRNA stability in breast and ovarian cancer cells. Biochim Biophys Acta Gene Regul Mech. 2019;1862:35–46.

    Article  CAS  PubMed  Google Scholar 

  38. Chen K, Shen D, Tan L, Lai D, Han Y, Gu Y, et al. A pan-cancer analysis reveals the prognostic and immunotherapeutic value of ALKBH7. Front Genet. 2022;13:822261.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Zhao Y, Zhao Q, Kaboli PJ, Shen J, Li M, Wu X, et al. m1A regulated genes modulate PI3K/AKT/mTOR and ErbB pathways in gastrointestinal cancer. Transl Oncol. 2019;12:1323–33.

    Article  PubMed  PubMed Central  Google Scholar 

  40. 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:6314.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. You Y, Wen D, Zeng L, Lu J, Xiao X, Chen Y, et al. ALKBH5/MAP3K8 axis regulates PD-L1+ macrophage infiltration and promotes hepatocellular carcinoma progression. Int J Biol Sci. 2022;18:5001–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Liu Z, Wang T, She Y, Wu K, Gu S, Li L, et al. N(6)-methyladenosine-modified circIGF2BP3 inhibits CD8(+) T-cell responses to facilitate tumor immune evasion by promoting the deubiquitination of PD-L1 in non-small cell lung cancer. Mol Cancer. 2021;20:105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Gao Y, Wang H, Li H, Ye X, Xia Y, Yuan S, et al. Integrated analyses of m(1)A regulator-mediated modification patterns in tumor microenvironment-infiltrating immune cells in colon cancer. Oncoimmunology. 2021;10:1936758.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Lu C, Rong D, Zhang B, Zheng W, Wang X, Chen Z, et al. Current perspectives on the immunosuppressive tumor microenvironment in hepatocellular carcinoma: challenges and opportunities. Mol Cancer. 2019;18:130.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Garnelo M, Tan A, Her Z, Yeong J, Lim CJ, Chen J, et al. Interaction between tumour-infiltrating B cells and T cells controls the progression of hepatocellular carcinoma. Gut. 2017;66:342–51.

    Article  CAS  PubMed  Google Scholar 

  46. Duffy AG, Ulahannan SV, Makorova-Rusher O, Rahma O, Wedemeyer H, Pratt D, et al. Tremelimumab in combination with ablation in patients with advanced hepatocellular carcinoma. J Hepatol. 2017;66:545–51.

    Article  CAS  PubMed  Google Scholar 

  47. Yau T, Hsu C, Kim TY, Choo SP, Kang YK, Hou MM, et al. Nivolumab in advanced hepatocellular carcinoma: Sorafenib-experienced Asian cohort analysis. J Hepatol. 2019;71:543–52.

    Article  CAS  PubMed  Google Scholar 

  48. Bonneville R, Krook MA, Kautto EA, Miya J, Wing MR, Chen HZ, et al. Landscape of microsatellite instability across 39 cancer types. JCO Precis Oncol. 2017;2017:PO.17.00073.

    PubMed  Google Scholar 

  49. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51:202–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We acknowledge the TCGA and ICGC databases for providing their platforms and contributors for uploading their meaningful datasets.

Funding

This research was funded by the National Clinical Key Specialty Construction Program, 2021; Fujian Research and Training Grants for Young and Middle-aged Leaders in Healthcare, 2022; Fujian Provincial Clinical Research Center for Cancer Radiotherapy and Immunotherapy, 2020Y2012; and Fujian Clinical Research Center for Radiation and Therapy of Digestive, Respiratory and Genitourinary Malignancies, 2021Y2014.

Author information

Authors and Affiliations

  1. Department of Radiation Oncology, Clinical Oncology School of Fujian Medical University, Fujian Cancer Hospital, Fuzhou, 350014, China

    Yu Xiao, Jinluan Li & Junxin Wu

Authors
  1. Yu Xiao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jinluan Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Junxin Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YX designed the study, YX and JL collected and analyzed the data, JL conducted the experiments, YX and JL prepared the figures and tables, YX wrote the manuscript, JL and JW revised the article critically.

Corresponding authors

Correspondence to Jinluan Li or Junxin Wu.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Ethics Committee of Shanghai Outdo Biotech Company. Ethics approval number: SHYJS-CP-1804023.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

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: Table S1.

Clinical data of TCGA hepatocellular cancer dataset. Table S2. Clinical data of ICGC hepatocellular cancer dataset. Table S3. The primers’ sequences for qRT PCR analysis.

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 http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) 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

Xiao, Y., Li, J. & Wu, J. Development and validation of a novel prognostic signature based on m6A/m5C/m1A-related genes in hepatocellular carcinoma. BMC Med Genomics 16, 177 (2023). https://doi.org/10.1186/s12920-023-01611-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-023-01611-x

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us