We analyzed 139 HBV-HCC samples from TCGA using non-negative matrix factorization (NMF) to perform unsupervised clustering. The cophenetic value and clustering heatmap were used to determine the optimal number of clusters, and we found that three clusters showed the greatest inter-group variability and the least intra-group variability. We then compared the overall survival (OS) rates of the three clusters using the “survival” and “survminer” packages in R software. The expression of genes related to programmed cell death (PCD) in the three clusters was visualized using a heatmap generated with the “pheatmap” package. We also analyzed the clinical characteristics of the three clusters using a stacked bar plot generated with the “ggplot2” package.

Analysis of immune cell infiltration was performed using TIMER2.0 (http://timer.cistrome.org), which utilizes gene expression data to estimate the abundance of various immune cell types in tumor tissues. Seven different methods were used to evaluate immune cell infiltration. Heatmap visualization of differentially expressed immune cells was generated using the “pheatmap” package. The microenvironment of each cluster was evaluated using the “estimate” package, and boxplots were generated using the “ggpubr” package to compare the stromal score, immune score, and ESTIMATE score between clusters. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu) was used to predict the response to immune checkpoint blockade therapy, and the results were visualized using violin plots and boxplots with the “ggpubr” package.

Gene set variation analysis (GSVA) was performed using the “GSVA” and “GSEABase” packages in R to calculate the pathway activity of HALLMARK gene sets in the three HBV-HCC subgroups. Heatmaps were generated to visualize the pathway activity using the “pheatmap” package. The HALLMARK gene sets were obtained from the Molecular Signatures Database (MSigDB) (http://www.gsea-msigdb.org/gsea/msigdb). Additionally, Copy number variations (CNVs) of PCD-related genes in HBV-HCC subgroups were analyzed using CNV data obtained from the UCSC Xena browser. Lollipop charts were generated to visualize the CNV variations using the “ggplot2” package.

RNA sequencing data from three subgroups of HBV-HCC were obtained and analyzed using R software. Differential expression analysis was performed using the “limma” package to identify differentially expressed genes (DEGs) between each subgroup. DEGs with an adjusted p-value < 0.05 and a log2 fold change > 1 were considered significant. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed on the DEGs for each subgroup using the “clusterProfiler” package in R. The enriched GO terms and KEGG pathways with a p-value < 0.05 were considered significant. To visualize the GO enrichment results, GO circle plots were generated using the “circlize” and “ComplexHeatmap” packages in R.

The univariate cox regression analysis was performed to identify genes significantly associated with OS with p-value<0.05. Two methods, LASSO regression and random survival forest (RSF), were used for further screening of survival-associated genes. The optimal lambda value was used to select genes in the LASSO regression, and the top 10 genes with the highest importance score based on the Gini coefficient were selected in the RSF analysis. The intersection of genes selected by the two methods was used for further analysis. A stepwise multiple cox regression was conducted to build a prognostic model using the selected genes, and hub genes were identified. Risk scores were calculated based on the expression levels and coefficients of the hub genes for distinguishing patients into high- and low-risk groups. Kaplan-Meier survival curves and receiver operating characteristic (ROC) curves were used to evaluate the performance of the model. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to explore the expression pattern of the hub genes and to visualize the clustering of patients in the high- and low-risk groups.

We used the same risk score formula derived from the stepwise multivariate Cox regression analysis in the training and validation sets. Patients were ranked according to their risk scores, and the risk score distribution and survival curves were plotted to evaluate the prognostic value of the risk score. In addition, to compare the gene expression patterns between the high- and low-risk groups, we selected the top five genes from the multivariate Cox regression analysis and compared their expression levels in the two groups. The expression levels were presented as a heatmap using the “pheatmap” package.

GSVA method was used to evaluate the relative enrichment score of 29 immune cell types and immune processes in HBV-HCC samples. The ssGSEA score of each immune cell type and immune process was calculated using the “GSVA” and “GSEABase” packages in R. The ImmunCellAI algorithm was used to evaluate the sensitivity of immune therapy for HBV-HCC patients. The immune cell score, immunotherapy exclusion score, and cytotoxic score were obtained through the ImmunCellAI web tool (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/analysis). The tumor immune dysfunction and exclusion (TIDE) algorithm was used to evaluate the MSI (Microsatellite instability, MSI) score of HBV-HCC samples. MSI is an important factor in the occurrence and development of tumors. The results were visualized using violin plots and boxplots with the “ggpubr” package.

Drug IC50 data were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database. Drug sensitivity analysis was conducted with “pRRophetic” and box plots were drawn by “ggplot2” in R software. IC50 values between the high- and low-risk groups were compared using the t-test. Drugs with significantly lower IC50 values in the low-risk group were considered potentially suitable for low-risk patients, while drugs with significantly lower IC50 values in the high-risk group were considered potentially suitable for high-risk patients.

Patients with complete clinical information and survival data were included in this study. Univariate cox regression analysis was performed to extract factors with p<0.05 with “graphics” package and construct a multivariate cox regression model with “StepReg” and “regplot” packages. The ROC curve was evaluated to assess the discrimination ability of the model using “timeROC” package. The calibration curve was plotted with “timeROC” package to evaluate the calibration of the model. The clinical prediction model was divided into high and low-risk groups based on the model with “survival” package.

Liver cancer cells were lysed using RIPA buffer (Cell Signal Technology, MA), centrifuged for the supernatant. The protein concentration was measured using bicinchoninic acid (BCA) assay (Cwbio, Beijing, China). The lysates were then diluted in loading buffer and denatured by heating at 100°C. Standard Western blot assay were performed using DLAT antibody (Proteintech, 68303) and GAPDH antibody (Abcam, ab77109) as the loading control.

Total RNA was extracted from the liver cancer cells using Trizol reagent (Invitrogen, USA) and cDNA was synthesized using the M-MLV Reverse Transcriptase Kit (Cwbio) following the manufacturer’s instructions. RT-PCR was performed using Real SYBR Mixture (Cwbio) on a Lightcycler 480 II instrument (Roche Applied Science, USA). GAPDH severed as the internal control.

Immunohistochemistry staining was performed using the streptavidin-peroxidase method according to the manufacturer’s instructions (Ultrasensitive; MaiXin, Fuzhou, China). The tissue microarray HLivH180Su09 which was related to HBV infection were incubated with an anti-DLAT antibody (mouse anti-human; dilution, 1:2000; HPA040786) at 4°C overnight, followed by the biotinylated anti-mouse IgG secondary antibody. The result of IHC were independently scored by two investigators who were blinded to the clinical data. The scores were obtained by evaluating the staining intensity and percentage of positive cells in representative areas. We used the following strategy to assess the results: intensity, 0 (no signal), 1 (weak), 2 (moderate), or 3 (high); percentage of cells, 0%-100%. We multiplied the scores of the staining intensity and percentage to obtain a final score (range 0–3). When the IHC score≥1.5, they had a high DLAT expression. When the IHC score<1.5, they were defined as low DLAT expression.

We collected 139 HBC-HCC samples from TCGA-LIHC which contained clinical information and survival outcomes. According to NMF unsupervised clustering (Figure 1A), the samples were divided into three subgroups, labeled Cluster 1, Cluster 2, and Cluster 3. Among them, Cluster 2 showed the worst Overall Survival (OS) probability (Figure 1B), and contained more overexpressed PCD-related genes in the heatmap (Figure 1C). Moreover, exploration of clinical characteristics revealed that Cluster 2 had the largest proportion of high-risk groups: more samples were at stage of G3/G4, III/IV and T3/T4 grade in histological grade, pathological stage, and T stage, respectively (Figures 1D-F).

We perform 7 algorithms to evaluate the immune infiltration of three clusters (Figure 2A). Cluster 2 displayed higher abundance of immune cells compared to the other clusters. Consistent with the above results, Cluster 2 had the highest stromal score, immune score and estimate score (Figures 2B-D), indicating that more active microenvironments existed in Cluster 2. The TIDE evaluation revealed that the exclusion score of Cluster 2 was also the highest (Figure 2E). Although Cluster 2 contained more immune cells, the cells were undergoing immune rejection and were unable to infiltrate. In the evaluation of dysfunctions, Cluster 1 received the lowest score (Figure 2F). This indicates that Cluster 1 was supposed to have the least immune rejection and dysfunction, which could be associated with better survival outcomes. Meanwhile, both Cluster 1 and Cluster 3 had a higher MSI score compared to Cluster 2, suggesting that immunotherapy was least effective in Cluster 2 (Figure 2G).

We further explored the pathway activation in three distinct clusters (Figure 3A). The heatmap revealed that metabolism-related pathways, such as fatty acid and bile acid metabolism, were significantly activated in Cluster 1. Cluster 2 exhibited activation of DNA replication pathways, including the G2M checkpoint, and tumor-related processes, such as p53 pathway, were significantly activated in Cluster 2, suggesting that Cluster 2 may have a closer relation to tumor progression. Analysis of CNVs in PCD-related genes revealed distinct patterns in the three HBV-HCC subgroups. In Cluster 1 (Figure 3B), 32 genes had more samples with amplifications in gene copy number compared to losses, while Cluster 2 and Cluster 3 had 22 and 24 genes, respectively (Figures 3C, D). These CNV variations may be correlated with the differences in pathway activation and survival outcomes among the subgroups.

We conducted GO and KEGG enrichment analyses on the differentially expressed genes within the three subgroups. The results of the GO enrichment analysis indicated that both Cluster 1 and Cluster 2 exhibited significant enrichment in GO terms related to metabolism (Figures 4A-C). However, the biological processes and signaling pathways related to metabolism were activated (GO z-scores > 0) in Cluster 1, while they were inhibited (GO z-scores < 0) in Cluster 2. KEGG pathway analysis showed that both Cluster 1 and Cluster 2 also shared significant enrichment in metabolic pathways (Figures 4D-F). However, the enriched pathways were mainly upregulated in Cluster 1, while they were downregulated in Cluster 2. These results suggest that the three subgroups have distinct metabolic states, with Cluster 1 showing activated metabolism, Cluster 2 showing inhibited metabolism, and Cluster 3 showing a different metabolic state.

We screened 20 PCD-related genes associated with OS using univariate Cox regression analysis. Two screening methods were used to identify potential genes: (1) According to LASSO regression, we selected 12 genes with the optimal lambda value (Figures 5A, B); (2) RSF analysis ranked the genes based on their importance, and we selected the top 10 genes (Figure 5C). The intersection of these two methods resulted in nine genes, which were further analyzed using multivariable Cox regression analysis. From this analysis, five genes (CHMP4C, DLAT, MMP1, NLRP6, and NOD2), were found to be associated with OS (Figure 5D). The risk score was calculated based on these five genes, and patients were divided into high-risk and low-risk groups. The KM survival curves showed that the high-risk group had significantly poorer OS than the low-risk group (p < 0.05) (Figure 5E). The ROC curves showed that the risk score had good accuracy in predicting 1-year (AUC: 0.766), 3-year (AUC: 0.804) and 5-year (AUC: 0.782) survival (Figure 5F). Additionally, PCA and tSNE analyses showed that the high-risk and low-risk groups were well separated based on their risk scores, indicating that the risk score represented the major differences in the patient samples (Figures 5G, H). These findings were consistent with those in the independent validation cohort (Figures 5I-L).

We ranked the patients according to the risk score of the training set, and found that patients with higher risk scores had significantly worse survival outcomes, indicating that the risk score was a reliable prognostic indicator (Figures 6A, B). Compared the expression levels of the five selected genes between the high-risk and low-risk groups, we found NLRP6 exhibited higher expression levels in the low-risk group while the other four genes expressed at higher levels in the high-risk group (Figure 6C). The same results were observed in the validation set (Figures 6D-F). These findings confirmed the prognostic value of the risk score and the potential clinical significance of the selected genes in HBV-HCC.

We evaluated the immune infiltration with ssGSEA analysis and discovered eight immune cell types were positively correlated with risk score, while two immune cell types were negatively correlated with the risk score (Figures 7A-J). In the high-risk group, most immune processes were significantly activated (Figure 7K), indicating a higher abundance of immune cells in this group. While low-risk group showed higher ImmunCellAI score, higher MSI score and lower exclusion score (Figures 7L-O) compared to high-risk group, indicating that the low-risk group was more likely to benefit from immunotherapy.

The IC50 values of 12 drugs were collected from GDSC. Among them, the IC50 of four drugs in the low-risk group was significantly lower than that in the high-risk group, indicating that these drugs may be more suitable for low-risk patients. On the other hand, the IC50 of 8 drugs was lower in the high-risk group, making them more suitable for high-risk patients (Supplementary Figure S1).

Combined with the risk score and clinical information (age, sex, T stage, N stage, M stage, and histological stage), univariate Cox regression analysis was performed yielding two factors T stage and risk score (p<0.05). These factors were closely associated with poor survival (Figure 8A). They were further used to construct a multivariate Cox regression model, visualized with a survival nomogram (Figure 8B). The ROC curve was drawn to assess the discrimination ability of the model (Figure 8C), with a larger AUC indicating better discrimination. The calibration curve was plotted to evaluate the calibration of the model, and the deviation between the actual curve and the ideal curve was small (Figure 8D). The clinical prediction model was further divided into high-risk and low-risk groups based on the model and significant survival differences were observed between the two groups (Figure 8E). Based on the risk score calculated by the 5 PCD-related genes and T stage of HBV-HCC tumors, we constructed a clinical model with good discrimination ability, calibration, and survival prediction.

We compared the expression of DLAT at both the mRNA and protein levels in Huh-7, and Huh-7 cells transfected with a plasmid containing the whole HBV genome (Huh-7/HBV). The results showed that DLAT levels increased after HBV transfection, both in terms of RNA and protein levels. (Figures 9A, B). Immunohistochemical staining was conducted on a tissue microarray of 76 HBV-HCC patients’ tumors and adjacent tissues. Three fields of view with high, medium, and low staining were chosen, revealing that the tumor exhibited stronger staining compared to the adjacent tissues. The immunohistochemical staining score of DLAT in cancer tissue was also significantly higher than in the adjacent tissues (p<0.001) (Figures 9C, D). Further analysis of the relationship between DLAT and the clinical characteristics of patients revealed that DLAT was associated with abnormal ALT and GGT levels (Figures 9E, F). It was speculated that DLAT is a gene associated with adverse effects in patients with HBV-HCC.