RNA-seq data for HCC samples were selected from the TCGA dataset (https://portal.gdc.cancer.gov/), and clinical data related to the patient were downloaded simultaneously. The study included 368 HCC and 50 normal liver tissue samples. The RNA-seq data of the verification set and the corresponding clinicopathological information were downloaded from ICGC dataset (https://dcc.icgc.org/), including 232 HCC samples and 202 normal samples adjacent to cancer. We obtained 33 LARs were obtained from previously published literature (Narita et al., 2019), of which 13 belonged to the acetyltransferase family, and 20 belonged to the deacetylase class. The 33 LARs information is provided in Supplementary Table S1.

The differential expression of LARs in HCC and normal liver tissue samples was analyzed with the “limma” R package and visualized by heatmap. The protein-protein interactions (PPI) network between differentially expressed LARs were analyzed with the STRING website (http://string-db.org/). In addition, the correlation between these LARs was calculated using Pearson correlation network analysis.

In the TCGA cohort, the expression of 33 LARs were analyzed to determine the HCC subtype by the “ConsensusClusterPlus” R package (Wilkerson and Hayes, 2010). The method was based on the classical K-Means algorithm, with the Euclidean distance, iterated 50 times, and 80% of tumor samples were taken in each iteration. The number of clusters was set to 2–9, and we determine the best cluster number through the clustering score for the cumulative distribution function (CDF). According to different subgroups, the survival curve was drawn with survival R packet. Combined with clinical characteristics, heat maps were used to analyze the expression and distribution of LARs in different subgroups.

In the TCGA dataset, the survival information of HCC patients and the expression data of LARs were combined, and Cox regression analysis was performed to obtain the hazard ratio (HR) with its 95% confidence interval (CI) for each LARs. According to HR to identify protective genes (HR < 1) and dangerous genes (HR > 1), LARs with p < 0.05, were selected for the next analysis. According to the LARs related to prognosis, nine risk genes were identified using the LASSO algorithm, and the risk coefficient of each gene was obtained. We constructed a risk-scoring equation based on LAR expression: R i s k   s c o r e = ∑ i = 1 n ( C o e f i ∗ x i ) where C o e f i   refers to the risk coefficient of the gene and x i is the expression value of the gene.

Taking the median risk score as the cutoff value, HCC patients were classified into high-risk and low-risk group. Using the “survival” R package to perform Kaplan–Meier survival curve analysis, we drew a risk curve from low to high patient risk values, and the package “timeROC” was used to depict the ROC curve, calculate the area under the curve (AUC), and evaluate the accuracy, the sensitivity and the specificity of the model. Principal component analysis (PCA) and t-distribution stochastic neighbor embedding (t-SNE) methods were used to study the distribution of different populations by using R packages stats and Rtsne.

Clinicopathological characteristics and risk scores were integrated with survival data to obtain independent prognostic analysis input files, and The R package “survival” was used to perform univariate Cox and multivariate Cox regression analyses on the input files to evaluate the role of predictive model risk scores in prognostic prediction. The results were a forest map for visualization. Then, we used the “rms” R package to construct these independent prognostic factors, and constructed 1-, 3-, and 5-years forecast nomograms and corresponding calibration charts. The concordance index (C-index) was applied to evaluate the performance of the nomogram, and the calibration curve of the nomogram was used to evaluate the consistency between actual and predicted value.

Based on the TCGA database, the “CIBERSORT” R package was used to calculate the infiltration proportions of 22 immune cells in each HCC patient and displayed with a barplot, heatmap, and violin chart. The correlation between the expression and distribution of 22 infiltrating immune cells in the two groups was analyzed, and the related heat map was drawn. Red denotes a positive correlation, blue indicates a negative correlation.

Gene expression data were analyzed by Gene Set Enrichment Analysis (GSEA) software (http://www.gsea-msigdb.org/gsea/index.jsp). The related pathways and molecular mechanisms of HCC patients were evaluated using the “Hallmark” and the “KEGG” gene sets. The minimum and maximum gene set sizes were 10 and 500, respectively. The nominal (NOM) p-value < 0.05 and false discovery ratio (FDR) q value < 0.25 were considered statistically significant.

Human hepatic cell line, THLE-2 was purchased from the American Type Culture Collection (MA, USA), and HCC cell line (MHCC-97H, HepG2, LM3) was purchased from the National Collection of Authenticated Cell Cultures (Shanghai, China). The cell lines were cultured in Dulbecco’s modified Eagle medium (Gibco, TX, USA), supplemented with 10% fetal bovine serum (Gibco) in a humidified incubator (37°C, 5% CO₂).

Total RNA was extracted from cells using TRIzol reagent (Invitrogen, MA, USA). cDNA was synthesized using the PrimeScript RT reagent kit (Takara, Shiga, Japan) following the manufacturer’s instructions. The samples were amplified by qPCR using SYBR Green qPCR Master Mix (Thermo Fisher Scientific, MA, USA). The primer sequence information is provided in Supplementary Table S2.

All statistical analyses were conducted using R. 4.0.4 (https://www.r- project. org/) and SPSS Statistics software (version 25, https://www.ibm.com/products/software). The student’s t-test was used to compare the differences between the two groups. A value of p < 0.05 indicated that a statistically significant difference.

The expression characteristics of 33 LARs were analyzed using TCGA and ICGC datasets. Compared with normal liver tissues, 23 same LARs were abnormally expressed in HCC in the ICGC and TCGA datasets (Figures 1A,B) Genes marked in red indicate that they were differentially expressed in the two databases. Considering the similarity of biological functions of LARs, we analyzed the interaction and correlation between 33 LARs differentially expressed in the TCGA dataset. A PPI network analysis showed that these LARs frequently interact with each other, and HDACs and KATs have a particularly high interaction with other LARs (Figure 1C). Figure 1D shows the correlation of 33 LARs in HCC. We found that there was a significant correlation between the gene expression patterns of LARs in the same functional category, and there was a high correlation between acetyltransferase and deacetylase. In correlation analysis, there was a high correlation among KAT2A, HDAC10, SIRT6 and SIRT7, while SIRT6 was negatively correlated with the expression of KAT2B, HDAC6 and HDAC7 in HCC.

In order to establish the prognostic characteristics based on LARs, Pearson correlation coefficient analysis was carried out. Consensus clustering analysis was used to divide 368 samples into k groups (k = 2–9) (Figures 2A,B). According to the cumulative distribution function (CDF) value, k = 3 is the best number of clusters for dividing HCC queues, that is, cluster 1, cluster 2 and cluster 3 (Figure 2C). In order to understand the differences among the three subgroups in more detail, we used principal component analysis (PCA) to evaluate the classification of mRNA expression profiles, and the analysis showed that there were significant differences among the three subgroups (Figure 2D). In addition, the results of survival analysis showed that there was a great difference in OS time among the three clusters, and the OS of patients with cluster 3 was significantly lower than that of cluster 1 (Figure 2E). Subsequently, we evaluated the relationship between clusters and clinicopathological features (Figure 2F). There are differences in pathological characteristics such as grade, age, sex, and status among the three groups, and there are also differences in expression levels of LARs between the three clusters.

In the TCGA dataset, we performed univariate Cox regression analysis to explore the prognostic values of LARs. We found that 10 LARs were significantly associated with OS in HCC patients (Figure 3A). Among the 10 genes, ESCO2, HAT1, HDAC1, HDAC2, HDAC4, and HDAC11 are risk factors for HCC, with a risk ratio >1, while KAT8, HDAC6, SIRT3, and SIRT5 are protective factors, with a risk ratio <1. LASSO cox regression analysis was then used to build a prognostic model from the 10 LARs for HCC patients in the TCGA dataset. Based on the minimum standard, the prognostic features based on 9 LARs are successfully developed (Figure 3B, Supplementary Figure S1A). Nine LARs coefficients were visualized from high to low using the LASSO algorithm, as shown in Figure 3C. Detailed information is provided in Supplementary Table S3.

According to the risk score, HCC patients in the training group were divided into high-risk group and low-risk group. Survival analysis showed that patients in the high-risk group had worse survival compared with those in the low-risk group (Figure 4A). The increase in LARs-related gene expression and mortality in HCC patients were parallel to the increase in risk score (Figure 4B). The AUC values of 1-year, 2-years and 3-years calculated by TCGA dataset are 0.732, 0.701 and 0.695 respectively (Figure 4C), indicating that our LARs prognostic characteristics have a good predictive ability. PCA and t-SNE analysis revealed that distribution patterns of patients were significantly different in the different risk groups (Figures 4D,E). Then, Kaplan-Meier survival analysis was performed to investigate the prognostic value of a single LAR in HCC patients. The results showed that the expression of ESCO2, HDAC2, HDAC11, HDAC1, HAT1, HDAC6, and SIRT3 were found to be linked to the prognosis of HCC (Figures 4F–N). To evaluate the predictive value of LARs signature obtained from the training set, we used the ICGC database as external verification dataset. Using the same formula, the risk scores of the patients in the validation set were calculated and grouped. As shown in Supplementary Figures S2A,B, the survival time of patients in the high-risk group was significantly shorter than that in the low-risk group (p < 0.001). Moreover, the verification results showed that the AUC expression of LAR-related mRNA was 0.752, 0.656, and 0.685 at 1, 2, and 3 years, respectively (Supplementary Figure S2C). The PCA and t-SNE analysis successfully separated and confirmed two patient subgroups (Supplementary Figures S2D–E). Kaplan-Meier survival curve revealed that patients with high expression of ESCO2, HDAC2, HDAC11, HAT1 and HDAC6 had poorer OS expression than patients with low expression. (Supplementary Figures S2F–N), which is consistent with results from training set. These results revealed that the LARs signature is excellent in predicting the prognosis of HCC in the validation set.

According to the heat map of risk and clinical correlation, the high-risk score was positively correlated with tumor classification, stage, grade, age, and status in the TCGA dataset (Figure 5A). In the ICGC dataset, the high-risk group was positively correlated with stage and status (Figure 5B). According to the risk and clinical correlation heat map of TCGA and ICGC datasets, HDAC1, HDAC2, HDAC4, HDAC11, HAT1, and ESCO2 were suggested to be high-risk LARs. In TCGA cohort, according to univariate Cox regression analysis, staging (HR = 1.657, 95% CI, 1.353–2.031, p < 0.001) and risk score (HR = 1.083, 95% CI, 1.060–1.107, p < 0.001) were significantly correlated with OS (Figure 5C). Multivariate analysis by Cox regression showed that, after correcting for other confounding factors, staging (HR = 1.515, 95% CI, 1.220–1.881, p < 0.001) and risk score (HR = 1.078, 95% CI, 1.051–1.106, p < 0.001) were still statistically significant (Figure 5D). The results were validated using the ICGC queue (Figures 5E,F). Based on multivariate analysis by Cox regression, sex, age, stage, and risk scores were introduced into the line chart to quantitatively predict the OS of HCC patients (Figure 5G). The calibration curves of the 1-, 3-, and 5-years line diagrams were very close to the best prediction curve, and the predicted OS rate and the actual observation value were highly consistent. (Figure 5H, Supplementary Figures S1B–F).

To further immune infiltration analysis, we employed CIBERSORT to calculate the infiltration abundance of 22 immune cell types between high-and low-risk HCC samples in the TCGA dataset. The infiltration ratio of immune cells was shown in Figure 6A. We also assessed the association between 22 types of immune cells, and the heat map revealed that the proportions of different subpopulations of tumor-infiltrating immune cells were weakly to moderately correlated (Figure 6B). Next, we used the heat map (Figure 6C) and the violin map (Figure 6D) to reveal possible differences in immune cell expression between groups. The results showed that memory activated CD4⁺ T cells and macrophages M0 were higher expressed in the high-risk group, while T cells gamma delta, NK cells regulatory, Macrophages M2, T cells follicular helper, and Resting mast cells were more infiltrated in the low-risk group.

GSEA was used to study tumor markers and signal pathways in patients with low-risk and high-risk HCC in TCGA dataset. We determined that four tumor characteristics, namely E2F targets, G2/M checkpoint, mitotic spindle and Wnt/β-catenin signaling were significantly enriched in high-risk HCC patients (Figure 7A). GSEA-based KEGG pathway analysis confirmed that cell cycle, notch signaling pathway, pathways in cancer and phosphatidylinositol signaling system were enriched in the high-risk subgroup (Figure 7B).

We used RT-qPCR to verify the expression level of these prognosis-related LARs genes.

Our results showed that compared with the HCC cell line (THLE-2), the expression of HAT1, HDAC1, HDAC2, HDAC4 and HDAC11 were up-regulated in the HCC cell line (97H, HepG2, LIM3) (Figure 7C).