Development and validation of m6A regulators' prognostic significance for endometrial cancer.

BACKGROUND: Endometrial cancer (EC) is the sixth most common cancer in women globally. It has been found that the expression levels of m6A regulators can be potentially used for prognostic stratification in some cancers, but the role of m6A regulators in EC prognosis remains unclear.
METHODS: The data of 584 EC samples were downloaded from The Cancer Genome Atlas and the mRNA expression profiles of 20 m6A regulators were analyzed, followed by functional enrichment analysis, immune infiltration analysis, and least absolute shrinkage and selection operator method-COX regression analysis.
RESULTS: The mRNA expression levels of 20 m6A regulators were significantly different between cancer samples across different grades. The 548 EC samples could be clearly divided into 2 clusters. Kaplan-Meier survival analysis proved that these two groups had highly different overall survival probabilities. Besides, the univariate regression analysis further reserved eight genes related to overall survival from the 20 m6A regulators. We established a prognostic signature including two genes, that is, IGF2BP1 and YTHDF3, that showed a strong ability for stratifying prognostically different EC patients. We identified 3239 differentially expressed genes between the high- and low-risk groups, involving in multiple biological processes and signaling pathways. Meanwhile, 6 differentially infiltrated immune cell types between the high- and low-risk groups could effectively distinguish the high- and low-risk EC groups. The expressions of immune checkpoints were different between high- and low-risk EC patients.
CONCLUSION: We first report the prognostic role of m6A regulators in EC, which should contribute to a better understanding of the underlying mechanisms of EC pathogenesis and progression.

Introduction

Endometrial cancer (EC) is the most common gynecological malignant disease and the sixth most common cancer in women globally.[ The incidence and disease-related mortality of EC continue to increase.[ Although early-stage EC patients have a favorable 5-year relative survival rate (96%), the rate is only 18% in patients with distant metastases.[ Unfortunately, outcomes for EC patients with systemic recurrence are horrible, with a median survival hardly exceeding 12 months.[ Furthermore, the investigation of EC has been lagging behind other cancers with only two known signals, leading to restricted biomarkers and targets.[ EC is one of the few human malignancies whose mortality is increasing, which underscores the urgent need to explore more predicted models for diagnosis and prognosis and develop effective treatment strategies for this disease.[

Emerging evidence has revealed that epigenetic regulation participates in the initiation and progression of multiple malignancies.[ And the epigenetic marks potentially serve as diagnostic, prognostic, and predictive biomarkers for cancer.[ Presently, with the deeper knowledge of epigenetics and the development of epigenome technology, >170 RNA modifications are discovered, and most RNA species contain 1 or multiple distinct chemical modifications and are widely linked to physiology and pathology.[ Several mRNA modification forms have been reported, including N6–2’-O-dimethyladenosine (m6A), inosine, pseudouridine, 5-methylcytidine (m5C), 5-hydroxymethylcytidine (hm5C), and N1-methyladenosine (m1A).[ Among these modification forms, the m6A modification is one of the most abundant internal modulations in mRNA.[ The effects of m6A on mRNA are mediated by expanding m6A regulators. The modification is installed by the m6A methyltransferases such as METTL3/14, WTAP, RBM15/15B, and KIAA1429 termed as “writers,” reverted by demethylases such as FTO and ALKBH5 termed as “erasers,” and recognized by m6A binding proteins such as YTHDF1/2/3, IGF2BP1 and HNRNPA2B1 termed as “readers.”[ Widespread genetic remodeling of m6A regulators and their expression levels are significantly correlated with the activity of cancer hallmark-related pathways and can be potentially useful for prognostic stratification.[ However, the comprehensive landscape of the expression of m6A regulators in EC and their roles in diagnosis and prognosis remain elusive.

Integrating multiple biomarkers into a single model would substantially improve the prognostic value.[ The least absolute shrinkage and selection operator method (LASSO) is a prevalent method for the regression of high-dimensional predictors.[ LASSO has been extended and broadly applied to the Cox proportional hazard regression model for survival analysis of high-dimensional data.[ Meanwhile, The Cancer Genome Atlas (TCGA) database contains a large amount of clinical, pathological, and biological data from patients with malignancies.[ We can more accurately predict the development tendency and dig deeper into the mechanism of tumors, providing a reliable research direction for treatment programs by a comprehensive analysis of the data.[ However, the application of the LASSO Cox model based on the TCGA database in assessing EC prognosis is still limited.

In this study, we aim to systematically analyze the expressions of 20 critical m6A regulators in 548 EC samples from TCGA database. We provided the expression information of each m6A regulator with different grade features. We revealed that the expression of m6A regulators was significantly associated with clinicopathological features and prognosis of EC, and a signature with 2 selected m6A regulators was designed to scale the prognosis of EC.

Materials and methods

Study population

The mRNA expression profiles of 584 EC samples were downloaded from TCGA, which consisted of 548 cancer and 36 adjacent samples. Among the cancer samples, 544 samples had complete survival information and were used for the analysis. Table 1 showed the detailed clinicopathological features of those samples.

Table 1

Clinicopathological characteristics of EC samples from TCGA database.

	EC Patients
Parameters	Training cohort (N = 416)	Validation cohort (N = 128)	χ2	P
Age (mean ± SD)	57.73 ± 9.62	56.89 ± 9.98	0.006156	.9375
Sex
Female	416 (100%)	128 (100%)	0	1
Pathologic stage
I	260 (62.50%)	78 (60.94%)	0.71838	.8689
ii	41 (9.86%)	10 (7.81%)
iii	95 (22.84%)	31 (24.22%)
iv	20 (4.80%)	9 (7.03%)
Grade
1	70 (16.83%)	28 (21.88%)	1.4661	.6901
2	95 (22.84%)	24 (18.75%)
3	246 (59.13%)	73 (57.03%)
High Grade	5 (1.20%)	3 (2.34%)
OS status
Dead	69 (16.59%)	36 (28.13%)	3.1996	.07366
Alive	347 (83.41%)	92 (71.87%)

OS = overall survival, TCGA = The Cancer Genome Atlas.

M6A regulators

Genes studied in this research were 20 m6A regulators collected from the study of Li et al,[ which were classified into 3 types according to their functions: RNA methyltransferases which added methyl groups to RNA (METTL3, METTL14, WTAP, VIRMA, RBM15, RBM15B, and ZC3H13), m6A binding proteins that identified and bound to m6A (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, HNRNPC, and HNRNPA2B1), and RNA demethylases which removed methyl groups from RNA (TOAL and KBS).

Analysis of immune infiltration

The immune infiltration difference of 22 immune cells in the samples was analyzed by using CIBERSORT software combined with the LM22 feature matrix.[ The sum of the proportions of all estimated immune cell types in each sample was equal to 1.

Differential gene expression analysis

Differentially expressed gene analysis was performed based on the limma[ in R language (version3.5.2). The |Log2FC| >1 and FDR ≤0.05 were served as the criteria to screen differentially expressed genes.

Functional enrichment analysis

The clusterProfiler package in R language was used to perform Gene ontology (GO, including Biological Process, Molecular Function, and Cellular Component) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis.[ The significantly enriched GO terms and KEGG pathways were screened using P value <.05 adjusted by the “BH” method as the threshold.

Construction of m6A regulator-based prognostic model

The 544 EC samples with complete survival information were randomly divided into training and validation sets with four-fifths and one-fifth of all samples, respectively. We used the univariate Cox regression method to analyze the association between these 20 m6A regulators and EC patients’ overall survival (OS) in the training set. The “glmnet” Bioconductor package in R was used to conduct LASSO-COX regression analysis for constructing EC prognostic model based on the significant regulators in univariate Cox analysis as follows:

Coef is the risk coefficient of each factor calculated by the LASSO-COX model, and X is the mRNA expression value.

Statistical analysis

We used the Wilcoxon rank test and the analysis of variance to determine the significance of gene expression differences between 2 groups and among multiple groups, respectively. Kaplan-Meier method was used for the estimation of OS probability of EC samples. Log-rank test was applied to the comparison of OS between different groups. P value <.05 was used as the significance threshold. All statistical analyses were performed in R software (version 3.4.1).

Results

Expression of m6A regulators was closely related to EC's grade

Figure 1A illustrated the mRNA expressions of the 20 m6A regulators across 548 EC tumor samples stratified by their grade as a heatmap, which showed that the expression levels of these regulators were variable in samples with different grades. Analysis of variance indicated significant differences of mRNA expressions of all 20 m6A regulators among EC samples with different grades (Figure 1B–F); specifically, almost all of those regulators exhibited gradually increasing expression with the increase of grade. The results indicated that m6A regulatory factors were highly correlated with the grade of EC patients.

Figure 1

The expression of m6A regulator is related to the clinicopathological features of EC. (A) Heat map of expression levels of m6A regulators in 548 samples across grade I to grade III. Red, upregulated. Blue, downregulated. x-Axis and y-axis refer to samples and genes, respectively. The box plots of expression levels of IGF2BP1, YTHDC1, YTHDC2, and IGF2BP2 (B); IGF2BP3, YTHDF1, YTHDF3, and HNRNPA2B1 (C); RBMX, HNRNPC, METTL3, and YTHDF2 (D); METTL14, WTAP, VIRMA and RBM15 (E); ALKBH5, FTO, RBM15B, and ZC3H13 (F) in EC samples with different grades. EC = endometrial cancer.

M6A regulators could separate EC samples with distinct prognosis

We performed clustering analysis for the EC tumor samples based on the mRNA expression levels of the 20 m6A regulators. Consensus clustering analysis using ConsensusClusterPlus Bioconductor package obtained the optimal number of sample clusters as two (Fig. 2A). Figure 2B illustrated the clustering result of the 548 EC tumor samples using the K-means method. To analyze the OS of these 2 groups, we used the survival function package to draw Kaplan-Meier survival curves (Fig. 2C), which showed a significant difference in the OS between the 2 groups (P < .05). Besides, according to the expression levels of all genes in different groups, the principal component analysis (PCA) could also distinctly distinguish samples within the 2 groups (Fig. 2D).

Figure 2

M6A regulators could separate EC samples with distinct prognosis. (A) Determination of the optimal number of consistent clusters. When the curve changes from steep to gentle, the corresponding coordinate on the horizontal axis is the optimal number of consistent clusters. (B) Clustering of EC samples. (C) The survival curve of EC samples based on Kaplan Meier method. The horizontal axis represents time in days, the vertical axis represents survival rate, and the color represents different groups. P value is calculated using log-rank test. (D) The PCA. The points of different colors represent samples of different groups. The distance between the points is closer; the expression of m6A regulators between 2 samples is more similar. EC = endometrial cancer, PCA = principal component analysis.

Prognostic value of m6A regulators in EC

Univariate Cox regression analysis identified IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF3, RBMX, VIRMA, and RBM15 as significant genes whose mRNA expression levels were closely associated with the OS of EC patients (Fig. 3). Besides, the high expression of all those 8 genes was correlated with inferior OS of EC patients according to their hazard ratio (>1).

Figure 3

The prognosis model constructed based on m6A regulators could well predict the OS probability of EC patients. (A) The forest plot of univariate Cox regression analysis based on 20 m6A regulators. HR represents the hazard ratio; 95% CI is the 95% confidence interval. (B) Determination of the optimal gene number based on LASSO Cox regression analysis. The horizontal axis is log(lambda), and the vertical axis is partial likelihood deviance. The optimal lambda value corresponds to the minimum value of partial likelihood deviance, and the optimal gene number is corresponded on the top of the optimal lambda. The survival curve of EC samples in the training set (C) and (D) validation set based on Kaplan Meier method. (E) The ROC curve analysis of the EC patients in the training set for 1-, 3- and 5-year survival. (F) The ROC curve analysis of the EC patients in the validation set for 1-, 3- and 5-year survival. EC = endometrial cancer, LASSO = least absolute shrinkage and selection operator method, ROC = receiver-operating characteristic.

Next, we used the training set to build a risk score model based on these 8 prognosis-related genes through the LASSO Cox regression algorithm and determined the best lambda value by cross-validation. As a result, 2 genes were selected for the model, that is, IGF2BP1 and YTHDF3 (Fig. 3B), and the model was constructed with the following equation: Risk score = 0.0904 × mRNA level of IGF2BP1 + 0.195 × mRNA level of YTHDF3. EC patients in the training set were assigned to a risk score and divided into low-risk and high-risk groups by the median risk score. Kaplan-Meier analysis showed that the prognosis of samples in the high-risk group was significantly worse than that in the low-risk group (Fig. 3C). The same grouping and calculation methods were applied to the samples of the validation set, which obtained consistent results with the training set (Fig. 3D). Moreover, the time-dependent receiver-operating characteristic curve analysis was performed. The area under curve (AUC) values of the training set for 1-, 3- and 5-year survival were 0.6552, 0.6408, and 0.6439 (Fig. 3E), and for the validation set, the corresponding AUC values were 0.8268, 0.8267, and 0.9062 (Fig. 3F), respectively. These results suggested that the established model could efficiently predict the prognosis of EC patients both in the training set and validation set.

Risk score was an independent prognostic factor for EC

We analyzed the risk score of the samples in the training set and validation set, which showed significant difference in the risk score among samples with different grades and stages (Fig. 4A and B). With the increase of Grade and Stage, the risk score showed an upward trend (Fig. 4D and E).

Figure 4

Risk score was an independent factor to predict the OS probability of EC patients. The box plots of risk scores of samples grouped by grade in training set (A) and validation set (D). The box plots of risk scores of samples grouped by stage in the training set (B) and validation set (E). The forest plots of multivariate Cox regression analysis in the training set (C) and the validation set (F). EC = endometrial cancer.

Next, we used the survival package in the R language incorporating grade and stage information to perform multivariate Cox regression analysis on the training and validation sets and determine whether the risk score was an independent prognostic indicator. It was found that the risk score was still significantly correlated with the OS of EC samples in the training and validation sets after adopting these factors into the multivariate Cox regression (Fig. 4C). More importantly, the risk of death was lower for samples with lower risk score compared to those with higher risk score (Fig. 4F). In summary, the risk score can independently predict the prognosis and clinicopathological characteristics of patients with EC.

Differential gene expression and functional enrichment analyses of EC patients in high- and low-risk groups

To further investigate the difference of gene expression between the high- and low-risk groups, and the potential molecular functions or biological processes through which these differentially expressed genes affected the prognosis of EC patients, we conducted differential gene expression and functional enrichment analyses of EC patients in the high- and low-risk groups. Significantly, a total of 3239 differentially expressed genes were identified in the high-risk group compared with the low-risk group, including 3121 upregulated genes and 118 downregulated genes (Fig. 5A). The result showed that the expression levels of differentially expressed genes were remarkably different between the high- and low-risk groups (Fig. 5B).

Figure 5

Differential gene expression and functional enrichment analyses of EC patients in high- and low-risk groups. (A and B). The differential gene analysis was performed using limma in R. (A) The differentially expressed genes were shown in the volcano plot. The red dots represented the down-regulated genes and the blue dots represented the up-regulated genes. (B) The differentially expressed genes were presented in the heatmap. (C and D) The GO and KEGG enrichment analysis based on the differentially expressed genes were conducted using clusterProfiler package in R. The top 20 significant GO terms (C) and the 3 significant KEGG pathways (D) were shown. EC = endometrial cancer, KEGG = Kyoto Encyclopedia of Genes and Genomes.

Then, we performed GO and KEGG enrichment analysis based on these 3239 genes, which were enriched in many GO terms, such as chromosome segregation, and 3 KEGG pathways, such as cell cycle. The top 20 GO terms and the 3 KEGG pathways were shown in Figures 5C and D, respectively. The detailed results of GO and KEGG enrichment analyses were shown in Table S1 (see Table, Supplemental Digital Content, functional enrichment analysis results).

Immune landscape of the low- and high-risk EC patients

Next, the immune infiltration difference of 22 immune cells in the high- and low-risk groups of EC patients was analyzed using CIBERSORT software combined with the LM22 feature matrix. The difference of immune cell infiltration of 544 EC patients was shown in Figure 6A. In addition, the correlation between the infiltration ratios of different immune cell types was weak (Fig. 6B), indicating that there was a large heterogeneity of the infiltration of different immune cells in EC patients. Meanwhile, 6 immune cell types presented significant infiltration difference between the high- and low-risk groups (Fig. 6C). Compared with the low-risk group, the activated dendritic cells, resting mast cells, and follicular helper T cells showed higher proportions of infiltration, whereas the monocytes, CD8 T cells, and regulatory T cells (Tregs) exhibited lower proportions of infiltration in the high-risk group, implying a potential association of these immune cells with EC patients’ prognosis. PCA analysis based on these 6 significantly different immune cell types could effectively divide the EC samples into high- and low-risk groups (Fig. 6D), suggesting that the infiltration difference of immune cells might be potentially correlated with EC patients’ prognosis.

Figure 6

Immune landscape of the low- and high-risk EC patients. (A) The immune infiltration difference of 22 immune cell types in high- and low-risk EC patients. (B) The correlation matrix of 22 immune cell infiltration. Red represented a positive correlation and blue represented a negative correlation. The darker the color, the greater the correlation. (C) The infiltration difference of 6 immune cell types was analyzed between the low-risk set and the high-risk set. (D) PCA analysis based on the 6 significantly different immune cell types between the high- and low-risk groups. EC = endometrial cancer, PCA = principal component analysis.

The expression of immune checkpoints has become a biomarker of immunotherapy for patients with EC.[ Significantly, we identified that the risk score was closely correlated with the key immune checkpoints, including CTLA4, PDL1, IDO1, TDO2, LAG3, and TIGIT (Fig. 7A). Meanwhile, the expressions of PDL1, IDO1, TDO2, LAG3, and TIGIT were significantly higher in the high-risk group than those in the low-risk group (Fig. 7B–F), implying that patients in the high-risk group might be more sensitive to the treatment of immunosuppressive drugs.

Figure 7

The correlation of the risk score with immune checkpoints in EC patients. (A) The chord diagram illustrating the correlation of the Risk Score with the expression of the 6 key immune checkpoints, including CTLA4, PDL1, IDO1, TDO2, LAG3, and TIGIT. (B–F) The expressions of immune checkpoints PDL1, LAG3, TIGIT, IDO1, and TDO2 in the low-risk set and high-risk set of EC patients. EC = endometrial cancer.

Discussion

EC carries a significant risk of systemic and locoregional recurrence.[ Although patients diagnosed and treated at an early stage possess relatively good survival rates, those women who are diagnosed with either advanced-stage disease or suffer a recurrence have a poor prognosis.[ M6A mRNA modification is involved in tumorigenesis and m6A regulators including m6A writers, erasers, and readers play key roles in cancer diagnosis, prognosis, and therapy.[ However, the effect of m6A regulators on EC development and prognosis remains unclear. In this study, we demonstrated that the expression of m6A regulators was closely associated with the malignancy and prognosis of EC. We identified 2 EC clusters by consensus clustering based on the expression of m6A regulators. These 2 subgroups were closely correlated with the prognosis and clinicopathological features of EC. Remarkably, we established a prognostic signature with two selected m6A regulators, which scaled the OS of EC patients into high- and low-risk categories.

It has been reported that m6A RNA methylation regulators affect multiple pathological processes in cancer progression. As an m6A writer, METTL3 is elevated in multiple cancers and promotes tumor proliferation.[ The reader IGF2BP1 promotes SRF-dependent transcription in cancer in an m6A- and miRNA-dependent manner.[ YTHDF1 regulates hypoxia adaptation and non-small cell lung cancer progression.[ The eraser ALKBH5 maintains the tumorigenicity of glioblastoma and promotes cell proliferation program.[ These findings suggest that the differential expression of specific m6A is linked to dysregulation of RNA in tumors. Considering the crucial biological functions of m6A regulators in tumorigenesis, we systematically investigated the relationships between the m6A RNA methylation regulator in each individual and the pathological features of EC. Interestingly, we observed that the expressions of most m6A RNA methylation regulators were significantly changed in EC samples with different grades. It provides new evidence that m6A regulators play a crucial role in EC progression, which is consistent with previous studies. Based on the expression similarity of m6A RNA methylation regulators in 548 EC samples, we identified 2 subgroups by consensus clustering. Surprisingly, we observed a significantly different OS in the 2 subgroups. And the PCA results also showed a clear distinction between them by comparing the expression levels. These data suggest that the expression of m6A RNA methylation regulators is closely associated with the clinicopathological features of EC. Moreover, these findings are also benefiting to develop novel therapeutic strategies by characterizing the expression of each individual m6A methylation regulator in EC, since targeting m6A methylation is considered as a new method for cancer therapy.[

Precision medicine needs accurate prognostic prediction.[ The LASSO method is a popular method for the regression of high-dimensional predictors.[ The LASSO models are widely used in the prognosis prediction of cancer.[ The previous study integrated the TCGA data analysis and the LASSO algorithm also showed that m6A regulators contributed to malignant progression and had clinical prognostic impact.[ Besides, LASSO models are used to identify the molecular biomarkers associated with EC progression and prognosis.[ In this study, we also sought to explore the prognostic value of m6A regulators in EC. The results of a univariate Cox regression analysis identified that 8 of 20 tested genes were significantly correlated with OS, in which IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF3, RBMX, VIRMA, and RBM15 were risky genes with HR >1. The previous study showed that IGF2BP1 and YTHDF1 could enhance ovarian cancer phenotype.[ IGF2BP2 serves as a potential prognostic biomarker of colorectal carcinoma.[ IGF2BP3 functions as a potential oncogene for many types of cancer and exhibits prognostic significance in endometrial clear cell carcinomas.[ The writer/reader VIRMA/YTHDF3 is upregulated in seminoma, constituting novel candidate biomarkers for patient management.[ RBMX is involved in the programmed cell death of breast cancer.[ RBM15 plays a crucial role in megakaryocytic leukemia.[ Our data further imply that these m6A regulators may serve as potential biomarkers and therapeutic targets of EC. To better predict the clinical outcomes of EC with m6A RNA methylation regulators, we applied the LASSO Cox regression algorithm to the 8 prognosis-associated regulators, in which IGF2BP1 and YTHDF3 were selected to build the risk signature and calculate the risk score by the coefficients obtained from the LASSO algorithm. Significantly, we observed a notable difference in OS between the 2 categories of TCGA database-derived EC patients that were separated into low- and high-risk groups based on the median risk score. It indicates that the consistency between our LASSO Cox model based on the m6A regulators and the clinical prognosis. Additionally, in the time-dependent receiver-operating characteristic curve analysis, the AUC values of the training and validation sets for 1-, 3-, and 5-year survival were relatively high, indicating that the established model could efficiently predict the prognosis of EC patients. Strikingly, we also observed that the risk scores of the above model were positively related to the clinical grade and stage of EC. These results indicate that the risk scores calculated by the signature can potentially predict EC patient outcomes and clinicopathological features. The multivariate Cox regression analysis was further performed to determine whether the risk signature was an independent prognostic indicator. Importantly, our results confirmed that the risk score originated from the m6A regulator could independently predict the prognosis of EC patients.

The differential gene expression is closely correlated with the prognosis of EC patients. For example, it has been reported that Akt1 is abnormally expressed in EC samples and associated with the prognosis of EC patients.[ The elevation of BP1 contributes to the malignant progression of EC and predicts poor prognosis of EC patients.[ Moreover, the progression of EC is complicated and multiple biological processes and signaling pathways are involved in the EC pathogenesis, such as chromosome segregation and cell cycle.[ In this study, a total of 3239 differentially expressed genes were found in the high-risk group compared with the low-risk group, which participated in several GO terms and KEGG signaling pathways, including chromosome segregation and cell cycle. These data are consistent with previous reports and provide valuable evidence for the association of abnormal gene expression with EC prognosis. Moreover, immune infiltration and immune checkpoints play crucial roles in the modulation of EC development and serve as potential targets for EC treatment.[ In the present study, we found that the infiltration of 6 immune cell types was significantly different between the high- and low-risk groups. Compared with the low-risk group, the activated dendritic cells, resting mast cells, and follicular helper T cells showed higher proportions of infiltration, while the monocytes, CD8 T cells, and regulatory T cells (Tregs) exhibited lower proportions of infiltration in the high risk group. Moreover, the clustering analysis based on these 6 immune cell types reliably distinguished the high- and low-risk groups. These results indicate that the difference in immune infiltration is significantly associated with the prognosis of EC patients. Dendritic cell infiltration is an independent factor related to the inferior prognosis of melanoma patients,[ and elevated CD8 T cell infiltration is associated with prolonged OS in hepatocellular carcinoma.[ These results are consistent with our study that the EC patients with high risk score have high infiltration proportion of activated dendritic cells, low proportion of infiltrating CD8 T cells and inferior prognosis. However, the infiltration of the remaining 4 immune cell types and their associations with prognosis in EC show distinct tendencies compared with other cancer types. For example, high infiltration proportion of resting mast cells presents a prognostic value for superior survival outcome in lung adenocarcinoma.[ Presence of follicular helper T cells is considered as a predictor of improved prognosis in breast carcinoma.[ Decreased regulatory T cells (Tregs) are related to improved prognosis in some cancers.[ Although there are no significant relationship between monocytes infiltration and the prognosis of renal cell carcinoma and breast cancer.[ These differences may be caused by the distinct functions and underlying mechanisms of the infiltrating immune cells in EC, which still needs further investigation.

Meanwhile, the expressions of immune checkpoints PDL1, IDO1, TDO2, LAG3, and TIGIT were obviously increased in the high-risk group compared with the low-risk group. Considering that the EC patients with high risk score had worse survival outcome than those with low risk score, our study also indicated a potential correlation of the increased expression of key immune checkpoints with the poor prognosis of cancer patients. These immune checkpoints and their association with the prognosis of cancer patients have been previously investigated. PDL1 plays a key role in the immune escape regulation of tumor cells, which shows obvious upregulation in multiple cancers, such as colorectal cancer, gastric cancer, non-small cell lung cancer, and breast cancer, and is associated with poor prognosis in patients with cancer.[ IDO1 is a crucial regulator of immunosuppression during the progression of cancer, and increased IDO1 expression is correlated with inferior survival outcome in glioblastoma.[ TDO2 is implicated in immune microenvironment and multiple processes associated with immune accommodation. Compared with the adjacent normal tissues, TDO2 is highly expressed in various cancer tissues, and its overexpression is associated with poor survival outcome in breast cancer.[ In patients with renal clear cell carcinoma and renal papillary cell carcinoma, a positive correlation is observed between the elevated LAG3 expression and poor prognosis.[ Similarly, gastric cancer patients with overexpressed TIGIT exhibit inferior prognosis.[ Our results are consistent with these researches. Furthermore, our study also implies that the EC patients in the high-risk group may be more available for anti-immune checkpoint therapy, which is of great clinical concern.

In conclusion, our finding is the first demonstration that systematically depicted the expression and prognostic value of m6A RNA methylation regulators in EC based on the TCGA database and LASSO Cox model. The expressions of m6A regulators are highly associated with the poor clinical features and prognosis of EC. Our study indicates new insights and basic evidence for future investigation on the role of m6A modification in EC, providing the potential of m6A regulators as biomarkers and therapeutic targets for EC.

Author contributions

Xuecheng Pang and Sumin Qian put forward the ideas of this article, written this article and analyzed the data. Xuecheng Pang, Xiang Zhang, Yue Huang and Sumin Qian helped for acquisition of data and analysis and interpretation of data. Xiang Zhang and Yue Huang helped for revising the manuscript. All authors read and approved the final manuscript.

Data curation: Xuecheng Pang, Sumin Qian.

Methodology: Xuecheng Pang, Sumin Qian.

Resources: Sumin Qian.

Software: Xuecheng Pang, Xiang Zhang, Yue Huang, Sumin Qian.

Supervision: Xuecheng Pang, Xiang Zhang, Yue Huang, Sumin Qian.

Validation: Xuecheng Pang, Xiang Zhang, Yue Huang, Sumin Qian.

Visualization: Xuecheng Pang, Sumin Qian.

Writing – original draft: Xuecheng Pang, Sumin Qian.

Writing – review & editing: Xuecheng Pang, Sumin Qian.