Identification of Hub Genes Associated With Immune Infiltration and Predict Prognosis in Hepatocellular Carcinoma via Bioinformatics Approaches

Abstract

Aims:

In the cancer-related research field, there is currently a major need for a greater number of valuable biomarkers to predict the prognosis of hepatocellular carcinoma (HCC). In this study, we aimed to screen hub genes related to immune cell infiltration and explore their prognostic value for HCC.

Methods:

We analyzed five datasets (GSE46408, GSE57957, GSE74656, GSE76427, and GSE87630) from the Gene Expression Omnibus database to screen the differentially expressed genes (DEGs). A protein–protein interaction network of the DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes; then, the hub genes were identified. Functional enrichment of the genes was performed on the Metascape website. Next, the expression of these hub genes was validated in several databases, including Oncomine, Gene Expression Profiling Interactive Analysis 2 (GEPIA2), and Human Protein Atlas. We explored the correlations between the hub genes and infiltrated immune cells in the TIMER2.0 database. The survival curves were generated in GEPIA2, and the univariate and multivariate Cox regression analyses were performed using TIMER2.0.

Results:

The top ten hub genes [DNA topoisomerase II alpha (TOP2A), cyclin B2 (CCNB2), protein regulator of cytokinesis 1 (PRC1), Rac GTPase-activating protein 1 (RACGAP1), aurora kinase A (AURKA), cyclin-dependent kinase inhibitor 3 (CDKN3), nucleolar and spindle-associated protein 1 (NUSAP1), cell division cycle-associated 5 (CDCA5), abnormal spindle microtubule assembly (ASPM), and non-SMC condensin I complex subunit G (NCAPG)] were identified in subsequent analysis. These genes are most markedly enriched in cell division, suggesting their close association with tumorigenesis. Multi-database analyses validated that the hub genes were upregulated in HCC tissues. All hub genes positively correlated with several types of immune infiltration, including B cells, CD4⁺ T cells, macrophages, and dendritic cells. Furthermore, these hub genes served as independent prognostic factors, and the expression of these hub genes combing with the macrophage levels could help predict an unfavorable prognosis of HCC.

Conclusion:

In sum, these hub genes (TOP2A, CCNB2, PRC1, RACGAP1, AURKA, CDKN3, NUSAP1, CDCA5, ASPM, and NCAPG) may be pivotal markers for prognostic prediction as well as potentially work as targets for immune-based intervention strategies in HCC.

Introduction

Hepatocellular carcinoma (HCC), the second leading cause of cancer-related death in the world, is a commonly fatal cancer with an unfavorable prognosis due to its complex genetics and clinical features (Siegel et al., 2015; Villanueva, 2019). To a certain extent, high heterogeneity contributes to the low survival rate of HCC despite excision, transplantation, transcatheter arterial chemoembolization, and radiofrequency ablation, among others, having been widely used in HCC treatment (Imamura et al., 2003; Forner et al., 2018; Caruso and Nault, 2019; Yin et al., 2019). Timely and effective intervention for HCC patients can improve not only their quality of life but also their survival rate (Deng et al., 2020; Kim et al., 2020). Therefore, the identification of new prognostic biomarkers and therapeutic targets plays a crucial role in HCC therapy.

Several prognostic biomarkers have been widely applied in HCC, such as alpha-fetoprotein and des-gamma-carboxyprothrombin (Fox et al., 2014; Abe et al., 2017). Nevertheless, these markers depend on the significant burthen of a tumor, which has resulted in their often limited application and inconsistent performance assessments (Park and Park, 2013). Numerous studies have demonstrated that valuable prognostic molecules can be identified from the bioinformatics analysis of high-throughput data, such as functional genes (Ma et al., 2020; Wang et al., 2020; Xie et al., 2020). From this, it can be determined that immune-associated genes may play a crucial role in HCC outcomes and targeted therapies on immune cells; thus, related genes have been developed for the reactivation of adaptive and innate immune systems and the creation of a strong antitumoral immune response. For instance, some researchers have found that inhibitors of programed death-1, programed death-ligand 1, and cytotoxic T lymphocyte-associated antigen 4 produce antitumoral effects on HCC cells (Langhans et al., 2019; Wang J. et al., 2019; Zhang Y. et al., 2020). Unfortunately, only sectional HCC patients with determinate tumor types react to the current immunotherapies; therefore, there is an urgent need to identify more underlying immune targets.

In the present study, differentially expressed genes (DEGs) were screened from the Gene Expression Omnibus (GEO) database. Then, we performed enrichment and protein–protein interaction (PPI) analyses of these genes to comprehend the functions of DEGs and identify the top ten hub genes in HCC. Next, we explored the potential correlations between each of the hub genes and infiltrated immune cells in the TIMER2.0 database. Furthermore, we visualized the prognostic landscape of candidate hub genes using several databases, including Oncomine (Rhodes et al., 2007), Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (Tang et al., 2019), Human Protein Atlas (HPA) (Ponten et al., 2008), and TIMER2.0 (Li et al., 2020).

Materials and Methods

Data Source

The gene expression datasets (GSE46408, GSE57957, GSE74656, GSE76427, and GSE87630) of HCC were obtained from GEO¹. All of the data included in the present study was available online. Information on these five datasets is summarized in Table 1.

DEGs Processing

GEO2R², an interactive online tool that can compare two different groups in a GEO dataset, was applied to screen the DEGs (Davis and Meltzer, 2007). The DEGs were defined as different expression genes between tumor and tumor-adjacent tissues in HCC patients with an adjusted p value < 0.05 and an absolute log fold-change (FC) > 1. Accordingly, to decrease the false discovery rate, the p values were adjusted using the Benjamini and Hochberg method. The overlapping up- and downregulated DEGs from these five datasets were identified using TBtools software (Chen et al., 2020).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analyses

Functional enrichment analyses played a crucial role in the identification of biological characteristics in transcriptome data. In this study, the Metascape database was used to conduct the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses in overlapping DEGs and identified hub genes (Zhou Y. et al., 2019). In order to choose the more remarkable terms within each cluster, a p value < 0.01 was considered statistically significant. We used bar charts and bubble diagrams, respectively, to visualize the results of the GO and KEGG analyses.

PPI Network Construction and Module Analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) is an interacting gene database designed to analyze PPI information (Szklarczyk et al., 2019). The overlapping DEGs were mapped in STRING to generate a network with functional interactions; then, this PPI network was visualized using Cytoscape software. Next, we employed the Molecular Complex Detection (MCODE) plugin to determine the most significant gene modules. Moreover, the cytoHubba plugin was applied to identify the hub genes using the Maximal Clique Centrality (MCC) method, and the top ten hub genes with the highest MCC scores were subjected to the subsequent analyses.

Validation for mRNA and Protein Levels of Hub Genes in Multi-Databases

The Oncomine database³ is a publicly available cancer database that facilitates the analysis of genome-wide expression in multifarious cancers. In the present study, the Oncomine was utilized to overview the mRNA expression of candidate genes with a p value < 0.0001 and |FC| > 1.5.

The GEPIA2 database⁴, which includes TCGA and GTEx data, was applied to analyze the differential expression of the hub genes in the HCC and normal groups, and the cutoff values were set as |log2FC| = 1.0 and p value = 0.01.

Furthermore, immunohistochemistry analysis was conducted online to observe the distribution and protein level of the candidate hub genes in the HPA database⁵.

Survival Analysis of Hub Genes

We employed GEPIA2 to perform overall survival (OS) and relapse-free survival (RFS) analyses, with a median group cutoff in 362 HCC patients. The survival curves with the calculated hazard rate (HR) and the log-rank p value were presented on the charts. Additionally, TIMER2.0 was utilized to conduct a univariate Cox analysis to validate the results from the GEPIA2 analyses.

TIMER2.0 Database Analysis

TIMER2.0⁶ is a comprehensive source that can explore the relationship between two genes, and the correlations of mRNA expression and immune infiltration. In this study, we analyzed the correlations of hub gene expression and several tumor-associated immune cells, including B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells (DCs). A p value less than 0.01 was considered statistically significant to identify the more prominent correlation between the hub genes and immune cells.

Further OS analyses were performed with macrophage and single hub gene expression. Moreover, we constructed ten multivariate Cox proportional hazard models, each of which comprised seven variables, including age, tumor stage, gender, race, tumor purity, macrophage level, and expression of the single candidate hub gene. The survival curves, featuring patterns of single gene expression, and macrophage levels were shown on the diagram. The association between each macrophage and OS was displayed as the low or high expression of a single hub gene.

Statistical Analysis

The Student’s t test or non-parametric Mann–Whitney test was utilized to compare the two independent groups, as appropriate. The correlations between the candidate hub genes as well as the relationship of these genes and immune cell infiltration were assessed using Spearman’s correlation. The log-rank test was used to calculate the HR and log-rank p value to compare the survival curves. Univariate and multivariate Cox regression models were applied to calculate the HR and Cox p value. If not specifically stated, p values < 0.05 were considered statistically significant.

Results

Identification and Enrichment Analysis of DEGs

After screening the DEGs according to the criteria, 1,414, 417, 454, 493, and 1,163 genes were identified from GSE46408, GSE57957, GSE74656, GSE76427, and GSE87630, respectively (Table 1). These genes shared 107 DEGs among these five datasets, of which, there were 18 upregulated and 89 downregulated genes (Figures 1A,B and Supplementary Figure 1).

FIGURE 1

As shown in Figure 2A, the results of the GO analysis suggested that the overlapping DEGs were principally enriched in the biological process (BP), especially in the detoxification of copper ion, the organic acid catabolic process, and the monocarboxylic acid metabolic process. In terms of cellular component (CC), the DEGs were enriched in plasma lipoprotein particles and extracellular matrix. In regard to molecular function (MF), the identified DEGs were significantly enriched in the aspects of peptide hormone binding, insulin-like growth factor binding, and carboxylic acid binding. Also, KEGG analysis showed that the overlapping DEGs were dramatically concentrated in mineral absorption and tryptophan metabolism (Figure 2B).

FIGURE 2

PPI Network Construction and Hub Gene Identification

To seek the interactions of the overlapping DEGs, a PPI network, which included 84 nodes and 191 edges, was constructed and visualized in the Cytoscape (Figure 3A). As Figures 3B–D shows, the three most prominent subnetworks were identified using the MCODE plugin according to the connective degrees. Moreover, the CytoHubba plugin was used to determine the top ten hub genes based on their MCC scores (Figure 3E). Interestingly, the hub-gene network from the CytoHubba analysis completely coincided with the highest score module from the MCODE analysis, both of which included the same 10 nodes and 45 edges. Notably, DNA topoisomerase II alpha (TOP2A) was the most significant gene, with the highest MCC score of 362,912, followed by cyclin B2 (CCNB2) (MCC score = 362,906), protein regulator of cytokinesis 1 (PRC1) (MCC score = 362,904), Rac GTPase-activating protein 1 (RACGAP1) (MCC score = 362,904), aurora kinase A, (AURKA) (MCC score = 362,888), cyclin-dependent kinase inhibitor 3 (CDKN3) (MCC score = 362,886), nucleolar and spindle-associated protein 1 (NUSAP1) (MCC score = 362,880), cell division cycle-associated 5 (CDCA5) (MCC score = 362,880), abnormal spindle microtubule assembly (ASPM) (MCC score = 362,880), and non-SMC condensin I complex subunit G (NCAPG) (MCC score = 362,880). All the hub genes were upregulated in HCC tissues. In addition, the analyses of the correlations between candidate hub genes on the mRNA level were conducted on the TIMER2.0. Of these ten genes, every two genes showed a significant correlation (Spearman’s rho value > 0.7; p < 0.05) with or without tumor purity adjustment (Supplementary Table 1).

FIGURE 3

The high scores of these candidate hub genes indicated that there would be potential biological effects in the hub genes network; thus, we further determined the functional enrichment of these genes. The BP analysis proved that the hub genes were dramatically enriched in terms of cell division, the positive regulation of mitotic nuclear division, female gamete generation, and mitotic cell cycle phase transition. Besides, these genes were significantly enriched in terms of the mitotic spindle in the CC analysis and in terms of protein kinase binding in the MP analysis (Supplementary Figure 2). These results implied that these hub genes are closely associated with tumorigenesis.

The Expression of Hub Genes Was Upregulated in Multi-Databases

To verify the dependability of the results from the bioinformatics analysis, we next determined the mRNA expression levels of the hub genes in the Oncomine and GEPIA2 databases. The results showed that the transcriptional levels of TOP2A, CCNB2, PRC1, RACGAP1, AURKA, CDKN3, NUSAP1, CDCA5, ASPM, and NCAPG were significantly overexpressed in HCC tissue when compared with the normal controls (Figures 4A–K), indicating their potential oncogenic effects.

FIGURE 4

In the HPA database analysis, we found that the protein levels of TOP2A, PRC1, RACGAP1, AURKA, NUSAP1, and CDCA5 were significantly higher in the HCC tissues when compared to normal liver tissues (Figures 5A–F).

FIGURE 5

The mRNA Levels of Hub Genes Are Positively Associated With Immune Infiltration in HCC

Numerous studies have demonstrated that the infiltration of tumor-associated immune cells and immune-related genes is correlated with the development and prognosis of HCC (Duan et al., 2019; Huang et al., 2020; Tang et al., 2020). Remarkably, targeting on these immune cells and/or genes has been a prospective workable approach in HCC therapy (Jayant et al., 2020). In this study, we attempted to explore the relationship between the mRNA expression of hub genes and immune cell infiltration using TIMER2.0. As Table 2 presents, each of the hub genes correlated with tumor purity in HCC tissues. Notably, we observed that these ten genes presented significant associations with infiltrating levels of B cells, CD4⁺ T cells, macrophages, and DCs, out of which these genes most strongly correlated with B cells [correlation coefficient (COR), 0.363 to 0.451; p < 0.01], macrophages (COR 0.262 to 0.384; p < 0.01), and DCs (COR, 0.449 to 0.548; p < 0.01), indicating that these hub genes were positively related to tumor-associated B cells, macrophages, and DCs in the HCC microenvironment.

The Overexpression of Hub Genes Predicts Poor Prognosis in HCC

In the GEPIA2 analysis, the KM plotter analyses showed that the upregulated mRNA levels of TOP2A, PRC1, RACGAP1, AURKA, NUSAP1, CDCA5, ASPM, and NCAPG were correlated with worse OS in HCC patients; however, there was no significant correlation between CCNB2/CDKN3 and HCC prognosis (Figures 6A–J). For the RFS analyses, the overexpression of all ten hub genes could predict unfavorable prognosis in HCC (Figures 7A–J). Additionally, we performed a univariate Cox regression analysis of the candidate hub genes in TIMER2.0. Table 3 summarizes the validation of the prognostic values of these hub genes in the TIMER2.0 database, suggesting that each of the hub genes may be an independent risk factor in HCC.

FIGURE 6

FIGURE 7

The Overexpression of Hub Genes Accompanied by a High Level of Infiltrated Macrophages Predicts Poor Prognosis in HCC

Tumor-associated macrophages (TAMs) majorly facilitate tumor angiogenesis, invasion, and metastasis and lead to a poor prognosis in HCC (Zhang et al., 2019; Zhao et al., 2020). There have been studies substantiating that targeting TAMs becomes a potential individualized precision single or combined therapy (Li et al., 2019). Therefore, the identification of TAM-related genes will contribute to providing more potential targets of the individualized precision treatment and improve the prognosis of HCC. In the present study, we evaluated the prognostic efficiency of the combination of TAMs and expression patterns for the single hub gene. As shown in Figures 8A–J, there was no significant relationship between the TAMs and prognosis under the low expression level of TOP2A/CCNB2/PRC1/RACGAP1/AURKA/CDKN3/NUSAP1/CD CA5/ASPM/NCAPG. However, under high CCNB2 expression, higher macrophage levels had a worse outcome in HCC [HR = 1.6, p = 0.0366]. Similarly, the high macrophage level predicted unfavorable prognosis under the high expression of RACGAP1 (HR = 1.79, p = 0.0131), AURKA (HR = 1.64, p = 0.0283), CDKN3 (HR = 1.94, p = 0.0066), ASPM (HR = 1.72, p = 0.0177), and NCAPG (HR = 1.66, p = 0.0304).

FIGURE 8

Moreover, we further established ten multivariate Cox regression analyses, each of which included seven variables: macrophage level, age, stage, gender, race, tumor purity, and expression of a single candidate gene (Figures 9A–J). The results showed that, after adjustments of age, stage, gender, race, and tumor purity, there was still no statistical correlation between TAM and prognosis with the low expression of TOP2A/CCNB2/PRC1/RACGAP1/AURKA/CDKN3/NUSAP1/CD CA5/ASPM/NCAPG; nevertheless, the lower level of TAM could predict favorable prognosis under the high expression of TOP2A (HR = 2.1, p = 0.0078)/CCNB2 (HR = 2.2, p = 0.0046)/PRC1 (HR = 1.92, p = 0.0221)/RACGAP1 (HR = 2.09, p = 0.0091)/AURKA (HR = 2.08, p = 0.0061)/ CDKN3 (HR = 2.7, p = 0.0009)/NUSAP1 (HR = 2.25, p = 0.0043)/CDCA5 (HR = 2.05, p = 0.0086)/ASPM (HR = 2.18, p = 0.0053)/NCAPG (HR = 2.25, p = 0.0043) (Figures 10A–J). These results suggested that each of the hub genes was an independent unfavorable prognostic biomarker and that combining their respective expression levels with the TAMs would help them play a more effective role in the prognosis prediction of HCC.

FIGURE 9

FIGURE 10

Discussion

Hepatocellular carcinoma, one of the malignant cancers with high heterogeneity, is frequently diagnosed in its middle and advanced stages (Friemel et al., 2015; Buczak et al., 2018). Surgical resection remains the most crucial technique for HCC treatment; however, its therapeutic effects are always unsatisfactory (Abe et al., 2017; Deng et al., 2020). Thus, there is a need to screen novel carcinoma biomarkers and therapeutic targets. In the present study, comprehensive and bioinformatics analyses of multi-databases were applied to determine ten hub genes that appeared to be correlated with infiltrated immune cells in HCC. These genes were identified as independent prognostic factors in HCC patients.

In our study, five GEO datasets (GSE46408, GSE57957, GSE74656, GSE76427, and GSE87630) shared 107 common DEGs in HCC tissues; following this, the PPI network was constructed based on these genes. The results revealed a most significant module using the MCODE analysis that completely coincided with the subnetwork identification from the CytoHubba analysis. Particularly, ten hub genes in the module were upregulated in HCC tissues: TOP2A, CCNB2, PRC1, RACGAP1, AURKA, CDKN3, NUSAP1, CDCA5, ASPM, and NCAPG, respectively. The enrichment analyses presented that these hub genes were dramatically enriched in several terms of the BP analysis, including cell division, the positive regulation of mitotic nuclear division, female gamete generation, and the mitotic cell cycle phase transition. These genes were also significantly enriched in terms of the mitotic spindle in the CC analysis and in terms of the protein kinase binding in the MP analysis, suggesting that there is a close association between the hub genes and tumorigenesis. The validation in Oncomine and GEPIA2 confirms that the mRNA levels of TOP2A, CCNB2, PRC1, RACGAP1, AURKA, CDKN3, NUSAP1, CDCA5, ASPM, and NCAPG were significantly overexpressed in the HCC tissues, and at this validation, the p < 0.0001 and p < 0.01 were set in Oncomine and GEPIA2, respectively to more accurately identify the expression pattern of hub genes between HCC and normal tissues. Further HPA analysis also demonstrated that, compared to normal liver tissues, HCC tissues had significantly higher protein levels of TOP2A, PRC1, RACGAP1, AURKA, NUSAP1, and CDCA5, while we could not obtain the protein expression of CCNB2, CDKN3, ASPM, and NCAPG of HCC from the HPA website. These hub genes were validated to be closely correlated with infiltrated immune cells using the TIMER2.0 database. Both survival curves and univariate Cox regression analyses suggested that these candidate hub genes have a strong predictive ability for HCC. Previous studies demonstrated that TAMs extremely facilitate tumor angiogenesis and lead to a detrimental prognosis in HCC (Zhang et al., 2019; Zhao et al., 2020). The identification of TAM-related genes will facilitate providing more potential targets of the individualized precision treatment and improve the prognosis of HCC. Thus, we explored the prognostic value for the combination of candidate gene and TAM expression patterns in HCC and found that there was no significant relationship between TAMs and prognosis under a low expression level of the single hub gene. Meanwhile, under the high expression of CCNB2/RACGAP1/AURKA/CDKN3/ASPM/NCAPG, high TAM levels predicted unfavorable prognosis. Furthermore, the multivariate Cox regression models indicated that all candidate hub genes were independent predictors and that combining their respective expression levels with TAM will help them play a more effective role in the prognosis prediction of HCC.

In the enrichment analysis of candidate hub genes, we observed that nine out of ten genes were significantly associated with cell division: TOP2A, CCNB2, PRC1, RACGAP1, AURKA, NUSAP1, CDCA5, ASPM, and NCAPG. The remaining gene, CDKN3, has a dual function in the regulation of the cell cycle. On one hand, CDKN3 serves as a cyclin-dependent kinase inhibitor while interacting with and dephosphorylating CDK2 kinase, thereby restraining its activation (Hannon et al., 1994; Johnson et al., 2002); on the other hand, CDKN3 can act as a MDM2-binding protein that forms a complex with MDM2 and P53, thus suppressing the production of P21, leading to the acceleration of cell cycle progression (Okamoto et al., 2006). These candidate hub genes have been demonstrated to work as oncogenes and are associated with clinical prognosis in numerous solid neoplasms, particularly in HCC (Roy et al., 2018; Gong et al., 2019; Lin J. et al., 2019). The TOP2A gene encoded a DNA topoisomerase, which controls and changes the topological status of DNA in the process of transcription and functions as a target for some antitumor agents (Delgado et al., 2018; Kitdumrongthum et al., 2020). Other bioinformatic analyses showed that TOP2A was related to development in cancers of the liver, esophagus, stomach, cervix, and lung, among others (Zeng et al., 2019; Kou et al., 2020; Zhang T. et al., 2020; Zhou et al., 2020). CCNB2 is an important element for the process of cell cycle regulation. U2AF homology motif kinase 1 facilitates the nuclear enrichment of MYB proto-oncogene like 2 by affecting the expression of CCNB2 to regulate cell cycle and proliferation (Wei et al., 2019), and reduced transmembrane protein 9 can contribute to decreased CCNB2 levels and then promote cell cycle arrest in HCC cells (Zhang et al., 2016). PRC1, RACGAP1, and CDCA5 were identified as the crucial genes in the pathological progression from cirrhosis to HCC, and their hypomethylation may drive the high expression of these genes (Lin Y. et al., 2019). AURKA can induce the metastasis of irradiated residuary HCC while promoting an epithelial-mesenchymal transition and cancer stem cell properties (Chen et al., 2017). MYC proto-oncogene and AURKA regulate the expression of each other at a mRNA level identified as a MYC-AURKA feedback loop (Lu et al., 2015). CDKN3 overexpression can shorten the survival of HCC cells and shift sensitivity to chemotherapeutic drugs across the AKT/P53/P21 signaling pathway; besides, CDKN3 has been shown to be downregulated in advanced tumor stages (Dai et al., 2016). On the contrary, Chunyang et al. presented that the upregulation of CDKN3 might facilitate cell proliferation via the stimulation of the G1-S transition (Xing et al., 2012). NUSAP1 is a target for mir193A-5p, and evidence has shown that mir193A-5p might block the tumorigenesis of HCC through reducing NUSAP1 (Roy et al., 2018). ASPM is associated with the development of HCV-related cirrhosis via the regulation of tumor-associated phosphorylation (Wang et al., 2017). ASPM is also considered a prognostic biomarker that predicts the increased possibility of invasive or metastatic HCC (Lin et al., 2008). NCAPG down-regulation indicates the suppression of HCC progression, possibly via the PI3K-AKT signaling pathway (Gong et al., 2019; Wang Y. et al., 2019). However, among these candidate hub genes, only PRC1 has been reported to be an immune-related gene in a weighted gene co-expression network analysis (Wang et al., 2020).

Previous bioinformatics analyses have revealed that some TOP2A, CCNB2, PRC1, RACGAP1, AURKA, NUSAP1, CDCA5, ASPM, and NCAPG can be identified as key genes basing on different screening rules (Cai et al., 2019; Wang M. et al., 2019; Zhou Z. et al., 2019; Song et al., 2020). Comparing these previous studies, ours has the following advantages: First, this study included five GEO datasets, while others included two or three gene expression microarrays. In general, a greater number of included samples indicate more credible results in integrated research. Second, we constructed ten gene-macrophage Cox regression models. Inevitably, there were still several limitations of the present study. The included datasets came from different platforms, which might lead to an uncertain systematic bias. Furthermore, TIMER2.0 is a visual website based on tumor tissue information from the Cancer Genome Atlas database (Li et al., 2020). Thus, although tumor purity adjustment was performed in the correlation analyses between the immune cell and candidate genes, there was still systematic bias. To overcome this issue, the application of single-cell RNA sequencing at a higher resolution should be conducted (Papalexi and Satija, 2018). Finally, future experiments in vivo/in vitro should be performed to verify the results of this bioinformatics analysis.

Numerous studies have revealed that immune cell infiltration, TAMs, for instance, can serve as a biomarker for the diagnosis and prognosis of various cancers (Ali et al., 2016; Zhou R. et al., 2019). Thus, we assessed the prognostic value of the combination of TAMs and expression patterns for each of the hub genes. The results showed that the high TAM level predicted unfavorable prognosis under the condition of the high expression of the hub genes, while there was no significant correlation between the TAM and prognosis under the condition of the low expression level of these genes. These results suggested that the combination of hub-gene expression and the TAM levels would play a more effective role in the prognosis prediction of HCC.

In summary, we identified ten genes with a positive correlation with infiltrated-immune cells. These candidate genes present the marked prognostic value of HCC and act as independent prognosis factors for patients with HCC. Moreover, these genes may function in the progression of HCC. Furthermore, we determined that combining the expression of these genes and TAMs can provide a more efficient HCC prognosis prediction. Overall, these findings suggest that these hub genes may be potential targets of immune therapy.

References

• 1

  AbeT.TashiroH.KobayashiT.HattoriM.KurodaS.OhdanH. (2017). Glasgow prognostic score and prognosis after hepatectomy for hepatocellular carcinoma.World J. Surg.411860–1870. 10.1007/s00268-017-3909-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AliH. R.ChlonL.PharoahP. D.MarkowetzF.CaldasC. (2016). Patterns of immune infiltration in breast cancer and their clinical implications: a gene-expression-based retrospective Study.PLoS Med.13:e1002194. 10.1371/journal.pmed.1002194

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BuczakK.OriA.KirkpatrickJ. M.HolzerK.DauchD.RoesslerS.et al (2018). Spatial tissue proteomics quantifies inter- and intratumor heterogeneity in hepatocellular carcinoma (HCC).Mol. Cell. Proteom. MCP17810–825. 10.1074/mcp.RA117.000189

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  CaiC.WangW.TuZ. (2019). Aberrantly DNA methylated-differentially expressed genes and pathways in hepatocellular carcinoma.J. Cancer10355–366. 10.7150/jca.27832

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  CarusoS.NaultJ. C. (2019). A dive into the deep heterogeneity of hepatocellular carcinoma.Gastroenterology1571477–1479. 10.1053/j.gastro.2019.10.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  ChenC.ChenH.ZhangY.ThomasH. R.FrankM. H.HeY.et al (2020). TBtools - an integrative toolkit developed for interactive analyses of big biological data.Mol. Plant[Epub ahead of print]. 10.1016/j.molp.2020.06.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  ChenC.SongG.XiangJ.ZhangH.ZhaoS.ZhanY. (2017). AURKA promotes cancer metastasis by regulating epithelial-mesenchymal transition and cancer stem cell properties in hepatocellular carcinoma.Biochem. Biophys. Res. Commun.486514–520. 10.1016/j.bbrc.2017.03.075

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  DaiW.MiaoH.FangS.FangT.ChenN.LiM. (2016). CDKN3 expression is negatively associated with pathological tumor stage and CDKN3 inhibition promotes cell survival in hepatocellular carcinoma.Mol. Med. Rep.141509–1514. 10.3892/mmr.2016.5410

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  DavisS.MeltzerP. S. (2007). GEOquery: a bridge between the gene expression omnibus (GEO) and BioConductor.Bioinformatics231846–1847. 10.1093/bioinformatics/btm254

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  DelgadoJ. L.HsiehC. M.ChanN. L.HiasaH. (2018). Topoisomerases as anticancer targets.Biochem. J.475373–398. 10.1042/BCJ20160583

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  DengY.ZhuJ.LiuZ.HuangM.ChangD. W.GuJ. (2020). Elevated systemic inflammatory responses, factors associated with physical and mental quality of life, and prognosis of hepatocellular carcinoma.Aging124357–4370. 10.18632/aging.102889

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  DuanM.GoswamiS.ShiJ. Y.WuL. J.WangX. Y.MaJ. Q.et al (2019). Activated and exhausted MAIT cells foster disease progression and indicate poor outcome in hepatocellular carcinoma.Clin. Cancer Res.253304–3316. 10.1158/1078-0432.CCR-18-3040

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  FornerA.ReigM.BruixJ. (2018). Hepatocellular carcinoma.Lancet3911301–1314. 10.1016/S0140-6736(18)30010-2

  • CrossRef
  • Google Scholar

• 14

  FoxR.BerhaneS.TengM.CoxT.TadaT.ToyodaH.et al (2014). Biomarker-based prognosis in hepatocellular carcinoma: validation and extension of the BALAD model.Br. J. Cancer1102090–2098. 10.1038/bjc.2014.130

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  FriemelJ.FrickL.WeberA. (2015). Intratumor heterogeneity in HCC.Aging7350–351. 10.18632/aging.100760

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  GongC.AiJ.FanY.GaoJ.LiuW.FengQ.et al (2019). NCAPG promotes the proliferation of hepatocellular carcinoma through PI3K/AKT signaling.OncoTargets Ther.128537–8552. 10.2147/OTT.S217916

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  HannonG. J.CassoD.BeachD. (1994). KAP: a dual specificity phosphatase that interacts with cyclin-dependent kinases.Proc. Natl. Acad. Sci. U.S.A.911731–1735. 10.1073/pnas.91.5.1731

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  HuangR.ChenZ.LiW.FanC.LiuJ. (2020). Immune systemassociated genes increase malignant progression and can be used to predict clinical outcome in patients with hepatocellular carcinoma.Int. J. Oncol.561199–1211. 10.3892/ijo.2020.4998

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  ImamuraH.MatsuyamaY.TanakaE.OhkuboT.HasegawaK.MiyagawaS.et al (2003). Risk factors contributing to early and late phase intrahepatic recurrence of hepatocellular carcinoma after hepatectomy.J. Hepatol.38200–207. 10.1016/s0168-8278(02)00360-4

  • CrossRef
  • Google Scholar

• 20

  JayantK.HabibN.HuangK. W.PoddaM.WarwickJ.ArasaradnamR. (2020). Immunological basis of genesis of hepatocellular carcinoma: unique challenges and potential opportunities through immunomodulation.Vaccines8:247. 10.3390/vaccines8020247

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  JohnsonL. N.De MolinerE.BrownN. R.SongH.BarfordD.EndicottJ. A.et al (2002). Structural studies with inhibitors of the cell cycle regulatory kinase cyclin-dependent protein kinase 2.Pharmacol. Ther.93113–124. 10.1016/s0163-7258(02)00181-x

  • CrossRef
  • Google Scholar

• 22

  KimS.RovgaliyevB.LeeJ. M.LeeK. W.HongS. K.ChoJ. H.et al (2020). Clinical significance of de novo malignancy after liver transplant: a single-center study.Transpl. Proc.[Epub ahead of print]. 10.1016/j.transproceed.2020.02.148

  • CrossRef
  • Google Scholar

• 23

  KitdumrongthumS.ReabroiS.SuksenK.TuchindaP.MunyooB.MahalapbutrP.et al (2020). Inhibition of topoisomerase IIalpha and induction of DNA damage in cholangiocarcinoma cells by altholactone and its halogenated benzoate derivatives.Biomed. Pharmacother. Biomed. Pharmacother.127:110149. 10.1016/j.biopha.2020.110149

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  KouF.SunH.WuL.LiB.ZhangB.WangX.et al (2020). TOP2A promotes lung adenocarcinoma cells’ malignant progression and predicts poor prognosis in lung adenocarcinoma.J. Cancer112496–2508. 10.7150/jca.41415

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  LanghansB.NischalkeH. D.KramerB.DoldL.LutzP.MohrR.et al (2019). Role of regulatory T cells and checkpoint inhibition in hepatocellular carcinoma.Cancer Immunol. Immunother. CII682055–2066. 10.1007/s00262-019-02427-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  LiT.FuJ.ZengZ.CohenD.LiJ.ChenQ.et al (2020). TIMER2.0 for analysis of tumor-infiltrating immune cells.Nucleic Acids Res.48W509–W514. 10.1093/nar/gkaa407

  • CrossRef
  • Google Scholar

• 27

  LiZ.WuT.ZhengB.ChenL. (2019). Individualized precision treatment: targeting TAM in HCC.Cancer Lett.45886–91. 10.1016/j.canlet.2019.05.019

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  LinJ.HouY.HuangS.WangZ.SunC.WangZ.et al (2019). Exportin-T promotes tumor proliferation and invasion in hepatocellular carcinoma.Mol. Carcinogenesis58293–304. 10.1002/mc.22928

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  LinS. Y.PanH. W.LiuS. H.JengY. M.HuF. C.PengS. Y.et al (2008). ASPM is a novel marker for vascular invasion, early recurrence, and poor prognosis of hepatocellular carcinoma.Clin. Cancer Res.144814–4820. 10.1158/1078-0432.CCR-07-5262

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  LinY.LiangR.YeJ.LiQ.LiuZ.GaoX.et al (2019). A twenty gene-based gene set variation score reflects the pathological progression from cirrhosis to hepatocellular carcinoma.Aging1111157–11169. 10.18632/aging.102518

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  LuL.HanH.TianY.LiW.ZhangJ.FengM.et al (2015). Aurora kinase A mediates c-Myc’s oncogenic effects in hepatocellular carcinoma.Mol. Carcinogenesis541467–1479. 10.1002/mc.22223

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  MaX.ZhouL.ZhengS. (2020). Transcriptome analysis revealed key prognostic genes and microRNAs in hepatocellular carcinoma.PeerJ8:e8930. 10.7717/peerj.8930

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  OkamotoK.KitabayashiI.TayaY. (2006). KAP1 dictates p53 response induced by chemotherapeutic agents via Mdm2 interaction.Biochem. Biophys. Res. Commun.351216–222. 10.1016/j.bbrc.2006.10.022

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  PapalexiE.SatijaR. (2018). Single-cell RNA sequencing to explore immune cell heterogeneity.Nat. Rev. Immunol.1835–45. 10.1038/nri.2017.76

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  ParkH.ParkJ. Y. (2013). Clinical significance of AFP and PIVKA-II responses for monitoring treatment outcomes and predicting prognosis in patients with hepatocellular carcinoma.BioMed Res. Int.2013:310427. 10.1155/2013/310427

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  PontenF.JirstromK.UhlenM. (2008). The Human Protein Atlas–a tool for pathology.J. Pathol.216387–393. 10.1002/path.2440

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  RhodesD. R.Kalyana-SundaramS.MahavisnoV.VaramballyR.YuJ.BriggsB. B.et al (2007). Oncomine 3.0: genes, pathways, and networks in a collection of 18,000 cancer gene expression profiles.Neoplasia9166–180. 10.1593/neo.07112

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  RoyS.HooiveldG. J.SeehawerM.CarusoS.HeinzmannF.SchneiderA. T.et al (2018). microRNA 193a-5p regulates levels of nucleolar- and spindle-associated protein 1 to suppress hepatocarcinogenesis.Gastroenterology1551951.e26–1966.e26. 10.1053/j.gastro.2018.08.032

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  SiegelR. L.MillerK. D.JemalA. (2015). Cancer statistics, 2015.CA Cancer J. Clin.655–29. 10.3322/caac.21254

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  SongX.DuR.GuiH.ZhouM.ZhongW.MaoC.et al (2020). Identification of potential hub genes related to the progression and prognosis of hepatocellular carcinoma through integrated bioinformatics analysis.Oncol. Rep.43133–146. 10.3892/or.2019.7400

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.Nucleic Acids Res.47D607–D613. 10.1093/nar/gky1131

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  TangY.ZengZ.WangJ.LiG.HuangC.DongX.et al (2020). Combined signature of nine immune-related genes: a novel risk score for predicting prognosis in hepatocellular carcinoma.Am. J. Transl. Res.121184–1202.

  • Google Scholar

• 43

  TangZ.KangB.LiC.ChenT.ZhangZ. (2019). GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis.Nucleic Acids Res.47W556–W560. 10.1093/nar/gkz430

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  VillanuevaA. (2019). Hepatocellular carcinoma.New Engl. J. Med.3801450–1462. 10.1056/NEJMra1713263

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  WangD.LiuJ.LiuS.LiW. (2020). Identification of crucial genes associated with immune cell infiltration in hepatocellular carcinoma by weighted gene co-expression network analysis.Front. Genet.11:342. 10.3389/fgene.2020.00342

  • CrossRef
  • Google Scholar

• 46

  WangF.ChangY.LiJ.WangH.ZhouR.QiJ.et al (2017). Strong correlation between ASPM gene expression and HCV cirrhosis progression identified by co-expression analysis.Digestive Liver Dis.4970–76. 10.1016/j.dld.2016.10.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  WangJ.MengJ.RanW.LeeR. J.TengL.ZhangP.et al (2019). Hepatocellular carcinoma growth retardation and PD-1 blockade therapy potentiation with synthetic high-density lipoprotein.Nano Lett.195266–5276. 10.1021/acs.nanolett.9b01717

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  WangM.WangL.WuS.ZhouD.WangX. (2019). Identification of key genes and prognostic value analysis in hepatocellular carcinoma by integrated bioinformatics analysis.Int. J. Genomics2019:3518378. 10.1155/2019/3518378

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  WangY.GaoB.TanP. Y.HandokoY. A.SekarK.DeivasigamaniA.et al (2019). Genome-wide CRISPR knockout screens identify NCAPG as an essential oncogene for hepatocellular carcinoma tumor growth.FASEB J.338759–8770. 10.1096/fj.201802213RR

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  WeiT.WeilerS. M. E.TothM.StichtC.LutzT.ThomannS.et al (2019). YAP-dependent induction of UHMK1 supports nuclear enrichment of the oncogene MYBL2 and proliferation in liver cancer cells.Oncogene385541–5550. 10.1038/s41388-019-0801-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  XieW.WangB.WangX.HouD.SuH.HuangH. (2020). Nine hub genes related to the prognosis of HBV-positive hepatocellular carcinoma identified by protein interaction analysis.Ann. Transl. Med.8:478. 10.21037/atm.2020.03.94

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  XingC.XieH.ZhouL.ZhouW.ZhangW.DingS.et al (2012). Cyclin-dependent kinase inhibitor 3 is overexpressed in hepatocellular carcinoma and promotes tumor cell proliferation.Biochem. Biophys. Res. Commun.42029–35. 10.1016/j.bbrc.2012.02.107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  YinZ.DongC.JiangK.XuZ.LiR.GuoK.et al (2019). Heterogeneity of cancer-associated fibroblasts and roles in the progression, prognosis, and therapy of hepatocellular carcinoma.J. Hematol. Oncol.12:101. 10.1186/s13045-019-0782-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  ZengL.FanX.WangX.DengH.ZhangK.ZhangX.et al (2019). Bioinformatics analysis based on multiple databases identifies hub genes associated with hepatocellular carcinoma.Curr. Genomics20349–361. 10.2174/1389202920666191011092410

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  ZhangJ.ChangL.ZhangX.ZhouZ.GaoY. (2019). Meta-analysis of the prognostic and clinical value of tumor-associated macrophages in hepatocellular carcinoma.J. Investig. Surg.[Epub ahead of print]. 10.1080/08941939.2019.1631411

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  ZhangT.WangB. F.WangX. Y.XiangL.ZhengP.LiH. Y.et al (2020). Key genes associated with prognosis and tumor infiltrating immune cells in gastric cancer patients identified by cross-database analysis.Cancer Biother. Radiopharm.35696–710. 10.1089/cbr.2019.3423

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  ZhangY.RanY.XiongY.ZhongZ. B.WangZ. H.FanX. L.et al (2016). Effects of TMEM9 gene on cell progression in hepatocellular carcinoma by RNA interference.Oncol. Rep.36299–305. 10.3892/or.2016.4821

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  ZhangY.ZhangL.XuY.WuX.ZhouY.MoJ. (2020). Immune-related long noncoding RNA signature for predicting survival and immune checkpoint blockade in hepatocellular carcinoma.J. Cell. Physiol.2359304–9316. 10.1002/jcp.29730

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  ZhaoX.WangX.YouY.WenD.FengZ.ZhouY.et al (2020). Nogo-B fosters HCC progression by enhancing Yap/Taz-mediated tumor-associated macrophages M2 polarization.Exp. Cell Res.391:111979. 10.1016/j.yexcr.2020.111979

  • CrossRef
  • Google Scholar

• 60

  ZhouR.ZhangJ.ZengD.SunH.RongX.ShiM.et al (2019). Immune cell infiltration as a biomarker for the diagnosis and prognosis of stage I-III colon cancer.Cancer Immunol. Immunother. CII68433–442. 10.1007/s00262-018-2289-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  ZhouW.WuJ.LiuX.NiM.MengZ.LiuS.et al (2020). Identification of crucial genes correlated with esophageal cancer by integrated high-throughput data analysis.Medicine99:e20340. 10.1097/MD.0000000000020340

  • CrossRef
  • Google Scholar

• 62

  ZhouY.ZhouB.PacheL.ChangM.KhodabakhshiA. H.TanaseichukO.et al (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets.Nat. Commun.10:1523. 10.1038/s41467-019-09234-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  ZhouZ.LiY.HaoH.WangY.ZhouZ.WangZ.et al (2019). Screening hub genes as prognostic biomarkers of hepatocellular carcinoma by bioinformatics analysis.Cell Transpl.2876S–86S. 10.1177/0963689719893950

  • Pubmed Abstract
  • CrossRef
  • Google Scholar