A novel immune signature to predict the prognosis of patients with hepatocellular carcinoma.

ABSTRACT: Aberrant immunity has been associated with the initiation and progression of cancers such as hepatocellular carcinoma (HCC). Here, we aim to develop a signature based on immune-related genes (IRGs) to predict the prognosis of HCC patients. The gene expression profiles of 891 HCC samples were derived from 4 publicly accessible datasets. A total of 1534 IRGs from Immunology Database and Analysis Portal website were obtained as candidate genes for prognostic assessment. Using least absolute shrinkage and selection operator (LASSO) regression analysis, 12 IRGs were selected as prognostic biomarkers and were then aggregated to generate an IRG score for each HCC sample. In the training dataset (n = 365), patients with high IRG scores showed a remarkably poorer overall survival than those with low IRG scores (log-rank P < .001). Similar results were documented in 3 independent testing datasets (n = 226, 221, 79, respectively). Multivariate Cox regression and stratified analyses indicated that the IRG score was an independent and robust signature to predict the overall survival in HCC patients. Patients with high IRG scores tended to be in advanced TNM stages, with increased risks of tumor recurrence and metastasis. More importantly, the IRG score was strongly associated with certain immune cell counts, gene expression of immune checkpoints, estimated immune score, and mutation of critical genes in HCC. In conclusion, the proposed IRG score can predict the prognosis and reflect the tumor immune microenvironment of HCC patients, which may facilitate the individualized treatment and provide potential immunotherapeutic targets.

Introduction

Hepatocellular carcinoma (HCC) is a dominant subtype of primary liver cancer and occurs commonly in patients with cirrhosis. This type of cancer has ranked the sixth most frequently-diagnosed malignancy and the second leading contributor to cancer-related deaths worldwide and, unfortunately, the disease burdens are stably increasing in Western populations.[ The therapeutic options for HCC include surgical resection, tumor ablation, liver transplantation, and multikinase inhibitors sorafenib.[ However, recurrence occurs highly in HCC patients even with treatment in early stage; and, for patients with advanced HCC, the survival rate remains unacceptably poor.[ In recent years, immune checkpoint inhibitors, such as those targeting programmed cell death protein 1 (PD-1) and cytotoxic T-lymphocyte protein 4 (CTLA-4), have emerged as potentially effective treatments for HCC patients in advanced stages.[ This emphasizes the functional importance of tumor immune microenvironment (TIM) in HCC.

TIM contains diverse of immune cell types that exhibit either immune promotive or suppressive roles, which has been considered to restrict the accumulation of cytotoxic T cells to the vicinity of cancer cells.[ In patients with HCC, the prognostic impacts have been demonstrated for tumor immune suppressors like regulatory T cells, tumor-associated macrophages, and myeloid-derived suppressor cells.[ To date, however, there has been no signature to predict the overall survival (OS) of HCC patients by systematically evaluating the TIM based on immune-related genes (IRGs). Chew et al[ previously developed an immune model that could determine the long-term survival in resectable HCC. Nevertheless, their study focused only on 14 IRGs with a limited number of included patients, which may not provide a comprehensive assessment on the TIM and patient's prognosis. Therefore, it is necessary to identify a signature that can reflect the status of TIM with prognostic capacity in HCC patients.

In this investigation, we aim to develop an immune signature to predict the OS of HCC patients based on the comprehensive list of IRGs in the Immunology Database and Analysis Portal database. The gene expression data of HCC samples profiled by RNA-sequencing or microarray from publicly accessible databases were used for analyses. We then established and validated the prognostic capacity of the IRG score, and examined its association with important clinicopathologic features. And last, the tumor immune-related characteristics involved by the IRG score were figured out. The expected results may enhance our understanding about the role of TIM in the development of HCC and promote the discovery of novel therapeutic targets.

Materials and methods

Acquisition of HCC datasets

In this study, the gene expression profiles and corresponding clinical information of HCC patients were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov), the International Cancer Genome Consortium (ICGC, https://icgc.org), and the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) databases. The RNA-sequencing data of 365 HCC samples in the TCGA-LIHC cohort from the TCGA database was used as the training set. The LIRI-JP cohort from the ICGC database was adopted as a testing set, which included 226 HCC patients and was profiled by RNA-sequencing. The microarray data of GSE14520 (n = 221) and GSE54236 (n = 79) from the GEO database were also selected as testing sets. The clinical features of the totaling 891 patients with HCC were summarized in Table 1. This study was a secondary analysis of publically accessible data, thus no ethical approval was required.

Table 1

Clinical data of patients with hepatocellular carcinoma in the 4 cohorts.

Variables	TCGA-LIHC (n = 365)	LIRI-JP (n = 226)	GSE14520 (n = 221)	GSE54236 (n = 79)
Age, yr
>60	192 (52.6)	177 (78.3)	40 (18.1)	–
≤60	173 (47.3)	49 (21.7)	181 (81.9)
Gender
Female	119 (32.6)	61 (27.0)	30 (13.6)	17 (21.5)
Male	246 (67.4)	165 (73.0)	191 (86.4)	62 (78.5)
Stage
I/II	254 (69.6)	138 (61.1)	170 (76.9)	–
III/IV	87 (23.8)	88 (38.9)	49 (22.2)
Unknown	24 (6.6)	0 (0.0)	2 (0.9)
AFP, U/L
>50	181 (49.6)	–	100 (45.2)	–
≤50	95 (26.0)		118 (53.4)
Unknown	89 (24.4)		3 (1.4)
Death
No	235 (64.4)	185 (81.9)	136 (61.5)	0 (0.0)
Yes	130 (35.6)	41 (18.1)	85 (38.5)	79 (100.0)

Data were presented as number (percentage).

– = not available, AFP = alpha fetoprotein.

Data preprocessing

The RNA-sequencing data were analyzed using fragments perkilobase of exon per million fragments mapped value with log2(x + 1) transformation. The Ensembel IDs in data expression matrix were matched to gene symbols using the GTF file from GENCODE (https://www.gencodegenes.org, version 23).[ For microarray data profiled by the Affymetrix platforms, the raw CEL files were downloaded and then processed with robust multichip average algorithm.[ For gene expression data measured by the Agilent microarray, we downloaded the normalized series matrix for subsequent analysis. The probe sets in the microarray data were then mapped to gene symbols according to the annotation information in corresponding platforms.

Derivation of the immune signature

The comprehensive list of IRGs was downloaded from the Immunology Database and Analysis Portal database (https://www.immport.org), which contained a total of 1534 genes (1024 non-duplicates). These genes were involved in 16 immune categories, such as antimicrobials, T-cell receptor signaling, B-cell receptor signaling, and chemokine.[ To construct the immune signature for prognosis prediction, the expression data of these IRGs were extracted from the training set, and those with zero expression values in >10% of samples were removed. Univariate Cox regression analyses were performed to determine the association between IRGs expression levels and patients’ OS, and those IRGs with Benjamini–Hochberg adjusted P-values < .05 were selected for further analyses. Then, least absolute shrinkage and selection operator (LASSO) penalized Cox regression model were utilized to determine the most significantly OS-related IRGs. In this model, the penalty parameter lambda.1se for prevention of overfitting was obtained using 10-fold cross validation.[ An IRG score for OS prediction was finally established by linear combination of the IRGs expression levels weighted by the coefficients from LASSO regression.

Statistic analyses

All statistical analyses were realized with R 3.5.0 software (The R Foundation for Statistical Computing, Vienna, Austria). To evaluate the prognostic ability of the IRG score, time-dependent receiver operating characteristic (ROC) curve were analyzed and the area under the curve (AUC) was calculated using the “survivalROC” package.[ For subsequent comparison, we separated patients into high- and low-risk groups according to the optimal cut-off IRG score derived from 3-year ROC curve. Kaplan–Meier survival analysis with log-rank test was used to examine the OS difference between the 2 groups. Multivariate Cox regression analysis was applied to test the independence of predictors for patients’ OS. Stratified analyses by TNM stage were carried out to evaluate prognostic ability of the IRG score in subpopulations. The association between IRG scores and important clinicopathologic features including TNM stage, metastasis risk, and tumor relapse were assessed by one-way ANOVA, Wilcoxon rank sum test, or log-rank test when appropriate.

Immune profiles and gene mutation

The immune infiltrations of HCC samples were estimated using the “MCPcounter” package in R, which robustly quantifies the absolute abundance of 6 immune cell types based on the transcriptomic data of bulk tumors.[ The overall immune scores of HCC samples were calculated using the “ESTIMATE” package in R, as proposed by Yoshihara et al.[ Pearson correlation test was then performed to explore the association of the IRG score with the estimated immune cell counts, gene expression of immune checkpoints, and immune scores. The gene-level copy number variation profile of TCGA-LIHC cohort was downloaded from the UCSC Xena database (https://xenabrowser.net, TCGA-hub). We only focused on the difference of commonly mutated genes in HCC,[ including TP53, CTNNB1, NCOR1, RB1, ERRFI1, CDKN2A, ALB, and APOB.

Results

Construction of IRG signature

The immune signature was constructed as mentioned above (see Figure S1, Supplemental Digital Content, which illustrates the process of data analysis). Univariate Cox regression analysis was conducted in 732 abundant IRGs from the TCGA-LIHC cohort, and we identified 59 genes with adjusted P-values < .05 for OS prediction. These genes were selected for LASSO penalized Cox regression analysis, and a total of 12 IRGs were screened out as the most informative predictors (Fig. 1A). A prognostic IRG score was then calculated by summarizing the expression levels of the 12 IRGs weighted by their coefficients from LASSO model (see Table S1, Supplemental Digital Content, which illustrates the results of univariate Cox regression and LASSO coefficients of the 12 IRG for OS prediction). In the training set, the AUC of the 3- and 5-year ROC curves achieved 0.70 and 0.72, respectively, suggesting a good performance of the IRG score in OS prediction of HCC patients (Fig. 1B). Patients in the training set were then assigned into the high- (n = 172) and low-risk (n = 193) groups according to the best cut-off IRG score (4.51) derived from the 3-year ROC curve. The distribution of OS status and expression patterns of the 12 IRGs were shown in Fig. 1C. In the high-risk group, as expected, 9 risky IRGs were upregulated and the remaining 3 protective IRGs tended to be downregulated. Patients in the high-risk group showed significantly poorer OS than those in the low-risk group (hazard ratio [HR]: 3.08, 95% confidence interval [CI]: 2.14–4.45, log-rank P < .001; Fig. 1D).

Figure 1

Construction of the immune signature in the training set. A, 10 cross-fold validation in the LASSO model. B, AUC of time-dependent ROC curves for predicting patients’ overall survival. C, Distribution of patients’ overall survival and the heatmap of the 12 IRGs expression. D, Kaplan–Meier curves for overall survival prediction by the immune signature. AUC = area under the curve, IRGs = immune-related genes, ROC = receiver operating characteristic.

For the 3 testing sets, the IRG scores were predicted by using the same formula developed in the training set, and patients were also divided into high- and low-risk groups according to the optimum cut-off value from the 3-year ROC curve. The AUC of the 3- and 5-year ROC curves were 0.69 and 0.78, respectively, for the LIRI-JP cohort and 0.65 and 0.67, respectively, for the GSE14520 dataset (Fig. 2A and B). For the GSE54236 dataset, the AUC of the 3-year ROC curve was 0.74 (Fig. 2C). Similar to the training cohort, there was a significant difference in OS between the high- and low-risk groups in the LIRI-JP cohort (HR: 5.44, 95% CI: 1.93–15.30, log-rank P < .001; Fig. 2D), the GSE14520 dataset (HR: 2.75, 95% CI: 1.73–4.40, log-rank P < .001; Fig. 2E), and the GSE54236 dataset (HR: 2.79, 95% CI: 1.74–4.47, log-rank P < .001; Fig. 2F).

Figure 2

Validation of the immune signature in the testing set. A–C, AUC of time-dependent ROC curves for predicting patients’ overall survival in the LIRI-JP cohort, the GSE14520 dataset, and the GSE54236 dataset. D–F, Kaplan–Meier curves for overall survival prediction in the LIRI-JP cohort, the GSE14520 dataset, and the GSE54236 dataset. AUC = area under the curve, ROC = receiver operating characteristic.

IRG score and clinicopathologic features

To explore the independence of the immune signature in predicting patients’ OS, we performed multivariate Cox regression by including age, sex, TNM stage, alpha fetoprotein, and IRG score as explanatory variables. The results demonstrated that in all of the datasets, the IRG score was an independent predictor of OS in patients with HCC (HR = 1.19, 1.12, 1.04, 1.02, respectively; Table 2). Besides, we found that TNM stage was also significantly associated with patients’ OS; thus, additional stratified analysis was performed to examine the prognostic ability of the IRG score in patients with different TNM stages. Due to the limited number of patients with stage III–IV, the GSE14520 and the GSE54326 datasets were not included in such analyses. As a result, the IRG score could predict the OS of HCC in patients with stages I–II or stages III–IV (see Figure S2, Supplemental Digital Content, which illustrates the results of stratified analysis by TNM stage).

Table 2

Multivariate Cox regression analysis in patients with hepatocellular carcinoma.

	TCGA-LIHC (n = 365)		LIRI-JP (n = 226)		GSE14520 (n = 221)		GSE54236 (n = 79)
Variables	HR (95% CI)	P-value	HR (95% CI)	P-value	HR (95% CI)	P-value	HR (95% CI)	P-value
Age∗	1.02 (1.00–1.04)	.051	1.00 (0.97–1.04)	.945	1.00 (0.98–1.02)	.751	–
Gender
Female	Reference		Reference		Reference		Reference
Male	0.78 (0.49–1.26)	.316	0.41 (0.21–0.79)	.008	1.21 (0.57–2.55)	.625	1.34 (0.75–2.39)	.318
Stage
I/II	Reference		Reference		Reference		–
III/IV	1.58 (0.96–2.61)	.074	2.32 (1.16–4.63)	.017	2.57 (1.47–4.20)	<.001
AFP
≤50 U/L	Reference		–		Reference		–
>50 U/L	0.75 (0.45–1.26)	.272			1.38 (0.88–2.14)	.157
IRG score∗	1.19 (1.12–1.27)	<.001	1.12 (1.04–1.21)	.002	1.04 (1.01–1.07)	.006	1.02 (1.01–1.03)	<.001

– = not available, AFP = alpha fetoprotein, CI = confidence interval, HR = hazard ratio.

Analyzed as continuous variable.

To further clarify the clinical implications of the immune signature, we investigated the association between IRG score and some important clinicopathologic features such as TNM stage, tumor metastasis, and tumor recurrence. In the TCGA-LIHC and the LIRI-JP cohorts, we identified a positive correlation between the IRG score and TNM stage (Fig. 3A). Also, in the GSE14520 dataset, patients with advanced TNM stage or increased metastasis risk had a higher IRG score (Fig. 3B). Besides, patients with high IRG score exhibited a shorter relapse-free survival in both the TCGA-LIHC cohort (HR: 1.50, 95% CI: 1.08–2.10, log-rank P = .016; Fig. 3C) and the GSE14520 dataset (HR: 1.91, 95% CI: 1.32–2.76, log-rank P < .001; Fig. 3D).

Figure 3

The associations of IRG score with important clinicopathologic features. A, The relationship between TNM stage and IRG score in the TCGA-LIHC and the LIRI-JP cohorts. B, The correlation of IRG score with TNM stage and metastasis risk in the GSE14520 dataset. C and D, Kaplan–Meier curves for relapse-free survival prediction in the TCGA-LIHC cohort and the GSE14520 dataset. IRGs = immune-related genes.

IRG score, immune profiles, and gene mutation

Tumor immune infiltration was estimated for the TCGA-LIHC cohort, the LIRI-JP cohort, and the GSE14520 dataset. As shown in Fig. 4A, the IRG score was significantly associated with the estimates of some immune cell types, such as overall T cells, CD8+ T cells, and monocytic lineage. The IRG score was also positively related to the gene expression of immune checkpoints like CD274, PDCD1, and CTLA4. Accordingly, in the 3 datasets, the IRG score was positively correlated with the immune score estimated by the R “ESTIMATE” package (Fig. 4B–D). To determine the mutations associated with the signature, we compared the IRG score between patients with and without mutation of common genes, including TP53, CTNNB1, NCOR1, RB1, ERRFI1, CDKN2A, ALB, and APOB. The results showed that for all of these genes (except ALB), the mutated-type group had a higher IRG score compared with the wild-type group (Fig. 4E).

Figure 4

The associations of IRG score with immune profiles and gene mutations. A, The correlation of IRG score with estimated immune cell counts and gene expressions of immune checkpoints. B–D, The relationship between IRG score and ESTIMATE immune score in the TCGA-LIHC cohort, the LIRI-JP cohort, and the GSE14520 dataset. E, The relationship between critical gene mutations and IRG score in the TCGA-LIHC cohort. IRGs = immune-related genes.

Discussion

The treatment opinions for HCC patients have experienced rapidly changes in recent years, especially the advances in immunotherapy. This highlights the critical role of TIM in the progression of HCC, and points towards the necessity to identify immune-related biomarkers for prediction of patients’ prognosis. In this study, we constructed a robust immune-related risk signature for HCC using the data from TCGA-LIHC cohort and validated its prognostic efficacy in another 3 independent datasets. The immune signature was also positively related to some important clinicopathologic features including TNM stage, metastasis risk, and tumor recurrence. More importantly, we found this signature was associated with estimated immune cell counts, immune score, gene levels of immune checkpoints, and mutation of critical genes in HCC.

The immune-related signature consisted of 12 IRGs with prognostic capacity for patients with HCC. In this signature, we observed that 10 of the 12 IRGs were cytokines or cytokine receptors, possibly functioning as important elements in the inflammatory process of tumorigenesis and tumor progression.[ This can be further elucidated by the concept that cytokines and their receptors can activate the oncogenic transcription factors of NFκB and STAT families to facilitate the development of cancer.[ Also, the tumor suppressive cytokines such as LECT2 may control the inflammatory phenotypes to constrain the growth of tumor.[ Therefore, the high-immune risk of our signature may represent an increased inflammatory status of the TIM, which may promote the progression of HCC and lead to poor OS of patients. Interestingly, our signature also contained CDK4, a fundamental driver of the cell cycle that is essential for the initiation and development of various tumors.[ Inhibition of CDK4/6 could augment antitumor immunity by promoting T-cell activation[; and, the selective CDK4/6 inhibitor palbociclib has shown encouraging results in preclinical models of HCC.[HSPA8 was another gene in the immune-signature that attracted our attention, which was associated with antigen processing and presentation.[ Inhibition of the ATPase activity of HSPA8 could enhance the immune response to protein antigens involved in cancer.[ These further shed light on the importance of our immune signature in the HCC microenvironment.

In order to understand the clinical values of the immune signature, we evaluated the relationship between IRG score and important clinicopathologic features. Patients with high IRG scores tended to have advanced TNM stage and increased risk of tumor metastasis and recurrence. These associations have also been documented in other studies exploring the prognostic immune signature in cancers such as renal papillary cell carcinoma and lung adenocarcinoma.[ It is likely that patients with high-immune risk are present with a TIM that can promote the cancer development, leading to advanced TNM stage and tumor relapse. Therefore, the proposed signature in our study can not only predict patients’ survival but reflect the probability of tumor progression, recurrence, and metastasis.

Additionally, we tried to classify the immune profiles and gene mutations associated with the IRG signature in HCC patients. Patients with high IRG score tended to have more T cells (e.g., CD8+ T cells) infiltration and higher estimated immune score, indicating an immunological microenvironment of HCC. Further analysis showed that the IRG score was positively correlated with the gene expression of PD-1, PD-L1, PD-L2, and CTLA-4. These results suggested that the function of T cells in the TIM of HCC was suppressed by the PD-1/PD-L1 or CTLA-4 signaling cascades, albeit with high infiltration of T cells. Thus, patients in the high-immune risk group may be more likely to benefit from the treatment of immune checkpoint inhibitors. Accordingly, Chen et al[ recently demonstrated that tumors with high PD-L1 expression and increased T cells infiltrating appeared to benefit more from blocking of immune checkpoints. Future studies are required to investigate the relationship between the immune signature and immunotherapy. Besides, we found that patients with mutation of critical genes in HCC had a higher IRG score. Apart from the well-known TP53, NCOR1 was also a mutated gene that had the mostly significant correlation with the immune signature. NCOR1 has been recognized as a new player on the field of T cell development,[ and it may exert a distinct role in liver regeneration and hepatocarcinogenesis.[ In our study, the mutation of NCOR1 caused a significant decrease in mRNA level (data not shown), which may contribute to the poor OS of HCC patients.[ These results further confirmed the prognostic ability of our signature.

Nevertheless, several limitations should be acknowledged in this study. First of all, the immune signature was developed using a retrospective design. Therefore, clinical validation of the IRG score is needed in prospective studies with large sample size. Secondly, due to the lacking of patients treated with immune checkpoint inhibitors, we cannot test the association of the immune signature with the response to immunotherapy. Thirdly, immunological studies focusing on the 12 IRGs individually or in combination should be conducted to explore their functions and support their clinical applications.

In summary, we developed an immune signature based on the expression data of IRGs from HCC patients. This signature showed robust ability in the prediction of patients’ OS and could reflect the risk of progression, metastasis, and recurrence of HCC. More importantly, the immune signature was associated with local immune profiles and mutations of critical genes in the TIM of HCC. These results may help us to understand the role of TIM in HCC and facilitate the discovery of novel biomarkers and therapeutic targets.

Author contributions

Jian Zhang conceived and designed the study. Qinghe Li and Bin Fan acquired and organized the gene expression data and relevant clinical information. Jun Ding and Xiaoxi Xiang performed the statistical analyses. Qinghe Li and Bin Fan wrote the manuscript, with key intellectual contents revised by Jian Zhang. All authors approved the final version of the manuscript.

Conceptualization: Jian Zhang.

Data curation: Jian Zhang.

Formal analysis: Qinghe Li, Bin Fan, Jian Zhang.

Investigation: Qinghe Li, Bin Fan.

Methodology: Qinghe Li, Bin Fan.

Resources: Qinghe Li, Bin Fan.

Software: Jun Ding, Xiaoxi Xiang.

Supervision: Jun Ding, Xiaoxi Xiang.

Validation: Jun Ding, Xiaoxi Xiang.

Visualization: Jian Zhang.

Writing – original draft: Qinghe Li.

Writing – review & editing: Jian Zhang.