Raw gene expression data from HCC and the corresponding clinical information were extracted from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Transcript data from 374 tumor specimens and 50 normal specimens were included. HTseq counts were normalized based on the transcripts per million (TPM) method.

The exclusion criteria were as follows: (i) histologic diagnosis ruled out HCC; (ii) extremely low gene expression values; and (iii) incomplete clinical data and a follow-up time<30 days. ID conversion was conducted by Perl (Perl Programming Language version 5.30.1)

To generate the m⁶A-related signature and quantify the m⁶A modification patterns of each patient, we aimed to construct a scoring system (called m6AScore) and assess all individuals with HCC.

Initially, we conducted consensus clustering analysis using the ConsensusClusterPlus package (21) in R statistical software (version 4.1.3, https://www.r-project.org/). Steps included data preparation with gene expression data from TCGA or other sources, applying preprocessing techniques, selecting the ConsensusClusterPlus algorithm, determining the optimal number of clusters using the scree plot, and interpreting cluster results through gene functional analysis, pathway enrichment, and correlation analysis with clinical features.

Then, the overlapped differentially expressed genes (DEGs) were obtained by Venn diagrams. Next, the prognostic value of protein-coding DEGs of the distinct m⁶A modification patterns was determined using Cox proportional hazards regression analysis. Moreover, the protein-coding DEGs with significant prognostic value were included in calculating the m6AScore. Given the prognostic protein-coding DEGs in HCC, principal component analysis (PCA) was then performed to establish the m6AScore. PCA is a dimension reduction method and has been extensively used in gene expression analysis (19, 22). Similar to a previous study, we added principal components (PC) 1 and 2 as the last gene signature scores.

Where i and j are the order and the total number of m⁶A-related prognostic protein coding in HCC. Finally, the Z-score of the m6AScore was used for further analysis.

For different clusters, the outstanding gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified and visualized by bubble diagrams.

In total, an independent cohort of 101 patients with HCC was enrolled at Huashan Hospital, Fudan University, between March 2015, and June 2020. None of the patients received radiotherapy or chemotherapy before surgery. Each patient provided informed consent. The study protocol was approved by the Human Research Ethics Committee of Huashan Hospital, Fudan University. Basic information about the samples is listed in the Supplementary Materials ( Tables S4, S5 ).

Copy number variation (CNV) data were downloaded from the TCGA database. Gene mutation data were also downloaded from the TCGA database, and gene mutations were identified in the high-/low-m6AScore groups. Moreover, the tumor mutation burden (TMB) of each sample was calculated via a Perl script (https://www.perl.org/). TMB was defined as the total number of mutations per megabase in the tumor tissue.

GSEA is a computational method that can execute GO and KEGG analyses with a given gene list. GSEA software version 4.2.3 (https://www.gsea-msigdb.org/) was employed to predict the potential functions of the m⁶A-related signature. Combined with the high- and low-risk groups determined by the m6AScore, GO and KEGG pathway enrichment analyses were conducted to visualize various genes involved in different pathways, biological functions, and their expression patterns. Data were corrected for multiple testing (number of permutations = 1000).

The R package of pRRophetic was used to predict the sensibility of common chemotherapeutic agents, as previously described (23). Sensibility indicates the effectiveness of a substance in inhibiting specific biological or biochemical functions. The group difference was tested by the Wilcoxon signed-rank test.

Gene signature analysis of T-cell dysfunction and prediction of cancer immunotherapy response by patients with cancer were performed using the tumor immune dysfunction and exclusion (TIDE) algorithm as previously described (24). TIDE utilized both T-cell dysfunction and exclusion signatures to model immune escape in tumors with different cytotoxic T lymphocyte levels, and the TIDE score is consistent with signatures of tumor immune evasion. A higher tumor TIDE score is associated with a worse immune checkpoint blockade response.

By performing single-sample gene set enrichment analysis (ssGSEA) (25), the immune function scores of each patient with HCC were calculated to further quantify the composition of tumor-infiltrating immune cells. According to the immune score, the degree of immune cell infiltration was quantified in HCC tissues. Immunocyte gene signatures enrolled in this study are listed in Table S6 .

Proteomic detection was performed in four steps: protein isolation from formalin-fixed paraffin-embedded tumor samples, protein digestion, peptide sequencing by MS, and data analysis to identify the proteins.

The details of the processing of the raw files generated by MS are as follows: MS raw files were searched in the NCBI Human Refseq database using the Mascot search engine (v2.3, Matrix Science Inc.; version 04-07-2013, total 32015 entries). The parent and daughter ion mass deviations were set at 20 and 50 ppm, respectively (QExactive HF) or 0.5 Da. The theoretical protein cleavage sites were arginine (R) and lysine (K), with a maximum of two missed cleavage sites allowed. Carbamidomethylation (C) for immobilization and acetyl (protein N-term) and oxidation (M) for dynamic modifications. All the identified peptides were obtained from the area under the MS1 peak calculation, and the peptide false discovery rate was controlled at 1%. At the protein characterization level, we only kept those proteins that contained at least one unique peptide and two high-quality peptides (strict peptide, i.e., Mascot Ion Score >20). For protein quantification, we used an intensity-based absolute quantification algorithm (i.e., iBAQ algorithm) and subsequently normalized each iBAQ value to FOT-iBAQ values. This was performed by dividing the iBAQ value by the sum of the iBAQ of all detected proteins in the corresponding sample. Then, to make this FOT size easier to read and write, this FOT value was multiplied by 105, resulting in the iFOT value.

The proteome sequencing of 101 HCC specimens obtained from the General Surgery Department of Huashan Hospital, Fudan University, was performed by the National Institute of Metrology (China) as described above. The sample quality control and the protein profile expression involved in this study have been shown in Figures. S1, S2 in Supplementary Material .

After the paraffin sections were dewaxed in water, the slides were immersed in citric acid antigen repair buffer and boiled in water for antigen repair. The sample was incubated with 3% H₂O₂ at room temperature before rinsing with phosphate-buffered saline (PBS). The solution was closed with 5% FBS at 37°C for 30 min and then incubated with the primary antibody at 4°C overnight. The sample was rinsed with PBS, incubated with Cy5/Sp Green-labeled secondary antibodies at 37°C, and then rinsed again with PBS. The solution was then incubated with Hoechst at room temperature for 15 min, and serial staining was performed repetitively by stripping off previous primary/secondary antibodies via microwave treatment.

The Kaplan–Meier analysis with the log-rank test was performed to compare the survival difference between the high-/low-risk HCC groups. Data were tested for normal distribution, and if normally distributed, the Student’s t-test was used or the nonparametric Mann–Whitney U-test was applied unless stated otherwise. Correlations were assessed according to the Pearson test for parametric data and the Spearman test for nonparametric data. All statistical analyses were conducted using Prism (version 8, https://www.graphpad-prism.cn/) and R statistical software (version 4.1.3, https://www.r-project.org/). A P-value of<0.05 denoted a statistically significant difference (*P< 0.05; **P< 0.01, and *** P< 0.001).

The TCGA database, which includes detailed clinical information on 33 cancer types, is regarded as the landmark of cancer genomics programs. To identify the m⁶A-related signature in HCC, we downloaded all the raw gene expression data of 374 patients with HCC from the TCGA database. After carefully checking the clinical information of each sample, 31 patients were excluded because of incomplete clinical data or a follow-up time of less than 30 days. Ultimately, 343 patients with reliable transcript data and detailed clinical information were enrolled in the study. Through literature reports, 15 well-established m⁶A regulatory genes (METTL3, METTL14, WTAP, ZC3H13, RBM15, YTHDC1, YTHDF1, YTHDF2, YTHDF3, IGFBP1, IGFBP2, IGFBP3, RBMX, FTO, and ALKBH5) were used for further analysis. The whole analysis process was summarized as a flowchart in Figure 1 .

Univariate Cox regression analysis of the gene expression level revealed that 11 of the 15 m⁶A regulators were significantly correlated with HCC prognosis, whereas METTL14, IGFBP1/2, and ALKBH5 failed to reach statistical association ( Figure 2A ). According to the consensus algorithm with the 11 m⁶A regulators, the slope of the scree plot sharply declined when the abscissa was k = 3, suggesting that the patients could be divided into three categories for better cluster presentation ( Figure 2B ). GO and KEGG analyses showed that each cluster displayed different features ( Figure 2C ). Cluster A was closely associated with RNA processing and splicing ( Figure 2C , left panel). Cluster B was involved in histone modification and ribosome assembly ( Figure 2C , middle panel). Cluster C was mainly enriched in extracellular structure, collagen matrix, and integrin binding ( Figure 2C , right panel). There were 85 DEGs overlapped in the three clusters ( Table S1 ), and 68 of these DEGs were significantly associated with HCC prognosis ( Figure 2D , left panel). DEGs are listed in Table S2 . Given that protein-coding genes are major executors of biological function, we performed proteomic detection using MS assays on 101 HCC samples at our center. In total, 7059 protein expression were successfully detected using this approach. Subsequently, 68 prognosis-related genes were intersected with MS-detected genes, and 20 protein-coding genes were eventually defined as m⁶A-related signatures for analysis ( Figure 2D , right panel). Their coefficients are listed in Table S3 .

The expression levels of the 20 m⁶A-related genes were compared between HCC tissues and adjacent normal counterparts. Most of them were upregulated, and five of them were downregulated ( Figure 2E ). The copy number variation (CNV) of loss and gain is shown in Figure 2F . Moreover, we found an abundance of high-confidential m⁶A clusters within these 20 genes via SRAMP online analysis (http://www.cuilab.cn/sramp) ( Figure S3 ), suggesting that these genes are m⁶A modification targets.

Taken together, we established an m⁶A-related signature composed of 20 protein-coding genes, which has a significant correlation with HCC prognosis.

According to the m⁶A-related signature, the m6AScore was calculated for each patient, as described in the method. Based on the Youden index (26), 343 HCC patients in TCGA were divided into low- and high-risk groups. As shown in Figures 3A, B . the death risk gradually increased with the elevation of m6AScore. Meanwhile, a high m6AScore was markedly correlated with advanced HCC stages ( Figure 3C ). Thus, patients with high m6AScores were stratified into the high-risk HCC group.

Then, Kaplan–Meier analysis for overall survival (OS) showed that the high-risk group frequently had a shorter survival time ( Figure 3D ). The 1-, 3-, and 5-year survival rates of high-risk vs. low-risk patients were 62.2% vs. 91.7%, 40.1% vs. 70.7%, and 21.9% vs. 57.0%, respectively. Receiver operating characteristic (ROC) curve analysis indicated that the m6AScore was a useful prognostic index for OS prediction of patients with HCC, as the areas under the ROC curve (AUCs) reached 0.746 and 0.669 in 1 and 3 years, respectively ( Figure 3E ). Moreover, the m6AScore model displayed good calibration, and the calibration curves for 1- and 3-year survival were all close to the ideal line ( Figure 3F ).

To further confirm the reliability of our model, we conducted validation using an external dataset of ICGC-JP. As shown in Figure S4 , the 1-, 3-, and 5-year survival rates of high-risk vs. low-risk patients were 84.8% vs. 94.9%, 71.2% vs. 84.5%, and 40.7% vs. 70.0%, respectively. The 1- and 3-year AUC was 0.707 and 0.739. These results indicated that our model had a favorable predictive activity.

Collectively, these findings suggest that the m⁶A-related coding-gene signature has a great capacity to predict the OS of patients with HCC.

To evaluate the reliability of the m⁶A-related protein-coding gene signature in prognostic evaluation, an independent cohort of 101 HCC specimens was used for further validation. The basic information about the validated HCC cohort in our center is listed in Table S4 . The protein levels of m⁶A-related signature were determined by MS assay, and the m6AScore was calculated to divide patients into high- and low-risk HCC groups as above mentioned.

In line with the above results in Figure 3 , the high m6AScore group had a higher death risk ( Figures 4A, B ) and advanced clinical stages ( Figure 4C ). In addition, the high-risk group had a remarkably shorter survival time ( Figure 4D ). The 1-, 3-, and 5-year survival rates of high-risk vs. low-risk patients were 87.5% vs. 92.2%, 44.1% vs. 71.8%, and 44.1% vs. 68.3%, respectively. ROC curve analysis further confirmed the accuracy and sensitivity of the m⁶A-related signature, as the AUC reached 0.722 and 0.604 in 1 and 3 years, respectively ( Figure 4E ). The calibration curves for 1 and 3 years were close to the ideal line and showed a satisfactory calibration effect ( Figure 4F ).

Therefore, these validated results confirm the reliability and robustness of our m⁶A-related signature that can readily stratify high-risk HCC patients and evaluate OS. The m⁶A-related signature may facilitate better clinical management of high-risk HCC at an early stage.

To exploit the underlying mechanisms of the m⁶A-related signature, GSEA was performed for GO and KEGG analyses to explore possible biological functions. As shown in Figures 5A, B , in the high-risk HCC group, aberrant DNA replication and recombination, epigenetic dysregulation, cell cycle chaos, and activation of oncogenic signaling pathways such as MAPK, mTOR, and VEGF were obviously enriched. In contrast, disordered amino acid and fatty acid metabolism as well as abnormalities of the cytochrome p450 pathway were markedly enriched in the low-risk HCC group.

Based on GSEA results, we probed potentially effective strategies for high-risk HCC treatments. Some Food and Drug Administration-approved inhibitors and chemotherapeutics that target the cell cycle (CGP.082996), repress DNA replication and synthesis (etoposide and mitomycin C), and inactivate MAPK signaling (PLK1, AKT, and p38 inhibitors) were selected for analysis. Then, half maximal inhibitory concentrations (IC₅₀) for these inhibitors were evaluated in low- and high-risk subgroups in the pRRophetic platform, which is an independent and authorized public resource for predicting inhibitor sensitivity and is based on gene expression microarray data (23). As shown in Figures 5C–K , high-risk HCC patients were more sensitive to these inhibitors, as low drug concentrations effectively inhibited cell proliferation.

These findings suggest that low- and high-risk groups may have different oncogenic mechanisms promoting HCC development. High-risk HCC patients may have elicited better responses to CDK and MAPK inhibitors and the chemotherapeutics etoposide and mitomycin C. Thus, our m⁶A-related signature was beneficial in selecting rational therapeutic strategies.

TMB is an important index to evaluate cancer progression and the efficacy of immunotherapy. We found that the gene mutation landscapes were markedly different between the low- and high-risk groups, as shown in the waterfall chart in Figures 6A, B . Among the low-risk patients, the CTNNB1 mutation was the most common event ( Figure 6A ), whereas in high-risk group, the frequency of the TP53 mutation appeared in more than half of the total mutant events ( Figure 6B ). Then, according to the combination of high/low TMB and m6AScore, HCC patients were divided into four subgroups as indicated in Figure 6C , L-TMB+L-m6AScore, L-TMB+H-m6AScore, H-TMB+L-m6AScore, and H-TMB+H-m6AScore. We found that the H-TMB+H-m6AScore group had the worst prognosis, whereas the L-TMB+L-m6AScore group had the best prognosis.

Recently, immunotherapy combined with immune checkpoint inhibitor approaches has generated encouraging clinical outcomes in different cancers, including HCC (27). TIDE is an algorithm that predicts immunotherapy responses; a high TIDE score frequently indicates poor responses (24). HCC patients with high m6AScores were inclined to have lower TIDE scores regardless of the TMB ( Figure 6D ), suggesting that high-risk groups benefit from immunotherapy. Furthermore, we analyzed the expression levels of several well-documented immune checkpoint molecules, including CTLA4, PD-1, PD-L1, LAG3, and TIM3, in TMB + m6AScore subgroups. Most immune checkpoint molecules had higher expression levels in the high m6AScore subgroup, indicating that high-risk HCC patients may suppress T-cell activation and evade antitumor immunity ( Figure 6E ). However, it also indicated that high-risk HCC patients may be sensitive to immune checkpoint therapy, such as antibodies directed against PD-1, CTLA4, TIM3, and LAG3.

Taken together, m6AScores combined with TMB could be used to stratify the poorest prognostic HCC subgroup. High-risk HCC patients are prone to immune evasion by upregulating immune checkpoint protein expression and may benefit from immune checkpoint therapy.

To systematically delineate the status of the tumor immune microenvironment, we compared the infiltration percentage of various immune cells via ssGSEA analysis. Notably, immunosuppressive cells such as M2 polarized macrophages and T-regulatory cell were obviously accumulated in the high-risk group ( Figure 7A ). Each of the 20 genes in the m⁶A signature individually contributed to immune cell infiltration ( Figure 7B ). It is well known that Foxp3 is an essential transcriptional factor of Treg; CD163, IRF4, and VSIG4 are the markers of M2-macrophages. Thus, the expression levels of these genes were analyzed in the TCGA dataset. The high-risk HCC group exhibited significant upregulation of Foxp3, CD163, IRF4, and VSIG4 ( Figures 7C–F ).

To confirm and validate these findings, multiplex immunofluorescence staining was performed to display tumor microenvironment features in part of our HCC cohort. In low-risk HCC samples, there were a few M2-macrophages (green) and Tregs (purple) observed ( Figure 8A ). Consistently, in high-risk HCC samples, there was obvious infiltration of M2-macrophages and Tregs within the tumor microenvironment ( Figure 8B ).

Hence, these results indicate that m⁶A signature can characterize the tumor immune microenvironment and that the high-risk HCC group may form an immunosuppressive microenvironment with plenty of M2-macrophages and Treg infiltration.