An Individualized Immune Prognostic Index is a Superior Predictor of Survival of Hepatocellular Carcinoma.

BACKGROUND The tumor microenvironment is largely orchestrated by the immune cells. Considerable evidence has shown their excellent clinicopathological application value in assessment of clinical outcomes and immunotherapy efficacy. Hence, a moderate, individualized prognostic signature based on immune cells that can estimate prognosis and reflect the immune microenvironment in hepatocellular carcinoma (HCC) patients is greatly needed. MATERIAL AND METHODS Here, we systematically analyzed the expression differences and survival prediction value of tumor infiltrating immune cells by analyzing 638 HCC patients from 3 public cohorts, including 2 microarray datasets and 1 RNA sequencing dataset. CIBERSORT software, a computational algorithm, was used to calculate the relative levels of immune cells. Three immune microenvironment subtypes were defined via ConsensuClusterPlus package. Univariate and multivariate survival analyses were used to develop an individualized immune prognostic index based on immune cell pairs. RESULTS Notably, HCC patients with higher immune signatures score, utterly appreciable, suffered inferior prognosis (hazard ratio=2.742; 95% confidence interval: 1.887-3.983; P.

Background

Primary liver cancers are the sixth most common and fourth lethal malignant tumors worldwide, causing people great mental and social pressure [1]. According to Global Epidemiological Statistics of Cancer in 2018, there are a total of 841 080 new cases and 781 631 death resulting from liver cancer [1]. Hepatocellular carcinoma (HCC) is the most common pathological subtype of liver cancer. HCC is mainly induced by chronic hepatitis B virus (HBV) and chronic hepatitis C virus (HCV) infection, alcohol abuse, autoimmune hepatitis, diabetes, obesity, and some metabolic diseases [2]. Hence, the complexity of the mechanisms of pathogenesis of liver disease has led to the heterogeneity of the disease in different populations in different countries and regions, which poses a huge challenge for personalized medicine. Many patients at diagnosis have reached an advanced stage of treatment and lost the surgical option as treatment. Characterization of the tumor heterogeneity, and analyses to identify effective personalized monitoring predictors is urgently needed to provide a wealth of information and clinical decision support in HCC.

Recently, cancer immunotherapy has revolutionized oncotherapy. Some immune checkpoint inhibitors, which are effective in reinvigorating antitumor immune responses, are being clinically used for therapy of cancers [2-5]. Previously, a checkpoint inhibitors tremelimumab, which is the first immunotherapy assessed in HCC, has been proven as an encouraging antitumor agent for the detection of other checkpoint inhibitors [6,7]. Encouragingly, nivolumab, a programmed cell death protein-1 (PD-1) immune checkpoint inhibitor, also has shown promising effects in patients with advanced HCC [8]. Despite the great advances made for HCC immunotherapy, there are many challenges to overcome to better utilize the benefits of immunotherapy. For instance, the objective response rate of immune checkpoint inhibitors is still less than 20%, which requires accurate appraisal of patients most likely to respond to immune-based therapies [3,9]. Furthermore, comprehensively analysis of the immune landscape and compositions of tumor cell types could provide precision monitoring strategies for HCC patients. Recently, Tian et al. developed an immune-clinical prognostic index by integrating 5-feature-based immune signatures and clinical data to offer good prediction capability for early/intermediate HCC stage [10]. However, the prognostic value of immune cells needs to be explored in a larger cohort.

Here, we identified and characterized the immune tumor microenvironment of HCC. Furthermore, associations between phenotypes of tumor cells infiltration level and molecular characteristics were further analyzed. More importantly, an individualized immune prognostic signature based on immune cell pairs were proposed and validated. These approaches will provide a novel insight in the field of immunogenomics.

Material and Methods

Data acquirement

RNA-Seq data of 374 HCC and 50 non-tumor tissues were downloaded from The Cancer Genome Atlas (TCGA) liver hepatocellular carcinoma (LIHC) cohort by using TCGAbiolinks package in R software [11]. Furthermore, clinical follow-up information, tumor pathological parameters, and tumor mutation status also been downloaded from TCGA pan-cancer database [12]. We also retrospectively collected the public HCC microarray datasets with clinical follow-up information from gene expression omnibus (GEO) dataset, including GSE14520 and GSE76427.

Immune profiles estimation

From the pan-cancer atlas project in The Genomic Data Commons (GDC) data portal, 22 types of immune cells infiltration levels (T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, B cells naïve, B cells memory, NK cells resting, NK cells activated, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, neutrophils and eosinophils) were acquired [13]. The proportion of immune cells in TCGA tumor tissues were estimated using CIBERSORT software, which is a computational algorithm used to calculate the relative levels of immune cells [14]. Beefily, relative proportions of immune cells were assessed by using the gene expression profiles matched to a set of 22 immune cell reference profiles.

Clinical significance of tumor infiltrating immune cells

Mann-Whitney U test was used to analyze differences in immune cell components between HCC and non-tumor tissues. To provide a more moderate result and removed the effect of perioperative period on prognosis. HCC patients with overall survival (OS) not less than 90 days were submitted to survival analysis. Based on the median value of each type of immune cells, patients were divided into high immune cell infiltrating group (immune cells >median value) and low immune cell infiltrating group (immune cells ≤median value). Univariate analysis of was evaluated using log-rank test to observe the survival status in the two group. Immune cells with P-value in survival analysis were identified as prognosis associated immune cells. Different immune cell infiltration patterns were distinguishing by using hierarchical clustering algorithm (based on Spearman correlation analysis).

To further estimated the performance of prognosis associated immune cells in separating patients with different OS. Unsupervised clustering K-means was conducted to classify HCC into different group based on immune cells infiltration levels by using the ConsensuClusterPlus package in R software [15]. The procedure was repeated 1000 times to make the stratification more stability. The clustering method was K-means algorithm with the Euclidean distance and the number of clusters was identified through cumulative distribution function.

Individualized prognostic signature based on immune cells pair

In order to estimate the prognosis prediction performance of immune cells, we generate a score named immune cells pair (ICP) for presenting the relative expression level between 2 kinds of immune cells. The ICP was scored 0 or 1 under the following system: 1, immune cell 1 was less than immune cell 2; 0 immune cell 1 was not less than immune cell 2. In this way, we could provide a relative immune cells expression profile and compare the results across different detection platforms without the need for normalization. We identified 231 ICPs based on 22 types of immune cells. Some ICPs with constant values (0 or 1) in were removed for further analysis.

Log-rank test was conducted to select the ICP that was significantly correlated to the OS of HCC patients in the TCGA database. Significant ICPs were candidate for prognostic signature development. Then, multivariate cox analysis was conducted to develop a prognostic index named ICP index (ICPI). ICPI was generated based on the coefficient of each ICP multiplied by the ICP score.

Validation of the ICPI

To further validate the survival prediction performance of ICPI. We also collected 2 microarray data, including GSE14520 [16] and GSE76427 [17], from the Gene Expression Omnibus (). The gene expression patterns of GSE54236 profiled by GPL3921 were used for analysis. After removing HCC patients with OS less than 90 days, a total 218 and 91 HCC patients in GSE54236 and GSE76427 were included for validating the proposed prognostic signature. CIBERSORT was used to estimate the relative levels of 22 immune infiltrating cells in the 2 datasets.

Development of the immune-clinical prognostic signature

TNM is the traditional clinically application risk stratification system. Hence, by using the multivariate COX analysis, we finally integrated TNM stage, and ICPI to generate the immune-clinical prognostic signature. The prognostic performance of continuous prognostic signatures we proposed was estimated with area under curve (AUC) of time-dependent receiver operating characteristic (ROC) by using “timeROC” packages in R software.

Results

Immune landscape of hepatocellular carcinoma

Tumor immune microenvironment consists of diverse immune cells orchestrates the antitumor response of immunotherapy. For 374 HCC tissues and 50 non-tumor tissues, we comprehensively analyzed the quantitative expression values of 22 immune cells. 12 out of 22 immune cells were differentially expressed between non-tumor and HCC tissues (Table 1). Then, we depicted the interactions between tumor-immune cell, and their effects on the OS of HCC. According the Spearman correlation of hierarchical clustering analysis, the 22 immune cells mainly separated into 5 cluster (Figure 1A). Immune cells in the same cluster interacted with each other significantly (Figure 1B).

Table 1

Expression profiles of immune cells in HCC and non-tumor samples.

Immune cells	HCC tissues (n=374) Median (IQR)	Non-tumor tissues (n=50) Median (IQR)	P-value
B cells naive	0.0140 (0.0017–0.0346)	0.0592 (0.0421–0.0856)	2.05E-13
B cells memory	0 (0–0.0182)	0 (0–0)	5.07E-07
Plasma cells	0.0148 (0–0.0375)	0.0220 (0.0106–0.0356)	0.0634
T cells CD8	0.1062 (0.0675–0.1614)	0.1125 (0.0762–0.1615)	0.5811
T cells CD4 naive	0 (0–0)	0 (0–0)	0.0267
T cells CD4 memory resting	0.1212 (0.0538–0.1840)	0.1099 (0.0707–0.1674)	0.6935
T cells CD4 memory activated	0 (0–0)	0 (0–0)	0.3751
T cells follicular helper	0.0345 (0.0096–0.0678)	0.0321 (0.0146–0.0511)	0.5309
T cells regulatory Tregs	0.0324 (0.0031–0.0704)	0.0017 (0–0.0116)	4.06E-09
T cells gamma delta	0 (0–0)	0 (0–0.0091)	1.29E-05
NK cells resting	0 (0–0.0040)	0 (0–0.0149)	0.4046
NK cells activated	0.0584 (0.0297–0.0831)	0.0563 (0.0317–0.0761)	0.8726
Monocytes	0.0403 (0.0205–0.0695)	0.0635 (0.0445–0.1048)	1.20E-06
Macrophages M0	0.0346 (0–0.0926)	0 (0–0.0152)	3.60E-07
Macrophages M1	0.0483 (0.0257–0.0766)	0.0484 (0.0266–0.0876)	0.1337
Macrophages M2	0.2524 (0.1887–0.3250)	0.3140 (0.2423–0.3663)	0.0012
Dendritic cells resting	0.0035 (0–0.0144)	0 (0–0.0007)	6.49E-06
Dendritic cells activated	0 (0–0)	0 (0–0)	0.0604
Mast cells resting	0.0467 (0.0143–0.1103)	0 (0–0.0537)	1.93E-08
Mast cells activated	0 (0–0)	0.0097 (0–0.0529)	1.36E-15
Eosinophils	0 (0–0)	0 (0–0)	0.0601
Neutrophils	0 (0–0.0048)	0.0075 (0.0022–0.0132)	2.44E-08

HCC – hepatocellular carcinoma; IQR – interquartile range.

Figure 1

Characteristics of immune cells relationships in hepatocellular carcinoma (HCC). (A) Unsupervised clustering of immune cells. Twenty-two types of immune cells are divided into 5 groups based on the expression correlation. (B) Cellular interaction of tumor microenvironment cell types. The size of each cell represents survival impact of each immune cell, calculated used the formula log10 (Log-rank test P value). The lines connecting immune cells represent cellular interactions. Red and blue lines represent the positive correlation and negative correlation, respectively.

The scored immune cells profiles were used to characterize the relationships between HCC intratumoral immune states and OS. For each type of immune cell, survival analyses between high- and low-infiltration level were conducted based on the median value of immune cells infiltration level. There were 5 immune cells (neutrophils, macrophages M0, macrophages M2, T cells CD4 memory resting, and B cells naive) that were markedly correlated to the OS of HCC (Figure 2A), to further estimate the ability of immune cells in distinguishing patients’ clinical outcome. Similarly, 329 HCC patients were separated into 3 subgroups based on immune cells infiltration levels using unsupervised clustering K-means and patients in different immune clusters suffered distinct clinical outcome (Figure 2B). Individual tumor grade, stage, gender, and HCC driver genes (TP53 and CTNNB1) varied substantially in their proportion of immune subtypes (Figure 3A). These findings indicated that immune cells could be useful in the survival stratification of HCC patients.

Figure 2

Immune phenotypes of HCC patients suffered distinct clinical outcome. (A) The relationships between immune cells infiltration levels and overall survival. Univariate Cox analysis was conducted to estimate the survival difference between high- and low- immune infiltration level of each type immune cells. (B) Identification of different subtypes of HCC patients based on immune cell infiltration levels by using consensus clustering analysis. Three subtypes of HCC patients had significant different clinical outcome.

Figure 3

Immune landscape of hepatocellular carcinoma. (A) Consensus clustering of 329 HCC patients from The Cancer Genome Atlas (TCGA) database based on immune cells infiltration. Mutation status of TP53 and CTNNB1, gender, histological grade, as well as stage were annotated in the lower panel. Patients are stratified by immune cell pair index (ICPI): (B) TCGA database; (C) GSE14520 database; and (D) GSE76427 database.

Construction and Definition of the ICPI

Considering that the prospective survival monitoring function of immune cells, we attempted to develop an individual prognostic signature for precision survival monitoring. We used 22 immune cells to construct 231 ICPs. The associations of the 276 ICPs with OS were analyzed in the TCGA dataset, resulting in 32 prognostic ICPs. Subsequently, multivariate Cox analysis was conducted to construct the immune prognostic signature by the 32 ICPs. Finally, we generated the ICPI based on 9 unique ICPs and their coefficients in order to predict patient survival (Table 2).

Table 2

Model information about ICPI.

Immune cell 1	Immune cell 2	Coefficient
Plasma.cells	Neutrophils	0.655
T.cells.CD8	Eosinophils	2.894
T.cells.CD4.memory.resting	Macrophages.M0	0.557
T.cells.CD4.memory.resting	Dendritic.cells.activated	1.896
T.cells.CD4.memory.activated	Neutrophils	−0.773
T.cells.regulatoryTregs	Neutrophils	0.857
NK.cells.resting	Neutrophils	0.801
Monocytes	Macrophages.M0	0.359
Macrophages.M1	Dendritic.cells.resting	0.514

ICIP – immune cells pair index.

We further estimate the discriminatory power and model calibration of the ICPI we proposed. Kaplan-Meier (K-M) survival plot showed that HCC patients were subdivided into 2 groups with distinct OS based on the median value of ICPI in the TCGA database (hazard ratio [HR]=2.742; 95% confidence interval [CI]: 1.887–3.983, P<0.001; Figure 3B). In the other cohorts GSE14520 and GSE76427, the ICPI also could separate patients into 2 groups with distinct OS (GSE14520: HR=1.631, 95% CI: 1.051–2.532, P=0.024, Figure 3C; GSE76427: HR=2.696, 95% CI: 1.059–6.861, P=0.038, Figure 3D).

To leverage the complementary value of ICPI and traditional clinical TNM stage, we integrated the ICPI with TNM stage to build a composite prognostic signature. TNM stage was coded as I=1, II=2, III=3, IV=4. Then, immune-clinical prognostic index was proposed based on multivariate: immune-clinical prognostic index=0.879*ICPI+0.445*TNM stage. According to the analysis results of time-dependent ROC curve, the ICPI and displayed excellent discriminatory power at different timeline (Figure 4A). And the AUC of ICPI and immune-clinical prognostic index consistently higher than the traditional TNM staging system. By TimeROC analysis, we also compared the survival predictive power of the ICPI (Figure 4B), TNM (Figure 4C), and immune-clinical prognostic index (Figure 4D) at 1, 3, and 5 years. ICPI model showed higher performance in HCC patients’ survival prediction than TNM stage.

Figure 4

Area under the curves (AUC) of the immune cell pair index (ICPI), clinical TNM stage and immune-clinical prognostic index survival models. (A) Survival prediction performance of the three survival models at different time points. At 1-year, 3-year and 5-year, the prognostic value of ICPI (B), TNM stage (C), and immune-clinical prognostic index (D) were also estimated. Subgroup analyses estimating clinical prognostic value of ICPI. Hazard ratios (HR) >1.0 indicate that high ICPI is a risky prognostic biomarker (E).

Subgroup analysis of the prediction performance of ICPI

To further observe the discriminatory power of ICPI in different subgroup of patients with HCC. In TCGA database, the ICPI showed excellent prognostic evaluation in different subgroups. Interestingly, the ICPI is more effective in G1 and G2 grade than that in G3 and G4 grade. The ICPI is more effective in early stage than that in advanced stage (Figure 4E).

Discussion

The heterogeneity of infiltrating immune cells is increasingly recognized as the pivotal constituents of tumors [18,19]. Characterization of the immune microenvironment is indispensable for survival prediction and immunotherapy strategy selection in HCC patients. In the present study, the immune cells infiltration landscape, HCC immune subtype, as well as prognostic value of immune cells in HCC patients were comprehensively investigated. The highlight of this manuscript is we proposed an ICPI by analyzing the relative ranking of immune cells infiltration and only involves pairwise comparison within the immune cell infiltration levels of a sample. The strategy increased the comparability of different platforms and eliminated the need of normalization. Therefore, the ICPI we proposed could serve as an individualized, superior predictor of survival of HCC. Based on the survival model developed based on tumor immune microenvironment, we could observe the roles of immune cells in HCC patients and guide the individualized medicine in HCC.

In our analysis, we elucidated the comprehensive immune landscape and immune cells interactions of HCC by quantifying 22 types of immune cells. The first objective of this study was to investigate the clinical significance of immune cells, including expression profiles and prognostic value. More than half immune cells differentially expressed between HCC and non-tumor tissues. More importantly, high infiltration levels of neutrophils, macrophages M0, macrophages M2 indicated dismal prognosis. And high infiltration levels of T cells CD4 memory resting and B cells naïve suggested the favorable prognosis. High T cells CD8 also tended to be the protective factor. Identification moderate prognostic indicators is always the major project of research in tumors. Recently, immunotherapy herald new era of HCC treatment. Survival monitoring indicators related to the tumor immune microenvironment possess application prospects for therapy targets identification and guide clinical management of HCC in the era of cancer immunotherapy [20,21]. Li et al. [22] found that the density of neutrophils was higher in the peritumoral compared with tumors and increased intratumoral neutrophils were markedly associated with decreased recurrence-free survival and OS after adjusting other clinical parameters. Our results also confirmed these findings. A recent meta-analysis included 3509 HCC patients suggested that high levels of intratumoral CD8+ lymphocytes were significantly correlated with favorable OS and disease-free survival of HCC [23]. Macrophages M2, has also been observed in HCC and may be associated with poor survival [24,25]. These studies focused on the prognostic value of immune cells have significant clinical practical value. Hence, the global immune landscape analysis for HCC is necessary. On the basis of previous frameworks, we systematically estimate the interactions between immune cells and proposed a novel effective prognostic signature for HCC patients.

As aforementioned, many studies have validated that an immune cell alone could be used as a biomarker for HCC diagnosis and prognosis evaluation. However, owing to the heterogeneity of the HCC, a single molecular event as the prognosis predictors tends to be unstable. Several previous studies have been trying to combine multiple molecular events that improve the sensitivity of diagnosis and prognosis [26,27]. Furthermore, in the present study, computational algorithm is bound to be some bias in the tumor cells infiltering. Here, we proposed an ICPI based on the relative ranking of immune cells infiltration. And only a 2–2 comparison of the immune cell infiltration level in each sample, thereby eliminating the need for data normalization. These processed also minimize the influence of bias. Similarly, Li et al. [28] constructed a prognostic signature based on immune-related genes. They constructed immune gene pair to estimate the heterogeneity of gene expression analysis method. Hence, the prognostic signature is effective in early-stage non squamous non-small cell lung cancer [28]. Subsequently, the approach was used in several cancers and displayed excellent performance [29,30]. The prognostic signature we proposed here is moderate in HCC patients’ survival prediction. More importantly, the prognostic efficacy was also validated by 2 independent microarray datasets. Cross-platform study also proves that the results are reliable. Interestingly, subgroup analysis showed that the risk stratification performance is particularly significant in early HCC either the training set or the validation set. These findings are meaningful for the early stage HCC patients.

Limitations of our study should be noted. First, the retrospective nature of the study limited its application, although we have used different datasets to validate its performance. Second, the great heterogeneity of immune microenvironment is difficult to be assessed accurately. Although we used relative ranking of immune cells infiltration to reduce certain bias, some bias may still remain. And future studies with biological analysis will provide a more comprehensive view in the survival monitoring.

Conclusions

In conclusion, we developed and validated an effective prognostic signature based on immune cell pairs, which is excellent in predicting HCC patients’ clinical outcome. This study provides a systematic view of the immune cells characteristics in HCC and suggests their superior survival monitoring performance.