Integrative Analysis of DNA Methylation and Gene Expression Profiling Data Reveals Candidate Methylation-Regulated Genes in Hepatoblastoma

Introduction

Hepatoblastoma (HB), a type of liver cancer in children and infants, occurs frequently and accounts for approximately 50% of pediatric liver cancers.¹ The incidence of HB is approximately 1/1,000,000 in children aged < 15 years, and nearly 20% of them have synchronous metastasis at the first diagnosis.^2,3 Although traditional approaches, including chemotherapy, surgical resection, and transplantation, remain the most common and effective treatment options for HB and have greatly progressed in the past 20 years, the prognosis of patients with clinically advanced HB has not improved.^2,3 Currently, the lack of effective biomarkers for early HB diagnosis greatly limits the treatment of HB.⁴ Therefore, it is essential to screen effective biomarkers for the early diagnosis and treatment of HB.

The development of molecular biology has made it possible to diagnose tumors based on biomarkers, which are more specific and sensitive than traditional diagnostic methods.^5–7 DNA methylation is a major epigenetic modification that can affect the expression and alternative splicing of genes and maintain the stability of genomes.⁸ It has been reported that aberrant DNA methylation in various cancers plays a critical role in the occurrence and development of tumors by regulating the expression of many tumor-associated genes.⁹ For example, some studies have indicated that hypomethylation of oncogenes and hypermethylation of tumor suppressor genes are associated with the occurrence of hepatocellular carcinoma.^10–12 In addition, aberrant DNA methylation affects the biological processes of bladder cancer.^13,14 Furthermore, DNA methylation of genes can also be selected as a potential biomarker for the diagnosis of tumorigenesis and the prediction of the prognosis of tumors based on their effective diagnostic and prognostic value.^15,16 Dulaimi et al revealed that hypermethylation of APC, RASSF1A, and ARF genes can be used for the early diagnosis of bladder cancer.¹⁷ Moreover, Zhou et al suggested that CDKN2A methylation has an important impact on the occurrence, development, and metastasis of head and neck squamous cell carcinoma and may serve as a promising diagnostic or prognostic marker.¹⁸ Hence, DNA methylation is crucial for the diagnosis and prognosis of tumors and is conducive to improving the treatment strategies for tumors.

Previous extensive studies have revealed that abnormal DNA methylation is related to the prognosis of patients with HB. For instance, Honda et al suggested that increasing GPR180 methylation is associated with a poor prognosis of HB.¹⁹ Moreover, Sakamoto et al also found that MT1G hypermethylation can be selected as a prognostic indicator of HB, and patients with HB may benefit from demethylating drug treatments by targeting MT1G.²⁰ In particular, by integrating analyses of DNA methylation, mutation, and gene expression data, Sekiguchi et al revealed the heterogeneity of HB and found that NQO1 and ODC1 can act as novel druggable targets.²¹ Moreover, combining genomic, transcriptomic, and epigenomic data, Carrillo-Reixach et al also identified choline kinase alpha as a promising therapeutic target for patients with intermediate- and high-risk HB.²² Thus, integrative analysis of DNA methylation and gene expression data for patients with HB may help elucidate the molecular mechanisms and identify novel therapeutic targets for HB.

Our study aimed to identify DNA methylation-regulated genes as novel potential diagnostic biomarkers and treatment targets for patients with HB by integrating methylation and gene expression data from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database,²³ which may contribute to the improvement of the molecular diagnosis and clinical treatment of HB.

Materials and Methods

Collection and Processing of Data

DNA methylation data of 19 HB tumor samples and 10 normal liver samples, generated from Illumina Human Methylation 450K BeadChip was obtained from the GSE78732 dataset in the GEO database. The methylation level of each methylated position was normalized using the β-value, which ranged from 0 to 1. The formula for calculating the β-value is M/(M + U + 100), where M represents the methylated allele frequency, and U represents the unmethylated allele frequency. A high β-value implies a high methylation level.

Gene expression profiling data of 53 HB tumor samples and 14 normal liver samples sequenced by an array was obtained from the GSE131329 dataset in the GEO database, and the clinical information of the HB samples is summarized in Table 1. Moreover, gene expression profiling data of 31 HB tumor samples and 32 normal liver samples sequenced by high-throughput sequencing were obtained from the GSE133039 dataset in the GEO database and were used to verify the expression of underlying genes. Raw data was normalized using the fragments per kilobase per million value.

Table 1 The Clinical Information of the 53 HB Patients the GSE131329 Dataset

Identification of Differentially Expressed Genes (DEGs) and Differentially Methylated Genes (DMGs)

Initially, principal component analysis (PCA) was conducted to determine the sample distribution using the gmodels R package. DEGs between the HB and normal liver samples from the GSE131329 dataset were identified using the Limma R package at a cutoff value of |log2 (FC) | ≥ 1 and a P value of < 0.05. Differential DNA methylation sites (DMSs) with a threshold of |log2 (FC) | ≥ 1 and an adjusted P value of ≤ 0.05 and differential DNA methylation regions (DMRs) with a threshold Fisher value of ≤ 0.05 and a |mean difference| of > 0.1 between the HB and normal liver samples were analyzed in the GSE78732 dataset using the methylationArrayAnalysis R package. Genes located in the DMRs were defined as DMGs. The ggplot R package was used to plot volcano plots and heatmaps of the DEGs and DMGs.

Integrative Analysis of DNA Methylation and Gene Expression

Integrative analysis of the DMGs and DEGs was performed to further screen for the methylation-regulated genes in HB. A Venn diagram was used to identify hypomethylated/upregulated genes and hypermethylated/downregulated genes from the DMGs and DEGs, respectively.

Survival Analysis

To screen for potential key genes in HB, we analyzed the correlation between the expression levels of the methylation-regulated genes and survival, including the overall survival (OS) and disease-free survival (DFS) of the patients with liver cancer, using the OncoLnc database (http://www.oncolnc.org/); genes with a P value of < 0.05 were reserved for further analysis.

Disease Association Analysis

To further investigate the association between the potential key genes and HB, we performed disease association analysis using the VarElect online tool. By importing “hepatoblastoma,” “liver cancer,” and “cancer,” we obtained the roles of the key genes in HB, liver cancer, or cancer. Moreover, based on the score, we selected the top three genes as the key DNA methylation-regulated genes.

Construction of the Protein–Protein Interaction (PPI) Network

The interactions among the key DNA methylation-regulated genes and other DEGs were analyzed using the Search Tool for the Retrieval of Interacting Genes (https://string-db.org/). Moreover, we further investigated the biological processes of these genes using the GeneMANIA online database (http://genemania.org/).

Biological Function Analysis

To further investigate the underlying molecular mechanisms of the key DNA methylation-regulated genes, Gene Ontology annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to explore the potential biological functions of the DEGs using the clusterProfiler R package at a cutoff value of P < 0.05.

Methylation Levels of the Key DNA Methylation-Regulated Genes

To further investigate the methylation levels of the key DNA methylation-regulated genes, we first analyzed the methylation sites and methylation regions of these genes. Next, a Circos diagram was plotted to show the DNA methylation and expression of the key DNA methylation-regulated genes using Circos v.0.69-8.

Expression Levels and Diagnostic Significance of the Key DNA Methylation-Regulated Genes

To further validate the expression levels of the key DNA methylation-regulated genes in HB, we first analyzed their expression levels in the HB samples from the GSE131329 dataset. We also validated their expression levels in the HB samples from the GSE133039 dataset. Furthermore, receiver operating characteristic (ROC) curve analysis was conducted to explore whether the key DNA methylation-regulated genes could differentiate the HB samples from the normal liver samples using the pROC R package. Finally, we compared the expression levels of the key DNA methylation-regulated genes among different clinical features.

Statistical Analysis

Statistical analysis was performed using the R software. Differences between the two groups were analyzed using Student’s t-test. Differences were considered statistically significant at P < 0.05.

Results

Identification of the DEGs in HB

To observe the distribution of the samples in the GSE131329 dataset, PCA was performed using the gene expression data of each sample. PCA showed that the HB samples and normal liver samples were distributed on the two flanks of the horizontal axis (Figure 1A), indicating that the samples in each group had good repeatability and could be used for subsequent analysis. A total of 1597 DEGs were identified between the HB and normal liver samples, including 850 upregulated genes and 747 downregulated genes in the HB samples compared with those in the normal liver samples (Figure 1B and C).

Figure 1 Identification of DEGs between HB and normal liver samples in the GSE131329 dataset. (A) PCA showed the distribution of HB and normal liver samples based on overall expression values. (B) volcano plot of DEGs; Red represents upregulated genes and blue represents downregulated genes. (C) heatmap of DEGs; Blue represents downregulated genes, and red represents upregulated genes.

Analysis of DNA Methylation in HB

Similarly, PCA was performed to observe the distribution of the samples in the GSE78732 dataset using M values of each sample. Clearly, PCA showed that the HB and normal liver samples exhibited distinct clustering (Figure 2A). Further, multidimensional scaling was conducted to evaluate the distribution of the M values in each sample. We found that the M values varied from −5 (hypomethylation) to 5 (hypermethylation) in the HB and normal liver samples (Figure 2B). A total of 24,460 DMSs were identified, including 4209 hypermethylated and 20,251 hypomethylated sites in the HB samples compared with those in the normal liver samples, at a cutoff value of |log2(FC)| ≥ 1 and a P value of ≤ 0.05 (Figure 2C and D). Next, we identified 451 DMRs, including 132 hypermethylated and 319 hypomethylated regions, at a threshold Fisher value of ≤ 0.05 and a |mean difference| of > 0.1. For the DMRs, we found 138 DMGs in the hypomethylated regions and 319 DMGs in the hypermethylated regions.

Figure 2 Analysis of DNA methylation for HB samples in the GSE78732 dataset. (A) PCA showed the distribution of HB and normal liver samples based on methylation M values. (B) The distribution of methylation density; The horizontal axis represents the methylation M value, and the vertical axis represents the methylation distribution density of gene. (C) volcano plot of DMSs between HB samples and normal liver samples; Red represents hypermethylation sites and blue represents demethylation sites. (D) Manhattan plot showed the distribution of DMSs in each chromosome; The vertical axis represents -log10(adj.P.Val), the horizontal axis represents chromosomes, each dot represents a methylation site, and different colors have no practical meaning.

Integrative Analysis of DNA Methylation and Gene Expression in HB

To identify the genes potentially regulated by DNA methylation, we overlapped the DMGs and DEGs and found 22 hypomethylated and upregulated genes and 19 hypermethylated and downregulated genes (Figure 3A and B). Thus, these genes were defined as the DNA methylation-regulated genes.

Figure 3 Identification of DNA methylation-regulated genes in HB. (A) hypomethylated and upregulated genes in HB samples compared with normal liver samples. (B) hypermethylated and downregulated genes in HB samples compared with normal samples.

Survival Analysis

To screen for the potential key genes from the DNA methylation-regulated genes in HB, we further examined the associations between these genes and the survival of the patients with liver cancer. We found that 13 of 41 DNA methylation-regulated genes, including ALDOB, CDC20, CXCL2, GC, GPD1, HSD17B13, SLC1A5, SPP1, TTC36, TTYH3, UHRF1, HEY1, and SKA2, were related to the OS or DFS of the patients with liver cancer (Figures S1 and S2). Therefore, these 13 genes were used for subsequent analyses.

Disease Association Analysis

We uploaded the 13 potential key DNA methylation-regulated genes to the VarElect online tool and searched for “hepatoblastoma,” “liver cancer,” and “cancer.” The top three genes were SPP1, UHRF1, and HEY1 (Table 2). In particular, we found that SPP1 and HEY1 were associated with HB (Table 2). Thus, SPP1, UHRF1, and HEY1 were regarded as the key DNA methylation-regulated genes.

Table 2 Correlation Between Genes and Disease

Construction of the PPI Network and Biological Function Analysis

To further explore the potential molecular mechanisms of SPP1, UHRF1, and HEY1, we first established the PPI network based on SPP1, UHRF1, HEY1, and other DEGs. As shown in Figure 4A, the PPI network included SPP1, UHRF1, HEY1, and 143 other DEGs. Notably, we found that all three key DNA methylation-regulated genes could affect TP53, while both UHRF1 and HEY1 could influence BMP4 (Figure 4A). Moreover, we found that the genes in the PPI network were mainly involved in cell division (Figure 4B).

Figure 4 PPI and GeneMANIA network analysis. (A) PPI network of key DNA methylation-regulated genes and DEGs. (B) GeneMANIA network showed the interactions and biological process of key DNA methylation-regulated genes and DEGs.

To fully analyze the roles of SPP1, UHRF1, and HEY1 in HB, we further investigated the biological functions of the DEGs. In terms of BPs, the DEGs were mainly involved in small molecule catabolic processes, carboxylic acid biosynthetic processes, and organic acid biosynthetic processes (Figure 5A). In terms of cellular components, the DEGs were mainly related to collagen-containing extracellular matrix, endoplasmic reticulum lumen, and cytoplasmic vesicle lumen (Figure 5A). In terms of molecular function, the DEGs were mainly associated with coenzyme binding, oxidoreductase activity, and iron ion binding (Figure 5A). In terms of the KEGG pathways, the DEGs were mainly related to the cancer-related signaling pathways, such as the cell cycle, PPAR signaling pathway, and p53 signaling pathway (Figure 5B). Notably, we found that both the cell cycle and p53 signaling pathway were related to TP53, which might be regulated by SPP1, UHRF1, and HEY1 (Figures 5B and 6A, B).

Figure 5 GO annotation and KEGG pathway analysis of DEGs. (A) the GO terms related to DEGs. (B) the enriched KEGG pathways in DEGs.

Figure 6 The regulatory mechanisms of cell cycle and P53 signaling pathway. (A) cell cycle signaling pathway. (B) P53 signaling pathway.

Methylation Levels of the Key DNA Methylation-Regulated Genes

We further investigated the methylation sites and regions of SPP1, UHRF1, and HEY1. As illustrated in Figure 7, SPP1, UHRF1, and HEY1 exhibited three, four, and five hypomethylated sites, respectively, in the HB samples when compared with those in the normal liver samples. Peculiarly, the methylation sites of UHRF1 and HEY1 were located in the body regions of the genes; the methylation site of SPP1 was located in the 1stExon or 5ʹUTR region; and the cg11701868 site (chr19: 4912006) of SPP1 was located in the transcription enhancer region (Figure 7).

Figure 7 The methylation sites and regions of key DNA methylation-regulated genes and the Circos plot showed the relationship between DNA methylation and gene expression. (A) the methylation sites and regions of SPP1. (B) the methylation sites and regions of UHRF1. (C) the methylation sites and regions of HEY1. (D) the Circos plot showed the relationship between DNA methylation and gene expression; The first layer (outermost) is the chromosome; The second layer is the heatmap of gene expression in HB and normal liver samples with the outer circle as normal liver samples and the inner circle as HB samples; The third layer is the histogram of the gene expression change (logFC) in HB samples compared with normal liver samples, and blue represents downregulation and red represents upregulation; The fourth layer is the histogram of the mean difference of DMRs in in HB samples compared with normal liver samples, and blue represents demethylation and red represents supermethylation; The fifth layer is the number of CpG in the DMR, and the horizontal axis is chromosome and the vertical axis is CpG number.

Expression Levels and Diagnostic Significance of the Key DNA Methylation-Regulated Genes

We found that SPP1, UHRF1, and HEY1 were upregulated in the HB samples compared with those in the normal liver samples from the GSE131329 dataset (Figure 8A). Moreover, we found that all of them were also upregulated in the HB samples compared with those in the normal liver samples from the GSE133039 dataset (Figure 8B). Next, we investigated whether SPP1, UHRF1, and HEY1 could differentiate the HB samples from the normal liver samples. Interestingly, both the ROC curves for SPP1, UHRF1, and HEY1 from the GSE131329 and GSE133039 datasets showed that all areas under the curves were greater than 0.7 (Figure 8C and D). Furthermore, we found that UHRF1 was upregulated in the HB samples with metastasis compared with that in the HB samples without metastasis (Figure S3).

Figure 8 Expression levels and diagnostic significance of the key DNA methylation-regulated genes. (A) expression levels of key DNA methylation-regulated genes in the GSE131329 dataset. (B) expression levels of key DNA methylation-regulated genes in the GSE133039 dataset. (C) ROC curves in the GSE131329 dataset. (D) ROC curves in the GSE133039 dataset. **P < 0.01, ***P < 0.001, ****P < 0.0001.

Discussion

HB is a markedly heterogeneous tumor because of its composition, including many tumors originating from different immature liver precursors, such as hepatocytes and biliary and other epithelial or mesenchymal cells.^24,25 It has been reported that HB patients with different histological subgroups can have a variety of clinical outcomes.^25–28 Moreover, other clinicopathological parameters, including tumor stage, distant metastasis, multifocality, patient age, and birth weight, can also affect the prognosis of HB.^29,30 However, the molecular mechanisms underlying HB have not been fully elucidated. Increasing evidence has revealed that comprehensive analysis of DNA methylation and gene expression can provide new insights into the molecular mechanisms of tumorigenesis and tumor development and identify novel biomarkers.^31–33 Therefore, we aimed to further analyze the molecular mechanism and screen for novel biomarkers of HB by integrating DNA methylation and gene expression data.

In our study, we first identified 24,460 DMSs, including 4209 hypermethylated and 20,251 hypomethylated sites, in the HB samples compared with those in the normal liver samples (Figure 2C and D), indicating that patients with HB exhibit global DNA hypomethylation, which is consistent with the results of current studies.^19,22,34,35 Moreover, 22 hypomethylated and upregulated genes and 19 hypermethylated and downregulated genes were regarded as the DNA methylation-regulated genes (Figure 3A and B). We also found that ALDOB, CDC20, CXCL2, GC, GPD1, HSD17B13, SLC1A5, SPP1, TTC36, TTYH3, UHRF1, HEY1, and SKA2 were associated with the survival of the patients with liver cancer. Thus, we speculated that these genes may play key roles in HB. Finally, we identified three genes (SPP1, UHRF1, and HEY1) as the key DNA methylation-regulated genes based on the correlation found between the genes and diseases. Although methylation and transcriptome data were derived from different samples, we found that SPP1, UHRF1, and HEY1 were upregulated in the HB samples compared with those in the normal liver samples from both the GSE131329 and GSE133039 datasets (Figure 8A and B); this indicates that SPP1, UHRF1, and HEY1 may play important roles in HB. In addition, the ROC curves indicated that SPP1, UHRF1, and HEY1 could be used as diagnostic biomarkers of HB (Figure 8C and D).

To our knowledge, SPP1 has rarely been studied in patients with HB; our study is the first to show that SPP1 might be associated with HB through abnormally methylated regulation. Consistently, it has been suggested to be upregulated in different malignancies, such as colorectal cancer, thyroid cancer, and hepatocellular carcinoma.^31,36,37 Interestingly, SPP1 is hypomethylated and upregulated in thyroid cancer and hepatocellular carcinoma.^31,37 It can also act as a diagnostic and prognostic biomarker of hepatocellular carcinoma.³¹ It has been revealed that UHRF1 is overexpressed in HB samples compared with that in control samples and can serve as a prognostic biomarker and potential molecular target for HB.^38,39 Similarly, HEY1 has been suggested to play a key role in glioblastoma and act as a potential target for future glioblastoma therapy.⁴⁰ It is also associated with metastasis of breast cancer.⁴¹ In summary, our study found that SPP1, UHRF1, and HEY1 may play key roles in HB, and their expression changes might be caused by methylation.

The PPI network showed that SPP1, UHRF1, and HEY1 could affect TP53 and that both UHRF1 and HEY1 could influence BMP4 (Figure 4A). TP53 has been reported to play a key role in HB.^42,43 Moreover, BMP4 has been reported to affect the progression of bladder cancer and breast cancer.^44,45 Therefore, we speculated that SPP1, UHRF1, and HEY1 may affect the development of HB by regulating TP53 and BMP4. Conversely, the functional enrichment analysis indicated that all three of them may affect the development of HB by regulating TP53-regulated pathways, such as the cell cycle and p53 signaling pathway (Figures 5B and 6A, B). Thus, SPP1, UHRF1, and HEY1 may affect the development of HB by regulating the cell cycle and p53 signaling pathway.

Conclusion

In conclusion, we identified three key DNA methylation-regulated genes, including SPP1, UHRF1, and HEY1, which may be regarded as novel and potential diagnostic biomarkers and therapeutic targets for HB. However, further studies are needed to concentrate on the regulatory mechanisms of these genes in HB occurrence and development.