Lipid metabolism–related prognostic signature identification

Uni- and multi-variate Cox regression methods were implemented on 558 lipid metabolism–related genes to identify prognosis–related ones in the TCGA HBV–related HCC dataset with disease–free interval (DFI) survival data. HBV–related HCC showed heterogeneity in OS. Within 2 years, HBV+ HCC patients had poorer survival than HBV– patients in the TCGA HBV–related live cancer cohort, whereas after 2 years, the survival trend was the opposite (Fig. S3). Considering lipid metabolism’s important function in tumorigenesis, this study developed a lipid metabolism-based prognostic model of HBV+ HCC. In total, 59 genes were related to prognosis in uni-variate Cox regression analysis. Next, these genes were analyzed using multi-variate Cox regression and 43 genes significantly (P  <   0.05) related to survival were selected. The coefficients of multi-variate Cox regression were set as the risk parameters. Finally, these 43 genes were screened using LASSO analysis (Fig. S1) to identify the most prognostic genes, and 11 genes were identified to include in the prognostic risk model.

The mRNA expression patterns and corresponding protein expression in HBV+ liver carcinoma and normal tissue adjacent to tumor (NAT) of these 11 genes were examined. In the TCGA HBV–related HCC cohort, eight genes showed significantly different expression levels between NAT and HBV+ HCC groups. Seven of the genes (MMP1 being the exception) showed significantly lower expression in HBV+ HCC (P < 0.05; Fig. S4A). In the Gao et al. cohort, nine genes showed significantly different expression levels between HBV+ HCC and NAT. Seven genes had lower expression in HBV+ live cancer compared with NAT samples, with LPCAT3 and MMP1 being the exceptions (P < 0.05; Fig. S4B). Protein expression data were unavailable for the 11 markers in the TCGA HBV–related HCC cohort but available for 5 of the 11 markers in the Gao et al. cohort. All five markers’ protein expression levels were significantly lower in HBV+ live cancer than in NAT samples (P < 0.001; Fig. S4C).

Prognostic risk model construction and validation of predictive performance

After identifying lipid metabolism–related prognosis signatures, the prognostic risk model based on these 11 markers was constructed. For each sample, the risk value was calculated based on the product of risk-score parameters and the expression value of the 11 genes (see Materials and methods). On the basis of risk value, samples were divided into two risk groups (named high– and low–risk groups). The clinical features of the two risk groups in both cohorts were listed in Tables 1 and 2. In the TCGA HBV–related HCC patients, tumor stage distribution was significantly different between the two groups. The high–risk patients harbored larger fractions of stage II and III tumors (Table 1, 31.8 and 54.5%, respectively, P = 0.0001216) and a larger proportion of stage I in low–risk tumors (Table 1, 60%, P = 0.0001216). In Gao et al. cohort, the distribution of Barcelona clinic liver cancer stage (BCLC) and the distribution of tumor size was significantly different in the two risk groups, with larger proportions of small tumors (Table 2, P = 0.002239, <= 5.5 cm; 75%) and stage A tumors of BCLC (Table 2, P = 0.009862, 61.1%) in the low–risk patients.

Table 1 Clinical pathological characteristics of patients and the correlation between those parameters and overall survival in two risk groups for the TCGA cohort. BCLC, Barcelona clinic liver cancer stageTable 2 Clinical pathological characteristics of patients and the correlation between those parameters and overall survival in two risk groups for the Gao et al. cohort. BCLC, Barcelona clinic liver cancer stage

Subsequently, survival analysis was utilized to determine whether this risk model could effectively stratify patient survival in the TCGA HBV–related HCC and another two independent cohorts, including Gao et al. cohort and Roessler et al. cohort. Kaplan–Meier curve showed a significantly worse prognosis of high–risk patients in the TCGA HBV–related HCC dataset (Fig. 1 A, B; OS: P < 0.001, DFI: P = 0.002). The Gao et al. cohort consistently showed the same trend (Fig. 1 C; OS: P = 0.005), in which low–risk patients exhibited better survival. This study obtained 55 active HBV replication chronic patients with overall survival (OS) and RFS information from the Roessler et al. cohort. The risk model divided the cohort into two risk groups, with survival analysis showing a significant difference in RFS between the two groups. In addition, the high-risk groups had poorer survival (Fig. S5A; P = 0.033, HR = 2.12) than low-risk group. The OS showed survival probability difference among different follow-up periods (Fig. S5B; P = 0.47, HR = 1.35), in which high-risk groups showed poorer OS after about 32 months but better OS within about 32 months.

Fig. 1

Predictive performance and independent prognostic capacity of the risk model. A, B Kaplan–Meier survival curves of OS and DFI stratified by the risk model of the TCGA HBV–related HCC tumors. C OS Kaplan–Meier survival curves for HBV–related HCC tumors stratified by the risk model in the Gao et al. cohort. D, E) Uni- and multi-variate Cox regression methods results based on OS of the TCGA and Gao et al. cohorts, respectively

The risk model can independently predict prognostic for HBV–related HCC

Survival analytic results showed the prognostic capacity of the risk model. Next, to explore whether the model can independently predict survival, uni- and multi-variate Cox regression analysis were carried out in the TCGA and Gao et al. datasets. The uni-variate Cox regression method showed the tumor stage and the risk score significantly associated with OS (Fig. 1 D; risk score: P = 0.00027, hazard ratio [HR] = 3.7; tumor stage: P = 1.3E–06, HR = 5.7) in the TCGA HBV–related HCC cohort. Risk score (Fig. 1 E; P = 0.0079, HR = 3.5, 95% CI = 1.4–9) and tumor stage (Fig. 1 E; P = 0.0083, HR = 2.2) were significantly related to OS in the Gao et al. also. Multi-variate Cox regression showed risk score and prognosis remained significantly related after adjustment of other clinic factors in both the TCGA (Fig. 1 D; P = 0.026, HR = 2.3) and Gao et al. (Fig. 1 E; P = 0.0073, HR = 3.6) cohorts. Uni- and multi-variate Cox regression were also conducted for DFI in the TGCA dataset and relapse–free survival (RFS) in the Gao et al. cohort, respectively (Fig. S6). The results show that risk score can predict prognosis for HBV–related HCC, independent of other clinicopathological factors, and that lower risk scores indicate longer survival times.

Decoding TME context in two risk groups

To decode the immune-cell landscape of two different groups, the infiltration of immune–cells was quantified (see Materials and methods). Results using the xCell methodology indicated that in the low–risk patients in each dataset, endothelial cells, macrophages, M1 and M2 macrophages all accounted for larger proportions (Fig. 2 A, TCGA HBV–related HCC cohort: P = 3.43E–05, P = 0.0061, P = 0.0030, and P = 0.03; Fig. 2 C, HBV–related HCC in the Gao et al. dataset: P = 0.016, P = 0.0059, P = 9.57E–05, and P = 0.031). In the low–risk patients, Plasmacytoid dendritic cells (pDCs) harbored higher infiltration (TCGA HBV–related HCC cohort: P = 1.50E–02; Gao et al. cohort: P = 6.31E–05).

Fig. 2

Immune–cell infiltration. Immune–cell infiltration differences between two groups of the TCGA HBV–related HCC dataset (A, B) and Gao et al. cohort (C, D) quantified by xCell and ESTIMATE. * a single asterisk means P < 0.05, ** two asterisks mean P < 0.01, and *** three asterisks mean P < 0.001

According to the cibersort method, the low-risk patients harbored significantly larger proportions of CD8 T cells and dendritic-resting cells (TCGA HBV–related HCC dataset: P = 0.039, P = 0.024, respectively; Gao et al. dataset: P = 0.00018, P = 0.0057). In the low–risk patients of the Gao et al. dataset, activated mast cells and M1 macrophages harbored significantly higher expression levels (P = 0.0053, P = 0.047, respectively). M0 macrophages exhibited a converse trend, with higher proportions in the high–risk patients (TCGA HBV–related HCC cohort: P = 0.29, Gao et al. cohort: P = 0.00058; Fig. S7).

Consistently, StromalScore, ImmuneScore, and ESTIMATEScore quantified by the ESTIMATE method showed a higher score in low–risk cases (Fig. 2 B, D; TCGA HBV–related HCC cohort: P = 0.0018, P = 0.013, P = 0.003, respectively; P = 2.2E–07, P = 6.3E–10, P = 5E–10, respectively in Gao et al. cohort; Wilcoxon rank–sum test). Overall, immune cell infiltration was significantly greater in the low–risk groups.

The intratumoral immune composition was explored more comprehensively by assessing 28 immune–cell subpopulations reported in a pan–cancer immunogenomic analysis [20]. In general, the high–risk groups showed “cold” tumors with lower immune cells expression levels, whereas the low–risk patients had “hot” tumors (Fig. 3; TCGA HBV–related HCC cohort: P = 0.049, Gao et al. cohort: P = 2.821E–08; Fisher’s exact test).

Fig. 3

Heatmap of 28 immune–cell infiltration. Infiltration Heatmap of immune cells in the TCGA HBV–related HCC dataset (A) and Gao et al. (B) cohorts. The row means the type of immune cell. The column corresponds to each sample. The upper right bar chart in each panel shows the proportions of “hot” and “cold” immune states in two risk groups

Function enrichment characterizes high–risk and low–risk patients

This study performed function analysis to explore the intrinsic biological mechanisms between two risk groups and the interaction of lipid-metabolism with the immune microenvironment. The overlapping pathways were exhibited, including lipid metabolism–related and immuno–related pathways among all the significant enrichment pathways in the two datasets (Fig. 4). The detailed and complete results are available in Fig. S8. Collectively, the results revealed that in the high–risk groups, “Homologous recombination”, “Cell cycle”, “DNA replication”, “Mismatch repair” and “Base excision repair” pathways significantly enriched (Fig. 4). “Antigen processing and presentation”, “Leukocyte transendothelial migration”, and “Th17 cell differentiation” pathways were shared among low–risk patients in each dataset (Fig. 4). In high–risk patients of both cohorts, results of Hallmark analysis also showed the “DNA repair” pathway significantly enriched. Metabolic pathways consisting of “Fatty acid biosynthesis”, “Fatty acid degradation”, “Fatty acid metabolism”, “PPAR signaling pathway”, and “Biosynthesis of unsaturated fatty acids” significantly enriched in the low–risk group of the TCGA HBV–related HCC cohort (Fig. 4). Additionally, pathways of immune, including “Leukocyte transendothelial migration”, “Antigen processing and presentation”, and “Th17-cell differentiation” showed significant enrichment in low–risk cases of the TCGA HBV–related HCC dataset (Fig. 4). In the Gao et al. dataset, low–risk enriched more immune pathways, involving “Leukocyte transendothelial migration”, “PDL1 expression and PD1 checkpoint pathway in cancer”, “Th1 and Th2 cell differentiation”, “B cell receptor”, “Natural killer cell mediated cytotoxicity”, “Th17 cell differentiation”, and “T cell receptor”, and “Antigen processing and presentation” (Fig. 4).

Fig. 4

Functional enrichment analysis. A, B Representative and significantly enriched KEGG pathways in two risk groups of the TCGA and Gao et al. cohorts. C, D Hallmark pathway analysis results in two risk groups of the TCGA and Gao et al. cohorts. Green and blue indicate enriched pathways in high– and low-risk patients with significant P

The immune genes’ expression [23] was also examined in two groups. A heatmap (Fig. 5 A) outlines the expression landscape of those genes in the two cohorts, and boxplots show the significantly expressed results between two risk groups (Fig. 5 B, P < 0.05, TCGA HBV–induced HCC cohort; Fig. S9, P < 0.05, Gao et al. cohort; Wilcoxon rank sum test). In line with the immune–related pathways enrichment results observed in low–risk groups, 17 immune–related genes, including a co–inhibitor (SLAMF7), ligands (CD40LG, CXCL9), a receptor (BTLA), a cell adhesion (SELP), molecules involved in antigen presentation, ENTPD1, GZMA, and PRF1 expressed significantly higher in low-risk patients of both cohorts (Fig. 5 B, Fig. S9). In low–risk groups, the immune checkpoint gene BTLA had significantly elevated expression in both cohorts. In the Gao et al. cohort, PDCD1 and CD274 had significantly increased expression (Fig. S9).

Fig. 5

Immune gene expression differences between two risk groups. A Average expression heatmap of 78 immune genes. The row and column represent a kind of gene and a group. B Box plots of genes with significantly different expression levels between two risk groups of the TCGA HBV–related HCC cohort

More genetic mutations occurred in high–risk groups

Further exploration of gene mutation levels was conducted in HBV+ HCC. As indicated in Fig. 6, the 20 most frequently altered genes were identified. Mutation of TP53 and LRP1B was identified in 40.6 and 9.9%, respectively, in the cohort of TCGA. In the high–risk patients, TP53 and LRP1B mutations rate had significantly higher levels (63.6 and 27.3%, respectively) (Fig. 6; P < 0.05). TTN and DNAH8 showed mutation frequencies of 28.7 and 8.9%, respectively, in the whole cohort. Similar to TP53 and LRP1B, the mutation rate of TTN and DNAH8 were significantly higher in high–risk patients. To further explore mutation numbers and genetic heterogeneity between the two groups, we examined TMB and MATH values and high–risk cases showed higher median TMB and MATH values (Fig. S10A–B; median TMB: 2.29 vs. 2.16, median MATH value: 86.5 vs. 83), but the differences were not significant.

Fig. 6

The gene mutation landscape of the TCGA HBV–related HCC cohort

Mutation landscape of two risk groups. The top bar chart shows the gene mutation counts of each sample. The table on the left indicates gene mutation frequency in two risk groups and the whole cohort. The heatmap portrays the gene mutation landscape, in which different mutation types were annotated with different colors. The bar plot next to the heatmap indicates the mutation type proportions in all samples for each gene. The bar graph on the far right exhibits mutation proportion in two risk groups for each gene. The bottom bar chart shows clinical characteristics. * a single asterisk means P < 0.05, ** two asterisks mean P < 0.01.