A novel immune-related gene signature predicts the prognosis of hepatocellular carcinoma.

Immune-related genes play a key role in regulating the cancer immune microenvironment, influencing the overall survival of patients with hepatocellular carcinoma (HCC). Along with the rapid development of immunotherapy, identifying immune-related genes with prognostic value in HCC has attracted increasing attention. Here, we aimed to develop a prognostic signature based on immune-related genes. By investigating the transcriptome landscape of 374 HCC and 160 non-HCC samples in silico, a total of 2251 differentially expressed genes were identified. Among which, 183 differentially expressed immune-related genes were subjected to a univariate Cox proportional hazard model to screen for genes with possible prognostic significance. A 10-gene prognostic signature, including HLA-G, S100A9, S100A10, DCK, CCL14, NRAS, EPO, IL1RN, GHR and RHOA, was generated employing a multivariate Cox proportional hazard model. Kaplan-Meier and Receiver Operator Characteristic (ROC) curves were used to evaluate the prognostic utility of the 10-gene signature. Moreover, the underlying mechanisms of these genes were analyzed via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment. According to the Tumor Immune Estimation Resource (TIMER) database, our prognostic signature was significantly associated with tumor-infiltrating B cells, CD4 T cells, dendritic cells, macrophages and neutrophils. Our study provides a novel prognostic signature based on immune-related genes associated with clinical outco mes of HCC.

Introduction

Hepatocellular carcinoma (HCC) is a highly heterogeneous and lethal malignancy which ranks as the second leading cause of cancer-related death worldwide (1, 2). Despite considerable efforts to improve clinical outcomes, the prognosis remains unsatisfactory. The Barcelona Clinic Liver Cancer (BCLC) remains the most widely used evaluation system for predicting HCC in clinical practice (3). It focuses primarily on a subset of patients who may benefit from ablative or intra-arterial treatment, which may not meet the predictable need for innovative treatment (4). Therefore, there is still a need for an optimized prediction system to guide clinical treatment.

Over the past decade, immunotherapy has become one of the most promising cancer treatments (5). Immune checkpoint inhibitors, which are essential components of immunotherapy, play a significant role in the treatment of advanced HCC (6). Programmed cell death protein 1 (PD-1) monoclonal antibodies, such as nivolumab and pembrolizumab, have already been approved by the Food and Drug Administration (FDA) for the treatment of advanced HCC (7). Nivolumab and pembrolizumab treatment resulted in a 2-year survival rate of 80% for patients with advanced HCC. Additionally, cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) blockade with tremelimumab induced sustained objective remission in patients with HCC who were also diagnosed with hepatitis C virus (HCV) (8). However, poor clinical responses to immunotherapy are still common due to limited infiltration of immune checkpoint inhibitors and the impaired anti-tumor effect of CD8+T cells in the tumor microenvironment (9). Immune-related genes that participate in immunotherapy resistance and tumor progression (10) play a central role in modulating the tumor immune microenvironment (11, 12). Therefore, developing a prognostic monitoring system based on immune-related genes would provide a precise prediction model for patients and insight into new therapeutic approaches.

Materials and methods

Dataset acquisition and pre-processing

The Level-3 gene expression data and the corresponding clinical characteristics of HCC patients were obtained from The Cancer Genome Atlas (TCGA) portal (https://portal.gdc.cancer.gov/). Additionally, a normal liver transcriptome dataset derived from the Genotype-Tissue Expression Project (GTEx) was obtained from the University of California, Santa Cruz (UCSC) Xena database (http://xena.ucsc.edu/). All raw data were annotated according to the human genome annotation file (GRCh38.p12) obtained from Ensembl (http://asia.ensembl.org/Homo_sapiens/Info/Index). Following normalization using the same criteria, log2(FPKM+1), the raw datasets from TCGA and GTEx were merged to form a study project with 374 HCC samples and 160 non-tumor samples. A subset of immune-related genes was downloaded from the ImmPort database (https://www.immport.org/shared/home). The characteristics of tumor-infiltrating immune cells were retrieved from the Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) database, which was compiled based on the TCGA data.

Differentially expressed gene analysis

Differentially expressed genes (DEGs) between HCC tumor and non-tumor samples were obtained using the EdgeR package (13) in R software (version 3.6.1) (14). Genes were considered DEGs if they met the following screening criteria: false discovery rate (FDR) < 0.05 and a log2 |fold change| > 1. Differentially expressed immune-related genes (DEIGs) were extracted from DEGs by intersecting them with the immune-related gene dataset.

Development of a prognostic model and analysis

The primary endpoint of the study was death. Patients with a survival time of less than 30 days or those lost to follow-up were censored (15). A total of 342 patients with HCC were subjected to univariate Cox analysis to identify genes associated with survival. The identified DEIGs were considered hub genes. Statistical analyses were conducted via the survival package (16) in R. A value of p < 0.01 was considered statistically significant.

Hub genes were subsequently subjected to a multivariate Cox regression analysis to construct a prognostic model. A receiver operator characteristic (ROC) curve was developed to assess the performance of the prognostic model. The prognostic model and clinical parameters (age, gender, the grade and stage of disease) were subjected to univariate and multivariate analyses to identify the independent prognostic factors. Genes that were eligible to construct the prognostic model were identified as key genes.

Survival analysis

All patients with HCC were allocated with a risk score according to the prognostic risk model formula, then the median value of the risk scores was employed as the cutoff. Patients with a risk score higher than the median value were considered the high-risk group and patients with a risk score lower than the median value were considered the low-risk group. We then compared the mortality rates of the two groups by R survival (16) and limma packages (17). Univariate survival analysis was performed for each of the key genes within the prognostic model. The median expression level was used to divide the patients into high- and low-expression groups.

Gene enrichment analysis

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were conducted to explore the possible molecular mechanisms of the DEIGs and the biological functions of the 10 key genes in the signature.

Correlation analysis

The relationship between the prognostic model and the infiltration level of immune cells in HCC samples was analyzed and visualized via R. The p-value was set to p < 0.05 as significance level and the absolute value of the correlation coefficient was set to < 0.3 as the cutoff values.

Results

Identification of DEIGs

A flow diagram of the present study is presented in . The tumor group consisted of 374 HCC samples obtained from TCGA, and the non-tumor group included 50 paracarcinoma samples from TCGA and 110 normal liver samples from GTEx. We identified 2251 DEGs by comparing the expression profiles of the tumor and non-tumor groups. There were 880 upregulated and 1371 downregulated DEGs observed in HCC. By intersecting immune-related genes with DEGs, 183 DEIGs were identified, of which 75 were upregulated in HCC ( ).

Figure 1

Flow diagram of the study.

Figure 2

DEIGs between 374 HCC and 160 non-tumor samples. DEIGs between tumor and non-tumor groups are depicted in the heatmap (A) and volcano plot (B). Among the 183 identified DEIGs, 75 genes were upregulated and 108 genes were downregulated in HCC. GO enrichment analysis revealed 180 enriched DEIGs (C, left). The top five GO items with the most enriched genes are listed (C, right). KEGG enrichment analysis revealed 130 enriched DEIGs (D, left) and the 10 KEGG items with the most enriched genes are listed (D, right). DEIGs, differentially expressed immune-related genes; HCC, hepatocellular carcinoma; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; N, non-tumor group; T, tumor group; BP, biological processes; CC, cellular components; MF, molecular functions.

Subsequent enrichment analyses indicated that 180 and 130 DEIGs were enriched by GO ( ) and KEGG ( ) analyses, respectively. The inflammatory related items were most frequently mentioned in the enrichment analyses. In the GO enrichment analysis, “leukocyte migration”, “side of membrane” and “receptor-ligand activity” were the most frequent biological terms among “biological processes (BP)”, “cellular components (CC)” and “molecular functions (MF)” categories, respectively ( ). In the KEGG enrichment analysis, “Epstein-Barr virus infection” was the most predominant item that significantly enriched 29 DEIGs ( ).

Construction of the prognostic model

Twenty-nine immune-related genes were identified as hub genes with clinical significance in HCC according to the univariate Cox proportional hazard model. Among these genes, 12 genes with a hazard ratio (HR) less than 1.0 were regarded as protective factors, whereas 17 genes with an HR greater than 1.0 were considered as risk factors ( ). We compared the expression disparity of the 29 hub genes between tumor and non-tumor tissues. The results showed that 13 hub genes were upregulated and 16 were downregulated in HCC ( ).

Figure 3

Immune-related hub genes in HCC. Twenty-nine DEIGs were identified that were associated with clinical outcomes in HCC (A, left). Twelve had an HR less than 1.0 and were considered protective factors, whereas 17 genes had an HR greater than 1.0 (A, right) and were regarded as risk factors. The heatmap delineating the expression profiles of the 29 hub genes (B). The volcano plot showing that 13 hub genes were upregulated and 16 genes were downregulated in HCC (C). DEIGs, differentially expressed immune-related genes; HCC, hepatocellular carcinoma; HR, hazard ratio; N, non-tumor group; T, tumor group.

Subsequently, the hub genes were subjected to multivariate Cox analysis to construct a prognostic risk score model. Ten hub genes were identified as key genes for developing the prognostic model. According to the model formula, each patient was labelled with a risk score. The risk model is as follows:

Risk score = [Expression level of HLA-G* (-0.38009)] + [Expression level of S100A9* 0.23467] + [Expression level of S100A10* 0.14180] + [Expression level of DCK* 0.81204] + [Expression level of CCL14* 0.52696] + [Expression level of NRAS* 0.59515] + [Expression level of EPO*0.24953] + [Expression level of IL1RN* (-0.18841)] + [Expression level of GHR* (-0.30230)] + [Expression level of RHOA* (-0.43270)].

Verification of the prognostic model

To investigate the prognostic value of the model, we divided the 342 enrolled patients with HCC into high- and low-risk groups based on the median value of the risk scores according to the prognostic risk model ( ). A higher risk score indicated a shorter survival time ( ). Additionally, patients in the low-risk group had significantly longer survival times than those in the high-risk group ( ). The five-year survival rate of the high- and low-risk groups were 33.6% and 62.4%, respectively ( ). In addition, we assessed the predictive ability of the prognostic model using the ROC curve. The area under the curve (AUC) was 0.808 ( ), indicating that the prediction model was highly accurate.

Figure 4

Evaluation of the prognostic model. (A) Distribution of 342 patients with HCC according to the prognostic model. (B) Survival status based on the risk score of the prognostic model. Patients with a low-risk score exhibited a longer survival time. (C) Based on the survival curves, the OS was significantly lower in the high-risk than in the low-risk groups. (D) The survival rate of high- and low-risk groups is listed. (E) According to the ROC curve, the clinical utility of the model is validated (AUC=0.808). Both univariate (F) and multivariate (G) Cox analyses revealed that the prognostic model was an independent predictive indicator for HCC. OS, overall survival; ROC, receiver operator characteristic; AUC, area under curve. TNM, Tumor-Node-Metastasis staging.

To determine whether the model is an independent predictor of prognosis, we examined the prognostic model along with commonly used clinical parameters, including age, sex, grade and cancer stage of the patients. A prognostic analysis was conducted on 221 patients with complete clinical data via Cox regression analysis. In conjunction with both univariate ( ) and multivariate ( ) Cox analyses, the prognostic model was demonstrated to be a reliable and independent indicator of HCC clinical outcomes.

Estimation of tumor infiltrating immune cells in HCC samples

The mechanisms and pathways involved in the prognostic model were depicted by performing 10 key gene-based enrichment analyses. GO enrichment analysis indicated that no genes were enriched in the GO-CC item. The most frequently occurring items in GO-BP and GO-MF were “leukocyte cell-cell adhesion” and “cytokine activity”, respectively ( ). In the GO-term enrichment analysis, RHOA, EPO and HLA-G were the most widespread genes, found in 119, 69 and 59 terms, respectively. In KEGG analysis, the pathway “cytokine-cytokine receptor interaction” was the most significant pathway ( ) and RHOA and NRAS were the most abundant key genes.

Figure 5

Enrichment analysis of the 10 key genes. Bubble plot (A, left) and list box (A, right) show the top 10 GO items for BP and MF. Significant KEGG pathways of the 10 key genes visualized by bubble plot (B, left) and list box (B, right). GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological processes, MF, molecular functions.

We speculated that the 10 key genes contribute to modulating the tumor immune microenvironment. Therefore, we examined the correlation between the prognostic model and the infiltration level of the six immune cells. Our results suggested a positive correlation between the prognostic model and B cells, CD4 T cells, dendritic cells, macrophages and neutrophils ( ) in HCC. There was no correlation between the prognostic model and the presence of tumor infiltrating CD8 T cells ( ).

Figure 6

Relationship between the prognostic model and the infiltration level of six immune cell types. B cells (A), CD4 T cells (B), dendritic cells (D), macrophages (E) and neutrophils (F) were significantly correlated with the prognostic model, but not CD8 T cells (C).

Signatures of the 10 key genes

Univariable survival analysis was carried out for all genes involved in the prognostic model to gain further insight into the components of the predictive model. Patients were divided into high- and low-expression groups based on the median values of normalized expression levels. Except for IL1RN, nine other genes were significantly associated with the overall survival time of patients with HCC ( ). In accordance with the results of the Cox analysis ( ), patients with high expression levels of CCL14, HLA-G and GHR had longer survival times ( ), whereas those with high expression levels of S100A9, S100A10, DCK, NRAS, EPO and RHOA had shorter survival times ( ).

Figure 7

Survival curves of the 10 key genes. Patients with high expression of CCL14 (A), HLA-G (B) and GHR (C) survived longer than those with low expression of these genes. High expression of S100A9 (D), S100A10 (E), DCK (F), NRAS (G), EPO (H) and RHOA (I) was associated with poor clinical outcomes.

Discussion

HCC is a highly heterogeneous disease with poor clinical outcomes. With the advent of immunotherapy in malignant diseases, immune-related genes have been of interest to the research community in recent years. In this study, we constructed an immune-related 10-gene model to predict HCC prognosis by mining public datasets from TCGA and GTEx.

We collected transcriptome data from 374 HCC and 50 para-tumor samples from TCGA. To reduce the bias in the analysis caused by the disparity in sample sizes, expression profiles of 110 normal liver samples from GTEx were added to the non-tumor group. As HCC often occurs in the context of cirrhosis and hepatitis (18), the transcriptome profiles of para-tumors may differ significantly from those of normal liver tissues. Adding normal liver samples to the non-tumor group provides a better analysis strategy with improved validity and reliability.

Following the acquisition of DEIGs by intersecting DEGs with immune-related genes, univariate Cox analysis was conducted to identify immune-related genes that were associated with clinical prognosis. A 10-gene prognostic model was developed and evaluated using univariate and multivariate Cox proportional hazard models, Kaplan–Meier estimates and ROC curves. According to the evaluation, the 10-gene model was capable of reliably predicting HCC outcomes.

The prognostic model was constructed based on expression signature of the10 genes, including CCL14, HLA-G, GHR, S100A9, S100A10, DCK, NRAS, EPO, RHOA and IL1RN, which was quite different from the previous studies that focused on DEGs (19, 20). Patients with HCC with high CCL14, HLA-G and GHR expression had significantly better outcomes. This is in accordance with other studies demonstrating that CCL14 (1, 21, 22), HLA-G (23, 24) and GHR (25) are potential tumor suppressors in HCC. High expression levels of S100A9 (26, 27), S100A10 (28), DCK (29), NRAS (30, 31), EPO (32, 33) and RHOA (34, 35) were associated with poor prognosis. The relationship between IL1RN and clinical outcomes is unclear and little is known about its biological function in HCC.

Furthermore, we validated the close correlation between the prognostic model and the level of immune-infiltrating cells. According to the TIMER database (36), B cells, CD4 T cells, dendritic cells, macrophages, and neutrophils had significant positive correlations with the prognostic model. Although the immune context of HCC is far more complex, the relationship between the prognostic model and the above five types of immune cells provides insight into the regulation of the tumor immune microenvironment.

It is noteworthy that the 10 key genes identified in this study were predominantly involved in ligand-receptor mechanisms for regulating the tumor immune microenvironment. Specifically, CCL14 (C-C motif chemokine ligand 14) is a ligand of CCR1 and CCR5 (37). HLA-G is a ligand of ILT2, ILT4 and KIR2DL4 (38). GHR is a well-known receptor for growth hormones (39). EPO is a ligand of EPOR. IL1RN is a natural interleukin 1 receptor antagonist that inhibits IL-1 activity (40). The ligand-receptor pattern appears to be a promising target for cancer immunotherapy. In summary, the prognostic effect of the immune-related predictors offers novel insights into the tumor immune microenvironment, however, a comprehensive study is still required.

This study has certain limitations that should be addressed. The prognostic risk model must be validated in an independent clinical cohort. Furthermore, the sensitivity and specificity of our prognostic immune biomarkers should be examined both in vitro and in vivo. Finally, this study concentrated on transcriptome expression profiles, therefore, it cannot reflect the full landscape of immune-related genes in HCC.

In conclusion, using transcriptome profiles of immune-related genes, we developed a reliable model for predicting survival outcomes in patients with HCC. Additionally, as an independent prognostic factor, the model predicted the abundance of immune-infiltrating cells in HCC.

Data availability statement

The original contributions presented in the study are included in the article/, further inquiries can be directed to the corresponding author/s.

Author contributions

L.Y designed, guided the study; SJ.H analyzed data and drafted the manuscript; JQ.Q and L.W acquired the data. All authors contributed to the article and approved the submitted version.

Funding

National Natural Science Foundation of China (82030076, 82070161, 81970151, 81670162 and 81870134). Shenzhen Science and Technology Foundation (JCYJ20190808163601776 and JCYJ20200109113810154). Shenzhen Key Laboratory Foundation (ZDSYS20200811143757022). Sanming Project of Shenzhen (SZSM202111004).

Acknowledgments

The data supporting this publication is available at UCSC Xena (http://xena.ucsc.edu/), Ensembl (http://asia.ensembl.org/Homo_sapiens/Info/Index), ImmPort (https://www.immport.org) and TIMER (https://cistrome.shinyapps.io/timer/).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.