Construction of a risk score prognosis model based on hepatocellular carcinoma microenvironment

Fa-Peng Zhang, Yi-Pei Huang, Wei-Xin Luo, Wan-Yu Deng, Chao-Qun Liu, Lei-Bo Xu, Chao Liu

Abstract BACKGROUND Hepatocellular carcinoma (HCC) is a common cancer with a poor prognosis.Previous studies revealed that the tumor microenvironment (TME) plays an important role in HCC progression, recurrence, and metastasis, leading to poor prognosis. However, the effects of genes involved in TME on the prognosis of HCC patients remain unclear. Here, we investigated the HCC microenvironment to identify prognostic genes for HCC.AIM To identify a robust gene signature associated with the HCC microenvironment to improve prognosis prediction of HCC.METHODS We computed the immune/stromal scores of HCC patients obtained from The Cancer Genome Atlas based on the ESTIMATE algorithm. Additionally, a risk score model was established based on Differentially Expressed Genes (DEGs)between high and lowimmune/stromal score patients.RESULTS The risk score model consisting of eight genes was constructed and validated in the HCC patients. The patients were divided into high- or low-risk groups. The genes (Disabled homolog 2, Musculin, C-X-C motif chemokine ligand 8, Galectin 3, B-cell-activating transcription factor, Killer cell lectin like receptor B1, Endoglin and adenomatosis polyposis coli tumor suppressor) involved in our risk score model were considered to be potential immunotherapy targets, and they may provide better performance in combination. Functional enrichment analysis showed that the immune response and T cell receptor signaling pathway represented the major function and pathway, respectively, related to the immune-related genes in the DEGs between high- and low-risk groups. The receiver operating characteristic (ROC) curve analysis confirmed the good potency of the risk score prognostic model. Moreover, we validated the risk score model using the International Cancer Genome Consortium and the Gene Expression Omnibus database. A nomogram was established to predict the overall survival of HCC patients.CONCLUSION The risk score model and the nomogram will benefit HCC patients through personalized immunotherapy.interest related to this article.

Key words: Hepatocellular carcinoma; Prognostic model; Immune related gene;Microenvironment; Risk score; Overall survival

INTRODUCTION

Hepatocellular carcinoma (HCC) is the sixth most common cancer worldwide and is the main cause of cancer-related death[1]. Although there are multiple methods to treat HCC, including surgical resection, liver transplantation, radiofrequency ablation and chemotherapy, the efficacy is limited by high recurrence rates and low rates of surgery and transplants because this cancer is usually diagnosed in a late stage[2-4].Recently, immunotherapy has emerged as a novel and effective therapy and is being applied in various tumors including HCC[5]. In particular, treatment targeting immune checkpoints has achieved success and improved patient survival[6-8].However, only a small number of patients receiving immunotherapy treatment responded to the treatment due to the immunosuppressive microenvironment. Hence,it is necessary to investigate biomarkers that enable us to predict the benefit of immunotherapy, which may help in clinical decision making for individualized treatment.

The HCC microenvironment includes various cells, such as hepatic stellate cells,cancer-associated fibroblasts, endothelial cells, neuroendocrine cells, immune cells,bone-marrow derived cells, and the extracellular matrix (ECM), which plays a crucial role in tumor initiation, progression, drug resistance and immune escaping[9-15].Immune cell infiltration is closely related to the survival of patients[16-20], indicating that understanding and reshaping the tumor microenvironment (TME) may improve the efficacy of cancer treatments in the future[21-23]. ESTIMATE was an algorithm designed by Yoshihara et al[24]to calculate immune and stromal scores based on the gene expression profile. Subsequent studies have applied the ESTIMATE algorithm in multiple cancers such as prostate cancer[25], breast cancer[26], and colon cancer[27],showing the capability of such big-data based algorithms, although the efficacy on HCC has not been verified.

In our study, we investigated the TME and the gene expression profile of HCC to construct a risk score prognostic model for HCC based on The Cancer Genome Atlas(TCGA) database. Furthermore, we have validated this model using the International Cancer Genome Consortium (ICGC) and Gene Expression Omnibus (GEO) databases.

MATERIALS AND METHODS

Gene expression datasets

The data of this study was mainly obtained from public databases. The transcriptional profiles and clinical materials from HCC patients were downloaded from the TCGA website (https://cancergenome.nih.gov/). HCC patients with R0 surgical resection were chosen, who did not receive other treatment for their disease and had a survival time of more than 1 mo. Among these patients, 361 HCC samples with complete transcriptional data and the corresponding clinical information were selected for consequent analyses. Immune and stromal scores were calculated by applying the ESTIMATE algorithm to the mRNA expression data[24]. For further verification, the code of LIRI-JP (n = 232) obtained from the ICGC database (https://icgc.org/) and the dataset GSE14520 (n = 221) downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database were selected for validation. The data downloaded from the TCGA, ICGC, and GEO databases were publicly available and accessible. The present study was conducted following pertinent guidelines and regulations approved by the TCGA, ICGC, and GEO.

Differentially expressed genes analysis

To select the intersection genes, 361 HCC patients obtained from the TCGA dataset were divided into high and low immune/stromal score groups according to the ESTIMATE results. The differentially expressed genes (DEGs) were identified using the package limma[28]in R software (Version 3.6.1; https://www.r-project.org/), and the cutoffs were fold change ＞ 1.5 and adjust P ＜ 0.05. Volcano plot and heatmaps were generated using the ggplot2 and pheatmap package in R software, respectively.

Overall survival curve

Kaplan-Meier (K-M) plots were generated to illustrate the correlation of immune/stromal scores with patients' overall survival (OS). The relationship was tested by the log-rank test.

Functional enrichment analysis of DEGs

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)pathway enrichment analyses of DEGs were performed using the Database for Annotation, Visualization and Integrated Discovery (Version: 6.8; https:// david.ncifcrf.gov/) to investigate the potential function of the DEGs[29]. Significant biological processes and pathways were presented using the ggplot2 R packages.

Construction and validation of the risk score prognostic model

Univariate, the Least Absolute Shrinkage and Selection Operator (LASSO) and multivariate Cox regression analyses were performed to explore the relationship between DEGs and patients' OS. In the univariate Cox regression analysis, P ＜ 0.05 was considered significant. To further narrow down correlated genes, the LASSO with L1-penalty was applied[30]. Based on the LASSO analysis, the pivotal genes were extracted from DEGs which were regarded as significant in the univariate Cox regression analysis. Then, a sub-selection of genes associated with patient prognosis was identified[31]. LASSO Cox regression analysis was performed using the glmnet R package (Version: 2.0). To evaluate the contribution of each gene to prognosis, the multivariate Cox regression analysis was performed. A stepwise method was used to further determine the best rprognostic model. Finally, eight genes were selected to construct a risk score prognosis model. HCC patients were divided into low- and high-risk groups based on the median risk score. The K-M survival curves for the cases with low or high risk were performed. The predictive ability of the risk score prognosis model was assessed by the survival receiver operating characteristic (ROC)package in R software. The concordance index (C-index) was calculated to investigate the risk score prognostic model performance[32]. Then, the risk score prognosis model was verified using the ICGC and GEO dataset, respectively.

Estimated immune cell type fractions

CIBERSORT is a gene expression-based deconvolution algorithm to describe the cell constitution of tissues[33]. LM22 is defined as “barcodes” with 547 gene expression signatures that distinguish 22 human hematopoietic cell phenotypes, including plasma cells, myeloid subsets, seven T cell types, naive and memory B cells and natural killer (NK) cells. We used CIBERSORT in combination with LM22 to sort the portions of 22 human immune cell types in HCC samples. For each sample, the sum of all estimate immune cell type fractions equals to 1[34].

Independence of the risk score prognostic model

Among 361 HCC patients with survival data, 289 patients with full clinical parameters, including gender, age, histologic grade, pathologic stage and vascular invasion, were subjected to consequent analyses. Univariate and multivariate Cox regression analyses were performed to assess the predictive ability of the risk score prognostic model for HCC patients. All statistical tests were two-sided, and P values＜ 0.05 were considered as statistically significant.

Construction and validation of the nomogram

The nomogram is widely utilized to predict the prognosis of cancer patients'prognosis[35]. In the present study, a nomogram was generated based on the independent prognostic factors identified by multivariate analysis to investigate the probability of 1-, 3-, and 5-year OS of HCC patients. The nomogram and calibration plots were generated using the rms R package (Version: 5.1-3). The calibration curve of the nomogram was drawn to evaluate the nomogram prediction possibilities against the observed rates[36].

Statistical analysis

The R software v3.6.1 (R Fo undation for Statistical Computing, Vienna, Austria) and GraphPad Prism v8.0 (GraphPad Software Inc., United States) were used for statistical analyses. All statistical tests were two-tailed with a statistical significance level set at 0.05.

RESULTS

Immune scores and stromal scores are significantly related to HCC prognosis

The flow chart of the study procedure is presented in Figure 1. In this study, 361 HCC patients with complete gene expression data and the corresponding clinical information were downloaded from the TCGA database for subsequent analysis. As shown in Supplementary Table 1, patients were 59.65 ± 13.33 years old, including 116 females (32.1%) and 245 females (67.9%). According to the ESTIMATE algorithm,immune scores ranged from -1209.16 and 2934.36, and stromal scores ranged from -1741.56 to 1195.07. Then, we investigated the relationship between immune/stromal scores and clinical characteristics. The result showed that the immune score was significantly negatively correlated with the pathologic stage (P = 0.032) (Figure 2B),but not with histologic grade (P = 0.968) (Figure 2A). In contrast, the stromal score was significantly negatively related to histologic stage (P = 0.008) (Figure 2C), but not with pathologic grade (P = 0.329) (Figure 2D).

To investigate the influence of immune/stromal scores on prognosis, we divided 316 HCC patients into high- and low-score groups and constructed K-M curves. We found that the immune scores and stromal scores were significantly positively correlated with OS (Figure 2E and F). Overall, we found that the immune and stromal scores were significantly correlated with poor prognosis.

Comparison of gene expression profiles with immune scores and stromal scores in HCC

To explore the difference of gene expression profiles between high- and lowimmune/stromal score groups, we compared 361 HCC cases obtained from the TCGA database. Volcano maps in Supplementary Figure 1 showed the differential gene results of the low vs high score group differential gene results. Heatmaps showed the most different 100 gene expression profiles of cases belonging to low and high immune scores/stromal scores groups. In comparison to the high immune scores group, 689 genes were downregulated and 106 genes were upregulated in the lowimmune score group (Supplementary Table 2). Similarly, there were 906 downregulated genes and 91 upregulated genes in the low stromal - score group(Supplementary Table 3). In addition, Venn diagrams (Figure 3A) showed that 13 genes (Supplementary Table 4) were commonly upregulated in the low score groups,and 473 genes were commonly downregulated (Supplementary Table 5).

Figure 1 Overall design of the present study. TCGA: The Cancer Genome Atlas database; HCC: Hepatocellular carcinoma; LASSO: Least absolute shrinkage and selection operator; GEO: Gene Expression Omnibus databases;ICGC: International Cancer Genome Consortium database.

Furthermore, the potential functions of the DEGs were evaluated by GO analysis and KEGG pathway. The top 5 GO terms identified included immune response,inflammatory response, cell adhesion, chemotaxis and extracellular matrix organization (Supplementary Figure 2). The KEGG pathway analysis revealed that the DEGs were enriched in the tumor necrosis factor signaling pathway, cytokinecytokine receptor interaction, complement and coagulation cascades, chemokine signaling pathway, and cell adhesion molecules (Figure 3B).

Construction of a risk score model and evaluation of its predictive ability in the TCGA HCC cohort

To explore the potential roles of DEGs in OS, a univariate analysis was performed.The results showed that 60 of the 486 intersection DEGs were significantly correlated with OS in HCC patients (Supplementary Table 6). To screen the most prognostic genes, LASSO analysis was performed. The contributions of 60 intersection DEGS were weighted by their relative coefficients (Supplementary Figure 3A and B). As a result, 13 genes were selected (Supplementary Table 7). Lastly, a multivariate analysis was performed, and 8 genes were chosen to establish a risk score model (Supplementary Figure 3C and Supplementary Table 8), and the final risk score formula was as follows: Risk score = [0.184 × expression level of Disabled homolog 2 (DAB2)] + [0.102× expression level of Musculin (MSC)] + [0.118 × expression level of C-X-C motif chemokine ligand 8 (CXCL8)] + [0.147 × expression level of Galectin 3 (LGALS3)] +[0.147 × expression level of B-cell-activating transcription factor (BATF)] + [ 0.48 ×expression level of killer cell lectin like receptor B1 (KLRB1)] + [ 0.393 × expression level of Endoglin (ENG)] + [ 0.113 × expression level of Adenomatosis polyposis coli tumor suppressor (APCS)]. We then calculated the risk score according to the formula and divided the patients into high- or low-risk groups based on the median risk score.According to the K-M analysis (Figure 4A), there was a significant difference in patients' OS between high- and low-risk groups, and patients in the high-risk group had significantly shorter OS than those in the low-risk group. Figure 4B shows the distribution of the risk score and gene expression data. The prediction capability of the risk score model was evaluated by calculating the area under the ROC curve(AUC) (Figure 4C). The AUC of the prognostic model for OS was 0.778, 0.754 and 0.75 for the first, third, and fifth years, respectively, suggesting that the risk score model had good performance.

Figure 2 Immune and stromal scores are associated with clinical characteristics and overall survival of hepatocellular carcinoma patients in the Cancer Genome Atlas database dataset. A-D: Correlation of the immune/stromal score with histologic grade and pathologic stage; E: Hepatocellular carcinoma (HCC) cases were divided into two groups based on their immune scores, as indicated by the log-rank test; F: Similarly, HCC cases were divided into two groups based on their stromal scores, as indicated by the log-rank test. TCGA: The Cancer Genome Atlas database; HCC: Hepatocellular carcinoma.

Validating the risk score model in GEO and ICGC dataset

To verify the robustness of our findings, this risk score model was further evaluated using the GEO and ICGC dataset, which included 222 and 233 HCC patients. The patients from the GEO and ICGC dataset were divided into high- and low-risk groups based on the previous formula. In agreement with the previous results, patients in the high-risk group had significantly worse OS than the low-risk group (Figure 4E and I).Figure 4F and J show the distribution of the risk score and gene expression data from the GEO and ICGC HCC cohort. Furthermore, the risk score model yielded an AUC of 0.661 at 1 year, 0.663 at 3 years and 0.676 at 5 years based on the GEO HCC cohort data (Figure 4G), and an AUC of 0.753 at 1 years, 0.65 at 3 years, and 0.715 at 5 years based on the ICGC HCC cohort data (Figure 4K). Recently, Long et al[36]constructed an immune prognostic model (IPM) including two genes (Triggering Receptor Expressed On Monocytes 1 and Exonuclease 1) to assess the prognosis of HCC patients. We calculated the C-indexes to compare the prognostic values of their model and our risk score model. As shown in Figure 4D, 4H and 4L, the concordance index of the risk score model for the first-, third-, and fifth-year OS was higher than the IPM in the TCGA, GEO and ICGC HCC cohorts, indicating that our risk score model had a better performance in evaluating prognosis. Above all, the risk score prognosis model is robust and efficient.

The difference of immune infiltration between the low- and high-risk HCC patients in the TCGA database

Figure 3 Comparison of gene expression profile with immune scores and stromal scores of hepatocellular carcinoma in the Cancer Genome Atlas database dataset. A: Venn diagrams showing the number of commonly upregulated or downregulated differentially expressed genes (DEGs) in stromal and immune score groups; B: Kyoto encyclopedia of genes and genomes (KEGG) analysis of DEGs, top 10 GO terms were displayed. False discovery rate of KEGG analysis was acquired from the database for annotation, visualization and integrated discovery functional annotation tool. DEGs: Differentially expressed genes; TNF: Tumor necrosis factor; KEGG: Kyoto encyclopedia of genes and genomes; DAVID: Database for annotation, visualization and integrated discovery.

We estimated the difference of immune infiltration between low- and high-risk HCC patients in 22 subpopulations of immune cells using the CIBERSORT algorithm. The proportion of immune cells in HCC varied significantly between the high- and lowrisk groups (Figure 5A and B). Figure 5C shows that a high fraction of M0 macrophages, T regulatory cells (Tregs) and T follicular helper cells (Tfh) mainly infiltrated in high-risk group patients. In contrast, a high fraction of CD8 T cells,resting CD4 memory T cells and monocytes mainly infiltrated in low-risk HCC patients. In addition, the proportions of different tumor-infiltrating immune cells(TIICs) showed a weak to moderate correlation (Figure 5D). Therefore, these results indicated that the different immune infiltration in patients with HCC could be used as a prognostic indicator and targets for immunotherapy.

Immunotherapy is increasingly applied to the clinical management of multiple tumors. However, only a part of patients with cancer showed a response to immunotherapy, and the efficacy of immunotherapy can be improved by identifying the type of immune infiltration. Hence, we investigated the relationship between patient risk scores and the expression of commonly immune checkpoints, and the results showed that the risk score was significantly associated with the expression of cytotoxic T-Lymphocyte associated protein 4 (CTLA-4), programmed cell death 1 (PD-1), and T-cell immunoglobulin mucin receptor 3 (TIM-3) (P ＜ 0.05) (Figure 6A, Supplementary Table 9). Furthermore, we explored the expression of the immune checkpoints in the high- and low-risk HCC patients. High-risk HCC patients had significantly higher expression of CTLA-4, PD-1 and TIM-3 than the low-risk HCC patients (P ＜ 0.05) (Figure 6B-D). Additionally, there was no significant difference in the expression of T Cell Immunoreceptor With Ig And ITIM Domains and Lymphocyte-activation gene 3 between the high- and low-risk HCC group (Supplementary Figure 4A and B). This suggests that the immunosuppressive microenvironment in high-risk patients may be responsible for their poor prognosis.

Alteration in immune-related pathways between high- and low-risk group patients

To investigate the difference of immune genes between high- and low-risk patients,we compared the gene expression profiles between the two groups. A total of 332 immune-related genes were extracted from the Molecular Signatures Database v6.2[20][(Immune system process M13664, Immune response M19817); http://www.broadinstitute.org/gsea/msigdb/index.jsp]. The results showed that there were 193 immunerelated genes in the DEGs between high- and low-risk groups(Supplementary Figure 4C). Furthermore, the GO and KEGG pathway analysis was performed to identify the potential functions of the 193 immune-related genes that were differentially expressed (Figure 6E and F, Supplementary Figure 4D and E)(Supplementary Table 10). The top 5 GO terms identified included immune response,inflammatory response, innate immune response, and cell-cell signaling (Figure 6E).The KEGG pathway analysis showed that these genes were enriched in the T cell receptor signaling pathway, cytokine-cytokine receptor interaction, rheumatoid arthritis, chemokine signaling pathway and tuberculosis (Figure 6F).

Figure 4 Prognostic analysis of the risk score model. A-C: Kaplan-Meier survival, risk score and time-dependent receiver operating characteristic (ROC) curves of the risk score model for the Cancer Genome Atlas database (TCGA) hepatocellular carcinoma (HCC) cohort; E-G: Kaplan-Meier survival, risk score and timedependent ROC curves of the risk score model for the Gene Expression Omnibus databases (GEO) HCC cohort; I-K: Kaplan-Meier survival, risk score and timedependent ROC curves of the risk score model for the International Cancer Genome Consortium database (ICGC) HCC cohort. A, E and I: OS was significantly higher in the low-risk score group than in the high-risk score group; B, F and J: Relationship between the risk score (upper) and the expression of eight prognostic immune genes (lower) is shown; C, G and K: Time-dependent ROC curve analysis of the risk score model; D, H and L: The concordance index (C-index) was used to evaluate prognostic performance for survival prediction. Performance was compared between the risk score model and immune prognostic model by calculating the C-index in the TCGA, GEO and ICGC HCC cohorts. TCGA: The Cancer Genome Atlas database; GEO: Gene Expression Omnibus databases; ICGC: International Cancer Genome Consortium database; IPM: Immune prognostic model; HCC: Hepatocellular carcinoma.

The risk score model is independent of conventional clinical characteristics

Figure 5 Landscape of immune infiltration in high- and low-risk hepatocellular carcinoma patients in the Cancer Genome Atlas database dataset. A:Relative proportion of immune infiltration in high- and low-risk patients; B: Heat map of the 22 immune cell proportions in high- and low-risk patients; C: Violin plots visualizing significantly different immune cells between high-risk and low-risk patients; D: Correlation matrix of all 22 immune cell proportions.

As shown in Figure 7, the risk score was significantly correlated with patient age,histologic grade, pathologic stage and vascular invasion (Figure 7A-Figure D, P ＜0.05). To explore the independent prediction of the risk score model, univariate and multivariate analyses were performed, the results showed that risk score and pathologic stage were independent prognostic indicators (Figure 7E). Taken together,these results indicated that the risk score model was an independent prognostic factor for OS.

Developing and validating a nomogram based on the risk score model

To establish a clinically applicable method to assess the prognosis of HCC patients,we developed a nomogram that included risk score and pathologic stage (Figure 7F).The concordance index of the nomogram was 0.71. The calibration plot for the possibility of 1-, 3- and 5-year survival showed good agreement between the prediction by risk score and actual observations (Figure 7G-I). The AUC was 0.77,0.799 and 0.773 for the 1, 3 and 5-year survival times, respectively. These results demonstrate that the nomogram might be a better model for predicting OS, and aid clinical management. Schematic diagram of the main altered pathway in the high- and low-risk patients is shown in Figure 8.

DISCUSSION

HCC remains a major challenge for public health worldwide[1,37]. Despite multiple therapeutic methods, such as surgical resection, liver transplantation, radiofrequency ablation and chemotherapy, the efficacy is limited. Moreover, effective prognostic indicators that can be utilized to guide cancer therapy is still lacking. Previous studies have shown that the TME is associated with tumor progression and patient prognosis[38-40]. Therefore, it is important to investigate the TME to identify biomarkers that can predict HCC patients' OS.

Figure 6 Enrichment analysis of the immune prognostic model. A: Correlation of the risk score with the expression of several prominent immune checkpoints; BD: Violin plots visualizing significantly different immune checkpoints between high-risk and low-risk patients; E-F: The Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis of immune related genes, top 10 GO terms were displayed. False Discovery Rate of GO and KEGG analysis was acquired from the DAVID functional annotation tool. CTLA-4: Cytotoxic T-Lymphocyte associated protein 4; PD-1: Programmed cell death 1; TIM-3: T-cell immunoglobulin mucin receptor 3; KEGG: Kyoto encyclopedia of genes and genomes; FDR: False discovery rate.

Figure 7 Relationship between the risk score model and other clinical information. A: Correlation of the risk score with the age; B: Correlation of the risk score with the histologic grade; C: Correlation of the risk score with the pathologic stage; D: Correlation of the risk score with the vascular invasion; E: Univariate and multivariate regression analysis of the relation between the risk score prognostic model and clinicopathological features regarding prognostic value; F: Nomogram for predicting the probability of 1-, 3-, and 5-year overall survival (OS) for hepatocellular carcinoma patients; G-I: Calibration plot of the nomogram for predicting the probability of OS at 1, 3, and 5 years; J-L: Time-dependent receiver operating characteristic curve analyses of the nomogram. ROC: Receiver operating characteristic curve; CI: Confidence interval.

Figure 8 Schematic diagram of the main altered pathway in the high- and low-risk patients. TME: tumor microenvironment; DAB2: Disabled homolog 2; MSC: Musculin; CXCL8: C-X-C Motif chemokine ligand 8; LGALS3:Galectin 3; BATF: B-Cell-activating transcription factor; KLRB1: Killer cell lectin like receptor B1; ENG: Endoglin;APCS: Adenomatosis polyposis coli tumor suppressor.

The ESTIMATE algorithm has been applied to multiple cancers, showing the effectiveness of the algorithms for big data analysis[41-44]. To investigate the TME of HCC, the ESTIMATE algorithm was applied to attain immune and stromal scores,where the results showed that high immune scores and stromal scores were significantly associated with worse OS, indicating that the TME was related to the prognosis of HCC patients. Subsequently, we divided HCC patients into high- and low-immune/stromal score groups to screen prognostic genes. We then analyzed the intersection DEGs yielded from the comparison of high vs low immune/stromal scores groups. GO analysis suggested that the DEGs were mainly involved in the TME, such as immune response, inflammatory response, cell adhesion, chemotaxis and extracellular matrix organization. The results were constitent with previous reports that immune cells and ECM molecules in the TME were interrelated[45-48],hence reshaping the immune microenvironment may improve the effect of cancer treatment[49-51]. The DEGs were identified, and univariate, LASSO and multivariate analyses were performed to construct a risk score model that could identify HCC patients with unfavorable outcomes. The AUC for the risk score model for predicting the 1-, 3-, and 5-year survival were 0.778, 0.754 and 0.75, respectively. Also, the risk score model was validated in the GEO and ICGC dataset. The results suggested that eight immune-related genes had a good performance for survival prediction,indicating that shaping the TME may inhibit or eliminate tumors, resulting in a favorable prognosis. The genes (DAB2, MSC, CXCL8, LGALS3, BATF, KLRB1, ENG and APCS) that compose our risk score model could be considered to be potential targets, and they may provide better performance in combination.

DAB2 is an adapter protein with signaling roles in the domain of endocytosis, cell differentiation, cell adhesion, angiogenesis, and homeostasis[52,53]. Previous studies have reported that DAB2 is a negative regulator of immune function[54]. DAB2 is expressed by macrophages, and it functions as a negative immune regulator of TLR4 endocytosis and signaling, controlling the inflammatory response to endotoxins[55]. In addition, a recent study revealed that low expression of DAB2 by dendritic cell (DC)enhanced IL-12 and IL-6 expression, besides improving the ability of DCs for antigen uptake, migration and T cell stimulation. DAB2-silenced DCs inhibited tumor growth[53]. MSC, a member of the basic helix-loop-helix transcription factors, enforced Foxp3 expression and promoted the unidirectional development of induced Treg cells(iTreg cells) by repressing the Th2 developmental program[56,57]. However, the role of MSC in HCC remains unclear. CXCL8 is a chemotactic factor secreted by malignant cells and various immune cells of multiple cancer types, which promotes tumor progression, recurrence and metastasis through shaping pro-tumoral vascularization,inflammation and immunity[58,59]. Accumulating studies demonstrate that CXCL8 is a prognostic marker and is a potential therapeutic target for HCC[60-62]. LGALS3, a glycan-binding protein secreted by cancer cells, has been regarded as an important regulator of multiple functions critical to cancer biology[63]. LGALS3, which is an important biomarker in prostate cancer[64,65], plays important roles in the progression and metastasis of colon cancer, acute myeloid leukemia, melanoma and pituitary tumors[66-69], and correlates with the infiltration of M2 macrophages[70]; BATF contains basic leucine zipper domains, is a member of the AP-1/ATF superfamily of transcription factors, plays a role in the growth and expansion of interleukin 17 (IL-17)-producing helper T cell (Th17) and iNKT cells expressing IL-17[71-73], the differentiation of Th17, Tfh and CD8+T cells[74], and controls the tumor formation and the progression of colorectal cancer[75]. KLRB1 (CD161) is a C-type lectin receptor expressed by most NK cells and subsets of T cells[76]. Previous studies indicated that KLRB1expressed by CD4+or CD8+T cells was associated with favorable outcomes in certain cancers such as lung cancer[77,78]. Recent studies demonstrated that the expression of LLT1 (ligand of KLRB1) by tumor cells may facilitate their escape from the immune system. Hence, the KLRB1/LLT1 receptor/ligand system appears to be a novel therapeutic target in the treatment of cancer[79,80]. ENG is a transforming growth factor beta 1 (TGFβ1) binding receptor. Teama et al[81]found that high expression of ENG/ TGFβ1 may contribute to carcinogenesis and the progression of HCC in cirrhotic patients. APCS (Amyloid P Component, Serum), a glycoprotein belonging to the pentraxin family of proteins, has a characteristic pentameric organization. Further functions of ENG and APCS remain unknown. In our study, increased expression of ENG and APCS correlated with a favorable outcome in HCC. To our knowledge, the eight gene signature related risk score prognostic model has not been previously reported and could be a novel prognostic factor of HCC.

In addition, multivariate Cox analyses demonstrated that the risk score and pathologic stage were independent prognostic indicators. Subsequently, a nomogram that included risk score and the pathologic stage was constructed. The calibration plot for the possibility of 1-, 3- and 5-year survival showed good agreement between the prediction by risk score and actual observation. The main feature of the nomogram is that it affords a personalized scoring system for patients and is feasible to predict prognosis. During the progression of tumor, the immune system plays a dual role in the complicated interactions between tumor and host; it conveys protective immunity by recognizing tumor-specific antigens to eliminate tumor cells, but also benefits tumor progression, either by altering tumor immunogenicity or by creating a microenvironment that can promote tumor outgrowth or aid in a subsequent metastatic cascade[82-85]. Therefore, tumor cell escape can occur through various immunosuppressive mechanisms, such as recruiting immunosuppressive cells (e.g.,Treg cells), increasing the expression of inhibitory ligands such as Programmed death ligand 1 (PD-L1), and decreasing the expression of major histocompatibility complex class I molecules, which results in immune tolerance[86-88]. Moreover, a potential antitumor immune response can be unleashed by blocking the function of immunosuppressive cells and immunosuppressive mechanisms[89]. We explored the immune mechanisms and the component of TIICs subpopulation between patients in the low- and high-risk groups. The results have shown that a high fraction of M0 macrophages, Tregs and Tfh were found in patients in the high-risk group. In contrast, a high fraction of CD8 T cell, resting memory CD4 T cells and monocytes was mainly found in low-risk group patients. We also explored the expression of the immune checkpoints between the high- and low-risk HCC patients. The high-risk HCC patients had significantly higher expression of CTLA-4, PD-1 and TIM-3 than the low-risk HCC patients (P ＜ 0.05). Previous studies revealed that resting memory CD4 T cells can be further differentiated and endowed with multiple functions,including restoration of immune tolerance to self-antigen or alloantigen and the promotion of CD8+ T cells to anti-tumor[90,91]. Tregs, which expressed CTLA-4, plays a vital role in the inhibition of anti-tumor immune responses. Treatment with an anti-CTLA-4 antibody has emerged as an effective therapy for the treatment of cancer[92-96].In our model, high-risk HCC patients had higher fraction of Tregs and a higher expression of CTLA-4, and a worse prognosis, indicating that the immunosuppressive environment and high expression of immune checkpoints may be responsible for the poor prognosis. Furthermore, these results suggest that treatment with antibodies against immune checkpoints will benefit the high-risk HCC patients more than the low-risk patients, hence resulting in a better prognosis.

In this study, we constructed a novel risk score model for prognostic prediction of HCC. One of the advantages of our risk score model is that it has high sensitivity and specificity in predicting the OS, being further validated using external databases. In addition, the risk score model is associated with the immunosuppressive environment and immune checkpoint expression, thus assisting clinicians in selecting personalized immunotherapy for HCC patients.

However, there are several limitations in our study. Firstly, the risk score model needs to be further validated in multicenter clinical trials and prospective studies.Secondly, the functional and mechanistic studies of the eight immune-related genes should be further carried out.

In summary, our research established and validated a risk score model that is based on eight immune-related genes to predict the OS of HCC, which may help in clinical decision making for individualized treatment. Notably, the risk score model provides an immunological viewpoint to clarify the mechanisms that determine the clinical outcome of HCC.

ARTICLE HIGHLIGHTS

Research background

Hepatocellular carcinoma (HCC) is a common malignant tumor with a poor prognosis. In recent years, immunotherapy has emerged as a novel and effective therapy and is being applied in various tumors including HCC. However, the influence of genes involved in the tumor microenvironment on the prognosis of HCC patients remains unclear. And the high-throughput studies that investigated the potential prognostic role of immune prognostic models in HCC are still lacking.

Research motivation

So far, only a small number of HCC patients receiving immunotherapy treatment exhibited responses due to the immunosuppressive microenvironment. Hence, it is necessary to

investigate the HCC microenvironment to identify prognostic genes that enable us to predict the benefit of immunotherapy, which may help in clinical decision making for individualized treatment.

Research objectives

To identify a robust gene signature associated with the HCC microenvironment to improve prognosis prediction and effectiveness of immunotherapy of HCC, we analyzed the data from The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO) and International Cancer Genome Consortium (ICGC) databases.

Research methods

We computed the immune/stromal scores of HCC patients obtained from TCGA based on the ESTIMATE algorithm. Univariate analysis, multivariate analysis and the least absolute shrinkage and selection operator, were utilized to construct our predictive model. This model was performed based on the significant differentially expressed genes screened established based on mRNA expression profiles from the TCGA database. The robustness of this model was validated using GEO and ICGC datasets.

Research results

The risk score model consisting of eight genes (Disabled homolog 2, Musculin, C-X-C motif chemokine ligand 8, Galectin 3, B-cell-activating transcription factor, Killer cell lectin like receptor B1, Endoglin, and Adenomatosis polyposis coli tumor suppressor) was constructed and validated based on HCC patients who were divided into high- or low-risk group. The receiver operating characteristic curve analysis confirmd the good potency of the risk score prognostic model. Moreover, we investigated the relationship between patient risk scores and the expression of common immune checkpoints, and the results showed that the risk score was significantly associated with the expression of Cytotoxic T-Lymphocyte associated protein 4,Programmed cell death 1, and T-cell immunoglobulin mucin receptor 3. To establish a clinically applicable method to assess the prognosis of HCC patients, a nomogram involving risk score and the pathologic stage was formulated.

Research conclusions

Our research established and validated a risk score model that is based on eight immune-related genes to predict the overall survival of HCC, which may help in clinical decision making for individualized treatment. The risk score model and the nomogram will benefit HCC patients through personalized immunotherapy.

Research perspectives

The risk score model provides an immunological viewpoint to clarify the mechanisms that determine the clinical outcome of HCC. Identifying effective molecular biomarkers and predictive markers of immunotherapy is a future direction for improving the effectiveness of immunotherapy.