m6A Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Pancreatic Cancer

Recent studies have shown that RNA N6-methyladenosine (m6A) modification plays an important part in tumorigenesis and immune-related biological processes. However, the comprehensive landscape of immune cell infiltration characteristics in the tumor microenvironment (TME) mediated by m6A methylation modification in pancreatic cancer has not yet been elucidated. Based on consensus clustering algorithm, we identified two m6A modification subtypes and then determined two m6A-related gene subtypes among 434 pancreatic cancer samples. The TME characteristics of the identified gene subtypes were highly consistent with the immune-hot phenotype and the immune-cold phenotype respectively. According to the m6A score extracted from the m6A-related signature genes, patients can be divided into high and low m6A score groups. The low score group displayed a better prognosis and relatively strong immune infiltration. Further analysis showed that low m6A score correlated with lower tumor mutation burden and PD-L1 expression, and indicated a better response to immunotherapy. In general, m6A methylation modification is closely related to the diversity and complexity of immune infiltration in TME. Evaluating the m6A modification pattern and immune infiltration characteristics of individual tumors can help deepen our understanding of the tumor microenvironment landscape and promote a more effective clinical practice of immunotherapy.


INTRODUCTION
More than 160 RNA modifications including N 7 -methylguanine (m 7 G), N 6 -methyladenosine (m 6 A), and N 5 -methylcytosine (m 5 C) have been identified. These modifications play a significant role in regulating RNA fate (1). In eukaryotes, m 6 A is regarded as the most important and abundant mRNA modification, accounting for more than 80% of all RNA methylation modifications (2). It is now clear that m 6 A methylation exists in almost all types of RNA, including coding RNA and non-coding RNA (3). The m 6 A modification is catalyzed by RNA methyltransferases such as METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B (writers), while the modification is removed by demethylases such as FTO and ALKBH5 (erasers). In addition, modifications can be recognized by m 6 A binding proteins, such as YTHDC1/2, YTHDF1/2/3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1/2/3, and RBMX (readers) (4)(5)(6). A growing body of evidence shows that m 6 A regulators are involved in vital biological processes and are dynamically regulated in many physiological and pathological processes (7)(8)(9). Abnormal expression and genetic alterations of m 6 A regulators are closely related to events such as developmental defects, metabolic disorders, abnormal immune regulation, and tumor progression (10,11). A comprehensive understanding of potential m 6 A regulators' expression perturbation and genetic variation under cancer heterogeneity will facilitate the identification of therapeutic targets based on RNA methylation (12,13).
Pancreatic cancer (PC) is a highly lethal malignant tumor, which is also one of the leading causes of cancer death worldwide (14). In the past 25 years, its global burden has more than doubled (15). With the in-depth comprehending of pathology, the diversity and complexity of tumor microenvironment have been increasingly understood, and immune cell subgroups which intimately involved in tumor genesis, metastasis and treatment are gradually recognized (16)(17)(18)(19). Tumor microenvironment (TME) has been found to play an important role in ineffective treatment and poor prognosis of pancreatic cancer. Many molecules and related signal transduction pathways in the microenvironment can promote cancer metastasis or immunosuppression. A variety of soluble immunosuppressive molecules and immunosuppressive cells can lead to the disorder of immune effector cells, thus forming a unique immunosuppressive environment for pancreatic cancer (20).
Increasing research focuses on whether the components of TME (including acellular matrix, pancreatic stellate cells, immune cells and soluble factors) can be used as effective targets for PC therapy (21). Recently, immune checkpoint inhibitors and chimeric antigen receptor T (CAR-T) cell therapy have become popular immunotherapies related to TME in pancreatic cancer (22,23). However, whether these therapies can bring clinical benefits remains to be fully studied. A comprehensive understanding of the characteristics of the tumor microenvironment in PC will make a beneficial contribution to the research of immunotherapy and provide new insights for basic and clinical applications.
An increasing number of studies have confirmed the close correlation between TME infiltrating immune cells and m 6 A modification, which could not be completely explained by the mechanism of RNA degradation. According to the study of Liu et al. (24), FTO enhances protein expression by regulating the m 6 A modification of JUNB and CEBPB genes, thereby promoting tumor glycolysis and inhibiting T cell effects. The FTO inhibitor Dac51 can inhibit FTO-mediated demethylation, inhibit the glycolytic ability of tumor cells, increase T cell infiltration, and have a synergistic effect with anti-PD-L1 therapy. The study of Han et al. (25) showed that YTHDF1 recognizes and binds to the transcript encoding lysosomal protein modified by m 6 A, increases the translation of lysosomal cathepsin in dendritic cells (DC), while inhibition of cathepsin can significantly increase the ability of crosspresenting antigen of dendritic cells. The absence of YTHDF1 in DC can enhance the cross-presentation of tumor antigens and the cross-priming of CD8+ T cells, thereby increasing the antitumor response of CD8+ T cells and enhancing the therapeutic effect of PD-L1 checkpoint blockade. It has been reported that cytotoxic tumor-infiltrating CD8+ T cell is increased in METTL3 or METTL14 deficient tumor. Depletion of METTL3 and METTL14 can inhibit m 6 A modification and enhance the response to anti-PD-1 therapy in colorectal cancer and melanoma (26). However, the above studies have focused on one or two m 6 A regulators and immune cell types, while tumor formation and suppression are the results of the highly synergistic effects of multiple regulatory factors. Therefore, the comprehensive analysis of the infiltration characteristics of tumor microenvironment mediated by m 6 A regulator is helpful to promote the cognition of tumor immune regulation.
At present, there are widely accepted molecular subtypes in pancreatic cancer, such as two tumor-specific subtypes and stromal subtypes identified by Moffitt et al. and purity Independent Subtyping of Tumors (PurIST) developed by Rashid et al. (27,28). These classifications may mainly be concerned with the components of the tumor at the pathological level and show good clinical value. We aimed to identify new subtypes from the m 6 A methylation modification direction and construct scores to supplement the existing clinical variable information, and correlate these analyses with the tumor microenvironment. Meng et al. have proposed an m 6 A-related mRNA signature, which can play a good prognostic predictive effect in pancreatic cancer (29). However, their study divided groups based on the existence of alterations (mutation and/or CNV) of m 6 A-related genes, and then identified differentially expressed genes for model construction. According to previous studies (30), data including gene expression profile, somatic mutation, and DNA methylation information can be used to identify the primary sites and origins of tumors, but the accuracy of gene expression data is the highest, especially in pancreatic cancer. Since gene expression data may provide more clinical value, this study identified subtypes based on m 6 A regulator gene expression to mine more accurate clinical subtypes and prognostic indicators.
In this study, we integrated the transcriptome information of 434 pancreatic cancer samples from five independent cohorts, comprehensively evaluated m 6 A modification patterns, and correlated the characteristics of immune cell infiltration in TME. Through the unsupervised clustering method, we identified two different m 6 A modification subtypes and defined two m 6 A-related gene subtypes. We found that distinct subgroups were accompanied with different immune cell infiltration characteristics. In addition, we constructed a scoring scheme to quantify individual m 6 A modification patterns, and predicted the prognosis and response to immunosuppressive therapy based on the score. Our findings indicate that m 6 A modification is closely related to TME immune cell infiltration, and can be used as a favorable predictor of prognosis and immunotherapy in pancreatic cancer.

Collection and Preprocessing of PC Public Datasets
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/).

Unsupervised Clustering of 23 m 6 A Regulators
We searched the relevant literature on m 6 A methylation modification to identify recognized m 6 A regulators for subsequent analysis. A total of 23 m 6 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 6 A methylation modification subtypes according to the expression of 23 m 6 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 (D(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).

Estimation of Immune Infiltrating Cells in TME
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.

Identification of Differentially Expressed Genes Between Different m 6 A Modified Phenotypes
Previous consensus clustering algorithms divided patients into two different m 6 A modification subtypes based on the expression of 23 m 6 A regulators. R package Limma was used to identify DEGs between the two m 6 A modification clusters with adjusted P value < 0.05.

Functional and Pathway Enrichment Analyses of DEGs
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.

Construction of m 6 A Score
To quantify the modification pattern of m 6 A in individual PC patients, we used principal component analysis (PCA) to construct the m 6 A scoring scheme. Firstly, univariate Cox regression model was performed on DEGs identified between different m 6 A modification clusters, and the genes with significant prognosis effect were selected for clustering samples and constructing m 6 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 6 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 6 A score: m6Ascore =

Verification of the m 6 A Score
To verify the reliability and clinical application value of m 6 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 6 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 6 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 6 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.

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

Obtain the Prediction Indicators of Immune Response
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.

Statistical 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 6 A regulators and prognosis. According to the correlation between m 6 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 6 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.

Construction of Genetic Variation, Immune Infiltration, and Prognostic Landscape of m 6 A Regulators
In this study, 23 m 6 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 6 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 6 A regulators. There were differences in the expression of multiple m 6 A regulators between the wild group and mutant group ( Supplementary  Figures 3, 4). Analysis of copy number variation of 23 m 6 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 6 A regulators were correlated with the prognosis of PC patients (Supplementary Figure 2). Univariate Cox regression model revealed the prognostic value of 23 m 6 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 6 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 6 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).  showed the expression of 23 m 6 A regulators in the two modified subtypes by heatmap ( Figure 2E). The expression levels of 23 m 6 A regulators between the two m 6 A clusters were also compared and shown in Figure 2F. To explore the internal biological changes under different m 6 A modification modes, we compared the composition of immune cells in TME. The result showed that m 6 A cluster A was characterized by higher infiltration of memory B cells and activated memory CD4+ T cells. In m 6 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 6 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-b signaling pathway. These results suggested that the m 6 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.

Identification of m 6 A-Related Gene Subtypes
Although the consensus clustering algorithm based on 23 m 6 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 6 A clusters, and obtained 719 survival-related genes which were named m 6 A-related signature genes (Supplementary Table 5). Based on representative m 6 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 6 A-related signature genes in two gene clusters ( Figure 3E). The expression levels of 23 m 6 A regulators between the two m 6 A-related gene clusters were also compared. Significant differences in the expression of m 6 A regulators were observed, which was consistent with the expected result of m 6 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.

Construction of m 6 A Score
Although the above results demonstrated the role of m 6 A methylation modification in the regulation of immune cell infiltration and prognosis, these analyses cannot accurately predict the m 6 A methylation modification pattern in a single tumor patient. To obtain a quantitative index of the m 6 A modification landscape of PC patients, we extracted the scores of principal component 1 and principal component 2 for calculating the final m 6 A score. Figure 4A showed that patients' m 6 A score in m 6 A cluster A was lower than m 6 A cluster B, and Figure 4B depicted that the m 6 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 6 A score construction ( Figure 4C). Subsequent analysis revealed the prognostic significance of m 6 A score. According to Figure 4D, the patients' survival rate (41% vs. 23%) in the m 6 A low score group was much higher than that in the high score group. The overall average m 6 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 6 A score group was significantly better ( Figure 4F, P < 0.001).

Validation of m 6 A Score and Its Application in Clinical Evaluation
To verify the m 6 A score, the 1-, 3-, and 5-year ROC curves were drawn, and the value of the area under curve (AUC) of the m 6 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 6 A score had the most considerable AUC value ( Figure 5C) 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 6 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.

Correlation Between m 6 A Score and Somatic Variation
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 6 A score to clarify the genetic imprint of each m 6 A score subgroup. Correlation analysis showed a significant and positive association between m 6 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 6 A score, and the low m 6 A score group always showed a survival advantage ( Figure 6C). In addition, we analyzed the somatic mutation landscape in the high and low m 6 A score groups, and found that the high m 6 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 6 A score group, which may be related to the poor prognosis of high m 6 A score group. These data can more comprehensively describe the impact of m 6 A score on genomic variation, and may provide new ideas for studying the potential interaction between m 6 A methylation modification and somatic mutation.

The Role of m 6 A Score in Predicting the Effect of Immunotherapy
The treatment of immune checkpoint inhibitors represented by CTLA-4/PD-1 inhibitors is undoubtedly a major progress in antitumor therapy. We compared the expression of common immune checkpoint genes between the high and low m 6 A score groups, and found that PD-L1 was highly expressed in the high m 6 A score group, PDCD1 was highly expressed in the low m 6 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 6 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 6 A methylation modification in pancreatic cancer may play an important role in mediating the immune response.

DISCUSSION
Recent studies have shown that m 6 A methylation modification plays an indispensable role in a variety of immune-related biological processes, including innate and acquired immune response, immune recognition, immune cell dynamic balance, and anti-tumor immune response (22). Since most studies mainly focus on the regulatory relationship between single m 6 A regulator and immune cell type, the comprehensive landscape of TME immune-infiltrating mediated by multiple m 6 A regulators in PC has not been fully understood. Therefore, clarifying the characteristics of immune cell infiltration in different m 6 A modification patterns will help us to improve our understanding of the anti-tumor immune response in TME, and provide new insights for the risk stratification of patients and the choice of clinical treatment strategies. Based on 23 m 6 A regulators, we first identified two m 6 A modification subtypes. Differences in mRNA transcriptomes between different m 6 A modification subtype were found to be closely related to tumor metabolism and immune-related biological processes, and differentially expressed genes (DEGs) were significantly enriched in glucose metabolism-related pathways, chemokine-related pathways, cytokine-related pathways, and TGF-b signaling pathways. DEGs correlated to the prognosis of PC were defined as m 6 A-related signature genes. Based on the m 6 A signature genes, we determined two m 6 Arelated gene subtypes. In the two gene clusters, we found that gene cluster A had lower immune score, stromal score and immune response-related T cell infiltration, which suggested an immune cold phenotype. In contrast, gene cluster B showed a relatively high immune score and T cell infiltration, which corresponded to the immune activation phenotype, namely hot  tumor. Survival analysis showed that gene cluster A characterized by immune cold phenotype was associated with poor prognosis, while gene cluster B characterized by anti-tumor immune response was associated with good prognosis. We speculate that patients in gene cluster B may benefit from immunotherapy. Our results are consistent with those of previous TME studies, which also indicates that m 6 A methylation modification is of great significance for shaping different TME immune characteristics.
In view of the individual heterogeneity of m 6 A methylation modification, it is necessary to quantify the m 6 A modification pattern of a single tumor sample. Scoring models based on specific biomarkers between m 6 A modified subtypes have been well established in gastric cancer and colorectal cancer to improve the choice of clinical treatment and prognosis of patients (12,13). Based on principal component analysis and a method similar to GGI, we established an m 6 A scoring scheme for PC patients. Gene cluster B with immune hot phenotype showed lower m 6 A score, while gene cluster A with immune cold phenotype indicated higher m 6 A score. Kaplan-Meier analysis showed that the m 6 A score had good prognostic predictive ability. The survival rate of patients in the low m 6 A score group was higher and the prognosis was better. These results suggest that the m 6 A score is a reliable index to comprehensively evaluate the m 6 A modification pattern of individual tumors, and can be used to further determine the characteristics of TME immune cell infiltration, namely tumor immunophenotype. Besides, verification in the TCGA-PAAD cohort showed that m 6 A score can be used as an independent prognostic indicator for PC patients. And the nomogram constructed by m 6 A score combined with other clinical variables can effectively predict the prognosis of patients. Evaluation of potential mutation driver genes in tumors is an important means to explore the potential mechanisms of cancer occurrence and development, which is conducive to cancer diagnosis and rational selection of treatment strategies. We found that the mutation rates of KRAS, TP53, CDKN2A, and AMAD4 were significantly increased in the high m 6 A score group. An important feature of KRAS mutant tumors is the immunosuppressive state (40). KRAS signaling induces the expression of immune regulatory factors and inflammatory cytokines in tumor cells, and subsequently recruits neutrophils and myeloid-derived suppressor cells (MDSC) to form an immunosuppressive tumor microenvironment (41). Mutated KRAS in pancreatic cancer plays a central role in tumor development and growth by regulating T-cell cytokines in TME. By acting on downstream effectors, KRAS leads to impaired T-cell recognition of tumor cells, which may mediate immunity escape (22,42). There are mutations of TP53 in most types of cancers. The deletion or mutation of TP53 in cancers will affect the recruitment and activity of T cells, leading to immune evasion and promoting cancer progression (43,44). The loss of P53 (encoded by TP53) in pancreatic cancer leads to increased infiltration of regulatory T cells (Tregs) in the peripheral and intratumoral tissues (45). CDKN2A is a multifunctional gene that prevents the cell cycle at the G1/S  checkpoint through the CKD4/6 regulatory mechanism. It is reported that approximately 60% of patients with pancreatic ductal adenocarcinoma carry the CDKN2A mutation, and this mutation is associated with a high risk of tumor development (46,47). SMAD4 is a member of the SMAD family and participates in the transforming growth factor-b (TGF-b) pathway, which inhibits the activity of normal immune cells and promotes the immune escape of cancer cells (48,49). These studies suggest that KRAS, TP53, CDKN2A, and AMAD4 mutations may be involved in the formation of immune suppression and immune escape in the high m 6 A score group. These m 6 A score-related gene mutations are closely related to the immune activity in TME, indicating that there may be a potential interaction between m 6 A methylation modification and tumor immune genomics. Because the mutation data of pancreatic cancer in the TCGA database is not sufficient, and only a few genes have obvious somatic mutations, it is necessary to verify the mutation oncoprint and explore the underlying mechanism in a larger data set.
In this study, we demonstrated that m 6 A modification patterns played an important role in the formation of different TME immune infiltration landscapes, which suggested that m 6 A methylation modification may affect the therapeutic effect of immune checkpoint inhibitors. We found that the expression of PD-L1 was higher in the high m 6 A score group. Previous studies suggested that pancreatic cancer had an immunosuppressive tumor microenvironment with high PD-L1 expression, which inhibited the cytotoxicity of activated T cells, and PD-L1 overexpression was associated with a poor prognosis (50,51). We also compared the IPS that predicted the efficacy of anti-PD-1/ CTLA-4 regimens in the high and low m 6 A score groups. The low score group had higher IPS, which indicated a relatively better immunotherapy effect. However, our results do not imply causal associations of m6A score and anti-tumor immunity in PC, more clinical evidence needs to be collected in future studies to verify the relationship between m 6 A score and immunotherapy. The above analysis suggests that m 6 A modification characteristics combined with TME status, tumor mutation burden, neoantigen load, PD-L1 expression, IPS and other biomarkers may be a more effective predictive strategy for immunotherapy.
Our research still has some shortcomings. Although we included 23 recognized m 6 A regulators through literature  review, newly identified regulators still need to be added into the model to improve the accuracy of the identification of m 6 A methylation modification patterns. In addition, since not all patients with low m 6 A scores can benefit from immunotherapy, more clinicopathological features need to be combined to improve the accuracy of prediction. Although we obtained 434 PC samples from different cohorts, the number of samples may be relatively insufficient, and our findings need to be further validated in a prospective cohort of PC patients receiving immunotherapy.
The study was done within tumor microenvironment in a whole, without distinguishing tumor component, immune component, and stromal component furthermore. This may cause some subtype information to be masked due to the mixture of the component, which is also a shortcoming of our research. We were more concerned about proposing molecular subtypes related to m 6 A methylation in the overall tumor microenvironment and further constructing scores. Subsequent clinical analysis showed that the m 6 A score could be used as an important supplement to existing clinical variables, and could effectively predict the prognosis of patients in combination with other clinical indicators. We may refine the differentiation of the various components of the tumor microenvironment in subsequent research work, and try to use single-cell analysis to distinguish cell types to obtain more information.
This study has provided some new insights into the clinical application of immunotherapy. By targeting m 6 A regulators or m 6 A-related signature genes to change the m 6 A modification pattern and further reverse the poor infiltration of immune cells in TME, that is, the transformation of immune cold tumors to hot tumors, may contribute to the future development of new immunotherapy drugs or combination therapy strategies. In addition, the combination of therapeutic strategies for KRAS, TP53, CDKN2A, and AMAD4 mutations with immunotherapy may open up a new way for the selection of treatment options and reverse the immunosuppressive state in tumors. These findings are conducive to the identification of different immunophenotypes, thereby improving the patient's response to immunotherapy, and can promote the clinical practice of personalized immunotherapy for cancer.
In conclusion, we evaluated 23 regulators-mediated m 6 A methylation modification landscapes based on 434 PC samples, and correlated m 6 A modification with the TME immune-infiltrating characteristics. And we constructed the m 6 A score, which can comprehensively evaluate the m 6 A modification pattern and immune-infiltrating characteristics of individual tumors, and further determine the tumor immunophenotype to guide clinical application.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.