Public RNA-seq expressed data and full clinical annotations are available from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. The exclusion criteria include removing all samples without clinical follow-up information, removing all samples with unknown survival time < 30 days and removing all samples without survival status. After excluding samples with incomplete survival data, a total of 3 eligible OC cohorts, including GSE9891, GSE29691 and TCGA-OV, were collected to be processed. RNA sequencing data of gene expression (FPKM values) and clinical information in TCGA dataset were accessed directly from GDC website (https://portal.gdc.cancer.gov/). FPKM data were converted to transcripts per kilobase million (TPM) data. Normalized RNA sequencing data and corresponded clinical information for GSE9891, GSE29691 microarray dataset were accessed directly from GEO website (http://www.ncbi.nlm.nih.gov/geo). Batch effects between these cohorts were removed using “sva” R package (34). OC patients with mutational data were obtained from TCGA database. A total of 16 OC specimens and 16 normal tissues were obtained from the first affiliated hospital of Nanjing medical university. We obtain all the written informed consent from patients.

Total 10 regulators were selected from TCGA and GEO datasets for determining m1A modification patterns. On the basis of expression of 10 m1A regulators, we used the ConsensuClusterPlus package to apply unsupervised clustering analysis in order to facilitate the identification of different m1A modification patterns and the classification of patients. The consensus clustering algorithm determines how many clusters there are, and the process is repeated a thousand times to ensure reliability (35, 36).

GSVA, a non-parametric unsupervised method, can be applied to assess the difference gene set enrichment between different m1A modification modes (37). Download the gene set “c2.cp.kegg.v6.2.symbols” from MSigDB database and use it for running GSVA. An adjusted P<0.05 was regarded as statistically significant.

CIBERSORT provides expression data for 22 common immune cells LM22. Based on these data, we calculated the association between m1Ascore and immune cell infiltration (38, 39). We utilized the ssGSEA (single sample gene set enrichment analysis) algorithm to assess the degree on immune infiltration. ESTIMATE algorithm was utilized to score stromal and immune gene sets and calculate tumor purity (40).

We applied the empirical Bayesian approach of the limma package to find out DEGs of the three m1A modification patterns (41). The adjusted p for these genes is less than 0.05.

In order to quantify m1A modification pattern of single OC patients, we built a score scheme to assess m1A modification patterns, which we termed as m1Ascore. The m1Ascore was built in the following steps: first, we normalized the DEGs extracted from the different m1Aclusters and extracted overlapping DEGs. overlapping DEGs were analyzed using unsupervised clustering to classify OC patients into several groups. Consensus clustering algorithms were used to determine the count of gene clusters and their stabilization. We conducted prognostic analysis for each overlapping DEG with Cox regression method and screened genes with P<0.05. Principal component analysis (PCA) was used for building m1A-associated gene signature. Principal components 1 and 2 were both chosen for feature scores. m1Ascore = ∑ (PC1i + PC2i). Where i is m1A phenotype-related genes’ expression (42, 43).

IPS is a representative gene associated with immunogenicity calculated using z-score. It uses PD-L1 expression of the four tumor-associated immune cells individually as an evaluation metric to distinguish the beneficiary population. Higher scores were linked to high immunogenicity (44). OC patients’ IPS comes from the Cancer Immunome Atlas (TCIA) (https://tcia.at/home).

CMap, as a gene expression profile database, allows the comparison of differentially expressed gene profiles with database reference datasets to obtain highly correlated agents with diseases (45). We utilized the CMap database to obtain the linkage between m1Ascore, OC and drugs. The 3D structures of the obtained agents can be accessed from Pubchem website.

The Human Protein Atlas (HPA) (https://www.proteinatlas.org) includes tissue and cellular protein distribution data from 44 different normal tissue categories and 17 major cancer categories. Immunohistochemical staining intensity, number, location and patient information are available online. After exploring 10 m1A-related genes in HPA database, seven m1A regulators expression (ALKBH1, ALKBH3, TRMT6, TRMT10C, TRMT61B, YTHDC1 and YTHDF2) in normal and OC tissues was obtained.

Trizol reagent (Thermo Fisher Scientific, USA) was utilized to isolate and extract RNA from tissue samples. NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) was designed to assess RNA quantity control and concentration. A high-capacity reverse transcription kit (Takara, Japan) was designed for reverse transcription of total RNA to cDNA. qRT-PCR was conducted in a Light Cycler 480II (Roche) using SYBR Green technology (Takara). Record the cycling threshold (Ct) for each gene and calculate the target gene mRNA expression with the 2-ΔΔCt method. All steps of the qRT PCR were performed according to the reagent instructions and all experiments were repeated 3 times. PCR primers are showed in Table S1.

Correlation coefficients between immune cells and m1A regulators expression were calculated by Spearman correlation analysis. Differences between three groups were compared with Kruskal-Wallis test, and associations among categorical covariates were tested with χ2 test (46). Based on the relevance of m1Ascore to patient prognosis, the optimal cut-off value for each dataset sub-group was defined using the survminer R package. This value dichotomizes the patients into high and low m1Ascore subgroups. Log-rank statistics were used for reducing the batch effect of calculations. Kaplan-Meier method was applied for drawing OS plots and log-rank test was utilized to identify statistics differences. Univariate Cox regression were used to compute risk ratios for m1A regulators and genes associated with m1A phenotype. Multivariate Cox regression was applied for identifying independent survival factors as well as the “forestplot” R package to visualize the results. The Maftools package and its “oncoplot” feature are used to present mutational differences. P < 0.05 is considered statistically significant. All data were processed in R 3.6.1 software.

The total workflow is as shown in the following figure (Supplementary Figure 1). According to previous literature reports, we identified 10 m1A methylated genes. First, we integrated somatic mutations and copy number variants (CNVs) of these genes to characterize the mutations. As can be seen in Supplementary Figure 2, the mean mutation rate of m1A regulators is very low, and only 7 of 436 specimens have m1A regulator mutations, with a frequency of 1.61%. Subsequently, we analyzed the change frequency of the CNV in m1A regulators. Figure 1A shows that copy number variation was present in all 10 regulators, of which, most were dominated by copy number amplification, with only YTHDF2 exhibiting a copy number deletion. Figure 1B further visualizes the location of CNV alterations in the m1A regulators on the chromosome. By analyzing the expression of these 10 m1A regulators, we identified two subgroups that did not cross over, which demonstrated that we could completely distinguish between OC and normal samples based on m1A regulators expression levels (Figure 1C).

For searching the link between m1A methylation regulators and OC, we compared m1A regulator mRNA expression in OC and normal tissues from two databases, TCGA and GEO, respectively. TCGA data showed aberrant expression of all m1A regulators in OC. The expression of all genes in OC was greater than that in normal (p<0.001, Figure 1D). However, the different expression of these m1A regulators was not the same as that in GSE27651 database. In GSE27651 dataset, there was no statistical difference in ALKBH3 and TRMT61A expression between OC and normal samples. And compared with normal tissues, ALKBH1 expression (P<0.01) was significantly decreased in OC, while the expression of YTHDC1 increased (P<0.05, Figure 1E). Furthermore, we compared the differences in m1A regulator expression at the protein level between two tissues. The results proved that ALKBH3 expression was lower in tumor compared to normal tissues, TRMT61A expression differences were not meaningful and the rest were expressed higher in OC (Supplementary Figure 3A). The above analysis indicates a high degree of heterogeneity in genetic variation and expression differences between normal and OC tissues, which indicates that the imbalance in the expression of m1A regulators makes a critical function in the development of OC and that genetic variation may interfere with their expression levels.

Next, we tried to clarify the relationship among the 10 m1A RNA methylation regulators. Analysis of the String database (Figure 1F) showed that ALKBH1 may be a central gene of the m1A RNA methylation regulator. However, further analysis did not display an expression correlation between ALKBH1 and other regulators. Interestingly, in TCGA cohort, we found that TRMT10C was positively related with seven m1A regulators, particularly YTHDF1 and YTHDF2. In contrast, YTHDC1 was negatively correlated with other eight regulators, particularly YTHDF1 and TRMT10C (Figure 1G). This suggests that TRMT10C and YTHDC1 may be key genes in m1A RNA methylation regulators affecting OC occurrence and development. However, in the GSE27691 cohort, we found that TRMT61A was negatively correlated with the expression of seven m1A regulators, particularly YTHDC1 and YTHDF1. In contrast, YTHDF1 was positively correlated with six other regulators, particularly YTHDC1 (Figure 1H). Broadly speaking, the correlation trends between the regulators in the TCGA and GEO databases were generally consistent.

We discuss the prognostic relevance of m1A regulators in patients with OC (Supplementary Figure 3B). We found that some m1A regulators, including ALKBH1, ALKBH3, TRMT6, TRMT10C, YTHDF1 and YTHDF2 are oncogenic and that overexpression of these genes results in worse prognosis for OC patients. The m1A regulator network diagram further describes the connections and interactions between m1A regulators and their prognostic significance for OV patients (Figure 1I). We revealed that m1A regulators expression showed significant correlations both within the same category and between different categories.

According to 10 m1A regulator expression, we used R program to classify OC patients with different m1A modification patterns and eventually identified three different modification patterns by using an unsupervised clustering approach, namely m1A cluster-A (227 patients), B (148 patients), and C (236 patients) (Figure 2A and Supplementary Figures 4A–C). To discover the prognostic worth of the three m1A modification patterns, we concluded Kaplan-Meier analysis, the survival curves showed that the m1A cluster-B modification pattern was more pronounced among the three m1A modification subtypes (P = 0.009, Figure 2B). We also plotted heat maps to discuss the correlation between clinicopathological characteristics of OC patients and m1A modification patterns (Figure 2C).

Subsequently, we performed GSVA enrichment analysis of the three m1A modification patterns to investigate the associated pathways and biological significance. As shown in Figures 2D, E, m1A cluster-A and m1A cluster-C are similar in that both are significantly enriched for a number of gene replication and repair-related pathways, including homologous recombination, mismatch repair, and DNA replication. Differently, ECM-receptor interaction and focal adhesion were heavily abundant in m1Acluster-B, which may be related to the cellular matrix and its role in TME.

To better understand the relationship between m1A modification patterns and immunity, we explored immune characteristics of different m1A modification patterns. Analysis of cellular infiltration in the TME showed that the three m1A modification patterns had distinctly different TME cellular infiltration characteristics, with the most abundant immune cellular infiltration in m1A cluster-B, including many CD4+ T cells, plasmacytoid dendritic cells, monocyte, NK cells, MDSC, etc. (Figure 3A). This immunological profile corresponds to a tumor immune-inflamed phenotype (47), and tumor patients with this profile are the most responsive to immunotherapy, and m1A cluster-B patients do show a matched survival advantage (Figure 2B).

Furthermore, the results of GSVA analysis showed (Figure 2E) that the cluster B modification pattern significantly correlated with intrastromal interactions. Therefore, we hypothesized that the matrix in cluster B activates the anti-tumor effects of immune cells. We then employed ESTIMATE algorithm to calculate the proportion of the immune matrix component of each m1A modification pattern in the TME and presented it as immune score, stromal score and estimate score, which were positively correlated with immunity, stroma and the sum of both. We found that clusters-A and B had significantly more immune and stromal components in TME than cluster-C (p<0.05, Figure 3B).

In clinical tumor immunotherapy, the main immune checkpoint inhibitors currently in use include CTLA-4 inhibitors and PD-1 inhibitors. In order to advance the study of ICIs, understanding immune checkpoints expression in three m1A modification patterns is necessary. Therefore, we explored immune checkpoints expression in different m1A modification patterns. As seen in Figure 3C, immune checkpoint expression was slightly different in m1A subtypes, with slightly higher immune checkpoint expression in cluster A compared to cluster C (p < 0.05).

To explore the biology of each m1A modification pattern, we identified 527 m1A-associated DEGs (Figure 4A) and conducted GO and KEGG analysis on these genes (Supplementary Figure 4D). Based on 527 m1A phenotype-related genes, an unsupervised cluster analysis was then performed, and the patients were divided into different genomic groups. Similarly, unsupervised clustering algorithm exhibited three m1A modification genomic phenotypes, which termed as m1A gene clusters A, B, C (Figure 4B and Supplementary Figures 4E–G). This demonstrates that there are indeed three distinct patterns of m1A modification in OC. We observed that the three different gene clusters have different clinical characteristics (Figure 4B). Compared with m1A gene cluster-A and cluster-C, OC patients in m1A gene cluster-B displayed poor differentiation and higher frequency of recurrence (Figure 4B); and of the three m1A gene clusters, patients in gene cluster-C had longer survival times (p=0.001, Figure 4C). Moreover, m1A regulators’ expression in the three m1A gene clusters also had significant differences (Figure 4D). In addition, we performed immune correlation analysis on the m1A gene cluster A-C. As shown in Figure 4E, the three m1A gene modification patterns have different TME cell infiltration characteristics. Unlike the m1A modification pattern, the m1A gene modification pattern has the most abundant immune cell infiltration in m1A gene cluster-A, such as activated CD8+T cell, MDSC, immature DC, monocyte, plasmacytoid DC, etc. After applying ESTIMATE algorithm, we found that gene cluster-A has the most immune and matrix components in TME, followed by gene cluster-B (P<0.05, Figure 4F). These findings demonstrate that m1A methylation modifications play an important regulatory effect on building of different TME landscapes.

Due to individual differences in m1A modification, based on 527 DEGs, we built a score system to accurately assess m1A modification pattern of single OC patients, and termed it the m1Ascore. First, the alluvial map visually displays the attribute changes of a single OC patient (Figure 5A). To better characterize the immunological profile of m1A feature, we analyzed the association between immune cells and m1Ascore (Figure 5B). Kruskal-Wallis test showed significant differences in m1Ascore between m1A gene clusters. Gene cluster-C scored highest median value and the lowest was for gene cluster-B (Figure 5C). In contrast, the highest median value was observed in m1A cluster B, while the lowest median value was observed in m1A cluster C (Figure 5D). Combining previous survival analysis in different clusters (Figures 2B and 4C), this indicates that the m1Ascore may be positively correlated with the survival of OC patients.

Until we determined the prognostic value of m1Ascore, we evaluated the significance of the m1Ascore in immunological terms. First, we examined the correlation between immune cells and m1Ascore, and found that M1 macrophages, CD8⁺ T cells, gamma delta T cells, T follicular helper cells decreased as m1Ascore increased, with a negative correlation between the two; in contrast DCs and resting CD4 memory T cells were positively related to m1Ascore, and cell infiltration increased with increasing score (Figure 5E). Subsequently, following the optimal cutoff values defined from survminer package, we dichotomized OC specimens into high- and low-m1Ascore groups and initially assessed the TME immune infiltration in these two groups. Figure 5F shows that intrinsic immune cells, including monocyte, NK cells and plasmacytoid DCs, were predominantly increased in high m1Ascore subgroup. On the contrary, there are more adaptive immune cells, consisting of B cells, CD4+ T cells and CD8+ T cells in low m1Ascore subgroup. Similarly, we apply the ESTIMATE algorithm between the two groups, and noticed no significant difference in infiltrating components in TME between the two groups (Supplementary Figure 4H). Furthermore, GSEA revealed that the low m1Ascore group mainly enriched some metabolic pathways and gene replication and repair-related pathways, including DNA replication and oxidative phosphorylation signaling. These are all tumor-related pathways, suggesting a poor prognosis for the low m1Ascore group (Figure 5G).

Next, we deeply analyze the value of m1Ascore in the prognosis of OC patients. High m1Ascore group patients showed better survival benefits (Figure 6A). Similarly, the high m1Ascore group consistently showed a marked survival advantage in patients stratified by different clinical characteristics (Figure 6B). Moreover, high m1Ascore group had a significantly greater proportion of early-stage patients and recurrence-free population (Figure 6C), and these patients had higher m1Ascore (P<0.05, Figure 6D). This means that high scoring populations characterized by m1Acluster-B modification modes and immune activation phenotypes have a more favorable prognostic outlook. The above outcomes indicate that m1Ascore can also be employed to assess several clinical features of OC patients, including clinical stage, grade, and survival status.

There is a statistics association between tumor burden mutations (TMB) and tumor stage, grade, residual tumor size, and immune infiltrating cells in TME, which indicates that TMB has a major role in predicting OC survival and guiding immunotherapy in OC patients (48). Given TMB’s clinically important nature, we attempted to probe the inherent relevance between TMB and m1Ascore to clarify genuine markings of the two groups. First, as shown in Figure 7A, although not statistically significant, low m1Ascore subgroup patients showed greater TMB than high m1Ascore group. Correlation analysis also confirmed that m1Ascore is significantly negatively associated with TMB (Spearman coefficient: R = -0.18, p = 0.0051; Figure 7B). Then, patients were dichotomized into high and low TMB groups. As shown in Figure 7C, we observed that OS was better in high-TMB patients (p=0.016). Collectively, these findings indicate that m1Ascore can be treated a predictor independent of TMB. Besides, we use maftools package to compare the distribution of somatic variation between low and high m1Ascore groups in TCGA. Figure 7D illustrates that low m1Ascore group exhibited a more extensive TMB than the high m1Ascore subgroup. These results provide new insights into the mechanisms underlying tumor m1Ascore composition and gene mutation, and indirectly confirm the value of m1Ascore in predicting the outcome of immunotherapy

In recent years, platinum-based drugs and paclitaxel have been used as representative drugs for chemotherapy of OC (49). For this reason, based on the two treatment cohorts of TCGA and GSE9891, we explored whether m1A modification characteristics can predict the response of patients to these two first-line chemotherapeutic agents, including cisplatin and paclitaxel. It can be seen from Figure 8A that patients with low m1Ascore showed significant treatment sensitivity to both groups of drugs, which indicates that the m1Ascore facilitates the forecasting of patient reaction to chemotherapy.

Moreover, we utilized CMap database to examine the potential role of m1Ascore in the development of new drugs for OC. First, we screened for differentially expressed genes between the groups with high and low m1Ascore (Supplementary Figures 5A, B) and conducted GO and KEGG enrichment analysis on these genes (Supplementary Figures 5C, D). Subsequently, we utilized CMap database to analyze these genes and find drugs with high association with OC. The top 10 small molecule drugs highly correlated with OC are shown in Table 1. Among these agents, we obtained chemical structure information from PubChem for 4 most important small molecule agents (Figures 9A–D).

ICI has emerged as a key drug for immunotherapy (23). However, there is still several patients who do not respond to immunotherapy, which to some extent limits the application of ICIs. Therefore, Charoentong et al. developed a quantitative scoring scheme called immune phenotype score (IPS) to ascertain the determinants that influence tumor immunogenicity. In this scoring scheme, IPS is an excellent predictor for detecting anti PD-1 and anti CTLA-4 antibody responses (44). Here, we completely analyzed the link that exists between IPS and immune characteristics. The IPS-CTLA4-neg-PD-1-neg, IPS-CTLA4-neg-PD-1-pos, IPS-CTLA4-pos-PD-1-neg, IPS-CTLA4-pos-PD-1-pos scores were designed to assess the likelihood of patients receiving ICIs therapies. The results proved that the score increased significantly in the high m1Ascore group (Figure 8B): IPS-CTLA4-neg-PD-1-neg, P = 0.023; IPS-CTLA4-neg-PD-1-pos, P = 0.036. This indicates that the high m1Ascore group shows higher IPS and appears to have more immunogenic phenotypes. These results suggest that the effectiveness of ICIs may be better in patients with high m1Ascore.

Besides, we also investigated some common immune molecules expressions in different m1Ascore groups, including CD274, PDCD1, CD40, CXCL9 and so on. We observed that most genes, except CD70, showed higher expression levels in the low m1Ascore group (Figure 8C). This may also be a major factor in the lower survival rate of low m1Ascore patients (Figure 6A). These results above imply that quantification of m1A modification patterns might be used as a prospective and stable biomarker for chemotherapeutic response and immunotherapy evaluation.

We used the HPA database to analyze m1A regulator protein expression. Since immunohistochemical data for TRMT61A, YTHDF and YTHDF3 were missing from the HPA data, we only analyzed the remaining seven m1A regulators. Compared to the normal tissue, TRMT10C, TRMT6 and YTHDF2 were relatively highly expressed in the tumor tissue, while YTHDC1 was relatively less expressed. TRMT61B was moderately expressed in both tumor and normal samples. ALKBH1 and ALKBH3 were both lowly expressed in tumor and normal samples (Figures 10A–G).

In terms of mRNA levels, TRMT6, TRMT61A, TRMT61B, TRMT10C, YTHDF1 and YTHDF2 YTHDF3 were overexpressed in tumor sample. While ALKBH1 and YTHDC1 were overexpressed in normal sample. ALKBH3 expression had no difference between the tumor and normal sample (Figures 11A–J). The mRNA expression detected by PCR was in general agreement with the TCGA database (Figure 1D).