N6-methyladenosine regulators in hepatocellular carcinoma: investigating the precise definition and clinical applications of biomarkers

Abstract

Background

Accurately identifying effective biomarkers and translating them into clinical practice have significant implications for improving clinical outcomes in hepatocellular carcinoma (HCC). In this study, our objective is to explore appropriate methods to improve the accuracy of biomarker identification and investigate their clinical value.

Methods

Concentrating on the N6-methyladenosine (m6A) modification regulators, we utilized dozens of multi-omics HCC datasets to analyze the expression patterns and genetic features of m6A regulators. Through the integration of big data analysis with function experiments, we have redefined the biological roles of m6A regulators in HCC. Based on the key regulators, we constructed m6A risk models and explored their clinical value in estimating prognosis and guiding personalized therapy for HCC.

Results

Most m6A regulators exhibit abnormal expression in HCC, and their expression is influenced by copy number variations (CNV) and DNA methylation. Large-scale data analysis has revealed the biological roles of many key m6A regulators, and these findings are well consistent with experimental results. The m6A risk models offer significant prognostic value. Moreover, they assist in reassessing the therapeutic potential of drugs such as sorafenib, gemcitabine, CTLA4 and PD1 blockers in HCC.

Conclusions

Our findings suggest that the mutual validation of big data analysis and functional experiments may facilitate the precise identification and definition of biomarkers, and our m6A risk models may have the potential to guide personalized chemotherapy, targeted treatment, and immunotherapy decisions in HCC.

Introduction

Liver cancer is a prevalent malignant tumor that significantly impacts public health, ranking sixth in incidence and third in mortality [1]. The primary pathological type is hepatocellular carcinoma (HCC), which accounts for 75–85% of cases. The widespread adoption of comprehensive, precise, and standardized treatments, along with the continuous development of various techniques such as surgery, interventional therapy, radio-chemotherapy, targeted therapy, and immunotherapy, has improved patient prognosis. However, the overall prognosis for HCC remains poor, particularly for patients with advanced HCC. The majority of these patients face challenges in effectively controlling disease progression, resulting in an approximate five-year overall survival rate of only 10% [2]. Therefore, the development of targeted drugs and establishment of accurate and reliable molecular subtypes for HCC are urgently needed.

Accurately determining the role of biomarkers in tumors is crucial for both targeted drug development and constructing molecular subtypes to guide individualized treatment decisions. Recent advances in high-throughput technologies and the accumulation of extensive omics data have significantly improved our understanding of cancer biology. An increasing number of studies recommend combining data from multiple centers and large datasets to achieve more reliable results [3]. Through mutual verification of data analysis and fundamental experiments, the biological functions of molecules can be more precisely delineated.

N6-methyladenosine (m6A) modification is a prevalent epigenetic phenomenon in eukaryotic cell biology, influencing both mRNA and non-coding RNA [4]. Current research highlights the critical roles of m6A regulators in processes such as RNA stability, splicing, maturation, subcellular localization, translation, protein interactions, and microRNA processing. These regulators govern various biological behaviors, encompassing stem cell self-renewal, differentiation, and tumor initiation and metastasis [5, 6]. For instance, the m6A demethylase FTO has been found to be upregulated in HCC and shown to promote HCC cell proliferation, migration, and invasion [7]. By targeting the m6A reader protein YTHDF2 with small interfering RNAs, it is possible to inhibit immune evasion and angiogenesis in HCC [8]. The m6A binding protein YTHDF1 can enhance stemness and therapeutic resistance in HCC through upregulating NOTCH1 expression [9]. Furthermore, the m6A methyltransferase METTL16 contributes to HCC pathogenesis and bolsters cancer stem cell (CSC) self-renewal by increasing the efficiency of mRNA translation [10]. Collectively, these studies underscore the significance of m6A regulators in HCC. Accordingly, in this study, we focus on m6A modifications, employing big data analysis and mutual verification with basic research to evaluate the characteristics and biological functions of m6A regulators in HCC, and aim to establish a foundation for clinical translation.

Materials and methods

Single-cell RNA-seq data analysis

We downloaded the single-cell RNA-seq dataset GSE125449 from the GEO database, wherein the set1 dataset, predominantly comprising HCC samples, was employed to examine gene expression across various cell types. Cell types were annotated according to the information furnished by the original authors. Data analysis was conducted using Seurat v4.3.0, and visualization was achieved through a UMAP plot generated by the DimPlot function at a resolution of 0.2.

Assessing correlations in multiple datasets

Investigating relationships between various factors across dozens of datasets is intricate and difficult to depict clearly. To tackle this, we evaluate the strength and direction of variable relationships by calculating the sum of correlation coefficients (Cumulative correlation). For example, when examining relationships among genes A, B, and C, we compute Pearson correlation coefficients between them across 31 transcriptome datasets and sum the results. Assuming the analysis shows a cumulative correlation of − 5 between genes A and C and 4 between genes B and C, it implies a negative correlation for A and C, a positive correlation for B and C, with a stronger association between A and C.

Construction and validation of the m6A models

To identify key m6A regulators in HCC and construct molecular models, we calculated the Pearson correlation between 29 m6A regulators and 12 HCC-related signatures (reflecting characteristics such as proliferation, metastasis, vascular invasion, tumor stemness, and prognosis) in each transcriptome dataset. By summing the Pearson correlations of each m6A regulator with the 5 beneficial HCC-related signatures or the 7 harmful signatures across 31 HCC datasets, we obtained cumulative correlation coefficients. Based on the absolute values of these cumulative correlation coefficients, we identified m6A regulators most closely associated with beneficial and harmful signatures. Subsequently, by analyzing the impact of these m6A regulators on HCC prognosis and their biological functions in HCC from previous studies, we compared and selected harmful m6A regulators and beneficial regulators. Based on these key factors, we constructed m6A.ES.benefit and m6A.ES.harm models using the GSVA method. To validate the clinical value of these two models in HCC, we used correlation analysis, t-tests, and non-parametric Mann-Whitney U-tests to analyze the associations between these two models and various HCC clinical and molecular features, such as HCC-related signatures, TNM staging, histological differentiation, tumor doubling time, proliferation, metastasis, and tumor stemness. In addition, Receiver operating characteristic (ROC) analysis were performed to validate the predictive ability of the m6A models. P-values less than 0.05 were considered statistically significant.

Drug sensitivities analysis

To explore the potential of m6A models in predicting chemotherapy drug response, we used the pRRophetic algorithm, suitable for untransformed log2 gene expression microarray data, to assess the sensitivity of over 200 drugs in 24 HCC microarray datasets. We applied the pRRopheticPredict function with tissueType and dataset parameters set to “all” and “cgp2016”, respectively. Next, we employed the calcPhenotype function of the oncoPredict method to examine drug sensitivity in 7 HCC sequencing datasets, using “GDSC2” as the drug screening dataset. Based on the results of pRRophetic, we calculated the cumulative correlation between drug sensitivity, m6A.ES.benefit, and m6A.ES.harm. By cross-validating these findings with the results of oncoPredict, we identified the drugs most closely related to m6A.ES.benefit and m6A.ES.harm.

Tumor immune microenvironment analysis

The CIBERSORT algorithm and GSVA’s ssgsea method were used to determine immune cell infiltration levels in HCC samples. Based on Yin He et al.‘s research, a gene set with 47 immune checkpoints was collected [11]. By analyzing the Person correlation between the m6A.ES.benefit and m6A.ES.harm models, immune cell infiltration, and immune checkpoint expression, the association between m6A models and the tumor immune microenvironment was evaluated.

Cell lines

Hep3B and huh7 cell lines were provided by Shanghai outdo biotech company (Shanghai, China). All cell lines were characterized by short tandem repeat (STR) analysis by Genetic Testing Biotechnology Corporation (Suzhou, China).

RNA interference, RNA extraction and real-time quantitative PCR

The small interfering RNA (siRNA) oligonucleotides targeting METTL4 and MARCKSL1 were synthesized by Shanghai Outdo Biotech Company (Shanghai, China). Total RNA was extracted by using Trizol reagent (Sigma-Aldrich, Germany). First-strand cDNA was synthesized by using the PrimeScript™ RT Master Mix (Takara, Shiga, Japan) with random primers. Real-time quantitative PCR was performed with a ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China). All reactions were run in triplicate on the LightCycler^® 480II Real-Time PCR Detection System (Roche, Indianapolis, IN, USA). The β-actin was applied as endogenous control. The relative expression was calculated using the comparative Ct (2^−ΔΔCT) method.

Sequences of siRNA against specific target in this study:

	Sequences
METTL4 siRNA1 sense	5’-GGAGUUCACUACUUCUGUUTT-3’
METTL4 siRNA1 anti-sense	5’-AACAGAAGUAGUGAACUCCTT-3’
METTL4 siRNA2 sense	5’-CCAGCUGUUCAUAAAGAAUTT-3’
METTL4 siRNA2 anti-sense	5’-AUUCUUUAUGAACAGCUGGAG-3’
METTL4 siRNA3 sense	5’-AAGCCCUACGAAGGUAUUAUATT-3’
METTL4 siRNA3 anti-sense	5’-UAUAAUACCUUCGUAGGGCUUTT-3’
MARCKSL1 siRNA1 sense	5’-GCCUGUGAAUCAACUGUGUCUTT-3’
MARCKSL1 siRNA1 anti-sense	5’-AGACACAGUUGAUUCACAGGCTT-3’
MARCKSL1 siRNA2 sense	5’-GGCUAGUGCAGCCUCAGAATT-3’
MARCKSL1 siRNA2 anti-sense	5’-UUCUGAGGCUGCACUAGCCTT-3’
MARCKSL1 siRNA3 sense	5’-GGCCAACGGCCAGGAGAAUTT-3’
MARCKSL1 siRNA3 anti-sense	5’-AUUCUCCUGGCCGUUGGCCTT-3’
Control siRNA sense	5’-UUCUCCGAACGAGUCACGUTT-3’
Control siRNA anti-sense	5’-ACGUGACUCGUUCGGAGAATT-3’

Primers used in this study:

Primer names	Sequences
METTL4 forward	5’-TGGTCTGTGGAGGTAGTTGC-3’
METTL4 reverse	5’-CGGTGGCTTATGTGAGTGAA-3’
MARCKSL1 forward	5’-CAGAAAGGAAGGGGCTATT-3’
MARCKSL1 reverse	5’-GGGGAAGAAAACCAACAACT-3’
β-actin forward	5’-GAAGAGCTACGAGCTGCCTGA-3’
β-actin reverse	5’-CAGACAGCACTGTGTTGGCG-3’

Cell transfection

1 × 10⁶ cells per well were transfected with siRNAs using X-tremeGENE siRNA transfection reagent (Roche, Indianapolis, IN, USA) according to the manufacturer’s protocol. The cells were harvested 24 h after transfection.

In vitro cell proliferation

Cells transfected with METTL4 siRNA-2, MARCKSL1 siRNA-1, or the control siRNA were seeded into 96-well plates at a density of 5,000 cells/well and incubated for 24–48 h to facilitate growth. Subsequently, 10 µl of Cell Counting Kit-8 (CCK-8) from Vazyme (Nanjing, China) was added to each well, followed by incubation at 37 °C for 1 h. The absorbance of each well at 450 nm was measured using the EonTM Microplate Reader from BioTek (VT, USA). Three independent experiments were performed.

In vitro invasion assays

The transwell chambers were coated with a mixture of 60 µl Matrigel matrix (Corning, USA) diluted in an 8:1 ratio with serum-free DMEM medium. The lower chamber was filled with 600 µl DMEM medium containing 10% FBS. In the upper chamber, 10,000 cells transfected with METTL4 siRNA-2, MARCKSL1 siRNA-1, or control siRNA were seeded. Following an incubation period of 168 h at 37 °C, the invaded cells were fixed using 4% paraformaldehyde and stained with 1% crystal violet (Solarbio, Beijing, China) for 20 min. The average number of invading cells across three fields was recorded using an inverted microscope. Three independent experiments were performed.

Statistical analysis

For analysis based on public datasets, R (https://www.r-project.org/, v 4.2.2) was used for data processing, analysis, and visualization. Student t-tests or non-parametric Mann-Whitney U-tests were applied to detect continuous variable differences. Pearson correlation coefficients were computed using the cor.test function. Survival analysis was conducted using survival and survminer packages, while X-tile (v3.6.1, Yale University School of Medicine) determined the optimal cut-off value. The P-value threshold for statistical significance was set at 0.05, unless otherwise specified in the manuscript.

The mean ± standard error of the mean (SEM) values were used to represent the experimental data, which derived from three independent experiments. GraphPad Prism 9 software was utilized to conduct the statistical analysis. Student’s t-test or a nonparametric test was used to detect differences between two groups. Statistical significance was determined by P values less than 0.05. * P < 0.05, ** P < 0.01, *** P < 0.001.

Results

m6A regulators are expressed abnormally in HCC

To enhance the reliability of our findings, we performed a thorough literature search and collected 31 transcriptome datasets to study the expression of the 29 m6A regulators in HCC. From these datasets, 20 included both HCC and non-tumor samples. Our analysis indicated that most m6A regulators exhibit high expression in HCC across these 20 datasets, with factors such as METTL4, METTL3, CBLL1, IGF2BP2, and IGF2BP3 showing strong consistency (Fig. 1A–D and Table S3). We validated the protein levels of 10 well-consistent and retrievable m6A regulators using the Human Protein Atlas (HPA). Except for IGF2BP3, which is undetectable in both HCC and normal hepatocytes, all other 9 molecules show high expression in HCC (Fig. 1E, F and Figure S1). In addition, Single-cell RNA sequencing data analysis from GSE125449 revealed that the 29 m6A regulators are mainly expressed in malignant cells, tumor-associated endothelial cells (TEC), cancer-associated fibroblasts (CAF), and T cells, within the HCC tumor microenvironment (Fig. 1G, H).

Fig. 1

Expression patterns of m6A regulators in HCC. A, Differences in m6A regulator expression between HCC and non-tumor samples across 20 datasets. B, Differences in m6A regulator expression in GSE14520 − GPL3921. C–D, Expression differences of METTL4 and CBLL1 in GSE14520 − GPL3921. E–F, METTL4 and CBLL1 proteins expression in the human protein atlas. G, Analysis of m6A regulators expression in different cells using the single-cell sequencing gene dataset GSE202642. H, Single-cell sequencing analysis of IGF2BP2 expression in various cells. Down, under-expressed in HCC tissues; Up, over-expressed in HCC tissues. eraser, demethylase; reader, m6A binding protein; writer, methyltransferase. Exp, normalized gene expression value. TEC, tumor-associated endothelial cells. CAF, cancer-associated fibroblasts. TAM, tumor-associated macrophages. HPC-like, cells with an unknown entity but express hepatic progenitor cell markers. *P < 0.05, ***P < 0.001

Genetic features of m6A regulators in HCC

To investigate the genetic features of m6A regulators, we analyzed DNA copy variations and DNA methylation changes using TCGA data and assessed their impact on gene expression. Our analysis revealed that 6 m6A regulators undergo DNA copy gains in HCC, while 9 molecules display copy losses (Fig. 2A–D). Except for IGF2BP1, HNRNPG, HNRNPA2B1, IGF2BP3, and HNRNPC, all other 24 m6A regulators’ expression is positively correlated with copy number changes (Fig. 2E–G). DNA methylation analysis showed lower levels in HCC for 11 m6A regulators compared to non-tumor tissues, while 6 regulators have higher levels (Fig. 2H–J). Overall, our study suggested that the expression of most m6A regulators inversely correlates with their DNA methylation, with a few exceptions such as CBLL1, YTHDF1, ZC3H13, and ALKBH5 (Fig. 2K, L). Additionally, SNP analysis across five datasets indicated each m6A regulator’s mutation burden is relatively low (≤ 2%) in United States, France, and Japan. However, in the Chinese population, ZC3H13, KIAA1429, and RBM15 show higher mutation frequencies (4-8%) (Fig. 2M and Figure S2A–D). Furthermore, we examined the impact of mutations on gene expression in three datasets with both mutation and transcription information: TCGA, LICA-FR, and LIRI-JP. KIAA1429, ZC3H13, and LRPPRC were suitable for analysis as they exhibited mutations in more than three HCC samples and had corresponding mRNA expression values. The results indicate that mutations do not affect the expression of these three genes (Figure S2E).

Fig. 2

Genetic features of m6A regulators in HCC. A, CNV features of m6A regulators in TCGA. B, Differences in seg.mean of m6A regulator between HCC and non-tumor samples. C–D, Seg.mean of IGF2BP3 and IGF2BP2 in HCC and non-tumor samples. E, Relationship between CNV and gene expression of m6A regulators. F–G, Relationship between CNV and gene expression of ALKBH5 and KIAA1429. H, Differences in DNA methylation of m6A regulator between HCC and non-tumor samples. I–J, DNA methylation of ALKBH5 and ALKBH3 in HCC and non-tumor samples. K, Relationship between DNA methylation and gene expression of m6A regulators. L, R Association between DNA methylation and gene expression of ALKBH3. M, SNP features of m6A regulators in TCGA. Seg.mean, normalized segmentation mean value. CNV, Copy number variation. SNP, Single nucleotide polymorphism. Gain, copy gain variations in HCC tissues. Loss, copy loss variations in HCC tissues. Down, down-regulated in HCC tissues; Up, up-regulated in HCC tissues. *P < 0.05, **P < 0.01, ***P < 0.001

m6A regulators affect the progression and prognosis of HCC

To examine the association of m6A regulators with the progression and prognosis of HCC at a data level, we leveraged the GSVA method to evaluate the 12 HCC-related signatures of each tumor sample across 31 datasets. These signatures represent various features of HCC such as proliferation, metastasis, vascular invasion, tumor stemness, and prognosis (Table S2).

We estimated the correlation between m6A regulators and HCC-related signatures in each dataset (Table S4). By calculating the cumulative correlation of each m6A regulator with the 5 beneficial or 7 harmful signatures (Fig. 3A, B, and Figure S3A, B) as well as the composition of these correlations (Fig. 3C and Figure S3C), we were able to assess and compare their roles. Among them, regulators like IGF2BP2, IGF2BP3, and IGF2BP1 demonstrated a positive correlation with harmful signatures and a negative correlation with beneficial ones. Conversely, factors such as METTL14, ZC3H13, and FMR1 showed an opposite trend. For instance, IGF2BP2, the most prominent among the regulators, exhibited a positive correlation with proliferation, vascular invasion, metastasis, stem cell, and poor survival signals in almost all analyzable datasets, while being negatively correlated with beneficial signatures (Fig. 3D–G).

Fig. 3

The influence of m6A regulators on tumor progression and prognosis. A, Visualization of m6A regulators negatively correlated with beneficial HCC-related signatures across 31 datasets. The gene node size signifies the negative sum of the correlation between the corresponding gene and beneficial HCC-related signatures. B, Representation of m6A regulators positively correlated with beneficial HCC-related signatures. The gene node size represents the cumulative correlation between the corresponding gene and signatures. C, Display of the composition ratio of the correlation between m6A regulators and beneficial HCC-related signatures among 31 datasets. D, Relationship between IGF2BP2 and HCC-related signatures among 30 datasets. E, Association between IGF2BP2 and HCC-related signatures in GSE14520-GPL3921. F–G, Relationship of IGF2BP2 with proliferation inhibition and metastasis signatures in GSE14520-GPL3921. H, Correlation between nine harmful m6A regulators and overall survival among 14 datasets. I–J, Kaplan-Meier analysis of the correlation between IGF2BP2/IGF2BP3 and overall survival. K, Relationship between four beneficial m6A regulators and overall survival among 12 datasets. L–M, Association between METTL14/ZC3H13 and overall survival. Negative, negative correlation. Positive, positive correlation. No, absence of correlation. Cor, correlation coefficient. *P < 0.2 and |Cor|>0.1, **P < 1e-6 and |Cor|>0.1, ***P < 1e-11 and |Cor|>0.1

Additionally, we classified the association strength of m6A regulators with HCC signatures into strong, moderate, and weak levels based on the tertiles of the absolute value of cumulative correlation. Through an extensive literature review, we summarized the biological functions of each m6A regulator in HCC, providing cross-validation between data analysis and experimental findings (Table S5). The predicted functions of 10 strong associated regulators were highly consistent with previous research, except for METTL4, whose role in HCC remains unreported. However, among the 19 moderate and weak associated m6A regulators, the functions of YTHDF2, FTO, and ALKBH5 remain controversial, and the roles of YTHDC2, FMR1, ALKBH3, and YTHDF3 appear inconsistent with our analysis results.

Further, we analyzed the relationship between m6A regulators and overall survival across multiple datasets (Fig. 3H–M and Figure S3D). We found that the 9 strong associated harmful m6A regulators generally indicated a worse prognosis when highly expressed (Fig. 3H). Conversely, the 4 strong or moderate associated beneficial factors primarily suggested a better prognosis when expressed at high levels (Fig. 3K). Given that these 9 harmful m6A regulators and 4 beneficial factors are significantly associated with HCC-related signatures such as proliferation, metastasis, vascular invasion, tumor stemness, and prognosis, this suggests that they play important roles in HCC. Therefore, we selected these regulators for further analysis.

Identification and functional analysis of target genes of m6A regulators

In previous studies, 701 target genes controlled by 27 m6A regulators were discovered in various cells using techniques such as RIP and MeRIP. These findings are documented in the RM2Target database (http://rm2target.canceromics.org). To explore if these genes display similar regulatory patterns in HCC, we computed the correlation coefficients between each target gene and its associated regulators across 31 datasets (Table S6). Based on cumulative correlation strength(|cumulative correlation|>1), we pinpointed 91 potential target genes for 9 harmful m6A regulators (IGF2BP3, IGF2BP2, METTL4, HNRNPC, HNRNPA1, YTHDF1, IGF2BP1, HNRNPG, METTL16) and 29 potential targets for 4 beneficial m6A regulators (FMR1, METTL14, ZC3H13, YTHDC2), with most target genes showing a positive correlation with corresponding regulators (Fig. 4A, Figure S4A). Since there is only one target gene for METTL4 in the RM2Target database, its absolute cumulative correlation in 31 HCC datasets is less than 1. Therefore, Fig. 4A does not include information about METTL4.

Fig. 4

Downstream target genes analysis of m6A regulators in HCC. A, Visualization of the downstream target genes correlated with harmful m6A regulators across 31 datasets. The gene node size represents the cumulative correlation between the corresponding target gene and harmful m6A regulator. The red line indicates negative correlation and the blue line signifies positive correlation. B, Main KEGG enrichment-pathways of the target genes associated with harmful m6A regulators. C, Main GO enrichment-pathways of the target genes associated with harmful m6A regulators. D, Relationship between CDK4 and 3 harmful m6A regulators among 29 datasets. E–G, Correlation of CDK4 with HNRNPC, YTHDF1 and IGF2BP3 in LICA-FR. H, Single-cell sequencing analysis of CDK4 expression in various cells in GSE202642. I, CDK4 protein expression in the human protein atlas. J, Relationship between CDK4 and HCC-related signatures among 29 datasets. K–M, Association of CDK4 with vascular invasion inhibition, proliferation and metastasis signatures in GSE112790. Cor, correlation coefficient. Negative, negative correlation. Positive, positive correlation. TEC, tumor-associated endothelial cells. CAF, cancer-associated fibroblasts. TAM, tumor-associated macrophages. HPC-like, cells with an unknown entity but express hepatic progenitor cell markers. *P < 0.2 and |Cor|>0.1, **P < 1e-6 and |Cor|>0.1, ***P < 1e-11 and |Cor|>0.1

Using GO and KEGG analyses, we found that these potential target genes primarily participate in the progression of malignant tumors like HCC, colorectal cancer, and gastric cancer. They are notably linked to Hepatitis B/C infection, p53, PI3K-Akt, and other cancer-related signaling pathways, and are involved in biological processes including cell proliferation, differentiation, cell cycle, and radiation response (Fig. 4B, C, Figure S4B, S4C, and Table S7–10).

Among all the target genes, CDK4 and MARCKSL1 have the highest absolute cumulative correlation with m6A regulators. Therefore, we further analyzed these two genes. Across multiple datasets, both genes almost always display a positive correlation with their regulators (Fig. 4D–G, Figure S4D–G) and are mainly upregulated in tumor samples (Figure S4H, Table S11). Single-cell analysis revealed that CDK4 and MARCKSL1 are expressed in immune cells, tumor cells, and tumor-associated cells, with notably higher expression levels observed in TECs and CAFs (Fig. 4H, Figure S4I). Further validation in the HPA database demonstrated that CDK4 and MARCKSL1 have higher expression levels in HCC cells compared to hepatocytes (Fig. 4I, Figure S4J).

Moreover, similar to their regulators, CDK4 and MARCKSL1 mainly exhibit positive correlations with harmful HCC-related signatures and negative correlations with beneficial ones (Fig. 4J–M, Figure S4K, and Table S12). Increased expression of both targets implied a worse prognosis (Figure S4L–N).

Functional validation of METTL4 and MARCKSL1

In the previous analysis, we conducted big data analysis to examine the functions of 10 strong associated regulators in HCC. Except for METTL4, whose role in HCC remains unexplored, the predicted functions of the other regulators were highly consistent with previous research (Table S5). Furthermore, our data analysis indicated that two major target genes of m6A regulators, CDK4 and MARCKSL1, may contribute to HCC progression (Fig. 4J–M, Figure S4K–N). The oncogenic function of CDK4 has been confirmed by Wan et al. [12], while the role of MARCKSL1 in HCC has not been reported. To further validate the predictive capability of our data analysis, we conducted cell experiments to investigate the effects of METTL4 and MARCKSL1 on the proliferation and invasion of HCC cells. Using RNA interference technology, we reduced the expression of METTL4 in Hep3B and Huh7 cells (Fig. 5A, B). CCK-8 revealed that downregulation of METTL4 could inhibited the proliferation of HCC cells (Fig. 5C, D), and transwell invasion assays demonstrated that downregulation of METTL4 could suppress the invasive ability (Fig. 5E). Similarly, we demonstrated that downregulation of MARCKSL1 could inhibit the proliferation and invasion of HCC cells (Fig. 5F–J). These experimental results confirm the oncogenic functions of METTL4 and MARCKSL1 in HCC, which are consistent with our big data analysis results.

Fig. 5

The effect of METTL4 and MARCKSL1 on tumor progression. A–B, SiRNA down-regulates the METTL4 expression in Hep3B and Huh7 cells. C-D, CCK-8 assays for Hep3B and Huh7 cells transfected with METTL4 siRNA-2 or the control. E, Transwell invasion assays for Hep3B and Huh7 cells transfected with METTL4 siRNA-2 or the control. Scale bars = 100 μm. F-G, SiRNA down-regulates the MARCKSL1 expression in Hep3B and Huh7 cells. H-I, CCK-8 assays for Hep3B and Huh7 cells transfected with MARCKSL1 siRNA-1 or the control. J, Transwell invasion assays for Hep3B and Huh7 cells transfected with MARCKSL1 siRNA-1 or the control. Scale bars = 100 μm. Data are presented as mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001. NS, Not Significant

Construction of risk models with m6A regulators and prognostic impact analysis

In the preceding analysis, we identified 9 harmful m6A regulators and 4 beneficial factors that are closely related to the progression and prognosis of HCC. To explore their potential clinical application value, we constructed m6A.ES.harm and m6A.ES.benefit models separately in the analyzable datasets through the GSVA method.

The correlation composition ratio indicates a robust association between the two m6A models and HCC-related signatures (Fig. 6A, Table S13). m6A.ES.harm is predominantly positively correlated with harmful signatures, and primarily displays a negative correlation with beneficial signatures (Fig. 6B–E). Conversely, m6A.ES.benefit exhibits the reverse trend (Fig. 6F–I). Survival analysis suggests that high m6A.ES.harm scores mainly predict a worse prognosis, while high m6A.ES.benefit scores typically point towards a better prognosis (Fig. 6J–L, Figure S5A–B). Additionally, both scores are strongly correlated with clinical characteristics such as tumor doubling time, tumor size, the grade of tumor malignancy, clinical pathological staging, vascular invasion, and Alpha-Fetoprotein (AFP) levels (Fig. 6M–P, Figure S5C–D). Furthermore, we employed ROC analysis to compare the prognostic value of two m6A models with the traditional HCC biomarker AFP (Figure S6). The results demonstrated that the AUC values of m6A.ES.harm at 2- and 5-year OS are superior to those of AFP (0.638 vs. 0.572 and 0.682 vs. 0.615, respectively). The 2-year AUC value of m6A.ES.benefit is higher than that of AFP (0.618 vs. 0.572), but its 5-year AUC value is lower than that of AFP (0.594 vs. 0.615).

Fig. 6

Correlation of m6A models with tumor progression and prognosis. A, Composition ratio of the correlation between m6A models and HCC-related signatures. B, Relationship between m6A.ES.harm and HCC-related signatures among 31 datasets. C-E, Association of m6A.ES.harm with proliferation, poor survival and vascular invasion inhibition signatures in GSE25097. F, Correlation between m6A.ES.benefit and HCC-related signatures among 29 datasets. G–I, Relationship of m6A.ES.benefit with proliferation inhibition, metastasis and good survival signatures in GSE25097. J. The impact of m6A.ES.harm and m6A.ES.benefit on overall survival among 11 datasets. K–L, Survival analysis among two groups stratified by the m6A.ES.benefit and m6A.ES.harm in GSE14520-GPL3921. M–N, Differences of m6A.ES.harm among fast/slow tumor growth groups and grade subgroups. O–P, Differences of m6A.ES.benefit among tumor size and TNM subgroups. Negative, negative correlation. Positive, positive correlation. No, absence of correlation. Cor, correlation coefficient. **P < 0.01, ***P < 0.001

Assessment of drug sensitivity based on m6A models

To explore the potential of m6A models in predicting drug response, we conducted a drug sensitivity analysis. In 24 HCC microarray datasets, we used the pRRophetic method to calculate the correlation between more than 200 drugs and the two m6A models (Table S14). The top 100 drugs with the highest absolute cumulative correlation with either m6A.ES.harm or m6A.ES.benefit are presented in Fig. 7A. Additionally, we used the oncoPredict method to compute the association between over 500 drugs and the two models in seven HCC sequencing datasets (Table S15).

Fig. 7

Drug sensitivities analysis based on m6A models. A, Network diagram displaying drugs with significant correlations between sensitivity and m6A.ES.harm or m6A.ES.benefit across 24 datasets through pRRophetic method. The gene node size represents the cumulative correlation between the corresponding drug sensitivity and m6A.ES.harm or m6A.ES.benefit. B, Correlation of m6A.ES.harm with sensitivity of 12 anti-cancer drugs through pRRophetic method. C–D, Relationship of m6A.ES.harm with gemcitabine and sorafenib sensitivity in GSE36376. E, Validation of the correlation between m6A.ES.harm and sensitivity of 12 anti-cancer drugs through oncoPredict method among 7 datasets. F–G, Association of m6A.ES.harm with gemcitabine and sorafenib sensitivity in TCGA. H, Relationship between m6A.ES.benefit and sensitivity of 4 anti-cancer drugs through pRRophetic method. I–J, Correlation of m6A.ES.benefit with paclitaxel and bortezomib sensitivity in GSE25097. K, Validation of the correlation between m6A.ES.benefit and sensitivity of 4 anti-cancer drugs through oncoPredict method. L–M, Association of m6A.ES.benefit with paclitaxel and bortezomib sensitivity in GSE114564. N–O, Differences of m6A.ES.harm and m6A.ES.benefit among TACE subgroups. Negative, negative correlation. Positive, positive correlation. Cor, correlation coefficient. TACE, trans-catheter arterial chemoembolization. ***P < 0.001

Through these two drug sensitivity analysis methods, we found that 12 drugs were significantly negatively correlated with m6A.ES.harm, including gemcitabine, sorafenib, and mitomycin C (Fig. 7B–G). This suggests that a lower m6A.ES.harm score may benefit from the treatment of these drugs. Besides, our analysis also indicates that m6A.ES.benefit was mainly positively correlated with drugs like paclitaxel, bortezomib, PI.103, and BI.2536 (Fig. 7H–M). This implies that patients with a higher m6A.ES.benefit score may derive therapeutic benefits from the four agents.

Furthermore, in the GSE104580 dataset, the m6A.ES.harm score was lower and the m6A.ES.benefit score higher in the TACE response group (Fig. 7N, O). Considering TACE in GSE104580 was based on doxorubicin, mitomycin C, carboplatin, and fluorouracil [13,14,15], and our analysis also showed a strong correlation of mitomycin C with m6A.ES.harm, this proves the value of our multi-dataset drug sensitivity analysis approach using pRRophetic and oncoPredict.

Analysis of immune microenvironment and immune response based on m6A models

Previous studies have demonstrated the involvement of m6A modification in tumor immunity [16]. To investigate the impact of m6A.ES.harm and m6A.ES.benefit on the tumor immune microenvironment (TIM), we analyzed their associations with immune cell infiltration and immune checkpoint expression (Table S16, S17, S18).

CIBERSORT analysis revealed that m6A.ES.harm is positively correlated with M0 macrophages and follicular helper T cells, and negatively correlated with resting mast cells and gamma delta (γδ) T cells (Fig. 8A–D). ssGSEA analysis showed that m6A.ES.harm is negatively correlated with effector memory CD8 T cells and type 1 T helper cells, and positively correlated with activated CD4 T cells and type 2 T helper cells (Fig. 8E–H). Moreover, m6A.ES.harm exhibits positive correlation with immune checkpoints TNFSF4 and CTLA4, and negative correlation with IDO2 and CD160 (Fig. 8I–L). Additional associations of m6A.ES.harm with immune cells and immune checkpoints can be found in the figures. m6A.ES.benefit also displays associations with immune cells, such as M0/M1 macrophages, γδT cells, regulatory T cells, and activated CD4 T cells, as well as immune checkpoint expression including CD274(PD-L1), CD40, and TNFSF14 (Figure S7A–F). These results imply that m6A.ES.harm and m6A.ES.benefit may affect immune cell infiltration and immune checkpoint expression in TIM of HCC.

Fig. 8

Immune infiltration and immunotherapy prediction based on m6A models. A, Correlation of m6A.ES.harm with 22 tumor-infiltrating immune cell types among 31 datasets through CIBERSORT method. B, Volcano plot of four immune cell types significantly associated with m6A.ES.harm. C–D, Correlation of m6A.ES.harm with M0 macrophage and gamma delta T cell in GSE104580. E, Relationship between m6A.ES.harm and 28 infiltrating immune cell types through ssGSEA method. F, Heat-map of eight immune cell types significantly associated with m6A.ES.harm. G–H, Correlation of m6A.ES.harm with activated CD4 T cell and CD56dim natural killer cell in GSE124751. I, Relationship between m6A.ES.harm and expression of immune check points. J, Heat-map of eight immune check points significantly associated with m6A.ES.harm. K–L, Association of m6A.ES.harm with the expression of TNFSF4 and IDO2 in TCGA. M–N, Differences of response to PD-1 or CTLA4 blocker among low and high m6A.ES.harm groups in TCGA. O–P, Differences of response to PD-1 or CTLA4 blocker among low and high m6A.ES.benefit groups. Negative, negative correlation. Positive, positive correlation. Cor, correlation coefficient. Strong, immunophenoscore (IPS) = 10. Moderate, 8 ≤ IPS ≤ 9. Weak, IPS ≤ 7. *P < 0.05, **P < 0.01, ***P < 0.001. NA, not available

To explore the influence of m6A.ES.harm and m6A.ES.benefit on the efficacy of immune checkpoint inhibitors, we divided them into high and low subgroups based on the optimal cut-off values from survival analysis and assessed immune response differences between the subgroups when treated with PD-1 or CTLA4 blocker. The results showed that the low m6A.ES.harm subgroup exhibited a more pronounced response to treatment with PD-1 blocker or CTLA4 blocker (Fig. 8M, N). The high m6A.ES.benefit subgroup showed a favorable response to treatment with CTLA4 blocker, while no significant differences were observed between different m6A.ES.benefit subgroups when treated with PD-1 blocker (Fig. 8O, P).

Discussion

HCC is a highly heterogeneous malignant tumor [17], which makes it difficult to identify robust and clinically meaningful biomarkers. Previous studies have demonstrated that m6A regulators play critical roles in cancer [6, 18, 19]. However, their functions in HCC are controversial. For example, the impact of the m6A reader YTHDF2 on HCC growth and metastasis is inconsistent across different studies [20,21,22]. By integrating big data analysis with basic experiments, we can accurately identify the function of key biomarkers, as demonstrated in our previous research [23]. In this study, we aim to explore the functional roles of m6A regulators and their potential clinical applicability by combining large data analysis with basic research.

Multi-omics data analysis reveals that m6A regulators such as METTL4, METTL3, CBLL1, IGF2BP2, and IGF2BP3 are highly expressed in HCC and demonstrate high consistency in different datasets, suggesting their potential for HCC diagnosis. Although certain molecules such as ZC3H13, METTL4, and YTHDC1 exhibit low expression in HCC in some datasets, they show opposite results in other datasets, indicating that the expression patterns of these molecules in HCC remain controversial (Fig. 1A). Most regulators show predominant expression in malignant cells, TEC, CAF, and T cells within the tumor microenvironment (Fig. 1G), which are closely related to the immune regulation, angiogenesis, invasion, and metastasis in HCC [24,25,26,27], implying that m6A regulators may play a role in managing these biological processes. CNV gain typically leads to an elevation in gene expression [28], aligning with our findings particularly for genes ALKBH5, WTAP, and KIAA1429 (Fig. 2E). DNA methylation usually inhibits gene transcription [29]. Although our results exhibit this trend, the correlation is not robust (Fig. 2K). This outcome might be attributed to our data source, Illumina HumanMethylation450, which does not capture all methylation sites. Additionally, besides CNV and DNA methylation, gene expression is influenced by multiple factors, including lncRNAs, acetylation, and phosphorylation [30]. This complexity might explain why, in the TCGA set, YTHDF2 displays CNV loss and increased DNA methylation in tumor samples, but its expression levels actually rise (Figs. 1A and 2B, and Fig. 2H). ZC3H13, KIAA1429, and RBM15 exhibit higher mutation frequencies in the Chinese HCC population (Figure S2A), and prior studies have noted their involvement in HCC proliferation and metastasis [18, 31, 32]. Whether these variations will induce therapeutic differences in Chinese HCC patients is a question that warrants further exploration.

To facilitate the description of correlations among quantitative data across multiple datasets, we introduce the concept of cumulative correlation. This approach examines quantitatively the association between m6A regulators and HCC-related signatures, such as proliferation, vascular invasion, and metastasis, thus identifying key factors. In our research, we discovered a strong correlation between 10 m6A regulators and HCC-related signatures, aligning well with the results of fundamental research (Figs. 3 and 5C–E, Figure S3, and Table S5). However, regulators with weak correlations presented disagreements in basic research. To some extent, this can be rationalized as follows: factors that have a strong correlation are more likely to be confirmed in vivo and in vitro studies, while factors with weak influence may be more susceptible to interference from various experimental variables, leading to inconsistent outcomes. In summary, calculating cumulative correlation through large-scale datasets and mutually verifying with basic experiments may be an effective approach to precisely define biomarker functions.

Upstream genes frequently exert biological effects by modulating the expression of downstream genes [30]. In this research, we pinpointed potential targets of m6A regulators using cumulative correlation. These targets are primarily involved in oncogenic pathways and biological processes in various cancers, including HCC (Fig. 4A-C, Figure S4A–C), further corroborating the role of m6A modification in HCC. The cell cycle regulatory protein CDK4, a well-known oncogene, has been confirmed to promote the progression of HCC in previous research [12]. Additionally, its inhibitors have demonstrated antitumor activity against HCC [33]. In our analysis, we discovered a strong correlation between CDK4 and three m6A readers: HNRNPC, YTHDF1, and IGF2BP3 (Fig. 4D–G). Yinmin Gu et al. also confirmed that IGF2BP3 can stabilize CDK4 expression in an m6A-dependent manner [34]. Hence, we infer that HNRNPC, YTHDF1, and IGF2BP3 could enhance CDK4 expression in HCC via m6A modification, thereby realizing their oncogenic role. This is consistent with our findings that all four display harmful properties (Figs. 3 and 4J–M, Figure S4L–M). As for the target gene MARCKSL1, which exhibits a strong correlation with IGF2BP1, IGF2BP2, IGF2BP3, previous studies have identified its role in the proliferation and metastasis process of lung adenocarcinoma and thyroid carcinoma [35, 36]. Our analysis suggests that MARCKSL1, similar to its three potential regulatory genes, plays an oncogenic role in HCC (Figs. 3 and 5H–J, Figure S4D–J, K, L, N). In summary, our target gene analysis provides further insight into the biological roles of m6A regulators and offers data support for their mechanistic exploration in HCC.

To further explore the clinical implications of m6A regulators in HCC, we constructed two molecular models: m6A.ES.harm and m6A.ES.benefit. Both models are significantly associated with HCC-related signatures and clinical features (Fig. 6, Figure S5). Drug sensitivity analysis suggests that patients with low m6A.ES.harm scores might benefit from medications such as gemcitabine, sorafenib, mitomycin C, while those with high m6A.ES.benefit scores might respond better to treatments with paclitaxel and botezomib (Fig. 7). Previous studies have shown that the disease-control rate is about 40% for advanced HCC treated with sorafenib [37], and the overall response rate ranges from 2.1 to 17.8% for patients treated with gemcitabine alone [38, 39]. Research by Ohnishi K et al. found that chemoembolization with mitomycin C microcapsules is a useful palliative measure for the treatment of inoperable HCC [40]. Paclitaxel, bortezomib, temozolomide, cytarabine, etoposide have also been shown to be effective in certain HCC patients [41,42,43,44,45]. However, except for sorafenib, prescribed as a frontline therapy for advanced HCC, other medications have been increasingly disregarded because of their limited efficacy. Pharmaceutical development is an expensive, lengthy, and complex process. Identifying suitable patients through molecular subtype for impactful drug intervention carries substantial clinical and financial importance. In this context, our m6A molecular models might provide new opportunities to repurpose these drugs for HCC treatment. Furthermore, our analysis suggests that the efficacy of innovative drugs such as BI2536 and PI-103 correlates with m6A.ES.harm and m6A.ES.benefit. Studies have demonstrated that BI2536 can reverse multidrug resistance in hepatoma cells [46], and PI-103 notably inhibits Huh7 proliferation [47]. Both compounds exhibit potential for HCC treatment. Therefore, our models might aid in designing personalized treatment approaches for these novel drugs, though clinical trial validation is necessary.

Our preliminary analysis indicated that the majority of m6A regulators are notably expressed in tumor-associated immune cells, including TEC, CAF, and T cells (Fig. 1G). Therefore, we explored the immune microenvironment based on the m6A models. The results indicate that m6A.ES.harm and m6A.ES.benefit are associated with the infiltration of immune cells such as M1 macrophages, natural killer T cells, MDSC, and regulatory T cells, and also with immune checkpoints including CTLA4 and PD-L1 (Fig. 8, Figure S7). Tumor immune regulation is a complex process involving immune cells and molecules. In HCC, the role of many immune cells and checkpoints remains unclear due to limited research, which makes data interpretation challenging. For instance, m6A.ES.benefit inversely correlates with regulatory T cells, MDSC, and natural killer T cells (Figure S7D), known mediators of immune suppression in HCC [48,49,50]. This may suggest a tumor microenvironment with high m6A.ES.benefit scores has less infiltration of these pro-carcinogenic immune cells, thus potentially improving patient prognosis. However, m6A.ES.benefit also correlates positively with M1 macrophages and PD-L1 (Figure S7B, F), both known to suppress immune responses in HCC and exhibit pro-tumor effects [51], which presents an interpretation conundrum. Therefore, we recommend a comprehensive analysis that considers the infiltration ratio and activity of different immune cells, as well as the expression of various immune checkpoints, to determine the intensity of the immune response. For example, in a high-m6A.ES.benefit environment with abundant immune-suppressive and activating cells, a predominance of immune activating cells may indicate stronger anti-tumor immunity. Of course, such integrated assessments necessitate further experiments to accurately define the functions of different immune cells and checkpoints in HCC. In addition, we assessed response differences to PD1/PDL1 and CTLA4 blockers among various m6A model subgroups. The results suggest that the low m6A.ES.harm and high m6A.ES.benefit subgroups exhibited stronger immune responses. Currently, the use of immune checkpoint inhibitors is a primary systemic approach for treating advanced HCC. PD1/PDL1 or CTLA4 inhibitors such as sintilimab, atezolizumab, tislelizumab, pembrolizumab, and tremelimumab have been shown to improve the efficacy and extend survival in patients with advanced HCC. However, the overall objective response rate is about 15-30%52–56, meaning most patients do not benefit from immunotherapy. The m6A models may help clinicians identify patients who are likely to respond to immunotherapy, thereby improving treatment efficacy and avoiding futile treatments.

In previous analysis, we combined big data with basic research and found that many m6A regulators are closely related to HCC proliferation, invasion, and prognosis. Based on these key factors, we constructed m6A models. Additionally, through drug sensitivity analysis and immune-related analysis, we discovered that these m6A models may have important clinical significance in guiding chemotherapy, targeted therapy, and immunotherapy for HCC. Indeed, beyond HCC, recent studies have widely demonstrated that m6A modifications are involved in the development, drug resistance, and immune regulation of various tumors. For instance, YTHDF1, IGF2BP1, and IGF2BP2 can regulate the growth and metastasis of head and neck squamous carcinoma, esophageal squamous cancer, ovarian cancer, or breast cancer [57,58,59,60]. Furthermore, FTO, YTHDF1, and METTL3 are associated with resistance to chemotherapy, targeted therapy, and immunotherapy in lung cancer, pancreatic cancer, colorectal cancer, breast cancer, or thyroid cancer [61,62,63,64,65,66]. Therefore, identifying truly meaningful and critical m6A regulators in different tumors using a ‘big data + basic research’ approach, constructing molecular subtypes based on this information, and thereby assessing the prognosis of different cancer patients and guiding treatment plans hold potential for future development.

Certainly, this study has some limitations. Our data analysis primarily relied on dozens of public datasets, each exhibiting variations in ethnicity, gender, age, treatment choices, and disease stages. These differences may introduce interference into our analytical results. Additionally, these datasets were sourced from different detection platforms and technologies. Currently, high-throughput detection technologies, including next-generation sequencing and gene chips, have inherent limitations in accuracy and reproducibility. Even when analyzing the same sample with the same technology, variations in results can occur. Moreover, numerous methods exist for data analysis, and different methods may lead to different outcomes. Furthermore, at the data level, this study analyzed the downstream regulatory genes of m6A regulators, as well as the relationship between m6A models and the tumor immune microenvironment and drug sensitivity. These findings require further validation through basic and clinical research.

Conclusions

In this study, we emphasized the importance of integrating big data analysis with experiment validation to discover effective biomarkers. Additionally, we proposed a novel method called cumulative correlation to quantify the relationships between disparate variables in multi-datasets. Moreover, by using m6A methylation modification as an entry point, we employed big data to analyze the expression patterns and genetic features of m6A regulators in HCC from various perspectives. In conjunction with research findings, we redefined the function of m6A regulators, and built m6A risk models to provide direction and data support for exploring individualized chemotherapy, targeted therapy, and immunotherapy strategies in HCC.