Identification of significant gene and pathways involved in HBV-related hepatocellular carcinoma by bioinformatics analysis

Introduction

Hepatocellular carcinoma (HCC) is the sixth most commonly diagnosed type of malignant tumor and is the second leading cause of cancer-related deaths worldwide. It is estimated that there were about 841,000 new cases and 782,000 deaths caused by liver cancer worldwide in 2018 (Bray et al., 2018), with Chinese patients making up more than half of the global HCC burden (Jemal et al., 2011; Bray et al., 2018). The high incidence of HCC in parts of Asia is mainly due to the prevalence of hepatitis B virus and C virus infections, especially the hepatitis B virus (Bray et al., 2018). Accumulating evidence has shown that carcinogenesis and progression of HCC are closely related to overexpression of various oncogenes and inactivation of tumor suppressor genes.

The poor prognosis associated with HCC is attributed to the lack of effective diagnostic and therapeutic methods in the early stage of the disease. In recent years, gene targeting therapy has been increasingly used for the treatment of advanced HCC, and significant progress has been made. The most commonly reported genetic alterations in HCC include mutations in the TERT promoter, TP53, CTNB1, AXIN1, ARID1A, CDKN2A, ARID2, RPS6KA3, and CCND1 (Khemlina, Ikeda & Kurzrock, 2017). Sorafenib targets multiple kinases and has been approved by the US Food and Drug Administration for the treatment of advanced HCC (Zhu et al., 2017). However, sorafenib has many shortcomings, such as low efficiency, high cost, and multiple side effects (Hu et al., 2013). Therefore, there is an urgent need to explore the relationship between the new gene function and the occurrence, development, and malignant characteristics of HCC, as well as to elucidate the precise molecular mechanisms underlying HCC, develop early screening methods, and discover novel and effective therapeutic strategies.

Recently, high-throughput technologies and gene chips have served as rapid methods for the identification of differentially expressed genes (DEGs) and functional pathways involved in the initiation and development of various diseases (Roh et al., 2010; Vogelstein et al., 2013). In these studies, a large number of tumor samples can be analyzed and thousands of genes can be identified; as a result, bioinformatics methods have become necessary for the analysis of gene expression profiles. However, obtaining reliable results from a single gene expression profile data is difficult, considering the potentially large number of differentially expressed genes, lack of stability and reproducibility, and high false-positive rates.

The RobustRankAggreg (RRA), which is based on a statistical model, is a biological analysis method for the integration and analysis of multiple gene lists (Kolde et al., 2012). RRA is a rank aggregation algorithm that assumes that all genes are arranged randomly in each dataset. A gene with a higher ranking in all datasets has a lower P-value and has a higher likelihood of being a DEG. Compared to other strategies used for the meta-analysis of datasets from multiple databases, the RRA method is more robust and easier to compute and facilitates better evaluation of the significance of the results. In addition, the RRA algorithm can handle the variable number of genes identified from different microarray platforms. More importantly, the RRA method does not strictly require the use of certain subset of problems or complete datasets to produce highly reliable results (Kolde et al., 2012). Therefore, the RRA algorithm is highly suitable for the integrated analysis of datasets from multiple databases.

In the present study, four GSE datasets GSE19665 (Deng et al., 2010), GSE55092 (Melis et al., 2014), GSE94660 (Yoo et al., 2017), and GSE121248 (Wang, Ooi & Hui, 2007) were downloaded from GEO; these datasets comprise a total of 299 samples, including 145 hepatitis B virus (HBV)-related HCC tissues and 154 non-cancerous tissues (chronic hepatitis B patients). The chip probe IDs were converted to their corresponding gene symbols. Bioinformatics analysis using R software and RRA method was then performed to obtain the integrated differentially expressed genes (DEGs). The Gene Ontology (GO; http://www.geneontology.org) is a public bioinformatics resource that provides information about gene product function using ontologies to represent biological knowledge (Gaudet & Dessimoz, 2017). KEGG (Kyoto Encyclopedia of Genes and Genomes) is a knowledgebase used for the systematic analysis of gene functions and for linking genomic information with higher-order functional information (Kanehisa & Goto, 2000). Herein, enriched GO terms and KEGG pathways were identified using the online tool DAVID 6.7. Then, the protein-protein interaction (PPI) network of the DEGs was constructed, and the hub genes were identified. We constructed the network using the hub genes and their co-expressed genes and analyzed the biological processes associated with the hub genes. Finally, survival analysis was performed based on the hub DEGs by generating the Kaplan–Meier curves in the cBioPortal. Therefore, the hub DEGs and the associated enriched pathways identified in this study can serve as reliable molecular markers for HBV-related HCC.

Material and methods

Results

Identification of DEGs in HCC

Four HBV-related HCC gene expression profiles were downloaded from the NCBI GEO database. Afterwards, the gene expression data was normalized and DEGs were identified with the limma R package (adjusted P < 0.05 and —log fold change (FC)— > 1), and the results are shown in Fig. S1. We screened out 648, 1,043, 1,171, and 580 DEGs respectively (Table 2, Fig. 1, Table S1). Through the integration and analysis of RRA, a total of integrated 341 DEGs were identified from the four datasets, including 117 upregulated genes and 224 downregulated genes (Table S2). The top 20 upregulated genes and the top 20 downregulated genes were charted on a heat map, as shown in Fig. 2.

Table 2:

Information of DEGs screened from each dataset.

GEO	Sample	Number of DEGs	Number of upregulated genes	Number of downregulated genes
GSE19665	hepatocellular carcinoma (HBV)	648	257	391
GSE55092	hepatocellular carcinoma (HBV)	1,034	409	634
GSE94660	hepatocellular carcinoma (HBV)	1,171	360	811
GSE121248	hepatocellular carcinoma (HBV)	580	167	413

DOI: 10.7717/peerj.7408/table-2

Notes:

GEO

gene expression omnibus

Functional enrichment analyses of DEGs

To further investigate the biological functions of the 314 DEGs, GO analysis was performed using online database DAVID 6.7. As shown in Figs. 3A and 3B and Table 3, GO analysis can be divided into three functional groups: biological process group (BP), the cellular component group (CC), and the molecular function group (MF).The results of GO analysis exhibited that the integrated DEGs were particularly enriched in the BP, including cell cycle, M phase, cell cycle phase, mitosis and nuclear division for the upregulated DEGs and oxidation reduction, innate immune response, complement activation, response to wounding and activation of plasma proteins involved in acute inflammatory response for the downregulated DEGs. For the CC, the upregulated DEGs were mainly enriched in mitotic spindle, microtubule cytoskeleton, cytoskeletal part, cytoskeleton and spindle pole and the downregulated DEGs were enriched in extracellular region, extracellular region part, extracellular space, microsome, vesicular fraction. In the MF, the upregulated DEGs were principally enriched in ATP binding, adenyl ribonucleotide binding, adenyl nucleotide binding, purine nucleoside binding, nucleoside binding, and the downregulated DEGs were enriched in electron carrier activity, oxygen binding, iron ion binding, heme binding, tetrapyrrole binding.

Table 3:

Top 15 GO enrichment terms of differentially expressed genes associated with hepatitis B-related hepatocellular carcinoma.

Expression	Category	Term	Count	%	P value
upregulated	BP	GO:0007049∼cell cycle	43	38.39286	7.63E-28
	BP	GO:0000279∼M phase	32	28.57143	2.48E-27
	BP	GO:0022403∼cell cycle phase	33	29.46429	1.59E-25
	BP	GO:0007067∼mitosis	27	24.10714	2.03E-25
	BP	GO:0000280∼nuclear division	27	24.10714	2.03E-25
	CC	GO:0005819∼ mitotic spindle	21	18.75	1.43E-21
	CC	GO:0015630∼microtubule cytoskeleton	24	21.42857	2.84E-13
	CC	GO:0044430∼cytoskeletal part	25	22.32143	3.19E-09
	CC	GO:0005856∼cytoskeleton	28	25	5.91E-08
	CC	GO:0000922∼spindle pole	7	6.25	6.75E-08
	MF	GO:0005524∼ATP binding	24	21.42857	3.45E-05
	MF	GO:0032559∼adenyl ribonucleotide binding	24	21.42857	4.27E-05
	MF	GO:0030554∼adenyl nucleotide binding	24	21.42857	9.64E-05
	MF	GO:0001883∼purine nucleoside binding	24	21.42857	1.22E-04
	MF	GO:0001882∼nucleoside binding	24	21.42857	1.35E-04
downregulated	BP	GO:0055114∼oxidation reduction	36	16.66667	7.03E-14
	BP	GO:0045087∼innate immune response	15	6.944444	1.67E-09
	BP	GO:0006956∼complement activation	10	4.62963	1.78E-09
	BP	GO:0009611∼response to wounding	27	12.5	1.84E-09
	BP	GO:0002541∼activation of plasma proteins involved in acute inflammatory response	10	4.62963	2.23E-09
	CC	GO:0005576∼extracellular region	83	38.42593	2.05E-21
	CC	GO:0044421∼extracellular region part	50	23.14815	9.40E-16
	CC	GO:0005615∼extracellular space	41	18.98148	8.87E-15
	CC	GO:0005792∼microsome	17	7.87037	2.50E-07
	CC	GO:0042598∼vesicular fraction	17	7.87037	3.70E-07
	MF	GO:0009055∼electron carrier activity	23	10.64815	1.17E-13
	MF	GO:0019825∼oxygen binding	12	5.555556	5.61E-12
	MF	GO:0005506∼iron ion binding	22	10.18519	5.88E-10
	MF	GO:0020037∼heme binding	14	6.481481	6.05E-09
	MF	GO:0046906∼tetrapyrrole binding	14	6.481481	1.33E-08

DOI: 10.7717/peerj.7408/table-3

Notes:

BP

biological process

CC

cellular component

MF

molecular function

GO

gene ontology

Integration of PPI network and module analysis

PPI network of the 341 DEGs was established by STRING. A total of 288 DEGs (99 upregulated genes and 189 downregulated genes) were filtered into the DEG PPI network complex, which contained 288 nodes and 2,259 edges (Fig. S2, Tables S3 and S4). Among the 288 nodes, 59 central node genes were selected as hub genes with the criteria of filtering degree ≥ 10 (ie, each node has more than 10 connections/interactions). The top 10 hub genes were as follows: SHCBP1, FOXM1, KIF4A, ANLN, KIF15, KIF18A, FANCI, NEK2, ECT2 and RAD51AP1 (Table 5).

Table 5:

The top 10 most degree values hub genes between HBV-related HCC and normal samples.

Up indicated that the gene was identifed as up-regulated in HCC; Down indicated that the gene was reported as down-regulated. UN suggested the gene has not been reported in current HCC associated studies.

Gene symbol	Gene description	logFC	Degree	Up/down
SHCBP1	SHC binding and spindle associated 1	1.211943105	43	Up
FOXM1	forkhead box M1	1.832558689	43	Up
KIF4A	kinesin family member 4A	1.850277648	43	Up
ANLN	anillin actin binding protein	1.982343654	43	Up
KIF15	kinesin family member 15	1.247277869	43	UN
KIF18A	kinesin family member 18A	1.26437819	43	Up
FANCI	FA complementation group I	1.269987824	43	UN
NEK2	NIMA related kinase 2	1.652396787	43	Down
ECT2	epithelial cell transforming 2	1.418475417	43	Up
RAD51AP1	RAD51 associated protein 1	1.585151756	42	UN

DOI: 10.7717/peerj.7408/table-5

Notes:

FC

fold change

Furthermore, the two most significant modules (Figs. 4A and 4B, Tables S5 and S6) of the PPI network were selected for KEGG pathway enrichment analysis. Results showed that the genes in module 1 were mainly enriched in cell cycle, oocyte meiosis, p53 signaling pathway, progesterone-mediated oocyte maturation, and the genes in module 2 were mainly enriched in complement and coagulation cascades, prion diseases and systemic lupus erythematosus (Table 6).

Discussion

HCC is a common malignant tumor affecting the digestive system. The incidence of HCC in developing countries is considerably higher than that in developed countries. Chronic infection with HBV or hepatitis C virus (HCV) is the primary etiological factor related to HCC in certain parts of Asia and sub-Saharan Africa (Bray et al., 2018). In western countries, the main risk factors include obesity, type 2 diabetes, heavy alcohol consumption, and smoking (El-Serag, 2011; Mittal & El-Serag, 2013). Although significant progress has been achieved in the diagnosis and treatment of HCC in recent years, the prognosis of HCC remains poor (Okajima et al., 2017). Thus, there is an urgent need to understand the detailed molecular mechanisms underlying HCC for early detection, diagnosis, and treatment of the disease. Recently, the wide application of microarray and high-throughput technologies has provided an effective way to screen thousands of genes that are likely to be involved in the occurrence and development of HCC. These genes could serve as potential targets for the diagnosis and treatment of HCC.

In the present study, we analyzed four mRNA expression datasets and identified a total of 341 DEGs, comprising 117 upregulated genes and 224 downregulated genes, in HBV-related HCC samples. GO and KEGG enrichment analyses showed that the upregulated genes were associated with the cell cycle (ontology: BP), mitotic spindle (ontology: CC), and adenosine triphosphate binding, adenyl ribonucleotide binding, adenyl nucleotide binding, purine nucleoside binding, and nucleoside binding (ontology: MF). The majority of the downregulated genes were enriched in oxidation reduction (ontology: BP), extracellular region (ontology: CC), and electron carrier activity (ontology: MF). The above results suggested that these DEGs are involved in the proliferation and cell cycle of chronic hepatitis B-induced HCC cells. KEGG pathway analysis revealed that the integrated DEGs were significantly enriched in retinol metabolism, drug metabolism, metabolism of xenobiotics by cytochrome P450, caffeine metabolism, and tryptophan metabolism. Previous studies have shown that dysregulation of the cell cycle and oxidation reduction play vital roles in the initiation or progression of tumors (Choi et al., 2001; Wang et al., 2017). Cytochrome P450 enzymes are involved in various types of cancer via several mechanisms, including the catalysis of the bioactivation of chemical procarcinogens, participation in the activation of cancer therapeutic drugs, as targets for cancer treatment, and as metabolic enzymes (Hrycay & Bandiera, 2015).

After screening out two modules with the most significant interactions, ten genes with the highest degrees of interaction were subsequently identified by constructing the PPI network. Results of survival analysis showed that the patients with alterations in FOXM1, NEK2, RAD51AP1, ANLN, and KIF18A expression patterns showed worse prognosis.

ANLN, an actin-binding protein, has been reported to be dysregulated in diverse human tumors (Hall et al., 2005). ANLN is not only considered as a proliferative marker, but also a prognostic and therapeutic indicator. ANLN knockdown was demonstrated to inhibit cell growth and migration in breast cancer (Zhou et al., 2015). In addition, another study showed that high ANLN levels in the nuclear fraction are involved in the cell cycle and are correlated with poor prognosis (Magnusson et al., 2016). Upregulation of ANLN expression plays an important role in non-small cell lung cancers (NSCLC) by activating RHOA and participating in the phosphoinositide 3-kinase/AKT pathway (Suzuki et al., 2005). A previous study showed that ANLN expression in gastric cancer (GC) is a molecular predictor of intestinal and proliferative type gastric tumors (Pandi et al., 2014). Furthermore, ANLN was identified as a biomarker for the prognosis of bladder urothelial carcinoma (Zeng et al., 2017), colorectal cancer (Wang et al., 2016), and lung adenocarcinoma (Long et al., 2018).

ANLN mRNA levels were upregulated in human HCC tissues compared to non-tumor liver tissues. ANLN knockdown was demonstrated to slow down the growth rates of tumors, reduce the number of tumors, and prolong the survival of mice (Zhang et al., 2018). Consistent with the findings reported in previous studies (Lian et al., 2018), our results showed that higher ANLN expression levels are associated with worse clinical outcomes and a shorter survival times of patients with HCC, thereby highlighting the potential use of ANLN as a prognostic biomarker.

As a member of the kinesin-8 family, KIF18A plays crucial roles in regulating microtubule dynamics, chromosome congression, and cell division (Mayr et al., 2007). In fact, elevated KIF18A expression was observed in several human cancers. KIF18A overexpression in human breast cancer has been closely associated with tumor grade, metastasis, and poor survival (Zhang et al., 2010; Kasahara et al., 2016). Some researchers even advocate measurement of KIF18A levels in patients with estrogen receptor positive (ER+) breast cancer (BC) prior to receiving endocrine therapy (Alfarsi et al., 2019). KIF18A expression levels were found to positively contribute to tumor stage, lymphatic invasion, lymph node metastasis venous invasion, and peritoneal dissemination in CRC (Nagahara et al., 2011). Proteomic analysis indicated that KIF18A is a promising biomarker for the early diagnosis of cholangiocarcinoma (CCA) (Rucksaken et al., 2012).

KIF18A expression levels were found to be markedly higher in liver cancer tissues compared to adjacent normal liver tissues (Liao et al., 2014). KIF18A has been suggested to promote proliferation, invasion, and metastasis of HCC cells by activating cell cycle signaling pathway and the Akt and MMP-7/MMP-9-related signaling pathways (Luo et al., 2018). KIF18A levels were found to be significantly related to clinicopathologic factors associated with alpha-fetoprotein (AFP) concentrations (≥200 ng/mL), tumor size (≥5 cm), clinical tumor-node-metastasis (TNM) stage, and portal vein tumor thrombus (PVTT). In survival analysis, TCGA, Provisional higher KIF18A expression had worse prognosis (shorter DFS and OS) (Liao et al., 2014). The above results indicated that KIF18A could serve as a novel biomarker for the diagnosis and treatment of HCC. Our findings are consistent with previous studies and demonstrated that previous experiments did not show whether these HCC were associated with HBV. Therefore, the role of KIF18A in HBV-related HCC types should be verified by further experiments.

As a member of the forkhead box (Fox) transcription factor family, FOXM1 is acts as an oncogene in many tumors, such as breast, cervix, and prostate cancers, and is known to play crucial roles in the prognosis and chemoresistance of tumors (Zhu et al., 2018). FOXM1 mRNA levels were upregulated in human HCC tissues and had positive relevance to the development, metastasis, recurrence, and worse clinical outcomes in HCC patients after orthotopic liver transplantation (Sun et al., 2011; Dai et al., 2015). SHCBP1, KIF4A, and ECT2 have been reported to mediate tumor initiation and progression of human HCC (Tao et al., 2013; Chen et al., 2015; Hou et al., 2017). NEK2 could serve as a useful predictor and potential therapeutic target in HCC (Fu et al., 2017). However, previous studies reported low NEK2 expression levels, which are inconsistent with our current findings. Other research groups reported that abnormal KIF15 levels were evidently associated with HCC progression and prognosis (Chen et al., 2017); however, these findings were not verified by cell or animal experiments. Similarly, these studies did not show that the occurrence of these HCCs is closely related to HBV.

FANCI and RAD51AP1 have been identified as new markers for HBV-related HCC, but have not been widely reported based on literature retrieval. Some findings provided new insights that RAD51AP1 is likely to mediate the molecular mechanisms underlying HCV-induced pathogenesis (Nguyen et al., 2018). However, further studies are required to verify the exact roles of these two genes.

In accordance with our findings, previous studies have also identified DEGs that participate in HBV-related liver cancer (Zhou et al., 2017). For example, Zhou et al. analyzed the gene expression profiles of GSE14520 and HCC samples from the Zhongshan Hospital affiliated with Fudan University, which comprised 63 paired HCC and non-tumor samples. All patients of these two cohorts were infected with hepatitis B virus. A total of 965 DEGs (389 upregulated genes and 576 downregulated genes) that were differentially regulated by at least two-fold with statistical significance. HSP90AB1, RPL8, NPM1, and MCM3 were selected as the hub genes from the PPI network. Nevertheless, the study by Zhou et al. comprised relatively fewer samples (66 primary HCC tumors and paired adjacent non-tumor tissues), and the main purpose of the study was to identify copy number variation (CNV)-driven DEGs. In another study, the DEGs that are common from four datasets were visualized using a Venn diagram (Chen et al., 2019). As a result, the number of DEGs obtained using this method was relatively small (84 upregulated and 46 downregulated). In the end, the following top ten hub genes were obtained: TOP2A, RFC4, CCNB1, CDC20, CDKN3, BUB1B, CCNB2, TPX2, PEN1, and MAD2L1. Compared to the previous two studies, we analyzed four GEO datasets comprising 299 samples, and 341 DEGs (117 upregulated and 224 downregulated) were identified. In addition, data analysis was conducted using the RobustRankAggreg (RRA) method, which is highly suitable for the analysis of datasets from multiple databases. Therefore, the hub genes identified in the present study are more reliable and comprehensive.

Conclusion

In summary, by conducting an integrated bioinformatics analysis using multiple datasets, we identified DEGs and the association pathways involved in HBV-induced HCCs. In addition, we identified several key candidate genes and biological pathways that can provide a deeper and more comprehensive understanding of the occurrence and development of HCC and its association with HBV. Our findings provided valuable insights for the identification of novel biomarkers for the diagnosis and treatment of HCC.

Supplemental Information