Introduction
Hepatocellular carcinoma (HCC) has become the second leading cause of cancer-associated death, which affects patients worldwide, and is associated with early recurrence and a poor response to chemotherapy (1,2). As our understanding of the roles of immune checkpoints in tumor cells and the surrounding non-tumor cells in the cancer microenvironment has advanced, novel technologies, such as Chimeric antigen receptor T-cell therapy and PD-1/PD-L1 checkpoint inhibition therapy, have been developed to target the immune environment of HCC to improve the prognosis of patients following HCC resection (3-6). However, the overall responses rates of patients treated with specific checkpoint blockers in HCC, such as those targeting PD-1 or CTLA-4, are not favorable, possibly due to unknown changes to the immune microenvironment (7).
The emergence of high-throughput nucleotide sequencing analysis provides new perspectives to understand the genomic changes in tumors, revealing the differentially expressed genetic signatures between tumor tissues and normal tissues (8,9). Several studies in different types of cancer, including breast cancer, thyroid cancer, non-squamous non-small cell lung cancer and colorectal carcinoma, have examined differences in the patterns of immune signatures to improve our understanding of the immune environment and the mechanisms underlying tumor development and progression (10-14). However, the specific immune genetic changes in HCC have not been extensively studied, although one study found that the levels of immune cells infiltration, specifically by T-cells, cytotoxic cells, Th2 cells and macrophages, in HCC were associated with improved survival in patients based on n silico analysis, suggesting that the type of immune cells present in HCC tissues were different from the immune cell profile of the normal liver (15). As the liver is now considered to a 'immune associated organ', the presence of immune cells in HCC should be taken into consideration as a leading factor for predicting prognosis following resection, and should not be restricted to specific types of immune cells (16-18).
In the present study, the changes in expression of immune related genes in HCC tissues were compared with the adjacent healthy matching tissues, using bioinformatics analysis. The immune-associated genes identified was derived from a comprehensive study of the immune landscape of 20 solid tumors, which allowed for evaluation of relevant immune functions and the immune status of solid tumors in a simplified manner (19,20). The aim of the present study was to identify the immune-related genetic changes between HCC tissues and normal liver tissues, to understand the effects of immune regulation of HCC, and the effect on progression of HCC. Additionally, an immune evaluating model for prognostic evaluation in HCC patients was constructed, with the aim of differentiating patients into sub-populations for more personalized clinical treatment to maximize the efficacy of therapies used, particularly for treatment with immune checkpoint inhibition.
Materials and methods
Datasets
Data on patients with HCC were obtained from The Cancer Genome Atlas (TCGA; cancer.gov/tcga) and ICGC (icgc.org), which are publicly available databases (21,22). The databases contained information on 370 (TCGA) and 232 cases (ICGC) of HCC, which included RNA sequencing information and the clinical characteristics (Table I). In the data obtained from TCGA, there were 249 men and 121 women with a median age of 61 (range, 16-85). In the ICGC dataset, there were 171 men and 61 women with a median age of 69 (range, 31-89). The list of immune-related genes for analysis was obtained from previous studies (19,20) which contained a total of 821 immune related genes.
Differential expression analysis
DEIGs between adjacent and HCC tissues were analyzed using the limma package on the cohort from TCGA (23). The raw data were normalized and log2(x+1) transformed. Genes with a fold change >1 and an adjusted P-value <0.05 were considered significant (based on false discovery rates using the Benjamini-Hochberg approach) (24). A heatmap of significantly up or downregu-lated immune-associated genes was plotted using the heatmap package version 1.8.0 (git.bioconductor.org/packages/heat-maps), and these genes were used for further prognostic analysis.
Gene ontology annotation and pathway enrichment
Immune genes determined to be significantly differentially expressed were functionally annotated using the clusterProfiler package (25), which stratified pathways according to one of the following categories: Cellular compartment, biological process or molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to enrich the pathways associated with the identified genes (26).
Construction of the immune risk score (IRS) model
The entire cohort of patients with HCC from TCGA were randomly divided into a training set and a testing set to construct and assess the prognostic model. The DEIGs were evaluated using a univariate Cox model for individual risk factors affecting OS of the training set (P<0.05), and the associated genes were analyzed together using a Lasso penalty linear regression model, which were subsequently used to construct a multivariate Cox model. In Lasso regression, the patients were subsampled 1,000 times, and the genes with an occurrence >900 instances were selected. In multivariate Cox regression analysis, a stepwise method is used, where all combinations of the identified genes are assessed to construct the best combination of the immune associated gene set. The prognostic value of the linear IRS model was validated in the testing set, the entire TCGA cohort and the independent cohort from the ICGC database separately, with patients divided into high- and low- risk sub-populations according to the median IRS. Kaplan-Meier survival curves and time-dependent receiver operating characteristic curves (ROC) were used to demonstrate the prognostic value of the 5-gene IRS model, using the R packages of survival (rdocumentation.org/packages/survival) and survivalROC (27).
Independent prognostic value of IRS and the construction of a nomogram
The independent prognostic value of IRS was further examined through univariate and multivariate Cox regression analysis in combination with the clinical characteristics, such as age, sex, tumor grade and tumor stages. Following evaluation of the risk effect of clinical characteristics and IRS, a nomogram model was constructed for prognostic prediction, which included the IRS and tumor stage. The predictive value of the nomogram was further confirmed using ROC curves for prediction of the 1, 3 and 5 year OS rates, in which the predictive value of the single risk factors were also assessed independently. The C-index of the nomogram was calculated with a bootstrap of 1,000 resamples, and the results ranged between 0.5-1.0, where 0.5 indicated a random chance and 1.0 indicated perfect separation of the outcomes. Calibration curves were also plotted to demonstrate the precision of the nomogram, contrasting the predictive probability with the actual incidence.
Tumor infiltrating immune cells and their correlation with IRS
The calculation of tumor infiltrating immune cells in patients with HCC from TCGA was performed using Tumor IMmune Estimation Resource (TIMER), an online tool which contains the reanalyzed genomic expression data across 32 types of cancer, including over 10,897 samples from TCGA (28,29). The online portal calculates the abundance of 6 types of infiltrating immune cells; B cells, CD4+ T cells, CD8+ T cells, neutrophils, dendritic cells and macrophages. The abundance of infiltrating immune cells was correlated with IRS, and significance was examined using a Pearson's correlation test. P<0.05 was considered to indicate a statistically significant difference. The degree of correlation between immune cell abundance and IRS was defined as follows: Very low, 0.0-0.2; low, 0.2-0.4; medium, 0.4-0.6; high, 0.6-0.8; and very high, 0.8-1.0.
Statistical analysis
Statistical analysis was performed in R (version 3.6.1), using R studio (version 1.2.1335) (30,31). DEIGs between adjacent and HCC tissues were analyzed using a Wilcoxon Signed-rank test. Univariate Cox regression, Lasso regression and multivariate Cox regression analysis were performed to construct the IRS model. The infiltration levels of different immune cells between HCC and para-tumor tissues were compared using Pearson's correlation coefficients. P<0.05 was considered to indicate a statistically significant difference.
Results
Analysis strategy and overview of the DEIGs
Analysis of the TCGA dataset identified 247 significantly differently expressed genes, of which 200 were upregulated and 47 were downregulated (Fig. 1A and B). Gene ontology of the DEIGs were primarily associated with immune cell activation, adhesion or responses to stimuli (Fig. 1C and E). The enriched pathways for those DEIGs were primarily enriched in cytokine and cytokine receptor interactions between cells, in which the cytokine-cytokine receptor interaction pathways had the highest counts of associated genes and significance. Additionally, the z-score of the cytokine-cytokine receptor interaction pathways was the highest ranked amongst all enriched pathways (Fig. 1D and F).
Relative risk effect of the immune associated genes and construction of the IRS based
Patients from TCGA were divided into two sets; a training set and an internal test set. In the training set, significant DEIGs (P<0.05) associated with OS in the Cox model (Fig. 2A), were further analyzed using a Lasso penalty linear regression model (Fig. 2B). The final multivariate Cox model was constructed using 5 genes [Enhancer of zest homology 2 (EZH2), ferritin light chain (FTL), complement factor H related 3 (CFHR3), isthmin 2 (ISM2), cyclin dependent kinase 5 (CDK5)], of which EZH2 and ISM2 still significant and had high risk effect following adjustment (Fig. 2C).
The final model was then used to calculate the IRS of patients for prognostic evaluation in four separate cohorts: The training cohort from TCGA (Fig. 3A-C), the test cohort from TCGA (Fig. 3D-F), the entire TCGA cohort (Fig. 4A-C) and the independent cohort from ICGC (Fig. 4D-F). The median IRSs of the four cohorts were used to stratify patients into high- and low-score groups.
Survival validation of the IRS model
The IRS model, based on the 5 significantly changed immune associated genes, was able to divide the patients with HCC into high- and low-risk sub-groups based on the corresponding score levels. Patients with higher scores had a worse prognosis in all four cohorts (Fig. 5A-D). The area under the curve (AUC) values of the model for 0.5, 1, 2, 3 and 5 year survival in all four cohorts were all ~0.7, with a lowest value of 0.643 for 5 year survival in the test cohort (Fig. 5E-H) Further analysis of patients with different tumor stages and grades, showed that a high IRS predicted worse survival in both datasets from TCGA and ICGC (P<0.05), demonstrating the independent prognostic value of IRS for clinical use. In patients with stage I & II, III & IV and grade 1 & 2 cancer in TCGA, the curves for high- and low-IRS showed notable differences in the 6 year survival, whereas after 6 years of follow-up, the curves nearly overlapped (Fig. 6A-C). In patients with grade 3 & 4 cancer from TCGA, and stage I & II, and III & IV patients from ICGC, the curves of patients with high- and low-IRS diverted with no overlap (Fig. 6D-F). Regarding disease free survival, patients with a high IRS also exhibited worse outcomes compared with patients with a low IRS, and this difference was significant in the entire TCGA cohort, in patients with stage I & II, and grade 1 & 2 cancer from TCGA (Fig. S1).
Nomogram of IRS and other associated clinical factors
To assess the clinical relevance and significance of the IRS model, the IRS model was combined with the clinical characteristics for prognostic prediction of the data obtained from TCGA and ICGC. In univariate and multivariate analysis of patients from TCGA, cancer stage and IRS were significantly associated with survival, with or without adjustment (Fig. 7A and B). In the ICGC dataset, in addition to cancer stage and IRS, the presence of a previous malignancy was also significantly associated with survival following adjustment (Fig. 7C and D). The correlations between the 5-gene model and clinical characteristics are presented in Fig. 7E and F.
For clinical use, a nomogram was constructed for all the significant factors, including clinical tumor stage and IRS, using the entire TCGA cohort. (Fig. 8A) The AUCs for 1, 3 and 5 year survival were higher in the nomogram compared with IRS or tumor stage (Fig. 8C) The C-index for the nomogram-predicted OS was 0.749, with 1,000 cycles of bootstrapping (Table II). Calibration graphs were drawn to evaluate the corresponding performance of the nomogram for predicting 1, 3 and 5 year OS, and the lines almost overlapped, suggesting its accuracy (Fig. 8B). These results all show the value of the nomogram for predicting OS in patients following resection, and was shown to be more accurate than tumor stage or IRS, for both short- and long-term.
Potential roles of immune infiltrating cells in prognostic prediction
Using the evaluation scores of six common types of immune infiltrating cells from TIMER, the abundance of immune-cells in HCC tissue was estimated. Patients in TCGA were used to calculate the abundance of infiltration, and a score was generated by the tool. The correlation between the scores of B cells, CD8+ T cells, dendritic cells, macrophages, neutrophils and IRS were significant and positive, suggesting an association between increased infiltration and IRS (Fig. 9). However, five types of infiltrating immune cells had very low correlation coefficients with IRS, and macrophages had a low correlation, suggesting that the IRS based on the five-gene model was primarily dependent on change sin expression of immune-associated genes in HCC tissues, as opposed to infiltration of immune cells.
Discussion
Malignant HCC is associated with poor outcomes and with high recurrence rates following resection (1,2). In addition, a growing body of evidence highlights the vital value of immune regulation in HCC, and poor responses to treatment with chemotherapy highlight the need for drugs with greater specificity for immune targets (3-6). The development of immune checkpoint blockers, such as nivolumab and pembrolizumab have not yielded optimistic results for patients with HCC, and this may be associated with the immune microenvironment of HCC tissues, as patients with a high degree of immune cell infiltration often exhibit more favorable outcomes (3). Understanding the immune environment and immune status of HCC may result in improved strategies for treatment of patients, resulting in improved prognosis.
In the present study, the differentially expressed immune associated genes between HCC and adjacent tissues were identified, and following Lasso regression and multivariate Cox analysis, 5 significantly differentially expressed OS-related immune associated genes (EZH2, FTL, CFHR3, ISM2, CDK5) were used to build a prognostic nomogram in combination with clinical characteristics. The prognostic values of the nomogram for 1, 3 and 5 year survival achieved was >0.8, and performed better compared with the individual clinical risk factors, and may thus be used to stratify patients with HCC in clinical practice, preventing early recurrence.
EZH2 is expression is low in the normal liver, and is associated with methylation of histone H3K27 and recruitment of methyltransferases, which are involved in DNA replication for cancer progression, and stem cell maintenance and differentiation of other cell lineages, such as immune cells (32,33). Several studies have confirmed its prognostic value and importance in various types of cancer, including lymphoma, glioma, head and neck carcinoma, and cervical neoplasia (34-38). Novel therapeutic methods have been developed to target EZH2, although high expression of EZH2 is not always correlated with the malignancy of a cancer, such as in colorectal cancer (39-41). In the present study, it was shown that high expression of EZH2 was significantly associated with poor OS in patients in the univariate analysis, including after adjustment, highlighting the prognostic value of high EZH2 expression in patients with HCC. EZH2 may promote the development and proliferation of HCC, and may thus result in recurrence following tumor resection.
Ferritin light chain (FTL) is the light subunit of the ferritin protein, which is involved in iron release and uptake in tissues. Several allergic diseases, inflammation status and oxidative damage are associated with the roles of FTL, including systemic lupus erythematosus, cataracts, hepatitis E virus infection and Alzheimer's disease (42-47). FTL also serves a role in cancer, where the dysfunction of iron metabolism is a hallmark of various types of cancer, and is involved drug resistance and malignant progression (48-50). Although the risk effect of FTL was small following adjustment in the multivariate Cox model, the potential role of changes to iron metabolism in the progression of HCC should not be ignored.
CFHR3 was demonstrated to exert a protective effect in patients with HCC, and physiologically, CFHR3 is exclusively expressed in normal liver (51). CFHR3 is associated with compliment factor H, which can bind to the C3d region of C3b, regulating the function of compliment system (52). The loss of CFHR3 results in age-related macular degeneration, and high expression levels of CFHR3 may result in atypical hemolytic-uremic syndrome (53-54). There are no published articles regarding the expression of CFHR3 in cancer, and the protective role of CFHR3 observed in the present study may highlight a potential change in the expression profile that may be used for improving the prognosis of patients with HCC.
In both univariate and multivariate analysis, ISM2 was significantly associated with poor outcomes and was considered a high risk factor in patients with HCC. ISM2 is a component of thrombospondin (THBS), which promotes the activity of mesenchymal and stromal cells through TGF-β, and regulates secretion of inflammatory cytokines through the NF-κB signaling pathway (55,56). THBS promotes epithelial-to-mesenchymal transition in melanoma, and exacerbates the progression of prostate cancer to more advanced stages (57,58). Additionally, THBS may serve a role in gastric carcinogenesis, and invasiveness of breast cancer and nodal metastasis (59-62). The role of ISM2 or THBS in HCC has not been explored to the best of our knowledge, and further analysis is required to understand their potential roles and effects.
CDK5 has been extensively studied as an important factor in tumor development and metastasis (63-65). Although CDK5 shares homologous structure with other CDKs, it is not cyclin-dependent and does not need to be phosphorylated in the T-loop for activation (66). CDK5 expression is upregulated in several types of cancer, and inhibition induces cancer cell death through a FOXO1-Bim pathway or mitochondrial dysfunction (67-73). Ehrlich et al (74) showed that expression of CDK5 was increased in HCC tissues, and was correlated with malignant phenotypes. Additionally, CDK5 was most active in the G2/M phase of cancer cells in the nucleus, and regulated DNA damage response through phosphorylation of ataxia telangiectasia mutated kinase, validating the prognostic role of CDK5 in the present study.
Infiltration of immune cells in HCC tissues was assessed, and 5 of the 6 common types of immune cells were significantly associated with IRS; however the correlation coefficients for all 6 types of cells were either low or very low. Thus, although patients with a high degree of immune cell infiltration may have a high IRS, the IRS based on the five immune associated genes primarily accounted for the functional status of the microenvironment in HCC tissues. Recent studies have focused on the roles of infiltrating immune cells in the tumor microenvironment, to explain the mechanisms underlying immune evasion and to predict drug response or prognosis (75-77). Further analysis of the sub-types of immune cells in HCC may improve our understanding of immune function in tumor, and with advances in technologies, changing the types of immune cells present may be considered as a potential treatment strategy, emphasizing the importance of restoring immune function in HCC (78-80).
The IRS model and the nomogram developed in the present study may exhibit value in clinical practice for prognostic prediction. Based on the IRS model and nomogram, it may be possible to tailor therapeutic regimens to each specific patient, or they may be useful for predicting/detecting early recurrence, and to evaluate immune function in HCC to optimize the benefits of monoclonal targeting therapies.
The present study has some limitations. First, the cohorts from TCGA and ICGA are primarily from several local populations, and thus may not be applicable to all ethnicities. Second, the DEIGs in this analysis may not reveal the holistic changes in the immune microenvironment in HCC. Third, the present study focused primarily on overall survival of patients following tumor resection, and disease-free survival was not assessed as this information was not contained in the data-sets. Furthermore, experimental validation of the prognostic signatures in HCC cell lines and human tissues is required to validate their relevance and improve our understanding of their respective roles, and will be performed in future experiments.
In conclusion, a 5-gene model was constructed from differentially expressed immune associated genes to evaluate the IRS of patients for independent prognostic prediction. By combining IRS with clinical tumor stage, a nomogram was constructed with efficient predictive value for 1, 3 and 5 year OS for patients with HCC. This nomogram may be used clinically for monitoring early recurrence and prognostic prediction.