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 RBMX were decreased in ovarian cancer (Figure 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) (Figure 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 (Figure S2A). In addition, except for IGF2BP1, there were significant positive correlations between m6A regulatory factors (Figure 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” YTHDF1 had the highest frequency of mutations in ovarian cancer (Figure S1B, Supplementary Table 2). 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 (Figure S1D, Supplementary Table 3). Furthermore, CNV has been proved that positively correlated with the expression of m6A regulators (Figure 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 (Figure 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 (Figure S2C). A clear separation was observed among them (Figure 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, Figure 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 (Figure 3C).
As shown in the Figure 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[28-32] (Figure S4A-B). We also found that m6Acluster C was associated with inflammation pathways like KRAS signaling, interferon alpha and interferon gamma response, reactive oxygen species pathway (Figure S4C). 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. We explored associations between 9 prognostic factors and 23 types of immune cells. The results suggested that most of the m6A were negatively correlated with immune cells, except for METTL14 and WTAP (Figure 4A). What’s more, we investigated the immune infiltration in different m6A modified pattern and found that 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 (Figure 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, Figure 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 (Figure 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.
Then, 97 overlapping DEGs were correlated with prognosis by univariate cox analysis (p < 0.05). The correlation between 97 genes and 9 prognostic m6A regulators was analyzed, and shown in Figure S5B. As shown in figure 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 (Figure 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 (Figure 5A-B). With the cutoff value 6.018, patients were divided into high or low-risk groups. Encouragingly, the results from m6A2Target database demonstrated that the seven genes are potential targets for 9 m6A regulators (Supply Table 4). There was a better prognosis in low-risk group (p < 0.001, Figure 5C). The area under ROC curve was 0.627 (Figure 5D). Moreover, PCA and t-SNE showed that the two groups were presented in two different patterns (Figure 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 (Figure 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 (Figure 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 (Figure S6D-F). In the prognostic model, the risk scores and the number of deaths were elevated in high-risk group (Figure 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 (Figure 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 (Figure 6B). The above results showed thatit 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 (Figure 6C-D). 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 (Figure 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 (Figure 6F-H).
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, CD200 were significantly upregulated in high-risk score, showing a consistent trend with poor prognosis (Figure 7A). In addition, IPS was higher in low-risk group, implicating a better response to ICBs therapy (Figure 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 (Figure 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 (Figure 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 (Figure 7C). In short, the burden of tumor mutations was more extensive in the low-risk group (95.88 vs. 94.17%) (Figure 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) (Figure 8A), and this association remained statistically significant (p < 0.001, 95%CI HR: 1.195-1.486) after adjusting for age, stage, and grade (Figure 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 (Figure 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) (Figure 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 (Figure 8C). Additionally, the representative immunohistochemical images from the HPA database were obtained. And the results were in line with qRT-PCR (Figure S7).