Integrated analysis of gene expression and DNA methylation profiles in ovarian cancer

Background Ovarian cancer is an epithelial malignancy that intrigues people for its poor outcome and lack of efficient treatment, while methylation is an important mechanism that have been recognized in many malignancies. In this study, we attempt to assess abnormally methylated gene markers and pathways in ovarian cancer by integrating three microarray datasets. Methods Three datasets including expression (GSE26712 and GSE66957) and methylation (GSE81224) datasets were accessed. GEO2R platform was used to detect abnormally methylated-differentially expressed genes. Protein-protein interaction (PPI) networks were built and analysed for hypermethylated and hypermethylated differentially expressed genes using Cytoscape software and Mcode app. GEPIA and cBioPortal platforms were used to validate the expression of the hub genes and the correlation between their mRNA expressions and methylation levels. Kaplan Meier-plotter platform were used to assess the prognostic significance of the hub genes. Results Six hundred eighty-one hypomethylated-upregulated genes were detected and involved in Rap1 signaling pathway, biosynthesis of amino acids, endocrine resistance, apoptosis, pathways in cancer. The hub genes were TNF, UBC, SRC, ESR1, CDK1, PECAM1, CXCR4, MUC1, IKBKG. Additionally, 337 hypermethylated-downregulated genes were detected and involved in pathways in cancer, focal adhesion, sphingolipid signaling pathway, EGFR tyrosine kinase inhibitor resistance, cellular senescence. The hub genes were BDNF, CDC42, CD44, PPP2R5C, PTEN, UBB, BMP2, FOXO1, KLHL2. TNF, ESR1, MUC1, CD44, PPP2R5C, PTEN, UBB and FOXO1 showed significant negative correlation between their mRNA expressions and methylation levels. TNF, ESR1 and FOXO1 showed prognostic significance. Conclusions Two novel gene networks were found for ovarian cancer. TNF, ESR1, MUC1 and FOXO1 are our candidate genes that might take part in ovarian cancer progression in an epigenetic approach, TNF, ESR1 and FOXO1 may serve as potential markers for ovarian cancer prognosis evaluation.


Background
Ovarian cancer, an epithelial malignancy, is reported to be the most common lethal malignancy among gynecologic cancer [1]. And the outcome of the patients is poor due to late diagnosis for lack of early signs and symptoms [2]. Ovarian cancer is prone to metastasis and recurrence [3], but the pathogenesis is still unclear, the suspected etiology include ovulation, hormones, genetics, environmental factors [4][5][6].
The standard treatment for ovarian cancer include surgery and chemotherapy, other treatments include radiation, hormone, immunotherapy, but the survival rate for it is still low (https://www.cancerresearchuk.org/ about-cancer/ovarian-cancer/survival) due to advanced stage when diagnosed and frequent recurrence (which often accompanied with increases chemoresistance) [5], thus more attention should be paid to ovarian cancer.
DNA methylation, the most studied epigenetic mechanism, is reported to be related to mRNA and miRNA expression regulation, thus contribute to cancer initiation or progression [7]. Recently, more and more studies indicated that abnormal gene methylation in promoter regions is involved with chemical therapy and targeting therapy of ovarian cancer [8][9][10].
In this study, DNA methylation datasets in ovarian cancer were screened. A series of bioinformatics tools were used for integrated analysis and detection of hub genes. Then the levels of hub genes, and the correlation between their expression and methylation level of them were confirmed by GEPIA [11], and cBioPortal [12,13] platforms, with the data from The Cancer Genome Atlas (TCGA) and GTEx (The Genotype-Tissue Expression). These confirmed genes are our candidate genes for a deeper study later on in ovarian cancer progress.

Gene datasets and differentially expressed genes identification
Gene expression datasets: GSE26712 (185 ovarian cancer samples, 10 normal samples), GSE66957 (57 ovarian cancer samples, 12 normal samples) and gene methylation datasets: GSE81224 (10 ovarian cancer samples, 5 normal samples) were screened from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/ geo/) for the biggest number of cases and the data provided are normalized that no further adjustment is needed. GEO2R platform were used to detect differentially expressed genes (DEGs) and abnormally methylated genes between normal and ovarian cancer samples. The parameter for DEGs and abnormally methylated genes were set with |t| > 2, P < 0.05. Then the DEGs and abnormally methylated genes were processed in Funrich software (http://www.funrich.org) for integrated analysis and Venn diagram visualization.

PPI network construction and modular analysis
STRING was used to build the protein-protein interaction (PPI) network (https://string-db.org/) [14], with minimum required interaction score set in 0.4. The result data was imported into Cytoscape [15] for visualization, subsequently, app Mcode (i.e., Molecular Complex Detection) was used to perform module analysis (Node score cut off > 3.5) for the PPI network complex. Hub genes were determined by connectivity degree > 27 for the hypomethylated-upregulated (HOUP) genes and connectivity degree > 12 for the hypermethylateddownregulated (HEDW) genes, the calculation of connection degrees of genes was performed using Microsoft Office Excel.

GO and Reactome, KEGG pathway enrichment analysis
The STRING website described above also helped in the integrated analysis of biological meaning of the proteins (or genes), Gene Ontology (GO) analysis (including the biological process (BP), cellular component (CC), and molecular function (MF)), Reactome pathway and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted for the selected genes, P value < 0.05 was considered statistically significant, and top 10 of BP, CC, MF items and top 5 of KEGG, Reactome pathways were illustrated using a web tool (http://www.ehbio.com/ImageGP/index.php/Home/ Index/index.html).
Validation of the expression of hub genes and correlations between methylation and mRNA levels in TCGA samples GEPIA platform was used to validate the expression levels of the hub genes and |log2FC| > 1 and P value < 0.05 was considered statistically significant. Methylation data TCGA 27 k methylation data were selected to assess the correlations between methylation levels and corresponding mRNA levels using cBioPortal platform (TCGA 450 k methylation data were not chosen because it contains only10 cases).

Prognostic analysis of hub genes in ovarian cancer tissue samples
The overall survival (OS), progression-free survival (PFS) and post progression survival (PPS) curves of each hub gene were drawn using the online platform, the Kaplan Meier-plotter (http://kmplot.com/analysis/index.php) [16]. Both logrank P value and Hazard Ratio (HR, and 95% confidence intervals) were evaluated, P < 0.05 was considered statistically significant.

Results
Detection of abnormal methylated-differentially expressed genes in ovarian cancer GSE26712 and GSE66957 (expression datasets) and GSE81224 (methylation datasets) were analyzed with GEO2R, two Veen diagrams were drawn and 681 HOUP genes and 337 HEDW genes were detected from these three datasets (Fig. 1).

Biological function and pathways enrichment analysis
The HOUP and HEDW genes were imported into STRING, respectively. Biological classification was evaluated for the two lists of genes including GO term enrichment analysis, KEGG pathways and Reactome pathways enrichment analysis.
The HOUP genes were enriched in biological process such as positive regulation of biological process, positive regulation of cellular process, cellular process, response to organic substance, response to stimulus. GO cell component analysis showed that the genes were significantly enriched in cytoplasmic part, cytoplasm, intracellular part, intracellular, cell. GO molecular function analysis showed significantly enrichment in.
protein binding, binding, enzyme binding, carbohydrate derivative binding, catalytic activity. KEGG and Reactome pathways analysis indicated enrichment in Rap1 signaling pathway, biosynthesis of amino acids, endocrine resistance, apoptosis, pathways in cancer, immune system, innate immune system, cytokine signaling in immune system, signal transduction, signaling by receptor tyrosine kinases (Fig. 2).
The HEDW genes were enriched in biological process such as regulation of anatomical structure morphogenesis, negative regulation of cellular process, negative regulation of biological process, regulation of biological quality, positive regulation of cellular process. GO cell component analysis showed that the genes were significantly enriched in intracellular, cytoplasm, intracellular part, cytosol, cytoplasmic part. GO molecular function analysis showed significantly enrichment in binding, protein binding, enzyme binding, actin binding, cytoskeletal protein binding. KEGG and Reactome pathways analysis indicated enrichment in pathways in cancer, focal The mRNA expression levels of the HOUP hub genes (analyzed by GEPIA platform). TNF, ESR1, CDK1, CXCR4 and MUC1 were confirmed significantly higher in ovarian cancer samples than that in normal samples adhesion, sphingolipid signaling pathway, EGFR tyrosine kinase inhibitor resistance, cellular senescence, negative regulation of the PI3K/AKT network, metabolism, metabolism of lipids, signal transduction, diseases of signal transduction (Fig. 2).
(MUC1), inhibitor of nuclear factor kappa B kinase regulatory subunit gamma (IKBKG). Three significant modules (Fig. 3b-d) were detected using app Mcode. Function analysis of the top three module genes using STRING showed enrichment in positive regulation of nitrogen compound metabolic process, positive regulation of macromolecule metabolic process, positive regulation of cellular metabolic process, positive regulation of cellular protein metabolic process, regulation of cellular process.
The PPI network graph for the HEDW genes was illustrated in Fig. 4a, 334 nodes and 772 edges were included in the network. Nine nodes with the highest connectivity degrees were judged as hub genes, including brain-derived neurotrophic factor (BDNF), cell division control protein 42 homolog (CDC42), CD44 antigen (CD44), Serine/threonine-protein phosphatase 2A 56 kDa regulatory subunit gamma isoform (PPP2R5C), phosphatase and Tensin homolog (PTEN), polyubiquitin-B (UBB), bone morphogenetic protein 2 (BMP2), forkhead box protein O1 (FOXO1), Kelchlike protein 2 (KLHL2). Two significant modules (Fig.  4b-c) were detected by app Mcode. Function analysis showed enrichment in protein modification by small protein conjugation, post-translational protein modification, protein ubiquitination, protein polyubiquitination, mitotic cell cycle. Fig. 6 The mRNA expression levels of the HEDW hub genes (analyzed by GEPIA platform). FOXO1 was significantly lower in ovarian cancer tissues than that in normal tissues

The validation of hub genes expression levels in TCGA ovarian cancer samples
To validate the expression levels of the HOUP and HEDW hub genes, GEPIA platform was used with the data origin from TCGA and GTEx. The results of the HOUP hub genes confirmed that the expression levels of TNF, ESR1, CDK1, CXCR4 and MUC1 (Fig. 5), indicating these 5 hypomethylated genes were activated in ovarian cancer development. For the HEDW hub genes, FOXO1 was significantly lower in ovarian cancer tissues than that in normal tissues (Fig. 6), indicating it might participate in ovarian cancer development with a different approach.

Methylation level and mRNA expression level correlation analysis of the hub genes in ovarian cancer
For further validation of the effect of methylation on mRNA expression, we performed a correlation analysis for the two lists of hub genes with cBioPortal, negative correlations were found in TNF, ESR1, MUC1, CD44, PPP2R5C, PTEN, UBB and FOXO1 (Table 1) between DNA methylation levels and corresponding mRNA levels.
Taken together, TNF, ESR1, MUC1 and FOXO1 are selected as our candidate genes for further research.

The prognostic significance of validated hub genes in ovarian cancer
To estimate the prognostic significance of abnormal expressed TNF, ESR1, MUC1 and FOXO1, the survival time (include OS, PFS and PPS) and gene expression levels were acquired from Kaplan Meier-plotter website. The analysis results showed that higher level of TNF is related to longer OS and PPS time (based on the survival curves and logrank P value), and higher level of ESR1 and lower level of FOXO1 are potential protective factors for ovarian cancer patients' survival (based on the HR and 95% confidence intervals value) (Fig. 7). These three key genes may serve as potential markers for ovarian cancer prognosis evaluation.

Discussion
Ovarian cancer is fatal, with high rate of metastasis, drug resistance and recurrence. But the progression and treatment for it remains unsolved. Recent studies focus on its proliferation, migration, and invasion showed relevance with miRNAs [17][18][19], lncRNAs [20,21], and EMT (epithelial-mesenchymal transition) [22]. And there are also studies focus on treatment for it, such as PARP inhibitors [23,24], hyperthermic intraperitoneal chemotherapy [25] and nanotherapy, but until now, survival rate of ovarian cancer remains low. Inspiringly, more and more microarray were applied for the investigation of its progression and drug resistance [26,27] presented us with massive information.
In our study, three datasets were used and 681 HOUP genes and 337 HEDW genes were identified. The HOUP genes were enriched in biological processes of positive regulation of biological process, positive regulation of cellular process, cellular process, response to organic substance, response to stimulus, it is reasonable. Cell component analysis for the HOUP genes showed enrichment in cytoplasmic part, cytoplasm, intracellular part, intracellular, cell. For molecular function, significant enrichment was found in protein binding, binding, enzyme Table 1 Methylation analysis of hub genes in TCGA ovarian  cancer. Negative correlations were found between TNF, ESR1,  MUC1, CD44,  binding, carbohydrate derivative binding, catalytic activity. The analysis is also reasonable because material and energy metabolism and protein binding are common activities in tumor cells including ovarian cancer. KEGG pathways analysis for the HOUP genes showed enrichment in Rap1 signaling pathway, biosynthesis of amino acids, endocrine resistance, apoptosis, pathways in cancer, immune system, innate immune system, cytokine signaling in immune system, signal transduction, signaling by receptor tyrosine kinases. It is also reasonable because immune, cytokine, biosynthesis, apoptosis in accordance with the nature of cancer. KEGG pathways analysis for the HEDW genes showed enrichment in pathways in cancer, focal adhesion, sphingolipid signaling pathway, EGFR tyrosine kinase inhibitor resistance, cellular senescence, negative regulation of the PI3K/ AKT network, metabolism, metabolism of lipids, signal transduction, diseases of signal transduction. Related studies about sphingolipid [28], EGFR tyrosine kinase inhibitor [29], cellular senescence [30], PI3K/AKT network [31] and metabolism of lipids [32] in ovarian cancer were found. In summary, our biological analysis results were logical and in accordance with previous researches. GEPIA platform were used to validate of the mRNA expression of 9 HOUP hub genes, named TNF, UBC, SRC, ESR1, CDK1, PECAM1, CXCR4, MUC1 and IKBKG, with data from TCGA and GTEx, TNF, ESR1, CDK1, CXCR4 and MUC1 were confirmed significantly higher in ovarian cancer samples than that in normal samples. To test whether DNA hypomethylation caused abnormal high expression, we explored the cBioPortal platform, and TNF, ESR1, MUC1 were confirmed with negative correlations with mRNA levels ( Table 1). The mRNA expression of 9 HEDW hub genes, i.e., BDNF, CDC42, CD44, PPP2R5C, PTEN, UBB, BMP2, FOXO1, KLHL2, were also validated using GEPIA. FOXO1 was Fig. 7 The prognostic significance of validated hub genes in ovarian cancer. The analysis results showed that higher level of TNF is related to longer OS and PPS time, and higher level of ESR1 and lower level of FOXO1 are potential protective factors for ovarian cancer patients' survival (judged by both P value and 95% confidence intervals of HR). The probe ID used were listed below: 207113_s_at (TNF), 205225_at (ESR1), 213693_s_at (MUC1), 202724_s_at (FOXO1) significantly lower in ovarian cancer samples than that in normal samples. Then cBioPortal platform was used again, and FOXO1 were confirmed with negative correlation with their mRNA levels (Table 1).
TNF, ESR1, MUC1 are HOUP genes, suggesting a potential role in ovarian cancer progression. According to the literature, TNF is widely studied in a variety of cancers, including ovarian cancer, it is also a major mediator of inflammation, there are also report trying to combine its effect in inflammation and carcinogenesis [33,34], while no report concerning methylation of TNF in ovarian cancer was found. ESR1, reported in many solid malignancies [35], is frequently methylated in ovarian cancer cell lines [36], its methylation was also detected in cfDNA of high-grade serous ovarian cancer patients [37]. High expression of MUC1 is associated with epithelial ovarian cancer progression [38], hypomethylated MUC1 in poorly-differentiated ovarian cancer indicated functions in tumorigenesis, metastasis, invasion and migration [39], and a combination of MUC1 vaccination and anti-PD-L1 blockade resulted improved survival in ovarian cancer patients [40].
FOXO1, a HEDW genes, effects in many carcinomas including ovarian cancer [41], it is also related to drug resistance in ovarian cancer [42,43], while it is not reported in ovarian cancer concerning methylation of FOXO1.

Conclusions
Taken together, with the datasets of gene expression and methylation, our study presented an integrated bioinformatics analysis of abnormally methylated DEGs for ovarian cancer. Hub genes including TNF, ESR1, MUC1 and FOXO1 might be potential targets for diagnosis or treatment of ovarian cancer in an epigenetic approach, TNF, ESR1 and FOXO1 may serve as potential markers for ovarian cancer prognosis evaluation.