Identification and Analysis of Immune-Related Gene Signature in Hepatocellular Carcinoma

Abstract

Background: Hepatocellular carcinoma (HCC) originates from the hepatocytes and accounts for 90% of liver cancer. The study intends to identify novel prognostic biomarkers for predicting the prognosis of HCC patients based on TCGA and GSE14520 cohorts.

Methods: Differential analysis was employed to obtain the DEGs (Differentially Expressed Genes) of the TCGA-LIHC-TPM cohort. The lasso regression analysis was applied to build the prognosis model through using the TCGA cohort as the training group and the GSE14520 cohort as the testing group. Next, based on the prognosis model, we performed the following analyses: the survival analysis, the independent prognosis analysis, the clinical feature analysis, the mutation analysis, the immune cell infiltration analysis, the tumor microenvironment analysis, and the drug sensitivity analysis. Finally, the survival time of HCC patients was predicted by constructing nomograms.

Results: Through the lasso regression analysis, we obtained a prognosis model of ten genes including BIRC5 (baculoviral IAP repeat containing 5), CDK4 (cyclin-dependent kinase 4), DCK (deoxycytidine kinase), HSPA4 (heat shock protein family A member 4), HSP90AA1 (heat shock protein 90 α family class A member 1), PSMD2 (Proteasome 26S Subunit Ubiquitin Receptor, Non-ATPase 2), IL1RN (interleukin 1 receptor antagonist), PGF (placental growth factor), SPP1 (secreted phosphoprotein 1), and STC2 (stanniocalcin 2). First, we found that the risk score is an independent prognosis factor and is related to the clinical features of HCC patients, covering AFP (α-fetoprotein) and stage. Second, we observed that the p53 mutation was the most obvious mutation between the high-risk and low-risk groups. Third, we also discovered that the risk score is related to some immune cells, covering B cells, T cells, dendritic, macrophages, neutrophils, etc. Fourth, the high-risk group possesses a lower TIDE score, a higher expression of immune checkpoints, and higher ESTIMATE score. Finally, nomograms include the clinical features and risk signatures, displaying the clinical utility of the signature in the survival prediction of HCC patients.

Conclusions: Through the comprehensive analysis, we constructed an immune-related prognosis model to predict the survival of HCC patients. In addition to predicting the survival time of HCC patients, this model significantly correlates with the tumor microenvironment. Furthermore, we concluded that these ten immune-related genes (BIRC5, CDK4, DCK, HSPA4, HSP90AA1, PSMD2, IL1RN, PGF, SPP1, and STC2) serve as novel targets for antitumor immunity. Therefore, this study plays a significant role in exploring the clinical application of immune-related genes.

Keywords: clinical; hepatocellular carcinoma; immune; mutation; prognosis.

Figure 1

The differential expression analysis in TCGA. (A) Heatmap. (B) Volcano plot. The green dot represents the low-expression genes, the red represents the high-expression genes and the black represents no differential genes in the volcano.

Figure 2

Lasso regression analysis, ROC curve, and Kaplan–Meier survival curves. (A) Univariate Cox analysis. (B) Distribution of the LASSO coefficients. (C) The 10-fold cross-verification of variable selection in the LASSO algorithm. (D) ROC curve in the TCGA cohort. (E) ROC curve in the GSE14520 cohort. (F) Survival curves in the TCGA cohort. (G) Survival curves in the GSE14520 cohort.

Figure 3

Prognosis of the risk model in two groups. (A) Univariate prognosis analysis in the TCGA cohort. (B) Multivariate prognosis analysis in the TCGA cohort. (C) Univariate prognosis analysis in the GSE14520 cohort. (D) Multivariate prognosis analysis in the GSE14520 cohort. (E) Heatmap of the ten genes’ expression in the TCGA cohort. (F) Heat maps of the ten genes’ expression in the GSE14520 cohort. (G) The overall survival risk scores, survival time, and survival status distribution in the TCGA cohort. (H) The distribution of the overall survival risk scores, survival time, and survival status in the GSE14520 cohort.

Figure 4

Clinical feature analysis and survival analysis. ** p < 0.01; *** p < 0.001. (A) Heatmap of the clinical correlation in the TCGA cohort. (B) Heatmap of the clinical correlation in the GSE14520 cohort. (C) Box plot of the clinical correlation in the TCGA cohort. (D) Box plot of the clinical correlation in the GSE14520 cohort. (E) The OS of the ten genes in the TCGA cohort. (F) The OS of ten genes in the GSE14520 cohort.

Figure 5

Tumor mutation analysis in TCGA cohort. (A) Waterfall plots in the high-risk group. (B) Waterfall plots in the low-risk group. (C) Survival curves of two groups. (D) Survival curves of four groups.

Figure 6

Immune cell infiltration analysis. * p < 0.05; ** p < 0.01; *** p < 0.001. (A) ssGSEA scores of immune cells and immune function in the risk group. (B) Immune cell bubble plot. (C–H) Correlation between the risk score and six immune cells. (I) The expression of immune checkpoints. (J–M) Violin graph of TIDE score, MSI, dysfunction, and exclusion scores between the low- and high-risk groups, respectively. (N–P) Box graphs of ESTIMATEScore, ImmuneScore, and StromalScore between the low- and high-risk groups, respectively.

Figure 7

Drug sensitivity analysis and nomogram. *** p < 0.001. (A–F) Drug sensitivity. 5-Fluorouracil (A), VX-11e (B), and sapitinib (C) were more effective in the high-risk group. Selumetinib (D), sorafenib (E), and gemcitabine (F) were more effective in the low-risk group. (G) Nomogram in the TCGA cohort. (H) Nomogram in the GSE14520 cohort. (I) Calibration curves in the TCGA cohort. (J) Calibration curves in the GSE14520 cohort.