Using Integrated Bioinformatics Analysis to Identify Abnormally Methylated Differentially Expressed Genes in Hepatocellular Carcinoma

Background

As the second main cause of death for both men and women, cancer also ranks first or second in causing deaths in each age group among women.¹ Hepatocellular carcinoma (HCC), the major form of the primary liver cancer, has been the second principal cause of deaths related to cancer, with annual cases of 800,000 globally,² and involves many risk factors. Substantial studies have indicated that the risk factors influencing HCC comprise liver cirrhosis, hepatitis B virus (HBV), hepatitis C virus (HCV), diabetes, aflatoxins, alcohol, non-alcoholic fatty liver disease (NAFLD), and certain immune-related factors.^3–5 HCC is treated differently depending on the tumor stage, as well as liver function and performance. Early-stage radiofrequency ablation (RFA), liver transplantation (LT), and surgical resection are the three common treatments for early liver cancer, which are considered potentially curative.^6,7 In addition, primary liver tumors have high metastatic potential and are chemotherapy-resistant, so the HCC prognosis is inaccurate, and the postoperative recurrence is high.⁸ According to the Global Cancer Statistics in 2018, liver cancer is predicted to be the fourth leading cause of cancer death, and approximately 782,000 patients die.⁹ Therefore, screening for early HCC and treatment for advanced HCC are crucial for reducing the mortality of HCC.

Mechanisms of epigenetic regulation are critical for coordinating cellular processes as well as functions, including methylation of deoxyribonucleic acid (DNA), non-coding ribonucleic acid (RNA), and histone acetylation.¹⁰ As a key factor influencing the epigenetic modification, DNA methylation plays a vital role in the regulation of gene expression as well as strengthening its connection with chromatin remodeling and histone modification.¹¹ At the precancerous stage, the epigenetic abnormalities accumulate in the liver and thereby lead to tumorigenesis by altering hepatocyte differentiation and survival.¹² Many studies have reported changes in DNA methylation specific in HCC, especially at the complete genome level.^12–14 At present, changes in DNA methylation can be easily detected quantitatively from fixed tissues; thus, using the DNA methylation spectrum as an alternative biomarker is still an active field in clinical cancer research. At the same time, the detection of abnormally methylated expression genes and the in-depth understanding of their characteristics may help us understand the molecular mechanism of HCC, thus, providing new ideas for the liver cancer treatment.

With the large-scale application of sequencing technology in recent decades and the reduction of its cost, more studies on genomic liver cancer have provided new ideas for the diagnosis and treatment of hepatic carcinoma.^15,16 Under this context, microarrays based on the high-throughput platforms have become a potential new tool to identify target genes and epigenetic variations, offering new approaches for the HCC treatment.¹⁷ The data for this study were obtained from the GEO database, an international public repository, freely available, with a wide variety of microarray data from different researchers, including different species regions.¹⁸ In this study, we used bioinformatics methods to synthesize data from the gene expression profiling microarrays (GSE57957, GSE76427) and gene methylation-spectral microarrays (GSE89852, GSE54503). Target genes were identified using R software to explore the interaction networks and signaling pathways associated with liver cancer. Our aim is to find new abnormally methylated genes that can serve as biomarkers for an early liver cancer diagnosis in the future. The in-depth study of key signaling pathways can provide innovative ideas for the precise treatment of liver tumors.

Methods

Data Information of Microarray

We obtained two gene expression profiling datasets (GSE57957, GSE76427) and two gene methylation profiling datasets (GSE89852, GSE54503) through GEO (https://www.ncbi.nlm.nih.gov/geo/). The statement to confirm that all methods were carried out in accordance with relevant guidelines and regulations, and that all experimental protocols were approved by a named institutional and/or licensing committee, and that informed consent was obtained from all subjects or, if subjects are under 18, from a parent and/or legal guardian can be obtained from GEO (SingHealth Centralized Institutional Review Board, PMID: 25093504; SingHealth Institutional Review Board, PMID: 29117471; the Ethics Committees of the National Cancer Center, the National Center for Global Health and Medicine, and Yotsuya Medical Cube, PMID: 28426876; the Institutional Review Board of Columbia University Medical Center (CUMC), PMID: 25861255). The GPL10558 platform (Illumina Human HT-12 V4.0 expression BeadChip) was used to register GSE57957 and GSE76427. GSE57957 obtained 39 tumor samples and 39 adjacent non-tumor samples; whereas, GSE76427 included 115 primary tissue samples and 52 adjacent non-tumor tissue samples derived from 115 HCC patients. GSE89852 and GSE54503 were registered on the GPL13534 platform (Illumina Human Methylation 450 BeadChip). GSE89852 included 37 paired normal and HCC samples from 37 patients; whereas, the samples of 66 HCC tumor and 66 adjacent non-tumor were obtained in GSE54503.

Data Acquisition and Processing

The methylation and expression profiling microarray datasets were analyzed to identify differential expression in the microarray data by using the limma package (v3.12) in R software (v3.6.1) (https://www.r-project.org/). The cut-off standard of P < 0.05 was set to analyze the differentially expressed genes (DEGs). The cut-off values of FDR < 0.05 and β > 0.2 were used to analyze the differentially methylated genes (DMGs). Finally, for the differential expression analysis results, the Venn Diagram package (v1.6.20) in R software was used to identify the overlapping DMGs and DEGs. Overlapping hypermethylated and down-regulated genes generated a Venn diagram of genes with hypermethylation and low expression. Overlapping hypomethylated and up-regulated genes generated a Venn diagram of hypomethylated and highly expressed genes.

Analysis of Functional and Pathway Enrichment

The functional profiles for genes and gene were analyzed and visualized by using the clusters R package software of clusterProfiler (v3.0.4). We investigated the functional profiles of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene ontology (GO). We set the P value Cutoff to 0.05 and the Q value Cutoff to 0.05.

Generation of PPI Network and Module Analysis

STRING database is a search tool for gene and protein interactions, which helps users easily obtain unique, wide-ranging experiments, and predictive interaction information, and it was used in this study to construct the PPI network of differential genes. The interaction score ≥0.4 was set as a threshold. We then imported the results of the STRING database into Cytoscape software (v3.7.2) and used cytoHubba (v3.18.0) to view the hub genes. To select the six hub genes, we used five ranking methods and overlapped the top ten genes. Betweenness, Closeness, Degree, EcCentricity and Edge percolated component (EPC) methods were used to analyze hypermethylation and down expression genes, while, Betweenness, Bottleneck, Closeness, Radiality and Stress methods were used for hypermethylation and down expression genes. In the Cytoscape software, the MCODE plugin was utilized for module analysis, and the standard values for defining modules were used (MCODE score >3; the number of nodes >5).

Hub Gene Verification

The online database GEPIA (http://gepia.cancer-pku.cn/) is used to interactively analyze gene expression profiles founded on the TCGA database. Customizable features, such as normal/tumor differential expression analysis, tumor type, pathology staging, survival analysis, dimensional analysis, correlation analysis, and similar gene detection, are provided by GEPIA. In this study, we identified differential mRNA expression levels of the hub genes between the liver cancer tumor and non-tumor samples using GEPIA. TCGA comprises comprehensive data of cancer genomic alterations, including mutation, copy number variation, mRNA expression, and methylation data. Besides, the HPA database was employed to validate the translational levels of the hub genes.

Survival Analysis and Genetic Alteration of Hub Genes

The multidimensional cancer genomics are explored, visualized, and analyzed through a resource network provided by cBioPortal (http://www.cbioportal.org/). In this database, large-scale cancer genomic datasets for a variety of cancers are allowed to be visualized, analyzed, and downloaded. In this study, the cBioPortal was used to explore the genetic alterations of the hub genes and the correlation between DNA methylation status and mRNA expression. The GEPIA online tool was employed to assess the prognostic value of the hub genes. We set P < 0.05 as a threshold to indicate the significant difference.

Results

Identification of Abnormally MDEGs in HCC

Figure 1 shows a study design flowchart. By R software, we separately analyzed the data from each microarray to acquire the DEGs or DMGs. In two datasets, 4264 genes were up-regulated (5379 in GSE57957; 6826 in GSE76427) while 2888 were down-regulated (4213 in GSE57957; 6873 in GSE76427). For microarrays of gene methylation, 1780 overlapping hypermethylated genes (7391 in GSE89852; 5116 in GSE54503) together with 4550 overlapping hypomethylated genes (41,046 in GSE89852; 23,579 in GSE54503) were identified. By comparing these two sets of genes, a number of 162 hypermethylated, low-expressing genes and 190 hypomethylated, highly expressed genes were identified (Figure 2). The clustered heat map of GSE57957 and GSE76427 (with the top 50 DEGs) and the top 50 DMGs of GSE89852 and GSE54503 are shown in Figure 3.

Figure 1 The flow chart of this study. Abbreviations: DEGs, differentially expressed genes; DMGs, differentially methylated genes; Hyper-LGs, hypermethylation and low expression genes; Hpo-HGs, hypomethylation and high expression genes.

Figure 2 Identification of abnormally methylated-differentially expressed genes (MDEGs) in mRNA expression profiling datasets (GSE57957, GSE76427) and gene methylation profiling datasets (GSE89852, GSE54503). (A) Hypermethylation and low expression genes. (B) Hypomethylation and high expression genes.

Figure 3 Clustered heat map of the top 50 differentially expressed genes (DEGs) and differentially methylated genes (DMGs). Red: upregulated genes; Blue: downregulated genes. (A and B) The heat map of top 50 DEGs in GSE57959 and GSE76427; (C and D) The heat map of top 50 DMGs in GSE89852 and GSE54503.

Analysis of GO Functional and KEGG Pathway Enrichment

We used the clusterProfiler R package to perform GO functional annotation and KEGG pathway enrichment analysis. The results of the significant enrichment analyses are shown in Table 1. For hypomethylated and highly expressed genes, biological process (BP) analysis results were mainly enriched in calcium homeostasis, keratinocyte proliferation, and collagen-activation signaling pathways. The analysis results of Molecular function (MF) recommended that these genes are mainly involved in growth factor binding, PDZ domain binding, insulin-like growth factor binding, histone serine kinase activity, and adenosine monophosphate (AMP)-activated protein kinase activity. In addition, cellular component (CC) analysis primarily indicated cell-matrix junctions, focal adhesions, cell-matrix adhesion junctions, and membrane anchoring components. For the hypermethylated and low-expressing genes, the BP analysis was mainly enriched in peptide hormone responses, multicellular biological processes, female pregnancy, and regulation of estradiol. MF enrichment indicated proximal promoter sequence-specific DNA binding, organic anion transmembrane transporter, RNA polymerase II proximal promoter sequence-specific DNA binding, and peptidoglycan receptor activity. In addition, CC analysis results were mainly enriched in the basolateral plasma membrane, multivesicular body, filamentous actin, and astrocyte projection.

Table 1 Gene Ontology Analysis of Aberrantly Methylated Differentially Expressed Genes in Hepatocellular Carcinoma

Analysis of KEGG Pathway

It can be seen from the KEGG pathway analysis that the hypomethylated, highly expressed genes were remarkably enriched in extracellular matrix (ECM)-receptor interaction, focal adhesion, and phosphatidylinositol 3-kinase (PI3K)-Akt signaling pathways. Table 2 shows that Hypermethylated and low-expressing genes were enriched in axon guidance, cyclic adenosine monophosphate (cAMP) signaling, and Wnt signaling pathway.

Table 2 KEGG Pathway Analysis of Aberrantly Methylated Differentially Expressed Genes in Hepatocellular Carcinoma

Construction of PPI Network and cytoHubba Analysis

We constructed 80 nodes and 106 edges from the hypomethylated, highly expressed genes, and 108 nodes and 166 edges from the hypermethylated, low-expressing genes. The PPI networks for hypomethylated, highly expressed genes and hypermethylated, low-expressing genes are shown in Figures 4 and 5, respectively. The hub genes were selected by using cytoHubba. By overlapping the top ten genes based on five ranking methods in cytoHubba, three down-regulated, hyper methylated genes (CFTR, ESR1, and CXCL12) and three up-regulated, hypomethylated genes (PTK2, VWF, and ITGA2) were identified (Table 3).

Figure 4 Protein–protein interaction (PPI) network of the hypomethylation and high expression genes.

Figure 5 Protein–protein interaction (PPI) network of the hypermethylation and low expression genes.

Table 3 Hub Genes Among the Aberrantly Methylated Differentially Expressed Genes Ranked in cytoHubba

Module Analysis

Figure 6 shows that two modules in the hypomethylated, highly expressed genes network, and one module in the hypermethylated, low-expressing genes network were statistically significant. In Table 4, we analyzed the GO terms and KEGG pathways. The results indicated that the hypomethylated, highly expressed gene was enriched in pathways significantly concerning ECM-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway. Hypermethylated, low-expressing genes showed enrichment in Estrogen signaling pathway and cytokine-cytokine receptor interaction.

Figure 6 Module analysis of abnormally methylated differentially expressed genes (MDEGs). (A) The hypomethylation and high expression genes in Module 1; (B) the hypomethylation and high expression genes in Module 2; (C) the hypermethylation and low expression genes in Module 1.

Table 4 Module Analysis of the Protein–Protein Interaction Network

Identification and Validation of Six Selected Genes

The GEPIA database was employed to investigate the six-hub gene expression, contributing to the further validation of the results. We found that the expression of the three up-regulated, hypomethylated genes and the three down-regulated, hypermethylated genes were dramatically different between tumor and normal tissues (Figure 7). However, the expression of CFTR and ITGA2 was not significantly different between tumor and normal samples. Further experiments are required to confirm this result. Nevertheless, data for the rest four genes were in agreement with our data. Next, TCGA was employed to confirm the results, and the outcome is shown in Table 5. It is found that agreed with our findings, the methylation and expression statuses were altered in TGCA data significantly. Furthermore, the HPA database was used to detect the hub gene expression at the translation level. For up-regulated, hypomethylated genes, ITGA2 and PTK2 were expressed higher in the liver cancer tissue (Figure 8). The difference in VWF between normal liver and cancer samples was not distinct. For down-regulated, hypermethylated genes, the lower expression in the liver cancer was demonstrated through the immuno-histochemical ESR1 results. The differences between ESR1 and CFTR in normal liver and cancer samples were not apparent. To determine if the identified genes are translationally relevant to the liver cancer, further verifications from more clinical samples are required.

Table 5 Validation of the Hub Genes in TCGA

Figure 7 The expression of the six genes in mRNA expression, using data from the TCGA database in GEPIA.

Figure 8 Validation of the six hub genes on a translational level using the Human Protein Atlas database.

Genetic Information of the Six Genes

The connection between six genes and prognosis was evaluated by conducting a survival analysis with GEPIA. We discovered that ITGA2 and ESR1 were tightly related to the overall survival (OS; Figure 9A). For the up-regulated, hypomethylated oncogene ITGA2, its upregulation led to shorter OS. However, downregulated expression of the down-regulated, hypermethylated oncogene ESR1 resulted in worse OS in HCC patients. As for the rest four genes, the trends were close to our prognosis, but not statistically significant. By using the cBioPortal tool, a hub gene network was constructed to find the corresponding drug and the potential therapeutic target. A network made up of six hub genes and 50 most frequently altered neighboring genes of them are illustrated in Figure 9B. The drug targets network showed that five genes, ESR1, CFTR, CXCL12, VWF, and PTK2, are the most popular drug targets nowadays (Figure 9C). We surmise that in the future, ITGA2 could find application as a new drug target. Figure 10A and B exhibits the information on the hub gene alteration. It can be seen that the six hub genes were varied in 277 (51%) in the 442 sequenced cases/patients. Meanwhile, PTK2 and ITGA were altered most often (46% and 8%), including deep deletion, amplification, and missense. The correlation between mRNA expression and DNA methylation of the hub genes in the HCC patient dataset of TCGA is shown in Figure 10C. A negative correlation, which indicates that the methylation regulates the genes’ mRNA expression (except for ESR1), was found and demonstrated the important role of methylation in expressing these genes.

Figure 9 Overall survival analyses of ESR1 and ITGA2 in HCC patient. (A) Downregulation of ESR1 was closely associated with worse overall survival; (B) downregulation of ITGA2 resulted in longer overall survival; (C) the network contained our six hub genes and the 50 most frequently altered neighbor genes. The relationship between hub genes and drugs is also exhibited.

Figure 10 Genetic alterations connected with the six genes and the relationship between mRNA expression and DNA methylation in the TCGA HCC study. (A) A visual summary of alteration based on a query of the six hub genes, which were altered in 277 (51%) in the 442 sequenced cases/patients; (B) alteration frequency of hub genes; (C) the correlation between mRNA expression and DNA methylation of the six hub genes in the HCC patient dataset of TCGA.

Discussion

In the development of HCC, abnormal methylation of genes related to tumor not only occurs in the terminal tumor, but is also a common early event. In addition, the frequency of promoter abnormal methylation develops from a precancerous lesion to HCC.¹⁹ Related studies have demonstrated that epigenetically silenced tumor suppressor genes can be re-expressed by the use of demethylation and histone modifiers.^20,21 It could be a routine treatment similar to other malignancies in the coming years. Therefore, epigenetic changes in tumors can not only serve as indicators or “biomarkers” for screening patients at high stake for HCC, but also provide innovative ideas for the treatment of HCC. In our research, we identified three hypomethylated, highly expressed genes, and three hypermethylated, low-expressing genes using bioinformatics analysis. We clearly verified the functional and enrichment analysis of the hub genes. The interaction network results revealed that, with the development of HCC, relevant genes may involve important pathways related to molecular regulation.

The GO enrichment analysis showed that the BP of genes with hypomethylation and high expression mainly focused on calcium ion homeostasis and keratinocyte proliferation. It is reported that intracellular calcium homeostasis can be adjusted from autophagy, which is closely linked to the emergence and development of tumor formation and autophagy.^22,23 Among the genes with high methylation and low expression, BP mainly concentrates in the response to peptide hormones, multi-multicellular organism processes, and female pregnancy. Previous studies have shown that atrial natriuretic peptides (ANP) can eliminate pancreatic adenocarcinoma in up to 80% of thymus-free mice.²⁴ Vasoactive intestinal peptide (VIP), as a regulator of the inflammatory response, is expressed in many types of tumors. Studies have shown that VIP can suppress the apoptosis of HCC cells through the cAMP/Bcl-xl pathway to prevent the progression of HCC.²⁵ In terms of CC, genes with hypomethylation and high expression are mainly concerned with cell-substrate junctions, focal adhesions, and cell-substrate adherences junctions; whereas, genes with hypermethylation and low expression are mainly concerned with the composition of the multivesicular body, basolateral plasma membrane, and filamentous actin. The adhesion and capture of tumor cells in the microcirculation are the primary conditions for tumor vascular metastasis. Therefore, overexpression of cell adhesion molecules can enhance the interaction between focal adhesion and ECM receptors, and accelerate the metastasis of liver cancer cells in the blood vessels and settlement at the metastasis site.^26–28 Leda-1/Pianp is a transmembrane protein of type I located in the basolateral membrane of polarized epithelial cells and enhances the tumor cells’ resistance to unstable junctions and the tightness and stability of intercellular junctions.²⁹ In terms of MF, the main function of genes with hypomethylation and high expression is binding to growth factors, PDZ domain binding, and insulin-like growth factor binding. However, the MF of genes with hypermethylation and low expression is mainly for RNA polymerase II proximal promoter sequence-specific DNA binding and organic anion transmembrane transporter activity. For example, the insulin-like growth factor-I receptor (IGF-IR) in most of the transformed cells exhibits potential cell-survival and antiapoptotic activities. Additionally, its expression is negatively regulated by p53, BRCA1, and WT1.³⁰ Studies based on animal models, as well as in vitro, have demonstrated that IGF plays a key role in the cell cycle, proliferation, survival, migration, apoptosis inhibition, protein synthesis, and cell growth of HCC.³¹

From the KEGG pathway enrichment analysis, hypomethylated, highly expressed genes were found to be significantly enriched in HPV infection, ECM-receptor interactions, focal adhesions, the PI3K-Akt signaling pathway, and the AMPK signaling pathway. Studies show that by inhibiting the PI3K/Akt pathway, petroleum ether (PE) may induce the apoptosis of hepatocellular carcinoma cells, thus producing an anti-tumor effect.³² Hou et al indicated that inhibition of the PI3K/Akt pathway could significantly reduce the phosphorylation level of FOXO3a, accompanied by an increase in the total level of FOXO3a, and significantly inhibited the proliferation of tumor cells.³³ Stearoyl-CoA desaturase 1 (SCD1) regulates autophagy through the AMPK pathway and has a great influence on HCC.³⁴ At the same time, we found that interaction of neuroactive ligand-receptor is the most abundant enrichment pathway with low expression of hypermethylated genes. In addition, concentration is also found in the cAMP pathway, oxytocin signaling pathway, axon guidance, and Wnt signaling pathway. Liu et al proposed a Bayesian method to demonstrate that mitogen-activated protein kinase (MAPK) signaling pathway and neuroactive ligand–receptor interaction are closely related to HCC.³⁵ Tao et al found that the upstream KEGG pathway of methylated genes in individual hepatocytes in HBV-associated hepatocellular carcinoma is a “neuroactive ligand–receptor interaction”.³⁶ These results show that the neuroactive ligand-receptor pathway greatly affects the pathogenesis of liver cancer.

The results based on PPI analysis showed that three hypomethylated and highly expressed hub genes were PTK2, VWF, and ITGA2. PTK2 activates Wnt/β-catenin signaling by fostering the beta-catenin nuclear accumulation in HCC cells, which is connected with OS, marker expression of liver cancer stem cell, and recurrence-free survival (RFS) in HCC patients.³⁷ In addition, increased PTK2 expression can improve the prognosis of chronic lymphocytic leukemia (CLL) patients treated with cyclophosphamide and rituximab-cyclophosphamide. PTK2 expression may be a useful biomarker in future clinical trials.³⁸ In our study, PTK2 expression in tumor samples was higher than that in normal tissues, and there was a statistically dramatic difference between these two groups. At the same time, our data showed a negative correlation between PTK2 mRNA expression and DNA methylation. The PTK2 mutation rate was found to be up to 46%, suggesting that abnormal methylation of DNA may result in mutations and can potentially be a target for future drug therapy. As the largest plasma protein in humans, VWF is an adhesive polysaccharide protein that intervenes the adhesion of platelets to the endothelial surface and subendothelial matrix. It is also the carrier of clotting factor VIII in the blood.³⁹ A review of elevated VWF levels concerned with increased risk of myocardial infarction, stroke, venous thrombosis, and diabetes suggests that the VWF gene plays a crucial role in the regulation of the blood clotting and bleeding balance.⁴⁰ In addition to its important role in hemostasis, an increasing amount of literature has shown that VWF has additional anti-tumor effects, mainly through negative regulation of angiogenesis and apoptosis.⁴¹ Liu et al demonstrated that VWF is connected to the pathogenesis of HBV infection and duplication, and is also related to the clinicopathological staging of HBV-infected HCC patients.⁴² The results showed that the VWF expression in normal and tumor tissues was dramatically different, and that the VWF expression in liver cancer samples was apparently higher than that in normal tissues. However, data obtained from HPA showed no obvious difference in VWF gene expression between normal liver tissues and tumor tissues. The mutation rate of VWF was up to 7%, and the expression of mRNA was negatively correlated with methylation. This may be the result of abnormal methylation, which needs to be confirmed by further experiments. It has been reported that ITGA2 is essential in cell survival, invasion, migration, and angiogenesis.^43,44 Overexpression of ITGA2, which is methylation-independent, is connected with inaccurate prognosis in acute myeloid leukemia (AML).⁴⁵ In addition, ADAR1 can enhance HCC metastasis by increasing ITGA2 expression and promoting tumor cell adhesion to the ECM.⁴⁴ In our study, we found that highly expressed ITGA2 is closely related to OS and RFS. In the verification of HPA data, the expression of ITGA2 in liver cancer tissues was obviously higher than that in normal samples. The mutation rate of ITGA2 was 8%, suggesting that an analogous mechanism might exist in HCC tumorigenesis. At the same time, we found that the expression of ITGA2 was negatively correlated with the methylation state. Therefore, we can speculate that ITGA2 gene mutation may lead to its hypomethylation; thus, contributing to the emergence and development of liver cancer. The three hypomethylated and highly expressed genes screened in this study were all highly expressed in the liver cancer tissues, and the mRNA expression levels were negatively correlated with the methylation state. This suggests that abnormal methylation of these three genes may accelerate the development of liver cancer.

With regard to the hypermethylated and low-expressing genes, three genes were identified: CFTR, ESR1, and CXCL12. CFTR is a large glycoprotein located on chromosome 7, belonging to the adenosine triphosphate (ATP)-binding superfamily of proteins. CFTR is expressed in various cell types and is concerned with epithelial cells in the airway, sinuses, gastrointestinal tract, sweat glands, and genital and nervous systems.⁴⁶ At present, there are no reports on the impact of CFTR on the occurrence and development of liver cancer. Further experiments are needed to confirm its function in HCC. In both males and females, the estrogen receptor (ER), ESR1, is the main receptor in the liver and mediates the liver’s response to estrogen.⁴⁷ As an important sex hormone, estrogen functions in various processes such as cell metabolism, cell differentiation, and tissue development.⁴⁸ Villa et al proved that the transcript status of ER in patients with HCC has a novel value in prognosis.⁴⁹ Moreover, in inoperable patients with HCC, ER transcription was the largest negative predictor of survival. Naugler et al found a protective effect of ER alpha in liver cancer.⁵⁰ In males, Kupffer cells (KCs) can produce estrogen-mediated inhibitions of interleukin (IL)-6 to reduce the risk of liver cancer, and these outcomes may be useful for the prevention of HCC. In this research, we discovered that the ESR1 expression in tumor tissues was lower than that in normal tissues, and the difference between the two tissues was statistically significant. This finding was consistent with the HPA validation data. Highly expressed ERS1 is associated with better OS and RFS. The mutation rate of the ESR1 gene was 6%, and the mRNA expression was positively correlated with the methylation state. This conclusion requires further experimental confirmation. Chemokines are potential targets for anti-tumor angiogenesis therapy, and the chemokine, CXCL12, promotes tumor metastasis and angiogenesis mainly through the receptor, CXCR4.⁵¹ A meta-analysis investigating the relationship between CXCL12 expression and various cancer survivals indicated that CXCL12 plays a crucial role in cancer occurrence.⁵² However, its role in liver cancer needs to be further verified. From our results, we found that the expression of CXCL12 was indeed low in HCC samples from multiple datasets. The gene mutation rate of CXCL12 was 5%, and its expression was negatively correlated with the methylation state. We speculate that the mutation of CXCL12 may lead to hypermethylation of the gene and decreased mRNA expression, thus, contributing to the emergence and development of liver cancer. Compared with non-tumor tissues, the three differentially expressed genes screened in this study were down-regulated in liver cancer; however, in liver cancer, their action mechanism of remains unclear. Therefore, we surmise that the abnormal methylation of these three genes may promote the development of liver cancer through other aspects besides expression.

Module analysis of PPI network for hypermethylated and low-expressing genes indicated that the estrogen signaling pathway, cAMP signaling pathway, and interaction of neuroactive ligand-receptor, might be involved in the HCC progression. As discussed above, the neuroactive ligand-receptor interaction pathway plays a critical role in liver cancer, which is highly consistent with our conclusion. Although the liver is not a classic target organ for estrogen, studies have shown that ER exhibits a rapid hepatocyte response to estrogen, suggesting that estrogen signaling contributes to rapid occurrence and progression of liver cancer. Meanwhile, studies have shown that estrogen can regulate the HCC malignant degree in vivo by depressing the invasion of tumor cells, impeding cell cycle progression, and facilitating apoptosis.⁵³ Cyclic Adenosine monophosphate (cAMP) is a key intracellular second messenger that can affect plenty of cellular processes.⁵⁴ Hara et al pointed out that VIP can impede the HCC progression through apoptosis via the cAMP/Bcl-xL pathway.²⁵ More research may be needed to verify our results. The PPI network module analysis for hypomethylated and highly expressed genes demonstrated that the focal adhesion, ECM-receptor interaction, human papillomavirus infection, and proteasome might be connected with HCC progression. The focal adhesion and ECM-receptor interaction are crucial cellular processes in the development of cancer. However, the roles of HPV infection and the proteasome in HCC, together with the impact of abnormal methylation, remain unknown. In HCC occurrence and progression, future investigations are needed to explain how abnormal methylation influences the functional roles of these pathways.

There are some limitations to be noted in our research. First, there is a lack of analysis of clinical data in our study. Secondly, there are inconsistencies in the conclusions obtained from data obtained from different platforms; and, further experiments are required to verify our conclusions. Thirdly, HCC is closely related to HBV, HBC, chronic alcoholism, smoking, aflatoxins, and other factors. No relevant data sets were selected for analysis in this study. In general, our conclusions are based on bioinformatics analysis, which has certain reference value; however, to further verify our results, complementary molecular experiments are encouraged.

Conclusions

In summary, new abnormally methylated genes and HCC-related pathways were discovered in this study. Moreover, we have verified the six selected genes with abnormal methylation from multiple aspects, and these outcomes could provide us with a greater insight into the development and progression of HCC. Based on abnormal methylation, hub genes, including PTK2, VWF, ITGA2, CFTR, ESR1, and CXCL12, could serve as biomarkers for an accurate HCC diagnosis and treatment. It is hoped that these findings can attract more attention to the molecular mechanisms of cancer in the future.