Bioinformatics analysis of aberrantly methylated-differentially expressed genes and pathways in hepatocellular carcinoma.

AIM: To discover methylated-differentially expressed genes (MDEGs) in hepatocellular carcinoma (HCC) and to explore relevant hub genes and potential pathways.
METHODS: The data of expression profiling GSE25097 and methylation profiling GSE57956 were gained from GEO Datasets. We analyzed the differentially methylated genes and differentially expressed genes online using GEO2R. Functional and enrichment analyses of MDEGs were conducted using the DAVID database. A protein-protein interaction (PPI) network was performed by STRING and then visualized in Cytoscape. Hub genes were ranked by cytoHubba, and a module analysis of the PPI network was conducted by MCODE in Cytoscape software.
RESULTS: In total, we categorized 266 genes as hypermethylated, lowly expressed genes (Hyper-LGs) referring to endogenous and hormone stimulus, cell surface receptor linked signal transduction and behavior. In addition, 161 genes were labelled as hypomethylated, highly expressed genes (Hypo-HGs) referring to DNA replication and metabolic process, cell cycle and division. Pathway analysis illustrated that Hyper-LGs were enriched in cancer, Wnt, and chemokine signalling pathways, while Hypo-HGs were related to cell cycle and steroid hormone biosynthesis pathways. Based on PPI networks, PTGS2, PIK3CD, CXCL1, ESR1, and MMP2 were identified as hub genes for Hyper-LGs, and CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10 were hub genes for Hypo-HGs by combining six ranked methods of cytoHubba.
CONCLUSION: In the study, we disclose numerous novel genetic and epigenetic regulations and offer a vital molecular groundwork to understand the pathogenesis of HCC. Hub genes, including PTGS2, PIK3CD, CXCL1, ESR1, MMP2, CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10, can be used as biomarkers based on aberrant methylation for the accurate diagnosis and treatment of HCC.

Core tip: We explored methylated-differentially expressed genes in hepatocellular carcinoma (HCC) using a series of bioinformatics databases and tools. In total, we categorized 266 genes as hypermethylated, lowly expressed genes (Hyper-LGs) referring to endogenous and hormone stimulus, as well as 161 hypomethylated, highly expressed genes (Hypo-HGs) referring to DNA replication and metabolic process. Pathway analysis showed Hyper- LGs were mainly enriched in cancer, while Hypo-HGs were essentially related to cell cycle. Finally, we identified hub genes that might be utilized as biomarkers based on aberrant methylation, which might be useful for the accurate diagnosis and treatment of HCC.

INTRODUCTION

Hepatocellular carcinoma (HCC), as the most frequent type of liver cancer, is one of the main aggressive malignant cancers worldwide and the third leading cause of cancer-related deaths[1,2]. HCC embodies a complicated, multi-step disease, and the processes involved are related to genomic amplifications, deletions, insertions, or mutations to induce a series of epigenetic and genetic alterations. Despite significant advances in early diagnosis and interventional therapies with the development of surgical and treatment approaches, most HCC patients are usually diagnosed at an advanced stage of cancer progression with a low 5-year survival rate and poor prognosis[3,4]. Therefore, a better understanding of the molecular mechanisms and functional pathways of HCC and the development of new critical gene targets for early HCC detection are urgently needed.

Tumor epigenetics, acknowledged as inherited modifications in gene expression, encompasses DNA methylation, noncoding RNA, and histone acetylation[5]. DNA methylation is the main epigenetic modification, affecting independent loci in gene transcriptional regulation and preserving genome stability. A variety of tumors have a special deregulation signature that is characterized by aberrant DNA methylation[6]. Altered methylation in DNA sequences, including hypomethylation of oncogenes and hypermethylation of tumor suppressor genes, are regarded as a key event in carcinogenesis, including in HCC[7-9]. Thus, the detection of methylated-differentially expressed genes (MDEGs) and a better understanding of their characteristics may be useful for discovering the molecular mechanism and pathogenesis of HCC.

Previous studies have shown by analyzing profiling arrays that the pathogenesis of HCC is a complicated biological process involving epigenetic and genetic changes[10-12]. However, most of the above studies mainly focused on either gene expression or methylation data and did not perform a conjoint analysis. Methylated expressed genes can be detected concurrently by joining gene expression and methylation microarray data, thus allowing us to identify more accurately biological characteristics of HCC[13,14]. In the present study, we explored the interaction network of differentially expressed genes (DEGs) and differentially methylated genes (DMGs) along with interrelated signalling pathways in HCC by analyzing the expression profile of gene expression microarray data (GSE25097) and gene methylation microarray data (GSE57956) using bioinformatics tools. We aimed to identify novel insights into the biological characteristics and pathways of MDEGs in HCC and make notional viewpoints available for the development and progression of HCC.

MATERIALS AND METHODS

Microarray data

We identified MDEGs between adjacent non-tumor samples and HCC samples by analyzing mRNA microarray and methylation profiling datasets. One gene expression profiling dataset, GSE25097, and another gene methylation dataset, GSE57956, were downloaded from Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/). In total, 243 adjacent non-tumor samples and 268 HCC tumor samples were registered in GSE25097 (platform: GPL10687 Rosetta/Merck Human RSTA Affymetrix 1.0 microarray, Custom CDF). For the gene methylation microarray data, GSE57956 was comprised entirely of 59 adjacent non-tumor tissues and 61 HCC tumor tissues [platform: GPL8490 Illumina HumanMethylation27 BeadChip (HumanMethylation27_270596_v.1.2)].

Data processing

We used an online tool, GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), to analyze the differential expression by comparing two groups of samples across setup conditions in a GEO series. In the study, we used P < 0.05 and |fold change| > 2 as the cut-off standard to define the DEGs and DMGs. “MATCH function” was performed to categorize overlapping MDEGs between the GSE25097 and GSE57956 data sets. Finally, overlapping down-regulated and hypermethylation genes were identified as hypermethylated, lowly expressed genes (Hyper-LGs); similarly, overlapping up-regulated and hypomethylation genes were considered hypomethylated, highly expressed genes (Hypo-HGs).

Functional and pathway enrichment analysis

DAVID (the database for annotation, visualization and integrated discovery, https://david.ncifcrf.gov/) is an online tool for functional annotation and enrichment analysis to reveal biological features related to large gene lists[15]. Gene ontology (GO) analysis, including biological process, cellular component, and molecular function, is a main bioinformatics analysis method for annotating genes[16]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database used to obtain high-level functions and utilities of the biological system originated from genome sequencing or high-throughput experimental technologies[17]. GO function and KEGG pathway enrichment analyses were performed for MDEGs using DAVID. A P-value < 0.05 was considered as statistically significant.

Protein-protein interaction network generation and module analysis

We built a protein-protein interaction (PPI) network of Hyper-LGs and Hypo-HGs using the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins, http://string-db.org/) database. STRING is an online database used to predict PPI[18], which is essential for recognizing the mechanisms of cell activities at the molecular level in cancer progressions. The cut-off standard was defined as an interaction score (median confidence) of 0.4. Consequently, the PPI network was visualized by Cytoscape (http://www.cytoscape.org/), and hub genes were ranked by cytoHubba. Molecular Complex Detection (MCODE) analysis was performed to screen modules within the PPI network in Cytoscape software. A MCODE score > 4 and number of nodes > 5 were taken as the criteria to define a module.

RESULTS

Screening of MDEGs in HCC

Online analysis was performed by GEO2R software to identify DEGs or DMGs. By comparing the 1873 DEGs (676 up-regulated genes and 1197 down-regulated genes) with the 7242 DMGs (2652 hypermethylated genes and 4590 hypomethylated genes), we categorized 266 Hyper-LGs and 161 Hypo-HGs in GO, KEGG, and PPI analyses. The flowchart is presented in Figure 1.

Figure 1

Flowchart of bioinformatics analysis. DMGs: Differentially methylated genes; DEGs: Differentially expressed genes; Hyper-LGs: Hypermethylated, lowly expressed genes; Hypo-HGs: Hypomethylated, highly expressed genes.

GO functional enrichment analysis

GO enrichment analysis was performed by DAVID, and the results are shown in Table 1. For Hyper-LGs, enriched biological processes (BP) included response to endogenous and hormone stimulus, cell surface receptor linked signal transduction, and behavior. Cell component (CC) mainly displayed extracellular region and plasma membrane part, intrinsic to plasma membrane. Additionally, molecular function (MF) enrichment indicated glycosaminoglycan, pattern and polysaccharide binding, and protein tyrosine kinase activity as important related processes. Hypo-HGs were enriched in BP of DNA replication and metabolic process, cell division and cycle, and chromosome organization. CC was mainly involved in chromosome, chromatin, and and extracellular region part. With regards to MF, enrichments were focused on peptidase and enzyme inhibitor activity, phosphorus-oxygen lyase and cyclase activity, as well as cytoskeletal protein binding.

Table 1

GO enrichment analysis of methylated-differentially expressed genes related with hepatocellular carcinoma

Category		Term	Count	%	P value
Hyper-LGs	GOTERM_BP_FAT	GO:0048545~response to steroid hormone stimulus	18	6.79	1.51E-08
	GOTERM_BP_FAT	GO:0009725~response to hormone stimulus	24	9.06	3.15E-08
	GOTERM_BP_FAT	GO:0007166~cell surface receptor linked signal transduction	60	22.64	1.75E-07
	GOTERM_BP_FAT	GO:0009719~response to endogenous stimulus	24	9.06	1.89E-07
	GOTERM_BP_FAT	GO:0007610~behavior	25	9.43	6.76E-07
	GOTERM_CC_FAT	GO:0044421~extracellular region part	48	18.11	1.97E-11
	GOTERM_CC_FAT	GO:0005886~plasma membrane	107	40.38	9.53E-10
	GOTERM_CC_FAT	GO:0005576~extracellular region	68	25.66	1.10E-08
	GOTERM_CC_FAT	GO:0044459~plasma membrane part	71	26.79	3.27E-08
	GOTERM_CC_FAT	GO:0031226~intrinsic to plasma membrane	47	17.74	1.22E-07
	GOTERM_MF_FAT	GO:0005539~glycosaminoglycan binding	11	4.15	7.78E-05
	GOTERM_MF_FAT	GO:0004714~transmembrane receptor protein tyrosine kinase activity	8	3.02	8.88E-05
	GOTERM_MF_FAT	GO:0001871~pattern binding	11	4.15	1.72E-04
	GOTERM_MF_FAT	GO:0030247~polysaccharide binding	11	4.15	1.72E-04
	GOTERM_MF_FAT	GO:0004713~protein tyrosine kinase activity	11	4.15	3.15E-04
Hypo-HGs	GOTERM_BP_FAT	GO:0006260~DNA replication	10	6.25	6.22E-05
	GOTERM_BP_FAT	GO:0051301~cell division	12	7.50	8.37E-05
	GOTERM_BP_FAT	GO:0051276~chromosome organization	15	9.38	1.41E-04
	GOTERM_BP_FAT	GO:0006259~DNA metabolic process	15	9.38	2.19E-04
	GOTERM_BP_FAT	GO:0007049~cell cycle	19	11.88	2.44E-04
	GOTERM_CC_FAT	GO:0005694~chromosome	14	8.75	3.77E-04
	GOTERM_CC_FAT	GO:0044427~chromosomal part	12	7.50	1.01E-03
	GOTERM_CC_FAT	GO:0000793~condensed chromosome	7	4.38	1.34E-03
	GOTERM_CC_FAT	GO:0000785~chromatin	8	5.00	2.77E-03
	GOTERM_CC_FAT	GO:0044421~extracellular region part	19	11.88	3.60E-03
	GOTERM_MF_FAT	GO:0030414~peptidase inhibitor activity	6	3.75	1.24E-02
	GOTERM_MF_FAT	GO:0016849~phosphorus-oxygen lyase activity	3	1.88	1.65E-02
	GOTERM_MF_FAT	GO:0009975~cyclase activity	3	1.88	1.80E-02
	GOTERM_MF_FAT	GO:0004857~enzyme inhibitor activity	7	4.38	3.53E-02
	GOTERM_MF_FAT	GO:0008092~cytoskeletal protein binding	10	6.25	3.83E-02

Top five terms were listed on the basis of P value if over five terms in the category, Hyper-LGs (hypermethylated, lowly expressed genes), Hypo-HGs (hypomethylated, highly expressed genes).

KEGG pathway analysis

The results of the KEGG pathway enrichment analysis implied that Hyper-LGs demonstrated enrichment in pathways of complement and coagulation cascades, dilated cardiomyopathy, cancer, Wnt, and chemokine signalling pathways. Hypo-HGs were significantly involved in cell cycle and steroid hormone biosynthesis pathways (Table 2).

Table 2

KEGG pathway analysis of methylated-differentially expressed genes related with hepatocellular carcinoma

Category		Term	Count	%	P value	Gene
Hyper-LGs	KEGG_PATHWAY	hsa04610:Complement and coagulation cascades	7	2.64	6.62E-03	C7, CR1, CD55, THBD, MASP1, SERPINE1, PLAUR
	KEGG_PATHWAY	hsa05414:Dilated cardiomyopathy	7	2.64	2.50E-02	LAMA2, ITGA9, ADCY1, ADRB1, ITGB8, ADCY5, TGFB3
	KEGG_PATHWAY	hsa05200:Pathways in cancer	15	5.66	2.94E-02	FGFR2, PTGS2, FLT3, PIK3CD, FZD1, TGFB3, MMP2, WNT2, LAMA2, RAC2, NKX3-1, LAMC2, WNT11, HHIP, GSTP1
	KEGG_PATHWAY	hsa04310:Wnt signaling pathway	9	3.40	3.16E-02	WNT2, SFRP5, NKD2, RAC2, PRICKLE1, SFRP1, FZD1, WNT11, FOSL1
	KEGG_PATHWAY	hsa04062:Chemokine signaling pathway	10	3.77	3.97E-02	CXCL1, ADCY1, DOCK2, CCL23, RAC2, TIAM1, ADCY5, PIK3CD, CCL19, CXCL6
Hypo-HGs	KEGG_PATHWAY	hsa04110:Cell cycle	7	4.38	2.32E-03	CCNE2, E2F2, PRKDC, CDC20, MCM2, SFN, PTTG1
	KEGG_PATHWAY	hsa00140:Steroid hormone biosynthesis	3	1.88	9.07E-02	CYP17A1, CYP7A1, UGT2B11

Top five terms were listed on the basis of P value if over five terms in the category, Hyper-LGs (hypermethylated, lowly expressed genes), Hypo-HGs (hypomethylated, highly expressed genes).

PPI network construction and cytoHubba analysis

MDEGs were analyzed by STRING. Ultimately, 264 nodes and 456 edges and 159 nodes and 290 edges were established in the Hyper-LGs and Hypo-HGs networks, respectively. The PPI networks for Hyper-LGs and Hypo-HGs, as shown in Figures 2 and 3, exhibited significantly more interactions than expected with a PPI enrichment P-value < 1.0e-16. We then visualized the Hyper-LGs and Hypo-HGs network in Cytoscape and selected hub genes using cytoHubba. A total of five and six hub genes were identified for Hyper-LGs and Hypo-HGs, respectively, by overlap of the top 10 genes according to six ranked methods in cytoHubba (Tables 3 and 4). Hyper-LGs were annotated as prostaglandin-endoperoxide synthase 2 (PTGS2), phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit delta (PIK3CD), C-X-C motif chemokine ligand 1 (CXCL1), estrogen receptor 1 (ESR1), and matrix metallopeptidase 2 (MMP2). Hypo-HGs were annotated as cell division cycle 45 (CDC45), denticleless E3 ubiquitin protein ligase homolog (DTL), Aurora kinase B (AURKB), cyclin dependent kinase inhibitor 3 (CDKN3), minichromosome maintenance complex component 2 (MCM2), and minichromosome maintenance 10 replication initiation factor (MCM10).

Figure 2

Protein-protein interaction network of hypermethylated, lowly expressed genes. Disconnected nodes were hid in the network.

Figure 3

Protein-protein interaction network of hypomethylated, highly expressed genes. Disconnected nodes were hid in the network.

Table 3

Hub genes for hypermethylated, lowly expressed genes ranked in cytoHubba

Catelogy	Rank methods in cytoHubba
	MNC	Degree	EPC	Closeness	Radiality	Stress
Gene symbol top 10	PTGS2	ADCY5	PTGS2	PIK3CD	PIK3CD	PTGS2
	PIK3CD	MMP2	PIK3CD	PTGS2	PTGS2	PIK3CD
	ADCY5	PTGS2	MMP2	MMP2	MMP2	MMP2
	ADCY1	PIK3CD	ADCY5	ESR1	ESR1	PRKG1
	CXCL1	PRKG1	ADCY1	PRKG1	TLR2	ESR1
	ESR1	ADCY1	ESR1	TLR2	SERPINE1	FYN
	MMP2	ESR1	CXCL1	FYN	PRKG1	RAC2
	FYN	FYN	TLR2	CXCL1	SNAI1	SERPINE1
	TLR2	CXCL1	CALCA	SERPINE1	CRP	CXCL1
	SERPINE1	TLR2	PTGER2	ADCY5	CXCL1	ADCY5

Bold gene symbols were the overlap hub genes in top 10 by six ranked methods respectively in cytoHubba. MNC: Maximum neighborhood component; Degree: Node connect degree; EPC: Edge percolated component.

Table 4

Hub genes for hypomethylated, highly expressed genes ranked in cytoHubba

Catelogy	Rank methods in cytoHubba
	MCC	MNC	Degree	EPC	Closeness	Radiality
Gene symbol top 10	CDC45	CDC45	CDC45	CDC45	CDKN3	CDKN3
	DTL	AURKB	AURKB	AURKB	CDC45	PRKDC
	RACGAP1	DTL	CDKN3	CDKN3	AURKB	CDC45
	AURKB	RACGAP1	DTL	DTL	PTTG1	PTTG1
	CDC20	CDC20	RACGAP1	RACGAP1	DTL	MCM10
	CDKN3	CDKN3	CDC20	CDC20	MCM2	BRCA1
	RRM2	MCM2	MCM2	MCM2	MCM10	PI3
	MCM2	PTTG1	PTTG1	PTTG1	RACGAP1	AURKB
	MCM10	RRM2	RRM2	RRM2	CDC20	MCM2
	MKI67	MCM10	MCM10	MCM10	RRM2	DTL

Bold gene symbols were the overlap hub genes in top 10 by six ranked methods respectively in cytoHubba. MCC: Maximal cilque centrality; MNC: Maximum neighborhood component; Degree: Node connect degree; EPC: Edge percolated component.

Module analysis

In total, four modules in the Hyper-LGs network and three modules in the Hypo-HGs network were established as statistically significant. The following GO and KEGG pathways were analyzed (Table 5). Enrichment analyses for the Hyper-LGs modules demonstrated that the pathways are mainly associated with neuroactive ligand and ECM-receptor interaction, axon guidance, and chemokine signalling pathway. For Hypo-HGs modules, enrichment analysis showed associations with the cell cycle and chemokine signalling pathway. The visualized genes of modules in the Hyper-LGs and Hypo-HGs network are shown in Figure 4A-D and Figure 5A-C.

Table 5

Modules analysis of the protein–protein interaction network

Category	Module	Score	Nodes	Enrichment and pathway description	Genes
Hyper-LGs	1	10.00	10	GO.0005886: plasma membrane	ADRB1, VIPR1, PTGDR, SCTR, CALCA, GPR83, ADCY1, ADCY5, PTGER2, PTGER4
				GO.0007187: G-protein coupled receptor signaling pathway
				GO.0004016: adenylate cyclase activity
				has04080: Neuroactive ligand-receptor interaction
	2	6.00	6	GO.0051953: negative regulation of amine transport	CCL19, ADRA2B, P2RY12, CXCL6, NPY5R, CXCL1
				GO.0008009: chemokine activity
				has04062: Chemokine signaling pathway
	3	5.68	20	GO.0005886: plasma membrane	OXT, SERPINE1, EDNRB, PTGS2, ADRA1A, IL18, CRP, SOCS3, EFNB3, MMP2, EPHB1, TIAM1, EFNA5, TBXA2R, EPHA2, TLR2, SNAI1, FYN, GNA14, PROK2
				GO.0005003: ephrin receptor activity
				GO.0051240: positive regulation of multicellular organismal process
				hsa04360: Axon guidance
	4	4.50	5	GO.0005605: basal lamina	ITGA9, ITGB8, COL6A2, LAMC2, LAMA2
				GO.0030198: extracellular matrix organization
				hsa04512: ECM-receptor interaction
Hypo-HGs	1	17.56	19	GO.0022402: cell cycle process	CDCA5, KIF14, BRCA1, CENPF, RACGAP1, NEK2, CDKN3, DTL, MCM2, MCM10, CDC45, PTTG1, ANLN, CDC20, RRM2, AURKB, MKI67, STIL, CCNE2
				GO.0015630: microtubule cytoskeleton
				GO.0003688: DNA replication origin binding
				hsa04110: Cell cycle
	2	5.00	5	GO.0007188: adenylate cyclase-modulating G-protein coupled receptor signaling pathway	CNR1, CCL20, ADCY6, HTR1D, CCL25
				hsa04062: Chemokine signaling pathway
	3	4.50	5	GO.0005615: extracellular space	REN, PLA2G1B, TIMP1, MMP9, VWF

Hyper-LGs: Hypermethylated, lowly expressed genes; Hypo-HGs: Hypomethylated, highly expressed genes.

Figure 4

Hypermethylated, lowly expressed genes modules.

Figure 5

Hypomethylated, highly expressed genes modules.

DISCUSSION

The occurrence and development of HCC is a complex and multistage process that involves multiple molecular changes of cumulative genetic and epigenetic disorders. As with many other tumors, epigenetic disturbances contribute significantly to the etiology of HCC, especially DNA methylation. Overall, identifying biomarkers in complex diseases, such as HCC, contributes to our understanding of the pathogenesis and diagnosis of diseases[12]. In this study, we identified 266 Hyper-LGs and 161 Hypo-HGs by utilizing public datasets and online bioinformatics tools to analyze microarray profiling data of gene expression (GSE25907) and gene methylation (GSE57956) in HCC. The findings of the interaction network disclosed that the related genes may be involved in molecular regulation of important pathways associated with the development and progression of HCC. Functional and enrichment analyses of the genes verified definite pathways, as well as hub genes associated with methylation, which may offer novel viewpoints for revealing the pathogenesis of HCC.

In view of the analysis in DAVID, Hyper-LGs in HCC, GO enrichment analysis demonstrated BP included response to endogenous and hormone stimulus, cell surface receptor linked signal transduction, and behavior. One fundamental endogenous genotoxic stimulus could cause DNA damage response in cancers[19]. MF enrichment indicated glycosaminoglycan, pattern and polysaccharide binding, and protein tyrosine kinase activity. Previous studies reported that receptor tyrosine kinases restrained tumor angiogenesis and proliferation[20]. In our study, KEGG enrichment analysis revealed the involvement of complement and coagulation cascades, dilated cardiomyopathy, cancer, Wnt, and chemokine signalling. These pathways can promote tumor cell proliferation and metastasis and alter the microenvironment in the pathogenesis of HCC[20,21].

Hypo-HGs in HCC were enriched in the BP of DNA replication, the metabolic process, cell division and cycle, and chromosome organization. MF of GO analysis largely showed enrichments in peptidase and enzyme inhibitor activity, phosphorus-oxygen lyase and cyclase activity, and cytoskeletal protein binding. Previous research showed that the cell cycle played a critical role in cancer by controlling cell division, and there are significant associations among cell proliferation, cell cycle deregulation, and cell cycle-related kinase with HCC incidence and metastasis[22,23]. In addition, cell cycle and steroid hormone biosynthesis pathways were disclosed by KEGG enrichment analysis in present study. It is plausible that metabolites involved in sterol and sphingolipid biosynthesis and phosphoinositides are related to the development of HCC[24]. In summary, understanding biological processes and the signalling pathways involved in MDEGs can help elucidate the pathogenesis of HCC and identify new therapeutic targets.

Based on the PPI network generated for MDEGs, significantly more interactions than expected were observed for Hyper-LGs and Hypo-HGs, with a PPI enrichment P-value < 1.0e-16, and a number of MDEGs appear to be involved in the development and progression of HCC. Finally, we visualized the networks in Cytoscape and identified hub genes for Hyper-LGs using cytoHubba in Cytoscape software: PTGS2, PIK3CD, CXCL1, ESR1, and MMP2. PTGS2 is a proinflammatory enzyme induced by prostaglandins involved in cell proliferation, tumorigenesis, progression, and metastasis[25] . The PTGS2 gene is commonly up-regulated and plays a role in apoptosis and proliferation of cells in numerous types of cancers[26,27]. It is noteworthy that one meta-analysis showed that the HCC susceptibility is associated with a PTGS2 variant[28]. PIK3CD is a protein coding gene related to transferase activity, transferring phosphorus-containing groups and kinase activity, and involves in lymphocyte activation, proliferation and differentiation[29]. One study found that PIK3CD activating mutations can cause immunodeficiency in patients with a complex phenotype combining defective B and T cell responses[30]. Immunodeficiency plays a crucial role in the progression of many cancers, and immunodeficient patients with HCC exhibit higher morbidity and mortality[31,32]. With receptor binding and chemokine activity, CXCL1 is a member of the CXC subfamily of chemokines and plays multifarious roles during HCC development, metastasis, and prognosis[33,34]. ESR1, as one of the potential tumor suppressor genes, is a ligand-activated transcription factor closely related with cancer progression. Promoter hypermethylation is a possible mechanism by which ESR1 is silenced in human HCC[35]. MMP2 is among the most well characterized MMPs, and this factor plays essential roles in degrading ECM components. Depending on activation by the PI3K signalling pathway or the ERK and JNK pathways, it has been shown to have critical functions in the progression, invasion, and metastasis of HCC[36,37].

Regarding Hypo-HGs, we identified six hub genes: CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10. CDC45 plays an essential role in the initiation of DNA replication, which is consistent with the biological progress of Hypo-HGs GO analysis. It may be related with the progression of cancers due to induced DNA damage mediating regulation of DNA replication. DTL, also named CRL4 (CDT2), is a ubiquitin-protein ligase complex that plays important roles in the cell cycle, DNA synthesis, and DNA damage[38]. Previous studies disclosed that DTL had an oncogenic function in cancers, including HCC[39,40]. AURKB is a member of the Aurora kinase subfamily of conserved Serine/Threonine kinases with higher expression in tumor cells than normal cells, and its overexpression has been associated with biological characteristics of cancer as well as diagnosis[41-43]. A member of a protein phosphatase family with dual function in cell cycling, aberrant expression of CDKN3 is associated with carcinogenesis in many cancers, including HCC[44,45]. MCM2 is a highly conserved and essential mini-chromosome maintenance protein involved in the initiation of DNA and eukaryotic genome replication. The expression level of MCM has been associated with outcomes in many cancers and is closely related to HCC recurrence[46]. MCM10 plays a key role in cell cycle progression by mediating DNA replication initiation and elongation as well as preventing DNA damage and protecting genome integrity[47]. Evidence suggests that MCM10 is associated with inherited diseases resulting from genome instability and abnormal proliferation, and the level of MCM10 expression has been correlated with cancer progression and aggressiveness[48]. These findings indicate that the MDEGs in HCC may have a regulatory function in these biological processes and molecular function, and they are reliable with functional enrichment analysis. However, as some genes and pathways identified in the present study have not been formally investigated as targets in the progression of HCC, further research is needed.

Module analysis of the PPI network for Hyper-LGs suggested that the neuroactive ligand and ECM-receptor interaction, axon guidance, and chemokine signalling pathway might be involved in HCC progression. ECM-receptor interaction and axon guidance are critical cellular processes during the development of cancer. In addition, we found the neuroactive ligand-receptor interaction pathway to be related to hypermethylation, potentially resulting in abnormal expression of genes in cancers; more studies are necessary to validate these findings. Module analysis of the PPI network for Hypo-HGs showed complex roles for cell cycle and chemokine signalling pathways during HCC development. The cell cycle is a vital process involving DNA replication and translation, with a tendency for dysregulation in cancer[22]. Interestingly, the chemokine signalling pathway, as an essential process, was disclosed in both Hyper-LGs and Hypo-HGs modules, and this pathway has been shown to influence pathogenesis and metastasis of HCC by altering the tumor microenvironment[20].

In the present study, several limitations should be mentioned. First, the study lacked further experimental verification of the effects of aberrant methylation on gene expression and functions in HCC. Second, we did not investigate clinical parameters and prognosis, owing to the accessibility of data by bioinformatics arrays and tools. Third, as only two microarray profiles were analyzed, the sample size was not sufficiently large; thus, large-sample studies are required to validate the findings. In addition, HCC is closely related to hepatitis B and C, chronic alcoholism, tobacco smoking, and aflatoxins, and etiological factors were not analyzed in our study. Therefore, supplementary molecular experiments should be encouraged to verify further the results of our investigation.

In conclusion, using a series of bioinformatics databases and tools, we found that interactions among MDEGs of different functions and signalling pathways are related to the pathogenesis of HCC. Hub genes for Hyper-LGs of HCC included PTGS2, PIK3CD, CXCL1, ESR1, and MMP2; such genes for Hypo-HGs included CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10. As special biomarkers based on aberrant methylation, these hub genes might be useful for accurate diagnosis and treatment of HCC. This study provides hypothetical and biological characteristic insight into the pathogenesis of HCC. Additional molecular-level studies are needed to confirm the identified genes and pathways in HCC and to elucidate potential mechanisms.

ARTICLE HIGHLIGHTS

Research background

Pathogenesis of hepatocellular carcinoma (HCC) is a complicated biological process involving epigenetic and genetic changes. Most prior studies, however, mainly focused on either gene expression or methylation data but not the association and did not perform a conjoint analysis. The detection of methylated-differentially expressed genes (MDEGs) and a better understanding of their characteristics may be useful for discovering the molecular mechanism and pathogenesis of HCC.

Research motivation

In view of the insights from previous studies that MDEGs can be detected concurrently by joining gene expression and methylation microarray data, we explored the interaction network of differentially expressed genes and differentially methylated genes along with interrelated signalling pathways to find novel insights into the biological characteristics and pathways of methylated-differentially expressed genes in HCC.

Research objectives

The objective was to discover MDEGs in HCC, and explore relevant hub genes and potential pathways to make notional viewpoints available for the development and progression of HCC.

Research methods

We analyzed differentially methylated genes and differentially expressed genes using a series of bioinformatics databases and tools including GEO Datasets, DAVID, STRING, and Cytoscape.

Research results

We categorized 266 hypermethylated, lowly expressed genes (Hyper-LGs) and 161 hypomethylated, highly expressed genes (Hypo-HGs) in GO, KEGG, and PPI analyses. Hyper-LGs mainly refer to endogenous and hormone stimulus, cell surface receptor linked signal transduction, and behavior, while Hypo-HGs refer to DNA replication, metabolic processes, cell cycle, and cell division. Pathway analysis showed that Hyper-LGs were enriched in cancer, Wnt, and chemokine signalling pathways, while Hypo-HGs were related to cell cycle and steroid hormone biosynthesis pathways. Based on PPI networks, PTGS2, PIK3CD, CXCL1, ESR1, and MMP2 were identified as hub genes for Hyper-LGs, and CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10 were identified for Hypo-HGs by combining six ranked methods of cytoHubba.

Research conclusions

We found that interactions among MDEGs of different functions and signalling pathways are related to the pathogenesis of HCC by a series of bioinformatics databases and tools. Hub genes for Hyper-LGs of HCC included PTGS2, PIK3CD, CXCL1, ESR1, and MMP2; such genes for Hypo-HGs included CDC45, DTL, AURKB, CDKN3, MCM2, and MCM10. As special biomarkers based on aberrant methylation, these hub genes might be useful for accurate diagnosis and treatment of HCC. This study provides hypothetical and biological characteristic insight into the pathogenesis of HCC.

Research perspectives

The present findings indicate that the MDEGs in HCC can have a regulatory function in biological processes and molecular function and that they are reliable with functional enrichment analysis. As some genes and pathways identified in the present study have not been formally investigated as targets in the progression of HCC, further research is needed.