Development of Potential Prognostic Biomarkers Based on DNA Methylation-Driven Genes for Patients with Endometrial Cancer

Introduction

Endometrial cancer (EC) is a common and lethal gynecologic malignancy in the world, as evidenced by the 320,000 new cases and 72,000 deaths reported in 2012¹ (EC is histopathologically classified into two types. Type I EC, attacking more than 70% of all cases, has generally a low grade and association with unopposed estrogen exposure. Type II EC has a high grade and often develops into papillary serous or clear cell cancer highly prone to relapse and metastasize. The first-line treatment is surgery if the patient shows indicative histopathology (tumor type, grade and stage). Chemotherapy, targeted drugs and hormone therapy are also used to increase survival and decrease recurrence, radiation therapy is another option for patients surgically inoperable. The application of these drugs delayed the progress of the disease and prolonged the survival of patients to a certain extent.² However, limited curative effect, serious toxic and side effects and easy to produce drug resistance led to poor prognosis of patients with advanced or recurrent EC, which is far from changing the current situation of EC treatment. With the deepening of the research on the molecular mechanism of EC, a series of new targets such as VEGF, PARP, PD-1/PD-L1 have been found, which brought hope to EC patients.³ For early-stage EC patients, initial primary surgical resection brings a good prognosis, the five-year survival rate up to 70%, but for advanced-stage patients, this rate falls to 20%.⁴ So, early diagnosis and management of EC as well as effective targeted drugs are necessary.

Malignancy is backed by complicate molecular mechanisms. DNA methylation is essential to the initiation and progression of malignancies.⁵ Current microarray and high-throughput sequencing technologies allow us to observe the expression and methylation levels of thousands of human genes simultaneously. Meanwhile, bioinformatic analysis offers an insight into the molecular mechanism of cancers on the whole genomic level. The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/publications/ publication guide lines) is an open database that provides cancer genetic and epigenetic profiles. MethylMix, an algorithm in R, can be used to screen out specific methylation-driven genes based on a β-mixed model.⁶ In the present study, the mRNA expression and methylation data of EC patients were downloaded from the TCGA database, and the methylation-driven genes were extracted by using the MethylMix algorithm. On genome level, we explored the correlation between the oncogenesis and DNA aberrant methylation in EC by MethylMix algorithm for the first time. The two key methylation-driven genes in EC can be used for new target drugs.

Materials and Methods

Data Collection and Analysis

The mRNA expression data of 587 samples (552 EC samples and 35 normal samples), methylation data of 485 samples (439 EC samples and 46 normal samples), and all their corresponding clinical information were downloaded from the TCGA database. The mRNA expression data and methylation data were normalized and analyzed in LIMMA package to obtain differentially expressed and aberrantly methylated genes.⁷ Then, we calculated the correlation between gene expression and gene methylation level by using MethylMix algorithm in R. After constructing a β-mixed model, we identified differentially expressed genes (DEGs) in EC and the EC-specific methylation-driven genes.

Functional Annotations for Methylation-Driven Genes

The functional annotations and related pathways of identified EC-specific methylation-driven genes were investigated. Gene Ontology (GO) analysis of methylation-driven genes covered three aspects, including molecular function (MF), biological processes (BP) and cellular composition (CC) were analysed in Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (http://david.abcc.ncifc rf.gov/).⁸ The results of GO enrichment analysis were plotted in the GOplot R package (https://cran.r-project.org/web/packages/GOplot/).⁹ In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed in ConsensusPathDB (http://cpdb.molgen.mpg.de) online software.¹⁰ The P-value <0.05 was set as the cutoff criterion in GO and KEGG analyses.

Construction of Protein–Protein Interaction (PPI) Network

To evaluate the interrelations between the products coded by these methylation-driven genes, these genes were uploaded to the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org/).¹¹

Risk Assessment Model and Risk Scoring

To screen out candidate methylation-driven genes with EC prognostic value, we constructed a linear risk assessment model by univariate and multivariate Cox analyses.¹² Firstly, the univariate Cox proportional hazards regression analysis was performed to investigate the overall-survival-associated driver genes. Those showing significant association were selected as candidate variables and then coped with multivariate Cox proportional hazards regression analysis. The prognostic index calculating formula was shown as follows:

Prognostic Index = expRNA1*βRNA1+expRNA2*βRNA2+expRNA3*βRNA3+ … expRNAn*βRNAn (expRNA was the expression level of each gene, and βRNA was the regression coefficient calculated by the multivariate Cox regression analysis).

According to the median prognostic value, EC samples were divided into high- and low-risk groups, and the overall survival rates of the low-risk group and high-risk group were calculated by Kaplan–Meier curve analysis in Survival R package, with the log-rank P < 0.05 was considered statistically significant. In addition, a time-dependent receiver operating characteristic (ROC) curve was plotted to estimate the predictive efficiency of the prognosis risk assessment model. Few expression profiling, sequencing, or array data have been published for different prognostic outcomes with EC. Therefore, according to the different expression levels of various subtypes, we identified data with different clinical subtypes to verify the key genes.

Survival Analysis and Mapping of Methylation Sites of Candidate Driver Genes

To identify key prognosis-associated genes in EC, the methylation level and methylation-expression levels of candidate driver genes in the risk assessment model were further analyzed via Survival R package. Moreover, we extracted relevant methylation sites of key genes from the downloaded EC methylation data. The correlation between the methylation and expression was determined, using |Cor| >0.35 as a cutoff criterion. Then, a joint time-dependent ROC curve was plotted to validate the predictive power of these candidate methylation-driven genes. Gene set enrichment analysis (GSEA) was introduced. In TCGA set validation, 546 EC samples were divided into two groups according to the median expression level of the candidate driver genes. To identify the potential function of these genes, GSEA (http://software.broadinstitute.org/gsea/index.jsp) was employed to determine the differentially expressed gene between the two groups.¹³ Annotated gene sets of HALLMARK in Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp) were chosen as the reference.¹⁴ The P-value <0.05 was set as the cutoff criterion. “joyplot” and “dotplot” functions in “Cytoscape” software were used to visualize the distribution.

Building the Nomogram

A nomogram and calibration curve was built by the “rms” package on R. The correctness was examined to check the consistency index between actual observation frequency and predicted probability. Then, we presented the predicted and observed results in the calibration curve to visualize the performance of the nomogram.

Results

Data of Methylation‑Driven Genes

The gene expression data of 587 samples (552 EC samples and 35 normal samples) and the methylation data of 485 samples (439 EC samples and 46 normal samples) extracted from TCGA database were put in the LIMMA software package. As a result, 127 methylation-driven genes of EC were screened out (|log Fold change (FC)| >0, adjusted-P < 0.05, |Cor| >0.45). The heatmaps of these aberrant methylation-driven genes are shown in Figure 1.

Figure 1 Heatmaps of EC-related aberrant methylation-driven genes. (The color from green to red shows a progression from hypomethylation to hypermethylation).

Biofunctions of Driver Genes

As to the 127 identified driver genes, the results of GO analysis showed that those in BP group were mainly involved in “Regulation of transcription, DNA-templated and Transcription, DNA-templated”; those in MF group in “Transcription factor activity, sequence-specific DNA binding, Metal ion binding, Nucleic acid binding and DNA binding”; those in CC group in “Intracellular” (Figure 2). The KEGG Pathway analysis revealed that these 127 methylation-driven genes were mainly enriched in “Generic Transcription Pathway, Gene expression (Transcription), RNA Polymerase II Transcription” (Figure 3).

Figure 2 Results of GO enrichment analysis of 144 methylation-driven genes. (A) GO analysis divided methylation-driven genes into three functional groups: MF, BP, and CC; (B) GO enrichment significance items of methylation-driven genes in different functional groups; (C) distribution of methylation-driven genes for different GO-enriched functions.

Figure 3 The significant enriched pathways of methylation-driven gene. Notes: Node size: the number of genes; node color: P-value; edge width: percentage of shared genes; Edge color: genes from input.

PPI Network of Methylation-Driven Genes

The PPI network obtained from STRING is shown in Supporting Figure 1. A total of 101 nodes and 21 edges formed in the network of 127 methylation driver genes.

Prognosis-Related Genes Pinpointed by Risk Assessment Model

By performing univariate and multivariate Cox regression analyses, we screened out prognosis-related methylation-driven genes of EC and constructed a linear risk assessment model (Supporting Table 1). As a result, the top ten significant genes (MSX1, ERMN, RP11-968O1.5, DCAF12L1, ALS2CR11, PDCL2, AQP5, SIT1, ZNF354C, RPS6KA6 and RP11-760D2.5) were screened out in the univariate Cox proportional hazards regression analysis. These genes were further coped with multivariate Cox proportional hazards regression analysis, and four genes (RP11-968O1.5, DCAF12L1, MSX1 and ALS2CR11) were identified as candidate prognosis-related genes in EC. Next, the risk assessment model was constructed using each of the four candidate genes as an independent indicator. The prognostic index = (0.162*RP11-968O1.5)+ (0.074*DCAF12L1) + (−0.104*MSX1) + (0.085*ALS2CR11) (Supporting Figure 2). According to the median prognostic value, expression data of 534 EC samples and their corresponding clinical information were divided into a high-risk group contained 267 samples and a low-risk group contained 267 samples, respectively. The results of Kaplan–Meier survival curve analysis showed that the overall survival rate of the high-risk group was significantly lower than that in the low-risk group (Figure 4A). The five-year overall survival ROC curve showed the area under curve (AUC) of the four-gene prognostic risk assessment model was 0.719 (Figure 4B), suggestive of its predictive power for the five-year overall survival of EC patients. The distribution of risk score, survival status, and the expression levels of four genes in each patient were also analyzed (Figure 4C-E). In addition, we analyzed the relationship between the clinical parameters of each patient and the risk score of the four genes. The univariate and multivariate Cox proportional hazards regression showed that the age, tumor status and risk score were independent prognostic indicators of EC (Supporting Figure 3). The heatmap showed the expression levels of the four genes in high- and low-risk groups based on the TCGA dataset. We observed significant between-group differences in tumor status, grade, histological type, age and stage (P < 0.001) (Supporting Figure 4).

Figure 4 A linear risk model based on four candidate methylation-driven genes. (A) Kaplan–Meier survival curve of overall survival between high‐risk group and low‐risk group; (B) five-year survival time dependent ROC curve; (C–E) the distributions of the four-gene signature, survival status, and expression profiles of the four genes of patients in the training data set.

Key Genes in Different Subtypes

There was a significant difference in the expression of the key genes in different subtypes based on TCGA data, using one-way analysis of variance (ANOVA). Among the 4 grade subtypes, MSX1 had the highest expression level compared with the RP11-968O1.5, DCAF12L1, ALS2CR11 (Figure 5A). The consistent trend was displayed in 4 different stage subtypes and 3 different histological type subtypes (Figure 5B and C).

Figure 5 Expression levels of Candidate biomarkers in different subtypes. (A) in different grade subtypes; (B) in different stage subtypes; (C) in different histological type subtypes. **P < 0.01, ***P < 0.001, ns, not significant.

Stratification Analysis of Hub Genes and Clinical Features

The stratification analysis was performed according to histological type, grade, age, tumor status and stage. Patients were stratified into endometrioid subgroups, grade G3&G4 subgroup, stage I & stage II subgroup, stage III & stage IV subgroup, with tumor and tumor-free subgroup, age >60 subgroup and age ≤60 subgroup. For the patients in endometrioid subgroup, the survival time of patients in the high-risk group was significantly shorter than that of patients in the low-risk group (Figure 6A), consistent with the trends for the mix and serous subgroups, the grade G3&G4 subgroup, stage I & stage II subgroup, stage III & stage IV subgroup, with tumor subgroup, age >60 subgroup and age ≤60 subgroup (Figure 6B-H).

Figure 6 Survival time of patients in high-risk and low-risk group of different subgroups. (A) Endometrioid subgroup, (B) mix and serous subgroup, (C) grade G3&G4 subgroup, (D) stage I & stage II subgroup stage, (E) stage III & stage IV subgroup, (F) with tumor subgroup, (G) age >60 subgroup, (H) age ≤60 subgroup.

Building Predictive Nomogram

To achieve the goal of establishing a clinical method to predict the survival probability with EC patients, a nomogram was created using TCGA cohort so as to assess the probability of the 1-, 3‐ and 5‐year OS. The predictors of the nomogram contained 6 independent prognostic factors including stage, age, histological type, grade, tumor status and risk score (Figure 7A). The 45° line represented the best prediction. Calibration plots suggested that the nomogram performed well (Figure 7B). ROC curve analysis showed that the risk score AUC value of the 1-, 3‐ and 5‐year model was 0.632, 0.650 and 0.655. The clinical factors AUC value of the 1-, 3‐ and 5‐year model was 0.603, 0.651 and 0.640. Combining risk score and clinical factors, the AUC value of the 1-, 3‐ and 5‐year model was 0.628, 0.685, and 0.674 (Figure 7C).

Figure 7 The nomogram to predict 1-, 3‐ or 5‐year OS and prognostic value of hub genes. (A) The nomogram for predicting proportion of patients with 1-, 3‐ or 5‐year OS. (B) The calibration plots for predicting patient 1-, 3‐ or 5‐ year OS. Nomogram‐predicted probability of survival is plotted on the x‐axis; actual survival is plotted on the y‐axis. (C) 1-, 3‐ or 5‐ year ROC curve analysis of the hub genes according to risk score and clinical factors.

Reciprocity Between Methylation and Expression of Driver Genes

Kaplan–Meier survival analysis further illustrated that the methylation levels and joint methylation-expression levels of DCAF12L1 and MSX1 were associated with the overall survival rate of EC patients. Hypermethylated (or hypermethylated and lowly expressed) DCAF12L1 linked with a higher survival rate, and hypomethylated (or hypomethylated and highly expressed) MSX1 linked with a higher survival rate (Figure 8). In addition, the methylated sites and the gene expression of DCAF12L1 and MSX1 are shown in Figure 9. According to the results of five-year dependent ROC analysis, the AUC of DCAF12L1 and MSX1 was 0.704 and 0.73, respectively. The joint ROC analysis revealed the AUC of DCAF12L1-MSX1 was 0.867, which suggested both have a good EC-diagnosing efficiency (Figure 10A). We then coped DCAF12L1 and MSX1 with GESA analysis (Figure 10B), finding both were mainly associated with the KRAS signaling pathway.

Figure 8 The survival analysis of candidate methylation-driven genes. (A) The survival analysis of gene DCAF12L1 methylation; (B) the combination of gene DCAF12L1 methylation and expression; (C) the survival analysis of gene MSX1 methylation; (D) the combination of gene MSX1 methylation and expression.

Figure 9 (A–G) The correlation of methylation sites and gene DCAF12L1 expression; (H–M) the correlation of methylation sites and gene MSX1 expression.

Figure 10 (A) The joint ROC analysis revealed the AUC of DCAF12L1-MSX1; (B) the GSEA analysis of genes DCAF12L1 and MSX1.

Discussion

Cancer is a product of multiple factors, such as somatic, genetic, epigenetic, and endocrine aberrations. It has been proved that DNA methylation can alter gene expression in various diseases, including malignancies (Kulis and Esteller 2010). Thus, the abnormal methylation in cancer-associated genes has been investigated for teasing out sensitive biomarkers and effective therapeutic targets of cancers.^15,16 In the present study, we identified 127 methylation-driven genes possibly regulating the initiation and development of EC. GO analysis showed the top three enriched terms were “Regulation of transcription, DNA-templated, Transcription factor activity, sequence-specific DNA binding, and Transcription, DNA-templated”, suggesting that these driver genes may be key components in the transcription associated process and function. The abnormal methylation and expression levels of these genes can dysregulate transcription of oncogenes or anti-oncogenes. Moreover, GO analysis also revealed that these driver genes were enriched in terms “Metal ion binding, Nucleic acid binding and DNA binding” in MF group, indicating that products coded by these genes might associate with the function of metal ion, nucleic acid and DNA binding, which can provoke mutations and may lead to cancer reversibly (eg, conformational changes) or irreversibly (eg, cleavage).¹⁷ As the results of KEGG pathways showed, the genes were mainly enriched in “Generic transcription pathway, Gene expression (Transcription) and RNA polymerase II transcription” pathways. These findings indicated that most of these driver genes might play in EC-associated as gene regulators, such as the transcription factors and cofactors. The mutations and epigenetic modification of these regulators can trigger tumorigenesis.¹⁸ In addition, the PPI network also verified the interactions between these driver genes.

Our univariate and multivariate Cox analyses identified four candidate methylation-driven genes (RP11-968O1.5, DCAF12L1, MSX1 and ALS2CR11) that were used to construct an EC-prognostic risk assessment model. The EC samples were divided into high-risk and low-risk groups based on the prognostic values calculated by this model. Survival analysis showed significant difference in overall survival rate between the two groups. The AUC of five-year overall survival ROC of these four driver genes was 0.719. These findings suggest that the four-gene risk prognosis model can predict the overall survival rate of EC patients. However, further studies are needed to validate the efficiency of this risk model in clinical applications.

We found that the methylation levels and joint methylation-expression levels of DCAF12L1 and MSX1 were associated with the prognosis. Patients with hypermethylated DCAF12L1, or hypermethylated and lowly expressed DCAF12L1, displayed a higher survival rate. This finding suggests that DCAF12L1 can act as a tumor suppressor in EC, and the hypermethylation and down-regulation of this gene may suppress tumorigenesis of EC. DCAF12L1, also known as KIAA1892L or WDR40B, is a member of Ddb1- and Cul4-associated factors (DCAFs). DCAFs can regulate DNA replication in cell cycle.¹⁹ Meanwhile, we found patients with hypomethylated (or hypomethylated and highly expressed) MSX1 showed worse prognosis, indicating that gene MSX1 can play as a tumor promoter, and raising the methylation and expression of MSX1 may curb the tumor development. MSX1 has been widely considered as a tumor suppressor. Tao et al found that MSX1 suppressed glioblastoma cell migration and invasion by inhibiting the Wnt/β-catenin signaling pathway.^19,20 Yue et al demonstrated that methylation of MSX1 promoter down-regulated MSX1 expression, triggering oncogenesis by mediating Notch signaling in cervical cancer cell lines.²¹ However, up to date, the expressions and functions of two genes in EC have not been validated. Our study suggests that the two genes can serve as diagnostic or prognostic biomarkers and potential therapeutic targets for EC.

Meanwhile, the ROC curve analysis further validated that DCAF12L1 and MSX1 were independent prognosis-related genes, and their combination showed a higher prognostic value. The GSEA analysis found that KRAS signaling pathway was significantly activated in high-expressed groups KRAS is a RAS protein that can activate RAF family by binding to GTP to promote cell proliferation.²² Thus, mutation, overexpression and amplification of RAS proteins are key steps in the development of cancer. RAS mutation appears in around 20% of human tumors, most frequently in KRAS (85%).²³

Conclusion

In conclusion, four candidate methylation-driven genes (RP11-968O1.5, DCAF12L1, MSX1 and ALS2CR11) were identified as EC-prognosis-related. DCAF12L1 and MSX1 can serve as EC prognostic markers and drug targets. However, their role in EC development should be clarified with future experimental evidences.