Journal of Oncology

Journal of Oncology / 2019 / Article

Research Article | Open Access

Volume 2019 |Article ID 2531932 | 15 pages |

Noninvasive Identification of Immune-Related Biomarkers in Hepatocellular Carcinoma

Academic Editor: Thomas R. Chauncey
Received31 Jan 2019
Accepted16 Jul 2019
Published18 Aug 2019


Primary liver carcinoma is one of the most common malignant tumors with a poor prognosis. This study aims to uncover the potential mechanisms and identify core biomarkers of hepatocellular carcinoma (HCC). The HCC gene expression profile GSE49515 was chosen to analyze the differentially expressed genes from purified RNA of peripheral blood mononuclear cells, including 10 HCC patients and 10 normal individuals. GO and KEGG pathway analysis and PPI network were used, and the enrichment of the core genes out of 15 hub genes was evaluated using GEPIA and GSEA, respectively. We employed flow cytometry to count mononuclear cells to verify our opinions. In this analysis, 344 DEGs were captured, including 188 upregulated genes and 156 downregulated genes; besides that, 15 hub genes were identified. GO analysis and KEGG analysis showed the DEGs were particularly involved in immune response, antigen processing and presentation, and infection and inflammation. The PPI network uncovered 2 modules were also mainly involved in immune response. In conclusion, our analysis disclosed the immune subversion was the major signature of HCC associated closely with JUN, VEGFA, TNFSF10, and TLR4, which could be novel noninvasive biomarkers in peripheral blood and targets for early diagnosis and therapy of HCC.

1. Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignancies, especially in the aged, which accounts for approximately 90% of all primary liver cancers severely threatening public health [1]. The mechanism of HCC is a complex process associated with the incremental accumulation of gene mutation, giving rise to abnormal immune subversion, cell cycle, and angiogenesis [24]. As for immune subversion, effector immune cells could execute immune control of HCC, which efficiently decrease malignant transformed cells. However, progression of HCC clearly certifies failure of tumor immune control suggesting inhibition of anticancer immune responses [5]. Especially, tumor-related mononuclear cells collaborate within an inflammatory network, which result in the immune privilege in the tumor environment [6]. Therefore, immunosuppressive mononuclear cells are equivalent to heterogeneous cell lines, including lymphocytes and monocytes cooperating by direct cell contact, secretion of cytokines, or production of extracellular matrix, which lead to the suppression of the immune response in the tumor milieu [7].

Currently, imageological examination and pathological biopsy are the conventional diagnostic methods of HCC [8]. However, imaging displays poor specificity, and pathological biopsy is an invasive method which may result in iatrogenic injury [9]. Therefore, serum biomarkers are routinely used for tumor diagnostic. For example, alpha-fetoprotein (AFP) has been widely used in clinical practice [10]. Although many studies have reported the accuracy of AFP for HCC, solely AFP still has some false-positive or false-negative rate [11]. Hence, the identification of specific and sensitive biomarkers is necessary in order to achieve accurate diagnosis and treatment of HCC as early as possible, especially noninvasive biomarkers.

High-throughput gene microarray is increasingly being widely used, which can analyze cancer and noncancer samples indicating us tumor-related genes at multiple levels from molecular diagnosis and pathological classification to therapeutic evaluation and prognosis prediction, as well as drug sensitivity and neoplasm recurrence [1214]. However, the use of microarrys in clinical application is restricted by countless number of genes identified by gene profiling, lack of both repeatability and independent verification, and requirement for complex statistical analyses. Moreover, most of the microarrys are based on the genes in tissues which are difficult to detect except by invasive methods [15]. Therefore, in order to put these expression profiles into clinical applications as soon as possible, it is necessary to identify an appropriate amount of serum genes and develop a suitable way that can be done by routine assay.

In this study, we downloaded the HCC gene expression profile GSE49515 in the Gene Expression Omnibus (GEO,, an online public collection database for microarray data and used GEO2R online software to compare gene expression profiles of tumor cells with normal liver cells to identify differentially expressed genes (DEGs). Then, we constructed the protein-protein interaction (PPI) network of the DEGs and selected 15 hub genes according to a high degree of connectivity. Following this, we analyzed gene ontology (GO) and pathway enrichment including the biological process (BP), molecular function (MF), cellular component (CC), and KEGG pathway of the DEGs. Moreover, we performed two modules and confirmed their enriched pathways. The core genes out of the 15 hub ones were found, and the interactions between any of them were detected with the help of GEPIA. After the analysis of the core status and biological function of any hub gene, we performed flow cytometry (FCM) to count mononuclear cells, confirming the findings.

2. Results

2.1. Identification of DEGs and Hub Genes

A comparison of 10 HCC samples with 10 normal samples in our study was performed by employing the GEO2R online analysis tool based on value < 0.05 and logFC ≤ −2 or logFC ≥ 2 criteria. A total of 344 DEGs were picked up after analyzing GSE49515, 188 of which were upregulated while 156 were downregulated (Figure 1(b)). The expression levels of the top 50 DEGs were displayed in a heat map to visualize the results (Figure 1(a)).

2.2. GO Function- and KEGG Pathway-Enrichment Analysis

To gain a more extensive and in-depth knowledge of those selected DEGs, we use DAVID to analyze significantly enriched GO function and KEGG pathways. After inputting all of the DEGs to DAVID online analysis tool, we obtained the GO analysis of these upregulated DEGs and downregulated DEGs. The results showed that these DEGs were mainly enriched in biological processes (BP), including apoptotic process, immune response, and inflammatory response, among which were positive regulation of NF-κB transcription factor activity and cell-cell signaling for downregulation, positive regulation of angiogenesis, negative regulation of cell proliferation, positive regulation of cell proliferation, mitotic spindle organization, and neutrophil chemotaxis for upregulation. Concerning molecular function (MF), the downregulated DEGs were particularly related to receptor binding, iron ion binding, haptoglobin binding, oxygen transporter activity, and peroxidase activity, while the upregulated DEGs were mainly implicated with nucleotide binding and ubiquitin protein ligase binding. Besides, GO cell component (CC) analysis indicated that the downregulated DEGs were mainly enriched in cytosol, extracellular exosome, Golgi membrane, blood microparticle, and nuclear chromosome (telomeric region) and the upregulated DEGs were principally enriched in nucleus, nucleoplasm, platelet alpha granule and extracellular space (Table 1).

CategoryTermCount% value

UpregulatedGOTERM_BP_DIRECTGO:0045766∼positive regulation of angiogenesis5<0.001<0.001
GOTERM_BP_DIRECTGO:0008285∼negative regulation of cell proliferation50.010.01
GOTERM_BP_DIRECTGO:0008284∼positive regulation of cell proliferation50.020.02
GOTERM_BP_DIRECTGO:0007052∼mitotic spindle organization3<0.001<0.001
GOTERM_BP_DIRECTGO:0030593∼neutrophil chemotaxis30.020.02
GOTERM_CC_DIRECTGO:0031091∼platelet alpha granule20.010.04
GOTERM_CC_DIRECTGO:0005615∼extracellular space90.060.05
GOTERM_MF_DIRECTGO:0000166∼nucleotide binding50.030.04
GOTERM_MF_DIRECTGO:0031625∼ubiquitin protein ligase binding30.020.05

DownregulatedGOTERM_BP_DIRECTGO:0006915∼apoptotic process90.060.02
GOTERM_BP_DIRECTGO:0006955∼immune response80.050.02
GOTERM_BP_DIRECTGO:0006954∼inflammatory response70.040.03
GOTERM_BP_DIRECTGO:0007267∼cell-cell signaling60.040.05
GOTERM_BP_DIRECTGO:0051092∼positive regulation of NF-kappaB transcription factor activity60.04<0.001
GOTERM_CC_DIRECTGO:0070062∼extracellular exosome260.160.04
GOTERM_CC_DIRECTGO:0000139∼golgi membrane90.060.03
GOTERM_CC_DIRECTGO:0072562∼blood microparticle60.04<0.001
GOTERM_CC_DIRECTGO:0000784∼nuclear chromosome, telomeric region40.030.05
GOTERM_MF_DIRECTGO:0005102∼receptor binding70.040.03
GOTERM_MF_DIRECTGO:0005506∼iron ion binding50.030.02
GOTERM_MF_DIRECTGO:0031720∼haptoglobin binding30.02<0.001
GOTERM_MF_DIRECTGO:0005344∼oxygen transporter activity30.02<0.001
GOTERM_MF_DIRECTGO:0004601∼peroxidase activity30.020.01

Afterwards, we analyzed the most significantly enriched KEGG pathway of the upregulated and downregulated DEGs, which is shown in Table 2. The downregulated DEGs were involved in measles, influenza A, rheumatoid arthritis, antigen processing and presentation, and legionellosis, while the upregulated DEGs were involved in bladder cancer, rheumatoid arthritis, malaria, herpes simplex infection, and osteoclast differentiation. The scatter plots in Supplement 1 (A, B, and C) show a GO and KEGG pathway-enrichment plot of HCC.

CategoryTermCount% valueGenes

Up regulatedssc05219: bladder cancer40.020.002VEGFA, CXCL8, HBEGF, THBS1
ssc05323: rheumatoid arthritis50.030.002JUN, IFNG, VEGFA, CXCL8, ITGB2
ssc05144: malaria40.020.004IFNG, CXCL8, ITGB2, THBS1
ssc05168: herpes simplex infection50.030.024SRSF3, JUN, IFNG, CYCS, CUL1
ssc04380: osteoclast differentiation40.020.047FOSL2, SQSTM1, JUN, IFNG

Down regulatedhsa05162: measles80.05<0.001CCNE2, TNFSF10, EIF2S1, TLR4, HSPA1A, HSPA1B, MSN, TLR7
hsa05164: influenza A80.050.001TNFSF10, EIF2S1, TLR4, HSPA1A, HSPA1B, TLR7, HLA-DQA1, HLA-DRA
hsa05323: rheumatoid arthritis60.040.001CXCL5, TLR4, LTB, HLA-DQA1, HLA-DRA, IL11
hsa04612: antigen processing and presentation50.030.004HSPA1A, HSPA1B, KLRD1, HLA-DQA1, HLA-DRA
hsa05134: legionellosis40.030.012NLRC4, TLR4, HSPA1A, HSPA1B

2.3. Hub Genes and Module Screening from PPI Network

Besides, 15 hub genes from Cytoscape software were identified in accordance with a high degree of connectivity selected (Table 3). We built the PPI network of the top 15 hub genes via the information of the STRING protein query from public databases (Figure 2(a)). The top 15 hub genes with a higher degree of connectivity are as follows: JUN, IL8, VEGFA, TLR4, IFNG, TNFSF10, EHHADH, ATF3, FUS, DUSP1, HSPA1A, CUL1, FPR2, POLR2H, and RHOB. Based on the GO function and KEGG pathway analysis of the profile, we uncovered JUN, VEGFA, TNFSF10, and TLR4 that were enriched in immune response-related pathway.

NameDegree valueLog FC


Moreover, we used MCODE plug-in to detect the highest modules in the PPI network. We chose the top 2 modules, and GO function- and KEGG pathway-enrichment analysis indicated that Module 1 was related to immune response, apoptotic process, and epithelial cell migration, while Module 2 was associated with a signaling pathway and cellular response to various substances (Figures 2(b) and 2(c), Table 4).

ModuleTerm valueFDRGenes

Module 1Immune response0.0222.28TNFSF10, JUN, CCR1
Positive regulation of cysteine-type endopeptidase activity involved in apoptotic process0.0224.07TNFSF10, CYCS
Cytokine-cytokine receptor interaction0.0114.95TNFSF10, CCR1, IFNG, VEGFA
Positive regulation of epithelial cell migration0.0221.94JUN, IFNG
HIF-1 signaling pathway0.0222.98IFNG, VEGFA, TLR4

Module 2Cellular response to corticotropin-releasing hormone stimulus<0.0011.74NR4A2, NR4A1
Intracellular receptor signaling pathway<0.0016.26NR4A2, NR4A1
Steroid hormone-mediated signaling pathway0.0123.88NR4A2, NR4A1
MAPK signaling pathway<0.0010.90DUSP1, NR4A1, HSPA1A
Negative regulation of cysteine-type endopeptidase activity involved in apoptotic process0.0222.82NR4A1, HSPA1A

2.4. The Kaplan–Meier Plotter and Expression Level of Hub Genes

The prognosis of the top 15 hub genes was analyzed in We found that the expression of JUN () was related to worse overall survival (OS) for HCC patients, as well as IL8 (called as CXCL8) (), VEGFA (), EHHADH (), FUS (), HSPA1A (), CUL1 (), and POLR2H (). The level of TLR4 (), IFNG (), TNFSF10 (), ATF3 (), DUSP1 (), FPR2 (), and RHOB () had no obvious difference in overall survival of HCC patients (Figure 3). However, the survival curves are analyzed with liver tissue, which can only indirectly explain the importance of hub genes in PBMC. These hub genes in PBMC which can be the biomarkers for early diagnosis may not be easy to detect in liver tissue.

Then, we employed DAVID to analyze the correlation of 15 hub genes. We found JUN, IFNG, VEGFA, TLR4, and TNFSF10 are the 5 high-degree-of-connectivity genes which are fully associated with immune response, inflammatory response, and HIF-1 signaling pathway (Table 5). Then, in order to confirm the most relevant hub genes, we used correlation analysis in GEPIA, and we detected JUN and VEGFA, JUN and ATF3, and JUN and RHOB are distinctly correlated ( value = 0, ; value = 0, ; value = 0, ), which means JUN may be the core gene of HCC (Figures 4(e)4(g)).

TermCount% valueGenesFDR

ssc05323: rheumatoid arthritis40.181.41E-04JUN, IFNG, VEGFA, TLR40.14
ssc05164: influenza A40.189.95E-04JUN, IFNG, TNFSF10, TLR41.04
ssc05321: inflammatory bowel disease (IBD)30.130.002JUN, IFNG, TLR42.66
ssc04066: HIF-1 signaling pathway30.130.007IFNG, VEGFA, TLR47.01
ssc04060: cytokine-cytokine receptor interaction30.130.031IFNG, VEGFA, TNFSF1027.91

2.5. Gene Set Enrichment Analysis

In order to make the further function of the hub gene clear, GSEA was used to map into GO analysis and KEGG pathways database. Under the cutoff criteria nominal value < 0.05, │enrichment score (ES)│ > 0.6, and gene size ≥ 100, six functional gene sets were enriched in total, which were particularly centralized on pathway associated with immune response and inflammatory response. Six pathways were “inflammatory response to antigenic stimulus,” “regulation of T cell migration,” “antigen processing and presentation via MHC class IB,” “negative regulation of B cell activation,” “negative regulation of NF-kB signaling,” and “positive regulation of interleukin-1 production” (Figure 5).

2.6. Identification of Biomarkers

To confirm the results we have stated above, we test the mRNA level of four key genes we predicted (JUN, TLR4, VEGFA, and TNFSF10). We found the mRNA level of JUN, TLR4, VEGFA, and TNFSF10 in the PBMC of HCC patients was significantly upregulated in the PBMC of HCC patients (), of which VEGFA increased obviously () (Figures 4(a)4(d)).

2.7. Immune Subversion in HCC

According to our prediction, the negative regulation of immune cells in HCC patients’ peripheral blood significantly occurred detected by FCM (Figures 6(a)6(h)). In detail, compared with the healthy subjects, the levels of T lymphocytes in peripheral blood, both helper T cells and cytotxic T cells, were significantly lower in HCC patients, which is mainly the immune mechanism of tumor patients () (Figure 6(k)). Meanwhile, the B lymphocytes as well as NK cells also decreased in HCC patients, especially NK cells () (Figures 6(i) and 6(j)). Taken together, our experiment of FCM demonstrated that T cell migration and B cell activation in adaptive immune and NK cells inhibition in innate immune were important mechanisms and could do further research in future.

3. Discussion

In recent decades, the morbidity and mortality of HCC have been increasing worldwide. Although the early diagnosis and treatment have developed a lot recently, the overall survival rates of HCC is still poor [16]. Therefore, the sensitive and specific biomarkers for HCC are urgently necessary. High-throughput studies can develop the thorough exploration of the vital mechanisms which lead to HCC. In our study, we identified DEGs between 10 HCC samples and 10 normal samples from the GEO database of GSE49515. In order to increase the statistical power of these DEGs, we defined that the absolute value of the logarithm (base 2) fold change (logFC) greater than 2 and a total of 344 DEGs were captured, including 188 upregulated genes and 156 downregulated genes. In order to have an in-depth detection of these DEGs, we employed GO function, KEGG pathway, PPI network, and connectivity analysis of these DEGs, via which we found that HCC-related genes and pathways have great importance in cancer initiation and progression.

There were plenty of mechanisms uncovered to contribute to the development of HCC, but the predominant mechanism implicated with tumorigenesis is still controversial, thereby causing difficulties to the diagnosis and treatment of HCC. GO term enrichment and PPI analysis in our study disclosed that downregulated DEGs were mainly associated with immune response. Our experimental results detecting the mononuclear cells in peripheral blood mononuclear cells of HCC patients and healthy individuals further confirmed the important role of immune subversion.

Immune response is well recognized to play a vital part in the initiation and progression of carcinogenesis because the development of HCC obviously records failure of tumor immune control which stands for immune subversion by the tumor environment [2]. To our knowledge, protective immune surveillance of tumor is mainly conducted by tumor-directed NK cells and lymphocytes, which can effectively identify and eradicate malignant cells [6]. On the one hand, for the reason that innate immune cells are capable of eliminating malignant cells and pertaining to the first-line defense to restrain tumor initiation and progression [17], they act as an essential player in the immunological surveillance. On the other hand, adaptive immune cells including B lymphocytes and T lymphocytes develop along with innate immune responses, which finally target tumor-associated antigens (TAAs) [18]. To be more specific, once the immune cells in peripheral blood are unable to generate or be recruited, specific immunity and nonspecific immunity function will degrade, eventually causing the rapid growth of carcinoma tissue. The quantity of immune cells is associated with many immune-related genes which may be differentially expressed during tumor manifestation and progression.

NK cells (CD3−/CD56+) reside in the liver and account for about 30% to 50% of the hepatic lymphocytes in humans [17]. Therefore, NK cells constitute a corresponding effector cell population of innate immune cells contributing to tumor surveillance within the liver. In our data from FCM, we identified that NK cells in HCC patients’ peripheral blood were significantly lower than those in the blood from normal people. Additionally, our GO analysis indicated that the downregulated DEGs were mainly related to immune response including immune cells manifestation and recruitment and release of cytokines. We found that there were 4 genes (JUN, VEGFA, TNFSF10, and TLR4) closely related to the downregulation of immune response which occurs in HCC. Mgrditchian et al. [19] found that the phosphorylation of JUN could induce NK cell infiltration into the tumor bed by inducing the transcription of CCL5, which eventually results in targeted autophagy. However, NK cells produce vascular endothelial growth factor A (VEGFA) in tumor tissues, which may enhance the formation of tumor by the way of angiogenesis. The presence of VEGFA-secreting NK cells is associated with a poor prognosis and depends on the type and stage of the tumor [20]. In our GO analysis, we detected that the upregulated DEGs were particularly associated with the positive regulation of angiogenesis, which indirectly verified that VEGFA-secreting NK cells might play an important part in the HCC. Wagner et al. reported that membrane-bound tumor necrosis factor ligand superfamily member 10 (TNFSF10) on the NK cells can supplement the perforin/granzyme pathway in a NK cell-mediated cytotoxicity, which can enhance the NK cell function of tumor elimination [21]. Besides, there were a few reports which showed that some antitumor substances can induce NK cells to proliferate and release IFN-γ via TLR4; thus, we predicted that TLR4 might be an important target to influence HCC development [2224]. Our result suggested that innate immune cells, especially the function of NK cells, are downregulated via some hub genes in HCC.

Adaptive immunity including humoral immunity and cellular immunity is highly specific elimination of transformed cells, which protects from tumor manifestation and progression. B lymphocytes, CD4+ T lymphocytes, and CD8+ T lymphocytes are vital members in adaptive immunity [25]. Equally, we employed FCM to detect the B lymphocytes and two T lymphocyte subsets which verified the reduction of adaptive immune function in HCC. We predicted the mechanisms of B lymphocytes and T lymphocyte declining may also be associated with the four hub genes (JUN, VEGFA, TNFSF10, and TLR4) via GO and KEGG analysis. c-JUN NH2-terminal kinase (JNK) signaling pathway was implicated in various T cell functions. The JNKs are synergistically activated by stimulation of the TCR with antibodies to its CD3 component and the CD28 auxiliary receptor, which was related to T cell activation [26]. Patterson et al. [27] discovered that Igα non-ITAM tyrosine 204 promoted T-independent B cell proliferation and differentiation via phosphorylation of JUN. In addition, tumor can produce VEGFA, which increases expression of inhibitory immune checkpoints mediating T cell exhaustion on intratumoral CD8+ T cells [28]. The previous report identified that VEGFA directly triggers regulatory T cells (Treg) proliferation resulting in tumors escape, which can become a therapeutic goal in the future [29]. TNFSF10 influences T cells function, which cannot only enhance the maximal suppressive function of Treg cells leading to the antitumor effect but also be involved in T cell-mediated killing of cancer cells [30, 31]. Furthermore, Fang et al. demonstrated that TLR4 appears to induce antitumor T cell response by activating dendritic cells [32]. Similarly, there is research confirming that TLR4 can promote tumor-specific cytotoxic T cell responses [33]. Taken together, from our perspective based on the GO and KEGG analyses and FCM, the regulation of immune-related genes including JUN, VEGFA, TNFSF10, and TLR4 can influence the immune response and be involved in the occurrence and progression of HCC. Drugs targeting immune response of peripheral mononuclear cells may be the potential candidates for therapy of HCC.

The PPI network could form a visible framework for a better understanding of the function of the proteome [34]. From the enriched pathways of top 2 modules, we discovered that the interactions among the proteins in HCC are particularly associated with pathways of immune response and cytokine-cytokine receptor interaction. It emphasizes that immune-related gene interaction can regulate the tumor-associated immune surveillance. Therefore, we utilized DAVID to uncover the correlation of 15 hub genes. Coincidentally, the 4 genes, JUN, VEGFA, TNFSF10, and TLR4, have the high degree of connectivity, especially JUN and VEGFA. We predicted that JUN can regulate mononuclear cells to release VEGFA, which may promote tumor angiogenesis, which was proved as the reason for tumor initiation and progression. Furthermore, we analyzed the gene set enrichment to make sure of the primary function of the 4 hub genes. We found these 4 hub genes were inextricably linked with the various processes of immune and inflammatory responses like the regulation of T and B cell migration and antigen processing via MHC class I B. Hence, monitoring the immune process of tumor immunity including immune-related genes and cytokines is of great importance for the diagnosis and treatment of HCC.

Moreover, HCC patients always used to suffer from hepatitis infected by hepatitis B virus (HBV) and hepatitis C virus (HCV). Patients who suffer cirrhosis or HCC from chronic virus hepatitis account for 90% [35]. It is consistent with our GO and KEGG analyses that the inflammation is an important process of HCC. JUN, VEGFA, TNFSF10, and TLR4 are connectivity genes involved in the inflammatory response. For instance, there was a report that stated JUN can act as an intermediate in antiviral immunity and TLR4 recognizes bacterial ligands to constrain the ability of antiviral immunity which combine bacterial and viral infections [36, 37]. As a result, we hypothesized that JUN, VEGFA, TNFSF10, and TLR4 in PBMC can also regulate virus infection of HCC patients. Detecting these genes can early protect patients from HCC occurrence and development as well as treat them early.

In conclusion, we provide a comprehensive and novel analysis of gene expression profiles to recognize DEGs which may play a core part in the development and prognosis in patients with HCC. Genes implicated with immune and inflammatory response were significantly altered in HCC patients. In order to acquire more precise correlation results, we plan to carry out subsequent authentication experiment later to prove these predictive results. Taken together, we found that JUN, VEGFA, TNFSF10, and TLR4 in PBMC play core roles in the immune response of HCC. Hence, these 4 genes may serve as potential serum biomarkers combining with AFP and targets of immunotherapy. We expect this analysis method will offer accurate and valuable information for future study on the molecular mechanisms of HCC and supply evidences for the detection of new diagnosis biomarkers and therapeutic strategies.

4. Materials and Methods

4.1. Microarray Data

We downloaded the gene expression profile of GSE49515 from GEO database. We chose a total of 20 samples, involving 10 cases of HCC and 10 healthy cases of purified RNA of peripheral blood mononuclear cells (PBMCs) in GSE49515, which was based on Agilent Gpl570 platform ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) by Shi et al. We also get the Series Matrix File of GSE49515 from GEO database.

4.2. Screen Genes of Differential Expression

The analysis of DEGs between HCC samples and normal liver tissue samples was performed by using GEO2R (, an online analysis method for GEO database on the basis of R language. We regulated DEGs as differentially expressed with logFC < −2 (upregulated genes) or logFC > 2 (downregulated genes), according to the criteria described in [38, 39]. The adjusted value < 0.05 was treated as statistically significant in order to reduce the false-positive rate. Thereafter, 344 DEGs were picked up, containing 188 upregulated genes and 156 downregulated genes. We utilized visual hierarchical cluster analysis to get the heatmap and volcano plot of two groups by ImageGP ( after the correlative raw data of TXT files were downloaded.

4.3. Gene Ontology and KEGG Pathway Analysis of DEGs

Gene ontology (GO) analysis is a common framework which can annotate genes and gene products including functions of cellular components, biological pathways, and molecular function [40]. Kyoto Encyclopedia of Genes and Genomes (KEGG) contains a set of genomes and biological pathways related to disease and drugs online database, which essentially is a resource for systematic understanding of biological system and certain high-level genome functional information [41]. The Database for Annotation, Visualization and Integrated Discovery (DAVID, is an online bioinformatics database. It has covered many biological data and relevant analysis tools, and then provided tools for the biological function annotation information for plenty of genes or proteins [42]. was considered as the cutoff criterion with significant difference. We could visualize the key biological processes, molecular functions, and cellular components and pathways of DEGs by using DAVID online database. And the scatter plot was performed by ImageGP according to the results of GO and KEGG pathway.

4.4. PPI Network and Module Analysis

Search Tool for the Retrieval of Interacting Genes (STRING) is an online tool for the assessment and integration of the protein-protein interaction (PPI) information, containing physical and functional associations. It covered 9,643,763 proteins from 2031 organisms in STRING version 10.0 [43]. We drew DEGs using STRING to evaluate the interactional associations among them, thereby utilizing the Cytoscape software to build a PPI network. In the meantime, we set maximum number of interactors = 0, confidence score ≥0.4 as the cutoff criterion. Moreover, the Molecular Complex Detection (MCODE) app was utilized to screen modules of the PPI network in Cytoscape in line with degree cutoff = 2, k-core = 2, node score cutoff = 0.2, and max. depth = 100. And we chose the top 15 genes with a high degree of connectivity as hub genes. The pathway analysis of genes in every module was worked out according to DAVID. Then, 15 hub genes were also mapped into STRING with maximum number of interactors ≤5 and confidence score ≥0.4. GO and KEGG pathway analyses were also used to explain the potential information.

4.5. Comparison of the Hub Genes Expression Level

GEPIA ( is a newly developed interactive web server which analyzes the RNA-sequencing expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx projects, employing a standard processing pipeline. GEPIA provides customizable functions such as tumor and normal differential expression analysis, profiling according to pathological stage or cancer types, patient survival analysis, similar gene detection, correlation analysis, and dimensionality reduction analysis [44]. In our study, we mainly used the correlation to reveal the relevance out of any two hub genes in HCC and normal people’s peripheral blood. Moreover, two suspicious genes that demonstrated a good manner in the scatter diagram were chosen.

4.6. Survival Analysis of Hub Genes

We used GEPIA database to analyze the relapse-free and overall survival information related to the hub genes. The hazard ratio (HR) with 95% confidence intervals and log rank value were computed and showed on the plot. was statistically significant.

4.7. Gene Set Enrichment Analysis (GSEA)

We divided 20 HCC samples from GSE49515 into two groups (high and low) on the basis of expression level of several key genes, and median expression value was regarded as the cutoff point. In order to explain the potential function of these key genes, GSEA ( was operated between the two groups. Annotated genes were selected as the reference gene sets. Nominal value < 0.05, gene size ≥ 100, and│enrichment score (ES)│ > 0.5 were regarded as the cutoff criteria.

4.8. Identification of Biomarkers

Based on the information in the individual MCODE modules, the node with the highest score was selected as the hub gene in GSE49515. According to our analysis of KEGG and PPI, we confirm four key genes (JUN, TLR4, VEGFA, and TNFSF10) which can be the noninvasive biomarkers in early HCC. To verify the vital role of the four key genes in HCC, the mRNA level of these four genes were measured in the PBMC of 20 HCC patients and 20 healthy individuals. Primers used for PCR are listed as follows: JUN-F: ACAGAGCATGACCCTGAACCT; JUN-R: TGTGCCCGTTGCTGGAC; TLR4-F: ATCAAGGACCAGAGGC; TLR4-R: CACTGAGGACCGACAC; VEGFA-F: GACGGACAGACAGACAGACACC; VEGFA-R: GAAGCGAGAACAGCCCAGA; TNFSF10-F: AGTGGCATTGCTTGTTT; TNFSF10-R: GAGCTGACGGAGTTGC.

4.9. Reidentification of Immune Subversion in HCC

Meantime, in order to confirm the key role of immune subversion in HCC, we screened peripheral blood mononuclear cells (PBMCs) using flow cytometry (FCM), which were obtain from 20 patients with HCC and 20 healthy individuals (ethical approval for the study was obtained from the Zhongnan Hospital of Wuhan University (no. 2019016)). Thereafter, we detected the level of some typical immune cells in PBMCs of HCC and healthy people via the molecules on the immune cell surface using BD FACSCanto II. The samples we detected were from anticoagulated blood with BD Multi TEST IMK Kit (Catalog No. 340503). The total level of T-lymphocytes was screened through CD3 molecules which includes helper lymphocyte T subsets detected by CD3 and CD4 molecules and inhibitory T lymphocyte subsets detected by CD3 and CD8 molecules. B lymphocytes and natural killer (NK) cells were determined by CD19, CD16, and CD56 molecules, respectively. This study was carried out according to legal requirements and supported by the Ethics Committee of the Zhongnan Hospital of Wuhan University. The HCC patients we chose are stage 0 and A of the Barcelona Clinic Liver Cancer (BCLC) stage as well as T1 and T2 of the TNM stage. It means the HCC patients are in the early stage of HCC. Besides, all patients did not take any invasive therapies, like tumor excision, interventional therapy, and hepatic artery embolism. The clinical characteristics of the participants we choose are summarized in Table 6.

GroupHepatocellular carcinoma (HCC) (n = 20)Healthy individuals

Age55 (42–76)43 (20–65)
Alpha-fetoprotein (AFP)
 >20 ng/ml200
 <20 ng/ml020
Tumor size5.57 (2.7–8.6)
 ≤3 cm2
 3–5 cm9
 >5 cm7
 Not available2
Barcelona Clinic Liver Cancer (BCLC) stage
 Stage 00
 Stage A9
 Stage B7
 Stage C0
 Stage D0
 Not available4
Viral infection
TNM stage
 Not available0

4.10. Statistical Analysis

All values were reported as means ± SD. Statistical significance was analyzed by SPSS 19.0 software. Differences were considered significant when .

Data Availability

The figure and table data used to support the findings of this study are included within the article.

Conflicts of Interest

All authors declare that they have no conflicts of interest to state.

Authors’ Contributions

Ling Li and Huijia Zhao contributed equally to this work.


This study was supported by the National Natural Science Foundation of China (nos. 81673503 and 30973582), the Open Research Fund Program of the State Key Laboratory of Virology of China (no. 2018KF005), and the Medical Science Advancement Program (Clinical Medicine) of Wuhan University (no. TFLC 2018003)

Supplementary Materials

Supplement 1: the scatter plots of GO function- and KEGG pathway-enrichment plot of HCC. (A) and (B) GO-pathway enrichment scatter plot of HCC. (C) KEGG pathway-enrichment scatter plots of HCC. (Supplementary Materials)


  1. R. X. Zhu, W. K. Seto, C. L. Lai, and M.-F. Yuen, “Epidemiology of hepatocellular carcinoma in the Asia-Pacific region,” Gut and Liver, vol. 10, no. 3, pp. 332–339, 2016. View at: Publisher Site | Google Scholar
  2. T. F. Greten, X. W. Wang, and F. Korangy, “Current concepts of immune based treatments for patients with HCC: from basic science to novel treatment approaches,” Gut, vol. 64, no. 5, pp. 842–848, 2015. View at: Publisher Site | Google Scholar
  3. H. S. Kim, K. S. Lee, H. J. Bae et al., “MicroRNA-31 functions as a tumor suppressor by regulating cell cycle and epithelial-mesenchymal transition regulatory proteins in liver cancer,” Oncotarget, vol. 6, no. 10, pp. 8089–8102, 2015. View at: Publisher Site | Google Scholar
  4. H. Yao, N. Liu, M. C. Lin, and J. Zheng, “Positive feedback loop between cancer stem cells and angiogenesis in hepatocellular carcinoma,” Cancer Letters, vol. 379, no. 2, pp. 213–219, 2016. View at: Publisher Site | Google Scholar
  5. T. Hato, L. Goyal, T. F. Greten, G. D. Duda, and X. A. Zhu, “Immune checkpoint blockade in hepatocellular carcinoma: current progress and future directions,” Hepatology, vol. 60, no. 5, pp. 1776–1782, 2014. View at: Publisher Site | Google Scholar
  6. M. Sprinzl and P. Galle, “Immune control in hepatocellular carcinoma development and progression: role of stromal cells,” Seminars in Liver Disease, vol. 34, no. 4, pp. 376–388, 2014. View at: Publisher Site | Google Scholar
  7. Y. Zhang, S. Petropoulos, J. Liu et al., “The signature of liver cancer in immune cells DNA methylation,” Clinical Epigenetics, vol. 10, no. 1, p. 8, 2018. View at: Publisher Site | Google Scholar
  8. T. Peeraphatdit, N. Naksuk, P. Phatharacharukul, B. J. Bell, and P. Ricci, “Adherence to diagnostic guidelines of hepatocellular carcinoma,” European Journal of Gastroenterology & Hepatology, vol. 27, no. 7, pp. 846–852, 2015. View at: Publisher Site | Google Scholar
  9. J.-Y. Huang, Q. Lu, and J.-B. Liu, “Delayed hepatic rupture post ultrasound-guided percutaneous liver biopsy,” Medicine, vol. 97, no. 9, p. e9955, 2018. View at: Publisher Site | Google Scholar
  10. B. I. Carr, V. Guerra, E. G. Giannini et al., “Significance of platelet and AFP levels and liver function parameters for HCC size and survival,” The International Journal of Biological Markers, vol. 29, no. 3, pp. 215–223, 2014. View at: Publisher Site | Google Scholar
  11. S. She, Y. Xiang, M. Yang et al., “C-reactive protein is a biomarker of AFP-negative HBV-related hepatocellular carcinoma,” International Journal of Oncology, vol. 47, no. 2, pp. 543–554, 2015. View at: Publisher Site | Google Scholar
  12. M. Liu, Z. Xu, Z. Du et al., “The identification of key genes and pathways in glioma by bioinformatics analysis,” Journal of Immunology Research, vol. 2017, Article ID 1278081, 9 pages, 2017. View at: Publisher Site | Google Scholar
  13. J. Wang, Y. Xie, X. Bai et al., “Targeting dual specificity protein kinase TTK attenuates tumorigenesis of glioblastoma,” Oncotarget, vol. 9, no. 3, pp. 3081–3088, 2018. View at: Publisher Site | Google Scholar
  14. J. C. Kasemeier-Kulesa, S. Schnell, T. Woolley et al., “Predicting neuroblastoma using developmental signals and a logic-based model,” Biophysical Chemistry, vol. 238, pp. 30–38, 2018. View at: Publisher Site | Google Scholar
  15. L. Li, L. Cai, L. Zheng et al., “Gefitinib inhibits bleomycin-induced pulmonary fibrosis via alleviating the oxidative damage in mice,” Oxidative Medicine and Cellular Longevity, vol. 2018, Article ID 8249693, 12 pages, 2018. View at: Publisher Site | Google Scholar
  16. W. Okajima, S. Komatsu, D. Ichikawa et al., “Liquid biopsy in patients with hepatocellular carcinoma: circulating tumor cells and cell-free nucleic acids,” World Journal of Gastroenterology, vol. 23, no. 31, pp. 5650–5668, 2017. View at: Publisher Site | Google Scholar
  17. K. Endo-Umeda, H. Nakashima, S. Komine-Aizawa et al., “Liver X receptors regulate hepatic F4/80+ CD11b+ Kupffer cells/macrophages and innate immune responses in mice,” Scientific Reports, vol. 8, no. 1, p. 9281, 2018. View at: Publisher Site | Google Scholar
  18. C. Schneider, A. Teufel, T. Yevsa et al., “Adaptive immunity suppresses formation and progression of diethylnitrosamine-induced liver cancer,” Gut, vol. 61, no. 12, pp. 1733–1743, 2012. View at: Publisher Site | Google Scholar
  19. T. Mgrditchian, T. Arakelian, J. Paggetti et al., “Targeting autophagy inhibits melanoma growth by enhancing NK cells infiltration in a CCL5-dependent manner,” Proceedings of the National Academy of Sciences, vol. 114, no. 44, pp. E9271–E9279, 2017. View at: Publisher Site | Google Scholar
  20. D. Gotthardt, E. M. Putz, E. Grundschober et al., “STAT5 is a key regulator in NK cells and acts as a molecular switch from tumor surveillance to tumor promotion,” Cancer Discovery, vol. 6, no. 4, pp. 414–429, 2016. View at: Publisher Site | Google Scholar
  21. J. Wagner, C. L. Kline, L. Zhou et al., “Dose intensification of TRAIL-inducing ONC201 inhibits metastasis and promotes intratumoral NK cell recruitment,” Journal of Clinical Investigation, vol. 128, no. 6, pp. 2325–2338, 2018. View at: Publisher Site | Google Scholar
  22. M. Ke, H. Wang, M. Zhang et al., “The anti-lung cancer activity of SEP is mediated by the activation and cytotoxicity of NK cells via TLR2/4 in vivo,” Biochemical Pharmacology, vol. 89, no. 1, pp. 119–130, 2014. View at: Publisher Site | Google Scholar
  23. M. F. Mian, N. M. Lauzon, D. W. Andrews, B. D. Lichty, and A. A. Ashkar, “FimH can directly activate human and murine natural killer cells via TLR4,” Molecular Therapy, vol. 18, no. 7, pp. 1379–1388, 2010. View at: Publisher Site | Google Scholar
  24. Y.-W. Tsao, Y.-C. Kuan, J.-L. Wang, and F. Sheu, “Characterization of a novel maitake (Grifola frondosa) protein that activates natural killer and dendritic cells and enhances antitumor immunity in mice,” Journal of Agricultural and Food Chemistry, vol. 61, no. 41, pp. 9828–9838, 2013. View at: Publisher Site | Google Scholar
  25. U. Pannicke, B. Baumann, S. Fuchs et al., “Deficiency of innate and acquired immunity caused by an IKBKB mutation,” New England Journal of Medicine, vol. 369, no. 26, pp. 2504–2514, 2013. View at: Publisher Site | Google Scholar
  26. K. Sabapathy, T. Kallunki, J. P. David, I. Graef, M. Karin, and E. F. Wagner, “C-Jun NH2-terminal kinase JNK1 and JNK2 have similar and stage-dependent roles in regulating T cell apoptosis and proliferation,” Journal of Experimental Medicine, vol. 193, no. 3, pp. 317–328, 2001. View at: Publisher Site | Google Scholar
  27. H. C. K. Patterson, M. Kraus, Y.-M. Kim, H. Ploegh, and K. Rajewsky, “The B cell receptor promotes B cell activation and proliferation through a non-ITAM tyrosine in the Igα cytoplasmic domain,” Immunity, vol. 25, no. 1, pp. 55–65, 2006. View at: Publisher Site | Google Scholar
  28. T. Voron, O. Colussi, E. Marcheteau et al., “VEGF-A modulates expression of inhibitory checkpoints on CD8+ T cells in tumors,” Journal of Experimental Medicine, vol. 212, no. 2, pp. 139–148, 2015. View at: Publisher Site | Google Scholar
  29. M. Terme, S. Pernot, E. Marcheteau et al., “VEGFA-VEGFR pathway blockade inhibits tumor-induced regulatory T-cell proliferation in colorectal cancer,” Cancer Research, vol. 73, no. 2, pp. 539–549, 2013. View at: Publisher Site | Google Scholar
  30. M. R. Pillai, L. W. Collison, X. Wang et al., “The plasticity of regulatory T cell function,” The Journal of Immunology, vol. 187, no. 10, pp. 4987–4997, 2011. View at: Publisher Site | Google Scholar
  31. Y.-T. Oh, J. Deng, P. Yue, T. K. Owonikoko, F. R. Khuri, and S.-Y. Sun, “Inhibition of B-Raf/MEK/ERK signaling suppresses DR5 expression and impairs response of cancer cells to DR5-mediated apoptosis and T cell-induced killing,” Oncogene, vol. 35, no. 4, pp. 459–467, 2016. View at: Publisher Site | Google Scholar
  32. H. Fang, B. Ang, X. Xu et al., “TLR4 is essential for dendritic cell activation and anti-tumor T-cell response enhancement by DAMPs released from chemically stressed cancer cells,” Cellular & Molecular Immunology, vol. 11, no. 2, pp. 150–159, 2014. View at: Publisher Site | Google Scholar
  33. A. Ahmed, J. H. Wang, and H. P. Redmond, “Silencing of TLR4 increases tumor progression and lung metastasis in a murine model of breast cancer,” Annals of Surgical Oncology, vol. 20, no. S3, pp. 389–396, 2013. View at: Publisher Site | Google Scholar
  34. N. Li, L. Li, and Y. Chen, “The identification of core gene expression signature in hepatocellular carcinoma,” Oxidative Medicine and Cellular Longevity, vol. 2018, Article ID 3478305, 15 pages, 2018. View at: Publisher Site | Google Scholar
  35. N. Coppola, L. Onorato, V. Iodice et al., “Occult HBV infection in HCC and cirrhotic tissue of HBsAg-negative patients: a virological and clinical study,” Oncotarget, vol. 7, no. 38, pp. 62706–62714, 2016. View at: Publisher Site | Google Scholar
  36. Y.-Q. Wang, X. Ma, L. Lu et al., “Defective antiviral CD8 T-cell response and viral clearance in the absence of c-Jun N-terminal kinases,” Immunology, vol. 142, no. 4, pp. 603–613, 2014. View at: Publisher Site | Google Scholar
  37. R. Mandraju, S. Murray, J. Forman et al., “Differential ability of surface and endosomal TLRs to induce CD8 T cell responses in vivo,” The Journal of Immunology, vol. 192, no. 9, pp. 4303–4315, 2014. View at: Publisher Site | Google Scholar
  38. A. Ward, A. Balwierz, J. D. Zhang et al., “Re-expression of microRNA-375 reverses both tamoxifen resistance and accompanying EMT-like properties in breast cancer,” Oncogene, vol. 32, no. 9, pp. 1173–1182, 2013. View at: Publisher Site | Google Scholar
  39. S. Luo, N. Cao, Y. Tang, and W. Gu, “Identification of key microRNAs and genes in preeclampsia by bioinformatics analysis,” PLoS One, vol. 12, no. 6, Article ID e0178549, 2017. View at: Publisher Site | Google Scholar
  40. P. Gaudet, N. Škunca, J. C. Hu, and C. Dessimoz, “Primer on the gene ontology,” Methods in Molecular Biology, vol. 1446, pp. 25–37, 2017. View at: Publisher Site | Google Scholar
  41. M. Kanehisa, “The KEGG database,” Novartis Foundation Symposium, vol. 247, pp. 91–101, 2002. View at: Google Scholar
  42. D. Huang, B. T. Sherman, Q. Tan et al., “The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists,” Genome Biology, vol. 8, no. 9, p. R183, 2007. View at: Publisher Site | Google Scholar
  43. D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, no. D1, pp. D447–D452, 2015. View at: Publisher Site | Google Scholar
  44. Z. Tang, C. Li, B. Kang, G. Gao, C. Li, and Z. Zhang, “GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses,” Nucleic Acids Research, vol. 45, no. W1, pp. W98–W102, 2017. View at: Publisher Site | Google Scholar

Copyright © 2019 Ling Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

1003 Views | 416 Downloads | 1 Citation
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly and safely as possible. Any author submitting a COVID-19 paper should notify us at to ensure their research is fast-tracked and made available on a preprint server as soon as possible. We will be providing unlimited waivers of publication charges for accepted articles related to COVID-19. Sign up here as a reviewer to help fast-track new submissions.