Identification of methylation-driven genes prognosis signature and immune microenvironment in uterus corpus endometrial cancer

Abstract

Background

Uterus corpus endometrial cancer (UCEC) is the main malignant tumor in gynecology, with a high degree of heterogeneity, especially in terms of prognosis and immunotherapy efficacy. DNA methylation is one of the most important epigenetic modifications. Studying DNA methylation can help predict the prognosis of cancer patients and provide help for clinical treatment. Our research aims to discover whether abnormal DNA methylation can predict the prognosis of UCEC and reflect the patient's tumor immune microenvironment.

Patients and methods

The clinical data, DNA methylation data, gene expression data and somatic mutation data of UCEC patients were all downloaded from the TCGA database. The MethylMix algorithm was used to integrate DNA methylation data and mRNA expression data. Univariate Cox regression analysis, Multivariate Cox regression analysis, and Lasso Cox regression analysis were used to determine prognostic DNA methylation-driven genes and to construct an independent prognostic index (MDS). ROC curve analysis and Kaplan–Meier survival curve analysis were used to evaluate the predictive ability of MDS. GSEA analysis was used to explore possible mechanisms that contribute to the heterogeneity of the prognosis of UCEC patients.

Results

3 differential methylation-driven genes (DMDGs) (PARVG, SYNE4 and CDO1) were considered as predictors of poor prognosis in UCEC. An independent prognostic index was finally established based on 3 DMDGs. From the results of ROC curve analysis and survival curve analysis, MDS showed excellent prognostic ability in TCGA-UCEC. A new nomogram based on MDS and other prognostic clinical indicators has also been successfully established. The C-index of the nomogram for OS prediction was 0.764 (95% CI = 0.702–0.826). GSEA analysis suggests that there were differences in immune-related pathways among patients with different prognosis. The abundance of M2 macrophages and M0 macrophages were significantly enhanced in the high-risk group while T cells CD8, Eosinophils and Neutrophils were markedly elevated in the low-risk group. Meanwhile, patients in the low-risk group had higher levels of immunosuppressant expression, higher tumor mutational burden and immunophenoscore (IPS) scores. Joint survival analysis revealed that 7 methylation-driven genes could be independent prognostic factors for overall survival for UCEC.

Conclusion

We have successfully established a risk model based on 3 DMDGs, which could accurately predict the prognosis of patients with UCEC and reflect the tumor immune microenvironment.

Introduction

Uterus corpus endometrial cancer (UCEC), as the sixth most common malignant tumor in women worldwide, has an increasing morbidity and mortality year by year [1,2,3]. At present, surgery combined with adjuvant therapy (radiotherapy and/or chemotherapy) is the standard treatment for UCEC patients [4, 5]. As an emerging adjuvant treatment, immunotherapy brings hope to patients who cannot choose standard treatment methods [6, 7]. Nevertheless, the prognosis of UCEC patients is still very different. Some patients with recurrent or advanced tumors (stage III and IV) still have a poor prognosis [8].

Clinical staging is now the most mainstream factor predicting the prognosis of UCEC patients. However, it can often be found in the clinic that patients with the same stage have completely different survival outcomes [9, 10]. Therefore, finding reliable molecular markers is an urgent need for personalized treatment and improving clinical efficacy.

In addition, in terms of treatment options, for patients who suffered from tumor progression after standard treatment, targeted therapy and immunotherapy may be a promising rescue method. However, only a small percentage of patients can really benefit from immunotherapy. Therefore, it is also very urgent to find biomarkers that can predict the response of patients to immunotherapy.

DNA methylation is an important epigenetic modification catalyzed by DNA methyltransferases which add a methyl group to position 5 of the cytosines present in the transcription regulatory regions of genomic DNA using S-adenosyl methionine as donor molecule. Generally, hypermethylation correlates with gene silencing while hypomethylation means gene activation. Hypermethylation in promoter regions of some important genes, such as tumor suppressor genes and DNA repair genes, leading to down-regulation of their expression, may result in abnormal cell differentiation and regulation and DNA damage that cannot be repaired. Which is believed to be closely related to the occurrence and development of cancer [11]. With the deepening of research, abnormal changes in DNA methylation are considered potential biomarker for cancer risk assessment, early diagnosis and prognostic monitoring. As early as 2016, Balázs et al. had confirmed that abnormal DNA methylation was closely related to the prognosis of ER+/HER2− breast cancer [12]. Subsequently, abnormal DNA methylation was found to be an independent prognostic biomarker for a variety of tumors (including esophageal cancer, lung cancer and head and neck squamous cell carcinoma) [13,14,15]. In UCEC, many studies have shown that the occurrence and development of tumors are closely related to abnormal DNA methylation. The genome-wide DNA methylation profile of endometrial cancer clarifies the importance of DNA methylation changes in the process of endometrial cancer genesis [16]. In 2018, abnormal methylation levels of KLF4 and HS3ST2 genes were revealed to be important in the development of endometrial hyperplasia and UCEC [17]. The abnormal methylation of PTEN, APC, CDO1, SOX11 was found to be an independent risk factor for UCEC [18,19,20].

In recent years, with the continuous advancement of high-throughput sequencing technology, the combined analysis of multiple omics has become more and more widely used in tumor research. As a new tumor marker with stability, reversibility and high frequency, abnormal DNA methylation changes have been used by many researchers to combine gene expression profiles and clinical data to develop prognostic tumors model with more accurate prediction capabilities. By integrating methylation and mRNA expression profile data, some molecular prognostic models with high prognostic ability have been established and verified in gastric cancer and colon cancer [21,22,23]. However, no one has made similar research in UCEC.

In addition, studies have shown that DNA methylation profiles can reflect the immune microenvironment of tumors, which is closely related to the efficacy of immunotherapy [24]. Therefore, the purpose of this article is to use the methylation data, gene expression data and clinical data from The Cancer Genome Atlas (TCGA) Data Portal to construct a molecular model that can effectively predict the prognosis of UCEC patients and reflect the tumor immune microenvironment.

Material and methods

Ethics committee approval

All data were downloaded from the public databases hence it was not required to obtain additional ethical approval for our study.

Data acquisition and processing

Three types of TCGA-UCEC data (including patient clinical information, methylation data based on Illumina Infinium HumanMethylation 450 platform and level-3 mRNA expression files were obtained from TCGA (https://tcga-data.nci.nih.gov/tcga/; accessed August 2020). Then, we retained patients with complete clinical data including age, stage, grade and histological type, methylation data or mRNA expression data for further study.

We obtained the complete clinical data of 544 patients from the TCGA database. The clinicopathological characteristics of the patients included in the study are shown in Additional file 12: Table S1. At the same time, we obtained DNA methylation data of 439 tumor tissues and 46 adjacent normal tissues and mRNA expression data of 552 tumor tissues and 23 adjacent normal tissues from the TCGA database. In the end, we retained 429 patients with complete survival data, methylation data and mRNA expression data for further research.Text: “Additional files: As per journal instructions, we have ignored to process the Original data files "different expression gene data.xls"; "Supplymentary data.xls"; "Gene expression data.txt" and "methylation data.txt" given in the submission package. Kindly check and confirm.we have checked it.Thanks.

Identification of differential methylation-driven genes (DMDGs)

By using the “Limma” package in the R statistics software, we obtained the difference in mRNA expression between tumor samples and normal samples. mRNAs satisfying both |fold change (FC)| > 1 and adjusted p < 0.05 were defined as significantly differentially expressed. The MethylMix algorithm was used to integrate the DNA methylation data and differentially expressed mRNA data into the same sample for correlation analysis. Differential methylation-driven genes were finally identified by using the “MethylMix” package.

Functional enrichment analysis

By using the ‘clusterProfiler’ package in the R statistics software, we performed gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis on the obtained DMDGs. We set an adjusted p < 0.05 as the cut-off value to screen for meaningful functional events.

Establishment of a methylation-driven signature (MDS) based on DMDGs

After matching the complete survival data with methylation and mRNA expression data, we acquired 429 patients. Then we acquired data from 429 patients (entire cohort) randomly assigned 216 patients as the training cohort with the help of “caret” R package. The remaining 213 patients were assigned as the testing cohort. Then, in the training cohort, the univariate Cox regression analysis and the multivariate Cox regression analysis were used to identify genes with prognostic value from differentially expressed DMDGs. We used the “survival” package in R software for Cox regression analysis. Next, “glmnt” package was used to perform the Lasso Cox regression to avoid overfitting.

We then utilized multivariate Cox model to verify the correlation between gene expression of patient’s survival status. Besides, multivariate Cox model can also assign the weight of each involved gene in reflect the prognosis information of patients. We constructed a risk methylation-driven signature (MDS) based on the correlation coefficient of the key genes from the multivariate Cox regression. The formula of MDS was defined as follows: β₁ × gene₁ expression + β₂ × gene₂ expression + $\cdots$ + β_n × gene_n expression, where β corresponded to the correlation coefficient. According to our prognostic model, each patient would get a risk score. We set the median risk score as the cutoff value for dividing UCEC patients into a high-risk group and a low-risk group.

Evaluation of the MDS

Firstly, we conducted overall survival (OS) curves and time-dependent receiver operating characteristic (ROC) curve analysis in the training cohort. OS curves were plotted by the Kaplan–Meier (K–M) method. ROC curve analysis was obtained through the “survivalROC” package. Next, we used the same methods to evaluate the prognostic value of MDS in the testing cohort. In order to more comprehensively evaluate the application ability of the prognostic MDS we constructed in UCEC patients, we have drawn the OS curves of the high-risk group and the low-risk group of patients with different clinical characteristics in the entire cohort.

To further evaluate whether our model could be used as an independent prognostic factor, we included age, stage, histological type, grade and MDS as independent variables. And then univariate Cox regression analysis and multivariate Cox regression analysis were carried out to identify the prognostic performance of these factors in the entire cohort. We also made multivariate ROC curves at different observation times to visualize the diagnostic performance of each factor. Then, the correlation between the prognostic model (MDS) and each clinical feature was analyzed to illustrate the reliability of the model we built in the entire cohort. Then, we explored the phenotypic differences between patients in the high-risk group and patients in the low-risk group in the entire cohort. Which includes the expression of risk genes, the methylation level of risk genes, the level of mRNAsi, Tumor mutational burden (TMB), the expression of m6A genes, immunophenoscore (IPS) and the infiltration of immune cells in the entire cohort. Among them, the results of mRNAsi in TCGA-UCEC were obtained from a previous study [25]. In order to obtain more reliable somatic mutation results, “VarScan2” software was used to identify tumor cell mutations. The “maftools” package was used to read and visually analyze the somatic variants of each sample. TMB was defined as the number of somatic variants per megabase of genome. The infiltration of immune cells was obtained through the CIBERSORT tool [26]. Scores of immune and stromal cells were also calculated by Estimation of Stromal and Immune cells in Malignant Tumours using Expression data (ESTIMATE) algorithm. What’more, IPS, which could predict the efficacy of immunotherapy, was calculated as described in the previous research [27].

Quantitative real-time RT-PCR (qRT-PCR) analysis

qRT-PCR was used to further verify the expression of prognostic genes in tissue samples. After obtaining the informed consent of each patient and the approval of the Institutional Review Committee of Nanjing Medical University, we collected 13 pairs of normal tissues and cancer tissues. These tissues were stored at − 80 °C before using.

Total RNA was extracted by using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). And the NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA) was used to measure the concentration of isolated RNA. High-capacity cDNA reverse transcription kits (Thermo Fisher Scientific) were used for reverse transcription of RNA extracted. The program of qRT-PCR was performed with the help of the SYBR Green PCR kit (Thermo Fisher Scientific) on Light Cycler 480 (Roche, Switzerland). The 2^–ΔΔCt method was used to assess the relative expression level of mRNAs. The primer sequence information is in Additional file 13: Table S2.Additional file: As per journal requirements, every additional file must have a corresponding caption. In this regard, please be informed that the caption of Additional file 13 was taken from the additional e-file itself. Please advise if action taken appropriate and amend if necessary.we have checked it.Thanks.

Construction of prognostic nomogram

After obtaining the results of Cox regression analysis in the independent evaluation of prognostic factors, we used the independent prognostic factors as parameters for constructing a nomogram with more accurate prognostic performance. The nomogram was carried out by means of “rms” package. We also made calibration curves and calculated C-index to evaluate the accuracy of the nomogram.

Gene-set enrichment analysis

The genome wide expression profiles of the UCEC in the high risk group and low risk group were employed to identify the biological functions and pathways by Gene-set enrichment analysis (GSEA) with GSEA tool (version 4.0.1) (http://software.broadinstitute.org/gsea/index.jsp) and annotated gene sets (h.all.v7.4.symbols.gmt) as previous research did [28]. The gene sets were filtered using the maximum and minimum gene set size of 500 and 15 genes, respectively. The enriched gene sets were obtained based on Nom.p < 0.05 and FDR < 5% after performing 1000 permutations.

Joint survival analysis of methylation and gene expression

By means of cox regression analysis and K–M methods, joint survival analysis based on gene expression levels and methylation was carried out.

Statistical analysis

Statistical analyses of all data utilized in this article were completed by R software (version 3.6.1, https://www.r-project.org/). When the difference met a joint satisfaction of adjust p < 0.05 and |fold changes (FC)| > 1, it was regarded to be statistically significant. Student’s t test was used for continuous variables, while categorical variables were compared with the chi-square (χ²) test. The Wilcoxon rank-sum test was utilized to compare ranked data with two categories. The Kruskal–Wallis test was utilized for comparisons among three or more groups. The univariate Cox regression analysis and multivariate Cox regression analysis were used to evaluate the relationship between MGs expression and survival data to establish a prognostic model. The Pearson Correlation Coefficient was used to investigate the correlation between MDS and other parameters. “rms” package of R software was used to create the nomogram. The receiver operating characteristic curves were created by the “survivalROC” package of R software and AUC values were also calculated by this package. If the AUC > 0.60, we would consider this model to have a certain predictive ability. If the AUC > 0.75, we would consider this prediction model to have excellent predictive value. All statistical tests were two-sided and p < 0.05 was considered to be statistically significant.

Results

Determination of DMDGs and functional enrichment analysis

The total workflow is shown in the following figure (Additional file 1: Figure S1). A total of 4774 genes were found to be abnormally expressed in tumor tissues, including 2724 up-regulated genes and 2050 down-regulated genes. We used the MethylMix algorithm to integrate the methylation data and gene expression data of each sample. Then, a mixed model including DNA methylation and mRNA expression was established. A Wilcoxon rank test was performed to determine the methylation-driven genes that were significantly differentially expressed. In the end, 48 methylation-driven genes were found to be significantly differentially expressed (|logFC| > 1; p < 0.05; Cor < − 0.3). Among these DMDGs, 12 were hypomethylated genes and the others were hypermethylated genes (Table 1). Figure 1 shows the distribution map of several genes with highest methylation levels. Figure 2 shows the correlation between gene expression and DNA methylation of top hypermethylated and hypomethylated genes. The methylation levels and gene expression levels of DMDGs in each sample we included in the study were displayed in the form of heat maps (Fig. 3A, B).

Table 1 Methylation-driven genes in TCGA-UCEC

Fig. 1

Summary of top methylated genes. The distribution map showing the methylation degree of methylated genes in TCGA-UCEC. The histogram demonstrates the distribution of methylation in the tumor samples (denoted as β-values, where higher β-values represent greater methylation)

Fig. 2

The correlation between gene expression and DNA methylation of top hypermethylated and hypomethylated genes. Average β-values are presented on the x-axis, log2 FPKM gene expression values are presented on y-axis

Fig. 3

Heatmap and functional enrichment analysis of methylation-driven genes in TCGA-UCEC. A Methylation values (β-values). The color change from green to red in the heatmap illustrates the trend from low to high methylation. B The expression of the methylation-regulated genes. The color change from green to red in the heatmap illustrates the trend from low to high expression. C The outer circle represents the methylation values (logFC) of methylated genes in each enriched GO (gene ontology) term: red dots on each GO term indicate upregulated genes, whereas blue dots indicate downregulated genes. The inner circle indicates the significance of GO terms (log10‐adjusted p values). D Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways

In order to explore the possible role of these DMDGs in the occurrence and development of UCEC, GO enrichment analysis and KEGG enrichment analysis of DMDGs were done. From Fig. 3C, we can see that these DMDGs were mainly enriched in eight GO terms, including benzene-containing compound metabolic process, xenobiotic catabolic process, cardiac chamber formation, DNA damage response, signal transduction resulting in transcription, catechol-containing compound biosynthetic process and so on. KEGG enrichment analysis indicated that these DMDGs may be related to Herpes simplex virus 1 infection, Glutathione metabolism, Drug metabolism—cytochrome P450, Platinum drug resistance, Metabolism of xenobiotics by cytochrome P450, Chemical carcinogenesis and so on (Fig. 3D).

Establishment of an independent prognostic index

Firstly, the patients included in the study were randomly divided into training cohort (n = 216) and testing cohort (n = 213) according to the composition of each clinical feature. Secondly, 3 DMDGs (PARVG, SYNE4 and CDO1) with prognostic significance were identified by univariate Cox regression analysis in the training group. The results of Lasso Cox regression analysis showed that the three prognostic DMDGs (PARVG, SYNE4 and CDO1) had no obvious overfitting (Additional file 2: Figure S2A, B). The relationship between the three DMDGs and the OS of UCEC patients was also confirmed in the multivariate Cox regression analysis. Then, we further verified the expression of these three genes in 13 pairs of normal tissues and tumor tissues. We found that SYNE4 and PARVG were significantly higher expressed in tumor tissues, while CDO1 showed no significant difference (Additional file 2: Figure S2C–E).

According to the 3 DMDGs and their corresponding coefficients in multivariate Cox regression analysis, an independent prognostic index (MDS) was established. The formula of MDS was defined as follows: 3.75767726426642 × PARVG expression + 5.91085180306775 × SYNE4 expression + 5.50757519645672 × CDO1 expression.

Evaluation of the prognostic model index

Based on the same formula, we calculated the risk score for each patient. Then, patients were divided into high-risk and low-risk groups according to the median risk score of MDS. After obtaining the prognostic model index for predicting the prognosis of UCEC patients, a series of measures were taken to evaluate the model.

Firstly, the survival status of patients in the high-risk group and the low-risk group is shown in Fig. 4. Figure 4A, B show the result of risk classification of patients in training cohort and in testing cohort according to MDS respectively. From Fig. 4C, D we found that no matter in training cohort or in testing cohort, as the risk score increases, the number of dead patients increases. Kaplan–Meier curves based on the log-rank test were created to visualize the prognostic value of MDS in training cohort and in testing cohort. From Fig. 4E, G, we found that whether in training cohort or testing cohort, patients with high-risk score have a poor prognosis. Figure 4F is a time-dependent ROC curve in training cohort. The AUC values of the ROC curves in 1, 3, and 5 year were all > 0.60. Among them, MDS showed an excellent performance in predicting the 1-year survival rate of UCEC patients. Figure 4H is a time-dependent ROC curve in testing cohort. The AUC values in 1, 3, and 5 year were 0.635, 0.635 and 0.62 respectively.

Fig. 4

Prognostic analysis of the training cohort and the testing cohort. A Risk score distribution of patients in the prognostic model in training cohort. B Risk score distribution of patients in the prognostic model in testing cohort. C Survival status scatter plots for patients in the training cohort. D Survival status scatter plots for patients in the testing cohort. E Kaplan–Meier curve analysis in the training cohort. F Time-dependent ROC curve analysis in the training cohort. G Kaplan–Meier curve analysis in the testing cohort. H Time-dependent ROC curve analysis in the testing cohort

Secondly, we further evaluated the prognostic value of MDS in patients with different clinical characteristics (Additional file 3: Figure S3). According to age, grade, pathological type and stage, we divided all patients into different groups. Kaplan–Meier curves were created in different groups. The prognostic model we constructed showed good prognostic ability in both patients ≤ 60 years and patients > 60 years (p = 0.002 and p = 0.022 respectively, Additional file 3: Figure S3A, B respectively). Additional file 3: Figure S3C, D show the survival curves of patients with grade 1–2 and patients with grade 3–4, respectively. Among them, MDS could better distinguish the survival in patients with grade 3 and grade 4 (p < 0.001). Additional file 3: Figure S3E reveals that MDS had a good prognostic ability in patients with endometrial cancer (p < 0.001), while seemed to be meaningless in other types of patients (p = 0.962, Additional file 3: Figure S3F). In patients with stage 1 to 2, MDS also showed good prognostic value (p = 0.020, Additional file 3: Figure S3G). However, in patients with stage 3 to 4, MDS was unable to distinguish patients with different survival outcomes (p = 0.090, Additional file 3: Figure S3H).

Thirdly, to further evaluate whether our model could be used as an independent prognostic factor, we included some key clinical characteristics containing age, stage, histological type, grade and MDS as independent variables. In training cohort, by means of univariate and multivariate Cox regression analysis, our established MDS remained significant (both P < 0.05, Table 2). While in entire cohort, stage, grade and MDS were regarded as independent prognostic factors.

Table 2 Cox regression analysis data

What’s more, we once again verified the differences of key prognostic DMDGs between the high-risk group and the low-risk group from the gene expression levels and methylation levels in entire TCGA cohort (Additional file 4: Figure S4). Then, the relationship between the prognostic model we constructed (including MDS and 3 risk genes) and clinical characteristics was also revealed (Additional file 5: Figure S5).

Figure 5A–F show the time-dependent ROC curves of the multi-prognostic signatures. Figure 5A, B show the multivariate ROC curves in 1 year. Figure 5D shows the multivariate ROC curves in 3 years. Figure 5E, F show the multivariate ROC curves in 5 years. We found that MDS, stage and grade have similar accuracy in predicting the prognosis of UCEC patients. In addition, it was interesting that after combining MDS and clinical information into a mixed model, more accurate prediction capabilities could be obtained (Fig. 5B, D, F).

Fig. 5

multivariate ROC curve analysis and nomogram with Calibration curves for the prediction of prognosis at 1, 3 and 5 years in the TCGA cohort. A, B ROC curve analysis with multiple variables for predicting 1-year OS. C, D ROC curve analysis with multiple variables for predicting 3-year OS. E, F ROC curve analysis with multiple variables for predicting 5-year OS. G Nomogram for OS. H Calibration curves at 1 year. I Calibration curves at 3 years. J Calibration curves at 5 years

Text: “Figure 5D and Figure 5D” has been changed as “Figure 5D”. Please check and amend if any.we have checked it.Thanks.

Construction of prognostic nomogram

Finally, to better predict the 1-year OS, 3-year OS and 5-year OS of UCEC patients, we constructed a new Nomogram based on the results of the multivariate Cox regression analysis of independent prognostic factors (Fig. 5G). Figure 5H–J show the Calibration curves of the nomogram for the probability of OS at 1, 3 and 5 year. The C-index of the nomogram for OS prediction was 0.764 (95% CI = 0.702–0.826). Additional file 6: Figure S6 shows the time-dependent C-index of MDS, stage, grade, and nomogram.

Gene-set enrichment analysis

To uncover the biological pathways most likely to be related to MDS, we conducted a GSEA analysis (Fig. 6). Among them, E2F TARGETS, G2M CHECKPOINT, MITOTIC SPINDLE, MYC TARGETS V2 and SPERMATOGENESIS were found to be significantly enriched in high-risk group. While in low-risk group, ALLOGRAFT REJECTION, IL2 STAT5 SIGNALING, IL6 JAK STAT3 SIGNALING, INFLAMMATORY RESPONSE and INTERFERON GAMMA RESPONSE were significantly enriched.

Fig. 6

GSEA enrichment analysis. A Gene-set enrichment analysis of genes that are differentially expressed in high-risk group. B Gene-set enrichment analysis of genes that are differentially expressed in low-risk group

TMB, mRNAsi and m6A regulators in two groups

Through GSEA analysis, we found an interesting phenomenon that some immune-related pathways may affect the prognosis of UCEC. Which caused us to think, was the prognosis of UCEC patients related to the difference in the immune microenvironment within the tumor? It has been reported that the stemness of tumors, the methylation modification mediated by m6A regulatory factors, and TMB are related to the immune microenvironment of tumors [29, 30]. Therefore, we decided to further explore whether there were differences between the two groups of patients in these aspects.

From Fig. 7A, B, we could find that patients in the low-risk group have higher TMB level than patients in the high-risk group. There was a significant negative correlation between TMB and MDS. The somatic mutations seen in the two groups are visually displayed in Additional file 7: Figure S7A, B. Additional file 7: Figure S7C suggests that TP53 has a higher mutation frequency in the high-risk group, while ARID1A, CTCF, KMT2B, MUC16 and PTEN have a higher mutation frequency in the low-risk group. In addition, the mRNAsi level of the high-risk group was significantly higher than that of the low-risk group (Fig. 7C). What’s more, there were significant differences in the expression of some m6A regulatory factors between patients in high-risk group and patients in low-risk group (Fig. 7D).

Fig. 7

The correlation between MDS and other phenotypic studies. A, B The correlation between MDS and tumor mutation burden (TMB). C The correlation between MDS and stem cell index (mRNAsi). D The correlation between MDS and gene expression of m6A regulators. Green represents the low-risk group, and red represents the high-risk group

IPS and the expression of immune check point in two groups

Through the ESTIMATE evaluation method, TumorPurity, ImmuneScore and StromalScore were calculated. Figure 8 indicates that patients in the low-risk group have lower TumorPurity and higher ImmuneScore and StromalScore. Figure 8E shows the correlation between MDS and immune cell infiltration. From which, we could find that compared with patients in high-risk group, patients in low-risk group had more T cells CD8, Eosinophils and Neutrophils, while the Macrophages M0 and Macrophages M2 was at low level. Furthermore, the score of macrophage M2 was found to be significantly positively correlated with MDS. While the score of T cells CD8 was found to be significantly negatively correlated with MDS (Additional file 8: Figure S8).

Fig. 8

The correlation between MDS and immune related indicators. A The correlation between MDS and tumor purity. B The correlation between MDS and ESTIMATEScore. C The correlation between MDS and ImmuneScore. D The correlation between MDS and StromalScore. E The correlation between MDS and immune cell infiltration

Then, the ssGSEA algorithm was used to quantify the level of immune infiltration. Comparison of the ssGSEA scores between different risk groups is shown in Additional file 9: Figure S9. From which we found that in addition to score of NK cells and score of type II IFN response, other parameters of immune cell infiltration indicate that patients in the low-risk group have higher immune activity.

To further explore whether MDS can be instructive in guiding immunotherapy, we compared the expression of immune checkpoints and immunophenoscore between two groups (Fig. 9). The results showed that 12 immune checkpoint molecules were highly expressed in patients in the low-risk group (Fig. 9A–C). Furthermore, the gene expression of these immune checkpoints was found to be significantly negatively correlated ted with MDS (Additional file 10: Figure S10). In addition, similar results were obtained in terms of immunogenicity. The IPS, IPS-CTLA4, IPS-PD1-PD-L1-PD-L2 and IPS-PD1-PD-L1-PD-L2-CTLA4 scores were higher in patients in the low-risk group (Fig. 9D).

Fig. 9

The correlation between MDS and gene expression of immune checkpoints. A–C The relationship between MDS and 12 immune checkpoint molecules. D The correlation between MDS and immunophenoscore (IPS)

Joint survival analysis

According to the results of correlation analysis, we finally determined 21 methylation sites (the methylation level of 11 methylation sites is negatively correlated with the expression level of CDO1, and the methylation level of 9 methylation sites is negatively correlated with the expression level of PARVG. and the level of methylation at a methylation site is negatively correlated with the expression level of SYNE4) (Additional file 11: Figure S11). The results of the univariate Cox regression analysis showed that 12 of methylation sites significantly related to CDO1, PARVG and SYNE4 were associated with positive survival outcomes (with all HR > 1) (Fig. 10A). After obtaining mixed models that integrating methylation level and gene expression, we conducted survival analysis on them. The results can be found in Fig. 10B–D. Which shows that hypermethylation and low gene expression PARVG was significantly related to the negative prognosis of UCEC patients. However, the expression of SYNE4 and CDO1 were not related to the prognosis of UCEC patients.

Fig. 10

Joint survival analysis. A Univariate Cox regression analysis assessing the relationship between methylation sites of 3 key MDMGs and overall survival in TCGA-UCEC. B–D Meaningful Kaplan–Meier survival curves for the methylation-driven genes and methylation joint survival analysis

Discussion

Uterus corpus endometrial cancer is one of the main gynecological cancers, its biological, clinical, histological and genetic characteristics are significantly heterogeneous. The traditional classification of endometrial cancer cannot fully reflect the heterogeneity of tumors. Therefore, there are certain limitations in predicting the prognosis of patients and the effect of treatment response. In recent years, more and more studies have provided the basis for prognostic judgment and clinical treatment of UCEC. In 2013, the American Cancer Genome Map Research Center for Endometrial Cancer provided revolutionary insights into the classification of endometrial cancer. According to the integrated genome characteristics, they reclassified endometrial cancer into 4 different types, each with different clinical outcomes. This classification provides more accurate information for clinical treatment and prognosis judgment. What’s more, it is recommended by NCCN guidelines to effectively predict the prognosis of UCEC [31]. Using the miRNA data of endometrial cancer in TCGA, Luo et al. created a prediction model for overall survival and recurrence-free survival, which can effectively predict the prognosis of endometrial cancer patients [32]. However, few evidence has suggested that alterations in DNA methylation could lead to tumor progression and local immunosuppression in the TME.

As we all know, abnormal DNA methylation is a major epigenetic event that promotes the development and progression of cancer. Abnormal DNA methylation is not only related to the shutdown of tumor suppressor genes, but also closely related to the tumor’s malignancy and prognosis [33]. Therefore, our article aims to discover new DNA methylation-driven biomarkers to more accurately predict the prognosis of UCEC, and to explore the possible mechanisms of DNA methylation-driven biomarkers differences between patients with different prognosis. In terms of methods, we used MethylMix to combine the methylation data with the mRNA data from TCGA-UCEC and perform correlation analysis [34]. After that 48 methylation-driven genes were found to be significantly differentially expressed. Then, TCGA-UCEC was randomly divided into training cohort and testing cohort according to the distribution of clinical characteristics, which were used to train and verify the prediction model. We further used univariate Cox regression analysis, multivariate Cox regression analysis and Lasso Cox regression analysis to filer prognosis related genes. 3 DMDGs were finally identified to construct an independent prognostic index. After that, we verified this model in the testing cohort, the results showed that this model could be used as an independent prognostic factor and this model had a good ability to predict the prognosis of UCEC patients. Furthermore, we combined MDS with prognostic clinical indicators to construct a new prognostic nomogram. On the one hand, the results of ROC curve analysis, survival curve analysis and Cox regression analysis all show that MDS can accurately distinguish UCEC patients with different survival outcomes. On the other hand, it is interesting that MDS is positively correlated with some known poor prognostic clinical features (old age, high grade, pathological type with high malignancy and high stage). Therefore, the prognostic model we constructed has a high degree of credibility to a certain extent.

PARVG, SYNE4 and CDO1 are the three risk DNA methylation-driven differential genes we identified. Among them, we verified the difference in expression of PARVG and SYNE4 in 13 pairs of normal tissues and cancer tissues. However, CDO1 was not differentially expressed between 13 pairs of normal tissues and cancer tissues. According to the large sample data showed in Additional file 4: Figure S4, the methylation level of CDO1 is lower mRNA expression of CDO1 is higher in the normal tissue than that in the tumor tissue. However, the range of the mRNA expression is much wide. When the number of the samples was limited like 13 pairs sample in our manuscript. The difference between groups will not be so significant. This might be the reason for the inconsistency between the qRT-PCR results.

PARVG, located in the q13.31 region of chromosome 22, has been proven not to be a tumor suppressor gene related to the occurrence and development of CRC and breast cancer [35]. The results of a methylation sequencing study on adenoid cystic carcinoma found that PARVG in tumor tissues showed lower methylation levels [36]. In addition, a report showed that the high expression of the PARVG gene caused by low-level DNA methylation is significantly related to the poor prognosis of renal cancer patients [37]. In our study, PARVG also showed consistent prognostic performance in UCEC (HR = 42.85, p < 0.05). SYNE4 is a member of the nesprin family of genes, that encode KASH (Klarsicht, Anc-1, Syne Homology) domain-containing proteins. SYNE4 is believed to be involved in the formation of the nucleus skeleton, but its role in tumors is still unknown. As for CDO1, it is the abbreviation of Cysteine Dioxygenase Type 1. CDO1 is a putative tumor suppressor gene in human cancers, and its abnormal DNA methylation level has been studied in a variety of tumors. In gastric cancer, researchers found that CDO1 has significantly DNA hypermethylation than normal tissues. More interestingly, the higher the DNA methylation level of CDO1 in gastric cancer patients, the better the survival prognosis of patients receiving chemotherapy [38]. In lung cancer, researchers also found a similar pattern that the promoter of CDO1 in tumor tissues is hypermethylated. Also, the higher the hypermethylation frequency of CDO1, the lower the degree of tumor differentiation [39]. Which may be related to the limit of the futile metabolism of cysteine and the consumption of cell NADPH caused by CDO1 being silenced, which promotes the proliferation of NSCLC [40]. In pancreatic ductal adenocarcinoma, the hypermethylation of the CDO1 promoter is highly specific in tumor tissues. This feature has been successfully used by researchers as a potential diagnostic biomarker [41]. In our study, the low expression of CDO1 driven by DNA hypermethylation is also significantly related to the poor prognosis of UCEC. Therefore, judging from the genetic makeup of MDS, the predictive prognosis model we constructed is highly reliable.

In addition to constructing a nomogram with high application prospects to accurately predict the prognosis of UCEC, we also explored the possible mechanisms that cause different survival outcomes among UCEC patients through GSEA. Which suggested that there were some changes in immune-related pathways between the high-risk group and the low-risk group. we could find that compared with patients in high-risk group, patients in low-risk group had more T cells CD8, Eosinophils and Neutrophils, while the Macrophages M0 and Macrophages M2 were at low level. Furthermore, the score of macrophages M2 was found to be significantly positively correlated with MDS. While the score of T cells CD8 was found to be significantly negatively correlated with MDS. As early as 2018, Zhang et al. proposed that macrophages induce PDPK1-mediated PGK1 phosphorylation in tumor cells by secreting IL6 to regulate the direction of PGK1 catalytic reaction, thereby enhancing the glycolysis of tumor cells and ultimately promoting the malignant progression of the tumor [42]. T cells CD8+ are the main effector cells that perform antigen-specific killing of tumor cells. Many reports have pointed out that targeted intervention in the glucose metabolism pathway of T cells CD8+ can effectively enhance their anti-tumor function. At the same time, we also explored the differences between the two groups in terms of immunity, RNA methylation modification (m6A), tumor stemness (mRNAsi), and tumor mutation burden (TMB). We found an interesting phenomenon that MDS has a clear relationship with immune-related events. Which was manifested in the significant negative correlation between MDS and TMB (cor = − 0.112, p < 0.05). The high tumor mutation burden (TMB) (TMB ≥ 10 Muts/Mb) has been approved by the US Food and Drug Administration (FDA) to guide biomarkers for immunotherapy. TMB was defined as the number of somatic variants per megabase of genome. These somatic mutations potentially result in new or fragmented proteins/peptides, which can be recognized as exogenous, and elicit an antitumor immune response. TMB can predict the efficacy of immunotherapy outcomes, and filtered candidate patients for immunotherapy treatment. However, in clinical practice, the TMB was assessed by whole-genome sequencing (WGS), whole-exome sequencing (WES) and targeted panel sequencing. Due to the high cost, relatively long turnaround time and the need for sufficient tissue samples. Evaluation of TME in clinical practice is not feasible. In this research, we found the closely relationship between MDS and TMB. With the help of this MDS, we can indirectly reflect the TMB status and provide valuable information on the results of immunotherapy response. The ssGSEA algorithm was also used by us to quantify the condition of immune infiltration. Interestingly, in addition to score of NK cells and score of type II IFN response, other parameters of immune cell infiltration indicate that patients in the low-risk group have higher immune activity. In addition, the expression of some immune checkpoints in patients with high MDS was significantly lower than that of patients with low MDS. The low expression of some immune checkpoints in high-risk group may cause by the low fraction of CD8+ T cells. Which was consist with the TMB results. Therefore, the MDS we constructed can not only accurately predict the prognosis of UCEC patients, but also accurately reflect their tumor immune microenvironment. What’s more, patients in the high-risk group had lower IPS. Which means that the results of our research also provide a reference method for guiding UCEC’s immunotherapy.

However, our study still has some limitations: the results of our study were only validated in the TCGA database. The prognostic model requires more data support from clinical patients. In addition, the mechanism by which DNA methylation-driven genes affect the prognosis of patients with UCEC needs to be further explored through in vivo and in vitro experiments.

Conclusion

All in all, we have successfully established a risk model (MDS) based on 3 DMDGs, which could accurately predict the prognosis of patients with UCEC. A nomogram combining clinical variables and MDS was also drawn to improve the accuracy of the prediction. In addition, our risk score model could connect with Methylation-Driven genes and tumor immune microenvironment, which provides a comprehensive perspective for clarifying the underlying mechanisms that determine the prognosis for UCEC. Our research may shed new light on UCEC patient’ prognosis and treatment management. However, further experiments are also required to validate our findings.