The RNA sequencing (RNA-seq) transcriptome data and corresponding clinicopathological features of PAAD samples were obtained from The Cancer Genome Atlas (TCGA) database in April 2021. Gene expression data (measured in fragments per kilobase of exon per million fragments mapped or FPKM) was transformed into transcripts per kilobase million (TPM). Simple nucleotide variation data was extracted from the TCGA database, while the copy number variation (CNV) data was obtained from the UCSC Xena Website (https://xena.ucsc.edu/). Sample differences in the tumor microenvironment (TME) were measured using Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) analysis with the “estimate” R package (Yoshihara et al., 2013). In addition, an eligible PAAD cohort (GSE21501) was downloaded from the Gene Expression Omnibus (GEO) database. In subsequent analysis, TCGA and GEO datasets were selected as training and validation sets, respectively.

A total of 22 m6A regulators were collected from relevant studies (Li et al., 2020a; Yi et al., 2020; Zheng et al., 2021); the regulators included 7 “writers” (WTAP, METTL16, VIRMA, RBM15B, METTL3, RBM15, and ZC3H13), 13 “readers” (YTHDC1, YTHDF1, YTHDC2, YTHDF3, IGF2BP2, LRPPRC, YTHDF2, HNRNPA2B1, HNRNPC, RBMX, EIF3A, G3BP1 and FXR1), and 2 “erasers” (ALKBH5 and FTO). The prognostic role of the m6A regulators was assessed using the Kaplan-Meier (KM) diagrams and Cox proportional hazards model.

Based on the expression matrix of 22 m6A regulators, patients in the TCGA cohort were classified into distinct clusters according to the best cutoff using the “ConsensusClusterPlus” R package (Wilkerson and Hayes, 2010). The number of clusters and their stability were confirmed by the consensus clustering algorithm. Survival analysis between distinct clusters was measured using the KM method. Differences in the biological processes between the distinct clusters were investigated through gene set variation analysis (GSVA) (Hänzelmann et al., 2013). The “c2.cp.kegg.v7.4.symbols” gene set was obtained from the Molecular Signatures Database (MSigDB). Adjusted p-value<0.05 was considered statistically significant.

The single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm was used to quantify the relative abundance of various immune cell subtypes in PAAD samples (Charoentong et al., 2017). Through enrichment score calculated by ssGSEA, the relative abundance of each immune cell type was represented in each sample. ESTIMATE analysis was performed to compare the differences in the TME with the “estimate” R package (Yoshihara et al., 2013). Furthermore, differences in the TIME and the expression of targeted immune checkpoint molecules between the distinct clusters were compared using the “limma” R package.

Principal component analysis (PCA) was used to test whether m6A regulators could separate distinct m6A modification patterns. Differentially expressed genes (DEGs) among the m6Aclusters were identified using the empirical Bayesian approach with the “limma” R package. The significance criterion of DEGs was set as p-value < 0.0001. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) biological processes analysis were performed to investigate the pathway signatures of the DEGs. A critical value of adjusted p-value = 0.05 was selected as the filter criteria. After identifying the DEGs, prognosis-related genes were filtrated from the DEGs by univariate Cox regression analysis. The significance criterion was set as p-value < 0.001.

Based on the expression of prognosis-related DEGs, samples in the TCGA cohort were classified into different subtypes according to the best cutoff using the “ConsensusClusterPlus” R package. The KM method was used to perform survival analysis between different subtypes. A heatmap revealed the expression of prognosis-related DEGs between different subtypes using the “pheatmap” R package. Furthermore, differences in immune cell infiltration, ESTIMATE score and the expression of targeted immune checkpoint molecules were compared with the “limma” R package.

Based on the expression of prognosis-related DEGs, PCA was used to score the samples in the TCGA and GEO cohorts. Principal components 1 and 2 were used to act as signature scores. The m6Ascore was defined using a method similar to Genomic Grade Index (GGI) (Sotiriou et al., 2006; Zeng et al., 2019): m 6 Ascore = ∑ ( PC 1   i + PC 2   i ) Where i is the expression of overlapping genes with significant prognosis-related DEGs the m6Aclusters.

According to the score, samples were divided into low- and high-m6Ascore groups. We then compared the biological differences between the low- and high-m6Ascore groups, including (O'Reilly et al., 2018) survival analysis (Owens et al., 2019), immuno-correlation analysis (Bilimoria et al., 2007), clinical-correlation analysis (Singhi et al., 2019), tumor mutation burden (TMB) (Mizrahi et al., 2020), targeted immune checkpoint molecules.

All statistical analyses in the study were performed using the R software (version 4.0.5). The Kruskal-Wallis test was used to perform difference comparison on three or more groups (Hazra and Gogtay, 2016). Continuous variables were dichotomized for patient survival using the optimal cutoff values determined by “survminer” R package. The KM and log-rank tests were used to evaluate the survival difference among different clusters. The CNV landscape of 22 m6A regulators in 23 pairs of chromosomes was plotted using the “RCircos” R package. The receiver operating characteristic (ROC) curves (R package “timeROC”) and the area under the curve (AUC) values were used to evaluate the prognostic ability of the m6Ascore for 1-, 2-, 3-, and 4-year overall survival (OS) (Blanche et al., 2013). Univariate and multivariate independent prognostic analyses were performed to assess whether the model was an independent prognostic factor of PAAD. All statistical p values were two-sided, with p < 0.05 deemed statistically significant.

A total of 22 m6A regulators (7 “writers,” 13 “readers,” and 2 “erasers”) were collected in this study. CNVs and somatic mutations were integrated to explore the prevalence of m6A regulator variations in PAAD. The CNV incidence of 22 m6A regulators are shown in Figure 1A. Most regulators focused on the deletion of copy number, while VIRMA, G3BP1, and other five m6A regulators had a prevalent frequency of CNV amplification. The CNV landscape of m6A regulators in 23 pairs of chromosomes are shown in Figure 1B. The overall average mutation frequency of m6A regulators was low, with only 6 (3.8%) of 158 samples having m6A regulator mutation (Figure 1C).

The prognostic value of the m6A regulators was evaluated by KM method and univariate Cox regression analysis (Supplementary Figure S1; Supplementary Table S1). The results showed that most m6A regulators were associated with survival. The network of m6A regulators comprehensively demonstrated the m6A regulators’ interactions, connection, and prognostic significance for PAAD patients (Figure 1D). The results showed that there was a distinct positive correlation between each other. Most regulators, such as IGF2BP2 and HNRNPC, presented tumorigenic characteristics, with higher gene expression levels correlating with poor prognosis. Conversely, several m6A regulators, such as ALKBH5 and METTL16, presented tumor-suppressing characteristics, with higher gene expression levels relating to favorable prognosis. These results suggested that the interrelations among regulators may have important effects on the development and progression of PAAD.

Based on the expression of 22 m6A regulators, model-based clustering was performed to classify PAAD patients using the “ConsensusClusterPlus” R package. Through unsupervised clustering, three different m6A modification patterns were uncovered ultimately (identified as m6Aclusters A-C), including 30 samples in cluster A, 42 samples in cluster B, and 105 samples in cluster C (Figure 2A). Prognostic analysis showed there was a survival disadvantage in m6Acluster B (Figure 2B). The heatmap showed m6Acluster A presented significantly low expression levels of all m6A regulators, while m6Acluster B was characterized by the high expression levels of all m6A regulators (Figure 2C). Moreover, GSVA showed different biological behaviors between the m6Aclusters (Figures 2D–F). The results suggested that different m6A modifications had significant correlation with biological behaviors and prognosis of PAAD.

Using ssGSEA, the study analyzed 23 different immune cell types in the m6Aclusters. The result revealed that m6Acluster B, which had a poor prognosis, was enriched in activated NK cells, mast cells, and T helper type 2 (Th2) cells. However, the abundance of CD56dim NK cells was enriched in m6Aclusters A and C (Figure 3A). These results indicated that m6A modification was associated with the infiltration of specific immune cell types and influenced the response to immunotherapy. In addition, the results of the ESTIMATE algorithm revealed that the stromal and ESTIMATE scores (p < 0.05) were higher in cluster B than in clusters A and C (Figures 3B–D). Combined with the heatmap, the study found that the expression level of m6A regulators showed a similar trend with the ESTIMATE score. Characterized by the higher expression levels of m6A regulators, m6Acluster B also had a higher ESTIMATE score. The results suggested that m6A regulators may play an important role in the regulation of the TME, thus affecting tumor progression and survival. Furthermore, the expression of targeted immune checkpoint molecules was different between the distinct clusters. As shown in the boxplots, the expression of the CTLA-4 gene was markedly high in m6Acluster B and the expression levels of the PD-1 and PD-L1 genes were markedly low in m6Acluster C (Figures 3E–G).

PCA showed that m6A regulators could separate distinct m6A modification patterns perfectly (Figure 4A). To further investigate the potential biological behavior of each m6Acluster, a total of 2457 DEGs among three m6Aclusters were extracted eventually (Figure 4B). Similarly, the “clusterProfiler” R package was used to implement GO enrichment analysis and KEGG pathway analysis for the DEGs. The results showed that the DEGs were enriched in biological processes related to tumorigenesis and tumor progression, such as FoxO signaling pathway and ErbB signaling pathway (Figures 4C,D). The results revealed that m6A modification played a significant role in the tumorigenesis and tumor progression of PAAD.

To further investigate this regulation mechanism, univariate Cox regression analysis was performed to extract prognosis-related genes among the DEGs (Supplementary Table S2). Based on the 53 prognosis-related DEGs obtained, patients were divided into three genomic subtypes through unsupervised clustering analyses (gene clusters A-C) (Figure 5A). Survival analysis showed that patients in gene cluster B had a worse outcome, while these in gene cluster C showed a prominent survival advantage (Figure 5B). The heatmap revealed that gene clusters A-C were characterized by different signature genes (Figure 5C). Prognosis-related DEGs were overexpressed in gene cluster B, and under-expressed in gene cluster C. Moreover, gene cluster B showed higher expression levels of all m6A regulators, while gene cluster C had lower expression levels of all m6A regulators (Figures 5D; Supplementary Figure S2A). Similarly, GSVA showed different biological behaviors between three gene clusters (Supplementary Figures S2B–D). In addition, ssGSVA showed gene cluster B had a high abundance of Th2 cells and NK cells (Figure 6A). The stromal and ESTIMATE scores (p < 0.05) were higher in cluster B compared to cluster C (Figures 6B–D). The results further demonstrated that m6A modification patterns were tightly associated with the TME of PAAD. Similarly, the expression of targeted immune checkpoint molecules was different between the distinct gene clusters (Figures 6E–G).

Considering the individual heterogeneity and complexity of m6A modification in PAAD, the study used PCA to quantify the m6A modification pattern of samples based on the prognosis-related DEGs and then divided them into low- and high-m6Ascore groups. The alluvial diagram was used to visualize the attribute changes of individual samples, which showed that m6Acluster B had a high proportion of gene cluster B and was linked to a low m6Ascore (Figure 7A). Furthermore, Kruskal-Wallis test indicated a difference in m6Ascore among the m6Aclusters. The results showed that m6Acluster B had the lowest median score, while m6Acluster A had the highest median score (Figure 7B). The similar results were obtained when analyzing the correlation between m6Ascore and gene clusters. Gene cluster B had the lowest median score, while gene cluster C had the highest median score (Figure 7C).

Patients in the high m6Ascore group demonstrated a prominent survival benefit in both the TCGA cohort and GEO cohort (GSE21501) (Figures 7D,E). Moreover, there was a negative correlation between m6Ascore and survival state. The low m6Ascore group had a high proportion of patients in dead (Figures 7F,G). As shown in Figure 7H, m6Ascore was markedly related to OS (hazard ratio (HR): 0.920, 95% confidence interval (CI): 0.881–0.961, p < 0.001). In addition, a multivariate Cox regression model including age, gender, tumor grade, m6Ascore, and tumor stage confirmed that m6Ascore was an independent prognostic factor of PAAD (HR: 0.927, 95% CI: 0.885–0.970, p = 0.001) (Figure 7G). ROC analysis revealed that m6Ascore had an acceptable prognostic value for PAAD patients (1-year AUC = 0.6671, 2-year AUC = 0.6657, 3-year AUC = 0.7171, 4-year AUC = 0.7708; respectively) (Figure 8A). These results indicated that m6Ascore had a robust and stable OS-predictive ability for PAAD.

As shown in Figures 8B,C, there was a negative correlation between the m6Ascore and the TMB. Patients with high TMB had a low m6Ascore. Survival analysis revealed that patients with high TMB had a poor prognosis; and the similar result was obtained from patients who had high TMB and low m6Ascore (Figures 8D,E). Of note, there were negative correlations between m6Ascore and some infiltrated immune cells, including activated Th2 cells, regulatory T cells, and NK cells et al. In addition, the expression levels of the PD-1and PD-L1 genes were higher in patients with low m6Ascore, while there was no difference in the CTLA4 gene expression (Figures 8G–I).