# Immunogenomic Landscape Analysis of Prognostic Immune-Related Genes in Hepatocellular Carcinoma

Authors: Peng Qin; Mengyu Zhang; Xue Liu; Ziming Dong
Journal: Journal of Healthcare Engineering
 (2021)

---

## Abstract

Hepatocellular carcinoma (HCC) is one of the leading causes of cancer-related death. HBV infection is an important risk factor for the tumorigenesis of HCC, given that the inflammatory environment is closely related to morbidity and prognosis. Consequently, it is of urgent importance to explore the immunogenomic landscape to supplement the prognosis of HCC. The expression profiles of immune‐related genes (IRGs) were integrated with 377 HCC patients to generate differentially expressed IRGs based on the Cancer Genome Atlas (TCGA) dataset. These IRGs were evaluated and assessed in terms of their diagnostic and prognostic values. A total of 32 differentially expressed immune‐related genes resulted as significantly correlated with the overall survival of HCC patients. The Gene Ontology functional enrichment analysis revealed that these genes were actively involved in cytokine‐cytokine receptor interaction. A prognostic signature based on IRGs (HSPA4, PSME3, PSMD14, FABP6, ISG20L2, TRAF3, NDRG1, NRAS, CSPG5, and IL17D) stratified patients into high-risk versus low-risk groups in terms of overall survival and remained as an independent prognostic factor in multivariate analyses after adjusting for clinical and pathologic factors. Several IRGs (HSPA4, PSME3, PSMD14, FABP6, ISG20L2, TRAF3, NDRG1, NRAS, CSPG5, and IL17D) of clinical significance were screened in the present study, revealing that the proposed clinical-immune signature is a promising risk score for predicting the prognosis of HCC.


---

## Body

## 1. Introduction

Hepatocellular carcinoma (HCC) ranks seventh among malignant tumors in terms of incidence. According to the latest Global Cancer Statistics, 841,080 new incidents of HCC and 781,631 deaths occurred during the year 2018 [1]. China is a country with a high incidence of HCC, which accounts for more than half of the world’s deaths [2]. With the development of modern medical science and technology, significant progress has been made in the treatment of HCC. As the clinical symptoms of early HCC are not typical, 70% to 80% of the patients have advanced disease at the time of diagnosis [3]. Existing treatment strategies are insufficient for patients with advanced HCC. Therefore, identifying novel and sensitive biomarkers is of critical importance for the early diagnosis of HCC.In previous studies, the therapeutic response of HCC patients was stratified based on the identification of molecular biomarkers, such as genes, microinterfering RNA (miRNA), circular RNA (circRNA), and long noncoding RNA (lncRNA). Chen et al. reported a four-gene (KPNA2, CDC20, SPP1, and TOP2A) based signature, which could be a candidate prognostic factor for patients with HCC [4]. The deregulation of miRNA-122 has been related to an increased risk of developing HCC [5]. Also, the upregulation of miRNA-372 has been associated with tumor progression and prognosis in HCC [6]. Several circular RNAs such as circRNA_0001955 [7] and circRNA_101505 [8] have been identified as potential biomarkers for HCC diagnosis and prognosis. Moreover, a five-long noncoding RNAs signature has been reported to improve survival prediction and be used as a prognostic biomarker for HCC patients [9].The liver is an essential organ for the proper functioning of the immune system, which is rich in various immune cells, and especially the cytotoxic T lymphocyte (CTL) that can recognize tumor antigens and eliminate the tumor cells from the tumor microenvironment. Over the last decade, cancer immunotherapy has proven to be a promising treatment protocol for various types of cancer [10, 11]. Certain studies have shown that HCC cells, which are in a highly immunosuppressive microenvironment, can induce host immunosuppression and avoid autoimmune response by downregulating major histocompatibility complex-1 (MHC-1), secreting immunosuppressive cytokines, and mediating negative costimulatory signals [12, 13]. Cancer immunotherapy can delay the progress of tumors by enhancing the immune response of the body, stimulating specific immunity of tumors, and breaking immune tolerance [14, 15]. Over recent decades, immunotherapy has been applied in the treatment of various types of tumors [11, 16, 17] and immune checkpoint inhibitors have become potential effective treatment in patients with advanced HCC [10]. In September 2017, nivolumab was approved by the FDA for liver cancer as a second-line treatment after failure of sorafenib based on the data of the multicohort phase 1/2 trial CheckMate 040 [18]. New immunotherapy technologies, such as chimeric antigen receptor T cells (CAR-T) [19], T cell receptor genetically engineered T cells (TCR-T) [20, 21], new antigen vaccines, and oncolytic viruses, have gradually found their application in clinic. These clinical results fully demonstrate the importance of immunology in liver cancer, so it is crucial to understand these molecular mechanisms, especially immune gene effects. The emergence of public, large-scale gene expression datasets has enabled researchers to identify responsible biomarkers for tumor monitoring and surveillance with much accuracy [22, 23]. The prognostic value of immune-related genes (IRGs) was explored to develop an individualized immune signature, which could improve prognostic estimation in patients with nonsquamous non-small cell lung cancer [24].The purpose of this research was to investigate whether IRGs have potential prognostic value for HCC and whether they can be used as biomarkers for immunotherapy. Initially, we combined the transcriptome RNA-sequencing data downloaded from TCGA to analyze the differentially expressed genes and differentially expressed immune-related genes in HCC. Then, we integrated IRGs expression profiles with clinical information, applying computational methods for the assessment of overall survival (OS) in HCC patients. Making the best of the complementary value of IRGs expression profiles and clinical characteristics, we investigated the potential clinical utility of IRGs on prognostic stratification and their implicational potential as biomarkers for targeted HCC therapy. Eventually, we build an individualized prognostic signature, which may support HCC prognosis. The study has the following contributions in this regard:(i)
The immune genes related differentially expressed genes (DEGs) were discovered, and an immune-related gene-based prognostic index (IRGPI) consisting of 10 genes (HSPA4, PSME3, PSMD14, FABP6, ISG20L2, TRAF3, NDRG1, NRAS, CSPG5, IL17D) was built using a Cox regression model.(ii)
The IRGPI might be employed as an independent prognostic factor, according to both univariate and multivariate Cox analyses. The IRGPI’s reliability was demonstrated through ROC analysis.(iii)
We investigated the link between this IRG risk profile and important clinical characteristics after establishing it (age, gender, survival state, and grade).


## 2. Materials and Methods

### 2.1. Transcriptome Expression Data and Clinical Information Acquisition

The transcriptome expression profiles and corresponding clinical information of hepatocellular carcinoma were downloaded from the Genomic Data Commons Data Portal of TCGA (https://cancergenome.nih.gov/), which contained data from 374 hepatocellular carcinoma and 50 noncancerous liver tissues. The IRGs list was derived from the Immunology Database and Analysis Portal (ImmPort) database [25].


### 2.2. Differential Gene Analysis

Differentially expressed genes (DEGs) between HCC and nontumor samples were screened by the R software edgeR package (http://bioconductor.org/packages/edgeR/) to select DEGs related to hepatocarcinogenesis [26]. The raw data were normalized by the Trimmed mean of M values (TMM) implemented in the edgeR Bioconductor package. Gene expression comparison was carried out by calculating the level of fold change (FC) in HCC versus noncancerous liver tissue with a false discovery rate (FDR) <0.05 and a log2 |fold change| >1 as the cutoff values. Differentially expressed IRGs were then extracted from all DEGs. The functional enrichment of Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of these differentially expressed IRGs was performed on the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) [27, 28].


### 2.3. Survival-Associated IRGs Analysis

The follow-up data of HCC patients were derived from TCGA’s Pan-Cancer Atlas. Differentially expressed IRGs, which were significantly correlated to overall survival (OS) in HCC patients, were selected via univariate COX analysis that was conducted using the R software survival package. These survival-associated IRGs were also used for functional enrichment analysis. Copy number alterations data of these IRGs were obtained from Cbioportal (http://www.cbioportal.org/) [26]. To clear up the potential molecular mechanisms of these survival-associated IRGs, we focused on the transcription factors (TFs), which are essential molecules that directly control the degree of gene expression. The expression profiles of 318 transcription factors (TFs) were downloaded from the Cistrome Cancer database, which is a valuable resource for experimental and computational cancer biology research [29]. Besides, a functional network between the TFs and these survival-associated IRGs was constructed via Cytoscape (version 2.8, http://cytoscape.org).


### 2.4. Construction and Validation of the Immune-Related Gene-Based Prognostic Index (IRGPI)

These selected survival-related IRGs were submitted for multivariate analyses, with integrated IRGs remaining as independent prognostic indicators for the development of the IRGPI. Prognostic IRGs with a false discovery rate of less than 0.05 were candidates to calculate the risk score value. Based on the results of the median risk score value, the IRGPI significantly stratified patients into high- and low-risk groups. The optimal IRGPI cutoff was determined by a time-dependent receiver operating characteristic (ROC) curve [30] at 5 years.


### 2.5. Statistical Analysis

Gene functional enrichment analyses were performed using R (version 3.6.1;https://www.r-project.org/) software cluster Profiler package [31]. AUC of the survival ROC curve was calculated by the survival ROC R software package to verify the reliability of the prognostic signature [30]. The differences in clinical parameters were tested by independent t-tests. Statistical significance was defined as P<0.05.


## 2.1. Transcriptome Expression Data and Clinical Information Acquisition

The transcriptome expression profiles and corresponding clinical information of hepatocellular carcinoma were downloaded from the Genomic Data Commons Data Portal of TCGA (https://cancergenome.nih.gov/), which contained data from 374 hepatocellular carcinoma and 50 noncancerous liver tissues. The IRGs list was derived from the Immunology Database and Analysis Portal (ImmPort) database [25].


## 2.2. Differential Gene Analysis

Differentially expressed genes (DEGs) between HCC and nontumor samples were screened by the R software edgeR package (http://bioconductor.org/packages/edgeR/) to select DEGs related to hepatocarcinogenesis [26]. The raw data were normalized by the Trimmed mean of M values (TMM) implemented in the edgeR Bioconductor package. Gene expression comparison was carried out by calculating the level of fold change (FC) in HCC versus noncancerous liver tissue with a false discovery rate (FDR) <0.05 and a log2 |fold change| >1 as the cutoff values. Differentially expressed IRGs were then extracted from all DEGs. The functional enrichment of Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of these differentially expressed IRGs was performed on the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) [27, 28].


## 2.3. Survival-Associated IRGs Analysis

The follow-up data of HCC patients were derived from TCGA’s Pan-Cancer Atlas. Differentially expressed IRGs, which were significantly correlated to overall survival (OS) in HCC patients, were selected via univariate COX analysis that was conducted using the R software survival package. These survival-associated IRGs were also used for functional enrichment analysis. Copy number alterations data of these IRGs were obtained from Cbioportal (http://www.cbioportal.org/) [26]. To clear up the potential molecular mechanisms of these survival-associated IRGs, we focused on the transcription factors (TFs), which are essential molecules that directly control the degree of gene expression. The expression profiles of 318 transcription factors (TFs) were downloaded from the Cistrome Cancer database, which is a valuable resource for experimental and computational cancer biology research [29]. Besides, a functional network between the TFs and these survival-associated IRGs was constructed via Cytoscape (version 2.8, http://cytoscape.org).


## 2.4. Construction and Validation of the Immune-Related Gene-Based Prognostic Index (IRGPI)

These selected survival-related IRGs were submitted for multivariate analyses, with integrated IRGs remaining as independent prognostic indicators for the development of the IRGPI. Prognostic IRGs with a false discovery rate of less than 0.05 were candidates to calculate the risk score value. Based on the results of the median risk score value, the IRGPI significantly stratified patients into high- and low-risk groups. The optimal IRGPI cutoff was determined by a time-dependent receiver operating characteristic (ROC) curve [30] at 5 years.


## 2.5. Statistical Analysis

Gene functional enrichment analyses were performed using R (version 3.6.1;https://www.r-project.org/) software cluster Profiler package [31]. AUC of the survival ROC curve was calculated by the survival ROC R software package to verify the reliability of the prognostic signature [30]. The differences in clinical parameters were tested by independent t-tests. Statistical significance was defined as P<0.05.


## 3. Results

### 3.1. Identification of Differentially Expressed IRGs

Transcriptional expression profiles and phenotype data of 377 HCC patients from the TCGA cohort were downloaded and integrated. Among them, there were 255 males and 122 females. The edgeR algorithm identified a total of 7,667 differentially expressed genes, 7,273 upregulated and 394 downregulated genes with the threshold of |log2FC| >1 and FDR <0.05 (Figures1(a) and 1(b)). From this set of genes, we extracted 329 differentially expressed IRGs, including 267 upregulated and 62 downregulated (Figures 1(c) and 1(d)). As expected, gene functional enrichment analysis revealed that inflammatory pathways were most frequently implicated. “Immune response,” “extracellular space,” and “growth factor activity” were the most frequent biological terms among biological processes, cellular components, and molecular functions, respectively (Figure 2(a)). For the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, cytokine-cytokine receptor interactions were most often enriched by differentially expressed IRGs (Figure 2(b)).Figure 1
Differentially expressed immune-related genes. Heat map (a) and a volcano plot (b) and a volcano plot (c) demonstrating differentially expressed genes between HCC and nontumor tissues. Heat map (d) demonstrating differentially expressed immune-related genes (IRGs) between HCC and nontumor tissues. Green dots represent downregulated genes, and red dots represent upregulated genes.
(a)(b)(c)(d)Figure 2
Gene functional enrichment of differentially expressed immune-related genes. (a) Gene ontology analysis; blue, red, and green bars represent a biological process, cellular component, and molecular function, respectively. (b) The top 10 most significant Kyoto Encyclopedia of Genes and Genomes Pathways.
(a)(b)


### 3.2. Identification and Characteristics of Survival-Associated IRGs

As predicting the prognosis is essential for clinical guidance, we focused on uncovering molecular biomarkers that could serve as viable prognostic indicators. In total, 32 IRGs that were significantly correlated with overall survival (OS) HCC patients (P<0.001; Table 1) were identified. Protein-protein interaction (PPI) network analysis demonstrated that HSP90AA1, PSMD10, and PSMD2 were the three hub genes among these datasets (Figure 3). A forest plot of expression profiles revealed that all of the 32 survival-associated IRGs were upregulated in HCC samples (Figure 4(a)). Given the important clinical significance of these IRGs, genetic alterations of these genes were examined, revealing that mRNA upregulation and amplification were the two most commonly occurring types of mutations (Figure 4(b)).Table 1
General characteristics of HCC survival-associated immune-related genes.
Gene symbolHRHR.95LHR.95HP valueHSPA41.0453663881.026195211.064895712.63E − 06HSP90AA11.0043132711.002007121.006624730.00024306PSMD21.020371921.009044311.03182670.00039904PSMD101.0516631351.023051881.081074550.0003444PSME31.0465705741.019774341.074070920.00058249PSMD141.1084183561.067902511.150471366.03E − 08IFI302.1248221911.392668493.241883740.00047121S100A101.0032152151.001536861.004896380.00017156S100A111.0015038161.000697371.002310910.00025609FABP61.1079410781.051312921.167619460.00012847ISG20L21.1506350341.091745411.212701211.65E − 07PPIA1.0127056741.00707131.018371579.19E − 06CACYBP1.0631697741.041890131.084884042.88E − 09TRAF31.3157934581.142871111.514879850.00013473DCK1.1299999021.063237311.200954648.37E − 05EED1.3709383461.138812521.650378720.0008582NDRG11.0074190931.004165661.010683077.51E − 06HDAC11.0430482361.026931441.059417971.13E − 07BIRC51.0286628181.014682991.042835255.17E − 05NRAS1.0760016331.047014131.105791681.46E − 07PLXNA11.1382905141.057028861.225799360.00060887PLXNA31.2391483571.111449391.381519180.00011147CSPG51.606267261.248640372.066323150.00022602GMFB1.1300501611.061485161.203044010.00012899IL17D1.1023301761.049876971.157404028.98E − 05KITLG1.2245507621.100158931.363007230.00021014STC21.0354484231.015675421.055606360.00039851BRD81.1574811071.066764791.255911820.00044469NR6A11.3001644481.12391051.50405890.00041304TNFRSF11A1.4165221671.151838611.742027950.00096886SHC11.0140173671.007070251.021012417.23E − 05CDK41.030718611.014076491.047633850.00026943Figure 3
Protein-protein interaction (PPI) network analysis.Figure 4
Prognostic values and mutation landscape of survival-associated immune‐related genes. Forest plot (a) of hazard ratios showing the prognostic values of genes. Mutation landscape (b) showed that ISG20L2 is the gene with the highest mutation frequency and there are 30 genes with a mutation rate ≥5%.
(a)(b)


### 3.3. TF Regulatory Network

To further understand the potential molecular mechanisms of these clinically related IRGs, we analyzed the regulatory mechanisms of these genes. The expression profiles of 318 transcription factors (TFs) were examined, and 117 were found to be differentially expressed between HCC and nontumor liver samples with the threshold of |log2FC| >2 and FDR <0.05 (Figures5(a) and 5(b)). The correlation analysis was constructed between these 117 TFs and the 32 survival-associated IRGs with a correlation score of more than 0.6 set as the cutoff values. The regulatory schematic acutely illustrated the regulatory relationships among these IRGs (Figure 5(c)).Figure 5
Transcription factor‐mediated regulatory network. Heat map (a) and a volcano plot (b) demonstrating differentially expressed TFs between hepatocellular carcinoma (HCC) and nontumor tissues, green dots represent downregulated TFs, and red dots represent upregulated TFs. (c) Regulatory network constructed based on differentially expressed TFs and survival-associated IRGs. Blue triangles represent TFs whereas orange circles represent high-risk IRGs.
(a)(b)(c)


### 3.4. Evaluation of Clinical Outcomes

In this study, we developed a prognostic signature based on the results of multivariate Cox regression analysis to divide the HCC patients into two groups with discrete clinical outcomes with regard to OS (Figure6). This index was based on the following formula: risk score = [expression level of HSPA4∗(0.029292)] + [expression level of PSME3∗(−0.04914)] + [expression level of PSMD14∗(0.054544)]+ [expression level of FABP6∗(0.054469)]] + [expression level of ISG20L2∗(0.096535)] + [expression level of TRAF3∗(0.153954) + [expression level of NDRG1∗(0.006116)] + [expression level of NRAS∗(0.040485)] + [expression level of CSPG5∗(0.238526)] + [expression level of IL17D∗(0.060712)]. The immune-based prognostic index (IRGPI) significantly stratified patients into low-risk (IRGPI < median value) and high-risk (IRGPI > median value) groups in terms of overall survival (Figure 7(a)). The area under the curve of the receiver operating characteristic (ROC) curve was 0.826, which suggested the moderate potential for the prognostic signature based on IRGs in survival monitoring (Figure 7(b)).Figure 6
Development of the prognostic index based on immune‐related genes. (a) The rank of prognostic index and distribution of groups. (b) Survival status of patients in different groups. (c) Heat map of expression profiles of these included genes.
(a)(b)(c)Figure 7
The prognostic value of the prognostic index based on immune‐related genes. (a) Patients in the high‐risk group suffered shorter overall survival. (b) Survival‐dependent receiver operating characteristic (ROC) curve validation of the prognostic value of the prognostic index. (c) Univariate and (d) multiple regression analysis of HCC.
(a)(b)(c)(d)The clinical data and risk scores were analyzed by univariate and multivariate regression analysis. TheP value of risk score was less than 0.001 in both univariate (Figure 7(c)) and multivariate regression (Figure 7(d)) analyses. These results indicated that the IRGPI obtained by our model could be used as an independent predictor after adjusting for other parameters, including age, gender, grade, tumor stage, tumor size, distant metastasis status, and the amounts of nodules (Table 2).Table 2
Univariate and multiple regression analysis of HCC.
VariablesUnivariate analysisMultivariate analysisHazard ratio (95% CI)P valueHazard ratio (95% CI)P valueAge0.995 (0.976–1.014)0.602599261.010 (0.990–1.031)0.323971151Gender (male/female)0.893 (0.532–1.500)0.669262011.555 (0.850–2.844)0.151666718Grade1.083 (0.772–1.520)0.644037011.057 (0.728–1.536)0.769568712Tumor stage2.084 (1.590–2.733)1.07E − 071.192 (0.404–3.521)0.750289146Tumor1.980 (1.541–2.543)8.94E − 081.515 (0.569–4.0340.405391572Distant metastasis4.769 (1.485–15.311)0.008669522.593 (0.637–10.553)0.183412547Node2.439 (0.593–10.035)0.216733452.311 (0.356–15.002)0.379943568Risk score1.264 (1.198–1.334)1.38E − 171.251 (1.172–1.336)2.10E − 11


### 3.5. Clinical Utility of IRGPI

To further assess the clinical value of the immune-related gene-based prognostic index (IRGPI), the relationship between this hub survival-associated IRGs and clinical characteristics including age, gender, survival state, grade, pathological stage, T stage, M stage, and N stage were analyzed (Table3). IRGPI showed a significant difference in survival state (Figure 8(a)), grade (Figure 8(b)), T stage (Figure 8(c)), and M stage (Figure 8(d)). However, no difference was observed between age, gender, and N stage.Table 3
Relationship between these survival-associated IRGs and clinical characteristics.
GenesAge (≥65/<65)Gender (male/female)State (dead/alive)Grade (IV-III/I-II)TPtPtPtPHSPA40.0240.981−1.2310.22−2.8330.006−3.4796.31E − 04PSME3−0.8870.377−0.3960.692−2.7110.008−2.3170.022PSMD14−1.0250.3070.3950.694−3.3180.001−2.0670.04FABP6−0.1430.886−1.2190.224−1.4130.162−1.7920.076ISG20L2−1.9470.0530.490.625−3.2470.002−2.8370.005TRAF30.2930.771.0480.297−2.1650.033−1.420.157NDRG10.4190.677−1.0920.276−2.5710.012−0.3560.723NRAS−1.7610.0810.0010.999−3.2140.002−2.2910.023CSPG5−0.3280.743−0.9020.368−2.3850.019−2.0120.046IL17D−1.3210.1881.4950.139−1.8960.062−2.0190.046Risk score−0.3160.753−0.2940.769−3.8612.60E − 04−2.1480.033Figure 8
The relationships between the immune-based prognostic index and (a) survival state, (b) grade, (c) T stage, (d) M stage.
(a)(b)(c)(d)


## 3.1. Identification of Differentially Expressed IRGs

Transcriptional expression profiles and phenotype data of 377 HCC patients from the TCGA cohort were downloaded and integrated. Among them, there were 255 males and 122 females. The edgeR algorithm identified a total of 7,667 differentially expressed genes, 7,273 upregulated and 394 downregulated genes with the threshold of |log2FC| >1 and FDR <0.05 (Figures1(a) and 1(b)). From this set of genes, we extracted 329 differentially expressed IRGs, including 267 upregulated and 62 downregulated (Figures 1(c) and 1(d)). As expected, gene functional enrichment analysis revealed that inflammatory pathways were most frequently implicated. “Immune response,” “extracellular space,” and “growth factor activity” were the most frequent biological terms among biological processes, cellular components, and molecular functions, respectively (Figure 2(a)). For the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, cytokine-cytokine receptor interactions were most often enriched by differentially expressed IRGs (Figure 2(b)).Figure 1
Differentially expressed immune-related genes. Heat map (a) and a volcano plot (b) and a volcano plot (c) demonstrating differentially expressed genes between HCC and nontumor tissues. Heat map (d) demonstrating differentially expressed immune-related genes (IRGs) between HCC and nontumor tissues. Green dots represent downregulated genes, and red dots represent upregulated genes.
(a)(b)(c)(d)Figure 2
Gene functional enrichment of differentially expressed immune-related genes. (a) Gene ontology analysis; blue, red, and green bars represent a biological process, cellular component, and molecular function, respectively. (b) The top 10 most significant Kyoto Encyclopedia of Genes and Genomes Pathways.
(a)(b)


## 3.2. Identification and Characteristics of Survival-Associated IRGs

As predicting the prognosis is essential for clinical guidance, we focused on uncovering molecular biomarkers that could serve as viable prognostic indicators. In total, 32 IRGs that were significantly correlated with overall survival (OS) HCC patients (P<0.001; Table 1) were identified. Protein-protein interaction (PPI) network analysis demonstrated that HSP90AA1, PSMD10, and PSMD2 were the three hub genes among these datasets (Figure 3). A forest plot of expression profiles revealed that all of the 32 survival-associated IRGs were upregulated in HCC samples (Figure 4(a)). Given the important clinical significance of these IRGs, genetic alterations of these genes were examined, revealing that mRNA upregulation and amplification were the two most commonly occurring types of mutations (Figure 4(b)).Table 1
General characteristics of HCC survival-associated immune-related genes.
Gene symbolHRHR.95LHR.95HP valueHSPA41.0453663881.026195211.064895712.63E − 06HSP90AA11.0043132711.002007121.006624730.00024306PSMD21.020371921.009044311.03182670.00039904PSMD101.0516631351.023051881.081074550.0003444PSME31.0465705741.019774341.074070920.00058249PSMD141.1084183561.067902511.150471366.03E − 08IFI302.1248221911.392668493.241883740.00047121S100A101.0032152151.001536861.004896380.00017156S100A111.0015038161.000697371.002310910.00025609FABP61.1079410781.051312921.167619460.00012847ISG20L21.1506350341.091745411.212701211.65E − 07PPIA1.0127056741.00707131.018371579.19E − 06CACYBP1.0631697741.041890131.084884042.88E − 09TRAF31.3157934581.142871111.514879850.00013473DCK1.1299999021.063237311.200954648.37E − 05EED1.3709383461.138812521.650378720.0008582NDRG11.0074190931.004165661.010683077.51E − 06HDAC11.0430482361.026931441.059417971.13E − 07BIRC51.0286628181.014682991.042835255.17E − 05NRAS1.0760016331.047014131.105791681.46E − 07PLXNA11.1382905141.057028861.225799360.00060887PLXNA31.2391483571.111449391.381519180.00011147CSPG51.606267261.248640372.066323150.00022602GMFB1.1300501611.061485161.203044010.00012899IL17D1.1023301761.049876971.157404028.98E − 05KITLG1.2245507621.100158931.363007230.00021014STC21.0354484231.015675421.055606360.00039851BRD81.1574811071.066764791.255911820.00044469NR6A11.3001644481.12391051.50405890.00041304TNFRSF11A1.4165221671.151838611.742027950.00096886SHC11.0140173671.007070251.021012417.23E − 05CDK41.030718611.014076491.047633850.00026943Figure 3
Protein-protein interaction (PPI) network analysis.Figure 4
Prognostic values and mutation landscape of survival-associated immune‐related genes. Forest plot (a) of hazard ratios showing the prognostic values of genes. Mutation landscape (b) showed that ISG20L2 is the gene with the highest mutation frequency and there are 30 genes with a mutation rate ≥5%.
(a)(b)


## 3.3. TF Regulatory Network

To further understand the potential molecular mechanisms of these clinically related IRGs, we analyzed the regulatory mechanisms of these genes. The expression profiles of 318 transcription factors (TFs) were examined, and 117 were found to be differentially expressed between HCC and nontumor liver samples with the threshold of |log2FC| >2 and FDR <0.05 (Figures5(a) and 5(b)). The correlation analysis was constructed between these 117 TFs and the 32 survival-associated IRGs with a correlation score of more than 0.6 set as the cutoff values. The regulatory schematic acutely illustrated the regulatory relationships among these IRGs (Figure 5(c)).Figure 5
Transcription factor‐mediated regulatory network. Heat map (a) and a volcano plot (b) demonstrating differentially expressed TFs between hepatocellular carcinoma (HCC) and nontumor tissues, green dots represent downregulated TFs, and red dots represent upregulated TFs. (c) Regulatory network constructed based on differentially expressed TFs and survival-associated IRGs. Blue triangles represent TFs whereas orange circles represent high-risk IRGs.
(a)(b)(c)


## 3.4. Evaluation of Clinical Outcomes

In this study, we developed a prognostic signature based on the results of multivariate Cox regression analysis to divide the HCC patients into two groups with discrete clinical outcomes with regard to OS (Figure6). This index was based on the following formula: risk score = [expression level of HSPA4∗(0.029292)] + [expression level of PSME3∗(−0.04914)] + [expression level of PSMD14∗(0.054544)]+ [expression level of FABP6∗(0.054469)]] + [expression level of ISG20L2∗(0.096535)] + [expression level of TRAF3∗(0.153954) + [expression level of NDRG1∗(0.006116)] + [expression level of NRAS∗(0.040485)] + [expression level of CSPG5∗(0.238526)] + [expression level of IL17D∗(0.060712)]. The immune-based prognostic index (IRGPI) significantly stratified patients into low-risk (IRGPI < median value) and high-risk (IRGPI > median value) groups in terms of overall survival (Figure 7(a)). The area under the curve of the receiver operating characteristic (ROC) curve was 0.826, which suggested the moderate potential for the prognostic signature based on IRGs in survival monitoring (Figure 7(b)).Figure 6
Development of the prognostic index based on immune‐related genes. (a) The rank of prognostic index and distribution of groups. (b) Survival status of patients in different groups. (c) Heat map of expression profiles of these included genes.
(a)(b)(c)Figure 7
The prognostic value of the prognostic index based on immune‐related genes. (a) Patients in the high‐risk group suffered shorter overall survival. (b) Survival‐dependent receiver operating characteristic (ROC) curve validation of the prognostic value of the prognostic index. (c) Univariate and (d) multiple regression analysis of HCC.
(a)(b)(c)(d)The clinical data and risk scores were analyzed by univariate and multivariate regression analysis. TheP value of risk score was less than 0.001 in both univariate (Figure 7(c)) and multivariate regression (Figure 7(d)) analyses. These results indicated that the IRGPI obtained by our model could be used as an independent predictor after adjusting for other parameters, including age, gender, grade, tumor stage, tumor size, distant metastasis status, and the amounts of nodules (Table 2).Table 2
Univariate and multiple regression analysis of HCC.
VariablesUnivariate analysisMultivariate analysisHazard ratio (95% CI)P valueHazard ratio (95% CI)P valueAge0.995 (0.976–1.014)0.602599261.010 (0.990–1.031)0.323971151Gender (male/female)0.893 (0.532–1.500)0.669262011.555 (0.850–2.844)0.151666718Grade1.083 (0.772–1.520)0.644037011.057 (0.728–1.536)0.769568712Tumor stage2.084 (1.590–2.733)1.07E − 071.192 (0.404–3.521)0.750289146Tumor1.980 (1.541–2.543)8.94E − 081.515 (0.569–4.0340.405391572Distant metastasis4.769 (1.485–15.311)0.008669522.593 (0.637–10.553)0.183412547Node2.439 (0.593–10.035)0.216733452.311 (0.356–15.002)0.379943568Risk score1.264 (1.198–1.334)1.38E − 171.251 (1.172–1.336)2.10E − 11


## 3.5. Clinical Utility of IRGPI

To further assess the clinical value of the immune-related gene-based prognostic index (IRGPI), the relationship between this hub survival-associated IRGs and clinical characteristics including age, gender, survival state, grade, pathological stage, T stage, M stage, and N stage were analyzed (Table3). IRGPI showed a significant difference in survival state (Figure 8(a)), grade (Figure 8(b)), T stage (Figure 8(c)), and M stage (Figure 8(d)). However, no difference was observed between age, gender, and N stage.Table 3
Relationship between these survival-associated IRGs and clinical characteristics.
GenesAge (≥65/<65)Gender (male/female)State (dead/alive)Grade (IV-III/I-II)TPtPtPtPHSPA40.0240.981−1.2310.22−2.8330.006−3.4796.31E − 04PSME3−0.8870.377−0.3960.692−2.7110.008−2.3170.022PSMD14−1.0250.3070.3950.694−3.3180.001−2.0670.04FABP6−0.1430.886−1.2190.224−1.4130.162−1.7920.076ISG20L2−1.9470.0530.490.625−3.2470.002−2.8370.005TRAF30.2930.771.0480.297−2.1650.033−1.420.157NDRG10.4190.677−1.0920.276−2.5710.012−0.3560.723NRAS−1.7610.0810.0010.999−3.2140.002−2.2910.023CSPG5−0.3280.743−0.9020.368−2.3850.019−2.0120.046IL17D−1.3210.1881.4950.139−1.8960.062−2.0190.046Risk score−0.3160.753−0.2940.769−3.8612.60E − 04−2.1480.033Figure 8
The relationships between the immune-based prognostic index and (a) survival state, (b) grade, (c) T stage, (d) M stage.
(a)(b)(c)(d)


## 4. Discussion

Hepatocellular carcinoma is also known as a clear example of inflammation-related cancer, given that more than 90% of HCCs arise in the context of hepatic injury and inflammation [32]. This fact highlights the importance of the differentially expressed IRGs. Previous studies have already addressed gene expression-based prognostic signatures in hepatocellular carcinoma [33, 34], thus providing a fundamental understanding of the pathogenesis of HCC at the genetic level. However, there is no comprehensive study that explored the characteristics of IRGs in HCC. Consequently, we conducted this comprehensive, genome-wide profiling study of IRGs to explore their clinical significance and verify reliable prognostic biomarkers that could be used to select patients at the highest risk for recurrence. Bioinformatic systems make it possible to explore their molecular mechanisms more deeply.To identify the prognostic biomarkers associated with immune genes of hepatocellular carcinoma, we first identified the immune genes associated DEGs and then constructed an immune-related gene-based prognostic index (IRGPI) composed of ten genes (HSPA4, PSME3, PSMD14, FABP6, ISG20L2, TRAF3, NDRG1, NRAS, CSPG5, IL17D) by Cox regression model and finally verified its correlation with OS. The IRGPI significantly stratified patients into high-risk and low-risk groups, and high-risk groups were associated with poor prognosis. Univariate and multivariate Cox analyses confirmed that this IRGPI could be used as an independent prognostic factor. The ROC analyses proved the reliability of this IRGPI.Immune characteristics in the tumor microenvironment are essential for the development of immunotherapies and the prediction of their clinical responses in cancers [35]. Our research focused on the comparison of immunogenomic profiles between hepatocellular carcinoma and healthy liver tissue, trying to identify some potential clinical implications. Gene functional enrichment analysis and KEGG suggested that these genes are mainly involved in growth factor activity and cytokine-cytokine receptor interactions, respectively.Hepatocyte growth factor (HGF) is the first factor to stimulate hepatocyte division and regeneration [36]. It also participates in enhancing angiogenesis, immune response, cell motility, and cell differentiation [37]. The interaction between HGF and hepatocytes can enhance HGF/c-Met signal transduction [38]. The expression of c-Met in HCC was higher than in surrounding tissues. Overexpression of c-met and other oncogenes have been identified as the causes of HCC invasiveness [39]. Our results showed that the change of the immune genome could affect the occurrence of hepatocellular carcinoma through the growth factor pathway. As hepatitis virus infection is the main cause of hepatocellular carcinoma [40], inflammatory response induced by cytokine-mediated immune response is considered the most important factor in the development of hepatocellular carcinoma [41]. Our KEGG analysis showed that the key immune regulatory molecules of hepatocellular carcinoma were mainly involved in the cytokine-cytokine receptor interaction pathway.Gene-regulatory networks modulate the entire process of gene expression and protein formation in living cells and therefore determine the fate of cells. TFs regulate gene expression by translating cis-regulatory codes into specific gene-regulatory events. In this study, we explored the main regulation consisting of transcription factors (TFs) and their impacting immune-related genes. Among HCC immune-related genes, we identified the potential targets of TFs. These datasets and their regulations were used to construct a comprehensive HCC immune-related genes TF mediated regulatory network. SIRT6, CENPA, and KDM1A are prominently featured in this network. It has been reported that SIRT6 overexpression in primary HCC tumors is correlated with tumor size and grade [42], while CEAPA, combined with KIF20A, PLK1, and NCAPG, form a 4-gene expression prognostic signature, which can be used to predict prognosis and to define a subgroup of high-risk HCC patients who could potentially benefit from JmjC inhibitor therapy [43]. GNPAT overexpression induced by c-myc/KDM1A complex transcriptional activation has been confirmed to be related to the progression of HCC [44]. However, the relationship between these TFs and IRGs has not yet been confirmed. Our network is conducive to a better understanding of the potential molecular mechanism of these IRGs.Previously, Zhao et al. performed genome-wide methylation profiling of the different stages of hepatitis B virus-related hepatocellular carcinoma [45]. Zucman et al. integrated signatures to study the genetic landscape and biomarkers of HCC [46]. Deng et al. and his team analyzed tumor microenvironment-related genes of prognostic value in hepatocellular carcinoma [47]. Although several HCC signatures based on immune-related genes have been developed recently [48–53], a more complete and reliable index that can predict both survival and immunotherapy success for HCC patients is urgently needed. In this study, we developed a prognostic signature based on ten immune-related genes for hepatocellular carcinoma. Our prognostic immune signature can be used to stratify clinically defined HCC patients into subgroups with different survival outcomes and can be clarified as an immune status indicator. Interestingly, our data showed that IRGPI performed moderately in prognostic predictions and was correlated with age, tumor stage, metastasis, number of lesions, and tumor burden. We further leveraged the complementary value of molecular and clinical characteristics and showed that combining both could provide a more accurate estimation of overall survival in HCC.


## 5. Conclusion

The present study identified the immune genes associated with DEGs that were then used to construct and validate the immune-related gene-based prognostic index for predicting the outcomes of HCC patients. Further study of these immune genes associated with DEGs will provide a new understanding of the potential relationship between immune genes and HCC prognosis.



---

*Source: 3761858-2021-10-29.xml*