Introduction
Liver cancer is the fourth leading cause of cancer deaths worldwide. Its incidence is increasing most rapidly in the USA, with 42,810 new cases estimated in 2020. Hepatocellular carcinoma (HCC), the most common form of liver cancer, is correlated with chronic infection with hepatitis B or C virus, aflatoxin-contaminated foodstuff consumption, and heavy alcohol intake.1,2 However, due to the universal late diagnosis, most HCC patients are unlikely to be amenable to potentially curative treatments.3 Some multi-kinase inhibitors are currently approved by the U.S. Food and Drug Administration for advanced HCCs. However, use of sorafenib or regorafenib has only improved the median overall survival (OS) duration by fewer than 3 months in advanced HCC patients, whose prognoses remain dismal.4,5 With an increasing number of HCC patients experiencing resistance to targeted therapy, efforts to develop novel therapy for HCC has drawn attention to the HCC microenvironment.6
Recent evidence has highlighted the importance of the tumor immune microenvironment for HCC outcome.7–9 Clinical trials of agents that target immune checkpoints such as programmed death 1 (PD-1, programmed death-ligand 1 (PD-L1), and cytotoxic T-lymphocyte-associate protein 4 (CTLA4) in HCC patients are ongoing, with an overall treatment response rate of less than 20%.10–12 The poor response to immunotherapy demonstrates that current patient-stratification methods are not sufficiently effective. The HCC immune microenvironment remains poorly understood and effective immune biomarkers for the evaluation of HCC prognosis are lacking. Identifying immune-related signatures for HCC patients that can serve to predicting the outcomes of and guide immunotherapy for HCC is needed.
We aimed to add to the knowledge of the association of immune-related genes with HCC prognosis and response to immunotherapy, which will set the groundwork for studies of the identified genes and provide a practical tool to optimize clinical HCC management. We studied 577 HCC cases from one cohort each in two public databases, identifying and validating immune-related gene set signature, IGSHCC as a robust and independent prognostic classifier for HCC with greater predictive accuracy than other prognostic models. We also established an IGSHCC-based nomogram to confirm the clinical applicability of our signature. We discovered a close association of IGSHCC score with tumor immunocyte infiltration status and response to immunotherapy and confirmed the utility of IGSHCC in prognostic risk-stratification and immunotherapy response prediction in another independent HCC cohort receiving anti-PD-1 therapy. Overall, the study finding support IGSHCC as a robust signature with potential clinical significance, given its value in individualized prognostic prediction and immunotherapy guidance for HCC.
Methods
Study populations and gene expression data preprocessing
The gene expression profiles and clinical information for HCC cases were collected from two public microarray datasets, The Cancer Genome Atlas (TCGA)13 and series GSE14520 in the Gene Expression Omnibus (GEO) from the Liver Cancer Institute (LCI) at Zhongshan Hospital, Fudan University.14 Only patients with available follow-up data, gene expression data, and clinical profiles were included. A total of 356 cases from TCGA and 221 cases from the GSE14520 were analyzed in our study.
Level 3 TCGA RNA sequencing data for HCC cases were downloaded and logarithmic-transformed. Gene expression data from GSE14520 were processed on the Affymetrix GeneChip Human Genome U133A 2.0 platform. The difference in the batch effect between the two datasets was adjusted by the ComBat method.15 To ensure consistency in gene names across the cohorts, expression profiles were collapsed from probe level to the corresponding gene symbols based on the annotation platform of each set.
The 356 patients in the TCGA dataset were randomly assigned to training (n=250) and internal validation (n=106) datasets. GSE14520 served as an external validation dataset (LCI validation set 1). Another independent cohort of HCC patient who underwent hepatectomy at Zhongshan Hospital in 2012 (n=56) was retrospectively evaluated as a polymerase chain reaction (PCR) array validation dataset (LCI validation set 2) to further test the reliability of IGSHCC. The clinical characteristics of the patients in the four datasets are summarized in Supplementary Table 1. This study was approved by the ethics committee of Zhongshan Hospital (approval number: B2021-143R), and written informed consent was obtained from all patients. The study was conducted following the ethical standards of the Declaration of Helsinki.
Immune-related gene set definition
A comprehensive list of immune-related genes was downloaded from the ImmPort database.16 The list consists of 1,534 genes that are involved in various immune pathways, including antimicrobial pathways, the T-cell receptor signaling pathway, interleukins, and the tumor necrosis factor family receptor signaling pathways. This list was used as a gene pool from which we screened out key genes to develop our immune-related prognostic signature. Among the genes in the list, 837 whose transcriptional expression levels could be obtained across all platforms, were selected for further analysis.
Development of the immune-related gene signature for HCC
LASSO, a method of regression analysis with high-dimensional factors, has been extended, and is used in Cox proportional hazards regression models for survival analysis.17,18 After application of univariate Cox regression analysis to the TCGA training dataset to identify candidate genes associated with prognosis for HCC (p<0.01) in the immune-related gene list, the LASSO Cox regression model was used with the glmnet package in the R computing language to screen the most powerful prognostic genes from the candidate genes. The best penalty parameter lambda was determined using 10-fold cross-validation. Using the formula below, a weighted mRNA level-based value, IGSHCC, was identified, and normalized to predict prognosis for HCC:
Prognostic score=∑i=1Nwixi
where N, w, and x represent the number of the selected genes, coefficient value, and gene expression level, respectively. The formula was used to calculate the prognostic scores of HCC patients. The information about and corresponding weight coefficients for 22 identified immune-related genes are listed in Supplementary Table 2. The cutoff value for the low- and high-IGSHCC groups was determined to be the median IGSHCC score in the training set and was used for further validation analysis.Validation of the immune-related gene signature for HCC
The IGSHCC score was further in the internal and the external validation sets, and patients were assigned to the high- or low-IGSHCC group based on the cutoff value described above. Univariate and multivariate Cox regression analyses were performed to evaluate the independent predictive ability of IGSHCC, age, sex, α-fetoprotein (AFP) level, tumor grade, and TNM stage for OS in the HCC patients. TCGA cases were placed in several subgroups according to TP53 genotype, CTNNB1 genotype, AFP level, age, tumor grade, and TNM stage. The performance of IGSHCC was then evaluated in the subgroups. The OS difference in the high- and low-IGSHCC groups was tested for significance by the Kaplan-Meier method and log-rank test.
Existing HCC prognostic signatures for comparison with IGSHCC
Three other published prognostic signatures were retrospectively collected for comparison with IGSHCC (Supplementary Table 3).19–21 The gene numbers in these signatures ranged from four to 65. The score for each signature was calculated by the corresponding formula, and a receiver operating characteristic (ROC) curve analysis was used to compare the predictive efficiency of different signatures using the pROC package in the R language. The area under the ROC curve (AUC) was calculated for all signatures, and the prognostic signature with the highest AUC value was considered the most efficient one.
Construction and evaluation of the IGSHCC-based nomogram
In cancer care, nomograms are models commonly used to predict prognosis. Nomograms generate an individual numerical probability of a clinical event by integrating diverse variables, fulfilling our desire for biological and clinical integration of models for prognosis prediction and assisting our drive toward personalized cancer medicine.22 Thus, IGSHCC was integrated with age and TNM stage, independent prognostic indicators for HCC identified via multivariate Cox regression analysis of the entire TCGA dataset, to establish a prognostic nomogram for predicting the likelihood of 3- and 5-year OS using the rms package in the R language.
Subsequently, to evaluate the predictive performance and clinical practicability of the nomogram, calibration plots and decision curve analysis (DCA) were performed in the entire TCGA and LCI datasets. Calibration plots were helpful in assessing the consistency of the actual prognoses with the predicted outcomes using the nomogram. In addition, time-dependent ROC analysis was carried out to evaluate the OS-predictive accuracy of the nomogram.
Single-sample gene set enrichment analysis (ssGSEA)
ssGSEA was conducted as described previously.23 Twenty-four immunocyte types that could be discriminated by gene sets specifically overexpressed in each immunocyte using a deconvolution approach were included in this analysis. The immunocytes are involved in both innate immunity (natural killer cells, CD56dim natural killer cells, CD56bright natural killer cells, plasmacytoid dendritic cells [DCs], immature DCs, activated DCs, neutrophils, mast cells, eosinophils, and macrophages) and adaptive immunity (B cells, central memory T cells, effector memory T cells, activated CD8+ T cells, T follicular helper cells, γδ T cells, T helper [Th] cells, type 1 (Th1), Th2, and Th17, and regulatory T cells). The T-cell infiltration score was defined as the average of the standardized ssGSEA values for CD8+ T, central memory T, effector memory T, Th1, Th2, Th17, and regulatory T cells. The cytotoxic cell infiltration score was calculated according to the geometrical mean of perforin (PRF1) and granzyme A (GZMA) expression levels.24 The DC infiltration score and macrophage infiltration score were defined as the standardized ssGSEA values for DC subsets (total, plasmacytoid, immature, and activated) and macrophages, respectively.
Molecular characterization and immunotherapy response prediction
To describe the difference in molecular and signaling pathway characterization between the high- and low-IGSHCC groups, gene set variation analysis25 was applied to our dataset with the C7 (immunological) and H (hallmark) gene sets downloaded from the Molecular Signatures Database.26 For the C7 gene set, the significant pathway sets were identified according to p-values <0.001 and fold-change values of at least 1.85. For the H gene set, the significant pathways were identified according to t-values of at least 2. In addition, the tumor immune dysfunction and exclusion (TIDE) algorithm and subclass mapping were used to predict the response of cancer to immunotherapy, such as anti-PD-1 and anti-CTLA4 therapy.27,28
Real-time PCR array analysis
An RT2 Profiler Custom PCR array (Qiagen, Hilden, Germany) was conducted to simultaneously examine the mRNA expression levels of 22 immune-related genes in our signature. The array was designed and produced by BioTNT and performed following the manufacturer’s instructions. Briefly, total RNA was extracted from HCC samples using TRIzol reagent (Invitrogen, Waltham, MA, USA), and cDNAs were prepared using an iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA). The PCR primers for the 22 immune-related genes are listed in Supplementary Table 4. Real-time PCR analysis was performed using RT2 SYBR Green qPCR Master Mix (Qiagen) with a ViiA 7 (384-well block) instrument (Thermo Fisher Scientific, Waltham, MA, USA). The relative expression levels of the target genes were calculated using the ΔΔCt method.
Histological and immunofluorescent analyses
HCC tissue sections were prepared and incubated with a rabbit anti-PD-L1 antibody (ab205921; Abcam, Cambridge, UK) and mouse anti-CD8 antibody (ab199016; Abcam). Next, the slides were incubated with goat antirabbit IgG antibody (Alexa Fluor 488; Abcam), goat antimouse IgG antibody (Alexa Fluor 647; Abcam), and 4′,6-diamidino-2-phenylindole (Invitrogen). Images of the sections were acquired using a confocal microscope with 40x objectives and analyzed using image processing software (Olympus Corp., Tokyo, Japan).
Statistical analysis
All statistical analyses were performed in the R language (version 3.6.0) and SPSS software (version 23.0; IBM Corp., Armonk, NY, USA). Continuous variables were reported as means ± standard deviation and were compared using Student t-test. Categorical variables were reported as frequencies and percentages (%) and compared using the Pearson chi-square test. A Cox proportional hazards model was used to conduct standard univariate and multivariate analyses. Prognostic differences between groups or subgroups were analyzed using the Kaplan-Meier method and log-rank tests. For all tests, two-tailed p-values <0.05 were considered significant.
Results
IGSHCC is a powerful prognostic signature with high predictive accuracy
The study design and workflow are shown in Figure 1. Univariate Cox regression analysis identified 78 genes from the immune-related gene list that were significantly associated with the OS of HCC patients (p<0.01) in the TCGA training set as the candidate genes for signature construction (Supplementary Table 5). Subsequently, by applying the LASSO Cox regression model to these genes, we constructed a prognostic signature consisting of 22 immune-related genes (Supplementary Fig. 1A, B). Most of the genes in the signature were cytokines, chemokines, and their receptors. To varying degrees, they played roles in the transforming growth factor-β signaling pathway, antigen processing and presentation, and adaptive immune response (Supplementary Fig. 1C). Detailed information on and corresponding coefficients for the 22 immune-related genes are presented in Supplementary Table 2 and Figure 1D. We calculated the IGSHCC score as the sum of the weighted expression levels for the 22 genes. To guarantee the consistency of our study findings, the cutoff value for high- and low-IGSHCC group was set as the median IGSHCC score in the TCGA training dataset and fixed as 1.291 during further validation processes (Supplementary Fig. 2).
Kaplan-Meier survival curve analysis (Fig. 2A) demonstrated an significant prognostic difference in OS between the high- and low-IGSHCC groups in the TCGA training dataset [hazard ratio (HR), 3.87, 95% confidence interval (CI), 2.44–6.14]; p<0.0001). We also confirmed the robust predictive performance of IGSHCC in the internal validation dataset (Fig. 2B), the entire TCGA dataset (Fig. 2C), and the LCI validation dataset (Fig. 2D). Of note, in the LCI validation dataset, the high- and low-IGSHCC groups also had markedly different clinical outcomes (HR: 1.99 95% CI: 1.30–3.05; p=0.0013). Additionally, we compared the expression levels of the 22 immune-related genes in the high- and low-IGSHCC groups in the four datasets (Supplementary Fig. 2). The heatmaps in Supplementary Figure 2 show that the expression levels of genes such as BMP6, GLP1R, and GMFB were higher in the high-IGSHCC group than in the low-IGSHCC group. In contrast, tumors with low IGSHCC scores had much higher expression of genes like CCR8, TNF, and RETN than did those with high scores. Taken together, the results demonstrate that the signature we constructed was a powerful and effective tool for prognostic risk prediction for HCC patients.
IGSHCC score is an independent risk factor for HCC prognosis
To confirm the independent prognostic value of our signature, we performed univariate and multivariate Cox regression analyses of the TCGA training dataset (Fig. 3A, B). In addition to IGSHCC, we included several common prognostic factors, like age, sex, AFP level, tumor grade, and TNM stage, in the analysis. The univariate regression model demonstrated that IGSHCC was significantly associated with OS (HR: 4.11, 95% CI: 3.03–5.57; p<0.0001). Furthermore, in the multivariate regression model, IGSHCC remained a strong independent prognostic predictor of OS (HR: 5.51, 95% CI: 3.40–8.96; p<0.0001). The results were consistent with those obtained from the internal validation dataset (Fig. 3C, D), LCI validation set (Fig. 3E, F), and the entire TCGA dataset (Supplementary Fig. 3A, B). In the LCI validation cohort, IGSHCC was still considered a significant prognostic variable (univariate p=0.001 and multivariate p=0.035). The results revealed that the IGSHCC score was a powerful independent risk factor for HCC prognosis.
To determine whether IGSHCC retained its prognostic value for OS in HCC patient subgroups with different molecular or clinical characteristics, we placed the patients in the entire TCGA dataset into different subgroups according to TP53 genotype, CTNN1B genotype, AFP level, age, tumor grade, and TNM stage. IGSHCC succeeded in distinguishing HCC patients with distinct OS durations in all of the subpopulations (Supplementary Fig. 4A–L). Again, these results demonstrated that IGSHCC is a promising predictive signature with high clinical relevance that can be used to further assess the prognostic risk in HCC patients on the basis of existing stratification systems.
IGSHCC outperforms previously published HCC prognostic signatures
The number of prognostic signatures for HCC reported in the literature is increasing. We chose representative signatures to compare their effectiveness with that of IGSHCC. The selection criteria were that the signature was identified in more than 100 HCC patients and validated in an independent cohort. The details of the selected studies are provided in Supplementary Table 3.19–21 Because ROC curves are useful for organizing classifiers and visualizing their performance, we performed ROC analyses to evaluate the predictive accuracy of these selected signatures and IGSHCC regarding OS within the TCGA dataset. As shown in Figure 4A, we drew an ROC curve for each signature according to their sensitivity and specificity in survival prediction. In each curve, we drew a 45-degree line in gray as the negative reference. IGSHCC had a greater AUC (0.67, 95% CI: 0.62–0.72) than those of the other prognostic signatures, indicating that it provided the most accurate prognosis prediction. Although the difference in AUC between the IGSHCC and TNM stage was not significant (0.67 vs. 0.65), the combination of IGSHCC and TNM stage had greater predictive efficacy than did the TNM stage alone (0.72 vs. 0.65; p=0.0046) and age alone (0.72 vs. 0.51; p<0.0001). These results implied that IGSHCC outperforms other published prognostic signatures and, when integrated with the TNM stage, increased the sensitivity and precision of the prediction of OS of HCC patients.
An IGSHCC-based nomogram is useful in clinical HCC management
Because IGSHCC proved to be a powerful prognostic signature model for HCC, we further attempted to provide clinicians with a quantitative method of predicting the 3- and 5-year OS rates in HCC patients. Therefore, we generated a nomogram by integrating IGSHCC with age and TNM stage to leverage the complementary value of these independent prognostic factors (Fig. 4B). Calibration curves were consistent with actual observations of OS duration and nomogram-based prediction of 5-year OS likelihood in the TCGA (Fig. 4C) and LCI validation (Fig. 4D) datasets. Similarly, regarding 3-year OS, the bias-corrected lines in calibration curves were still close to the ideal 45-degree line for both datasets (Supplementary Fig. 5A, B), indicating good predictive performance of our nomogram for HCC prognosis. In addition, we plotted ROC curves to evaluate the predictive accuracy of the nomogram. In the entire TCGA dataset, the AUCs for 3- and 5-year OS were both 0.72 (Fig. 4E). We further verified the predictive efficacy of IGSHCC in the LCI validation dataset, with AUCs of 0.71 and 0.69 for 3- and 5-year OS, respectively (Fig. 4F).
DCA can determine whether application of the nomogram to predicting cancer outcomes increases net clinical benefits when compared with the treat-none or treat-all scheme. The results of DCA in the CGA and LCI validation datasets showed that the nomogram provided HCC patients with greater net clinical benefits in predicting 5-year OS likelihood than did IGSHCC alone (Fig. 4G, H). The nomogram also exhibited greater prediction threshold probabilities than IGSHCC. Consistently, the DCA of 3-year OS also demonstrated increased net clinical benefits and threshold probabilities of the nomogram (Supplementary Fig. 5C, D). The results confirmed the clinical practicability of IGSHCC in that, when integrated into a nomogram, it can serve as a quantitative method for predicting OS rates in HCC patients.
IGSHCC can distinguish HCC patients with contrasting tumor immune microenvironments
An increasing number of patients with cancer have benefited from immune checkpoint blockade.7 Considering that we constructed the IGSHCC based on 22 immune-related genes, we wanted to determine whether the signature can provide some clues about tumor immune microenvironments, which would contribute to treatment responder identification and immunotherapy guidance among HCC patients. To that end, we applied ssGSEA to tumor samples from the CGA to deconvolve the relative abundance of 24 immunocyte types based on bulk tumor gene expression profiles. We visualized the abundance of immunocytes in each tumor sample using a heatmap and then categorized tumor immune infiltration into three subgroups (low, moderate, and high) using unsupervised hierarchical clustering (Fig. 5A). We also assessed the associations of immune cell infiltration with other clinical features, such as the TP53 mutation status. HCCs with high IGSHCC scores tended to exhibit relatively low tumor immune cell infiltration statuses and vice versa. This demonstrated that high-IGSHCC patients may have had “cold” tumors and were unlikely to respond to immunotherapy and that low IGSHCC patients may have benefited from immunotherapy. These results demonstrated that IGSHCC reveals, to some extent, the immune phenotypes for HCCs.
Furthermore, considering the concomitant infiltration of tumor promoting and suppressing immune cells in the tumor microenvironment, we compared the cytotoxic function of T cells in the low- and high-IGSHCC groups to see if the prognostic signature was correlated with the cytotoxic potential of immunocytes. As expected, the T-cell infiltration score (p<0.0001), cytotoxic cell infiltration score (p=0.0001), and DC infiltration score (p<0.0001) were increased in the low-IGSHCC group (Fig. 5B–D), patients in which tended to have high immunocyte infiltration statuses. In contrast, the macrophage infiltration score (p=0.0015) was higher in HCC patients with high IGSHCC scores than in those with low IGSHCC scores (Fig. 5E). The data demonstrated that HCC patients with low IGSHCC scores may have immune-activated phenotypes, characterized by increased immune cell infiltration and high cytotoxic potential, whereas those with high IGSHCC scores suffer from relatively low levels of immune activation.
Low- and high-IGSHCC groups have different sensitivities to immunotherapy
Because HCC immunotherapy response is closely related to immune activation in the local microenvironment, we then performed gene set variation analysis to explore the association between immune functional condition and IGSHCC score. As depicted in the heatmap in Figure 6A, distinct functional pathways were enriched in the low- and high-IGSHCC groups. Interleukin-6 signaling, interferon-α response, interferon-γ response, and T-cell activation were among the most significantly enriched biological processes in the low-IGSHCC group. Next, we characterized the low- and high-IGSHCC groups regarding the hallmark gene sets. Several tumorigenesis pathways, such as unfolded protein response, G2/M checkpoint, the phosphoinositide 3-kinase/Akt/mammalian target of rapamycin signaling pathway, the transforming growth factor-β signaling pathway, and DNA repair pathways, were enriched in the high-IGSHCC group (Fig. 6B). These results again showed that HCC patients in the low- and high-IGSHCC groups have distinct immune conditions in the tumor microenvironment.
To determine whether our prognostic signature can predict HCC patient responses to immunotherapy, we used the TIDE algorithm to calculate the probability of tumor response to immunotherapy in the TCGA dataset and LCI validation set. The four-fold contingency tables in Figure 6C, D showed that low-IGSHCC patients were more likely than high-IGSHCC patients to have immunotherapy-sensitive disease (TCGA dataset, p<0.001; LCI validation set, p=0.032). Moreover, we applied subclass mapping to compare the expression profiles for the HCCs in these two groups with a published dataset for 47 patients with melanoma that responded to anti-PD-1 and anti-CTLA4 immunotherapy.29 This demonstrated that compared with high-IGSHCC patients, low-IGSHCC patients were more likely to have responses to anti-PD-1 therapy (Bonferroni corrected p=0.008; Fig. 6E). The findings proved that our prognostic signature can effectively predict sensitivity of HCC to immunotherapy and thus contribute to identification of suitable candidates for anti-PD-1 therapy. IGSHCC score can serve as a useful reference for determining the immunotherapeutic strategy for HCC.
IGSHCC successfully stratifies HCC patients in an independent PCR array cohort
To validate the robust power of IGSHCC in prognostic risk-stratification and immunotherapy response prediction, we collected and analyzed data on an independent cohort of HCC patients in our institution. The detailed clinical characteristics of the HCC patients in this PCR array validation dataset (LCI validation set 2) are shown in Supplementary Table 1. Based on the formula described in Materials and Methods, the IGSHCC score for each individual was calculated, and patients were assigned to the high-IGSHCC (n=33) or low-IGSHCC (n=23) groups according to the cutoff value (Fig. 7A). Kaplan-Meier analysis demonstrated that IGSHCC retained the ability to predict prognosis for HCC patients in this cohort (HR: 2.49, 95% CI, 1.19–5.22; p=0.012; Fig. 7B).
To confirm the potential value of IGSHCC in identifying HCC patients who need immunotherapy, we examined the infiltration of CD8+ T cells and expression of PD-L1 in the tumors in patients who received nivolumab (anti-PD-1 antibody) in our PCR array validation dataset using immunofluorescence. Compared with the nonresponders, nivolumab responders tended to have greater CD8+ T-cell infiltration and lower expression of PD-L1 in tumors (Fig. 7C). Importantly, we found that HCC patients who responded well to treatment with nivolumab had lower IGSHCC scores than did nonresponders (Fig. 7D), which greatly supported our previous conclusion that the signature can help identify suitable HCC candidates for immunotherapy. Collectively, the observations highlighted the clinical significance of IGSHCC in predicting survival and immunotherapy response in HCC patients.
Discussion
The heterogeneity of HCC leads to greatly varied clinical outcomes among patients, and effective prognostic predictors are lacking. Furthermore, although the achievements of immunotherapy have revolutionized cancer treatment and provided significant survival benefits to patients,7,30 accumulating evidence from preclinical and clinical studies demonstrates that traditional immune-related markers such as PD-L1 are neither consistent nor reliable predictors of immunotherapy outcomes.31 In view of that, we designed this study to develop a signature based on immune-related genes that is associated with the tumor immune microenvironment and can provide an individualized prediction of both prognosis and immunotherapy sensitivity for HCC patients. A novel 22-gene signature named IGSHCC was developed, and proved to be a reliable tool to predict HCC survival and immunotherapy response, which was successfully validated in multiple public datasets, as well as an independent clinical cohort of HCC patients in our own institution. In addition, the IGSHCC-based nomogram had satisfying accuracy and sensitivity in predicting the likelihood of 3- and 5-year OS of HCC and provided added clinical benefits to HCC patients. The study results will serve to improve prognostic prediction and personalized therapy in HCC patients by measuring the expression of genes in the IGSHCC panel in tumors.
Based on large public datasets, researchers have developed several prognostic signatures for HCC,19–21 but they have yet to be incorporated into clinical practice because of limitations such as the overfitting risks arising from screening of a large gene pool. In addition, Cox proportional hazards regression analysis is unsuitable for high-dimensional microarray data in which the number of covariates is close to or greater than the number of observations.32 In this study, we used the LASSO Cox regression method, which can select the most valuable variables with robust prognostic value and low correlations to prevent overfitting.17 As expected, IGSHCC had better predictive performance than did other published prognostic models, and there was no reduction in stability and reliability of our signature when applied across independent datasets. In addition, there is no significant gene overlap between our signature and other published models, even if the same construction method was used. The main reason is that the distribution ratios we and other researchers used to randomly divided dataset into training set and internal verification set were different. Therefore, a unified modeling standard and procedure is warranted to promote the consistent use of predictive signatures in clinical practice.
Most of the genes in our prognostic signature were cytokines, chemokines, and their receptors that play key roles in immune and inflammatory response. Several of the genes, including STC2, TNF, CCR8, and BMP6, are associated with human cancers. For example, increased STC2 expression is correlated with aggressive disease and shortened OS in renal cell carcinoma patients.32 Also, BMP6 has been shown to inhibit mast cell recruitment in melanoma, which indicates the role of this gene in modulation of the tumor immune microenvironment.33 Some genes included in the signature may be therapeutic targets in HCC patients, but their biological functions have yet to be fully explored in HCC cases, which is well deserving of further investigation.
Some researchers have proposed that stratifying cancer patients based on the characteristics and quality of tumor immune infiltration is a critical next step in predicting the responsiveness of their disease to immunotherapy.34 In this study, the immune landscape of HCCs depicted by ssGSEA sheds light on the negative correlation between immune infiltration status and IGSHCC score. The complexity and diversity of tumor immune components drove us to explore the association of the IGSHCC score with T-cell activation and cytotoxic function. The decreased T-cell infiltration score, cytotoxic cell infiltration score, and DC infiltration score, and increased macrophage infiltration in patients with high IGSHCC scores demonstrated that their tumors are more likely than those in patients with low scores to be infiltrated-excluded and poorly immunogenic, with a lack of effective antigen-presenting cells and adaptive antitumor cytotoxic T cells. Moreover, several immune-related pathways were enriched in tumors in the low-IGSHCC group, whereas some tumorigenesis pathways were enriched in those in the high-IGSHCC group, further supporting the distinct immune biological processes in low- and high-IGSHCC tumors. In addition, given the close association between IGSHCC score and tumor immune phenotype, we performed an analysis using the TIDE algorithm and subclass mapping that demonstrated that HCC patients with low IGSHCC scores were more likely than those with high scores to benefit from anti-PD-1 immunotherapy.
Most important, we validated the robust prognostic risk-stratification ability and immunotherapy response prediction of IGSHCC in an independent cohort of HCC patients from our institution. External verification of the signature in a real-world cohort strongly demonstrated that prognosis risk and immunotherapy response in HCC patients can be readily evaluated by measuring the expression of 22 immune-related genes in tumor samples using PCR. Thus, our study, for the first time, shows the feasibility of translating IGSHCC into a clinical assay and set the groundwork for future prospective cohort studies to assess the overall translational potential of this signature.
The study had some limitations. First, the retrospective study design and single-center validation require future prospective validation of the results in multicenter cohorts. Second, the patients’ clinical information was incomplete, which impeded comprehensive analyses and decreased statistical reliability and validity. Moreover, the function of some genes in the signature in HCC cases is not clearly understood, so basic experimental studies are warranted to reveal the biological mechanisms underlying our prognostic signature. Finally, although we verified the predictive value of IGSHCC in HCC patients undergoing surgical resection, for major liver cancer patients who cannot receive surgery, whether the expression level of these genes in the tumor can be obtained through minimally invasive methods like needle biopsy and liquid biopsy needs further research to explore whether the signature can be applied to the majority of HCC cases.
In conclusion, IGSHCC is the first immune-related prognostic signature for HCC patients and is a robust and independent biomarker for predicting the risk of OS, even in specific clinical subgroups of patients. Furthermore, this immune-based signature provides a broad view of tumor immune phenotypes, resulting in improved identification of HCC patients who are likely to benefit from immunotherapy. The clinical feasibility of IGSHCC and the IGSHCC-based nomogram should be studied further prospectively in multicenter cohorts.
Supporting information
Supplementary Fig. 1
Development of the immune-related gene signature for HCC.
(A) Partial likelihood deviance of covariates in 78 OS-related candidate genes. (B) LASSO coefficient profiles of the 22 immune-related genes. The coefficient value was chosen using 10-fold cross-validation. (C) The clustered immune function of the 22 immune-related genes included in IGSHCC. (D) The coefficient of each gene in IGSHCC ranked from high to low value. The yellow bars represent the positive coefficients, whereas the blue ones represent the negative coefficients.
(DOCX)
Supplementary Fig. 2
Distribution of IGSHCC scores in HCC patients.
(A) TCGA training data set, (B) TCGA validation set, (C) entire TCGA data set, and (D) LCI validation data set. In each figure, the top panel represents the low- and high-IGSHCC groups defined based on the cutoff value for TCGA training data set, the middle panel represents the survival statuses and durations in the HCC patients in the low- and high-IGSHCC groups, and the bottom panel is a heatmap of the expression of the 22 immune-related genes in IGSHCC in each tumor.
(DOCX)
Supplementary Fig. 3
Cox regression analysis of OS in the HCC patients in the entire TCGA dataset based on the IGSHCC score and clinicopathological risk factors, including age, sex, AFP level, tumor grade, and TNM stage.
(A) Univariate and (B) multivariate. The forest plots show the HRs for death (solid blue squares) and their 95% CIs (horizontal red lines) for each variable.
JCTH-D-21-00283-s003.docx
(DOCX)
Supplementary Fig. 4
Kaplan-Meier analyses of OS in the HCC patients in subgroups with different molecular and clinical characteristics.
Patients in the entire TCGA dataset were stratified into subgroups including (A and B) TP53 genotype, (C and D) CTNNB1 genotype, (E and F) AFP level, (G and H) age, (I and J) tumor grade, and (K and L) TNM stage. The survival difference between the low- and high-IGSHCC patients in these subgroups was analyzed using the log-rank test.
JCTH-D-21-00283-s004.docx
(DOCX)
Supplementary Fig. 5
Evaluation of the predictive accuracy and net benefit of IGSHCC-based nomogram for 3-year OS.
(A and B) Calibration plots for the IGSHCC-based nomogram used to evaluate the consistency of predicted and actual 3-year OS durations in the entire TCGA dataset and LCI validation data set. The 45-degree line represents ideal performance of the nomogram. (C and D) DCA of the nomogram (red line) and IGSHCC (blue line) in the TCGA and LCI validation dataset. The black line at the bottom and gray line represents the assumptions that none and all HCC patients survived at 3 years, respectively. The x-axis represents the threshold probability of death, and the y-axis represents the corresponding net benefit of the model.
JCTH-D-21-00283-s005.docx
(DOCX)
Supplementary Table 1
Clinical characteristics of the HCC patients in the four data sets.
(DOCX)
Supplementary Table 2
Details of the 22 immune-related genes in IGSHCC.
(DOCX)
Supplementary Table 3
Summary of the three reported HCC prognostic signatures.
(DOCX)
Supplementary Table 4
Primer sequences for the 22 immune-related genes in IGSHCC used in the PCR array.
(DOCX)
Supplementary Table 5
Seventy-eight candidate immune-related genes identified by univariate Cox regression analysis.
(DOCX)
Abbreviations
- AFP:
α-fetoprotein
- AUC:
area under the curve
- BCLC:
Barcelona Clinic Liver Cancer
- CTLA4:
cytotoxic T-lymphocyte-associate protein 4
- DCA:
decision curve analysis
- DC:
dendritic cell
- GEO:
Gene Expression Omnibus
- GZMA:
granzyme A
- HCC:
hepatocellular carcinoma
- LASSO:
least absolute shrinkage and selection operator
- MVI:
microvascular invasion
- OS:
overall survival
- PCR:
polymerase chain reaction
- PD-1:
programmed death 1
- PD-L1:
programmed death-ligand 1
- PRF1:
perforin
- ROC:
receiver operating characteristic
- ssGSEA:
single-sample gene set enrichment analysis
- TCGA:
The Cancer Genome Atlas
- TIDE:
tumor immune dysfunction and exclusion
- TNM:
tumor-nodes-metastasis
Declarations
Acknowledgement
We thank the China Scholarship Council for supporting Dr. Chenhao Zhou and Dr. Yuan Gao as a visiting PhD candidate at The University of Texas MD Anderson Cancer Center. We also thank Donald R. Norwood at the Research Medical Library at the University of Texas MD Anderson Cancer Center for editing this manuscript.
Data sharing statement
All data in this study are available from the Cancer Genome Atlas website (http://cancergenome.nih.gov/) and series GSE14520 in the GEO from the Liver Cancer Institute (LCI) at Zhongshan Hospital, Fudan University (https://www.ncbi.nlm.nih.gov/geo/). Other data analyzed in this study are included in this article and its supplementary materials.
Funding
The work was supported by the National Natural Science Foundation of China (82103521, 82073208); the Shanghai Sailing Program (21YF1407500); the China Postdoctoral Science Foundation (2021M690674); the Special Foundation for Science and Technology Basic Research Program (2019FY101103); the Shanghai Shen Kang Hospital Development Center new frontier technology joint project (SHDC12021109).
Conflict of interest
QD has been an editorial board member of Journal of Clinical and Translational Hepatology since 2021. The other authors have no conflicts of interest related to this publication.
Authors’ contributions
Study concept, design, and supervision (NR), analysis and interpretation of data (CZ, JW), drafting of the manuscript (CZ, JW, YG), acquisition of data (CL, XZ, QZ, JS), preparation of figures and tables (CWL, MA), critical revision of the manuscript for important intellectual content (MCH, YS, QD, YL), provision of patient tissue samples (YY, QY).