Identification of regulatory role of DNA methylation in colon cancer gene expression via systematic bioinformatics analysis.

Colon cancer arises from the accumulations of genetic and epigenetic changes. Currently, profiles of DNA methylation and gene expression of colon cancer have not been elucidated clearly. This articles aims to characterize the profile of DNA methylation and gene expression of colon cancer systemically, and acquire candidate genes potentially regulated by altered methylation for this disease.Data were downloaded from The Cancer Genome Atlas database. Differentially methylated CpG sites (DMCs) and differentially methylated regions (DMRs) were calculated via COHCAP. Differentially expressed genes (DEGs) were identified by DESeq2. Weighted gene co-expression network analysis (WGCNA) package in R was applied for WGCNA.Data of 275 solid tumor tissues and 19 adjacent tumor tissues of colon cancer were obtained. A total of 1828 DMCs, including 1390 hypermethylated and 438 hypomethylated CpG sites, were identified between tumor and normal groups. A total of 789 DEGs, containing 435 upregulated genes and 354 downregulated genes were observed. It revealed that 8 DMRs-DEGs and 95 DMCs-DEGs pairs were significantly correlated. Furthermore, genes of yellow and brown modules from WGCNA were significantly correlated with tumor/normal status, and significantly enriched in peroxisome proliferator activated receptor signaling pathway, glutamatergic synapse, and neuroactive ligand-receptor interaction. Genes in the above 2 modules were also significantly enriched in DMCs or DMRs-associated genes. Specifically, ADHFE1, HAND2, and GNAO1 were hypermethylated and downregulated in colon cancer, suggesting that the low expression levels of these genes may be regulated by DNA hypermethylation. In addition, the 3 genes were involved in brown module of WGCNA, indicating their important roles in colon cancer.The investigation of the relationship between DNA methylation and gene expression may help to understand the effect of DNA methylation alteration on genes expression, especially gene co-expression network in the development of colon cancer. Genes such as ADHFE1, HAND2, and GNAO1 may be served as potential candidates for diagnosis and therapy targets in colon cancer.

Introduction

Colon cancer is a common malignant cancer worldwide and considered as the third most commonly diagnosed cancer in both men and women.[ Genetic and epigenetic alterations are involves in the development and progression of colon cancer. Approximately 75% of colorectal cancer arises from long-term accumulations of epigenetic alterations.[ In terms of genetic alterations, numerous chromosomal gains and losses are involved in colon cancer.[

In the past decades, epigenetic modification, especially DNA methylation regulation in carcinogenesis has recently become a focus for cancer researches. For example, the altered hypermethylation status of CpG islands (CGIs) in colon cancer is a common signature in colon cancer.[ Silencing of multiple tumor suppressor genes through DNA methylation pattern alterations plays a key role in the carcinogenesis of colon cancer.[ LINE-1 hypomethylation in early stage colon cancer is usually associated with a worse overall survival.[ However, the mechanism of DNA methylation alterations and their effect on gene expression has not been systematically elucidated.

Systematic networks of gene expression data can help to identify driver genes or pathway for the development of diseases. Moreover, correlation network analysis of genes is increasingly used in bioinformatics analysis. Based on weighted gene co-expression network analysis (WGCNA)[ that is a tool for systematic genetic data analysis, Horvath et al[ identified a novel gene named abnormal spindle microtubule assembly as from gene co-expression modules in glioblastoma, highlighting the use of WGCNA in recognizing crucial networks and molecular targets in cancer. Currently, little is known about the feature of co-expression network of colon cancer.

In this article, we comprehensively assessed genome-wide methylation and gene expression profiles associated with colon cancer with the large data provided by The Cancer Genome Atlas (TCGA) database, and tried to identify clusters (modules) highly correlated with colon cancer using the expression data of all samples via WGCNA. Finally, we focused the genes in the identified modules, and noted these genes in the DNA methylation and gene expression change to explain their role in colon cancer, which may be considered as diagnosis and therapy targets for colon cancer.

Materials and methods

Data collection

We downloaded DNA methylation (Illumina Infinium HumanMethylation450 BeadChip), mRNA expression (IlluminaGA_RNASeqV2.1.0.0), and clinical information of colon cancer from TCGA database.[ Approximately 275 solid tumor tissues (tumor group) and 19 corresponding adjacent tumor tissues (normal group) of colon cancer were obtained, and both DNA methylation and mRNA expression data were available for each of the 294 samples. This integrated study was a second analysis of TCGA data, so an ethical statement is not required from patients.

Identification of differentially methylated CpG sites and differentially methylated regions

The raw data of DNA methylation was preprocessed, we filtered the CpG sites (CpGs) for which DNA methylation data were missing in 80% samples, and finally 395,530 CpGs were remained. COHCAP package in R[ was used as a tool to calculate differentially methylated CpG sites (DMCs) and differentially methylated regions (DMRs). The threshold of DMCs was set as beta value >0.6 or beta value <0.3, delta value >0.2, false discover rate (FDR) <0.01. Furthermore, the threshold of DMRs was set as number of each region ≥2, FDR <0.05. Unsupervised hierarchical clustering analysis was performed to explore the characterization of DNA methylation pattern.

Identification of differentially expressed genes

Genes with count value showing 0 in >20% samples were filtered, and 16,853 genes were obtained for further analysis. Differentially expressed genes (DEGs) between tumor and normal group was calculated by DESeq2 package in R.[ The threshold for DEGs was set as |log2 (fold change)| >2.5 and FDR <0.01.

Construction of weighted gene co-expression network

The expression data of >16,853 genes were employed as input data to construct an unsigned weighted gene co-expression network (WGCN) via WGCNA package of R.[ The soft threshold was set as 10 to get better Scale-Free Topology Model Fit. The correlation between eigenvector of each module and tumor/normal status was investigated. Those with |coefficient of correlation| >0.2 and P value <.01 was considered as significant modules. Furthermore, genes that were most correlated with the eigenvector of their corresponding module were considered as hubs.

Integrated analysis of DNA methylation and gene expression

The DNA methylation and gene expression profiling was visualized by ClicO FS.[ In order to further explore the relationship between DNA methylation and gene expression, we mapped the DMRs to the nearest DEGs based on official 450 k array data annotation file via shell to obtain DMRs-DEGs pairs. Moreover, according to the official 450 k array data annotation file, the correlation analysis between DMCs and the nearest DEGs was analyzed via R. The intersections between above DEGs and selected modules were visualized by venny 2.1.0. The visualization of modules was performed by Cytoscape3.1.0 (http://clicofs.codoncloud.com).[

Genomic features and function enrichment analysis

To uncover the genomic features including CGI context and gene context of selected CpGs, Fisher exact test was performed. CGI context is composed of CGI, shore and shelf. Gene context contains TSS200, TSS1500, 5’UTR, 3’UTR, 1st Exon, and gene body according to the official 450 k array annotation file. Those with P value <.01 were considered as significant functional enrichments. Furthermore, enrichment analysis of Gene Ontology[ and Kyoto Encyclopedia of Genes and Genomes[ pathway was performed to explore the biological functions of selected genes. The online tool GeneCodis3 was used for function annotation analysis (http://genecodis.cnb.csic.es/analysis).[

Results

Differentially methylated CpG sites and differentially methylated region between tumor and normal groups in colon cancer

Based on the differential methylation analysis of colon cancer, we obtained 1828 DMCs, including 1390 hypermethylated and 438 hypomethylated CpGs. The top 10 DMCs ranked according to FDR values were cg20912169, cg20295442, cg18065361, cg04025964, cg16300300, cg09383816, cg09248054, cg17872757, cg22871668, and cg21900495, which were all hypermethylated. Specifically, the top 3 DMCs were all located in ADHFE1. Unsupervised hierarchical clustering analysis showed that the DNA methylation pattern of the tumor group was distinct from the normal group, and exhibited higher DNA methylation levels than did normal group (Fig. 1). Enrichment analysis of genomic features revealed that DMCs were significantly enriched in 5’UTR, body, 1st exon, TSS200, and TSS1500 for gene context, whereas for CGI context were island, shelf, and shore (Supplementary Table 1a–b). In addition, we further performed differential analysis for methylation regions to provide more robust results. Briefly, we identified 377 DMRs, including 359 hypermethylated regions and 18 hypomethylated regions.

Figure 1

Heat map of top 200 differentially methylated CpG sites (DMCs) [by false discover rate (FDR) rank] based on the unsupervised hierarchical clustering analysis of colon cancer.

Identification of differentially expressed genes between tumor and normal groups in colon cancer

After mRNA expression data were preprocessed, we performed differential expression between tumor and normal group. We obtained 789 DEGs, containing 435 upregulated genes and 354 downregulated genes.

Construction of weighted gene expression network

WGCNA was performed after preprocessing mRNA data. Finally, 19 modules were identified. To explore the characterization of those modules, we performed correlation analysis between modules and tumor/normal status (Supplementary Figure 1a). It showed that 9 modules were significantly correlated with tumor/normal status, among which yellow and brown modules were the most correlated with tumor/normal status (coefficient = −0.86, −0.61, respectively; P value = 6.50E − 86, 9.60E − 31, respectively) (Supplementary Table 2). The gene significance versus module membership plot for both yellow and brown modules was displayed in Supplementary Figure 1b, c. Functional enrichment analysis showed that genes in yellow module were significantly enriched in mineral absorption, nitrogen metabolism, proximal tubule bicarbonate reclamation, and peroxisome proliferator activated receptor (PPAR) signaling pathway. Meanwhile, genes in brown module were significantly enriched in vascular smooth muscle contraction, glutamatergic synapse, dilated cardiomyopathy, hypertrophic cardiomyopathy, cell adhesion molecules, calcium signaling pathway, and neuroactive ligand-receptor interaction. Hubs in each module were selected and displayed in Table 1. The genes network connections in yellow and brown modules were also constructed, and are displayed in Figure 2A–B.

Table 1

Hubs in yellow and brown modules.

Figure 2

: A, Visualization of weighted gene co-expression network analysis (WGCNA) network connections of the yellow module. It showed network connections whose topological overlap was above the threshold of 0.08, including 464 edges and 105 nodes. B, Visualization of WGCNA network connections of brown module. It showed network connections whose topological overlap was above the threshold of 0.3, including 534 edges and 82 nodes. The hubs genes were in the center of network, with bigger icons.

To better understand the potential effect of DNA methylation on gene expression, the profiles of DNA methylation and gene expression in colon cancer were visualized in Figure 3A. It was commonly known that the presence of methyl moieties in promoter region could inhibit corresponding gene expression, whereas in the gene body gene, could increase corresponding gene expression. In our study, DMRs-DEGs pairs changing in opposite direction was not found, and just obtained 8 genes associated with differential expression or DNA methylation according to the DMR mapping as described in the method section, consisting of 5 hypermethylated and downregulated genes, and 3 hypermethylated and upregulated genes (Table 2). According to the correlation analysis between DMCs and DEGs, 95 significantly correlated DMCs-DEGs (involving 39 DEGs) pairs were obtained (|coefficient of correlation| >0.2, P value <.05), including 63 negatively and 32 positively correlated pairs. In addition, from above DMRs-DEGs and DMCs-DEGs pairs, 40 differentially methylated DEGs were obtained. Functional enrichment analysis demonstrated that those 40 DEGs significantly enriched in regulation of transcription DNA-dependent, multicellular organismal development, noradrenergic neuron differentiation, positive regulation of cell proliferation, and synaptic transmission in Gene Ontology terms of Biological process, whereas for Kyoto Encyclopedia of Genes and Genomes pathways, Glutamatergic synapse, and Adipocytokine signaling pathway, Glutamatergic synapse, PPAR signaling pathway, and Neuroactive ligand-receptor interaction were significantly enriched.

Figure 3

Interaction of DNA methylation and gene expression of cancer. A, Circos plot of differentially methylated CpG sites (DMCs) and differentially expressed genes (DEGs). B, Venny plot of DMCs-associated, differentially methylated regions (DMRs)-associated genes, and genes in yellow and brown modules.

Table 2

The details of DMRs-differentially expressed genes pairs.

Next we further tried to uncover the relation between DNA methylation and gene expression in colon cancer. We integrated the data from DEGs from above DMRs-DEGs and DMCs-DEGs pairs and modules from WGCN.

We found that genes in yellow and brown modules were significantly enriched in DEGs or differentially methylated DEGs with P value <.01 (Supplementary Table 3a–b). The overlapping of DMCs-associated, DMRs-associated DEGs and genes in yellow and brown modules were performed (Fig. 3B).

It showed that there were 14 genes (MAB21L1, DCLK1, SFRP1, CBLN2, NOVA1, GRIK1, HAND2, RXRG, SORCS1, GNAO1, NBLA00301, PHOX2A, LONRF2, CHODL) involved in DMCs-DEGs pairs and brown module, 2 genes (ADHFE1, TNXB) involved in DMCs-DEGs, DMRs-DEGs pairs and brown module, and 4 genes (UNC5C, DMRTA1, DHRS9, CASR) involved in DMCs-DEGs pairs and yellow module. Specifically, HAND2 and GNAO1 were hubs of brown modules among the above referred 14 genes.

Discussion

In this article, we characterized the DNA methylation and gene expression profiling of colon cancer, and inferred the potential genes regulated by aberrant DNA methylation in colon cancer. We obtained 1828 DMCs (including 1390 hypermethylated and 438 hypomethylated CpGs in colon cancer) and 377 DMRs (including 359 hypermethylated regions and 18 hypomethylated regions), and the DNA methylation level was relatively higher in colon cancer than that in the normal tissues, which was in line with the previous study that reported that there was an increasing number of hypermethylated CGIs in colon cancer.[ The hypermethylation status in both CpGs and CGI level may lead to the dysregulation of specific genes, especially tumor suppression genes.[ To obtain more effective genes potentially regulated by DNA methylation, we performed comprehensive analysis of DMCs, DMRs, and genes in WGCN. A series of genes were implicated with the process of tumor formation and progression of colon cancer. In this study, we used differential analysis to compare DNA methylation or gene expression pattern between colon cancer and normal tissues, and performed functional analysis of DMCs or DEGs to connect the relationship between these changes and cancer development. However, differential analysis cannot assess the correlation among the DMCs or DEGs. A change of a diver gene in DNA methylation or gene expression level could lead to a series of changes from downstream or other related genes. Although WGCNA could determine clusters (modules) of highly correlated genes, and identify modules which are significantly correlated with the cancer status. In this study, we both used differential analysis and WGCNA to explore the crucial genes for the process of tumor formation or progression of colon cancer.

Alcohol alcoholism is confirmed to be a risk factor for colorectal cancer. In this article, alcohol dehydrogenase, iron containing 1 (ADHFE1) was significantly hypermethylated and downregulated in colon cancer. A previous study has also showed that ADHFE1 was hypermethylated in colorectal cancer.[ Also, ADHFE1 acted as one of the potential biomarkers for the rectal cancer.[ Moon et al[ confirmed that downregulation of ADHFE1 induced by alcohol-associated promoter hypermethylation, accelerated cell proliferation in colorectal cancer cells. In this study, ADHFE1 was also involved in brown module of WGCN, indicating the roles of ADHFE1 in colon cancer in a systematical aspect.

Heart and neural crest derivatives expressed 2 (HAND2) was a basic helix-loop-helix transcription factor, and played a very important role in the development and differentiation of heart and nervous system.[ In our study, HAND2 was significantly hypermethylated, downregulated in colon cancer. Moreover, HAND2 was one of the hubs in brown modules in WGCN. A recent study revealed that HAND2 was hypermethylated and downregulated in rectal adenocarcinoma.[ In addition, it was reported that a panel of 9 genes including HAND2 discriminated different types of cervical cancer based on the promoter methylation level, providing a satisfactory clues for case classification.[ The continuous proliferation of the endometrium was observed in mice with knockdown of HAND2.[ Here, bioinformatic results showed that the downregulation of HAND2 was potentially regulated by hypermethylation, which was first reported in colon cancer.

Guanine nucleotide–binding protein, α-activating activity polypeptide O (GNAO1) was hypermethylated and downregulated, and among hubs in brown modules of WGCN. GNAO1 is a member of subfamily of heterotrimeric G proteins, which switches signal transduction and their aberrant regulation can promote oncogenesis.[ Mutations of GNAO1 were involved in many diseases.[ It was reported that GNAO1 was downregulated, which increased cell proliferation in hepatocellular carcinoma.[ This is the first study to report the deregulation of GNAO1 in colon cancer. For the other genes of yellow and brown modules that were associated with differential gene expression or DNA methylation in colon cancer, such as MAB21L1, CBLN2, NOVA1, GRIK1, RXRG, SORCS1, NBLA00301, PHOX2A, LONRF2, and CHODL, their role was not reported in colon cancer, which may provide some clues for the further basic research in colon cancer.

Taken together, we found that a series of genes, such as ADHFE1, HAND2, and GNAO1, regulated by aberrant methylation plays a crucial role in colon cancer. These genes might serve as candidate diagnostic markers and therapy targets in colon cancer, which needs to be validated in more large cohorts to evaluate their potential as biomarkers in colon cancer. Moreover, the molecular mechanism of these genes needs to be uncovered by future studies.