Data on RNA-sequencing, including 424 cases (50 normal tissues and 374 tumor samples, FPKM value), somatic mutation and copy number variation (CNV) was obtained from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) database. In addition, RNA-sequence data were also obtained from GSE76427 (52 normal tissues and 115 tumor tissues). The FPKM values were converted to transcripts per kilobase million (TPM) values. Batch effect was corrected using “combat” algorithm based on SVA R package. This study combined RNA-sequencing data and clinical information from the GSE76427 and the International Cancer Genome Consortium (ICGC). We extracted 720 epigenetic-related genes from EpiFactors, a database for epigenetic factors, corresponding genes and products (14) Table S1 . In total, 200 inflammatory response-related genes (IRRGs) were downloaded from the gene set enrichment analysis (GSEA) hallmark inflammatory response gene set Table S2 . Project Achilles uses genome-scale CRISPR-Cas9 tool for gene knock out to identify candidate genes critical for cancer survival (15). Essential genes in HCC were identified using genome-wide CRISPR screening from the dependency map (DepMap) portal. Dependency scores for approximately 17,000 candidate genes were calculated using the CERES algorithm (15). Candidate genes were determined as crucial genes with a CERES score of <−1 across 75% of HCC cell lines. Further, the link between the identified candidate genes, the ERGs and the IRRGs was explored. Differential expression of genes (DEGs) was determined at FDR<0.01 and |log₂FC|≥1.5 using R software, limma package. RNA-seq analyses based on the adequate tumor purity and the tumor purity estimates by single sample gene set enrichment analysis (ssGSEA) for tumor tissue samples from HCC.

Unsupervised cluster analysis was performed based on DEGs to identify different epigenetic modification patterns and classify patients for further analysis. The consensus clustering algorithm determines the number and stability of the clusters. Consensus clustering was performed using the ConsenSuClusterPlus R package and repeated 1000 times to ensure the stability of the classification (16). Further, Gene Set Variation Analysis (GSVA) was performed using the “GSVA” R package to determine differences in biological processes between the epigenetic modification models. The gene set “c2.cp.kegg. V6.4. symbols” was obtained from MSigDB database for GSVA analysis. GSVA is a nonparametric, unsupervised method for estimating gene set enrichment variation from samples of expressed datasets.

The epigenetic score was used to evaluate individual epigenetic modifications with HCC outcomes. The epigenetic signature was constructed as follows: First, DEGs were standardized in all HCC samples, and the overlapping genes were identified. The patients were divided into several groups for subsequent analysis using unsupervised cluster analysis. The number and stability of the gene clusters was determined using consensus clustering algorithm. Second, univariate Cox regression analysis was used to determine the prognostic genes. Finally, principal component analysis (PCA) was used to determine the principal components. The epigenetic score was then determined using the formula (16): epigenetic score =Σ (PC1i + PC2i). Where i refers to the expression value of phenotype ERGs.

Univariate and multivariate Cox regression analyses were used to evaluate the prognostic value of ERGs. Significant genes in the univariate and multivariate Cox regression analyses were selected as the characteristic genes. The risk score for individual patients was calculated as follows: e^{sum (gene's expression×coefficient)}. The patients were then divided into high- or low-risk groups based on the median risk score cut-off value. Receiver operating characteristic (ROC) curves were used to predict the accuracy of the prognostic signatures. Kaplan–Meier survival curves were used to examine the effect of the signature on survival. The TCGA and GSE76427 databases were used in the training set, while ICGC data sets were used in the testing set. The prognostic signature for IRRGs was constructed similarly to the ERG signature.

Drug sensitivity of the high and low-risk groups was determined using the “pRRophetic” and “ISOpureR” packages (17). Expression profile data and somatic mutation data of human cancer cell lines (CCLs) were obtained from the Broad Institute-Cancer Cell Line Encyclopedia (CCLE) project based on the DepMap portal. Drug sensitivity data of the CCLs were obtained from the Cancer Therapeutics Response Portal (CTRP) and PRISM Repurposing dataset. The CTRP and PRISM datasets provide the area under the dose-response curve (area under the curve—AUC) values as a measure of drug sensitivity, with lower AUC values suggesting increased sensitivity. The detailed data processing flow and algorithm reference Chen et al.’s study (17). In addition, heatmaps were used to depict principal component differences of immune cells between the high- and low-risk groups based on XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT, and CIBERSORT-ABS algorithms.

All statistical analyses were performed in R statistical software version 3.6.2. Normally and not-normally distributed data were analyzed using the unpaired student's t-test and Wilcoxon rank-sum test, respectively. The relationship between the molecular signature and the clinicopathological parameters was assessed using the chi-square test. Logistic regression analyses were employed to determine whether the signature was an independent prognostic factor. The “surv-cutpoint” function was used to repeatedly test all possible cut points to find the maximum rank statistic, classified as the epigenetic score. The patients were divided into high and low-risk groups based on the maximum selected log-rank statistic to reduce the batch effect. The maftools package in R was used to analyze mutations in between the high- and low-risk groups. A P-value <0.05 was considered statistically significant.

Figure 1 shows the study flow chart. The clinicopathological characteristics of the patients are presented in Table 1 . The enriched ERGs signaling pathways were explored using GSVA analysis. The results showed that the ERGs were mainly enriched in cancer and immune-related pathways such as PI3K/AKT/mTOR, Wnt/beta-catenin, mTORC1, IL2/STAT5, TGF-beta, and TNF-alpha signaling pathways Figure 2A . Further, transcriptome data obtained from the TCGA and GEO databases were merged to obtain the DEGs of the 720 ERGs. Finally, using the unsupervised clustering method, three different clusters were identified, consisting of 162 cases of cluster-A, 215 cases of cluster-B, and 161 cases of cluster-C. GSVA enrichment analysis was then used to investigate the biological function of the different clusters. Cluster-A was mainly enriched in cancer-related pathways such as mTOR and ERBB signaling, and cell cycle progression Figure 2B . In contrast, Cluster-B was associated with metabolism pathways such as linoleic acid, arachidonic acid, and ether butanoate metabolism Figure 2C . Kaplan–Meier survival analysis showed that cluster-A was associated with poor survival Figure 2D . Therefore, we hypothesized that cluster A could affect epigenetics and other cancer-related processes, resulting in a poor prognosis. We then applied the ssGSEA to contrast immune cell differences between the three clusters Figure 2E . The results showed significant differences between cluster A and the other clusters in dendritic cells (DCs), neutrophils, MHC I, and regulatory T cells (Treg), suggesting that cluster A could play a key role in tumor immunity. Principal component analysis of the transcriptome profiles in the three clusters showed significant transcriptome differences Figure 3A . Further, the DEGs in each cluster were identified, and a Venn diagram was constructed to identify co-expressed genes. Finally, 5087 co-expressed genes in the three clusters were identified using the limma package Figure 3B . The KEGG analysis showed that the co-expressed genes were mainly enriched in pathways related to histone, cell cycle, DNA, and RNA functions Figures 3C, D . These findings indicate that epigenetic modifications play a key role in tumor development.

Unsupervised cluster analysis based on the acquired 5087 epigenetic-related genes was used to classify patients into different groups. Three distinct clusters 1/2/3 were obtained. The analysis showed that the three gene clusters had different characteristic genes. Cluster 2 was associated with the worst prognosis Figure 4A . Considering the individual heterogeneity, an epigenetic score was established to quantify the epigenetic pattern in each patient based on the phenotypic related genes. A ggalluvial diagram was used to visualize the relationship between the epigenetic score, epigenetic cluster and gene cluster Figure 4B . The relationship between immune-cell infiltration and the epigenetic score is presented in Figure 4C . The patients were divided into low or high-risk groups based on the median cut-off value. A low epigenetic score was associated with poor survival Figure 4D . The Kruskal-Wallis test was used to determine significant differences between the epigenetic score, gene clusters, and epigenetic clusters. A low epigenetic score, gene cluster 2 and epigenetic cluster A were associated with poor survival Figures 4E, F . Previous studies have reported that tumor mutation burden (TMB) plays an important role in tumor prognosis. Thus, we explored the relationship between the high/low epigenetic score groups and TMB. The results showed that the high TMB group had a poorer prognosis. In addition, low/high TMB combined with a low epigenetic score was also associated with poor outcomes Figures 4G, H .

Further, the maftools package was used to analyze differences in somatic mutations between the low and high epigenetic score groups. The results indicated that the low epigenetic score group had increased TMB than the high epigenetic score group. TP53 was the most common mutated gene Figures 5A, B . In addition, the relationship between the clinical-pathological characteristics and the epigenetic score was also explored. The results showed that a low score was associated with advanced tumor Figure S1 . Finally, the immunophenoscore (IPS) based on The Cancer Immunome Atlas (TCIA) database was used to predict responsiveness to CTLA-4 and PD-1. The results showed that the high- and low-score groups had a good response to PD-1 and CTLA-4 Figures 5C–F , which could provide useful insights for further exploration.

Univariate Cox analysis identified 57 genes that were significantly associated with survival. A correlation coefficient for each gene in the model was calculated to determine possible collinearity between the genes Table S3 . As a result, nine differentially expressed genes (ANP32B, ASF1A, BCORL1, BMI1, BUB1, CBX2, CBX3, CDK1, and CDK5) were selected as independent prognostic factors based on lasso regression in the TCGA-LIHC and ICGC databases Table S4 . Moreover, the risk score for each patient was calculated, and the cohort was divided into two groups (high and low risk) based on the median cut-off value. Kaplan-Meier survival curves showed that the high-risk group had poorer survival than the low-risk group Figure S2A . In addition, ROC curves were used to determine whether the expression pattern could be used as an early predictor of HCC. The ROC curves showed AUC of 0.788, 0.699 and 0.722 for one, three- and five-years survival, respectively Figure S2B . Further, we established a risk survival status chart of the patients. The chart showed that the number of patients dying increased as the risk score increased Figure S2C . Moreover, the risk score was significantly associated with the tumor stage and T stage Figure S2D .

Analyses of the CTRP and PRISM-derived drug response data yielded four CTRP-derived compounds (SB-743921, GSK461364, gemcitabine, and paclitaxel) and two PRISM-derived compounds (dolastatin-10 and LY2606368). The high-risk score group had lower estimated AUC values Figures 6A, B . The relationship between the nine genes and drug sensitivity was analyzed using the cellminer database [a relational database and query tool for the NCI-60 cancer cell lines (18)]. The results suggested that CBX2 has the connection with most drugs such as dasatinib, acrichine, and nelarabine Figure S3 .

The relationship between ERGs groups and other immune and molecular subtypes were also explored. A total of 355 immune samples were further classified according to a pan-patient immune subtype (19) Table S5 . As shown in Figure 6C , there were more C3 and C4 subtypes in the low ERGs subgroup. However, there were more C4 subtypes in the high ERGs subgroup (P=0.001, chi-square test). Further, Tumor Immune Dysfunction and Exclusion (TIDE) was used to evaluate the clinical efficacy of immunotherapy in different ERGs subgroups. A higher TIDE prediction score represents a higher likelihood of immune evasion, indicating that patients are less likely to benefit from immune checkpoint (ICI) therapy (20). In this study, patients with high ERGs benefited more from ICI treatment than patients with low ERGs Figure S4A . In addition, we found that the high ERGs subgroup had a higher T cell exclusion score. However, the low ERGs subgroup had a higher T cell dysfunction, with no differences in microsatellite instability (MSI) score between the two subgroups Figures S4B-D . Meanwhile, the ROC analysis suggested that the predictive value of the risk model was better than that of the 18-gene T-cell-inflamed signature (TIS) and TIDE models (21, 22) Figure S4E . Further, we explored differences in the expression of ICs between the two groups. The results showed significant differences in expression of PDCD-1 (PD-1) and CTLA4 between the two groups Figure S5A . In addition, there were significant differences in the expression of RBM15, YTHDC1, and TDHDC2 between the high and low-risk groups Figure S5B . The heatmap showing immune responses based on different algorithms is shown in Figure 7 . The results demonstrated that most immune cells showed a trend of high expression in the high-risk group when compared to low-risk group.

Univariate Cox analysis identified 50 genes that were significantly associated with survival Table S6 . Further, multivariate Cox regression analyses identified 18 independent prognostic genes based on lasso regression Table S7 . After that, the IRRGs signature was constructed. Kaplan-Meier survival curves showed that the high-risk group had poorer survival and an AUC value of 0.817, 0.769, and 0.754 for one, three, and five years survival, respectively. Univariate and multivariate Cox analysis revealed that the signature (HR: 3.797, 95CI: 2.744-5.255) and tumor stage (HR: 1.399, 95CI: 1.125-1.739) were independent prognosis factors of OS in HCC patients Figure S6 . The results of the training set in the present study were consistent with results of the testing set in ICGC Figure S7 . These results suggest that the IRRGs signature was stable and accurate and could be applied in the clinical management of HCC patients.

To determine key genes responsible for HCC malignancy, we explored genome-wide CRISPR-based loss-of function screens derived from DepMap. In total, 640 genes were identified as essential in the survival of HCC cell lines Table S8 . The heatmap showing the top 12 genes (CDC20, TOP2A, BIRC5, RRM2, CCNA2, MCM2, BOP1, SPC24, CCT3, MCM6, CENPW, and CDK1) is shown in Figure 8A . Finally, we evaluated the relationship between the ERGs, IRRGs, and CRISPR genes. The protein-protein interaction (PPI) network demonstrated that IRRGs PSEN1 was linked to most ERGs and CRISPR genes such as CDK1, TOP2A, CBX2, and CBX3 Figure 8B . Most genes were also shown to have a strong correlation with the immune cells Figure 8C . In addition, PSEN1 was over expression in HCC and high expression of PSEN1 indicated bad prognosis based on the GEPIA database Figures 9A, B . Immunohistochemistry results from the HPA database to illustrate that PSEN1 were significantly increased in tumor tissue Figures 9C, D .