Immune system‑associated genes increase malignant progression and can be used to predict clinical outcome in patients with hepatocellular carcinoma.

Hepatocellular carcinoma (HCC) is one of the most malignant types of cancer, and is associated with high recurrence rates and a poor response to chemotherapy. Immune signatures in the microenvironment of HCC have not been well explored systematically. The aim of the present study was to identify prognostic immune signatures and build a nomogram for use in clinical evaluation. Using bioinformatics analysis, RNA‑seq data and overall survival (OS) information on 370 HCC cases from TCGA and 232 HCC cases from ICGC were analyzed. The differential expression of select immune genes, based on previously published studies, between HCC and adjacent tissue were analyzed using the limma package in R. Enrichment of pathways and gene ontology analysis was performed using clusterProfiler. Subsequently, univariate Cox regression analysis, Lasso penalty linear regression and multivariate Cox regression models were used to construct a model for immune risk score (IRS). The R packages, survival and survivalROC, were used to plot survival and the associated receiver operating characteristic curves. Infiltration of immune cells was calculated using Tumor IMmune Estimation Resource, with significance examined using a Pearson's correlation test. P<0.05 was considered significant. Based on the analysis, expression of 200 immune genes were upregulated and 47 immune genes were downregulated immune genes. In the multivariate Cox model, 5 genes (enhancer of zest homology 2, ferritin light chain, complement factor H related 3, isthmin 2, cyclin dependent kinase 5) were used to generate the IRS. By stratifying according to the median IRS, it was shown that patients with a high IRS had poor OS rates after 1, 2, 3 and 5 years, and this result was consistent across the testing, training and independent validation cohorts. Additionally, the IRS was correlated with the abundance of infiltrating immune cells. The nomogram built using IRS and clinical characteristics, was able to predict 1, 3 and 5 year OS with area under the curve values of >0.8. These results suggest that the model developed to calculate the IRS may be used to monitor the effectiveness of treatment strategies and for prognostic prediction.

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.

Table I

Clinicopathological characteristics of patients in TCGA and ICGC.

TCGA. n=370
Clinical characteristics	n	%
Survival status
Survived	244	6555
Died	126	34.05
Age
≤65 years	232	62.70
>65 years	138	37.30
Sex
Male	249	67.30
Female	121	32.70
Histological grade
G1	55	14.86
G2	177	47.84
G3	121	32.70
G4	12	3.24
Stage
I	171	46.22
II	85	22.57
III	85	22.57
IV	5	1.35
T classification
T1	181	4852
T2	93	25.14
T3	80	21.62
T4	13	3.51
TX	1	0.27
M classification
M0	266	71.89
Ml	4	1.08
MX	100	27.03
N classification
N0	252	68.11
NI	4	1.08
NX	113	30.54
ICGC, n=232
Clinical characteristics	n	%
Survival status
Survived	189	81.47
Died	43	18.53
Age
≤65 years	90	38.79
>65 years	142	61.21
Sex
Male	171	73.71
Female	61	26.29
Clinical characteristics	n	%
Stage
I	36	15.52
II	106	45.69
III	71	30.60
IV	19	8.19
Prior malignancy
No	202	87.07
Yes	30	12.93

TCGA, The Cancer Genome Atlas; ICGC, International Cancer Gene Consortium.

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).

Figure 1

Differentially expressed immune genes between HCC tissues and matching normal tissues. (A) Differentially expressed genes. Green, downregu-lated; red, upregulated. Genes with a logFC>1 and a P-value <0.05 were considered significant. (B) Heatmap of 50 top-ranked immune genes between tumor tissue and normal tissue. (C) Gene ontology of the significantly differentially expressed genes. (D) Enriched pathways for the differentially expressed immune associated genes. (E) Ranks of gene ontology for the significantly differentially expressed immune associated genes. Significantly enriched gene ontologies are presented in ranked order based on the number of genes and the Z score. (F) Ranks of enriched pathways for the significantly differentially expressed immune associated genes. Significantly enriched pathways are presented in ranked order based on the number of genes and the Z score. HCC, hepatocellular carcinoma; FC, fold change; fdr, false discovery rate.

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).

Figure 2

Construction of immune risk score models using the significantly differentially expressed immune associated genes. (A) Univariate analysis to determine the risk effects of each immune associated gene alone in the training set. (B) Lasso penalty regression analysis using 8 significant risk genes. (C) Multivariate analysis of core-related immune associated genes. EZH2, enhancer of zest homology 2; FTL, ferritin light chain; CFHR3, complement factor H related 3; ISM2, isthmin 2; CDK5, cyclin dependent kinase 5.

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.

Figure 3

Stratification of patients in the training and test sets according to the median immune risk score. (A) Distribution of risk score, (B) survival overview and (C) mRNA expression levels in the internal training set. (D) Distribution of risk score, (E) survival overview and (F) mRNA expression levels in the internal test set. The dotted lines represent the median risk score cut-off dividing patients into low- and high-risk groups. The red dots and lines represent the patients in the high-risk group. The green dots and lines represent the patients in the low-risk group.

Figure 4

Stratification of patients based on the median immune risk score. (A) Distribution of risk score, (B) survival overview and (C) mRNA expression levels in the entire The Cancer Genome Atlas cohort. (D) Distribution of risk score, (E) survival overview and (F) mRNA expression levels in the International Cancer Gene Consortium cohort. The dotted lines represent the median risk score cut-off dividing patients into low- and high-risk groups. The red dots and lines represent the patients in the high-risk group. The green dots and lines represent the patients in the low-risk group.

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).

Figure 5

Survival curves of patients with different IRS and the predictive value of the IRS for 0.5, 1, 2, 3 and 5 year survival. Survival curves of patients stratified by the IRS in the (A) training set, (B) testing set, (C) entire TCGA cohort and the (D) external International Cancer Gene Consortium cohort. Survival was compared using a log-rank test; P<0.05 was considered to indicate a statistically significant difference. ROC curves for 0.5, 1, 2, 3 and 5 year survival in the (E) training set, (F) testing set, (G) entire TCGA cohort and (H) ICGC cohort. IRS, immune risk score; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Gene Consortium; ROC, receiver operating characteristic curve; AUC, area under the curve; CI, confidence interval.

Figure 6

Survival curves of patients stratified by IRS based on stage or grade. Differences in survival in patients with hepatocellular carcinoma. Patients were stratified according to IRS, and the survival of patients with (A) stage I&II, (B) stage III&IV, (C) grade 1 & 2 and (D) grade 3&4 cancer was compared in the patients from TCGA. (E) Survival of patients from ICGC stratified according IRS. Differences in survival with (E) stage I & II or (F) stage III & IV. Survival was compared using a log-rank test; P<0.05 was considered to indicate a statistically significant difference. IRS, immune risk score; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Gene Consortium; CI, confidence interval.

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.

Figure 7

Risk factors for overall survival in patients from TCGA and ICGC. (A) Univariate and (B) multivariate analysis of risk factors for overall survival in the data obtained from TCGA. (C) Univariate and (D) multivariate analysis of risk factors for overall survival in the data obtained from ICGC. Expression of five significant immune associated genes in correlation with clinical characteristics in (E) TCGA and (F) ICGC cohorts. TCGA, The Cancer Genome Atlas; ICGC, International Cancer Gene Consortium.

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.

Figure 8

Combined nomogram for IRS and clinical characteristics of 186 patients with hepatocellular carcinoma from TCGA. (A) Nomogram of clinical characteristics and IRS. (B) Validation of the nomogram for 1, 3 and 5 year survival. (C) ROC curves based on the nomogram, IRS or tumor stage alone for 1, 3 and 5 year survival prediction. IRS, immune risk score; ROC, receiver operating characteristic curve; AUC, area under the curve; OS, overall survival; FPR, False Positive Rate; TPR, True Positive Rate.

Table II

C-index analysis of models.

Model	C-index
Stage model	0.654
Prognostic model	0.746
Nomogram model	0.749

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.

Figure 9

Relationships between IRS and infiltrated immune cells in the hepatocellular carcinoma microenvironment. Comparisons were performed using a Pearson's correlation coefficient; P<0.05 was considered to indicate a statistically significant difference. Correlation between IRS and infiltration of (A) B cells, (B) CD4+ T cells, (C) CD8+ T cells, (D) dendritic cells, (E) macrophages, and (F) neutrophils. Survival was compared using a log-rank test; P<0.05 was considered to indicate a statistically significant difference. Cor, correlation coefficient.

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.