Common gene expression data and complete clinical records of BC patients were retrieved from TCGA and the GEO database (http://www.ncbi.nlm.nih.gov/geo/) (Xin et al., 2020). The BC patients without survival information were excluded. In this study, a total of six eligible BC GEO cohorts [GSE21653 (Sabatier et al., 2011), GSE24450 (Heikkinen et al., 2011), GSE42568 (Clarke et al., 2013), GSE45255 (Nagalla et al., 2013), GSE51783 (Zhao et al., 2014), and GSE61304 (Grinchuk et al., 2015)] with complete prognostic information and mRNA sequence information and The Cancer Genome Atlas-BC (TCGA-BC) were included together for further analysis (Supplementary Table S1). The averaging method of Affy and simpleafty was used to adjust background and quantile normalization. Our research flowchart is shown in Supplementary Figure S1. We downloaded the RNA sequencing data (FPKM value) from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) and used the R package TCGAbiolinks for comprehensive analysis. TCGA bio links package was a software package specially developed for GDC data analysis (Colaprico et al., 2016). The somatic mutation dataset was also from the TCGA database, and the R (version 3.6.3) and R biological conductor packages were used to analyze the data.

A total of 23 m6A RNA methylation regulators were selected from the published literature (Yang et al., 2018; Shi et al., 2019). Next, based on the expression levels of 23 m6A regulators, we used unsupervised cluster analysis to classify the patients into three m6A modification modes with optimal k value, used the ConsensClusterPlus R software package to cluster the patients, and performed cycle computation 1,000 for guaranteeing the stability of classification (Wilkerson and Hayes, 2010). The overall survival (OS) between different clusters was calculated by the Kaplan–Meier statistical method.

GSVA is a non-parametric and unsupervised method commonly used to estimate the path variation and the activity change of biological processes in samples of expression datasets (Hänzelmann et al., 2013). Then, the “GSVA” R package was used to analyze the GSVA enrichment of genes and studied the enrichment of biological processes between different m6A modification modes. We downloaded the gene sets of c2. cp. KEGG. v6.2.—symbols from the MSigDB database for GSVA analysis after correction. When the corrected p value was less than 0.05, it was considered statistically significant. The clusterProfiler R package was used to annotate the functions of m6A-related genes, and the critical value of false discovery rate (FDR) was <0.05.

To identify genes associated with m6A modification patterns, patients were divided into different m6A clusters based on the expression levels of 23 m6A regulators. Then, the empirical Bayesian method of the limma R package was used to identify DEGs in three different m6A modification patterns (Ritchie et al., 2015). The significance standard of DEGs was |logFC|> 1, p < 0.05.

To understand the degree of immune cell infiltration in the three subgroups, we quantified the relative abundance of each cell infiltration in the TME of the BC samples using single-sample gene-set enrichment analysis (ssGSEA). Then, we obtained the gene sets for each type of TME infiltrating immune cells and stored a relatively complete subset of human immune cells, including activated CD8 T cells, natural killer T cells, activated DCs, macrophages, and regulatory T cells (Supplementary Table S2). The relative abundance of different infiltrating cells in the TME in each sample was evaluated using the enrichment scores calculated by ssGSEA analysis.

In this study, we applied a method to quantify the m6A modification pattern of individual BC patients. Thus, we constructed a scoring system to evaluate the m6A modification pattern of BC patients, namely, m6Ascore. The specific procedures were as follows.

Firstly, the overlapping DEGs in different m6Aclusters were extracted from all BC samples. These overlapping DEGs were analyzed using unsupervised clustering, and patients were divided into several groups for further analysis. The consistency clustering algorithm was used to determine the number and stability of gene clustering. Next, we used the univariate Cox regression model to analyze the prognosis of each overlapping DEG and selected genes with significant prognostic differences for principal component analysis (PCA) to construct the m6A gene signature. Principal component 1 and principal component 2 were selected as signature scores, and the scores were concentrated on the set of gene blocks with the greatest correlation. In contrast, the contributions of other genes unrelated to other set members were weighted. Finally, m6Ascore was determined in a similar way to the Genome Grading Index (GGI) (Sotiriou et al., 2006; Zhang et al., 2020): m 6 A s c o r e = ∑ ( P C 1 i + P C 2 i ) , where i is the expression of overlapping DEGs with a significant difference in prognosis among m6Aclusters.

We used the gene set to assess the association between the m6Ascore and important biological processes (Şenbabaoğlu et al., 2016; Zeng et al., 2019), including 1) epithelial-mesenchymal transition (EMT) markers, EMT1, EMT2, and EMT3; 2) immune checkpoint; 3) CD8 T-effector signature; 4) mismatch repair; 5) antigen processing and presentation; 6) antigen processing machinery (APM); 7) Wnt targets; 8) DNA damage repair; 9) nucleotide excision repair; 10) pan-fibroblast transforming growth factor-β response signature (Pan-FTBRS); 11) DNA replication; and 12) angiogenesis signature.

The Kruskal–Wallis test and one-way ANOVA were used to evaluate the statistical differences among three or more groups (Mariathasan et al., 2018). Spearman’s correlation analysis was used to calculate the correlation coefficient between the expression of m6A regulators and TME infiltrating immune cells. Based on the correlation between the m6Ascore and the patient prognosis, the optimal cut-off point of the data subgroup was determined by the survminer R software package, and the patients were divided into high and low m6Ascore. The Kaplan–Meier method was used to plot the survival curve for the prognosis analysis. The log-rank test was used to determine the statistical significance of the difference. The multivariate Cox regression analysis was used to evaluate the independent prognostic factors. The mutation landscape of the high and low m6Ascore group was presented by the waterfall function of the R software maftools package (Hazra and Gogtay, 2016). The RCircos R package was used to describe the CNV of 23 m6A regulators in 23 pairs of chromosomes. In this study, all tests were bilateral, and p < 0.05 was considered statistically significant. All data processing was carried out in R software V3.6.3.

Firstly, a total of 23 genes that regulate RNA methylation were identified, including eight writers (METTL14, METTL3, RBM15/RBM15B, CBLL1, WTAP, ZC3H13, and VIRMA), two erasers (ALKBH5 and FTO), and 13 readers (FMR1, YTHDF1/2/3, HNRNPA2B1, YTHDC1/2, RBMX, HNRNPC, LRPPPRC, and IGF2BP1/2/3). The dynamic process and mechanism of m6A methylation in RNA modification are described in Figure 1A. Next, we integrated somatic mutation and CNV to explore the mutation rate of 23 m6A regulators in BC. The overall mutation frequency of the m6A regulator was relatively low, and only 55 in 983 samples caused m6A regulator mutation (Figure 1B). Among 23 m6A regulators, LRPPRC, YTHDF1, FMR1, WTAP, YTHDC1, YTHDF3, and RNPA2B1 were mutated in BC, while the other regulatory factors were not mutated. Furthermore, it was found that the CNV alteration frequency of 23 m6A regulators was prevalent in BC, most of which were with the copy number amplification. However, WTAP, RBM15, ZC3H13, YTHDC2, YTHDF2, and RBM15B were with the high frequency of CNV deletion (Figure 1C). The CNV alteration location of the m6A regulators on the chromosome is presented in Figure 1D. Based on the expression levels of these 23 m6A methylation-related regulators, BC samples could be markedly distinguished from normal samples (Figure 1E). Then, we investigated the mRNA expression levels of m6A regulators in BC and normal tissues to determine whether the abnormal expressions were associated with BC malignancy. Compared with normal breast tissue, 17 of the 23 m6A regulators were differentially expressed (Figures 1C,F). These results indicated that the mRNA expression levels of m6A regulators were highly heterogeneous in BC and normal tissues, demonstrating the underlying roles of the abnormal expression pattern of m6A regulators in the oncogenesis and progress of BC.

Next, the clinical significance of the m6A regulators in BC patients was explored. The network described the interaction between the m6A regulators and their impact on the BC prognosis (Figure 2A, Supplementary Table S3). This result denoted that some m6A regulators (such as YTHDF1 and FTO) were associated with a poor prognosis. However, other regulatory factors (such as HNRNPC and IGFBR3) were associated with a good prognosis of BC patients (Supplementary Figures S2–S4). Moreover, the expression of m6A regulators of the same functional category also showed a significant correlation. We then classified the patients with different m6A modification patterns based on the expression levels of 23 m6A regulators. By model-based cluster analysis, three different methylation modification patterns were identified, including 685 cases in m6ACluster A, 617 cases in m6ACluster B, and 498 cases in m6ACluster C (Figure 2B, Supplementary Table S4). Prognostic analysis of three major subtypes of m6A modification showed that the m6ACluster B modification pattern had a particularly significant survival advantage (Figure 2C). Moreover, the relationship between clinical factors and gene expression of three different m6A methylation patterns was analyzed and presented in Figure 2D.

Through GSVA enrichment analysis, we then intended to investigate the enrichment of biological processes in these different m6A modification patterns. m6ACluster A was significantly correlated with immunosuppression and other biological processes. m6ACluster B was significantly enriched in the pathways related to complete immune activation, including Toll_like_receptor_signaling_pathway, T_cell_receptor_signaling_pathway, B_cell_receptor_signaling_pathway, and Chemokine_signaling_pathway. Moreover, m6Acluster C showed significant enrichment in carcinogenic activation and matrix pathway, such as TGF_beta_signaling_pathway, Adherens_junction, and ECM_receptor_interaction (Figures 3A,B, Supplementary Table S5). PCA was used to analyze the transcriptome profiles of the three m6A modification patterns, showing significant differences in the transcriptome profiles of different modification patterns (Figure 3C). In the subsequent analysis of TME cell infiltration, it was found that the m6ACluster C was enriched in the infiltration of innate immune cells, such as macrophages, myeloid-derived suppressor cells, natural killer cells, monocytes, regulatory T cells, CD8 T cells, and DCs (Figure 3D, Supplementary Table S6). However, according to the results of our survival analysis, the patients with m6ACluster C did not possess the corresponding survival advantages. It has been proven that tumors with immune rejection phenotype were characterized by the presence of a large number of immune cells, which stayed in the matrix surrounding the tumor cell nest and did not penetrate the central zone (Zhang et al., 2020).

In order to explore the potential biological regulatory pathways in the different m6A modification patterns, we successfully identified 3429 DEGs associated with the m6A phenotype (Figure 4A). The clusterProfiler software package was used to analyze the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) function enrichment of the DEGs. The results are shown in Figures 4B–E. The KEGG analysis showed that the DEGs were enriched in the cell cycle, EGFR tyrosine kinase inhibitor resistance, and mTOR signaling pathway (Figures 4B,C). The GO analysis of the biological process showed that these DEGs were enriched in DNA replication, regulation of DNA metabolism, and histone modification (Figures 4D,E). Furthermore, the cellular component analysis showed that DEGs were abundant in histone methyltransferase complex and transferase complex and transferring phosphorus-containing groups. Molecular function analysis indicated that DEGs were mainly located in the helicase activity and transcription coactivator activity (Figures 4D,E).

To further understand this regulatory mechanism, we divided the patients into different genotypes by unsupervised cluster analysis of the 3429 genes related to the phenotype of m6A, using the R limma package from the transcriptomic profile of TCGA and GEO. By model-based cluster analysis, we finally identified three different methylation modification patterns, including 685 cases in gene cluster A, 617 cases in gene cluster B, and 498 cases in gene cluster C (Figure 5A). The results were consistent with the cluster grouping of m6A modification patterns, revealing three different m6A modified genomic phenotypes, named m6A gene cluster A, cluster B, and cluster C (Supplementary Figures S4A–D, Figure 5B, Supplementary Table S4). This indicates that these three different methylation patterns of m6A did exist in BC. There were significant differences in the expression of m6A regulators in the three m6A gene clusters, which was consistent with the expected results of previous methylation modification patterns of m6A (Figure 5C). Furthermore, there was an observable good prognosis in gene cluster B and a poor prognosis in gene cluster C (Figure 5D). However, the above could only be based on the analysis of the patient population and could not accurately predict the methylation pattern of m6A in individual patients.

Due to the individual differences and complexity of m6A methylation modification, we constructed a set of m6A modification models that could quantify the individual BC, namely, the m6Ascore. The changes in individual patient attributes are shown by the alluvial chart (Figure 6C). The optimal cut-off value of 7.795803 was calculated using the survminer software package, and patients were successfully divided into a high group and low group according to m6Ascore (Figure 6A). These results indicated that patients with a high score of the m6A group had a significant survival benefit. Then, the correlation between the m6Ascore and biological process was analyzed to clarify the genetic characteristics of m6A (Figure 6B). Besides, the Kruskal–Wallis test revealed a significant difference in m6Ascore among m6A gene clusters (Figure 6D). Gene cluster C showed the lowest median score, while gene cluster B showed the highest median score, indicating that the low m6Ascore group might be closely related to immune activation-related signatures, whereas the high m6Ascore group could be related to stromal activation-related signatures (Figure 6D). The score of m6Acluster B was higher than that of m6Aclusters A and C, and the m6Ascore of m6Acluster C was the lowest (Figure 6E).

Then, we used the maftools software package to analyze the distribution of somatic mutations between high and low m6Ascore of BC patients in the TCGA database (Figures 7A,B). The results showed that the mutation rate of PIK3CA was relatively high in the high m6Ascore cohort compared to the low m6Ascore cohort (37 vs. 27%). Subsequently, quantitative analysis of tumor mutational burden (TMB) confirmed that the m6Ascore was significantly negatively correlated with TMB (Figure 7C). The patients with low m6Ascore were significantly associated with the higher TMB (Figure 7D). At the same time, the m6Ascore analysis showed that the low M6Score group had a higher proportion of patient death (Figure 7E) and the average m6Ascore was higher in alive patients than dead patients (Figure 7F). Survival analysis of patients showed that patients with low TMB had a more significant prognostic advantage than high TMB patients (Figure 7G). To predict the molecular subtypes of BC samples from TCGA, we explored the proportion of basal, her2, luminal A, luminal B, and normal patients in the low and high m6Ascore groups. It was found that the luminal A type accounted for a large proportion in the high m6Ascore group, while the basal type accounted for a large proportion in the low m6Ascore group (Supplementary Figure S6A). Then, according to our established score, the basal subtype significantly exhibited the lowest m6Ascore among these subtypes, consistent with the overall prognosis of BC patients (Supplementary Figure S6B). The above results were consistent with the previous conclusion (Figure 6A).

Clinically, the PD-L1 expression of immune cells was considered separately for PD-L1 immune therapy, and the immune cell proportion score (IPS) was introduced as an indicator to distinguish beneficiaries. Furthermore, we also explored the association between IPS and the m6Ascore in BC. The IPS-PD1, IPS-CTLA4, and IPS-PD1/CTLA4 scores were significantly increased in the high m6Ascore group (Figures 8A–D). The treatment benefits from antibodies against immune checkpoints of PD-1 and CTLA-4 might differ between the high and low m6Ascore groups. There is evidence that patients with high TMB status have a relatively durable clinical response to PD-1/PD-L1 immunotherapy (Luchini et al., 2019). The above analyses proved that the m6A modification mode is intensively associated with TMB and PD-1/PD-L1 immunotherapy.