Mining Prognostic Biomarkers of Hepatocellular Carcinoma Based on Immune-Associated Genes

Abstract

This research aims to investigate the immune-associated gene signature from databases to improve the prognostic value in hepatocellular carcinoma (HCC) by multidimensional methods using various bioinformatic methods. Fifty-one immune-associated genes were mined out, which were associated with clinical characters through univariate and multivariate Cox regression analyses, and 51 immune-associated genes could be well-divided HCC samples into high-risk and low-risk clusters. Next, we performed least absolute shrinkage and selection operator (LASSO) Cox regression method to reveal 18 immune-associated genes' signature and calculate risk score of each gene for receiver operating characteristic (ROC) analysis. Comparing with low-risk cluster, high-risk cluster had higher risk score with unfavorable prognosis. Then, multivariate Cox regression analysis showed that risk score of 18 immune-associated genes' signature was associated with tumor invasion and tumor-node-metastasis (TNM) stage. ROC analysis indicated combined TNM stage, and risk score performed more sensitive and specific than single TNM stage or risk score in survival prediction. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis found that the pathways enriched in tumorigenesis were related to risk score, and those pathways could separate HCC samples into high and low clusters. In addition, the survival prediction of 18 immune-associated genes' signature was well validated in independent test data set, external data set, and Real-time Quantitative PCR (RT-qPCR) experiment. The 18 immune-associated genes' signature was constructed, which could be used in effective prediction of HCC prognosis.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors globally. In China, its morbidity and mortality are ranked first among all types of malignant tumors (Chen et al., 2016). HCC is characterized by high incidences, rapid progress, strong invasiveness, and poor prognosis. Although the diagnostic and therapeutic techniques for HCC are developing rapidly, the benefits to patients are limited (Dimitroulis et al., 2017; Hartke et al., 2017). This is because early clinical manifestations of this disease are not obvious and most of the patients in the middle or late stages of their first diagnosis lack the best opportunity for treatment. Therefore, the search for early molecular diagnostic markers of HCC and their use for accurate molecular classification are beneficial for individualized and accurate diagnosis and treatment of patients.

Pathological staging and clinical grading are still the most widely used guidance in clinical diagnosis and treatment of HCC despite their inadequate capabilities (Bruix and Sherman, 2011). However, molecular marker-based HCC phenotype and diagnosis markers, such as gene, DNA methylation, noncoding RNA, intestinal flora, exocrine among many more, are of great interest. In the current study, knowledge of HCC biology was translated into clinical practice to ensure that the HCC phenotype was closely related to gene mutation and transcriptomic classification (Calderaro et al., 2017). Tan et al. (2019) identified a potential therapeutic target and prognostic biomarker of the transcription factor, GATA6, which is capable of transforming the phenotypic reprogramming of HCC cells, is associated with metabolism and boost tumorigenesis, self-renewal, and metastasis (Tan et al., 2019). MiR-345 a novel biomarker in diagnosis and treatment is useful in HCC prognosis due to its function in the inhibition of tumor metastasis and epithelial–mesenchymal transition (Yu et al., 2017). Wang et al. (2019) developed a novel four-methylated mRNA signature-based risk score system, which was used in predicting prognosis and a useful guide for individualized treatment in HCC patients.

Recent studies have shown that immunological mechanisms play an important role in triggering various physiological and pathophysiological processes in tumors (Chen and Mellman, 2013). The chronic inflammatory immune response of the liver is the primary risk factor for tumorigenesis, and it is closely related to the development, growth, and treatment of the immune mechanism on HCC. In-depth research on how immune-associated genes promote HCC growth has led to the development of more effective therapeutic options (Makarova-Rusher et al., 2015). However, there is very limited research focused on the gene of immune phenotype in the HCC and its association with diagnosis, treatment, and prognosis. The aim of this study was to establish a genetic signature based on immune-associated genes to predict the prognosis of HCC. This was achieved through the multistep screening using various bioinformatic tools. The 18 immune-associated genes' signature for use in the prediction of prognosis in HCC patients was established. Multidimensional calculation demonstrated the predictive efficiency of the 18 immune-associated genes' signature, the accuracy of molecular phenotype in 18 immune-associated genes' signature, the correlation between 18 immune-associated genes' signature and clinical characteristics, and the systematic involvement of 18 immune-associated genes' signature in HCC functional pathways. Also, the prognostic value and accuracy of these 18 immune-associated genes' signature in the test data set were well validated and assessed.

Materials and Methods

Patients' information, data preparation, and processing

The HCC patients' information and corresponding clinical information were downloaded from the International Cancer Genome Consortium (ICGC) database using the ICGC data portal (Alaei-Mahabadi et al., 2016). The HCC immune-associated genes were collected from the ImmPort database, which is a public immunology database providing raw data and schemes from clinical trials, mechanism studies, and new methods of cellular molecular measurements (Zhou et al., 2018). Besides, four sets of data (HCCDB6: 225 HCC and 220 paracancer; HCCDB7: 80 HCC and 82 paracancer; HCCDB15: 356 HCC and 49 paracancer; HCCDB17: 115 HCC and 52 paracancer) were downloaded from HCCDB database, which were used to validate the survival prediction of 18 immune-associated genes' signature (Lian et al., 2018). Screening of 389 cases of RNA sequence samples from clinical information was done and these included 212 HCC cases and 177 paracancer cases. Of these cases, 212 HCC cases were selected and divided into train (106 cases) and test (106 cases) sets based on the following criteria: (1) HCC cases were randomly assigned to the train and the test set; and (2) patient age distribution, clinical staging, follow-up time, and mortality were similar between the two sets. One thousand three hundred thirty-four immune-associated genes of 212 HCC cases were obtained and were subjected to low-abundance gene filtration. To reduce the influence of low-abundance genes on the model, all genes with expression level less than 1 in HCC cases were eliminated, and finally, 1258 HCC immune-associated genes with high expression abundance were collected. This study was approved by the Clinical Research Ethics Committee of Shenzhen People's Hospital and HCC patients undergoing liver cancer surgery provided informed consent.

The selection of immune-associated genes for signature construction

Univariate Cox regression analysis model was used for each immune-associated gene and survival data analysis with p < 0.05 via R package survival coxph function. Univariate Cox regression analysis model was performed into clinical characters to reveal immune-associated gene-related risk signature. To identify immune-associated genes associated with survival prediction and to establish immune-associated gene risk score, multivariable Cox regression analysis was used to further refine these selected immune-associated genes through the addition of clinical risk signatures as covariables. A hierarchical cluster analysis of immune-associated genes with significant univariate and multivariable Cox regression factors was performed. Euclidean squared distance metric was used to identify consensus cluster based on immune-associated genes. The chi-square test was used to verify the significant differences in clinical risk signatures between consensus clusters (p < 0.01). However, the number of immune-associated genes that could be used to discriminate the consensus clusters was quite large and hence not useful for clinical testing. Therefore, there was a need to further narrow the range of immune-associated genes while maintaining high accuracy. Least absolute shrinkage and selection operator (LASSO) is a compression estimation, which constructs a penalty function to get a more refined model so that it compresses some coefficients and sets some to zero. It retains the advantage of subset contraction and it is a biased estimation for data with complex collinearity. It can realize the selection of variables while estimating parameters and better solving the multicollinearity problem in regression analysis (Li and Sillanpaa, 2012). Therefore, The R glmnet package was used for the LASSO Cox regression analysis to decrease the number of immune-associated genes.

Statistical analyses

The risk score for each patient was calculated based on the immune-associated genes as the following formula: Risk score = expression_{gene 1} × β _{gene 1} + expression_{gene 2} × β _{gene 2} + … + expression_{gene x} × β _{gene x}, where x was the number of immune-associated genes and β was the value of the coefficient for each immune-associated gene. The test set patients were divided into low risk and high risk based on the median risk score. For survival analysis, overall survival (OS) of the low-risk, high-risk, and different consensus clusters was subjected to Kaplan–Meier (KM) survival curve analysis (p < 0.01). The accuracy and sensitivity of survival prediction were based on the risk score (topography, lymph node, metastasis) tumor-node-metastasis (TNM) stage, and risk + TNM stage, which were verified by receiver operating characteristic (ROC) curve analysis. Multivariate Cox regression analysis on prediction relationship between risk score and clinical features to verify immune-associated genes was independently correlated with the clinical features. Furthermore, the risk score distribution in different clinical TNM stages and different clusters was analyzed by Student's t-test. Based on the results of multivariate Cox regression analysis, we used rms R package to construct a nomogram to provide a quantitative method for clinicians to predict the prognosis of HCC patients. In the meantime, calibration plots were used to assess the performance of the nomogram, and the 45° line was regarded as the best prediction. All the statistical analyses of this study were performed by R language, version 3.3.1.

Functional enrichment analysis of the immune-associated genes

To further understand the immune-associated gene function in HCC, GSVA package R of single-sample gene set enrichment analysis (ssGSEA) was used to analyze the Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment scores for each sample of the train set. Enrichment scores in each sample were calculated based on the normalized difference of empirical cumulative distribution functions of gene expression ranks. According to the risk scores, the functional pathway could be used to divide HCC cases into different clusters. In addition, Pearson's correlation analysis was used to investigate the association between risk scores and pathways.

Quantification of genes by quantitative real-time PCR

We collected surgical resection of cancer and paracancerous tissues from 37 patients with HCC who were newly diagnosed at Shenzhen People's Hospital between August and November 2019. Total RNA was separated using the TRIzol reagent (Invitrogen) according to the manufacturer's instructions. After extraction of RNA, the RNA was reverse transcribed to cDNA using PrimeScript RT polymerase (Takara, Dalian, China). SYBR Green Master Mix (Life Technologies) on a LightCycler 96 Detection System (Roche) was applied for RT-qPCR. GAPDH was used as an internal control for mRNA. Primers used in this study are listed in Supplementary Table S1. A comparison CT method was used to analyze the relative expression of RNA.

Results

Identification of immune-associated genes for HCC prognosis

Univariate Cox regression analysis for each immune-associated gene and for survival data was performed, and we obtained 220 immune-associated genes with significant prognostic difference with p < 0.05. Univariate Cox regression analysis was also conducted on the 220 immune-associated genes with clinical signatures, including tumor invasion, TNM stage, fibrosis, smoking, gender, alcohol, and age. The result showed that portal vein invasion, TNM stage, and vein invasion showed a significant difference in prognosis (p < 0.05) (Table 1). Furthermore, 220 immune-associated genes were subjected to multivariate Cox regression analysis in which the clinical signatures of portal vein invasion, TNM stage, and vein invasion were used as covariates. A total of 51 immune-associated genes displayed a significant prognostic potential and these were selected for further assessment. In addition, we determined whether the 51 immune-associated genes could divide cases into different cluster characteristics based on the clinical signature of portal vein invasion, TNM stage, and vein invasion. Thus, the 51 immune-associated genes were used for hierarchical clustering analysis, and HCC cases were divided into two clusters (Fig. 1A). The proportion of deaths, early stage, late stage, and invasion in cluster 1 and cluster 2 was 22% and 0%, 52% and 71%, 49% and 29%, and 41% and 12%, respectively (p < 0.001) (Fig. 1B). KM analysis showed that cluster 1 had significantly shorter OS than cluster 2 (p < 0.001) (Fig. 1C). These results demonstrated that the 51 immune-associated genes were related to HCC clinical signature and could be used to screen for prognostic markers of HCC.

FIG. 1.

The association between immune-associated genes' signature and clinical characteristics. (A) Hierarchical cluster analysis of expression profiles of immune genes with significant prognostic potential. (B) Differences between cluster 1 and 2 genes of the clinical signature determined by the chi-square test. (C) KM curve of cluster 1 and 2 based on hierarchical cluster analysis classification. KM, Kaplan–Meier.

Table 1.

Univariate Cox Regression Analysis of Clinical Signatures for Overall Survival

	β	HR	95% CI for HR	Wald test	p
Portal vein invasion	0.85	2.3	1.4–3.8	12	0.00056
TNM stage	0.89	2.4	1.3–4.5	7.9	0.0049
Vein invasion	1	2.7	1.3–6	6.4	0.011
Fibrosis	0.58	1.8	0.98–3.2	3.6	0.057
Smoking	−0.89	0.41	0.15–1.2	2.8	0.092
Gender	−0.87	0.42	0.14–1.2	2.6	0.11
Alcohol	−0.31	0.73	0.45–1.2	1.6	0.21
Age	0.0015	1	0.95–1.1	0	0.95
Bile duct invasion	−17	4.00E-08	0–Inf	0	1

CI, confidence interval; HR, hazards ratio; TNM, tumor-node-metastasis.

To maintain the accuracy and specificity of prognosis prediction of immune-associated genes, an LASSO Cox regression model was performed to narrow the scope of 51 immune-associated genes. First, the trajectory of each independent variable was analyzed, as shown in Figure 2A. It was evident that as the lambda increased, the number of independent coefficients was tending toward zero. The 10-fold cross-validation model construction and analyzing of the confidence intervals under each lambda are shown in Figure 2B. Based on Figure 2B, the model was optimal when γ = 29.39386 was used. Therefore, γ = 29.39386 was used as the final model and included 18 immune-associated genes (Table 2).

FIG. 2.

LASSO Cox regression analysis of immune-associated genes. (A) The trajectory of each independent variable. Lateral axis and vertical axis represent the lambda value and independent variable coefficient, respectively. (B) The confidence interval under each lambda. The red dotted line and solid vertical lines stood for the cross-validation curve and mean-squared error, respectively. The vertical dotted line represents the optimal values, which were plotted according to the minimum criteria and mean-squared error. LASSO, least absolute shrinkage and selection operator.

Table 2.

The Detail Information for 18 Immune-Associated Genes

Entrez ID	Gene ID	Gene name
3597	IL13RA1	Interleukin 13 receptor subunit alpha 1
4856	NOV	Nephroblastoma overexpressed
54567	DLL4	Delta_like canonical Notch ligand 4
55315	SLC29A3	Solute carrier family 29 member 3
353500	BMP8A	Bone morphogenetic protein 8a
2690	GHR	Growth hormone receptor
867	CBL	Cbl proto-oncogene
2264	fgfr4	Fibroblast growth factor receptor 4
7349	UCN	Urocortin
10397	NDRG1	N-myc downstream regulated 1
2649	NR6A1	Nuclear receptor subfamily 6 group A member 1
3480	IGF1R	Insulin-like growth factor 1 receptor
355	FAS	Fas cell surface death receptor
1524	CX3CR1	C-X3-C motif chemokine receptor 1
183	AGT	Angiotensinogen
23452	angptl2	Angiopoietin-like 2
7098	TLR3	Toll-like receptor 3
2215	FCGR3B	Fc fragment of IgG receptor IIIb

Application of the 18 immune-associated genes' risk signature

A risk score formula to build a risk signature for HCC based on the 18 immune-associated genes' expression was performed, which then depicted the risk score distribution of cluster 1 and cluster 2 (Fig. 3A). Based on this analysis, cluster 2 was distributed in the middle and back of risk score values, whereas cluster 1 was distributed in the front of risk score values, which showed a higher risk score. Moreover, compared with cluster 2, HCC cases of cluster 1 had shorter OS (Fig. 3B). These results suggested that the high-risk score HCC cases based on 18 immune-associated genes have a worse prognosis. Furthermore, ROC analysis was used to verify the prognostic classification of the risk score. All the area under the curve (AUC) of 1, 3, 5, and 10 years was greater than 0.7 (Fig. 3C), indicating the well prognostic classification based on the 18 immune-associated genes' risk score.

FIG. 3.

The risk signature of 18 immune-associated genes. (A) The association between 18 immune-associated genes and risk score distribution. (B) The KM survival analysis of cluster 1 and 2 signatures. (C) ROC curves of the risk score based on 18 immune-associated genes for 1, 3, 5, and 10 years. (D) The distribution of risk scores of tumor invasion for TNM-T stage. ROC, receiver operating characteristic.

To observe the relationship between the prediction accuracy of 18 immune-associated genes' risk score and their clinical signatures, multivariate Cox analysis was conducted. This further validated the 18 immune-associated genes used as independent robust signatures regardless of the adjusting of the signatures, including tumor invasion, TNM stage, fibrosis, smoking, gender, alcohol, and age (p = 0.00119) (Table 3). Multivariate Cox analysis showed that portal vein invasion (p = 0.03021), TNM stage (p = 0.04537), and vein invasion (p = 0.04195) were significantly associated with 18 immune-associated genes' risk score, an indication that those clinical signatures contributed to poor prognosis of HCC. Furthermore, the distribution of the risk score of tumor invasion (portal vein invasion, vein invasion, bile duct invasion) and the indifferent clinical TNM-T stages was analyzed. As shown in Figure 3D, there were significant differences in risk score among different stages, and the higher the stage, the higher the risk score, which indicated that a relationship existed between the risk score and clinical TNM stage.

Table 3.

Multivariate Cox Regression Analysis of Clinical Signatures for Overall Survival

	β	HR	95% CI for HR	p
Risk	1.18	3.25	1.59–6.64	0.00119
Alcohol	−0.79	0.45	0.21–1.00	0.04876
TNM stage	0.79	2.21	0.87–5.59	0.04537
Portal vein invasion	0.64	1.91	0.88–4.14	0.03021
Fibrosis	0.37	1.45	0.75–2.78	0.27012
Gender	−0.43	0.65	0.12–3.62	0.62256
Smoking	0.35	1.42	0.32–6.30	0.64428
Vein invasion	−0.25	0.78	0.25–2.37	0.04195
Age	0.01	1.01	0.95–1.07	0.72683
Bile duct invasion	−17.18	0	0.00–Inf	0.99789
Artery invasion	0	1	1.00–1.00	NA

In addition, ROC curve analysis was used to reveal the accuracy and sensitivity of prognostic prediction by making comparison among the risk scores, TNM stage, and risk score + TNM stage. AUC curve demonstrated that both risk score and TNM stage had a robust efficiency to predict the HCC prognosis. However, the efficiency of prognostic prediction of a combination of TNM stage and risk score was more significant than either a single risk score or TNM stage (Fig. 4A). The median risk scores of the three models (TNM stage, risk score, and risk score + TNM stage) were used as the threshold to divide the HCC cases into high- and low-risk groups, and then, OS analysis was performed. KM analysis represented that HCC cases of risk score + TNM stage in the high-risk group had shorter OS than the low-risk group, which showed the most significant difference in the prognosis of stage + risk score model with p < 0.0001 (Fig. 4B). Consistent results were found in HCC cases of risk score (p = 0.00071) (Fig. 4C) or TNM stage (p = 0.016) (Fig. 4D). It was evident that the efficiency of prognostic prediction of risk score + TNM-stage model was the most significant difference. All the results indicated that the risk score based on 18 immune-associated genes were associated with clinical signatures, and 18 immune-associated genes were an independent and efficient biomarker to predict HCC prognosis.

FIG. 4.

The accuracy and sensitivity of the prognostic prediction of TNM stage, risk score, and risk score + TNM stage. (A) ROC curves of TNM stage, risk score, and risk score + TNM stage. The KM survival analysis of (B) risk score + TNM stage, (C) risk score, and (D) TNM stage.

To provide a quantitative method that could be applied by clinicians for predicting the prognosis of HCC, a nomogram, which combined 18 immune-associated genes' risk score and clinical signature (gender, TNM stage, portal vein invasion, and vein invasion) (Fig. 5A), was used. Zero to 100 was treated as the value of the point assigned to each variable, and the sum of the values of the points of each variable was recorded as the total number of points. The 3-, 5-, and 10-year survival rates of HCC cases were determined by the vertical line of total point values and intersected at the 3-, 5-, and 10-year survival lateral axis. As shown in Figure 5, the risk score based on 18 immune-associated genes made more contribution to risk point compared with other clinical signature ranges from 0 to 100. Whereas gender contributed least to the risk point, it was consistent with the result of no significant association between risk score and gender by multivariate Cox analysis (p = 0.11; Table 3). As the total point values increased, the OS of 3, 5, and 10 years decreased. With the same total point values, the 10-year OS was less than 5 years, and the 5-year OS was less than 3 years. Calibration plots were used to assess the performance of the nomogram, which showed that the line was about the same as the 45° line (Fig. 5B), thereby indicating a good prediction of HCC prognosis by nomogram.

FIG. 5.

Construction of a nomogram for the 18 immune-associated genes. (A) The nomogram prediction of the 3-, 5-, and 10-year survival rate of HCC patients. (B) The calibration plot of the nomogram prediction. HCC, hepatocellular carcinoma.

Eighteen immune-associated genes related to pathway enrichment analysis

To explore the potential functional characteristics of 18 immune-associated genes, ssGSEA was carried out to elucidate the biological function based on enrichment scores for each sample. The most significant pathways are shown in Figure 6A, which revealed that most of the pathways were associated with tumorigenesis, such as JAK-STAT, MAPK, P53, cell receptor, and apoptosis, among others. Also, some pathways played a major role in promoting or suppressing the tumor through immune mechanisms, such as the toll-like receptor, natural killer cell-mediated cytotoxicity, T cell receptor, TGF beta, and chemokine. When 18 immune-associated genes' risk score was used for classification, the 20 pathways were applied to divide HCC cases into cluster 1 and cluster 2. The risk score distribution of the two clusters showed that the risk score in cluster 1 was significantly higher than cluster 2 (Student's t-test p < 0.05; Fig. 6B). Moreover, the relationship between the risk score and pathways was investigated by Pearson's correlation analysis and the pathways were associated with the risk score (Fig. 6C).

FIG. 6.

Analysis of the association between 18 immune-associated genes and pathway enrichment. (A) Hierarchical cluster analysis of functional characteristics of the 18 immune-associated genes. (B) The distribution of risk score based on 20 pathways of 18 immune-associated genes. (C) Relationship between pathways and risk score as determined by Pearson's correlation analysis.

Validation of the 18 immune-associated genes' signature in the test set

To verify the stability and reliability of the 18 immune-associated genes' signature, an 18 immune-associated genes' signature was performed into an independent test set. Similarly, the risk score of each HCC case was calculated, and the HCC cases of the test set divided into high-risk and low-risk groups based on the same risk score classification. The KM analysis showed that the OS of the high-risk group was shorter than the low-risk group, which was consistent with the train test (p = 0.00037; Fig. 7A). Furthermore, the ROC curve analysis revealed that the classification based on 18 immune-associated genes' risk score was efficient with an AUC value of 0.777 (Fig. 7B). The validation demonstrated the stability and reliability of 18 immune-associated genes' signature and indicated its potential for HCC prognostic prediction.

FIG. 7.

Validation of the 18 immune-associated genes' signature in the test set. (A) KM curve of clusters grouped by 18 immune-associated genes. (B) ROC curves of the test set based on 18 immune-associated genes.

Validation of the 18 immune-associated genes' signature by RT-qPCR experiments and external data set

Real-time Quantitative PCR (RT-qPCR) experiments of 18 immune-associated genes' signature on 37 HCC tissues and matched paracancerous tissues showed that most of the gene expressions were higher in HCC than paracancerous tissues (Supplementary Fig. S1). Although the expression of some genes both in cancer and paracancerous tissues was very low, the trend of the gene expression of cancer is higher than that of paracancerous tissues. Overall, the expression of genes in cancer was higher than paracancerous tissues, which was consistent with that high-risk score that HCC cases based on 18 immune-associated genes have a worse prognosis. To further validate the effectiveness of 18 immune-associated genes for HCC prognosis, the 18 immune-associated genes' signature was used to analyze in external big database. HCCDB database contains 15 public HCC expression data sets that cover around 4000 clinical HCC samples and serves as a one-stop online resource for exploring HCC gene expression. HCCDB also provides the results from computational analyses, including differential expression analysis, tissue-specific and tumor-specific expression analysis, survival analysis, and coexpression analysis (Lian et al., 2018). We can see that the high-risk group had shorter OS than the low-risk group in all four sets of data (HCCDB6, HCCDB7, HCCDB15, and HCCDB17) (Supplementary Fig. S2A). Furthermore, the ROC curve analysis validated the effectiveness of 18 immune-associated genes risk score with AUC value of 0.83, 0.65, 0.79, and 0.76 for HCCDB6, HCCDB7, HCCDB15, and HCCDB17, respectively (Supplementary Fig. S2B). All the results confirmed that the 18 immune-associated genes' risk score was a potential biomarker for the prognostic prediction of HCC.

Discussion

HCC has been associated with immunological reactions, especially chronic inflammation and fibrosis caused by different factors. These immunological reactions contribute to initiation and development of tumorigenesis (Nakamoto, 2017). Therefore, it is important to investigate the pathogenesis of HCC based on immunologic mechanisms. Recent studies have shown that immunotherapy has been successful in treating HCC patients (Tagliamonte et al., 2016). However, the mechanisms by which immunologic factors regulate HCC occurrence and whether immune genes are effective in diagnosis, treatment, and prognosis prediction are unclear. In this study, we constructed the 18 immune-associated genes' signature to predict the prognosis of HCC patients by multidimensional methods based on bioinformatic tools. The results showed that the 18 immune-associated genes' signature was robust, capable of evaluating the survival outcome of patients and was significantly associated with clinical features. The prognosis of most HCC patients was poor, mainly due to late diagnosis and improper prognosis assessment. This incentivized the need to develop better strategies for early diagnosis and prognostic prediction of HCC (Tsuchiya et al., 2015). Thus, this study presents effective prognostic biomarkers of 18 immune-associated genes and demonstrates their potential for application in prognostic prediction of HCC.

The signature in this research was based on immune-related genes and not only displayed better prognosis prediction of HCC patients, but also was significantly associated with clinical features and functional biology. The expression profile of immune genes involved in tumors has been identified, and their potential clinical correlations have been revealed in several cancers (Cheng et al., 2016; Li et al., 2017). Immune system genes are those that regulate the immune reactivity, and participate in the generation of immune responses in the wake of physical damage or other immune-related diseases (Wu et al., 2017). Previous research revealed that the characteristics of immune genes are highly conserved in many solid tumors, and they regulate various immune functions via diverse immune cell types. On this basis, tumor immune-specific genes can be used as effective biomarkers to evaluate the prognosis of cancer patients or to predict the response to immunoregulatory therapy (Chifman et al., 2016; Rhee et al., 2018). This research mined 18 immune-associated genes' signature and used them as molecular classifiers to evaluate the prognosis of HCC. The immune-associated 18 genes revealed that HCC patients with high-risk score had a worse prognosis compared with those with low-risk score. Although a large number of genes have been identified as cancer biomarkers, they are still not comprehensive enough to accurately achieve a single assessment of the prognosis of complex pathological diseases. In particular, the expression of a single gene could be modulated by a series of other factors. Owing to the identification of several molecular markers in tumor research, the molecular classification of tumors has become mainstream. Therefore, a bioinformatic statistical model was used to construct a gene signature comprising immune-associated genes. The model has been widely used in prognosis prediction and is superior to the prognosis based on a single-disease biomarker (Naoi and Noguchi, 2016; Montalbo et al., 2019). In addition, molecular-based tumor typing can complement the pathological classification of tumors in clinical prognosis. The molecular biomarkers provide more accurate clinical information for the diagnosis, treatment, and prognosis of patients (Li et al., 2016; Choi et al., 2017). Some studies suggest that the combination of molecular typing and pathological typing of tumors is superior to pathological typing alone in guiding clinical diagnosis and treatment. Here, we performed univariate Cox regression, multivariable Cox regression, and TNM stage + risk score ROC analysis, which confirmed the association between the immune genes and pathological classification (TNM stage, portal vein invasion, vein invasion, and bile duct invasion). From this analysis, we established the 18 immune-associated genes' signature. This signature was not only a robust prognostic biomarker but also was associated with clinical features, including tumor TNM-stage pathological classification. This research indicated that the 18 immune-associated genes' signature may complement HCC pathological classification, which is likely to promote personalized diagnosis and therapy of HCC patients.

Since 18 immune-associated genes' signature can distinguish patients into the high- and low-risk clusters, it was hypothesized that 18 immune-associated genes were involved in the biological function that affected the accuracy and specificity of prognostic prediction. Therefore, the 18 immune-associated gene-related pathway enrichment analysis based on enrichment scores was carried out by ssGSEA, and the result showed that the pathways were mainly enriched in two parts that included tumorigenesis and immune mechanisms. There is growing evidence that inflammation and the immune system are associated with cancer, and the tumor immune microenvironment has proven to be an important driver of biological behavior and a factor affecting chemotherapy (Pico de Coana et al., 2015; Andrejeva and Rathmell, 2017). Also, several studies have found that HCC does not adopt immune privilege, but coordinates powerful immune responses, involving innate and adaptive immune systems, which contribute to the occurrence and development of HCC (Zhou et al., 2017; Inada et al., 2019). The current research indicated the role of 18 immune-associated genes in the regulation of immune pathways and involved in tumorigenesis, which highlighted the value of 18 immune-associated genes in the study of HCC immune mechanism. The study also found that the most significant 20 pathways could be applied to divide HCC cases into cluster 1 with high risk and cluster 2 with low risk according to the gene risk score. Moreover, Pearson's correlation analysis found that the pathways were associated with the risk score. These results demonstrated that the 18 immune-associated genes' risk score was a robust biomarker for HCC prognosis and provided some in-depth immune mechanisms based on genetic molecular processes and progression of tumorigenesis in HCC.

The prognostic signature consisted of 18 immune-associated genes with strong prognostic potential. The study findings revealed that 9 of the 18 genes belonged to the immune factor receptor, which plays an important role in the inflammatory process of tumorigenesis and progression (Seebacher et al., 2019). Some of the nine genes related to immune factor receptor have been reported in previous research. Previous studies proved that five (IL13RA1, GHR, FGFR4, IGF1R, and FAS) of the nine genes related to immune factor receptor were associated with liver cirrhosis and progression into HCC, and thus, they were viewed as novel biomarkers for diagnosis, therapy, and prognosis. IL13RA1 plays a central role in a variety of inflammation, and it is found to be overexpressed in disease-related fibroblasts, neoplastic cells, and cancer (Kioi et al., 2006). GHR is the receptor for growth hormone, which is mediated by the STAT 5 pathway. The damage of liver GHR-STAT5 signal is related to abnormal glucose metabolism and fatty liver disease, and it contributes to hepatocytes sensitive to fibrosis and tumorigenicity transformation (Mueller et al., 2012). The well-known FGFR4 is involved in the regulation of tumor pathways, including cell proliferation, and cell differentiation, which has been revealed as a carcinogenic driver in HCC coupled with liver cirrhosis (Sheu et al., 2015; Joshi et al., 2017). Kim et al. (1996) revealed HBV may play a key role in the process of chronic hepatitis B by activating the expression of IGF1R, leading to liver fibrosis and progressing into HCC. FAS functions to modulate hepatic mitochondrial metabolism and fatty acid oxidation contributing to nonalcoholic fatty liver disease, which generally develops to liver fibrosis and HCC (Cullen and Martin, 2015; Item et al., 2017). This study explains the activation of potential carcinogenic signals by immune factor receptors (Pitt et al., 2016). Therefore, HCC patients with a high-risk score in 18 immune-associated genes' signature can reflect the increase of tumor inflammatory microenvironment, which promotes the progression of liver cirrhosis and leads to HCC with poor prognosis.

Conclusion

In conclusion, this study established the 18 immune-associated genes' signature with the ability to accurately and specifically act as independent predictors of prognosis in HCC patients. The 18 immune-related genes' signature provides new significance for immune mechanism and immunotherapy for HCC to guide the clinical diagnosis, treatment, and prognosis of HCC. In addition, the prognosis of HCC based on 18 immune-related genes' signature can be used as a complement pathological subtyping of HCC, which is not only accurate but also effective in identifying different immune phenotypes as biomarkers for better treatment of HCC.

Footnotes

Disclosure Statement

No competing financial interests exist.

Funding Information

This project was supported by the National Natural Science Foundation of China (No. 81671596) and the National Science Foundation for Young Scientists of China (No. 31700795).

Supplementary Material

Supplementary Figure S1

Supplementary Figure S2

Supplementary Table S1

References

Alaei-Mahabadi

, Bhadury

, Karlsson

J.W.

, Nilsson

J.A.

, and Larsson

(2016). Global analysis of somatic structural genomic alterations and their impact on gene expression in diverse human cancers. Proc Natl Acad Sci U S A, 113, 13768–13773.

Andrejeva

, and Rathmell

J.C.

(2017). Similarities and distinctions of cancer and immune metabolism in inflammation and tumors. Cell Metab, 26, 49–70.

Bruix

, and Sherman

(2011). Management of hepatocellular carcinoma: an update. Hepatology (Baltimore, MD), 53, 1020–1022.

Calderaro

, Couchy

, Imbeaud

, Amaddeo

, Letouze

, Blanc

J.F.

, et al. (2017). Histological subtypes of hepatocellular carcinoma are related to gene mutations and molecular tumour classification. J Hepatol, 67, 727–738.

Chen

D.S.

, and Mellman

(2013). Oncology meets immunology: the cancer-immunity cycle. Immunity, 39, 1–10.

Chen

, Zheng

, Baade

P.D.

, Zhang

, Zeng

, Bray

, et al. (2016). Cancer statistics in China, 2015. Cancer J Clin, 66, 115–132.

Cheng

, Ren

, Zhang

, Cai

, Liu

, Han

, et al. (2016). Bioinformatic profiling identifies an immune-related risk signature for glioblastoma. Neurology, 86, 2226–2234.

Chifman

, Pullikuth

, Chou

J.W.

, Bedognetti

, and Miller

L.D.

(2016). Conservation of immune gene signatures in solid tumors and prognostic implications. BMC Cancer, 16, 911.

Choi

, Ochoa

, McConkey

D.J.

, Aine

, Hoglund

, Kim

W.Y.

, et al. (2017). Genetic alterations in the molecular subtypes of bladder cancer: illustration in the Cancer Genome Atlas Dataset. Eur Urol, 72, 354–365.

10.

Cullen

S.P.

, and Martin

S.J.

(2015). Fas and TRAIL ‘death receptors' as initiators of inflammation: implications for cancer. Semin Cell Dev Biol, 39, 26–34.

11.

Dimitroulis

, Damaskos

, Valsami

, Davakis

, Garmpis

, Spartalis

, et al. (2017). From diagnosis to treatment of hepatocellular carcinoma: an epidemic problem for both developed and developing world. World J Gastroenterol, 23, 5282–5294.

12.

Hartke

, Johnson

, and Ghabril

(2017). The diagnosis and treatment of hepatocellular carcinoma. Semin Diagn Pathol, 34, 153–159.

13.

Inada

, Mizukoshi

, Seike

, Tamai

, Iida

, Kitahara

, et al. (2019). Characteristics of immune response to tumor-associated antigens and immune cell profile in patients with hepatocellular carcinoma. Hepatology (Baltimore, MD), 69, 653–665.

14.

Item

, Wueest

, Lemos

, Stein

, Lucchini

F.C.

, Denzler

, et al. (2017). Fas cell surface death receptor controls hepatic lipid metabolism by regulating mitochondrial function. Nat Commun, 8, 480.

15.

Joshi

J.J.

, Coffey

, Corcoran

, Tsai

, Huang

C.L.

, Ichikawa

, et al. (2017). H3B-6527 is a potent and selective inhibitor of FGFR4 in FGF19-driven hepatocellular carcinoma. Cancer Res, 77, 6999–7013.

16.

Kim

S.O.

, Park

J.G.

, and Lee

Y.I.

(1996). Increased expression of the insulin-like growth factor I (IGF-I) receptor gene in hepatocellular carcinoma cell lines: implications of IGF-I receptor gene activation by hepatitis B virus X gene product. Cancer Res, 56, 3831–3836.

17.

Kioi

, Seetharam

, and Puri

R.K.

(2006). N-linked glycosylation of IL-13R alpha2 is essential for optimal IL-13 inhibitory activity. FASEB J, 20, 2378–2380.

18.

, Cui

, Diehn

, and Li

(2017). Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol, 3, 1529–1537.

19.

, Wu

W.K.

, Xing

, Wong

S.H.

, Liu

, Fang

, et al. (2016). Distinct subtypes of gastric cancer defined by molecular characterization include novel mutational signatures with prognostic capability. Cancer Res, 76, 1724–1732.

20.

, and Sillanpaa

M.J.

(2012). Overview of LASSO-related penalized regression methods for quantitative trait mapping and genomic selection. Theor Appl Genet, 125, 419–435.

21.

Lian

, Wang

, Zhang

, Wang

, Luo

, Tang

, et al. (2018). HCCDB: A database of hepatocellular carcinoma expression atlas. Genomics Proteomics Bioinformatics, 16, 269–275.

22.

Makarova-Rusher

O.V.

, Medina-Echeverz

, Duffy

A.G.

, and Greten

T.F.

(2015). The yin and yang of evasion and immune activation in HCC. J Hepatol, 62, 1420–1429.

23.

Montalbo

, Lozano

J.J.

, Izquierdo

, Ingelmo-Torres

, BaNos

, Palou

, et al. (2019). Ability of a urine gene expression classifier to reduce the number of follow-up cystoscopies in bladder cancer patients. Transl Res, 208, 73–84.

24.

Mueller

K.M.

, Themanns

, Friedbichler

, Kornfeld

J.W.

, Esterbauer

, Tuckermann

J.P.

, et al. (2012). Hepatic growth hormone and glucocorticoid receptor signaling in body growth, steatosis and metabolic liver cancer development. Mol Cell Endocrinol, 361, 1–11.

25.

Nakamoto

(2017). Promising new strategies for hepatocellular carcinoma. Hepatol Res, 47, 251–265.

26.

Naoi

, and Noguchi

(2016). Multi-gene classifiers for prediction of recurrence in breast cancer patients. Breast Cancer (Tokyo, Japan), 23, 12–18.

27.

Pico de Coana

, Choudhury

, and Kiessling

(2015). Checkpoint blockade for cancer therapy: revitalizing a suppressed immune system. Trends Mol Med, 21, 482–491.

28.

Pitt

J.M.

, Vetizou

, Daillere

, Roberti

M.P.

, Yamazaki

, Routy

, et al. (2016). Resistance mechanisms to immune-checkpoint blockade in cancer: tumor-intrinsic and -extrinsic factors. Immunity, 44, 1255–1269.

29.

Rhee

J.K.

, Jung

Y.C.

, Kim

K.R.

, Yoo

, Kim

, Lee

Y.J.

, et al. (2018). Impact of tumor purity on immune gene expression and clustering analyses across multiple cancer types. Cancer Immunol Res, 6, 87–97.

30.

Seebacher

N.A.

, Stacy

A.E.

, Porter

G.M.

, and Merlot

A.M.

(2019). Clinical development of targeted and immune based anti-cancer therapies. J Exp Clin Cancer Res, 38, 156.

31.

Sheu

M.J.

, Hsieh

M.J.

, Chiang

W.L.

, Yang

S.F.

, Lee

H.L.

, Lee

L.M.

, et al. (2015). Fibroblast growth factor receptor 4 polymorphism is associated with liver cirrhosis in hepatocarcinoma. PLoS One, 10, e0122961.

32.

Tagliamonte

, Petrizzo

, Tornesello

M.L.

, Ciliberto

, Buonaguro

F.M.

, and Buonaguro

(2016). Combinatorial immunotherapy strategies for hepatocellular carcinoma. Curr Opin Immunol, 39, 103–113.

33.

Tan

H.W.

, Leung

C.O.

, Chan

K.K.

, Ho

D.W.

, Leung

M.S.

, Wong

C.M.

, et al. (2019). Deregulated GATA6 modulates stem cell-like properties and metabolic phenotype in hepatocellular carcinoma. Int J Cancer, 145, 1860.

34.

Tsuchiya

, Sawada

, Endo

, Saito

, Uemura

, and Nakatsura

(2015). Biomarkers for the early diagnosis of hepatocellular carcinoma. World J Gastroenterol, 21, 10573–10583.

35.

Wang

, Ruan

, Yu

, Tian

, Liang

, Jing

, et al. (2019). A four-methylated mRNA signature-based risk score system predicts survival in patients with hepatocellular carcinoma. Aging, 11, 160–173.

36.

C.E.

, Yu

C.W.

, Chang

K.W.

, Chou

W.H.

, Lu

C.Y.

, Ghelfi

, et al. (2017). Comparative global immune-related gene profiling of somatic cells, human pluripotent stem cells and their derivatives: implication for human lymphocyte proliferation. Exp Mol Med, 49, e376.

37.

, Xue

, Wang

, Shen

, Jiang

, Zhang

, et al. (2017). miR-345 inhibits tumor metastasis and EMT by targeting IRF1-mediated mTOR/STAT3/AKT pathway in hepatocellular carcinoma. Int J Oncol, 50, 975–983.

38.

Zhou

, Bi

, Yuan

, Yang

, and Sun

(2018). Gain of UBE2D1 facilitates hepatocellular carcinoma progression and is associated with DNA damage caused by continuous IL-6. J Exp Clin Cancer Res, 37, 290.

39.

Zhou

, Sprengers

, Boor

P.P.C.

, Doukas

, Schutz

, Mancham

, et al. (2017). Antibodies against immune checkpoint molecules restore functions of tumor-infiltrating T cells in hepatocellular carcinomas. Gastroenterology, 153, 1107–1119.e10.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.05 MB

0.22 MB

0.05 MB