Gene signature of m6A-related targets to predict prognosis and immunotherapy response in ovarian cancer

The aim of the study was to construct a risk score model based on m6A-related targets to predict overall survival and immunotherapy response in ovarian cancer. The gene expression profiles of 24 m6A regulators were extracted. Survival analysis screened 9 prognostic m6A regulators. Next, consensus clustering analysis was applied to identify clusters of ovarian cancer patients. Furthermore, 47 phenotype-related differentially expressed genes, strongly correlated with 9 prognostic m6A regulators, were screened and subjected to univariate and the least absolute shrinkage and selection operator (LASSO) Cox regression. Ultimately, a nomogram was constructed which presented a strong ability to predict overall survival in ovarian cancer. CBLL1, FTO, HNRNPC, METTL3, METTL14, WTAP, ZC3H13, RBM15B and YTHDC2 were associated with worse overall survival (OS) in ovarian cancer. Three m6A clusters were identified, which were highly consistent with the three immune phenotypes. What is more, a risk model based on seven m6A-related targets was constructed with distinct prognosis. In addition, the low-risk group is the best candidate population for immunotherapy. We comprehensively analyzed the m6A modification landscape of ovarian cancer and detected seven m6A-related targets as an independent prognostic biomarker for predicting survival. Furthermore, we divided patients into high- and low-risk groups with distinct prognosis and select the optimum population which may benefit from immunotherapy and constructed a nomogram to precisely predict ovarian cancer patients’ survival time and visualize the prediction results.


Introduction
is the most common modification of eukaryotic mRNA and lncRNA, playing a key role in epigenetics (Xu et al. 2022). The methylation modification of m6A is reversible, and the regulatory factors include methyltransferases (METTL3/14/16, WTAP, RBM15/15B, Generally, different m6A regulators might have distinct functions. However, even the identical m6A regulators may have distinct functions under distinct conditions, mainly due to distinct downstream-targeted genes. Liu et al. reported that METTL3 could interact with CUL4B, and then active PI3K/AKT signaling pathway, which plays a pivotal role in promoting gastric cancer progression (Liu et al. 2021a, b). Xia et al. discovered that METTL3 is involved in apoptosis mediated by Keap1/Nrf2 signaling pathway (Xia et al. 2021). Therefore, identification of the targets of m6A regulators is particularly important to elucidate the dynamic variation.
Sequencing technology has developed tremendously over the past decades, followed by genetic and epigenetic therapies began to emerge, such as targeted therapy and immunotherapy. Disappointingly, for the unscreened population, most patients benefit little or nothing from immunotherapy, far from meeting clinical needs b, c). In general, immunotherapy is closely related to the tumor microenvironment (TME), such as T cell abundance and tumor mutation burden (TMB) (Zhang et al. 2020a, b, c;Galon and Bruni 2019). M6A also played a crucial role in TME (Xia et al. 2021). In this paper, our bold conjecture is that m6A can influence the efficacy of immunotherapy by acting on downstream targets and TME.
In the current work, we screened 9 prognostic m6A regulators, and identified three modification patterns. Next, we identified 7 m6A-related mRNAs and disclosed that they may be potential target of m6A regulators using the m6A2Target database. We also constructed a risk score model to predict overall survival and immunotherapy response. Overall, the low-risk group was more suitable for immunotherapy. Finally, seven identified m6Arelated targets were validated with quantitative PCR and the Human Protein Atlas (HPA) database. The flow of the study is shown in Fig. 1.

Data collection and preprocessing
The expression profiles and clinical information of ovarian cancer were obtained from The Cancer Genome Atlas (TCGA, https:// portal. gdc. cancer. gov/), GSE140082 (GPL14951), and GSE13876 (GPL7759). The GSE13876, including 157 ovarian tumor samples, was used for validation. On the other hand, the expression matrix of normal ovarian tissue was acquired from the Genotype-Tissue Expression (GTEx, http:// commo nfund. nih. gov/ GTEx/ data). Next, we used "ComBat" function of sva package to correct non-biological batch effects. The somatic mutation data and Copy Number Variation (CNV) were obtained from TCGA dataset. Then, the correlation between m6A expression and CNV patterns was analyzed using the Kruskal-Wallis tests. In addition, the 'oncoplot' function of R package 'maftools' was utilized to depict the mutation landscape. The RCircos package in R software was Fig. 1 The flow chart of the study conducted chosen to visualize the m6A regulators in chromosomes (Mayakonda et al. 2018).

Unsupervised clustering for nine m6A regulators and differential expression analysis
We selected 24 m6A methylation regulators from the previous publications for this study (Fan et al. 2020;Wang et al. 2021). Kaplan-Meier (KM) survival analysis was chosen to compare the OS of ovarian cancer patients in distinct m6A groups. In addition, nine m6A regulators were screened. Unsupervised clustering analysis was performed to identify different m6A clusters based on nine prognostic m6A regulators. The R package "ConsensuClusterPlus" was adopted to perform the steps (Wilkerson and Hayes 2010).

Functional annotation between different m6A clusters
Overlapping differentially expressed genes among different m6A clusters were screened using R "limma" package. p < 0.001 and |logFC|> 0.5 were set as the inclusion criteria. The Search Tool for the Retrieval of Interacting Genes (STRING, https:// cn. string-db. org/) protein database and cytoscape was utilized to construct the Protein-protein interactions (PPI) and co-expression networks (Szklarczyk et al. 2015). GSVA algorithm was performed using the R "GSVA" package (Hanzelmann et al. 2013). Finally, based on overlapped differentially expressed genes (DEGs), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted. The significance criteria were set as an adjusted p < 0.05 and |logFC |> 0.1.

Estimation of immune cell infiltration
The single sample gene set enrichment analysis (ssGSEA) and CIBERSORT algorithm were employed to quantify the level of immune cell infiltration (Subramanian et al. 2005;Barbie et al. 2009;Newman et al. 2015). The gene set that labeled each immune cell type was derived from Charoentong's study (Charoentong et al. 2017). Immunophenotypic scores (IPS) of ovarian cancer patients were obtained from The Cancer Immunome (TCIA, https:// tcia. at/ home) database, which can reflect the efficacy of immune checkpoint blockade (ICB) immunotherapy (Charoentong et al. 2017).

Construction of the risk score model
Univariate Cox regression analysis was applied to screen prognostic genes using R "Survival" package. p < 0.05 was set as the inclusion criteria. 98 prognostic genes were screened. Next, 'pacman' package in R software was applied to conduct 'pearson' correlation analysis. The inclusion criteria were set as p < 0.05 and |cor|> 0.3. Finally, the least absolute shrinkage and selection operator (LASSO) regression was employed to further screen for potential prognostic risk characteristics: The R "Survminer" package was utilized to determine cutoff points, and then divided the patients into the highand low-risk groups according to cutoff points. The transcriptional profiles were compared by principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE).

Establishment and validation of a prognostic nomogram
The prognostic value of the risk score was accessed by the univariate and multivariable Cox regression analyses. The R packages "forestplot" were employed to visualize the results. The clinical information, including age, stage and grade were adopted for model correction. Next, a nomogram was established by multivariate Cox regression. Then, the calibration plots and the C-index were calculated for the nomogram models. These analyses were executed by the R package "rms." Receiver-operating characteristic (ROC) (R package "survivalROC") curves were constructed to evaluate the prognostic ability of the nomogram for 1/3/5-year OS and the area under the curve (AUC) values.

Validation of the expression of the genes in risk signature
All the specimens were from the Gynecology in Renmin Hospital of Wuhan University. All the patients provided informed consent and were approved by the Ethics Committee of Renmin Hospital of Wuhan University. We collected contralateral normal ovarian tissue from 3 cases with ovarian cysts and ovarian cancer tissue of 5 cases. Total RNA was isolated with TRIzon reagent (Invitrogen). Primers used in the qRT-PCR are present in Supply Table 1. Reverse transcription was performed using mRNA Reverse-Transcription Kit (Yeasen Biotech Co., Ltd.) according to the manufacturer's instruction. Random primers were added to initiate cDNA synthesis. qRT-PCR was executed using the SYBR Green PCR Mix (Yeasen Biotech Co., Ltd.) on a Bio-Rad CFX96 system. The 2 −ΔΔCt method was performed to calculate the gene relative expression. In addition, to further determine the expression level of the m6A-related genes, immunohistochemical images were obtained from the HPA database (https:// www. prote inatl as. org/).

Potential targets of m6A regulators
The validated and potential targets of m6A regulators were obtained from the m6A2Target databases (Deng et al. 2021).

Statistics
When the samples are divided into two groups, if the samples meet the parameter conditions (normal distribution and homogeneity of variance), t test is employed, otherwise the non-parametric two-sided Wilcoxon-rank sum test was chosen. Furthermore, when the samples were divided into multiple groups, ANOVA test was performed if the samples met parameter conditions, otherwise Kruskal-Wallis tests was employed. The R package "limma" was employed to identify differential transcription factors in distinct group (Ritchie et al. 2015). The threshold for differential transcription factor was as follows: p < 0.05 ("*"), p < 0.01 ("**"), and p < 0.001 ("***"). p < 0.05 was regarded as significant. All analyses were performed with R software V4.1.0 and GraphPad Prism.
The somatic mutations of 24 m6A regulators are 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. Figure 3A shows 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 ( 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, 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).

Consensus clustering of m6A regulators in ovarian cancer
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 (Batlle and Massague 2019; Meurette and Mehlen 2018;Doheny et al. 2020;Katoh and Katoh 2020; (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 has 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 activated CD8 + T cell, activated dendritic cell, MDSC, type 17 and 1 T helper cell, except for METTL14 and WTAP (Fig. 4A). What is 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 suggested that the tumor microenvironment of m6A cluster B may be immunosuppressive (Fig. 4B). Regarding immunosuppressive cells, such as Treg, and TAM, there was no significance among three clusters. However, contrary to our expectations, the level of MDSC was downregulated in m6A cluster B (Fig. 4B). Thus, further exploration of the tumor microenvironment is needed to explain this phenomenon in the future. 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 p < 0.05 ("*"), p < 0.01 ("**"), p < 0.001 ("***"). B Immune infiltration analysis in three clusters. Cluster A, blue; cluster B, orange; cluster C, red tumor-associated signaling pathway, KEGG analysis suggested that overlapping DEGs were enriched in PI3K-Akt signing pathway, focal adhesion, and ECM-receptor interaction.
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 lowrisk groups. Encouragingly, the results from m6A2Target database demonstrated that the seven genes are potential targets for nine m6A regulators ) (Supply Table 4). 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 ( 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). In addition, 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 highrisk 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. 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 Fig. 6 The immune landscape in risk groups. A Differences of stromal score, immune score and estimate score between risk high-score and low-score. B The abundance of TME infiltrating cell. The expres-sion of 13 immune-related functions (C), 9 prognostic m6A (D), cell proliferation-related gene (E), DNA repair-related gene (F), HLArelated gene (G), and EMT-related gene (H) between two risk groups 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-D). Next, we compared the expression of nine 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. 6F-H).
The discovery of immune checkpoints and the development of immune checkpoint inhibitors have opened a new chapter in the war on cancer. The binding of immune checkpoint associated receptors with ligands negatively regulates T cell activation, allowing tumor cells to escape T cell pursuit and inactivation (Pardoll 2012). On the other hand, with the advent of immune checkpoint inhibitor drugs, a very small number of patients can have significantly improved outcomes (de Miguel and Calvo 2020). We observed that a lot of Immune checkpoints, including CD80, CD276, CD86, LGALS9, PDCD1LG2, CD44, LAIR1, NRP1, HAVCR2 and CD200 were significantly upregulated and the expression of several immune co-stimulatory receptor or ligand, such as CD28, TNFSF4 and TNFSF9 were increased in high-risk group (Fig. 7A). Notably, the LGALS9, IDO1, and VTCN1 were elevated in low-risk group. To evaluate which groups are more suitable candidates for immunotherapy, we calculated the IPS and TMB between two group. Surprisingly, IPS was higher in low-risk group, implicating a better response to ICBs therapy (Fig. 7B). TMB was also regard as a specific marker to guide the ICBs. Higher TMB was higher neo-antigenicity, which subsequently improved the response of ICBs in tumors. We observed 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, which are warrant to be verified in more by more experiment (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).
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 Fig. 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). In addition, the representative immunohistochemical images from the HPA database were obtained. In addition, the results were in line with qRT-PCR ( Figure S7). Regarding the relationship between 7 genes and tumor-infiltrating lymphocytes, we obtained the expression profile of 157 patients with advanced stage serous ovarian cancers from GSE13876. The result of correlation analysis suggested that most of the seven genes were positively correlated with Tregs, and negatively related with activate memory CD4 + T cells. In addition, the seven genes tended to be negatively correlated with the tendency of CD8 + T cells, but were not statistically significant ( Figure S8).

Discussion
Ovarian cancer is a heterogeneous disease. Even among patients with similar characteristics and clinical stage, the prognosis may vary greatly. Accumulated evidence showed that m6A modification by targeting downstream genes played a considerable role in tumor proliferation, differentiation, development, invasion, and metastasis (Xu et al. 2021;Hua et al. 2018;Liu et al. 2020Liu et al. , 2018. Accumulated evidence reveled that m6A regulators and their targets may play a crucial role in the benefit of immunotherapy (Cong et al. 2021;Song et al. 2021;Liu et al. 2021a, b). Yet, the studies to date have tended to focus on m6A regulators rather than their targets and far too little attention has been paid to the potential role of m6A regulators and their targets influencing ovarian cancer immunotherapy. Hence, this study was undertaken to investigate the potential functions of m6A regulators and their targets and identify the candidate crowd that would benefit from ICB therapy.
In this work, we compared the different expressions of normal tissues with 24 m6A regulators of ovarian cancer patients, and then screened out nine prognostic m6A which revealed three m6A modification patterns with different TME cell infiltration in the landscape. What is more, we researched the mutation of m6A and CNV in ovarian cancer. The results indicated that CNV rather than mutation was closely related to m6A expression. Significantly, compared to the other two groups, m6Acluster B had the worst prognosis, with an infiltrative feature of immunosuppression. In addition, the result of GSVA demonstrated that malignant functional features of the tumor were significantly enriched in m6Acluster B. Considering the heterogeneity of m6A modification, we selected 97 overlapping DEGs which were correlated with prognosis of three m6A clusters to quantify the pattern of m6A modification in individual tumors. The signature genes and signal pathway in current work can provide ideas for laboratory experimental design to elucidate the molecular mechanism of ovarian cancer. Next, results from LASSO analysis revealed seven key genes, namely ELK3, AFF4, RCOR2, PHACTR2, CIZ1, NUAK1, and MCC, which were associated with m6A regulators and immune landscape. ELK3 has been reported that plays a key role in the migration and apoptosis of cancers and is associated with infiltration of M2 macrophages, CD4 + and CD8 + T cell Dazhi et al. 2020;Meng et al. 2020;Kim et al. 2020). As for RCOR2, it is known as an immunomodulatory gene and is related enriched in with low CD8 + T cell infiltration tumors (Routh et al. 2020). In addition, Xiong proved that Rcor1, Rcor2, and Lsd1 impaired Foxp3 + Treg function and promoted antitumor immunity (Xiong et al. D A nomogram to quantitatively predict survival based on risk score and clinical parameters. Calibration curves (E) of the nomogram and receiver-operating characteristic (ROC) curves (F) to estimate the accuracy and performance of the predictive nomogram 2020). However, except for ELK3 and RCOR2, there was no relevant study on five others and tumor immune-infiltrating cells. AFF4 is a validated downstream target of m6A regulators and is downregulated in ovarian cancer Liang et al. 2018). It has been extensively studied in the proliferation, migration, and invasion of cancers Han et al. 2019;Deng et al. 2018). Similar to AFF4, NUAK1, a member of ARK5 family, is negatively correlated with prognosis in ovarian cancer and implicated in tumor cell apoptosis and tumor cell invasion and metastasis (Monteverde et al. 2018;Phippen et al. 2016;Zhang et al. 2015;Riester et al. 2014). There is also abundance of evidence that CIZ1 and MCC were closely related with oncogenesis and tumor progression (Fukuyama et al. 2008;Xiao et al. 2017;Kinzler et al. 1991;Higgins et al. 2012;Pauzaite et al. 2016). There are limited studies on the role of PHACTR2 in tumorigenesis compared with others. PHACTR2 is identified as a candidate hub gene for lung cancer and esophageal squamous cell carcinoma (Musolf et al. 2021;Wang et al. 2013;Chen et al. 2020a, b). In this article, our results suggested that there may be m6A regulators related targets and may connect with different immune landscape, which provide a new insight into the mechanism of those genes in the pathogenesis of ovarian cancer.
To further explore the prognostic value and biological mechanisms of key targets, we constructed a m6A-related prognostic model to provide novel insights for the treatment of ovarian cancer. According to the best cutoff value of the risk score, the patients were stratified into high-and low-risk subgroups with different clinical outcomes. Moreover, the results of PCA and t-SNE showed that the risk model is discriminative. Next, we investigated the relationship between the risk score and clinical features and discovered that it was strongly associated with prognosis.
The tumor microenvironment continually changes over the course of cancer progression and the composition of the tumor microenvironment influences the response of the ICB. Usually, the responses to ICB therapy depend on the TMB and tumor-infiltrating lymphocytes (Li et al. 2020a, b). To explore whether there is a possible connection between ICB therapy and risk score, we researched the immune-infiltrating and TMB in high-and low-risk groups. There was no statistically significant difference in immune scores. In addition, patients in high-risk groups 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, activated CD4 + and CD8 + T cell, and type 17 helper cell had been discovered to be mostly enriched in low-risk group. Moreover, cytolytic activity, HLA, inflammation promoting, T cell co-stimulation and type I IFN response were significantly elevated in low-risk group. Chang, et al. proposed that the orientation and quality of the immune response influenced the likelihood of immunotherapy response. The results showed that the orientation of immune response was ponded to immune stimulatory and the quality tended to be "hot" T cell infiltration in low-risk group compared with high-risk group, which prompted that the likelihood of immunotherapy was higher. In addition, the cell proliferation, DNA repair, and EMT-related signaling pathways were enhanced in high-risk group, which may be associated with poor prognosis. Regarding immune checkpoint, the expression levels of IDO1, VTCN1 and LGALS9 in low-risk group were higher. Higher TMB has higher neoantigenicity, which subsequently improved the response of ICBs in tumors. We combined TMB and risk score to divide the population into four groups for prognostic survival analysis. The results elucidated that there was a best prognosis in high-TMB and low-risk subgroup. Moreover, IPS can effectively predict the effect of tumor immunotherapy and there was higher IPS in the low-risk group. Collectively, the above results suggested that the low-risk group might be the best candidate to receive immunotherapy.
Importantly, we demonstrated that the risk score has a good clinical prognostic stratification effect in advanced ovarian cancer. Univariate and multivariate analyses demonstrated that risk score was an independent prognostic factor for ovarian cancer. What's more, integrating other clinical features, we constructed a nomogram to precisely predict ovarian cancer patients' survival time and visualize the prediction results. Emerging sequencing technologies make it possible to individualize risk stratification and treatment at the molecular level. In this work, we innovatively combined risk score to construct a clinical powerful model for prediction of prognosis and ICBs.
Undeniably, there are several limitations in this study. First, we only use m6A2Target database to verify the m6A targets. Further experimental studies are needed to confirm the regulation of m6A modification on these targets and elucidate its role in ovarian cancer progression and immune escape. Secondly, all data came from public databases, undermining the future applicability of the results to the clinic, and we only determined the cutoff value based on the data of TCGA and GSE140082 to distinguish the two groups of high and low scores. Hence, large sample multicenter prospective is warranted to verify our risk model grouping and judgment of immune clustering effect and enhance clinical applicability, which may be our followup experimental plan. Finally, about the threshold of risk grouping, there may be some bias, and the best cutoff value still needs more research to explore and confirm. Even so, the rapid development of biological technologies will hopefully exploit the way for its realization.
Despite these limitations, our study disclosed that the CNVs may contribute to the dysregulation of m6A regulators, and identified three m6A modification patterns 1 3 with distinct immune landscape. We then identified the overlapped DEGs in three m6A modification patterns and screened m6A-related targets using correlation analysis and m6A2Target database. Overall, we detected seven m6Arelated targets as an independent prognostic biomarker for predicting survival. We comprehensively analyzed immune-infiltrating, immune-related pathways, tumor biology-related signaling pathways, immune checkpoints, IPS and TMB between high-and low-risk groups and revealed that the low-risk group is the best candidate population for immunotherapy.