A novel signature incorporating lipid metabolism- and immune-related genes to predict the prognosis and immune landscape in hepatocellular carcinoma

Abstract

Background: Liver hepatocellular carcinoma (LIHC) is a highly malignant tumor with high metastasis and recurrence rates. Due to the relation between lipid metabolism and the tumor immune microenvironment is constantly being elucidated, this work is carried out to produce a new prognostic gene signature that incorporates immune profiles and lipid metabolism of LIHC patients.

Methods: We used the "DEseq2" R package and the "Venn" R package to identify differentially expressed genes related to lipid metabolism (LRDGs) in LIHC. Additionally, we performed unsupervised clustering of LIHC patients based on LRDGs to identify their subgroups and immuno-infiltration and Gene Ontology (GO) enrichment analysis on the subgroups. Next, we employed multivariate, LASSO and univariate Cox regression analyses to determine variables and to create a prognostic profile on the basis of immune- and lipid metabolism-related differential genes (IRDGs and LRDGs). We separated patients into low- and high-risk groups in accordance with the best cut-off value of risk score. We conducted Decision Curve Analysis (DCA), Receiver Operating Characteristic curve analysis as a function of time as well as Survival Analysis to evaluate this signature's prognostic value. We incorporated the clinical characteristics of patients into the risk model to obtain a nomogram prognostic model. GEO14520 and ICGC-LIRI JP datasets were employed to externally confirm the accuracy and robustness of signature. The gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) were applied for investigating the underlying mechanisms. Immune infiltration analysis was implemented to examine the differences in immune between both risk groups. Single-cell RNA sequencing (scRNA-SEQ) was utilized to characterize the genes that were involved in the distribution of signature and expression characteristics of different LIHC cell types. The patients' sensitivity in both risk groups to commonly used chemotherapeutic agents and semi-inhibitory concentrations (IC50) of the drugs was assessed using the GDSC database. On the basis of the differentially expressed genes (DEGs) in the two groups, the CMAP database was adopted for the prediction of potential small-molecule compounds. Small-molecule compounds were molecularly docked with prognostic markers. Lastly, we investigated the prognostic gene expression levels in normal and LIHC tissues with immunohistochemistry (IHC) and quantitative reverse transcription polymerase chain reaction(qRT-PCR).

Results: We built and verified a prognostic signature with seven genes that incorporated immune profiles and lipid metabolism. Patients were classified as low- and high-risk groups depending on their prognostic profiles. The overall survival (OS) was markedly lower in the high-risk group as compared to low-risk group. Time-dependent ROC curves more precisely predicted patients' survival at 1, 3 and 5 years; the area under the ROC curve was 0.81 (1 year), 0.75 (3 years) and 0.77 (5 years). The DCA curves showed the value of the prognostic genes in this signature for clinical applications. We included the patients' clinical characteristics in the risk model for both multivariate and univariate Cox regression analyses, and the findings revealed that the risk model represents an independent factor that influences OS in LIHC patients. With immune analysis, GSVA and GSEA, we identified that there are remarkable differences between the two risk groups in immune pathways, lipid metabolism, tumor development, immune cell infiltration and immune microenvironment, response to immunotherapy, and sensitivity to chemotherapy. Moreover, those with higher risk scores presented greater sensitivity to the chemotherapeutic agents. Experiments in vitro further elucidated the roles of SPP1 and FLT3 in the LIHC immune microenvironment. Furthermore, four small-molecule drugs that could target LIHC were screened. In vitro qRT-PCR , IHC revealed that the SPP1,KIF18A expressions were raised in LIHC in tumor samples, whereas FLT3,SOCS2 showed the opposite trend.

Conclusions: We developed and verified a new signature comprising immune- and lipid metabolism-associated markers and to assess the prognosis and the immune status of LIHC patients. This signature can be applied to survival prediction, individualized chemotherapy, and immunotherapeutic guidance for patients with liver cancer. This study also provides potential targeted therapeutics and novel ideas for the immune evasion and progression of LIHC.

Keywords: individualized chemotherapy; lipid metabolism; liver hepatocellular carcinoma; prognostic gene signature; survival prediction; targeted chemotherapy; tumor immune microenvironment.

Figure 1

The flow diagram of the present study.

Figure 2

Exploration of lipid metabolism-related DEGs (LRDGs). (A) TCGA-LIHC volcano map of differentially expressed genes. (B) The Venn diagram displays the intersection of common genes among LIHC-related DEGs and lipid metabolism-related genes. (C) Kaplan-Meier curves for overall survival of LIHC patients in different clusters. (D) unsupervised consensus clustering heatmap. (E) The plot of the relative area changes from k=2 to 9 under the cumulative distribution function (CDF) curve. (F) Consistent CDF plot. (G) Tracing plot of clustered samples. (H, I) A bar plot (H) and chord diagram (I) showing immune-related biological processes enriched to LRDGs by GO analysis.

Figure 3

The landscape of LIHC immunity. (A) Comparison of immune cell infiltration patterns between tumor tissue and normal tissue by ssGSEA algorithm. (B) Comparison of immune cell infiltration patterns between different clusters. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 4

Construction of a prognostic signature for LIHC patients based on immune-related and lipid metabolism-related DEGs. (A) The Venn diagram displays the intersection of common genes among LIHC-related DEGs and lipid metabolism-related and immune-related genes. (B) The forestplot shows the results of hazard ratios and 95% confidence intervals of signature genes from the univariate Cox regression analysis. (C) The LASSO regression algorithm was used to select the optimal variable (l) with a 10-fold cross-validation method. (D) The solution path was plotted according to coefficients against the L1 norm. (E) The distribution of risk score, survival status, and the expression levels of coefficients in the prognostic signature. (F) The overall survival curves of LIHC patients between high-risk and low-risk groups were plotted based on the prognostic signature.

Figure 5

(A) The lollipop graph displays the variables and corresponding coefficients in the prognostic signature. (B) The forest plot shows the results of hazard ratios and 95% confidence intervals of signature genes from the multivariate Cox regression analysis. (C) Representative immunohistochemical staining images of SPP1 (antibody HPA074922, 10×), KIF18A (antibody HPA039312, 10×), SOCS2 (antibody CAB010356, 10×), and STC2 (antibody HPA045372, 10×) in normal and LIHC tissues are retrieved from The Human Protein Atlas database (HPA, https://www.proteinatlas.org/, accession date: January 2023). It should be noted that the immunohistochemistry staining of FLT3, GAL, and DGAT2L6 was absent from the HPA database. (D) The time-dependent ROC curves for different clinical characteristics in the TCGA cohort. (E) The time-dependent ROC curves for the prognostic signature in the TCGA cohort.

Figure 6

(A) Survival analysis of genes involved in the prognostic signature. (B) The genes involved in the prognostic signature expression among different tumor stages. DGAT2L6 was absent in the UALCAN database.

Figure 7

Survival curves and time-dependent ROC curves of patients with different gender (A, B), ages (C, D), and tumor stages (E, F) between high-risk and low-risk groups.

Figure 8

(A) The multivariate Cox regression model with clinical features included. (B) A nomogram model was constructed to predict the 1-year, 3-year, and 5-year overall survival of LIHC patients. (C) Calibration curves of the nomogram model for 1-year, 3-year, and 5-year overall survival. (D-F) Decision curve analysis for 1-year (D), 3-year (E), and 5-year (F) overall survival of the nomogram model.

Figure 9

GSE14520 and ICGC-LIRI JP was used for validation (A) The overall survival curves of GSE14520 patients between high-risk and low-risk groups were plotted based on the prognostic signature. (B) The time-dependent ROC curves for different clinical characteristics in the GSE14520 cohort. (C) Decision curve analysis for 3-year overall survival of the prognostic model. (D) The overall survival curves of ICGC-LIRI JP patients between high-risk and low-risk groups were plotted based on the prognostic signature. (E) The time-dependent ROC curves for different clinical characteristics in the ICGC-LIRI JP cohort. (F) Decision curve analysis for 3-year overall survival of the prognostic model.

Figure 10

(A) Gene set variation analysis (GSVA) between high and low-risk groups. (B) Gene set enrichment analysis (GSEA) between high and low-risk groups. (C, D) Gene mutation analysis between high-risk and low-risk groups.

Figure 11

Correlation analysis of the risk score and immune infiltration in LIHC patients. (A-D) Comparison of the TIDE score (A). exclusion score (B). estimate score (C). stromal score (D). between the high-risk and low-risk groups. (E) The heatmap diagram displays the Immune infiltration difference between the high-risk and low-risk groups via Cibersort Algorithm. (F) The box diagram displays the Immune infiltration difference between the high-risk and low-risk groups via Cibersort Algorithm. **p < 0.01, and ***p < 0.001, ****p < 0.0001.

Figure 12

Correlation of the prognostic signature with single-cell clusters. (A) UMAP plot of ten major cell clusters in the LIHC tumor microenvironment. (B) The distribution of the prognostic genes in cell clusters. (C) Violin plot of the prognostic signature expression at the single-cell level.

Figure 13

Sensitivity analysis of high-risk and low-risk patients to commonly used chemotherapy drugs.

Figure 14

(A-D) Molecular docking pattern of key pharmacodynamic substances and core targets. (A) Irinotecan-SPP1. (B) Idarubicin-GAL. (C) Idarubicin -STC2. (D) Irinotecan-KIF18A.: (E-H) 3D structures of small molecule drugs predicted by the PubChem open chemical database,including irinotecan (E), apilimod (F), idarubicin (G), and methoxsalen (H).

Figure 15

Expression of the prognostic genes in human. (A) IHC images of SPP1, KIF18A , FLT3 and SOCS2, in LIHC tissue and paracancerous tissue (magnification ×20). Scale bars: 100µm for 20×. N represents paracancerous tissues, and T represents LIHC tissues. (B) mRNA expression of SPP1, KIF18A , FLT3 and SOCS2 in LIHC tissues and paracancerous tissues. N represents paracancerous tissue, and T represents LIHC tissue. ****p< 0.0001.