The workflow of this study was shown in Figure S1 . The expression profile data and clinical information of pancreatic cancer samples were obtained from the Cancer Genome Atlas (TCGA, RRID : SCR_003193, https://tcga-data.nci.nih.gov/tcga/) and Gene-Expression Omnibus (GEO, RRID : SCR_005012, https://www.ncbi.nlm.nih.gov/geo/) database. This study collected 5 independent PC cohorts (TCGA-PAAD, GSE28735, GSE57495, GSE62452, GSE85916) for further analysis. RNA sequencing data in TCGA (FPKM format) were downloaded and transformed into TPMs (transcripts per kilobase million). We download the normalized matrix file in GEO, and applied ComBat algorithm in the R package SVA to eliminate batch effects between different GEO data sets. The survival status and survival time of the samples were extracted from the clinical information of the 5 PC cohorts. Data with follow-up time less than 31 days and duplicate data were excluded. Somatic mutation data was collected from the TCGA database. The copy number variation data (CNV) of TCGA-PAAD was downloaded from the UCSC Xena database (http://xena.ucsc.edu/).

We searched the relevant literature on m⁶A methylation modification to identify recognized m⁶A regulators for subsequent analysis. A total of 23 m⁶A regulators were included, including 8 writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B), 13 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, and RBMX), and 2 erasers (FTO and ALKBH5). Unsupervised clustering analysis was performed to identify different m⁶A methylation modification subtypes according to the expression of 23 m⁶A regulators, and patients were classified for further analysis. The number of clusters (K) and their stability were determined by the consensus clustering algorithm. Li et al. (31) tested four methods for finding K: the cumulative distribution function (CDF), the proportional change in the area under the CDF curve upon an increase of K (Δ(K)), GAP-PC (32) and CLEST (33). They found that CDF was able to reveal the correct K, as the CDF curve was flat only for the true K, reflecting a perfectly or near-perfectly stable partitioning of the samples at the correct K. So we took the K corresponding to the flattest CDF curve as the determined number of clusters. The R package consusclusterplus was utilized to perform the above steps (34).

The R package CIBERSORT was used to quantify the infiltration of different immune cells in PC samples from five cohorts. Leukocyte signature matrix (LM22) contains 547 reference genes, which can be used to distinguish 22 human immune cell phenotypes, including various types of T cells, B cells, NK cells, plasma cells, and myeloid subgroups. CIBERSORT is a deconvolution algorithm, which can calculate the proportion of different types of cells in the sample based on LM22 (35). The ESTIMATE algorithm infers the cell density and tumor purity of the tumor based on the transcriptome profile of the sample (36). Tumor tissue with rich immune cell infiltration indicates a higher immune score and lower tumor purity. The R package ESTIMATE was used to evaluate the immune and stromal content (immune and stromal score) in each sample.

Previous consensus clustering algorithms divided patients into two different m⁶A modification subtypes based on the expression of 23 m⁶A regulators. R package Limma was used to identify DEGs between the two m⁶A modification clusters with adjusted P value < 0.05.

GO (Gene Ontology) is a crucial bioinformatics tool for annotating and analyzing the biological functions of genes, including MF (molecular function), BP (biological processes), and CC (cellular components). As a database resource, KEGG (Kyoto Encyclopedia of Genes and Genomes) is mainly used to understand the high-level functions and values of biological systems from molecular-level information. To get annotation information and explore the biological functions of the above DEGs, GO and KEGG enrichment analyses were accomplished using clusterProfiler package with a cutoff of p-value < 0.05 and q-value < 0.05.

To quantify the modification pattern of m⁶A in individual PC patients, we used principal component analysis (PCA) to construct the m⁶A scoring scheme. Firstly, univariate Cox regression model was performed on DEGs identified between different m⁶A modification clusters, and the genes with significant prognosis effect were selected for clustering samples and constructing m⁶A score. The patients were divided into several groups for further analysis. The number and stability of gene clusters were determined by consensus clustering algorithm. Subsequently, principal component analysis was used to construct the m⁶A-related gene signature, and principal component 1 and principal component 2 were extracted as signature scores. This method focuses the score on the set with the largest block of strongly related or anti-related genes, while reducing the contribution of genes that are not tracked with other set members. Besides, the PCA algorithm can effectively achieve data dimensionality reduction and largely retain the information of the original data. We used a method similar to GGI (37, 38) to define the m⁶A score: m 6 A s c o r e = ∑ ( P C 1 i + P C 2 i ) where i is the expression of final genes related to the m⁶A phenotype.

To verify the reliability and clinical application value of m⁶A score, the receiver operating characteristic (ROC) curves of 1-, 3-, and 5-year were drawn. We first drew the ROC curve based on all samples. Then, the ROC curve was drawn solely in the TCGA-PAAD cohort, and the prognostic prediction performance of m⁶A score and other clinical indicators were compared. Univariate and multivariate Cox regression analyses were used to evaluate the correlation between the patient’s m⁶A score, clinical variables, and prognosis to determine whether the score can be used as an independent prognostic indicator of pancreatic cancer. P < 0.05 indicated that the difference was statistically significant. The results were shown in the forest diagram. Next, 8 indicators (age, gender, grade, stage, stage_T, stage_N, stage_M, and m⁶A score) were used to construct a nomogram to personally predict the 1-year, 3-year, and 5-year survival rates of patients. The ROC curve was drawn to show the predictive performance of the nomogram. The R packages survival, survminer, timeROC, rms, and regplot are used for calculation and graph drawing.

We calculated the copy number increase or loss frequency of 23 m⁶A regulators in the TCGA-PAAD cohort, and used the R package Rcircos to draw a copy number variation map of m⁶A regulators on human chromosomes. To determine the tumor mutation burden (TMB), we counted the total number of non-synonymous mutations in the TCGA-PAAD cohort. R package maftools was used to plot the oncoprint of gene mutation.

Immunophenoscore (IPS) is a favorable factor to predict the efficacy of anti-CTLA-4 and anti-PD-1 regimens, which can quantify the determinants of tumor immunogenicity and show the characteristics of tumor immune landscape (39). The score is calculated based on four categories of immune-related genes, including MHC molecules (MHC), immunomodulators (CP), effector cells (EC), and suppressor cells (SC). IPS of samples in TCGA-PAAD were downloaded from the online platform The Cancer Immunome Atlas (TCIA, RRID : SCR_014508, https://tcia.at/home) for further analysis.

All statistical analysis and graph drawing were completed by R-4.0.3. The Wilcoxon rank-sum test was used for comparison between two groups. Kaplan-Meier survival analysis and univariate Cox regression model were performed to calculate the relationship between m⁶A regulators and prognosis. According to the correlation between m⁶A score and patient survival, R package Survminer was used to repeatedly test all possible cut-off points to obtain the largest rank statistic, and the patients were divided into high and low m⁶A score groups based on the largest selected log-rank statistic. The Kaplan-Meier method was utilized to draw the survival curve for prognostic analysis, and the log-rank test was used to determine the significance of the difference. Spearman correlation analysis and distance correlation analysis were applied for the correlation test. All heat maps were generated by R package pheatmap. All statistical P value were two-tailed, and P < 0.05 was statistically significant.

In this study, 23 m⁶A regulators were identified, including 8 writers, 13 readers, and 2 erasers. We first calculated the incidence of somatic mutations in PC. Among 158 tumor samples, a total of 120 cases (75.95%) had genetic alterations, of which TP53 and KRAS mutations were the most frequent with more than 50% ( Figure 1A ). However, the mutation frequency of 23 m⁶A regulators was pretty low, genetic alterations occurred in only 5 (3.16%) samples ( Figure 1B ). We further analyzed the relationship between TP53 and KRAS mutations and the expression of m⁶A regulators. There were differences in the expression of multiple m⁶A regulators between the wild group and mutant group ( Supplementary Figures 3, 4 ). Analysis of copy number variation of 23 m⁶A regulators showed that CNV mutations were common in PC. VIRMA, HNRNPA2B1, IGFBP1, IGFBP3, YTHDF3, and YTHDF1 had widespread CNV amplification. However, METTL16, WTAP, ALKBH5, YTHDF2, and RBM15B showed extensive CNV deletions ( Figure 1C ). The location of CNV changes on the chromosome was illustrated in Figure 1D . Kaplan-Meier survival analysis showed that 15 m⁶A regulators were correlated with the prognosis of PC patients ( Supplementary Figure 2 ). Univariate Cox regression model revealed the prognostic value of 23 m⁶A regulators in PC patients ( Supplementary Table 1 ). As shown in Figure 1E , the network presents a comprehensive landscape of the interactions, connection, and prognostic significance of m⁶A regulators in PC. The results indicated that writers, readers, and erasers have a significant correlation in expression. The cross-talk among them probably plays an essential role in the formation of different m⁶A modification patterns, and might be related to the occurrence and development of cancer. In addition, we implemented CIBERSORT and ESTIMATE algorithms to quantify the activity or enrichment level of immune cells in pancreatic cancer tissues. The correlation coefficient heatmap was used to display a general landscape of immune cell interactions in the tumor microenvironment ( Figure 1F ).

Based on the expression of 23 m⁶A regulators, the R package ConensusClusterPlus was used to qualitatively classify patients with different m⁶A modification subtypes. Through the consensus clustering algorithm, two different m⁶A modification subtypes were finally identified, including 294 cases in subtype A and 140 cases in subtype B ( Figures 2A–D and Supplementary Table 2 ). We named these two subtypes m⁶A cluster A and m⁶A cluster B, showed the expression of 23 m⁶A regulators in the two modified subtypes by heatmap ( Figure 2E ). The expression levels of 23 m⁶A regulators between the two m⁶A clusters were also compared and shown in Figure 2F . To explore the internal biological changes under different m⁶A modification modes, we compared the composition of immune cells in TME. The result showed that m⁶A cluster A was characterized by higher infiltration of memory B cells and activated memory CD4+ T cells. In m⁶A cluster B, the infiltration of activated NK cells, M0 macrophages, and Neutrophils were significantly increased ( Figure 2G ). To reveal the potential biomolecular characteristics of different m⁶A modified phenotypes, R package LIMMA was used for differential expression analysis to determine the transcriptome differences between two subtypes. We identified 1159 differentially expressed genes and annotated DEG with R package clusterProfiler. Figures 2H, I summarized the significant biological processes of DEG enrichment, such as glucose metabolism, glycolysis, HIF-1 signaling pathway, Hippo signaling pathway, and TGF-β signaling pathway. These results suggested that the m⁶A methylation modification may involve in tumor metabolism and immune regulation, and was closely related to tumor genesis and progression. Supplementary Tables 3, 4 provides detailed descriptions.

Although the consensus clustering algorithm based on 23 m⁶A regulators classified PC patients into two subtypes, potential genetic changes and prognostic correlations in these phenotypes were not very clear. We performed univariate COX regression analysis on the 1159 DEGs between the previously identified m⁶A clusters, and obtained 719 survival-related genes which were named m⁶A-related signature genes ( Supplementary Table 5 ). Based on representative m⁶A-related signature genes, we adopted unsupervised cluster analysis and identified two stable transcriptome phenotypes, which were defined as gene cluster A and gene cluster B ( Figures 3A–C ). In addition, we explored the prognostic significance of gene subtypes by integrating transcriptome and survival information. Through Kaplan-Meier analysis and log-rank test, gene cluster B showed a better prognosis (P < 0.001, Figure 3D ). The heatmap showed the transcriptome profile of 719 m⁶A-related signature genes in two gene clusters ( Figure 3E ). The expression levels of 23 m⁶A regulators between the two m⁶A-related gene clusters were also compared. Significant differences in the expression of m⁶A regulators were observed, which was consistent with the expected result of m⁶A modification patterns ( Figure 3F ). Previous studies have shown that the immune system may produce favorable or unfavorable consequences, which may manifest as pro-tumor or anti-tumor activity. Monocytes and anti-tumor lymphocyte subsets such as CD8+ T cells, memory CD4+ T cells, and naive B cells had a higher level of infiltration in gene cluster B, while activated NK cells and mast cells infiltrated more abundantly in gene cluster A ( Figure 3G ). The immune and stromal score based on the ESTIMATE algorithm indicated that the infiltration of immune cells and stromal components in gene cluster B was higher. Therefore, we speculate that the abundant immune cell infiltration in gene cluster B forms an effective anti-tumor immune response.

Although the above results demonstrated the role of m⁶A methylation modification in the regulation of immune cell infiltration and prognosis, these analyses cannot accurately predict the m⁶A methylation modification pattern in a single tumor patient. To obtain a quantitative index of the m⁶A modification landscape of PC patients, we extracted the scores of principal component 1 and principal component 2 for calculating the final m⁶A score. Figure 4A showed that patients’ m⁶A score in m⁶A cluster A was lower than m⁶A cluster B, and Figure 4B depicted that the m⁶A score in gene cluster B was lower than that in gene cluster A. We drew an alluvial diagram to display the process of m⁶A score construction ( Figure 4C ). Subsequent analysis revealed the prognostic significance of m⁶A score. According to Figure 4D , the patients’ survival rate (41% vs. 23%) in the m⁶A low score group was much higher than that in the high score group. The overall average m⁶A score of the surviving patients was lower than that of the dead patients ( Figure 4E , P = 0.0025). Kaplan-Meier survival analysis showed that the prognosis of patients in the low m⁶A score group was significantly better ( Figure 4F , P < 0.001).

To verify the m⁶A score, the 1-, 3-, and 5-year ROC curves were drawn, and the value of the area under curve (AUC) of the m⁶A score was calculated. The results showed that the AUC values of all three curves were around 0.65 both in total samples ( Figure 5A ) and TCGA-PAAD cohort ( Figure 5B ). We also compared the 1-year ROC curve with other clinical characteristics in TCGA-PAAD cohort, and m⁶A score had the most considerable AUC value ( Figure 5C ). Univariate Cox regression analysis showed that age (p = 0.012, HR = 1.027, 95%CI [1.006-1.049]), grade (p = 0.026, HR = 1.392, 95%CI [1.041-1.862]) and m⁶A score (p = 0.002, HR = 1.022, 95%CI [1.008-1.037]) were considered statistically significant ( Figure 5D ). Multivariate Cox regression analysis showed age (p = 0.012, HR = 1.028, 95%CI [1.006-1.050]) and m⁶A score (p = 0.005, HR = 1.021, 95%CI [1.006-1.036]) were independent prognostic predictors ( Figure 5E ). By integrating multiple clinical indicators, the nomogram can be an effective tool for quantitatively assessing individual risks in the clinical environment. We constructed a nomogram to predict patients’ OS at 1-, 3-, and 5-year. Taking a random sample as an example, the total score of the patient was 340, the probability of survival time less than 1-, 3-, and 5-year were 0.123, 0.452, and 0.563, respectively ( Figure 5F ). The ROC curve was used to evaluate the predictive performance of the nomogram. The AUC values of 1-, 3-, and 5-year ROC curves were 0.718, 0.800, and 0.792, respectively ( Figure 5G ). The results showed that m⁶A score could be used as a new effective clinical predictor and can be combined with other clinical variables to improve the prognosis of patients with pancreatic cancer.

Previous studies have shown that tumor mutation burden may be an emerging and potential tumor marker, which can assist in the selection of patients for immune checkpoint therapy. In view of the important clinical significance of TMB, we tried to explore the inner link between TMB and m⁶A score to clarify the genetic imprint of each m⁶A score subgroup. Correlation analysis showed a significant and positive association between m⁶A score and TBM (Spearman coefficient: R = 0.18, P = 0.032; Figure 6A ). Next, we divided patients into two subgroups based on TBM. As shown in Figure 6B , we found that patients with low TMB showed better overall survival than patients with high TMB. Next, we evaluated the synergy of these scores in the prognostic stratification of PC. Survival analysis demonstrated that TBM status did not affect predictions based on m⁶A score, and the low m⁶A score group always showed a survival advantage ( Figure 6C ). In addition, we analyzed the somatic mutation landscape in the high and low m⁶A score groups, and found that the high m⁶A score group had a higher mutation rate (93.48%) than the low score group (70.65%). According to the results ( Figures 6D, E ), both KRAS (74% vs. 46%) and TP53 (78% vs. 43%) had a higher somatic mutation rate in the high m⁶A score group, which may be related to the poor prognosis of high m⁶A score group. These data can more comprehensively describe the impact of m⁶A score on genomic variation, and may provide new ideas for studying the potential interaction between m⁶A methylation modification and somatic mutation.

The treatment of immune checkpoint inhibitors represented by CTLA-4/PD-1 inhibitors is undoubtedly a major progress in anti-tumor therapy. We compared the expression of common immune checkpoint genes between the high and low m⁶A score groups, and found that PD-L1 was highly expressed in the high m⁶A score group, PDCD1 was highly expressed in the low m⁶A score group, while the expression of CTLA4 and IDO1 had no significant difference between the two groups ( Figures 7A–D ). As a new predictor of the immune response, IPS is widely used and recommended for evaluating the immune response of patients. Our analysis showed that the IPS of the low m⁶A score group was higher no matter in the case of anti-PD-1/CTLA-4 therapy alone, or combination therapy ( Figures 7E–H ). These results indicated that m⁶A methylation modification in pancreatic cancer may play an important role in mediating the immune response.