DNA methylation regulator-based molecular subtyping and tumor microenvironment characterization in hepatocellular carcinoma

Backgroud: Although recent studies have reported the regulation of the immune response in hepatocellular carcinoma (HCC) through DNA methylation, the comprehensive impact methylation modifications on tumor microenvironment characteristics and immunotherapy efficacy has not been fully elucidated.

Methods: In this research, we conducted a comprehensive assessment of the patterns of DNA methylation regulators and the profiles of the tumor microenvironment (TME) in HCC, focusing on 21 specific DNA methylation regulators. We subsequently developed a unique scoring system, a DNA methylation score (DMscore), to assess the individual DNA methylation modifications among the three distinct methylation patterns for differentially expressed genes (DEGs).

Results: Three distinct methylation modification patterns were identified with distinct TME infiltration characteristics. We demonstrated that the DMscore could predict patient subtype, TME infiltration, and patient prognosis. A low DMscore, characterized by an elevated tumor mutation burden (TMB), hepatitis B virus (HBV)/hepatitis C virus (HCV) infection, and immune activation, indicates an inflamed tumor microenvironment phenotype with a 5-year survival rate of 7.8%. Moreover, a low DMscore appeared to increase the efficacy of immunotherapy in the anti-CTLA-4/PD-1/PD-L1 cohort.

Conclusions: In brief, this research has enhanced our understanding of the correlation between modifications in DNA methylation patterns and the profile of the tumor microenvironment in individuals diagnosed with HCC. The DMscore may serve as an alternative biomarker for survival and efficacy of immunotherapy in patients with HCC.

Background

DNA methylation is known to involve in multiple biological processes, including cancer progression (1), through the modification of chromosomal proteins, which alter the 3-dimensional conformation of the genome and/or protein-DNA interactions (2). DNA methylation is a type of dynamic reversible process in mammalian cells and is regulated by transmethylases, demethylases and recognizing proteins, which are defined as “writers”, “erasers” and “readers”, respectively (3). A recent study suggested that DNA methylation is a dynamic marker correlated with tumor immune escape and T-cell exhaustion (4).

Hepatocellular carcinoma (HCC) has a high incidence of malignant tumors and is the third leading cause of death among all types of cancers throughout the world (5). Approximately 80% of HCC cases are linked to persistent infections caused by either hepatitis B virus (HBV) or hepatitis C virus (HCV) (6). In addition, viral infection represents an important contributor to the epigenetic and tumor immune changes observed in HCC (7). Immune checkpoint inhibitors (ICIs) have been demonstrated to improve survival but only in a limited fraction of multiple tumor types, including HCC (8–10). There is growing evidence implying a link between DNA methylation and tumor immunity/immunotherapy. The innate immune response, which serves as a tumor suppressor, can be suppressed by the DNA methyltransferase 1 (DNMT1), which maintains the silencing of retrotransposable elements (11, 12). In addition to PD-1/PD-L1 expression and tumor mutational burden (TMB), several methylation regulators, such as ubiquitin-like with plant homeodomain (PHD) and ring finger domains 1 (UHRF1) and Ten-eleven translocation 1 (TET1), have been reported to be potential biomarkers for ICIs therapy (13, 14). Xu et al. found that the IFN-γ/JAK/STAT/TET2 signaling pathway is involved in influencing tumor immunity, and that stimulating TET2 activity increases the efficacy of anti-PD-L1 drugs in solid tumors (15).

However, these previous studies have been limited by technical constraints and have focused only on a small number of methylation regulators and cell types. Nevertheless, because the antitumor effect of these agents involves multiple tumor suppressors working together in a well-coordinated manner. It is essential to gain a comprehensive understanding of how different DNA methylation regulators impact the tumor microenvironment (TME) to optimize the efficacy of immunotherapy for HCC. In this research, we utilized an unsupervised clustering technique to analyze gene expression data from 21 DNA methylation regulators in various cohorts to detect distinct patterns of modifications in DNA methylation regulators. We then evaluated the TME profiles represented by three distinct methylation regulator patterns. Finally, a DNA methylation score (DMscore) model was established to quantify the methylation status of the individuals. We validated the potential predictive value of the DMscore as an alternative biomarker for survival and immunotherapy efficacy.

Materials and methods

Data mining

The workflow of this study is shown in Supplementary Figure 1A. We analyzed 33 histologically confirmed HCC tissues and matched noncancerous liver tissue as normal controls. We named this cohort the “discovery cohort”, after which the preparation of mRNA libraries was performed followed by sequencing, and the gene expression profiles were analyzed as previously described (16, 17). The study received approval from the Research Ethics Committee at the First Affiliated Hospital of Zhejiang University, and written informed consent was acquired from all participating patients. Among the 33 patiens with HCC, 1 case was alcohol-induced cirrhosis and progressed to HCC, and the remaining 32 cases were HBV-related HCC. There were 27 males, compared with only six female patients. The lowest positive rate of eight serum markers was alkaline phosphatase (ALP), which was higher than the normal range in only two patients, and the positive rate was only 6% (2/33). The positive rate of alpha fetal protein (AFP) was the highest, which was 60% (20/33). All patients had grade A Child-Pugh scores. The clinical data of these patients were summarized in Supplementary Table 7.

Other gene expression data and full clinical annotations were downloaded from The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC) and Gene Expression Omnibus (GEO) databases. For the TCGA-LIHC (liver hepatocellular carcinoma) cohort, RNA sequencing (FPKM) of gene expression data were acquired using the R package TCGAbiolinks (18) and subsequently transformed into the more comparable transcripts per kilobase million (TPM) format. Somatic mutation and methylation 450K data were downloaded from the UCSE Xena browser. HBV/HCV infection was determined to positive if the participant of interest met the following criteria: HBV surface antigen; HBV DNA; HBV Core antibody; Hepatitis C Antibody; Hepatitis C Viral RNA; HCV Genotype (19). According to these criteria, 164 patients in the TCGA-LIHC cohort were HBV/HCV positive. For the HCC datasets in ICGC, RNA-seq data (raw counts), somatic mutation data and clinical information were obtained from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). These specimens were predominantly obtained from a cohort of individuals in Japan who had been diagnosed with either HBV or HCV infection (20). The raw count format of the RNA-seq data was also transformed into TPM values. Batch effects between the TCGA-LIHC cohort and the LIRI-JP RNA-seq dataset were corrected using the “ComBat” function of the sva package (Supplementary Figures 1B, C). All the baseline information of the patients in the eligible HCC datasets is summarized in Supplementary Table 1.

Tumor mutation burden analysis

To calculate the TMB of each patient, the total number of nonsynonymous mutations counted was divided by the exome size (38 Mb was utilized as the exome size). To calculate the TMB of each sample in the TCGA-LIHC cohort, we selected somatic mutation data processed by the VarScan platform and visualized them using the R package “maftools”.

The details of the mutations in the LIRI-JP cohort were visualized via a waterfall plot generated with the “GenVisR” package in R Studio software.

GSVA analysis and functional annotation

To assess the level of activity in each biological pathway, the R package GSVA was used to compute the ssGSEA score for every sample (21). The gene set used to identify different immune cell types within the TME was obtained from Charoentong’s study and included 23 distinct subtypes of human immune cells such as activated CD8+ T cells, dendritic cells in an activated state, macrophages, natural killer T cells, and regulatory T cells, and so on (Supplementary Table 2) (22). The relative abundance of each immune cell type in TME was determined by an enrichment score acquired from ssGSEA analysis.

To investigate the disparities in biological mechanisms associated with DNA methylation, we conducted GSVA enrichment analysis employing the GSVA package. The “c2.cp.kegg.v7.4.symbols” gene sets were downloaded from the MSigDB database for GSVA.

Unsupervised clustering of 21 DNA methylation regulators

Unsupervised clustering methods were used to identify different DNA methylation modification patterns and classify patients for further analysis. A set of 21 regulators (Supplementary Table 3) obtained from published literature was extracted from either the discovery cohort or the data-mining cohort to detect distinct DNA modification patterns facilitated by DNA methylation modifiers. We employed the ConsensuClusterPlus package to execute the aforementioned procedures, ensuring classification stability through 1000 repetitions (23).

Identification of differentially expressed genes and enrichment analysis

We utilized the R package limma to detect gene expression related to DNA methylation modifications and identified pathways that were enriched with relevant associations. We also conducted KEGG analysis using the R package clusterProfiler, employing a significance threshold of p < 0.05 and FDR < 0.05.

Generation of the DNA methylation score

Univariate Cox model analysis was used to assess the association between DEGs and overall survival. A total of 468 DEGs associated with a significant prognosis (p-value < 0.05) were identified for further analysis. Principal component analysis (PCA) was subsequently preformed to generate a signature relevant to DNA methylation. Both PC1 and PC2 were selected as signature scores. This approach offered the benefit of emphasizing the score on the subset containing a significant cluster of highly correlated (or anticorrelated) genes while reducing the impact of genes that do not exhibit similar patterns as other members within the subset. We defined the DMscore using a method similar to that reported previously (24–26):

$$\text{DMscore }={\sum }^{}(PC1i+PC2i)$$

where ‘i’ represents the expression level of genes associated with methylation regulation and their impact on prognosis.

Collection of genomic and clinical information on immune-checkpoint blockade

The present investigation employed five pretreatment tumor expression profiles obtained from cohorts receiving immune checkpoint blockade therapy to assess the agreement of the immunotherapy response. Data on anti-PD-L1 therapy efficacy in patients with metastatic urothelial tumors were obtained by the R package IMvigor210CoreBiologies (27). The sequencing raw count data were normalized and transformed into TPM values. We investigated the immunotherapeutic features of TCGA-LIHC patients by utilizing the publicly accessible Cancer Immunome Database (TCIA), which comprises relevant clinical pathology data.

Statistical analysis

Statistical analyses were conducted using R software (version 4.0.3). The Wilcoxon test was used to compare differences between two groups. One-way ANOVA and Kruskal-Wallis tests of variance were used as parametric and nonparametric methods, respectively. Principal component analysis (PCA) was also conducted to investigate the distributions of the different groups using “prcomp”. Correlation coefficients between the TME-infiltrating immune cells and the expression of regulators were computed by Spearman and distance correlation analyses. We conducted both univariate and multivariate Cox regression analyses to identify potential risk factors. Additionally, we assessed the statistical significance of the differences insurvival rates across various risk groups using Kaplan-Meier (K-M) analysis with log-rank p-values. Venn diagrams, heatmapss, boxplots, forest plots, and alluvial diagrams were drawn using R. The immune cell profiles in TME of each sample in the data-mining cohort were analyzed using CIBERSORT (28). Next, MethylMix was utilized to analyze DNA methylation data and paired gene expression data to detect noteworthy DNA methylation occurrences that influence the expression of corresponding genes, thereby revealing their classification as genes driven by DNA methylation (29).

Results

Expression of DNA methylation regulators in the discovery cohort

As shown in Figure 1A, 21 DNA methylation regulators were included in this study. As shown in Figure 1B, the protein-protein interaction (PPI) network revealed extensive protein interactions among methylation modifiers, with the exception of QSER1, which collaborates with TET1 to inhibit de novo DNA methylation mediated by the enzyme DNMT3 (30). Next, we performed a comparative analysis of the expression levels of DNA methylation regulators in paracancerous and HCC tissues within the discovery cohort (Figure 1C). Notably, HCC tissues presented significantly increased expression of QSER1, MBD3, UNG, DNMT3B, DNMT3A, SMUG1, MBD4, DNMT1, MBD2, TET1, TDG, MBD1, UHRF2, MECP2, ZBTB38, ZBTB33, TET3 and UHRF1. Only NTHL1 was expressed at significantly lower levels in HCC tumor tissue. Based on the expression levels of 21 DNA methylation regulators, we could distinguish HCC samples from normal liver samples (Figure 1D).

Figure 1

Figure 1 Landscape of DNA methylation regulators in HCC. (A) Summary of 21 DNA methylation regulators and their potential biological functions. (B) The PPI network downloaded from the STRING database indicates the interactions among the regulators. (C) The expression of 21 m6A regulators in normal and HCC tissues in the discovery cohort (*p-value < 0.05, **p-value < 0.01, ***p-value < 0.001). (D) Principal component analysis (PCA) of the expression profiles of 21 regulators to distinguish tumors from normal samples in the discovery cohort. (E) Spearman’s correlation heatmap between 21 DNA methylation regulators and immune cells in the discovery cohort (*p-value < 0.05, **p-value < 0.01, ***p-value < 0.001).

The ssGSEA algorithm was used to explore the association between the expression of these regulators and tumor immune cell infiltration in the TME. The expression of MECP2, TET2, ZBTB4, and DNMT1 was positively correlated with more than half of the immune cells (Figure 1E). The immune cell profiles most strongly associated with these marks are immature B cells, MDSCs, Macrophages (Kupffer cells), Mast cells, T regulatory cells, Th17, and T follicular cells. These cell types are largely immunosuppressive and involved in immune evasion, and thus often associated with a poor prognosis in patients with HCC (31). Poorly correlations were found between most methylation regulators and key effectors of anticancer immunity include CD8+ T cells, Eosinophils, Nature killer cells (NK), and Dendritic cells (32–34).

Landscape of DNA methylation regulators in data-mining cohort

By examining the expression patterns of DNA methylation regulators and their correlation with TME cell infiltration in the discovery cohort, we proceeded to investigate the genetic modifications of these regulators in a data-mining cohort. Within the HCC genome of the TCGA-LIHC cohort, the overall mutation frequency among all the regulators was determined (Figure 2A). Among these regulators, the eraser TET1 had the highest mutation rate (2%), and its mutation rate is a favorable prognostic marker of immunotherapy (14). Among the HCC patients, no mutations were observed in twelve other regulators, namely DNMT3B, MBD3, MBD4, ZBTB4, UHRF1, UHRF2, UNG, NTHL1, SMUG1, MBD2, ZBTB33 and QSER1 (Figure 2A). Analysis of copy number variation (CNV) alteration frequency revealed that regulators exhibited amplification in terms of copy number. Notably, MDB1/2/3, TET2, UHRF1, DNMT1, and ZBTB4 exhibited widespread CNV deletion (Figure 2B). The chromosomal location of CNV among these regulators are shown in Figure 2C.

Figure 2

Figure 2 Multiomics landscape of DNA methylation regulators in the TCGA-LIHC cohort. (A) The mutation frequency of regulators in the TCGA-LIHC cohort. Each column represents the TMB. The number on the right shows the mutation frequency of each variant type. The lower bar represents the sample annotations. (B) The CNV frequency of the regulators in the TCGA-LIHC cohort. (C) The chromosomal locations of the regulators. (D) Differences in the gene expression levels of the regulators between normal and tumor patients in the TCGA-LIHC cohort. *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

In the LIRI-JP cohort, the mutation rate of these regulators was similar to that in the TCGA-LIHC cohort (Supplementary Figure 2A). The reader ZBTB38, known as a tumor suppressor, had the most mutations, followed by UHRF2. However, no mutation in the eraser TET1 was found in these patients. The chromosomal locations of these mutations are shown in Supplementary Figure 2B.

To ascertain whether the expression patterns of methylation regulators found in the discovery cohort also occurred in the data-mining cohort, we investigated the alterations in expression in the data-mining cohort. The results showed that most of these regulators were highly expressed in HCC tissues, while MBD4 and TET2 were expressed at significantly lower levels in HCC tissues (Figure 2D, Supplementary Figure 2C). As a prognostic risk factor and most mutated gene, TET1 strongly inversely correlated with cytotoxic lymphocyte infiltration, which include CD8+ T cells and CD56dim NK cells, and positively correlated with Th2, which mostly related to tumor-promoting actives (35). In summary, these analyses revealed that variations in the genetics and expression of these DNA methylation regulators might play pivotal roles in the initiation, progression and heterogeneity of HCC.

Methylation modification patterns mediated by 21 regulators

Gene expression and clinical information from 608 patients with HCC from the TCGA-LIHC and LIRI-JP cohorts were collected for analysis. Univariate Cox regression analysis revealed associations between the expression of 21 DNA methylation regulators and overall survival in patients with HCC (Supplementary Figure 3A). The results revealed that TET1 and QSER1 had the highest hazard ratios (HRs) of 1.431 and 1.209, respectively, while the absence of a regulator was found to be a significant favorable factor (p-value < 0.05). We created a network of methylation regulators to provide a comprehensive overview of the interactions and connections between DNA methylation and prognosis in patients with HCC (Figure 3A, Supplementary Table 4). We found that the expression of DNA methylation regulators was correlated not only with in the same functional category, but also among writers, erasers, and readers.

Figure 3

Figure 3 DNA methylation modification pattern and relevant biological characterization. (A) Correlations and prognosis of DNA methylation regulators in HCC patients. The red line represents a positive correlation with a p-value < 0.0001, and the blue line represents a negative correlation with a p-value < 0.0001. The size of the node represents the p-value value of the log-rank test. The right semicircle represents the prognostic factors for HCC: green represents favorable factors for OS, and purple represents risk factors for OS. (B) Consensus clustering of 21 regulators in 608 data-mining cohort samples. (C) Survival analysis of patients in the data-mining cohort according to distinct DNA methylation patterns. (D, E) GSVA enrichment analysis showing the activation states of biological pathways associated with distinct methylation modification patterns. A heatmap was generated to visualize these biological processes, in which red represents activated pathways and blue represented inhibited pathways. The data-mining cohort was used for sample annotation. (D) Cluster A vs. Cluster B; (E) Cluster B vs. Cluster C.

To explore additional insights into the potential patterns of DNA methylation regulation, an analysis was conducted using the expression profiles of the data-mining cohort. Three distinct expression patterns were divided using unsupervised clustering (Supplementary Figures 3B-E) and were termed Clusters A, B, and C, respectively (Figure 3B). Among these, 130 patients had pattern A, 320 had pattern B, and 158 had pattern C. An analysis of prognostic factors indicated a strong correlation between pattern A of DNA methylation regulators and intermediate survival duration (Figure 3C).

Immune cell infiltration characteristics associated with distinct methylation modification patterns

We explored the differences in signaling pathways among the three patterns by GSVA enrichment analysis. As shown in Figure 3D and Supplementary Table 5, pattern A was enriched in nucleic acid biological processes, including RNA polymerase activity, ate biosynthesis, and pyrimidine metabolism. Pattern B was enriched for carcinogenic and immune fully activation, included the TGF-beta signaling pathway, focal adhesion and adheres junction, JAK-STAT signaling pathway, T cell receptor signaling pathway and cytokine-cytokine receptor interaction (Figures 3D, E). However, pattern C was strongly related to DNA replication and repair mechanisms, which included nucleotide excision repair, base excision repair, DNA replication, and mismatch repair (Figure 3E). The results of this research indicate that disruption of DNA methylation regulators can influence overall survival by affecting methylation levels, thus contributing to the significant variability observed in HCC.

We utilized the ssGSEA algorithm to perform an extensive evaluation of the infiltration of TME cells. Consistent with favorable survival outcomes, cluster A exhibited the lowest level of activated CD4+ T cell infiltration and the highest level of activated CD8+ T cell infiltration (Figure 4A; Supplementary Table 5), could be classified as immune-inflamed phenotype. Surprisingly, TME cell infiltration analysis indicated that cluster B contained the most infiltrating immune cells, including dendritic cell, natural killer cell, eosinophils, B cell, MDSC, macrophage, mast cell, monocyte, neutrophil, Tfh (follicular helper T) cell, and so on (Figure 4A, Supplementary Table 5). However, patients in Pattern B did not show a matching survival advantage (Figure 3C). Previous research has indicated that tumors with an immune-excluded phenotype exhibit a significant abundance of immune cells in the stromal compartment encircling tumor cells (36). Therefore, cluster B might be identified as immune-excluded phenotype. Additionally, cluster C was similar immune-desert phenotype, distinguished by immunosuppression (Figures 3D, E, Figure 4A).

Figure 4

Figure 4 TME characteristics and transcriptome traits associated with DNA methylation modification patterns. (A) Boxplot of the relative immune cell abundances for the three DNA methylation patterns. (*p-value < 0.05, **p-value < 0.01, ***p-value < 0.001). (B) Spearman’s correlation heatmap between 21 DNA methylation regulators and immune cells in the data-mining cohort. Red indicates a positive correlation; blue indicates a negative correlation. (*p-value < 0.05, **p-value < 0.01, ***p-value < 0.001). (C) Functional annotation of DEGs among methylation clusters using KEGG enrichment analysis. The color depth of the bar plots represents the q-value of the enrichment. (D) The proportions of patients with the three modification patterns according to HBV/HCV viral infection status. (E) The difference in the number of infiltrating cells in each TME according to HBV/HCV viral infection status. The asterisks represent the statistical p-value: *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

As DNA methylation regulators correlated with immune cell infiltration in the TME in the discovery cohort (Figure 1E), we further examined the correlation in the data-mining cohort. The expression of two writers, DNMT3A and DNMT3B, and three erasers, TET1/3 and QSER1, was negatively correlated with most immune cells (Figure 4B). However, the expression of ZBTB4, MBD2 and DNMT1 was positively correlated with that of most immune cell types (Figure 4B), which was consistent with the findings in the discovery cohort.

A total of 1105 DEGs among the three DNA methylation patterns were identified, and KEGG pathway analysis revealed that those DEGs were enriched in HBV infection events (Figure 4C). The distinct etiological mechanisms of HCC, both viral and non-viral, exhibit a significant association with a unique TME (10). The level of immunosuppression observed in microenvironments associated with viruses was found to be considerably greater than that observed in unrelated microenvironments (37). Uninfected HCC exhibited a relatively reduced response to ICIs in comparison to infected patient with HCC (38). Further investigation revealed that patients with viral-infected HCC predominantly exhibited Pattern C methylation modifications, whereas those without viral infection were characterized by Pattern B modification patterns (Figure 4D). The expression of DNMT1, DNMT3A, DNMT3B, MBD1, QSER1, TDG, TET1, TET2, TET3, UHRF1, UHRF2, UNG and ZBTB33 was significantly high in the virus-positive group (Supplementary Figure 4A).

The above results suggest that viral infection has the potential to impact the function of regulators responsible for DNA methylation, and might result in alterations in patterns of DNA methylation modifications. Additionally, we found the activated B cell, activated CD8+ T cell, eosinophil, immature B cell, neutrophil and Th17 were significantly down-infiltrated in viral-positive patients compared with viral-negative patients, while activated dendritic cell, gamma delta T cell and natural killer T cell was remarkably upregulated (Figure 4E). These results hold promise for advancing our understanding of the mechanisms underlying variations in methylation modification patterns observed in tumors.

After the completion of the aforementioned analyses, notable discrepancies were identified in the infiltration characteristics of TME cells with regard to patterns of DNA methylation modification. We subsequently employed the CIBERSORT technique to compare variations in immune cell components across the three patterns. The CIBERSORT analysis also showed significant correlation among DNA methylation regulators and immune infiltration (Supplementary Figure 4B). Conversely, macrophages M2 exhibited a significant enrichment of cluster A, while activated CD4+ memory T cells showed a notable enrichment of cluster B (Supplementary Figure 4C). Besides the stronger correlation to CD8+ T cell, and the differential expression of QSER1 in virus infection status, we found the low expression of QSER1 significantly prolonged OS (Supplementary Figure 4D). Taken together, these results indicated that DNA methylation modification did changes TME infiltration.

Generation of DNA methylation gene signatures and functional annotation

As mentioned above, we identified 1105 DEGs associated with DNA methylation phenotypes using the limma package (Supplementary Figure 5A). The DEGs were related to the cell cycle, HBV infection, DNA replication and transcriptional events according to KEGG analysis (Figure 4C), revealing that the dissimilar clinical and transcriptomic features observed in the DNA methylation regulatory patterns could arise from variations in the distinct genes associated with the DNA methylation signatures. We identified 468 prognostic genes (p-value < 0.05) via univariate Cox model analysis (Supplementary Table 6). According to the expression levels of 468 genes, an unsupervised clustering analysis was conducted to categorize patients with HCC into three distinct clusters. These clusters, namely, DNA methylation geneCluster A/B/C, exhibit varying prognosis in terms of survival and immune-infiltrating features. (Supplementary Figures 5B-J). These findings provide additional support for the proposition that each HCC subtype exhibits distinct clinical and immune characteristics.

However, these studies were limited by their population size and lacked the ability to accurately predict the methylation modification patterns in individuals. To assess DNA methylation status at the individual level, we developed a risk score system called the DMscore, which is based on 468 signature genes associated with DNA methylation. Kruskal-Wallis test revealed significant discrepancy on DMscore between methylation regulator patterns. From this point of view, the median DMscore for DNA methylation regulator pattern C was found to be the highest, whereas methylation modification Pattern A exhibited the lowest median DMscore (Figures 5A, C). Gene cluster A showed the highest median DMscore (Figures 5B, C). An alluvial diagram was generated to illustrate alterations in attributes among individual patients (Figure 5C). Furthermore, we also tested the correlation between the TME ssGSEA score and the DMscore. As shown in Figure 5D, the DMscore was significantly correlated with the infiltration of most of the 23 immune cell types.

Figure 5

Figure 5 Construction of the DMscore. (A) Differences in the DMscore among the three methylation clusters. (B) Differences in DMscore among three geneClusters. (C) Alluvial diagram showing the changes in methylation clusters, geneCluster, DMscore, HBV/HCV viral infection status, disease stage and survival status. (D) Correlations between the DMscore and TME infiltration in the data-mining cohort. The asterisks represent the statistical p-value < 0.05 (E) Survival analysis of the high- and low-DMscore groups in the data-mining cohort. Log-rank test, p-value < 0.001. (F) Kaplan-Meier curves for patients in the data-mining cohort stratified by both TMB and DMscore. Log-rank test, p-value < 0.001. (G) Stacked bar plot depicting different proportions of patients with HBV/HCV infection in the high- and low-DMscore groups in the data-mining cohort. (H) Boxplot of the DMscores for distinct HBV/HCV infection status groups in the data-mining cohort, p-value = 2.7e-06. The asterisks in A, B, and H represent the statistical p-value: *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

Subsequently, we assessed the prognostic value of the DMscore in patients with HCC. With the cutoff value of 10.87888 identified by the survminer package, HCC patients were separated into high and low DMscore groups. All the gene Cluster B samples belongs to the low DMscore group (Figure 5C). Patients in the low DMscore subgroup exhibited an extended duration of survival, with a 5-year survival rate that was twice as high as that in the high DMscore subgroup (7.8% compared to 3.0%) (Figure 5E). When TMB and DMscore were combined, the patients in the group with a low DMscore and high TMB exhibited superior overall survival compared to that of the remaining groups (Figure 5F, Supplementary Figure 8). These results suggested that the DMscore is a prognostic biomarker that could effectively predict survival of patient with HCC. In addition, we found that the DMscore was significantly correlated with viral infection status (Figures 5G, H).

Additionally, we analyzed TCGA-LIHC methylation data. Unsupervised clustering based on the methylation level of methylation-driven genes classified patients into two groups (Supplementary Figure 6A), which we called MethCluster C1 and C2. Similarly, the infiltration of most of the TME immune cells was significantly different between MethCluster C1 and C2 (Supplementary Figure 6B). In addition, MethCluster C1 had a significantly greater DMscore (Supplementary Figure 6C).

The present study demonstrated that the DMscore could reflect genomic methylation modifications, and effectively predict the survival of patients with HCC.

Predictive value of the DMscore in immunotherapy

The individuals who will experience the most significant advantages from immunotherapy have been identified, as immnotherapy has improved clinical outcomes in the treatment of diverse tumor types. Considering the correlation between DMscore and immune infiltration, we further explored whether the DMscore could be a prognostic factor for immunotherapy efficacy by analyzing four ICI treatment cohorts.

In the TCIA-LIHC cohort (anti-CTLA-4/PD-1 therapy), patients in the DMscore-Low group exhibited a significantly higher response rate (Figures 6A–D), indicating the immunotherapeutic benefits of CTLA-4/PD-1 antibody therapy in patients with a low DMscore.

Figure 6

Figure 6 The effect of the DMscore in the immune checkpoint treatment cohorts. (A-D) Immunotherapeutic benefits of anti-CTLA-4/PD-1 therapy in the TCGA-LIHC cohort. (E) Survival analyses for patients with low- (32 patients) or high- (202 patients) DMscore. (E) Patient groups in the anti-PD-L1 immunotherapy cohort were evaluated via Kaplan-Meier curves (IMvigor210 cohort). Log-rank test, p-value < 0.001. (F) The proportion of patients who responded to PD-L1 antibody immunotherapy in the low- or high-DMscore groups. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. (G) Distribution of the DMscore in distinct anti-PD-L1 clinical response groups. Wilcoxon test, p-value = 0.00067. (H) Correlations between the DMscore, TMB, and anti-PD-L1 immunotherapy response in the IMvigor210 cohort. (I) Survival analyses of patients receiving anti-PD-L1 immunotherapy stratified by both the DMscore and TMB using Kaplan-Meier curves. Log-rank test, p-value = 0.003. *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

Consistent with the above findings, a similar result was also obtained in the IMvigor210 cohort (anti-PD-L1 therapy). According to the K-M survival analysis of the anti-PD-L1 cohort (IMvigor210), patients with a low DMscore exhibited significantly improved overall survival (p-value < 0.001) (Figure 6E) and significant immunotherapeutic benefit (Figures 6E-G). Additional investigations into the measurement of TMB, a biomarker for immunotherapy efficacy, demonstrated a significant association between a decreased DMscore and increased TMB (Supplementary Figure 7A). A significant negative correlation was found between the DMscore and TMB (Figure 6H). Survival analysis revealed that patients with a combination of a high DMscore and high TMB had a great survival advantage (Figure 6I).

The detection of immune phenotypes in the IMvigor210 cohort was feasible because we could explore the variation in the DMscore across different phenotypes. Our results suggest that there is a notable association between decreased DMscore and immune phenotypes prevalent in desert regions, which may hinder the effectiveness of checkpoint inhibitors in combating tumors within this specific phenotype (Supplementary Figure 7B).

In the studies suggested a notable association between DNA methylation alterations and tumor immune characteristics as well as the efficacy of ICI immunotherapy. The utilization of the DMscore could aid in predicting the response to ICI immunotherapy and potentially enhance its predictive accuracy when used in combination with the TMB.

Discussion

Ongoing research has consistently confirmed the significant role of abnormal DNA methylation in promoting genome instability and its impact on regulating antitumor immune responses and immunotherapy outcomes (14, 39, 40). However, there is still a need for a comprehensive understanding of the overall modulation of DNA methylation modification and the immune microenvironment among patients with HCC.

The initial phase of the research involved identifying 21 enzymes responsible for regulating DNA methylation that exhibited differential expression in both HCC carcinomas and paracancerous tissues. Subsequently, we validated the changes in methylation regulatory enzyme levels during HCC development by analyzing the expression of these enzymes in two publicly available datasets (TCGA-LIHC and LIRI-JP). In contrast to Sun et al.’s findings of TET1 underexpression in HCC (41), our study consistently demonstrated overexpression of TET1 in tumor tissues across all three datasets. Although the differential expression profiles of some genes were not entirely consistent across the three datasets and some genes showed differential expression in opposite directions, the 21 DNA methylation regulatory enzymes overwhelmingly showed significant differences in expression levels between cancer and paracancerous tissues in all three datasets. The variation in the orientation of disparities in expression could be attributed to the considerable diversity observed in HCC (42). The distinct expression patterns of these regulatory enzymes in HCC and adjacent noncancerous tissues suggest the significant involvement of DNA methylation regulatory enzymes in the progression of HCC.

In this research, we classified three unique patterns of DNA methylation modification by analyzing the expression levels of 21 regulators involved in DNA 5-mc methylation. Cluster A exhibit activated adaptive immunity, consistent with the immune-inflamed phenotype; Cluster B featured activation of innate immunity and stroma, corresponding to the immune-excluded phenotype; and Cluster C had immunosuppressive characteristics, similar to the immune-desert phenotype. The immune-excluded and immune-desert types could be classified as noninflamed tumor, namely, “cold tumors”. A “hot tumor”, characterized by significant infiltration of immune cells in the TME, exhibits an immune-inflamed phenotype (36, 43, 44). Consistently, patients with pattern A had the longest overall survival. Although pattern B exhibited significantly abundant immune cells, the activated TGF-β pathway resulted in reduced infiltration of T cells into tumors and weakened their ability to eliminate tumor cells. This leads to an immune-excluded phenotype and a less favorable prognosis (27, 45). The presence of immune-desert phenotypes featured with immune tolerance and ignorance, due to a lack of activated and primed T cells. Therefore, it is not surprising that patients in Cluster C had poorer survival outcomes than patients in other clusters (46).

In addition, it has been demonstrated that the variations in the mRNA transcriptome among the three DNA modification patterns are significantly linked to biological pathways associated with DNA replication and viral immunity. Similarly, similar to the clustering outcomes of the phenotypes resulting from DNA methylation modifications, we identified three genomic subtypes based on the expression levels of the DEGs in the mRNA. These subtypes also exhibit distinct survival prognoses and are significantly correlated with immune status. This once again highlighted the crucial role of 5-mc methylation modification in influencing diverse immune environments.

Considering the heterogeneity of methylation modification between individuals, we established a DNA methylation scoring (DMscore) system to calculate the methylation modification features of individuals with HCC. As expected, Cluster A was associated with geneCluster A/B and was characterized by immune activation and a lower DMscore, which corresponded to longer survival time (Figures 3C, 4A, 5A-E). The objective response rates (ORRs) observed in the CheckMate 040 trial indicated that patients treated with nivolumab, a PD-1 inhibitor, exhibited lower efficacy in patients infected with HBV/HCV than in those without viral infections (47). Consistent with these findings, we found that patients with HBV/HCV infection had relatively high DMscores (Figures 5G-H), while patients who benefited from anti-CTLA-4/PD-1 therapy tended to have relatively low DMscores (Figures 6A-G). The presence of a virus-induced immune-suppressing tumor microenvironment may explain the limited effectiveness of immunotherapy in treating HCC. These findings also reinforce the significance of using the DMscore as a predictive tool in guiding immunotherapeutic approaches. Moreover, when combined with the TMB, a well-established biomarker for assessing response to immunotherapy (48), the DMscore demonstrated enhanced accuracy in predicting both patient survival (Figure 5F)and treatment outcomes (Figure 6I). Another interesting result is that a low-DMscore intent to be higher response rate to atezolizumab but presented a “desert” immunophenotype in IMvigor210 cohort. One possible explanation for this result is that DMscore was negatively correlated with immune cells infiltrating (Figure 5D). The desert phenotype in IMvigor210 cohort was identified by histologically CD8+ T cell, and demonstrated a low response to atezolizumab. However, CD8+ T cell only explained 4% variance of response to atezolizumab (27). In lung adenocarcinoma (49) and acute myeloid leukemia (50), DNA methylation regulators defined low-risk group was preferentially associated with TME and the sensitivity to immune response. The inter-connectedness of immune factors and the gaps in our understanding of the mechanisms of response to immune-oncology agents means that CD8+ T cell alone will be not sufficient to deliver precision medicine across all the indications for ICIs. In this study, DMscore could be an independent variable to predict the response of ICIs. In brief, the DMscore could be an important biomarker for evaluating methylation modification patterns and predicting patient response to anti-CTLA-4/PD-1 immunotherapy.

Even though we analyzed multiple datasets, this study has a few limitations. First, we included DNA 5-mc methylation-related regulators in this study; future research could integrate additional m6A regulators to further advance the understanding of epigenetic regulation in HCC. Second, there is no accessible transcriptome or clinical data for ICI therapy patients with HCC. The predictive value of the DMscore for immunotherapy efficacy in HCC needs to be validated in further studies.

Conclusions

In summary, this work investigated the comprehensive regulatory mechanisms of DNA 5-mc methylation modification in the HCC microenvironment and constructed a comprehensive scoring system for individual DNA methylation modification patterns. The DMscore serves as a valuable tool for predicting immune infiltration within the TME and refining the accuracy of immunotherapy prognosis.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA database, accession number: PRJNA762641.

Ethics statement

The study received approval from the Research Ethics Committee at the First Affiliated Hospital of Zhejiang University, and written informed consent was acquired from all participating patients.

Author contributions

JZ: Data curation, Formal Analysis, Funding acquisition, Writing – original draft. ZL: Methodology, Resources, Writing – original draft. KY: Resources, Validation, Writing – original draft. SS: Data curation, Methodology, Resources, Writing – original draft. JP: Conceptualization, Funding acquisition, Investigation, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Startup fund for Advanced Talents from Zhejiang Shuren University (grant no. 2023R038), the Shanghai Natural Science Foundation (grant no. 19ZR1431100), the Multidisciplinary Cross Research Foundation of Shanghai Jiao Tong University (grant no. ZH2018QNA42) and the Incubating Program for National Natural Science Foundation of China of Renji Hospital (grant no.RJTJ-PY-033).

Acknowledgments

The authors would like to thank The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC) and The Cancer Immunome Database (TCIA) for the availability of the data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer LJ declared a shared parent affiliation with the author JP to the handling editor at the time of the review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1333923/full#supplementary-material

Abbreviations

HCC, hepatocellular carcinoma; HBV, Hhepatitis B virus; HCV, hepatitis C virus; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; CNV, copy number variation; TMB, tumor mutation burden; DEGs, differentially expressed genes; OS, overall survival; KEGG, Kyoto Encyclopedia of Genes and Genomes; ssGSEA, single-sample gene set enrichment analysis; PCA, principal component analysis; HR, hazard ratio; TGF-β, transforming growth factor beta; TME, tumor microenvironment; Tfh, T follicular helper cell; MDSC, myeloid-derived suppressor cell; ICIs, immune checkpoint inhibitors; ORRs, objective response rate.

References

1. Klutstein M, Nejman D, Greenfield R, Cedar H. DNA methylation in cancer and aging. Cancer Res. (2016) 76:3446–50. doi: 10.1158/0008-5472.CAN-15-3278

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Kondo Y. Epigenetic cross-talk between DNA methylation and histone modifications in human cancers. Yonsei Med J. (2009) 50:455–63. doi: 10.3349/ymj.2009.50.4.455

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Nicholson TB, Veland N, Chen T. Chapter 3 - Writers, Readers, and Erasers of Epigenetic Marks. In: Gray SG, editor. Epigenetic Cancer Therapy. Academic Press, Boston (2015). p. 31–66.

Google Scholar

4. Xiao Q, Nobre A, Pineiro P, Berciano-Guerrero MA, Alba E, Cobo M, et al. Genetic and epigenetic biomarkers of immune checkpoint blockade response. J Clin Med. (2020) 9:286. doi: 10.3390/jcm9010286

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, et al. Integrated proteogenomic characterization of HBV-related hepatocellular carcinoma. Cell. (2019) 179:561–77. doi: 10.1016/j.cell.2019.08.052

PubMed Abstract | CrossRef Full Text | Google Scholar

6. El-Serag HB. Epidemiology of viral hepatitis and hepatocellular carcinoma. Gastroenterology. (2012) 142:1264–73. doi: 10.1053/j.gastro.2011.12.061

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Song MA, Kwee SA, Tiirikainen M, Hernandez BY, Okimoto G, Tsai NC, et al. Comparison of genome-scale DNA methylation profiles in hepatocellular carcinoma by viral status. Epigenetics. (2016) 11:464–74. doi: 10.1080/15592294.2016.1151586

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Brahmer J, Reckamp KL, Baas P, Crinò L, Eberhardt WE, Poddubskaya E, et al. Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer. N Engl J Med. (2015) 373:123–35. doi: 10.1056/NEJMoa1504627

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell. (2016) 165:35–44. doi: 10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Llovet JM, Castet F, Heikenwalder M, Maini MK, Mazzaferro V, Pinato DJ, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol. (2022) 19:151–72. doi: 10.1038/s41571-021-00573-2

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhao Y, Oreskovic E, Zhang Q, Lu Q, Gilman A, Lin YS, et al. Transposon-triggered innate immune response confers cancer resistance to the blind mole rat. Nat Immunol. (2021) 22:1219–30. doi: 10.1038/s41590-021-01027-8

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hong YK, Li Y, Pandit H, Li S, Pulliam Z, Zheng Q, et al. Epigenetic modulation enhances immunotherapy for hepatocellular carcinoma. Cell Immunol. (2019) 336:66–74. doi: 10.1016/j.cellimm.2018.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Li D, Chen B, Zeng Y, Wang H. UHRF1 could be a prognostic biomarker and correlated with immune cell infiltration in hepatocellular carcinoma. Int J Gen Med. (2021) 14:6769–76. doi: 10.2147/IJGM.S335016

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Wu HX, Chen YX, Wang ZX, Zhao Q, He MM, Wang YN, et al. Alteration in TET1 as potential biomarker for immune checkpoint blockade in multiple cancers. J Immunother Cancer. (2019) 7:264. doi: 10.1186/s40425-019-0737-3

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu YP, Lv L, Liu Y, Smith MD, Li WC, Tan XM, et al. Tumor suppressor TET2 promotes cancer immunity and immunotherapy efficacy. J Clin Invest. (2019) 129:4316–31. doi: 10.1172/JCI129317

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Huang P, Xu M, Han H, Zhao X, Li MD, Yang Z. Integrative analysis of epigenome and transcriptome data reveals aberrantly methylated promoters and enhancers in hepatocellular carcinoma. Front Oncol. (2021) 11:769390. doi: 10.3389/fonc.2021.769390

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Huang P, Zhang B, Zhao J, Li MD. Integrating the epigenome and transcriptome of hepatocellular carcinoma to identify systematic enhancer aberrations and establish an aberrant enhancer-related prognostic signature. Front Cell Dev Biol. (2022) 10:827657. doi: 10.3389/fcell.2022.827657

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. (2016) 44:e71. doi: 10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zhang F, Xue M, Jiang X, Yu H, Qiu Y, Yu J, et al. Identifying SLC27A5 as a potential prognostic marker of hepatocellular carcinoma by weighted gene co-expression network analysis and in vitro assays. Cancer Cell Int. (2021) 21:174. doi: 10.1186/s12935-021-01871-6

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Fujimoto A, Furuta M, Totoki Y, Tsunoda T, Kato M, Shiraishi Y, et al. Whole-genome mutational landscape and characterization of noncoding and structural mutations in liver cancer. Nat Genet. (2016) 48:500–09. doi: 10.1038/ng.3547

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

22. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–73. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. (2006) 98:262–72. doi: 10.1093/jnci/djj052

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, et al. Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures. Cancer Immunol Res. (2019) 7:737–50. doi: 10.1158/2326-6066.CIR-18-0436

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Huang X, Qiu Z, Li L, Chen B, Huang P. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in hepatocellular carcinoma. Aging (Albany NY). (2021) 13:20698–715. doi: 10.18632/aging.203456

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. (2018) 554:544–48. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–57. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Cedoz P, Prunello M, Brennan K, Gevaert O. MethylMix 2.0: an R package for identifying DNA methylation genes. Bioinformatics. (2018) 34:3044–46. doi: 10.1093/bioinformatics/bty156

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Dixon G, Pan H, Yang D, Rosen BP, Jashari T, Verma N, et al. QSER1 protects DNA methylation valleys from de novo methylation. Science. (2021) 372:eabd0875. doi: 10.1126/science.abd0875

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Llovet JM, Pinyol R, Yarchoan M, Singal AG, Marron TU, Schwartz M, et al. Adjuvant and neoadjuvant immunotherapies in hepatocellular carcinoma. Nat Rev Clin Oncol. (2024) 21:294–311. doi: 10.1038/s41571-024-00868-0

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Donne R, Lujambio A. The liver cancer immune microenvironment: Therapeutic implications for hepatocellular carcinoma. Hepatology. (2023) 77:1773–96. doi: 10.1002/hep.32740

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Varricchi G, Galdiero MR, Loffredo S, Lucarini V, Marone G, Mattei F, et al. Eosinophils: The unsung heroes in cancer? Oncoimmunology. (2018) 7:e1393134. doi: 10.1080/2162402X.2017.1393134

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Yang X, Yang C, Zhang S, Geng H, Zhu AX, Bernards R, et al. Precision treatment in advanced hepatocellular carcinoma. Cancer Cell. (2024) 42:180–97. doi: 10.1016/j.ccell.2024.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Chen J, Gong C, Mao H, Li Z, Fang Z, Chen Q, et al. E2F1/SP3/STAT6 axis is required for IL-4-induced epithelial-mesenchymal transition of colorectal cancer cells. Int J Oncol. (2018) 53:567–78. doi: 10.3892/ijo.2018.4429

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. (2017) 541:321–30. doi: 10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Li B, Yan C, Zhu J, Chen X, Fu Q, Zhang H, et al. Anti-PD-1/PD-L1 blockade immunotherapy employed in treating hepatitis B virus infection-related advanced hepatocellular carcinoma: A literature review. Front Immunol. (2020) 11:1037. doi: 10.3389/fimmu.2020.01037

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Pfister D, Núñez NG, Pinyol R, Govaere O, Pinter M, Szydlowska M, et al. NASH limits anti-tumour surveillance in immunotherapy-treated HCC. Nature. (2021) 592:450–56. doi: 10.1038/s41586-021-03362-0

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Pan Y, Liu G, Zhou F, Su B, Li Y. DNA methylation profiles in cancer diagnosis and therapeutics. Clin Exp Med. (2018) 18:1–14. doi: 10.1007/s10238-017-0467-0

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chiappinelli KB, Zahnow CA, Ahuja N, Baylin SB. Combining epigenetic and immunotherapy to combat cancer. Cancer Res. (2016) 76:1683–89. doi: 10.1158/0008-5472.CAN-15-2125

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Sun X, Zhu H, Cao R, Zhang J, Wang X. BACH1 is transcriptionally inhibited by TET1 in hepatocellular carcinoma in a microRNA-34a-dependent manner to regulate autophagy and inflammation. Pharmacol Res. (2021) 169:105611. doi: 10.1016/j.phrs.2021.105611

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zhu S, Hoshida Y. Molecular heterogeneity in hepatocellular carcinoma. Hepat Oncol. (2018) 5:10. doi: 10.2217/hep-2018-0005

CrossRef Full Text | Google Scholar

43. Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol. (2015) 15:669–82. doi: 10.1038/nri3902

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Gajewski TF, Woo SR, Zha Y, Spaapen R, Zheng Y, Corrales L, et al. Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr Opin Immunol. (2013) 25:268–76. doi: 10.1016/j.coi.2013.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gajewski TF. The next hurdle in cancer immunotherapy: overcoming the non–T-cell–inflamed tumor microenvironment. Semin Oncol. (2015) 42:663–71. doi: 10.1053/j.seminoncol.2015.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Kim JM, Chen DS. Immune escape to PD-L1/PD-1 blockade: seven steps to success (or failure). Ann Oncol. (2016) 27:1492–504. doi: 10.1093/annonc/mdw217

PubMed Abstract | CrossRef Full Text | Google Scholar

47. El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet. (2017) 389:2492–502. doi: 10.1016/S0140-6736(17)31046-2

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. (2019) 30:44–56. doi: 10.1093/annonc/mdy495

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Huang J, Huang C, Huang C, Xiang Z, Ni Y, Zeng J, et al. Comprehensive analysis reveals the prognostic and immunogenic characteristics of DNA methylation regulators in lung adenocarcinoma. Respir Res. (2024) 25:74. doi: 10.1186/s12931-024-02695-4

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Wen Q, Wang S, Hong L, Shen S, He Y, Sheng X, et al. m5C regulator-mediated methylation modification patterns and tumor microenvironment infiltration characteristics in acute myeloid leukemia. Immunity Inflammation Dis. (2024) 12:e1150. doi: 10.1002/iid3.1150

CrossRef Full Text | Google Scholar