Identification of prognostic immune cells and potential immune-related markers in hepatocellular carcinoma.

Background: Hepatocellular carcinoma (HCC) is one of the most common and deadly tumors worldwide. Immunotherapy has emerged as a promising strategy for HCC treatment, and understanding the immune microenvironment of HCC provides a theoretical basis for identifying new immune targets. However, the roles of immune components and their regulatory mechanisms in HCC require further clarified.
Methods: By analyzing HCC expression profiles from The Cancer Genome Atlas (TCGA) database, we depicted the proportion profile of immune cells for each sample using the software CIBERSORTx. Using R packages, we also characterized the distribution of differentially expressed genes (DEGs) in immune cells, calculated the correlation coefficient between immune cells and common DEGs, and analyzed their biology function by Gene-Ontology analysis.
Results: We found that seven immune cell types were related to the overall survival of HCC patients, and identified 3,692 differentially expressed immune-related genes, predominantly functioning in nucleic acid processing and metabolism. Moreover, 14 DEGs were identified as common candidates related to immune cells and overall survival. Conclusions: Our study not only presents an overview of the immune features of the microenvironment of HCC, but also provides potential targets related to immune components. 2022 Translational Cancer Research. All rights reserved.

Introduction

Primary liver cancer, one of the most common cancers, contributes to a large number of cancer-related deaths worldwide (1,2). Hepatocellular carcinoma (HCC) accounts for the majority of cases among all types of primary liver cancers (3). The main risk factors for HCC development include liver cirrhosis, infection with hepatitis B virus (HBV)/hepatitis C virus (HCV), alcohol abuse, and metabolic syndrome (4-7). In addition, other factors, such as smoking and aflatoxin B1 intake, have been identified as contributors to HCC (8,9). HCC is frequently characterized by mutations in TERT, TP53, and ARID1A, and dysregulation of the WNT, RAS, and mTOR signaling pathways (10,11). However, not all of these onco-drivers are considered druggable targets, and pharmaceuticals for HCC treatment currently available pharmaceuticals are currently limited. Tyrosine kinase inhibitors (TKIs) such as Sorafenib, Regorafenib, and Lenvatinib have been demonstrated to have limited survival advantages for HCC (12-14). However, identification of novel druggable targets for more effective treatment of HCC has inevitably become a priority.

Recently, there has been an increase in research focusing on the tumor microenvironment (TME), in particular, its immune component. It is widely accepted that the immune components of the TME significantly influence the pathophysiological progress of tumors, therapeutic responses, and clinical outcomes (15-17). Furthermore, immunotherapies that modulate the immune response against tumors are making promising clinical progress. Antibodies blocking immune-checkpoint molecules, such as PD-1 and CTLA-4, are being tested and put into clinical application with inspiring results (18). Thus, a deeper understanding of the immune components of the TME is indispensable for developing novel strategies for HCC treatment.

To explore the role of immune components in HCC, we used TCGA database and screened immune cells that were related to overall survival in HCC patients. Next, we identified common genes that were related to candidate immune cells and found that 14 candidates were correlated with overall survival. To study the biological role of the differentially expressed genes (DEGs) related to immune components, we used gene ontology (GO) analysis to identify the enriched pathways and functions. Our study characterizes the immune features of HCC, and provides potential clues on HCC therapies. We present the following article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-19/rc).

Methods

Samples collection and data processing

Expression data with clinical information from HCC patients were downloaded from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). These data included fifty files from normal tissues and 369 individual files from tumor tissues. The Fragments Per Kilobase per Million (FPKM) values were used in subsequent analyses. The clinicopathologic features of HCC patients were shown in and website: https://cdn.amegroups.cn/static/public/tcr-22-19-01.zip. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Since the data from the TCGA were de-identified and publicly available, no institutional review board approval was necessary and no informed consent was signed for this study.

Table 1

Clinicopathologic parameters of HCC patients in TCGA cohorts

Clinicopathologic parameters	TCGA
	Total (n=369)	%
Age, years
<60	169	45.80%
≥60	199	53.93%
Not reported	1	0.27%
Gender
Male	249	67.48%
Female	120	32.52%
T
T1	181	49.05%
T2	94	25.47%
T3	78	21.14%
T4	13	3.52%
Tx	1	0.27%
Not reported	2	0.54%
N
N0	250	67.75%
N1	4	1.08%
Nx	114	30.89%
Not reported	1	0.27%
M
M0	265	71.82%
M1	4	1.08%
Mx	100	27.10%
Stage
I	171	46.34%
II	86	23.31%
III	83	22.49%
IV	5	1.36%
Not reported	24	6.50%

HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.

Identifying immune cells related to survival

We calculated the proportions of 22 types of immune cells using software CIBERSORTx. First, we uploaded the matrix of FPKM values for all samples to the CIBERSORTx website (https://cibersortx.stanford.edu/index.php) as the Mixture file. Next, we selected the LM22 matrix (default setting) within the software as the signature gene file, and ran the software in the relative mode, with batch correction and quantile normalization both disabled. We set the permutations for significance analysis to 100. We then identified survival-related immune cells using the survival R package based on the proportion of each cell type. We used the log-rank test to analyze the survival data and set the medians of the proportions of cell types as the cut-off values. We analyzed the relationship between the proportions of immune cells and clinical stages using the one-way ANOVA.

Identification of immune-related genes and common genes

The specific genes for each type of survival-related immune cell were identified by calculating the DEGs between the high- and low-expressing samples using the median of proportions as the cut-off value. DEGs were analyzed using the limma R package. Benjamini & Hochberg (BH) method for multiple hypothesis testing on the P values were performed in differential expression analysis. Significantly altered genes for each cell type included those with |log2(FoldChange (FC))| >2, and P value <0.05. The results were visualized using the ggplot2 R package. Common genes were screened and visualized by UpSet R package.

Enrichment analysis of immune-related genes and common genes

Significant Gene Ontology (GO) enrichment analyses were performed using the cluster profiler R package. A P value <0.05, and |ave_log2FC| >2, were chosen as the thresholds.

Relationship between immune components and common genes

We calculated the Pearson correlation between the common genes and 22 types of immune components, and the result was visualized by the ggplot2 R package. Additionally, we calculated the correlation between the common genes and the overall survival by the Log-Rank test.

Statistical analysis

Method of one-way analysis of variance (ANOVA) was performed to compare differences in immune cell components among stages, while Tukey’s multiple comparisons test was performed to compare differences in components of immune cells between any two of stages.

Results

Identifying immune cells related to HCC overall survival

It is widely reported that immune infiltration is related to the occurrence, progression, and prognosis of HCC. To further explore the relationship between immune components and prognosis in HCC, the proportions of 22 types of immune cells were calculated using the online software CIBERSORTx () (19). The proportion of macrophages was 39.9%±11.8% and T cells was 36.28%±12.8%, which accounted for the majority of immune cells ().

Figure 1

The relationship between the proportion of immune cells and overall survival in HCC. (A) The profile of the proportion of immune cells in the HCC samples. Each sample is represented by one column, and the proportion of the specific immune cell types is indicated by the column height in that sample. (B) The proportion of 22 types of immune cells in the HCC samples. Medians and quartiles are shown in black lines. (C,D) Survival analysis for the proportion of (C) follicular helper T (Tfh) cells, M0 macrophages, eosinophils, (D) naïve B cells, resting memory CD4+ T cells, γδT cells, and resting mast cells. HCC, hepatocellular carcinoma.

By analyzing the survival data of HCC patients in TCGA cohorts, we identified the proportion of M0 macrophages (OS: HR =0.4688, 95% CI: 0.3284 to 0.6692, logrank P<0.0001), naïve B cells (OS: HR =1.579, 95% CI: 1.107 to 2.254, logrank P=0.011), resting mast cells (OS: HR =1.577, 95% CI: 1.109 to 2.242, logrank P=0.011), CD4+ T cells (memory resting) (OS: HR =1.494, 95% CI: 1.055 to 2.116, logrank P=0.023), follicular helper T (Tfh) cells (OS: HR =0.6425, 95% CI: 0.4536 to 0.9101, logrank P=0.012), γδT cells (OS: HR =1.648, 95% CI: 1.146 to 2.369, logrank P=0.0064), and eosinophils (OS: HR =0.5744, 95% CI: 0.3284 to 1.005, logrank P=0.049) as candidate immune cells related to overall survival (). Higher proportions of Tfh cells, M0 macrophages and eosinophils correlated with worse prognosis, while higher proportions of naive B cells, CD4+ T cells (memory resting), γδT cells, and resting mast cells correlated with better prognosis (). Thus, we considered these six types of immune cells as survival-related cells, and as immune candidates for further investigation.

Exploring the relationship between immune cells and clinical stages of HCC

Next, we investigated the relationship between the proportions of each type of immune cells at different clinical stages. The proportions of memory resting CD4 T cells were statistically different in distinct clinical stages (P<0.05), while the proportion of M0 macrophages was significantly lower in stage I than in stages II (P<0.05) and III (P<0.05), and the proportion of M1 macrophages was significantly higher in stage I than in stage III (P=0.0431). However, the proportions of most types of immune cells showed no statistical difference in the distinct clinical stages of HCC (). These results suggest that there is no clear relationship between the proportion of immune cells, including survival-related immune cells, and clinical stages.

Figure 2

The relationship between immune cells and HCC stages. The proportions of 22 types of immune cells in different clinical stages of HCC are shown in violin plots. Medians are shown in black lines. Quartiles are shown in dotted lines. *, P<0.05; **, P<0.01; ***, P<0.001. HCC, hepatocellular carcinoma.

Identification of DEGs related to immune cells and pathway analysis

Gene expression relating to the specific immune cell candidates was assessed as described in the ‘Methods’ section, revealing a total of 3,692 DEGs. Among these genes, 1,531 were related to naïve B cells, 2,321 to M0 macrophages, 1,277 to Tfh cells, 1,447 to memory resting CD4+ T cells, 797 to γδT cells, 1,683 to resting mast cells, and 1,283 to eosinophils ().

Figure 3

Identification of DEGs related to immune cells and GO pathway analysis. (A-G) Gene expression profiles related to naive B cells, M0 macrophages, Tfh cells, memory resting CD4+ T cells, γδT cells, and resting mast cells. Volcano plots are used to present the data. The upregulated/downregulated genes are presented as red/blue dots according to the criteria: |log2Fold Change| >2 and P value <0.05. (H-I) GO pathway enrichment analysis on genes related to each type of immune cells. The significance of enrichment results is indicated by the color, and the count of genes enriched for each result is indicated by the dot size. DEGs, differentially expressed genes; GO, Gene Ontology.

Next, to explore the biological function of the DEGs, we performed Gene Ontology (GO) pathway enrichment analysis for each group of immune-related DEGs. The results of top 30 pathways are shown in . The DEGs were mainly enriched in the regulation of RNA-related processes (like metabolism and splicing), immune responses (like T cell activation, neutrophils activation, and immune-related cell surface receptor signal pathway), and mitosis-related processes ( and website: https://cdn.amegroups.cn/static/public/tcr-22-19-01.zip). The results of the GO pathway enrichment analysis therefore suggest that RNA-related signal pathway may play a pivotal role in regulating immune components of HCC.

Identification and analysis of common genes related to immune cells

We identified the distribution of DEGs in the overall survival-related immune cells to explore the common genes which may play important roles in regulating immune cells (). Sixteen genes were differentially expressed in all seven types of immune cells, including MFSD10, DBN1, ALDOA, ELOVL1, CBX3, HLF, HMGA1, TPD52L2, WDR1, SLC25A6, SYNGR2, ENO1, G6PD, COLGALT1, MANBAL, and JPT2. The expression pattern of 16 genes is shown in , and it was shown that the expression pattern of HLF seemed to be negatively correlated to the other 15 genes. Thus, 16 genes were considered as candidate genes of immune cells in HCC.

Figure 4

Identification of common genes related immune cells. (A) The distribution profile of DEGs in the seven types of immune cells. (B) The expression profiles of common genes in HCC-TCGA cohort. DEGs, differentially expressed genes; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.

Exploring the relationship between common genes and clinical features of HCC

Relationship between overall survival and 16 candidates were explored, and 14 out of 16 genes (ALDOA, CBX3, COLGALT1, ELOVL1, ENO1, G6PD, HLF, HMGA1, JPT2, MANBAL, MFSD10, SLC25A6, SYNGR2, and TPD52L2) were shown to be correlated to overall survival significantly (). Higher expression of ALDOA (OS: HR =0.5686, 95% CI: 0.4008 to 0.8066, logrank P=0.0013), CBX3 (OS: HR =0.5425, 95% CI: 0.3811 to 0.7723, logrank P=0.00057), COLGALT1 (OS: HR =0.6135, 95% CI: 0.4323 to 0.8706, logrank P=0.0057), ELOVL1 (OS: HR =0.5433, 95% CI: 0.3827 to 0.7715, logrank P=0.00054), ENO1 (OS: HR =0.4941, 95% CI: 0.3465 to 0.7047, logrank P<0.0001), G6PD (OS: HR =0.5214, 95% CI: 0.3672 to 0.7404, logrank P=0.00021), HMGA1 (OS: HR =0.5283, 95% CI: 0.3717 to 0.7509, logrank P=0.0003), JPT2 (OS: HR =0.6893, 95% CI: 0.4878 to 0.9739, logrank P=0.034), MANBAL (OS: HR =0.6095, 95% CI: 0.4304 to 0.8630, logrank P=0.0048), MFSD10 (OS: HR =0.6126, 95% CI: 0.4325 to 0.8676, logrank P=0.0053), SLC25A6 (OS: HR =0.7061, 95% CI: 0.4984 to 1.000, logrank P=0.049), SYNGR2 (OS: HR =0.5931, 95% CI: 0.4173 to 0.8428, logrank P=0.0032), and TPD52L2 (OS: HR =0.5947, 95% CI: 0.4193 to 0.8435, logrank P=0.0032) correlated with worse prognosis, while higher expression of HLF (OS: HR =1.428, 95% CI: 1.007 to 2.025, logrank P=0.045) correlated with better prognosis.

Figure 5

The relationship between common genes and overall survival. The survival analysis for the common gene ALDOA, CBX3, COLGALT1, ELOVL1, ENO1, G6PD, HLF, HMGA1, JPT2 MANBAL, MFSD10, SLC25A6, SYNGR2, and TPD52L2.

Also, the relationship between prognostic candidates and clinical stages were explored. We found that the expression levels of ALDOA and MFSD10 increased with the increase of clinical stages. Meanwhile, the levels of HLF decreased with the increase of clinical stages. Also, the level of CBX3, COLGALT1, JPT2, MANBAL, SYNGR2, and TPD52L2 increased with clinical stage I to III. The levels of SLC25A6 increased in stage IV compared to stage I. The level of ELOVL1, ENO1, G6PD, and HMGA1 showed no difference in distinct clinical stages ().

Figure 6

The relationship between common genes and HCC traits. The expression level of common genes in different clinical stages of HCC were shown in box plots. Medians and quartiles were shown in black lines. HCC, hepatocellular carcinoma.

Exploring the relationship between candidate genes and immune cells

Additionally, we explored the correlation between candidate genes and the proportions of 22 types of immune cells, and found that candidate genes were mainly correlated to the proportions of naïve B cells, Tregs, γδT cells, resting NK cells, M0 macrophages, and M1 macrophages (). For example, the proportions of naïve B cells were negatively correlated to expression levels of 15 candidate genes (except for HLF), while the proportions of M0 macrophages were positively correlated to expression levels of 15 candidate genes (except for HLF). In summary, it is indicated that there were potential relationship between candidate genes and the function of immune cells to further exploration.

Figure 7

The relationship between common gene and immune components. The heatmap represents the Pearson correlation coefficients between common genes and immune components. *, P<0.05.

Discussion

Since the molecular and immune features of HCC have been comprehensively characterized (20,21), the prognostic value of immune-based markers has emerged. In this study, we aimed to identify the prognostic markers basing on key immune components of HCC. The key immune components of HCC were analyzed by calculating the proportions of 22 types of immune cells and the relationship with overall survival of HCC patients. The results showed a diversification of the proportions of immune cells among patients, indicating an important role and the complexity of the immune component of HCC.

In our study, seven types of immune cells were found to be related to HCC survival, including naive B cells, Tfh, CD4+ T Cells (memory resting), γδT cells, M0 macrophages, eosinophils, and resting mast cells. An increasing number of studies have identified specific immune cells as prognostic markers in HCC, including B cells, T cells, macrophages, and mast cells. Zhang et al. analyzed the distribution of B cells in HCC patients by multiplexed IHC and cytometry, and found that a high density of naïve B cells was a potential prognosticator for survival (22). Dysfunction of circulating Tfh cells has been reported to lead to HCC progression, and Tfh cells have potential as prognostic markers and a novel therapeutic targets (23,24). These findings, which were based on circulating Tfh, contradict with our results that a higher proportion of Tfh cells were associated with worse prognosis, therefore the role of infiltrating Tfh cells needs further investigation. These results suggest that these seven types of immune cells could be predictors and potential targets of HCC prognosis.

We also analyzed the DEGs related to immune cells and their biological functions, and found that DEGs were mostly enriched in RNA-related processes (like RNA splicing and RNA metabolism), immune responses (like T cell activation, neutrophils activation, and regulation of immune-related cell surface receptor signal pathway), and mitosis processes. It is implied that multiple processes especially RNA-related processes may play important roles in immune microenvironment of HCC, and deserve to explore their therapeutic potential. For instance, it has widely reported that RNA N6-methyladenosine, which is a modification pattern of mRNA and involves in mRNA stability splicing and translation, plays a crucial part in HCC carcinogenesis and drug resistance (25,26).

Additionally, 16 common genes were identified from immune-related DEGs, and 14 of them (ALDOA, CBX3, COLGALT1, ELOVL1, ENO1, G6PD, HLF, HMGA1, JPT2, MANBAL, MFSD10, SLC25A6, SYNGR2, and TPD52L2) are significantly correlated with overall survival of HCC patients. HMGA1, a nuclear protein associated with chromatin remodeling, is overexpressed in multiple cancers, including thyroid cancer, liver cancer, and pancreatic cancer (27-29). Wei et al. reported that HMGA1, activated by PD-L1, promoted colorectal cancer stem cells expansion (30). In addition, several studies have suggested that ALDOA, ENO1, HLF, and G6PD are associated with cancer progression and immune response (31-34). It is generally accepted that macrophages are the main component in microenvironment of HCC, and play important roles in promoting or inhibiting hepatocarcinogenesis under particular circumstances (35). Interestingly, the proportions of M0 macrophages were significantly correlated to 16 candidate genes, indicating that M0 macrophage are associated with multiple biological aspects including energy metabolism, chromatin processing and RNA processing, and their potential relationships deserve further investigations. These results suggest that candidate genes may play a significant role in regulating immune components of HCC, and deserve further exploration the underlying mechanisms.

Our study has several limitations. First, the proportions of immune cells were calculated based on the TCGA dataset, so further validation will be required using other resources. Second, the relationship between immune cells and prognostic candidates were calculated using mathematical methods. These results should be confirmed by molecular biology experiments such as RNA-interfering and proliferation assay. Third, the proportions of Tfh cells and γδT cells were zero in many cases, suggesting that the results relating to these two immune cell-types may rely on outliers. Therefore, conclusions regarding Tfh cells and γδT cells should be made with caution. Fourth, by the reason of limited information of viral load, the relationship between hepatitis virus infection and immune markers or common candidates were not explored further in our study.

In summary, we identified seven types of survival-related immune cell candidates and 16 common DEGs, and 14 of which were correlated with overall survival in HCC patients. These cells and candidate genes could be considered as prognostic markers, or potential targets for HCC treatment. Our research not only provides insights into immune cells of HCC, but also identifies potential therapeutic targets worth further exploration.

Conclusions

Our study not only presents an overview of the immune features of the microenvironment of HCC, but also provides potential targets related to immune components.

The article’s supplementary files as