Analysis of RNA m6A methylation regulators and tumour immune cell infiltration characterization in prostate cancer

Abstract Potential roles of RNA N6-methyladenosine (m6A) modification in tumour microenvironment (TME) cell infiltration has been demonstrated in recent studies. Nonetheless, the mechanism of its regulation remains unknown and immunotherapy has been marginal in prostate cancer. We demonstrated the expression of different m6A regulators within prostate cancer related to genetic variation, alternative splicing (AS), tumour mutational burden (TMB) and TME. Unsupervised clustering and risk prediction model constructed by 24 m6A regulators could predict scores of TME and prostate cancer patients prognosis. T cells CD8 was the intersection of immune cells which are related to multiple biological processes, and the fraction of T cells CD8 strongly correlates with immune associated gene sets. m6A methylation modification and immune cells infiltration played a nonnegligible role in prostate cancer. Our study represents a step towards personalized immunotherapy for prostate cancer patients.


Introduction
Prostate cancer is one of the most common malignancies worldwide and the first common cancer of the genitourinary tract in old man [1][2][3][4]. Immunotherapeutic agents offer new therapy options for patients with many cancers [5,6]. To date, responses to immune checkpoint inhibitor monotherapy have been marginal in prostate cancer and the majority of patients without clinical benefit. Tumour microenvironment (TME) plays a crucial role in the tumour progression, and numerous studies have shown that TME not only contained cancer cells but also stromal cells [7][8][9][10]. Therapeutic response of many tumours, including prostate cancer, correlates with tumour mutational burden (TMB). The higher score of TMB, the more tumour-related oncogenic mutations and the less normal-like cells, which made it easier to be detected by immune cells and became the target of tumour immunity, and more likely to be effective for immunotherapy [11,12].
N6-methyladenosine (m 6 A) is the most abundant modification among more than 100 patterns in mRNA [13,14]. m 6 A modification at different levels regulates mRNA (including structure, maturation, stability, splicing, export, translation and decay). The modification process of m 6 A is dynamic and reversible, which is mainly regulated by three kinds of enzymes: methyltransferases ("writers"), demethylases ("erasers") as well as m 6 A-binding proteins ("readers"). Under the joint efforts of these enzymes, the level of m 6 A maintains a dynamic equilibrium. The finding would help improve the abilities of predicting immunotherapeutic responsiveness, and reveal the role and mechanism of m 6 A modification [15].
Nowadays, more and more researches have revealed the special correlation between m 6 A methylation and TME immune infiltration, which cannot be explained via RNA degradation mechanism [16]. Therefore, comprehensive recognizing of the TME cell infiltration characterizations mediated by multiple m 6 A regulators will contribute to enhancing our understanding of TME immune regulation. In our study, we integrated the genomic information of prostate cancer and normal samples to comprehensively evaluate the roles of m 6 A regulators, and correlated the m 6 A regulators with the TME cell-infiltrating characteristics. We aimed to prove that m 6 A modification played a nonnegligible role in shaping individual TME characterizations.
were also downloaded from the TCGA database and ICGC database. The splicing event data for prostate cancer were downloaded from the TCGA SpliceSeq database. Regarding the copy number variations (CNVs) and single nucleotide polymorphism (SNP) analyses, we used TCGA data. R software (3.6.2) (Vienna, Austria) was used for data extraction and sorting to obtain the gene expression matrices and clinical data. Circos plots were performed by the R package "Rciorcos" in the catalogue of somatic mutations in cancer (COSMIC) database (https://cancer.sanger.ac.uk/cosmic), we get the m 6 A related genes mutation information of prostate cancer.

Prognostic analysis of alternative splicing for 24 m 6 A regulators
The splicing factor (SF) is an important regulator in the alternative splicing (AS) process. SF data of prostate cancer patients were downloaded from http://hapmap.ncbi.nlm.nih. gov/downloads/genotypes. The percent spliced in (PSI) value was used in quantifying AS events and calculating for seven types of AS events: exon skip (ES), mutually exclusive (ME) exons, retained intron (RI), alternate promoter (AP), alternate terminator (AT), alternate donor (AD) site and alternate acceptor (AA) site. The correlation between the expression of SF and the PSI value of the most significant survivalassociated AS events was analysed by Pearson's test.

Differential expression analysis and functional enrichment analysis
To identify differentially expressed genes (DEGs) between normal and tumour groups, we used the R language "limma" package to screen the DEGs from GEO and TCGA transcriptome data. Mann-Whitney's test was performed to determine differential expression levels of genes. Adjusted p value <.05 was considered to be statistically significant. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted with R language "clusterProfiler" and "enrichplot" package, false discovery rate (FDR)<0.05 was set as the threshold. Bar chart was used to visualize the biological process (BP), cellular component (CC) and molecular function (MF) of GO enrichment. The bar chart was also used to visualize the pathways of KEGG.

Estimation of TME, TMB and immune cell infiltration
We analysed 29 immune-associated gene sets which represented diverse immune cell types, functions and pathways. The immune infiltration levels were quantified using enrichment scores calculated by single-sample gene set enrichment analysis (ssGSEA) [17]. We used the ssGSEA score to classify patients into two (GEO: high, low) or three (ICGC and TCGA: high, medium, low) different groups. ESTIMATE was used to evaluate the immune cell infiltration level (immune score), tumour purity and stromal content (stromal score) for each prostate sample. CIBERSORT was used to calculate the proportions of 22 human immune cell subsets. To quantify the proportions of immune cells, the CIBERSORT algorithm was used in our study. We used the R package "CIBERSORT" to estimate the fraction of immune cells in prostate cancer samples. TMB is calculated by the total number of somatic mutations/total covered bases [18]. We constructed a model based on the TMB status and classified prostate cancer patients into high TMB and low TMB groups by the mean of population TMB. We tried to analyse whether there are differentially expressed m 6 A regulators between the high and low TMB subtypes.

Construction of risk prediction model
In order to build a prostate cancer prognostic model use m 6 A regulators, all 24 genes were used to conduct a univariate Cox survival regression, and p value <.1 was used to create a least absolute shrinkage and selection operator (LASSO) regression model. This method was applied to reduce potential over-fit and implemented through the "glmnet" R package. Univariate and multivariate Cox analyses were conducted with R language "forest" and "survival" package. The risk scores of each patient were calculated through the prediction formula of the risk prediction model. We calculated the cut-off value used to determine whether the patient is at high or low risk. The formula used for this model was: risk score ¼ P n i¼1 genei Á coefi.

Verification of the validity of the risk prediction model
Overall survival (OS) was calculated from the date of diagnosis to the date of death or last follow-up. Survival curves were estimated using the Kaplan-Meier method, and the logrank test was used to test for differences between groups. The time-dependent receiver operating characteristic (tdROC) curve analysis was first applied to evaluate the predictive accuracy of the model for cancer-specific death or biochemical recurrence based on the risk scores, with the help of the "survivalROC" R package. Chi-squared test or Fisher's exact test was used to investigate the correlation between risk model and clinicopathologic variables, and draw the heatmap through the "pheatmap" R package, p < .05 was considered statistically significant for all tests.

Identification of molecular subtypes by using consensus clustering
Unsupervised clustering was performed using the "ConsensusClusterPlus" package in R to identify subgroups based on the m 6 A regulators. This algorithm determined consensus clustering by measuring the stability of clustering results from the application of a given clustering method to random subsets of data. After executing ConsensusClusterPlus, we obtained the cluster consensus and item-consensus results. Graphical output results included heatmaps of the consensus matrices, which displayed the clustering results, consensus cumulative distribution function (CDF) plots and delta area plots, and which allowed us to determine an approximate number of clusters. Numbers of  relatively high consistency within the cluster, relatively low coefficient of variation, and no appreciable increase in the area under the CDF curve. Associations between both clinicopathologic characteristics and clustering were analysed using the v 2 test or Fisher's exact test, and draw the heatmap through the "pheatmap" R package, p < .05 was considered statistically significant for all tests.

Landscape of variation of m 6 A regulators in prostate cancer
A total of 24 m 6 A regulated genes including nine writers, two erasers and 13 readers were identified in our research. We obtained the expression levels of 24 m 6 A regulators from TCGA and GEO databases (652 prostate tumour samples and 101 normal samples). The genes that met the adjusted p value <.05 were considered DEGs. DEGs were shown in heatmap and boxplot to intuitively illustrate the differences between tumour and normal samples of TCGA ( Figure 1(A,B)) and GEO databases ( Figure S1A,B). The investigation of CNV alteration frequency showed a prevalent CNV alteration in 10 regulators and the location of CNV alteration of m 6 A regulators on chromosomes is shown in Figure 1(C) ( Table S1). Among 24 m 6 A regulators, 17 of them were experienced mutations. It was found that the ZC3H13 exhibited the highest mutation frequency, with frequency 16.67%. Further analyses revealed a significant mutation co-occurrence relationship between IGF2BP1, YTHDF2, along with RBM15B and YTHDC2 (Figure 1(D)). The COSMIC database contains 1404 prostate cancer samples and 392 metastasis samples, and the chromosomal mutation sites with differences in primary and metastasis of prostate cancer are shown in Figure  1(E). We identified differential somatic mutations between cancer and metastasis samples using the COSMIC, and found that the m 6 A regulators (KIAA1429, YTHDF2, RBM15, RBM15B, ELAVL1, IGF2BP1, LRPPRC) which experienced mutations could be the prominent factors resulting in tumour metastasis (Table S2).

m 6 A methylation modification regulatory network and enrichment analysis
The comprehensive landscape of m 6 A regulators interactions, regulator connection for prostate cancer patients was depicted with the m 6 A regulators network, and their correlations were also analysed using "corrplot" package in R software ( Figure 2(A)). In this network, the relationship between m 6 A regulators was shown in pink lines, and seven nodes have higher degrees. In TCGA and GEO databases, the top three groups with the highest correlation were KIAA1429/ YTHDF3, KIAA1429/LRPPRC and ZC3H13/FTO ( Figure S2A,B). They also have the higher correlation in ICGC database (p < .05, r > 0.4, Figure S2C). To explore the biological functions of m 6 A regulators, they were categorized into BP, CC and MF. Under stringent threshold conditions (p-adjust <.05), we identified 150 specific BP, 13 CC and 26 MF of GO terms that were enriched in these genes, and top 10 terms of each part are in Figure 2(B). Additionally, analysis using "clusterProfiler" indicated that these genes were significantly enriched in spliceosome KEGG pathway as shown in Figure 2(C).
We demonstrated that whether tumours with a high writer gene expression exhibit a low eraser gene expression actually depending on the different writer and eraser genes ( Figure S3A,B). It was found that tumours with a low expression of FTO showed a high expression of writer genes (KIAA1429, METTL14, RBM15, RBM15B, ZC3H13, RBMX, CBLL1), while it did not affect the expression of WTAP. Tumours with a high expression of eraser gene FTO showed a high expression of writer gene METTL14. In addition, it was also found that tumours with a low expression of ALKBH5 showed a high expression of writer genes (METTL14, RBM15B and ZC3H13). The above results indicated that cross-talk among the regulators of writers, readers and erasers may play critical roles in the formation of different m 6 A modification patterns.

Survival associated as events mediated by 24 m 6 A regulators
Forty-four thousand and nineteen AS events of 56,755 genes were identified in 499 prostate cancer patients, and ES was the most frequent splice signatures among the seven AS types ( Figure S4A). Notably, one gene could have multiple survival associated AS events in prostate cancer. Thus, a subset of overlapping AS events among the seven types of AS in prostate cancer were illustrated by UpSet plot diagram ( Figure 3(A)). Prognosis-related AS were shown in Volcano plot (Figure 3(B)), and the top 20 significant survival associated AS events of the seven types were presented in Figure S4B. We then screened the m 6 A genes as SFs, and the relationship between m 6 A SFs and the prognosis-related AS events is shown in Figure 3

TME cell infiltration characteristics identified by ssGSEA
In the TCGA-PRAD (prostate adenocarcinoma) cohort as well as GEO and ICGC cohorts, patients were divided into two or three immune groups based on ssGSEA method ( Figure  S5A-C). ssGSEA analysis was used to assess the immune infiltration level of every gene set. The patients in high immune infiltration cluster had a higher significant immune scores, estimate scores, stromal scores and a lower tumour purity (Figure 4(A,B), Figure S5D). The difference of TME scores (immune scores, estimate scores, stromal scores, tumour purity) between high and low immune infiltration groups was statistically significant in TCGA database (Figure 4(C)). A similar stratification of cases was confirmed in GEO and ICGC independent validation cohorts, all p < .05 ( Figure 4(D), Figure S5E). In addition, we assessed the association between immune infiltration clusters and the clinical prognosis features. This analysis clearly revealed immune infiltration clusters were related to stage N (p¼ .021) and OS. However, there was no relationship between other clinical characteristics and immune infiltration clusters (Figure 4(E)). Prognostic analysis also revealed low immune infiltration cluster to be markedly related to prolonged survival, while high immune infiltration cluster was characterized by poorer survival, all p < .05 (Figure 4(F,G)). One-way ANOVA test also confirmed the remarkable differences on m 6 A regulators expression between immune infiltration clusters in all databases ( Figure 4(H,I), Figure S5F). Venn diagram summarizes the common differentially expressed m 6 A regulators of immune infiltration clusters for TCGA, ICGC and GEO cohorts ( Figure  4(J)). Eventually, HNRNPA2B1 and METTL3 as the representative m 6 A genes played a non-negligible role in the immune regulation in TME. The expression levels of HNRNPA2B1 and METTL3 in the high immunity group were significantly lower than those in the low immunity group ( Figure S6A-C). To better illustrate the characteristics of 29 immune-associated gene sets, we also tested the correlation between the gene sets for TCGA and ICGC cohorts ( Figure 4(K,L)). The above results strongly suggested that the expression of m 6 A regulators was significantly correlated with immune infiltration levels.

Mutational genomic landscape and m 6 A genes expression in prostate cancer
The summary of mutational genomic landscape of prostate cancer is shown in Figure S7A. Missense mutation was the highest mutation in the nine variant classifications and SNP had the highest frequency in the total mutation frequency compared to INS and DEL. The results of single nucleotide   variants (SNVs) indicated C > T had the highest incidence in six base substitutions (C > A, C > G, C > T, T > A, T > C and T > G). The waterfall map summarized the high mutated genes and their mutation classification of prostate cancer ( Figure S7B). In addition, the proportion of more mutated genes was visualized using the genecloud ( Figure S7C). The co-occurrence and exclusive relationships between these mutant genes are shown in Figure S7D. The co-occurrence correlation between RYR1 and KMT2D was the most significant (p < .001). Kaplan-Meier's analysis indicated that there was significant difference of OS between the high group and low TMB group (p¼ .033, Figure S7E). Patients with less TMB had a better prognosis, suggesting that TMB as a risk-independent prognostic factor. The differences of m 6 A genes expression between the high TMB and low TMB groups were calculated and visualized by "limma" package. The result of heatmap showed that these m 6 A regulators (HNRNPA2B1, HNRNPC, KIAA1429, YTHDF2, YTHDF1, METTL3, ZC3H13, ELAVL1, RBMX, RBM15, FTO, RBM15B) presented differential expression profiles between high group and low TMB groups ( Figure S7F, Table S3).

Unsupervised clustering of 24 m 6 A regulators in the prostate cancer cohort
To further explore the characteristics of these m 6 A modification phenotypes in the different clinical traits and biological behaviours, we fixed attention on the TCGA-PRAD cohort. We selected the m 6 A regulators to identify subgroups of prostate cancer. Numbers of clusters were determined according to the following criteria: relatively high consistency within the cluster, relatively low coefficient of variation, and no appreciable increase in the area under the CDF curve. We calculated average cluster consensus and the coefficient of variation among clusters depending on category number ( Figure 5(A)). Hence, when k ¼ 4 ( Figure 5(B)) in prostate cancer could be the optimal choice with clustering increasing from k ¼ 2-9, and we named these three clusters as m 6 A gene cluster 1-4, respectively. Principal component analysis (PCA) could better distinguish cluster 4 and other clusters ( Figure 5(C)), indicating a good predictive value. Kaplan-Meier's survival analysis of Figure 5(D) revealed significant differences in prognosis among the four clusters (p¼ .012), cluster 1 had the best survival, while cluster 4 had the worst survival. However, there was no relationship between clinicopathologic characteristics and clusters ( Figure  S8A). The 29 immune-related gene signature in different clusters was presented in the heatmap with TME scores as the annotations ( Figure 5(E)). The result of Figure 5(F) revealed unsupervised clusters were related to TME, patients in cluster 4 had a highest significant immune scores, estimate scores, stromal scores and a lowest tumour purity. There were significant differences between group 1 and group 4 in comparison between immune scores, estimate scores, stromal scores and tumour purity. The differences of m 6 A genes expression between the four clusters were calculated and visualized by "limma" package. The results showed that 21 m 6 A regulators presented differential expression profiles ( Figure 5(G)). The above results strongly suggested that the unsupervised clustering of 24 m 6 A regulators was significantly correlated with TME.

Constructing and evaluating the risk prediction model by 24 m 6 A regulators
In order to investigate the prognostic value of these 24 m 6 A genes in prostate cancer, univariate cox survival regression analysis was performed based on the expression levels of the m 6 A regulators from TCGA-PRAD cohort ( Figure 6(A)). We selected the survival associated eight m 6 A regulators (p < .1) as candidates to create a lasso regression model. The lasso analysis results showed that four genes (HNRNPA2B1, ELVAL1, METTL14 and YTHDF2) were the powerful composed factors of prostate cancer risk prediction model ( Figure S8B). All samples in TCGA-PRAD cohort were divided into two groups according to the median value of risk score in our model. In the m 6 A risk prediction model, the risk scores and the number of deaths were significantly higher in the high risk group compared to the low risk group, and the heatmap of composed genes expression is shown in Figure 6(B). The expression levels of the m 6 A regulators in high-risk and lowrisk group were presented in the heatmap with clinicopathologic variables as the annotations (Figure 6(C)). In this risk prediction model, the results showed that there were significant differences between the high-risk and low-risk groups in term of stage T (p < .05), and stage N (p < .01). Univariate and multivariate Cox regression analyses were performed to test whether the risk model was an independent prognostic factor. As a result, only the risk score was significantly related to OS of prostate cancer patients in both univariate (p¼ .022, Figure S8C) and multivariate Cox regression analyses (p ¼ .044, Figure S7D). The ROC curve was used to assess the sensitivity and specificity of the risk model and the result showed that AUC values were 0.78, suggesting well-prediction performances (Figure 6(D)). Kaplan-Meier's survival analysis revealed that the high-risk group had significantly shorter survival time compared to low-risk group (Figure 6(E)). Among the four m 6 A regulators that was built of the risk prediction model, the expression levels of risk score were related to stage T (p < .001) and stage N (p¼ .002, Figure 6(F)), HNRNPA2B1 expression levels were related to stage T (p < .001) and stage N (p < .001, Figure S9A), ELVAL1 expression levels were related to stage T (p < .001) and stage N (p¼ .011, Figure S9B), YTHDF2 expression levels were only related to stage N (p¼ .011, Figure S9C).
The 29 immune-related gene signature in risk prediction model was presented in the heatmap with TME scores as the annotations (Figure 6(G)). The results of Figure 6(H) revealed  ARTIFICIAL CELLS, NANOMEDICINE, AND BIOTECHNOLOGY risk prediction model was related to TME, patients in highrisk group had a higher immune scores (p < .05), estimate scores (p < .01), stromal scores (p < .001) and a lower tumour purity (p < .01) compared to low-risk group. There were significant differences between high-and low-risk group in comparison between immune scores, estimate scores, stromal scores and tumour purity. The above results strongly suggested that the risk prediction model made by 24 m 6 A regulators was significantly correlated with TME.

Immune cell infiltration analysis in multiple types
In our study, we investigated the difference of immune infiltration in prostate cancer and normal tissues in 22 subpopulations of immune cells by CIBERSORT algorithm, and the result is shown by violin plot (Figure 7(A)). The differences of 22 immune cell phenotypes in three immune infiltration clusters (high, medium, low) were calculated and visualized by "limma" package. The result of boxplot showed that 16 immune cells presented differential expression profiles between three clusters (Figure 7(B)). We then used the CIBERSORT method to compare the component differences of immune cells among the two TMB groups (Figure 7(C)). T cells CD8, NK cells activated, monocytes, macrophages M2 were considered to be related to the TMB grouping. To better depict the function of immune cell phenotypes, we examined the immune infiltration signatures in unsupervised clustering and lasso risk prediction model of m 6 A regulators with prostate cancer patients. The result of boxplot showed that six immune cells (T cells CD8, T cells CD4 memory resting, T cells gamma delta, NK cells resting, macrophages M0, mast cells resting) presented differential expression profiles between four unsupervised clusters (Figure 7(D)), and four immune cells (T cells CD8, T cells CD4 naive, T cells gamma delta, dendritic cells activated) were considered to be related to risk prediction grouping (Figure 7(E)). We were interested in common immune cells which are related to multiple BPs, so we focussed on the T cells CD8 in the result of the Venn diagram analysis (Figure 7(F)).

Difference in the immune-associated gene sets among two T cell CD8 groups
To better illustrate the fractions of 29 immune-associated gene sets in low and high T cell CD8 groups, we tested the differences between the two T cells CD8 groups. The results showed that there was a significant difference in the level of T cells CD8 in five immune-associated gene sets (iDCs, inflammation-promoting, Th2 cells, B cells, T cell co-stimulation) ( Figure S10A). We then studied the genes expression in above five immune-associated gene sets, and screening the DEGs between high and low T cells CD8 groups. The results confirmed that T cells CD8 was significantly related to the genes expression in inflammation-promoting, Th2 cells, B cells and T cell co-stimulation sets, and not related to the genes expression in iDCs set ( Figure S10B-F). Discussion m 6 A plays a key role in the fate of mRNAs, including their nuclear export, translation, decay and AS [19,20]. Recent  evidence suggests that m 6 A modification is closely connected to prostate cancer. VIRMA and YTHDF3 are increased in stages III/IV than in stage II in prostate cancer. Li et al. [21] found that YTHDF2 is highly expressed in prostate cancer and the downregulation of YTHDF2 considerably increases the level of m 6 A and inhibits the proliferation and migration of prostate cancer cell lines. Cai et al. [22] reported that the silencing of METTL3 by shRNA inhibits cell proliferation, colony formation and invasion of cells and inhibits the growth of tumours in vivo. However, little is known about the role of different m 6 A regulators in the development and progression of prostate cancer. To the best of our knowledge, this study constitutes the first systemic analysis the genetic variations (including SNP, CNV and somatic mutations) of m 6 A regulators in prostate cancer. We found that different genetic variations could be the prominent factors related to the different m 6 A regulators expression. The investigation of genetic alteration frequency showed a prevalent CNV alteration in 10 m 6 A regulators, and m 6 A regulators which experienced somatic mutations could be the prominent factors resulting in prostate cancer metastasis.
Several studies in the exploration of AS signature have previously revealed that the aberrant expression of m 6 A is also associated with AS [23,24]. We identified m 6 A genes and AS events which are related to OS in prostate cancer through the analysis of TCGA program to gain comprehensive insight into differential m 6 A regulators splicing patterns. Observations in this study provided solid evidence that one m 6 A gene could yield several mRNAs via AS, resulting in multiple transcripts and various proteins, some exert antagonistic functions and some exert coordinating functions. However, more m 6 A genes play antagonistic role in survival in our research. This finding did not align with the previous study that reported most survival associated AS events in ovarian cancer were favourable prognostic factors [25]. Conversely, prostate cancer exhibits higher levels of adverse prognostic factors, but lower levels of favourable prognostic factors. Therefore, AS primarily regulates disease progression through the antagonizing functions of splice variants. One variant from AS might affect the regular function of another variant.
ssGSEA was implemented to calculate the immune cell infiltration levels for each sample in order to determine the immune infiltration groups. Prognostic analysis revealed low immune infiltration cluster to be markedly related to prolonged survival, while high immune infiltration cluster was characterized by poorer survival. The patients in high immune infiltration group had a higher significant immune scores, estimate scores, stromal scores and a lower tumour purity. The lower the tumour purity, the higher the degree of heterogeneity and malignancy. This result may be the main reason for the poor survival to the patients in high immune infiltration group. As most studies focus on single TME cell type or single regulator, the overall TME infiltration characterizations mediated by integrated roles of multiple m 6 A regulators are not comprehensively recognized. Identifying the role of distinct m 6 A modification patterns in the immune infiltration will contribute to enhancing our understanding of TME antitumor immune response, and guiding more effective immunotherapy strategies. From the three databases, we analysed the differential m 6 A regulators in the high and low groups. At the same time, HNRNPA2B1 and METTL3 of the above three clusters were identified as the immune strongly related genes associated with prostate cancer. These results may illustrate that HNRNPA2B1 and METTL3 expression may affect the immune microenvironment of prostate cancer to a certain extent, this is consistent with the conclusion of the recent study [26,27], and thus supports the reliability of our conclusions to a large extent.
TMB is associated with the abundance of neoantigens and increased immunogenicity. In this study, T cells CD8 and macrophages M2 cell infiltration in the high TMB group was significantly higher than that in the low TMB group for prostate cancer patients. However, it should be noted that the correlation between the infiltration of different immune cells and the level of TMB was not consistent in our research. NK cells activated and monocytes infiltration in the high TMB group was significantly lower than that in the low TMB group for prostate cancer. Patients in the low TMB group had a better prognosis in our study. However, higher TMB was associated with better prognosis for the patients likely with immunotherapy such as melanoma. The reason behind the discrepancy to our results could be that most of the TCGA patients were not treated with immunotherapy [28]. Further, the differences between expression of m 6 A regulators have been proved to be significantly associated with TMB and immune related biological pathways. These DEGs were considered as m 6 A-related signature genes. The above results strongly suggested that m 6 A regulators was significantly correlated with immune infiltration, TMB, SNP, CNV and AS. Our study not only focuses on immune-related m 6 A genes, and uses strict standard layer screening to obtain genes which may be potential prognosis targets of prostate cancer.  (G) The heatmap shows the unsupervised clustering of the 29 immune cell types from TCGA cohort among four m 6 A gene clusters. The distribution of tumour microenvironment (TME) was compared between the high-risk and low-risk groups. (H) The TME scores (immune scores, estimate scores, stromal scores, tumour purity) were compared between the high-risk and low-risk groups. The asterisks represented the statistical p value ( Ã p < .05, ÃÃ p < .01 and ÃÃÃ p < .001). Further, in this study, four m 6 A gene clusters were identified based on m 6 A regulators, which were significantly correlated with survival. The mRNA transcriptome differences between distinct m 6 A gene clusters have been proved to be also significantly associated with m 6 A and TME. This demonstrated again that the m 6 A clusters were of great significance in shaping different TME landscapes. Therefore, a comprehensive assessment of the m 6 A clusters will enhance our understanding of TME cell-infiltrating characterization. Murphy et al. [29] used multiple omics platforms to integrate biomarkers to improve the stratification of patients with aggressive and indolent prostate cancer. Univariate Cox regression followed by lasso regression was used to validate m 6 A regulators and construct a risk prediction model for prostate cancer patients in our study. The low-risk group had significantly better prognosis than the high-risk group in risk model. The results of ROC also showed that the model we constructed exhibited good classifier performance. Although the well predictive value of model was available, the significance of these m 6 A regulators in relation to tumour classification, survival time and prognosis need to be confirmed in more groups of patients. These findings emphasize the important role of m 6 A regulators in constructing tumour risk models and evaluating prognosis.
Therapeutic response of many tumours, correlates with amounts of tumour-infiltrating lymphocytes (TILs), tertiary lymphoid structures (TLSs), the number of neoantigens, and TMB [30]. The distribution of T cells within the TME can be suggestive of an efficacious anti-tumour T cell response, with several groups showing enrichment for T cells at the periphery or centre related to improved outcome [31]. Simoni et al. [32] found that human lung and colorectal cancer CD8þ TILs can not only be specific for tumour antigens, but also found that these bystander CD8þ TILs have diverse phenotypes that overlap with tumour-specific cells. In this study, the differences of T cells CD8 fraction were related to immune infiltration, m 6 A gene clusters, TMB and risk scores in prostate cancer. Our data also revealed the importance of T cells CD8 in the research of prostate cancer, which is a new discovery and also consistent with the results of other tumour before. However, a large number of experiments are still needed for further study.
Our study had some limitations. First, although we found that m 6 A regulators play different roles in the multiple BP, such as mutational landscape, epigenome, TMB and TME. However, the potential molecular mechanisms were not evaluated, and the detailed mechanism of m 6 A regulators in prostate cancer warrants further investigation. Second, m 6 A gene clusters and risk model might improve prognosis of prostate cancer patients at a higher risk group for immunotherapy. Hence, further research is necessary to explore whether m 6 A regulators could be used as diagnostic markers or therapeutic targets in prostate cancer and guiding more effective immunotherapy strategies.

Conclusions
In summary, our data provide strong evidence that m 6 A methylation modification in the mutational landscape, epigenome, TMB and TME. m 6 A gene clusters and risk model  might improve prognosis of prostate cancer patients at a higher risk group for immunotherapy. The difference of m 6 A regulators plays different roles in multiple BP. The comprehensive evaluation of individual tumour m 6 A modification pattern will contribute to enhancing our understanding of tumour characterization and guiding more effective immunotherapy strategies.

Data availability statement
The dataset supporting the conclusions of this article is included within the article.