Genetic alteration of m6A regulators in ovarian cancer
In TCGA cohort, nine out of 21 m6A regulators (RBM15, RBM15B, LRPPRC, ZC3H13, IGF2BP2, CBLL1, YTHDF1/2/3, and ELAVL1) had significantly high expression level in ovarian cancer. While the expression of METTL3/14/16, WTAP, FTO, HNRNPA2B1, HNRNPC, YTHDC1/2 and RBMXwere decreased in ovarian cancer (Fig. S1A, S1C). Furthermore, K-M analysis was carried out to evaluate the clinical significance of the 24 m6A regulators in ovarian cancer, nine of them were correlated with the overall survival (OS) (Fig. 2). The patients with over-expressions of METTL3, HNRNPC and YTHDC2 had longer OS. The PPI analysis revealed that WTAP, RBM15, METTL3, METTL14, FTO and HNRNPA2B1 had more internal connections with others, while IGF2BP2, IGF2BP1, IGF2BP3, LRPPRC and EIF3A had less interaction with others (Fig. S2A). In addition, except for IGF2BP1, there were significant positive correlations between m6A regulatory factors (Fig. S2B).
The somatic mutations of 24 m6A regulators were mapped in Figure S1B, which showed that the mutations in m6A regulators were rarely observed in ovarian cancer. Only 11.7% of the 436 patients had somatic mutations, and m6A “reader” YTHDF1had the highest frequency of mutations in ovarian cancer (Fig. S1B, Supplementary file 8). By the contrast, a larger number of broad CNVs were existed in m6A regulators, and most of which leads to amplification, but deletion of YTHDF2, YTHDC2, RBM15B and WTAP (Fig. S1D, Supplementary file 8). Furthermore, CNV has been proved that positively correlated with the expression of m6A regulators (Fig. S3). In summary, CNVs may contribute to the abnormal expression of m6A regulators, and further effected the progression and clinical outcomes of ovarian cancer.
Consensus clustering of m6A regulators in ovarian cancer
Figure 3A showed the crosslink among the 24 m6A regulators and their clinical significance (Fig. 3A). We observed the significant crosslink not only existed in the m6A regulators with same function, but also between the “writers”, “erasers” and “readers”. The above results suggested that the interconnection of the “writers”, “erasers” and “readers” might form a variety of modification patterns in ovarian cancer. The cohort was divided into three subgroups (m6Acluster A: N = 266; m6Acluster B: N = 183; m6Acluster C: N = 307) by applying the unsupervised clustering of nine prognostic m6A regulators. The optimal number of stable gene clusters were three (k = 3) by the consistent clustering algorithm (Fig. S2C). A clear separation was observed among them (Fig. S2D). Next, we compared OS between each cluster and proved that m6Acluster B had decreased OS than m6Acluster A and m6Acluster C (p = 0.015, Fig. 3B). Nevertheless, there was no statistical difference in OS between m6Acluster A and m6Acluster C. Notably, ovarian cancer patients in m6Acluster B had more expression of m6A (Fig. 3C).
As shown in the Fig. S4, metabolic pathways like inositol phosphate, glutathione metabolism and oxidative phosphorylation were enriched in m6Acluster A. While m6Acluster B were significantly related to carcinogenesis, including TGF-β signaling, Notch signaling, Hedgehog signaling, WNT signaling and ubiquitin mediated proteolysis (Fig. S4A-B) (Batlle and Massague 2019, Meurette and Mehlen 2018, Doheny et al. 2020, Katoh and Katoh 2020, Jiang et al. 2020). The above results suggested that m6A regulators were involved in many crucial cellular activities, such as carcinogenesis, and inflammation of ovarian cancer.
Immune landscape was significantly associated with the m6A modification pattern
Infiltration of immune cells within the TME have a pivotal role on the tumor progression and the efficacy of immunotherapy. The infiltration of B cells, CD4+ T cells, CD8+ T cells, dendritic cells, NK cells and T helper cell was low in m6A cluster B which presented immune-excluded patterns (Fig. 4A). Additionally, we explored associations between 9 prognostic factors and 23 types of immune cells and found that most of the m6A were negatively correlated with immune cells, except for METTL14 and WTAP (Fig. 4B).
Construction and verification of the m6A-score model
357 overlapping differentially expressed genes among the m6Acluster A, m6Acluster B and m6Acluster C were identified in this study (p < 0.001, |logFC| > 0.5, Fig. S5A). The result of GO analyses showed that overlapping differentially expressional genes (DEGs) were related to cell surface receptor, protein kinase activity, extracellular matrix, and cytoplasmic region (Fig. S5C). In the tumor-associated signaling pathway, KEGG analysis suggested that overlapping DEGs were enriched in PI3K-Akt signing pathway, focal adhesion, and ECM-receptor interaction (Supplementary file 8).
Then, 97 overlapping DEGs were correlated with prognosis by univariate cox analysis (p < 0.05) (Supplementary file 8). The correlation between 97 genes and 9 prognostic m6A regulators was analyzed, and shown in Figure S5B. As shown in Fig. 5B, 47 genes were strongly correlated with 9 prognostic m6A regulators (p < 0.05, |cor| > 0.3). Furthermore, co-expression network of the m6A-related genes and m6A regulators was plotted (Fig. S5D).
To avoid over-fitting of the prognostic model, 47 differentially expressed m6A-related genes with prognostic values were analyzed by LASSO regression analysis, and 7 genes were screened for prognostic model (Fig. 5A-B). With the cutoff value 6.018, patients were divided into high or low-risk groups (Supplementary file 8). Encouragingly, the results from m6A2Target database demonstrated that the seven genes are potential targets for 9 m6A regulators (Supplementary file 8) (Cheng et al. 2019). There was a better prognosis in low-risk group (p < 0.001, Fig. 5C). The area under ROC curve was 0.627 (Fig. 5D). Moreover, PCA and t-SNE showed that the two groups were presented in two different patterns (Fig. 5E-F).
We likewise examined the correlation between risk score and clinical characteristics. There was differential expression of risk scores between different grades. The risk score of G1-2 is higher than that of G3-4 (Fig. S6A). Due to the difficulty in early detection of ovarian cancer, most clinical samples were found to be in the middle and advanced stages, so the number of patients with G1-2 was slight, and the experimental results may be biased. Although the box plot showed risk scores at an advanced stage as well as age > 65 had higher levels, but this was not statistically significant (Fig. S6B-C). Additionally, we respectively studied the influence of risk group on prognosis under different clinical characteristics, and the results proved that, in all clinical subgroups, the prognosis was better in low-risk group (Fig. S6D-F). In the prognostic model, the risk scores and the number of deaths were elevated in high-risk group (Fig. S6G).
Evaluation of the correlation between prognostic risk score and ICBs
Immunotherapy represented by ICBs has shown impressive therapeutic efficacy with long-lasting responses in a minor subset of patients. In general, immunotherapy is closely related to TME, like T cell abundance and TMB. Thus, we compared the scores related to TME between two groups and hoped to screen the optimal groups for immunotherapy.
The stromal score was significantly elevated in high-risk group (Fig. 6A). However, there was no significant difference in the immune score. We next investigated the immune infiltration in two groups. The resulted indicated that patients in high-risk group had enhanced adaptive infiltrating immune cells, including eosinophil, macrophages, mast cell, dendritic cell, and regulatory T cell et al. On the other hand, activated B cell, CD4+ and CD8+ T cell, and type 17 helper cell had been discovered to be mostly enriched in low-risk group (Fig. 6B). The above results showed that it was the proportion of immune cells, rather than the content of immune cells, that differed between the two groups. Furthermore, the results of immune-related pathways elucidated that cytolytic activity, HLA, inflammation promoting, T cell co-stimulation and type I interferon (IFN) response were significantly elevated in low-risk group. High-risk group had a higher enrichment of type II IFN response, which plays critical roles in orchestrating both innate and adaptive immune responses (Fig. 6C). Next, we compared the expression of 9 prognostic m6A regulators between two groups. The results demonstrated that the high-risk group had higher expression of m6A regulators (Fig. 6E). Ultimately, the tumor biology-related signaling pathways were analyzed and the results indicated that cell proliferation, DNA repair, and EMT related signaling pathways were enriched in high-risk group (Fig. 6D-G).
The discovery of immune checkpoints and the development of immune checkpoint inhibitors have opened a new chapter in the war on cancer. Immune checkpoints CD80, CD276, CD86, LGALS9, PDCD1LG2, CD44, LAIR1, TNFSF4, NRP1, TNFSF9, CD28, HAVCR2, CD200were significantly upregulated in high-risk score, showing a consistent trend with poor prognosis (Fig. 7A). In addition, IPS was higher in low-risk group, implicating a better response to ICBs therapy (Fig. 7B).
TMB was regard as a specific marker to guide the ICBs. Higher TMB was higher neo-antigenicity, which subsequently improved the response of ICBs in tumors. The R package “maftools” was performed to analyze the effect of TMB on the survival prognosis of ovarian cancer, our results revealed that the higher TMB was correlated with better prognosis (Fig. 7D). TMB was integrated with risk score for prognosis analysis, and the results demonstrated that patients in the group with low risk-score and high TMB had the best survival prognosis compared with the other three groups (Fig. 7E). Subsequently, we investigated the association between TMB and risk-score and compared the differences in the distribution of somatic mutations between high- and low-risk groups and observed a negative correlation between TMB and risk score, and the TMB was higher in the low-risk group, but there was no statistic difference (Fig. 7C). In short, the burden of tumor mutations was more extensive in the low-risk group (95.88 vs. 94.17%) (Fig. 7G-H). The correlation between risk score and TMB is not significant, which may be related to insufficient samples.
Validation of the m6A signature
We referred risk score as independent protective factor (p < 0.001, 95% CI HR: 1.215-1.516) (Fig. 8A), and this association remained statistically significant (p < 0.001, 95%CI HR: 1.195-1.486) after adjusting for age, stage, and grade (Fig. 8B).
To provide a quantitative method to predict the probability of survival, we constructed a nomogram that integrated both the risk score and clinical characteristics (Fig. 8D). The calibration plots depicted in Figure 8E showed that the derived nomogram performed well. Similarly, nomogram model showed predictive ability at 1, 3 and 5-year overall survivals (1-year AUC = 0.679, 3-year AUC = 0.661, 5-year AUC = 0.624) (Fig. 8F). Integration of the risk score with clinical characteristics in the nomogram showed high specificity and sensitivity of survival prediction. To further confirm the expression of the seven prognostic genes, we obtained ovarian cancer and normal ovary tissues from patients including 5 cancer samples and 3 normal tissues. The results demonstrated that almost genes were significantly downregulated in ovarian cancer (Fig. 8C). Additionally, the representative immunohistochemical images from the HPA database were obtained. And the results were in line with qRT-PCR (Figure S7).