The TCGA data included TCGA-LIHC and other TCGA digestive cancers (TCGA-COAD, TCGA-READ, TCGA-STAD, TCGA-CHOL, and TCGA-PAAD) with clinical information were obtained from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/). ICGC-LIRI-JP datasets were downloaded from the International Cancer Genome Consortium (ICGC) Data Portal (https://dcc.icgc.org/). GSE14520 datasets were obtained from Gene Expression Omnibus databases (GEO). The IMvigor210 cohort (20) data were obtained from the “IMvigor210CoreBiologies” R package. Riaz et al. Cell 2017 (GSE91061) (21) and Lauss et al. Nat Commun 2017 (GSE100797) (22) were gathered from the TIDE website (http://tide.dfci.harvard.edu/) with a detailed clinical information and gene expression data. The transcripts per kilobase million (TPM) values of gene expression were used for further analysis.

Sum up to 18 m⁵C regulators were selected to identify different m⁵C clusters by m⁵C regulators. Based on the gene expression data of these m⁵C regulators in the training cohort (TCGA-LIHC cohort), we used an unsupervised cluster algorithm via the R package “ConsensuClusterPlus”. The R package “NbClust” and Silhouette algorithm were applied to verify the optimal cluster number further. Principal component analysis (PCA) and tSNE were used to visualize different clusters in two-dimensional.

We distinguished the DEGs between two m⁵C clusters via the “DESeq2” R package. The DEGs were defined as adjusted p-value< 0.05 and |log2FoldChange| > 1. A total of 1267 DEGs between two m⁵C clusters were extracted for further analysis.

To explore the functions of DEGs, enrichment analysis of these genes was conducted by R package “clusterprofiler” R package, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases. False discovery rate (FDR)< 0.25 and p-value< 0.05 were considered statistically significant.

To investigate the differential pathway of m⁵C clusters, the gene set variation analysis “GSVA” R package was used to perform GSVA. The gene sets of hallmarks were downloaded from the MSigDB database for GSVA analysis. The adjusted p-value< 0.05 was considered statistically significant.

We performed the prognostic analysis for DEGs in two m⁵C clusters using a univariate Cox regression. As a result, 142 DEGs with a significant prognosis while a p-value ≤ 0.0001 were extracted for further analysis. The prognostic risk model was constructed by using the “glmnet” and “My.stepwise” R package based on LASSO and stepwise Cox method. Finally, we obtained a 6-gene-based prognostic risk model to calculate each patient’s risk score by weighting the Cox regression coefficients. The formula is as follows:

RiskScore = coefficient of each modeled gene * expression of each modeled gene, where “ i “ represents the modeled gene, “ Coef “ represents the coefficients of regression, and “ expr “ represents the expression of gene.We verified the prognostic risk model in validation cohorts. The time-dependent receiver operating characteristic (ROC) curves were used to evaluate the prognostic prediction accuracy of the risk model and the area under the curve (AUC) was measured by the R package “survivalROC”.

To explore the correlation between the prognostic risk model and the immune cell infiltration level, we use the CIBERSORT algorithm to assess the infiltration abundance of 22 immune cells. And we obtained collected 20 inhibitory immune checkpoints from Auslander’s study (23). The Tumor Immune Dysfunction and Exclusion (TIDE) and IPS algorithms were used to evaluate tumor immune escape status. To compare the TME characterization of two different clusters, first, the R package “IOBR” was used to analyze the TME pathway via the single-sample gene set enrichment analysis (ssGSEA) (24, 25) algorithm. Meanwhile, a set of gene signatures that represent a non-inflamed TME, and tumor therapy-associated response were collected from Jiao Hu’s study (26).

Broad Institute-Cancer Cell Line Encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) contained 1019 cancer cell line RNA expression profile data. Drug sensitivity data of cancer cell lines were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC2, https://www.cancerrxgene.org/) and the PRISM Repurposing dataset (19Q4, released December 2019, https://depmap.org/portal/prism/). The GDSC2 contains 809 cell lines and sensitivity data for 198 compounds, and the secondary PRISM contains 499 cell lines and drug susceptibility data for 1448 compounds. All three of these drug datasets provide the area under the dose-response curve (area under the curve — AUC) values as a measure of drug susceptibility, and lower AUC values indicate increased susceptibility to compound responses. The drug susceptibility of the TCGA-LIHC cohort was estimated via the R package “oncoPredict”.

Connectivity Map (CMap) analysis was a complement to further investigate the therapeutic potential of drugs in HCC. We first analyzed differential expression genes between the TCGA-LIHC tumor and normal samples. And 150 up-regulated genes and 150 down-regulated genes were submitted to the CMap website (https://clue.io/query). The CMap analysis yielded a connectivity score for each perturbation, a negative score represents a gene expression pattern of a certain perturbation that is oppositional to the disease-specific expression pattern, suggesting a potential therapeutic effect of this perturbation in this disease.

55 fresh hepatocellular carcinoma tissues and part of paired normal liver tissues were collected from Xiangya hospital and immediately frozen in liquid nitrogen. A retrospective cohort was created, with a total of 55 patients enrolled as of September 1, 2021. After the follow-up period, the clinical characteristics were obtained, and survival analysis and multivariate Cox regression were performed. This study was reviewed and approved by the Xiangya Hospital Medical Ethics Committee of Central South University (No.201806928), and got consent from all participants.

The AG RNAex Pro Reagent (AG21102, Accurate Biology, Changsha, China) was used to extract total RNA from 55 fresh hepatocellular carcinoma tissues and part of paired normal liver tissues according to the manufacturer’s instructions. Before reverse transcription to cDNA, genomic DNA is eliminated by treatment for 2 min at 42°C with gDNA Clean Reagent. Evo M-MLV RT Kit (AG11705, Accurate Biology, Changsha, China) was used to synthesize the complementary RNA. The SYBRR Premix Pro Taq HS qPCR Kit (AG11708, Accurate Biology, Changsha, China) was utilized to achieve real-time quantification. The 2^-△△CT approach was applied to assess the relative expression levels of target genes, which were normalized by GAPDH. The PCR primers are listed in Supplementary Table S6 .

The protein expression levels of Programmed death-ligand 1 (PD-L1; also called B7-H1 or CD274) and cytotoxic T-lymphocyte antigen-4 (CTLA4) were estimated via immunohistochemical (IHC) staining. Briefly, after deparaffinization and rehydration, antigen retrieval was performed by heating the slides in EDTA buffer (G1203; Servicebio, Wuhan, China). Then the slides were treated with 3% hydrogen peroxide (Annjet, Shandong, China) for 25 mins to eliminate endogenous peroxidase and blocked with 3% BSA (G5001; Servicebio) for 30 minutes at room temperature to decrease nonspecific binding. The slides were then incubated overnight at 4°C with rabbit anti-CD274 (PDL1) and anti-CTLA4 primary antibodies at a dilution of 1:200 and 1:300, respectively (CD274: catalog no. Ab205921, Abcam, UK; CTLA4: catalog no. ab237712, Abcam, UK). Next, the slides were incubated with Goat Anti-rabbit IgG/HRP secondary antibody (G1215; Servicebio) for 50 minutes at 37°C. The DAB solution was used for coloration, and hematoxylin was used for counterstaining.

Immunostaining of CD274 (PDL1) and CTLA4 were defined on a scale of 0 to 3 in which 0 means no staining, 1 means mild staining, 2 means medium staining, and 3 means intense staining. The percentage score of stained cells were also calculated on a scale of 1 to 4 in which 1 represents (0–25%), 2 = (26–50%), 3 = (51–75%) and 4 = (76– 100%). To obtain the final score, the intensity score and percentage score were added to reach the final score ranging from 0 to 7.

Correlations between variables were explored with Pearson correlation coefficients. Continuous variables were compared between two groups through the Wilcoxon rank-sum test. Survival analysis including Kaplan-Meier and Cox regression analysis was performed by “survival” R package. The optimal cut-off value in training cohort TCGA-LIHC and validation cohorts were determined by the “surv_cutpoint” function in the “survminer” R package. Unless otherwise stated, P-value< 0.05 was regarded as statistically significant. The whole work was conducted in R 4.1.2 software.

We collected 18 m⁵C regulator genes from previous studies (27–31). The multi-omics landscape of these 18 m⁵C regulator genes was summarized from the TCGA-LIHC cohort ( Figure 2 ). Most of the m⁵C regulator genes were significantly differentially expressed in tumor and normal tissues. For instance, ALKBH1, ALYREF, DNMT1, DNMT3A, DNMT3B, NOP2, NSUN2, NSUN3, NSUN4, NSUN5, TET1, TET2, TET3, TRDMT1, and YBX1 were up-regulated in HCC, while DNMT3L, NSUN6, and NSUN7 were down-regulated ( Figure 2A ). We then explored the comprehensive correlations and prognostic value of m⁵C regulators in the TCGA-LIHC cohort ( Figure 2B ). The univariate Cox regression analysis further demonstrated that most of these m⁵C regulator genes played a risk factor role in HCC ( Figure 2C ). Then we found out that these genes had widespread copy number variations (CNVs) ( Figure S1 ) but infrequent mutations ( Figure 2D ). For example, ALYREF, DNMT3A, DNMT3B, DNMT3L, NOP2, NSUN2, NSUN4, NSUN5, NSUN6, NSUN7, TET1, and TRDMT1 focused on copy number amplification, whereas ALKBH1, DNMT1, TET2, and YBX1 preferred deletion, which suggests that CNVs’ dominant role in the regulation of m⁵C relative to mutation. Metascape (32) is a web tool for gene annotation and analysis, the enrichment analyses of 18 m⁵C regulators are displayed in Figure 2E .

A total of 371 samples from TCGA-LIHC were defined as the training cohort, and further divided into k clusters (k = 2~6) via the “ConsensusClusterPlus” R package. We found out that k = 2 was the best number of clusters according to the CDF curve of the consensus score ( Figures S2A, B ). The Silhouette algorithm and “NbClust” R package further confirmed the result ( Figures S2C, D ). The principal component analysis and tSNE method of 18 m⁵C regulator gene expression showed significant separation between two clusters ( Figures S2E, F ). Finally, 248 patients were classified as Cluster 1, and the rest 123 patients were classified as Cluster 2. Compared to Cluster 1, the Kaplan-Meier analysis showed a worse OS and progression-free survival (PFS) in Cluster 2 ( Figures S2G–H , log-rank test, p< 0.001).

To further explore the differences between the two clusters, we compared the clinical characteristics between the two clusters. The result obtained shows that showed Cluster 2 had a higher proportion of advanced disease stage, tumor stage, and histologic grade ( Figures 3A–C , Fisher’s exact test, p< 0.01). We then analyzed the differentially expressed genes (DEGs) by using the “DESeq2” R package, a total of 1267 DEGs between two m⁵C clusters were obtained according to the adjusted p-value less than 0.05, and the absolute value of log2 fold change less than 1. Next, these DEGs were used for biological functional enrichment analysis. These DEGs were primarily enriched in stromal and catabolic as well as metabolic pathways ( Figures 3D, E ), such as small molecule catabolic process (GO term), response to xenobiotic stimulus (GO term), organic acid catabolic process (GO term), complement and coagulation cascades (KEGG pathway), metabolism of xenobiotics by cytochrome P450 (KEGG pathway), and drug metabolism - cytochrome P450 (KEGG pathway). Further GSVA analysis on hallmark pathway analysis ( Figure 3G ) revealed that patients in Cluster 2 exhibited an obvious enrichment of pathways involved in cell cycle, DNA repair, and carcinogenic activation pathways for example WNT-β-Catenin signaling, TGF-β signaling, PI3K-AKT-MTOR signaling, MYC targets, etc. Moreover, we found out that Cluster 2 had a higher expression of immune checkpoints than Cluster 1 ( Figure 3F ). Our analysis also revealed that Cluster 2 had a significantly increased TIDE score and Exclusion score but decreased Dysfunction score ( Figure 3H ), which suggests that these patients in Cluster 2 may have an immune evasion mechanism.

We constructed a prognostic risk model to quantify the risk for each patient based on our m⁵C clustering. After univariate Cox regression, 142 DEGs with a significant prognosis at a p-value ≤ 0.0001 were analyzed by the LASSO regression algorithm in the TCGA-LIHC cohort ( Figures 4A, B ). In stepwise variable selection procedures, the Cox regression algorithm was used to filter the variables further. Finally, a total of 6 candidate genes including CBX2, SOCS2, LCAT, PBK, KRT17, and FTCD ( Figures 4C, D ) were selected to construct the prognostic risk model. The patients in the training cohort TCGA-LIHC were separated into low- and high-m⁵C score groups by the optimal cut-off value calculated by the “surv_cutpoint” function in the “survminer” R package. In the training cohort, TCGA-LIHC ( Table S1 ), the Kaplan–Meier analysis revealed that the high-m⁵C score group had a significantly worse OS and PFS than the low-m⁵C score group ( Figures 4E, G ). Multivariate Cox analysis of the training cohort revealed that the m⁵C score can be an independent prognostic factor ( Figure S3 ). The 1-, 3-, and 5-year AUC of m⁵C scores for OS was 0.781/0.762/0.711 ( Figure 4F ), and for PFS, the 1-, 3-, and 5-year AUC was 0.724/0.629/0.657 ( Figure 4H ). Likewise, in the validation cohorts ICGC-LIRI-JP and GSE14520, our risk model still had a stable prognostic capability, the high-m⁵C score patients manifested a significantly worse OS than the low-m⁵C score patients ( Figure 4I , Figure S4A ), and the 1-, 3-, and 5-year AUC were 0.717/0.694/0.755 and 0.619/0.622/0.627 ( Figure 4J and Figure S4B ).

To investigate the role of our risk model in the immune microenvironment of HCC, we looked into the correlation among the m⁵C score and expression level of immune checkpoints, immune cells infiltration, TIDE, and IPS scores in the training cohort TCGA-LIHC. The relative expression of most of the 20 immune checkpoints was significantly elevated in the high- m⁵C score group ( Figure 5A ). To further analyze the correlations between the immune microenvironment and our prognostic risk model, we used the CIBERSORT algorithm to calculate 22 types of immune cells, Figure 5B shows the correlation between the m⁵C score and immune infiltration cells. The m⁵C score displayed a significantly positive correlation with macrophages (M0), activated memory CD4 T cells, follicular helper T cells, and Eosinophils (p< 0.01), and a significantly negative correlation with resting memory CD4 T cells, resting NK cells, monocytes, macrophages (M2), and resting mast cells (p< 0.01). Not only CIBERSORT, but we also used TIMER and MCPcounter to analyze immune cell infiltration in high and low-risk groups, the result were displayed in Figure S5 . Additionally, we also analyzed the correlation between the expression of immune checkpoints and our risk model score. Most of the immune checkpoints had a significantly positive correlation with the m⁵C score ( Figure 5C ). We hypothesized that the m⁵C score might be related to immune evasion mechanisms, so we compared the TIDE and IPS scores between the low- and high- m⁵C score groups. The patient in the high-m5C score group had a higher inhibitory immune checkpoints expression level, TIDE and Exclusion score, and a lower Dysfunction score ( Figuress 5D–G , p< 0.001). Myeloid-derived suppressor cells (MDSC) (33) are a group of cells defined by their T cell immunosuppressive functions. We found out that the high-m⁵C score group also presented a higher MDSC score ( Figure 5H , p< 0.001), which further suggests that the high-m⁵C score group may be in an immunosuppressed state. IPS was significantly elevated in the low-m⁵C score group ( Figure S5 , p< 0.001). These findings indirectly demonstrate that our risk model had a valuable relevance in the immune microenvironment of HCC.

To further explore the correlation between TME and our risk model, we also analyzed the TME pathway and several therapeutic signatures with our m⁵C score. The result obtained shows a significant positive correlation between the m⁵C score and most TME and therapeutic pathways ( Figures 6A, B ). Growing evidence has displayed an association between somatic mutations in tumor genomes and response to immunotherapy (34–36). We found that though there was no notable difference in the total mutation frequency between the two m⁵C score groups, the high-m⁵C score group exhibited more TP53 mutation than the low-m⁵C score group (42% versus 15%, Fisher’s exact test, p< 0.001) ( Figures 6C, D ). Recent research has revealed that TMB may act as a biomarker of response to immune checkpoint inhibitors (37–40). We then compared the TMB between the two m⁵C score groups and found out that the high-m⁵C score group had a higher TMB than the low-m⁵C score group (p< 0.05) ( Figure 6E ).

To extend our model for clinical translation, we further explored the possible association of the m⁵C score with potential drug sensitivity. We performed analyses of GDSC2 and PRISM drug sensitivity databases in an attempt to find new potential compounds associated with the m⁵C score ( Figure 7A ). By correlation analysis, we screened the compounds in these two drug sensitivity databases that were negatively correlated with m⁵C scores, and the absolute value of the Pearson’s correlation coefficient obtained > 0.3. Also a significantly different AUC (Wilcoxon test, p< 0.05) was observed between the two groups of high- and low-m⁵C scores ( Tables S2 , S3 ). The top 10 associated compounds with m⁵C scores in two drug databases are presented in Figure 7B . Although these 20 candidate compounds showed increased drug sensitivity in patients with higher m⁵C scores, the above analyses alone can not lead to the conclusion that these chemicals had a therapeutic effect on HCC. Thus, we used the CMap analysis to find compounds in which gene expression patterns were oppositional to the HCC-specific expression patterns. If CMap scores< −1, it represents that these compounds might bring therapeutic benefits to HCC. Secondly, we calculated fold-change differences in the mRNA expression levels of candidates’ drug target genes between tumor and normal tissue, and a higher fold change value indicated a greater potential for a candidate agent for HCC treatment. Finally, a comprehensive literature search was performed in PubMed to find out the experimental and clinical evidence of candidate compounds for treating HCC. All results were presented in Figure 7C ( Table S4 ). Four compounds, including Sepantronium bromide (YM-155), axitinib, vinblastine, and docetaxel, which had strong in vitro and silico evidence, were thought to be the most promising therapeutic candidates in HCC patients with high m⁵C scores. Differences in the sensitivity of these four drugs between the high- and low-m⁵C score groups were displayed in Figure 7D-G . Our results can provide some reference for further clinical translation.

To validate the feasibility of our prognostic risk model in immunotherapy and other GI tumors, we applied our model to three immunotherapy cohorts that included IMvigor210, Riaz, et al. Cell 2017 (GSE91061), and Lauss et al. Nat Commun 2017 (GSE100797), and the other five additional TCGA digestive cancer cohorts that also included TCGA-CHOL, TCGA-PAAD, TCGA-STAD, TCGA-COAD, and TCGA-READ. We then found out that our prognostic risk model can distinguish between the high and low-risk groups with a completely different OS or PFS in three immunotherapy cohorts, and TCGA-CHOL, and TCGA-PAAD cohort ( Figures S6A–E , log-rank test, p< 0.05), but showed no significant difference in TCGA-STAD, TCGA-COAD, and TCGA-READ ( Figures S6F–H ) cohorts.

To increase the probability of early detection of HCC, we developed a tumor diagnostic model by these six genes of our prognostic risk model. 50 HCC tumor samples and 50 paired normal tissues from TGCA-LIHC as the training cohort, 247 HCC tumor samples and 241 normal tissues from GSE14520 as the validation cohort to verify the reliability of the model. The tumor diagnostic model was formulated as follows: logit (P-HCC) = 6.906 - 0.494 * SOCS2 expression level - 1.250 * LCAT expression level – 0.493 * FTCD expression level + 0.602 * KRT17 expression level + 1.950 * PBK expression level + 7.243 * CBX2 expression level. The results obtained are shown in Figures 8A–F . The AUCs of the training cohort and the validation cohort are 0.99 ( Figure 8B ) and 0.960 ( Figure 8E ), the diagnostic model achieved 94% sensitivity and 100% specificity in the training cohort ( Figure 8C ) and 93.9% sensitivity and 90.9% specificity in the validation cohort ( Figure 8F ). Meanwhile, we also constructed 2 diagnostic models of HCC tumors and liver cirrhosis. The training cohort GSE63898 contained 228 HCC samples and 168 liver cirrhosis samples, and the validation cohort GSE25097 contained 268 HCC samples and 40 liver cirrhosis samples. The tumor diagnostic model was formulated as follows: logit (P-HCC) = -2.534 – 1.240 * SOCS2 expression level - 1.320 * LCAT expression level + 0.335 * FTCD expression level – 1.911 * KRT17 expression level + 4.422 * PBK expression level + 2.006 * CBX2 expression level. The results obtained are shown in Figures 8G–L . The AUCs of the training cohort and the validation cohort are 0.973 ( Figure 8H ) and 0.929 ( Figure 8K ), the diagnostic model achieved 91.7% sensitivity and 97.0% specificity in the training cohort ( Figure 8I ) and 86.6% sensitivity and 92.5% specificity in the validation cohort ( Figure 8L ).

In the Xiangya HCC cohort ( Table S5 ), the prognoses were approximately the same as those seen in the training set. Survival analysis verification testified that the prognostic risk score model was a great independent prognostic factor in HCC ( Figures 9A–C ). The expression of 6 candidate genes in HCC and paired normal liver tissues, CBX2, PBK, and KRT17 were higher in tumors than in normal tissues. On the contrary, FTCD, LCAT, and SOCS2 were lower in tumors than in normal tissues ( Figure S7 , p< 0.05).

To evaluate the predictive value of the prognostic risk score model on immunotherapy response, we selected CD274 and CTLA4, two immune checkpoints approved for clinical treatment, for further verification of the model. We collected pathological tissues from the 55 patients in the Xiangya HCC cohort. The tissues were used to make paraffin sections and for immunohistochemical staining. The results obtained show that higher CD274 and CTLA4 expression in the high prognostic risk score group compared to that of the low prognostic risk score group, a finding in line with the results of our bioinformatics analysis ( Figures 9D–G ). Our research provides further verifications that patients with a low prognostic risk score could benefit from immunotherapy to a higher degree than those with a high prognostic risk score. Hence, a prognostic risk score can be utilized as a potential biomarker for immunotherapy response.