Significance of tumor-infiltrating immunocytes for predicting prognosis of hepatitis B virus-related hepatocellular carcinoma

Qi-Feng Chen, Wang Li, Pei-Hong Wu, Lu-Jun Shen, Zi-Lin Huang

Abstract

Key words: Immune risk score; Hepatitis B virus; Hepatocellular carcinoma; Prognostic signature; Cell type identification by estimating relative subsets of RNA transcripts

INTRODUCTION

Although the screening, diagnostic, and therapeutic techniques for hepatocellular carcinoma (HCC) have progressed, it remains a major cause of cancer-related death worldwide, with dismal clinical outcomes[1]. Hepatitis viruses have been epidemiologically related to HCC development[2,3], of which hepatitis B virus (HBV)has been recognized as a major factor for HCC in Asian countries[4]. Typically, HBV infection is more prevalent in Eastern Asia, such as China, where over two-thirds of patients are detected[5].

Although the etiological link between HBV and HCC has been well established, the effect of immunocytes on tumor-related microenvironment remains largely unclear.The tumor-infiltrating lymphocytes can serve as prognostic factors for HCC, and are a research hotspot[6,7]. Nonetheless, the non-lymphocyte immunocytes have also been recognized to play crucial roles in prognosis of various types of cancers[8]. Therefore, it is critical to comprehensively predict the immune-related prognosis.

The traditional approaches to measure immune infiltrate in tumors, like immunohistochemistry (IHC) and flow cytometry, cannot thoroughly evaluate the immune functions of various types of cells or effectively distinguish closely related cell populations, since few immune markers can be simultaneously detected using the existing technology. In recent years, the algorithm cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT), known as the metagene means, is recognized as the method with the highest accuracy. Typically, CIBERSORT can accurately estimate 22 human immunocyte phenotypes[9].

This study explored the proportions of immunocytes by appropriately selecting HBV-related HCC cases, so as to identify candidate immune biomarkers using The Cancer Genome Atlas data. The least absolute shrinkage and selection operator(LASSO) algorithm was employed to determine key immunocytes. Thereafter, a HCC immunoscore system was formulated, aiming to establish an immune-related marker that could be used to predict the survival of HBV-related HCC patients.

MATERIALS AND METHODS

Patient datasets and processing

Data were downloaded from the TCGA database. The data transfer tool of the Genomic Data Commons applications was employed to download gene expression profiles as well as clinical information of HCC patients (https://tcgadata.nci.nih.gov/). HCC patients with seroprevalence of hepatitis B surface antigen(HBsAg) were included. Among them, 58 patients were HBV-related, but patients with CIBERSORT P ≥ 0.05 were excluded (n = 5), leaving 53 HCC cases in the final cohort for analysis. The flow chart for this analysis is shown in Figure 1. A total of 53 HBV-related HCC cases were enrolled for comprehensive integrated and functional analyses. The mRNA expression in HCC samples was derived from the Illumina HiSeq RNASeq platform and normalized based on TCGA. All data were publicly available and open access. Hence, approval from the ethics committee was not required. All data were processed in accordance with the guidelines from data access policy as well as the NIH TCGA human subject protection (http://cancergenome.nih.gov/publications/publicationguidelines).

Immunocyte type fraction estimation

The CIBERSORT method was used to quantify immunocyte fractions[9], which could sensitively and specifically distinguish 22 phenotypes of human hematopoietic cells.Pornpimol et al[10,11]had used the deconvolution method through the CIBERSORT tool,with the threshold set at P ＜ 0.05; their results suggested that CIBERSORT could accurately predict the immunocyte population fractions. Hence, only those possessing the CIBERSORT P ＜ 0.05 were deemed as eligible in subsequent analysis. Moreover,the immunocyte proportions were separately estimated for 76 samples. For all samples, each estimate of immunocyte type fractions added up to be 1.

Immunoscore signature construction

The optimal threshold for each immunocyte fraction was generated using the survminer package. Then, the value of immunocyte fraction was assigned as 0 or 1;specifically, the value 0 was given to the cell type fraction below the corresponding threshold, otherwise the value 1 was assigned. The correlation of each immunocyte expression with the overall survival (OS) of patients was calculated using the univariate Cox model. Then, the immunocytes were screened, followed by verification through LASSO regression using “glmnet” package in R project. Finally, an immunoscore was constructed based on immunocytes by linearly combining the multiplication of regression model (β) by expression level as follows: Immunoscore =β1 × immunocyte 1 expression level + β2 × immunocyte 2 expression level + ····· + βn× immunocyte n expression level. In addition, the time dependent receiver operating characteristic (ROC) curves were employed to evaluate the accuracy of prognosis prediction in five years, which was achieved by the comparisons of sensitivity as well as specificity in predicting the survival using the immunoscore.

Confirmation of the immunoscore signature

The patients were distributed based on the immunoscore and their survival information. Furthermore, all cases were classified into a high- or low-risk group according to the immunoscore threshold. In addition, related Kaplan-Meier survival curves for patients with high or low predicted risk were plotted. Thereafter,univariate as well as multivariate Cox proportional hazard regression analyses were performed. We developed a nomogram with “rms” package in R project. The prognosis prediction accuracy of the nomogram was verified for discrimination and calibration. Each test was two-sided, and a difference of P ＜ 0.05 was deemed statistically significant. The R (version 3.4.1; R Foundation) software was adopted for all analyses.

Functional analysis of the prognostic immunoscore

Figure 1 General design of this study.

Inhibitory checkpoint molecules can be potentially used in numerous cancer types,which are thereby recognized as the immunotherapy targets for cancer. The correlations of the immunoscore with corresponding genes were further investigated.In addition, the correlation of the immunoscore with the determined immune markers was analyzed through Pearson’s correlation analysis based on their expression levels.The mechanisms underlying “Molecular Signatures Database” for the c2.cp.kegg.v6.2.symbols were investigated according to gene set enrichment analysis(GSEA)[12]with JAVA procedure (http://software.broadinstitute.org/gsea/index.jsp).The random sample permutation number was deemed at 1000 and a significant level was P ＜ 0.05.

RESULTS

Patient characteristics and composition of immunocytes

Detailed patient characteristics are listed in Table 1. The immunocyte compositions of cancer and normal tissues were analyzed. Figure 2A summarizes immunocyte composition in the entire cohort. Resting CD4 memory T cells ranked the first,accounting for 25% on average; while CD4 naive T cells and gamma delta T cells were not detected. Heatmap summarizing the distribution of immunocytes was also plotted (Figure 2B). Various immunocyte subpopulation fractions showed a low to moderate correlation (Figure 2C). Relative to para-carcinoma tissues, higher proportions of regulatory T cells as well as resting dendritic cells could be detected in cancer tissues, along with lower proportions of plasma cells, resting NK cells, M2 macrophages, and neutrophils (Figure 2D, P ＜ 0.05).

Immunoscore derivation

Table 2 displays the thresholds of all cell types. Figure 3A presents the forest plot of the correlation of 18 immunocyte subsets with OS. Through the LASSO algorithm(Figure 3B and C ), eight types of immunocytes (including plasma cells, naive B cells,CD4 memory resting T cells, M1 macrophages, activated NK cells, resting dendritic cells, M2 macrophages, and neutrophils) were selected to build the prognostic immunoscore model. The immunoscore was calculated according to the formula:(2.764 × fraction level of naive B cells) - (0.371 × fraction level of plasma cells) + (1.694× fraction level of resting CD4 memory T cells) + (21.45 × fraction level of activated NK cells) + (1.669 × fraction level of M1 macrophages) - (1.649 × fraction level of M2 macrophages) + (0.504 × fraction level of resting dendritic cells) + (3.647 × fraction level of neutrophils). Moreover, the immunoscore prognosis performance, as the continuous variable, was also studied within the cohort through time-dependent ROC analysis. The areas under the curves of the immunoscore biomarker prognostic model were 0.971, 0.912, and 0.975 for 1-, 3-, and 5-year OS, respectively (Figure 3D).

Table 1 Baseline patient characteristics

Confirmation of the immunoscore

The patients were distributed based on the immunoscore and their survival information (Figure 4A). Subsequently, the patients were classified into high vs low immunoscore groups according to the threshold (36.47) acquired using survminer package. The Kaplan-Meier curves for both groups are presented in Figure 4B. The patients with a high immunoscore were associated with worse OS than those with a low immunoscore [hazard ratio (HR) = 66.007, 95% confidence interval (CI): 8.361-521.105; P ＜ 0.0001]. The HR of the immunoscore was 2.010 (95%CI: 1.487-2.716) based on the univariate Cox proportional hazard regression model, and similar results were obtained from the multivariate Cox proportional hazard regression model (HR =2.997, 95%CI: 1.737-5.170) after adjustment for clinical covariables (Table 3).

Additionally, a prognostic nomogram incorporating the above prognostic factors for visualization and prediction was built to identify a composite factor to predict the prognosis for HCC cases (Table 3). Notably, the nomogram was constituted by six factors (Figure 5A), with a C-index of 0.757 (95%CI: 0.648-0.866). Meanwhile, the 1-, 3-,and 5-year calibration curves of the recurrence probability were consistent between predicted values by the nomogram and measured values (Figure 5B-D).

Functional analysis of the prognostic immunoscore

Gene expression profiles were examined for exploring the potential immunoscore model-related biological phenotypes. First, the correlation of the immunoscore with the screened immune-related gene expression was particularly emphasized. Figure 6 displays the significant negative correlations of the immunoscore with programmed death (PD)-1 (P = 0.024), PD-L1 (P = 0.026), PD-L2 (P = 0.029), and CD27 (P = 0.033)expression. Second, GSEA was performed for illustrating the biological roles for the prepared immunoscore model (Supplemental Table 1). The highly expressed genes showed significant enrichment in pathways including glycosaminoglycan chondroitin sulfate biosynthesis, proximal tubule bicarbonate reclamation, focal adhesion,autoimmune thyroid disease, the Notch signaling pathway, basal transcription factors, DNA replication, and homologous recombination (Figure 7).

DISCUSSION

The CIBERSORT was utilized for the quantification of composition of 22 immunocytes for HBV-related HCC. Furthermore, the immunoscore was developed, which was based on eight immunocyte fractions and served as a new tool to predict the prognosis of HBV-related HCC survival. The prognostic significance of the immunoscore signature was confirmed by the time-dependent ROC curve, Kaplan-Meier, and Cox regression analyses. The results suggested that the established immunoscore system displayed favorable accuracy for estimating OS of HBV-related HCC by nomogram analysis. Further functional analysis uncovered several immune checkpoints and pathways that played crucial roles in the immunoscore system during the progression of HCC.

There are few available immunoscore-based models for the prognosis of HCC. In 2015, Sun et al[13]confirmed that center tumor CD8+ T cells could predict the prognosisfor HCC patients, and were a valuable marker for predicting patient survival and tumor recurrence. However, only CD3+ and CD8+ T cells were measured by IHC in their research. Similarly, in 2016, Gabrielson et al[14]concluded that immunoscore as well as PD-L1 expression could be used as valuable markers to predict the prognosis for HCC patients undergoing resection of primary tumor. However, regarding etiology of HCC, only 16.9% (11/65) of cases suffered from HBV infection in their study. In 2016, Petrizzo demonstrated that immunoscore was a meaningful marker to predict the prognosis for HCC patients undergoing resection for primary tumor[15].However, 65 HCC patients with unknown virus condition were reviewed in that study, and IHC was only performed with anti-CD3, CD8, and PD-L1 monoclonal antibodies. In 2017, Yao et al[16]conducted a cohort study by recruiting 92 consecutive HBV-HCC patients, and concluded that the immunoscore classification showed close relation to patient outcome. However, only CD3+, CD8+, and CD45RO+ T cells were measured by IHC in that study.

Table 2 Thresholds for immunocyte fractions

The current study was superior to previous studies examining the prognosis value of immunoscore for HCC in several aspects. First, assessing CD3+ and CD8+lymphocytes alone could not comprehensively represent local immune status. Prior research was limited by its narrow view of immunological reaction due to technical limitations, and immunocytes were evaluated through IHC based on a single surface marker, which was ineffective for distinguishing the tightly correlated cell types, since numerous marker proteins would be expressed among various types of cells.Technically, IHC is limited to existing phenotypic markers, and may lead to inconsistent findings in clinical studies. However, using transcriptomic data for computationally describing tumor microenvironment can serve as a promising tool to overcome the technical limitations related to IHC, and can be more readily employed for the further characterization of 22 different immune populations than IHC. Second,homogeneity in the case cohort was largely enhanced through pre-analysis case selection. Viral and non-viral HCC were suggested to be associated with different prognoses, and different treatment strategies were recommended; therefore, they must be distinctively analyzed to develop the models for prognosis[17-20]. In contrast to previous studies that combined viral and non-viral related HCC, the current study was conducted in HBV-positive patients only. Third, in terms of the methodology,LASSO penalized regression could increase the accuracy of bioinformatic analysis,and allowed for simultaneous analysis of each independent variable to select the most valuable parameters[21]. Therefore, such approach would show much higher accuracy than stepwise regression of the multivariate Cox model[22]. Fourth, an easy-to-use nomogram model based on Cox analysis was constructed in our meticulously chosenpatient cohort with clinicopathological characteristics that might facilitate individualized prediction of patients.

Figure 3

Figure 3 Immunoscore model construction. A: Forest plots presenting the relationships of various immunocyte subpopulations with overall survival, in which the unadjusted hazard ratios as well as the corresponding 95% confidence intervals are displayed; B: Coefficient profiles of 20 immunocyte type fractions using the least absolute shrinkage and selection operator (LASSO); C: 10-fold cross-validation to select parameters for LASSO model, which can determine the cell type numbers adopted for LASSO model; D: Immunoscore determined using the time dependent receiver operating characteristic curves. For the 1-, 3-, and 5-year OS, the AUCs for the immunoscore were 0.971, 0.912, and 0.975, respectively. aP ＜ 0.05. LASSO: Least absolute shrinkage and selection operator; AUC: Areas under the curve.

Immune checkpoint inhibitors have shown significant therapeutic responses to tumors. Clinical trials of antibodies targeting the immune checkpoint inhibitors to treat advanced HCC are ongoing[23]. Our results showed that the immunoscore was negatively correlated with PD-1, PD-L1, PD-L2, and CD27 (Figure 6). There are two ligands for PD-1 receptor, namely, PD-L1 and PD-L2[24,25]. Our previous studies found that excessive PD-1-positive intratumoral lymphocytes could predict the cytokineinduced killer cell survival benefits for HCC patients[26], and serum soluble PD-1 and PD-L1 can independently predict the prognosis for HCC patients in distinct ways[27].CD27 is a cytokine that facilitates naive T cell expansion specific to antigen, which is of vital importance to generate T cell memory[28]. Moreover, pathway enrichment analysis suggested that such immunocytes could affect HCC genesis as well as progression by eight pathways whose biological functions have been previously reported in HCC. Most of these pathways are canonical and crucial, which participate in HCC genesis as well as progression.

This study had several limitations. First, information, such as acute infection,immune system diseases, as well as anti-inflammatory drugs, could not be obtained since this study was conducted based on the publicly available datasets, and this could potentially result in error or bias. Second, each case was retrospectively chosen,so an underlying bias may exist related to the imbalanced clinicopathological characteristics with treatment heterogeneity. Third, all data were downloaded from Western countries to establish the models, so the findings should be cautiously applied to the Asian population. Finally, further large-sample research is needed for external validation of our results.

In summary, this study uncovered the immunocyte signature in HBV-related HCC.The proposed immunoscore model may comprehensively provide the required clinical data to improve the personalized treatment for HBV-related HCC cases.

Table 3 Cox regression analysis for overall survival-related clinicopathological features

Figure 4 Verification of the immunoscore signature to predict the prognosis for hepatitis B virus-related hepatocellular carcinoma. A: The vital status and overall survival were distributed based on the immunoscore; B: Kaplan-Meier curve regarding the overall survival of high- vs low- immunoscore groups.

Figure 5 Construction of the nomogram. A: Nomogram to predict the 1-, 3-, and 5-year overall survival of hepatitis B virus-related hepatocellular carcinoma patients based on the immunoscore and clinicopathological parameters; Calibration curves for predicting patient survival at 1 (B), 3 (C), and 5 (D) years.

Figure 6 Scatter plots describing the relationship of immunoscore with the expression of immune checkpoint regulators, including PD-1 (A), PD-L1 (B),PD-L2 (C), and CD27 (D).

Figure 7

Figure 7 Gene set enrichment analysis delineating the biological pathways related to the immunoscore values using the “c2.cp.kegg.v6.2.symbols” gene sets.

ARTICLE HIGHLIGHTS