Pan‐cancer analysis of the expression of m6A regulators
m6A RNA modification is dynamically regulated by methylases and demethylases, whereas readers mediate m6A-dependent functions. In this study, we selected seven mRNA m6A regulators including three “writers” (METTL3, WTAP, and ZC3H13), two “erasers” (FTO and ALKBH5) and two “readers” (YTHDF1 and YTHDF2). The overall expression levels of these seven m6A regulators in pan-cancer are shown in Figure 1A. Our results indicate that ALKBH5 had the highest expression in pan-cancers, while FTO expression was lowest. Further analysis revealed that YTHDF1 had the highest expression in CHOL and READ, while ZC3H13 had the highest expression in STAD, WTAP was expressed at the highest levels in GBM, and FTO was most strongly expressed in CHOL (Fig. 1B). Correlations among the seven m6a regulators are presented in Figure 1C.
We also analyzed the differential expression of these m6A regulators between tumor and adjacent normal tissues in 18 TCGA cancer datasets (cancers with > 5 normal tissue samples were selected for this analysis). Our results show that the expression of ALKBH5 did not differ in ESCA, KIPP, READ, or STAD compared with adjacent normal tissue, but was clearly expressed at higher levels in CHOL, HNSC, COAD, KIRC, LIHC, BLCA, and BRCA, and expressed at significantly lower levels in GBM, KICH, LUAD, LUSC, PRAD, THCA, UCEC (Fig. 2D). The expression levels of the other six m6A regulators in pan-cancer are presented in Figure 2E–J.
Prognostic value of m6A regulators in pan-cancer
Subsequently, we evaluated the prognostic value for overall survival (OS) of the seven m6A regulators in pan-cancer. K-M curves showed that m6A regulators exhibited striking prognostic heterogeneity among cancers. Compared with the low expression group, the ALKBH5 high expression group had worse OS in BLCA and LAMA, while in ESCA, OV, PAAD, high ALKBH5 expression was related to good OS. High FTO expression was associated with worse OS in BLCA and THCA, while it was related to better OS in KIRC and LGG. The prognostic data for all seven m6A regulators in pan-cancer are presented in Figure 2. Overall, these data suggest that m6A regulators are key factors affecting the prognosis of patients with cancer, with differing roles across various cancers, depending on the type of tumor.
Correlations of m6A regulator expression with TMB and MSI in various cancers
We calculated TMB and MSI values for each pan-cancer patient and analyzed correlations of these values with the expression levels of the seven m6A regulators. Our data showed that ALKBH5 levels are positively correlated with TMB in COAD, UCEC, SARC, and LGG, and negatively correlated with TMB in KIRC, BRCA, THYM, THCA, PRAD, PAAD, ESCA, and DLBC (Fig. 3A). Correlations of the other six m6A regulators with TMB in pan-cancer are presented in Figure 3B–G.
ALKBH5 levels were positively correlated with MSI in COAD, UCS, UCEC, STAD, SARC, PCPG, and LUAD and negatively correlated with MSI in THCA and LAML (Fig. 3H). Correlations between the other six m6A regulators and MSI in pan-cancer are shown in Figure 3I–N.
These results demonstrate that the correlations of m6A regulators with TMB and MSI across various cancers are heterogeneous.
Correlations of m6A regulators with tumor stemness
Stemness indices are used to describe the proportions of tumor cells resembling stem cells. Here, we evaluated the associations of seven m6A regulators with stemness scores (RNAss and DNAss) (Fig. 4A–B). Expression levels of YTHDF1 and YTHDF2 were positively correlated with RNAss and DNAss in the majority of 33 TCGA tumor types tested, with higher expression levels associated with greater tumor stemness scores. In contrast, FTO and ZC3H13 expression levels were negatively correlated with tumor stemness scores in most TCGA tumors.
Correlation of m6A regulator expression with immune subtype and immune genes
Correlation coefficients between immune characteristics have been applied to cluster non-hematological tumors in TCGA into six immune subtypes (C1–C6). Here, we analyzed potential correlations between the expression levels of seven m6A regulators and immune subtypes in pan-cancer. The expression levels of the seven m6A regulators in pan-cancer data C1–C6 immune subtypes are shown in Figure 4C, and demonstrate that expression levels of the m6A regulators differed significantly among immune subtypes. METTL3, ZC3H13, FTO, and YTHDF1 were highly expressed in C5 subtype tumors; however, ALKBH5, YTHDF2, and WTAP were expressed at lower levels in C5 subtype.
Immune genes, particularly immune checkpoint genes, have a major impact on tumor immune infiltration, immune escape, and response to immunotherapy [26, 27]. Exploration of associations between expression of the seven m6A regulators and that of immune genes was conducted to better comprehend the potential influence of m6A regulators on immunotherapy. We evaluated the associations of seven m6A regulators with 47 immune genes in pan-gynecologic cancers. The results showed that expression levels of m6A regulators were mostly positively correlated with those of immune genes (Fig. 4D). In particular, in UCEC 31 of 47 immune genes were associated with WTAP expression. These findings suggest that these m6A regulators may interact with immune genes in various signal transduction pathways, and may represent ideal future immunotherapy targets.
Correlations of expression of m6A regulators with immune, stromal, and ESTIMATE scores
The ESTIMATE algorithm can use expression profile data to estimate the proportions of stromal and immune cells in malignant tumor tissues, and calculate corresponding stromal and immune scores, as well as an estimate score representing tumor purity. We calculated immune, stromal, and Estimate scores for each sample in pan-gynecological cancer using the Estimate algorithm, and analyzed correlations of expression levels of the seven m6A regulators with each score in gynecological pan-cancer. The results are presented in Figure 5 and show that expression levels of the seven m6A regulators were significantly positively or negatively correlated with stromal and immune scores in pan-gynecological cancers.
Correlations of m6A regulator expression and immune infiltration in pan-gynecological cancer
Tumor-infiltrating immune cells are closely associated with tumor metastasis and cancer progression, and may serve as therapeutic targets for improving patient survival [28, 29]. Here, we determined the percentages of immune cells in each pan-gynecological cancer (including CESC, OV, UCEC, etc.) and explored the correlations between levels of the seven m6A regulators and immune infiltration level using the CIBERSORT algorithm. We found that ZC3H13 expression was clearly related to B cells memory, Macrophages M2, T cells CD4 memory resting in CESC, B cells memory, B cells naïve, T cells CD4 memory resting in OV, and NK cells activated, Plasma cells, T cells CD4 memory resting in UCEC. Further FTO expression was with T cells CD4 memory resting in CESC, B cells memory in OV, NK cells activated, B cells memory, T cells CD4 memory resting in UCEC. The relationships between the other five m6A regulators and immune infiltration cells are illustrated in Figure 6. These results reveal that, compared with the other six m6A regulators, ZC3H13 expression has a strong correlation with the level of immune infiltration in pan-gynecological cancers.
Enrichment analysis of m6A regulator expression
We conducted single-gene GSEA analysis for each m6a regulator in pan-gynecological cancers, to explore which regulatory pathways and biological processes are enriched for these molecules. In KEGG pathway GSEA, among immune-related pathways in OV, YTHDF2 and WTAP associated genes were enriched for antigen processing and presentation, the TOLL-like receptor signaling pathway, and natural killer cell mediated cytotoxicity. Genes negatively correlated with ALKBH5 in UCEC were enriched in pathways including autoimmune thyroid disease and antigen processing and presentation, among others. YTHDF2 is a positive regulator of the autophagy pathway in OV, while in UCEC, ALKBH5 was a negative regulator of this pathway (figure 7). The results of GO analysis of m6a in pan-gynecological cancer are presented in Figure 8. FTO was enriched in olfactory receptor activity and sensory perception of smell. In OV, WTAP expression was positively correlated with immune response pathways, such as immune response regulating cell surface receptor signaling and regulation of lymphocyte activation. In CESC, YTHDF2 expression was positively correlated with response to interleukin 12.
We next explored potential therapeutic methods using the Cellminer database, with m6A regulators as targets. Correlations between NCI-60 drug z score and expression of the corresponding NCI-60 cell line RNA-seq/composite were analyzed, where higher z score for a cell line indicated that it was more sensitive to the corresponding drug. Here, we selected drugs approved by the US Food and Drug Agency and drugs undergoing clinical trials. Drugs sensitive to the seven m6A regulators are listed in Figure 9. ZC3H13 expression was positively correlated with sensitivity to the drugs selumetinib, trametinib, hypothemycin, vemurafenib, dabrafenib, and cobimetinib. Further, YTHDF2 expression was negatively correlated with drug sensitivity to dasatinib.