Introduction

Gastric cancer (GC) is considered a common cancer type [1]. Genetic alterations include gene mutations and SCNAs that inactivate tumor suppressors and activate oncogenes related to WNT, EMT, PI3K/AKT/mTOR, RTK/PIK, and RTK/RAS/MAPK pathways [2]. Aberrant promoter DNA methylation and histone modification are major epigenetic alterations in GC [2]. RNA modifications regulate gene expression and are considered a promising therapeutic target for GC [3].

RNA modifications are identified on mRNA, miRNA, and lncRNAs [3]. m6A is the most important RNA alteration in eukaryotic cells. “Writers”, “erasers”, and “readers” of m6A add, remove, or recognize m6A-modified sites and change biological processes. m6A is related to various aspects of RNA metabolism and plays a key role in GC development [4]. m6A regulator expression is linked to the diversity of TME [5]. The important m6A “writer” METTL3 facilitates the Epithelial-mesenchymal transition (EMT) program through METTL3/ZMYM1/E-cadherin signaling and promotes angiogenesis and glycolysis in GC through the enhancement of HDGF mRNA stability [6, 7]. The known m1A modification “writers” include TRMT10C, TRMT6, TRMT61A, and TRMT61B [3, 8]. m1A could regulate the structure and stability of tRNA and rRNA. Recent studies have also revealed m1A in eukaryote mRNAs, promoting translation [9, 10]. The m1A “writers” TRMT6 and TRMT61A increase m1A methylation in a subset of tRNAs and increase PPARδ translation for the activation of Hedgehog signaling and driving of self-renewal of liver CSCs and tumorigenesis [11]. The TRMT6/61A complex is overexpressed in bladder cancer, causing increased m1A modification on tRFs and dysregulation of the tRF targetome [12]. APA generates transcript isoforms in coding regions or 3’UTRs [13]. CPSF, CSTF, CFI, PCF11, CLP1, and NUDT21 form the core pre-mRNA 3’end processing complex. PABPN1 controls RNA transcripts’ poly(A) tail length through binding at proximal poly(A) sites [13]. APA events are related to a range of cancers, including GC. A deep sequencing-based approach revealed APA-mediated 3’UTR shortening of 513 genes across the GC genome. The 3’UTR shortening of the oncogene NET1 enhances the activity of transcription and increases GC cell migration and invasion in vitro [14]. APA factor NUDT21 is upregulated in GC and promotes tumor growth and metastasis through the upregulation of SGPP2 [15]. Editing of A-to-I RNA is mediated by ADAR enzymes. It is observed that ADAR can convert adenosines to inosines in double-stranded RNA substrates, which leads to codon changes and the diversification of protein functions [16]. The A→I deamination reaction is catalyzed by ADARs, including ADAR, ADARB1, and ADARB2 [16]. Primary GC displays an RNA mis-editing phenotype with ADAR/ADARB1 dysregulation from the genomic gain and loss of the ADAR and ADARB1 genes. This ADAR/ADARB1 imbalance is linked to poor prognosis in GC patients [17].

Emerging evidence highlights the crosstalk between numerous types of RNA modifications. FTO can bind to pre-mRNAs in intronic regions in the proximity of alternatively spliced (AS) exons and poly(A) sites [18]. m6A “writers” TRMT6/61A catalyse the formation of m1A and m6A in tRNAs [19]. Depleting m6A “writers” KIAA1429 and METTL3 induces 3’UTR lengthening several hundreds of mRNAs with over 50% targets [20]. m6A RNA modifications and A-to-I RNA editing negatively correlate [21]. Previous studies have been confined to single RNA modifications in cancer studies.

Immunotherapy using ICIs has rejuvenated tumor immunology. However, tumor efficacy is variable and limited to a subset of cancer patients [22]. Immune infiltrates in the TME play a central role in tumor development [22]. ImmAPA score pipelines have identified tumor-specific ImmAPAs [23]. Investigation of the integrated cross-talk of multiple RNA modifications during TME cell infiltration and immunotherapy is required to further improve current immunotherapeutic strategies. In this study, expression profiles and clinical information of 1051 GC patients were integrated to identify RNA modification patterns.

Results

Genetic and epigenetic variations of RNA modification “writers” in GC

Twenty-six RNA modification “writers” were investigated. The mutation frequency of the 26 “writers” were first analyzed across pan-cancer in TCGA (Supplementary Figure 2A). Amongst the 433 STAD samples, 113 (26.1%) had mutations in one or more “writers” (Figure 1A). We found ZC3H13 showing the highest mutation rate (8%). Missense mutations and frame-shift deletions were the predominant mutational types. The co-occurrence of mutations was identified in various pairs of genes, including CSTF3 and RBM15, ADARB1 and ZC3H13, and CFI and PCF11 (Supplementary Figure 2B and Supplementary Table 2). Patients with mutations in “writers” showed superior overall survival compared with patients lacking mutations (Figure 1B), suggesting they play a functional role in GC. Using GSVA analysis, the underlying mechanisms were related to the activation of tumor suppressor pathways, including P53 signaling, base excision repair, and mismatch repair pathways (Supplementary Figure 2C and Supplementary Table 3). Copy number variation (CNV) frequency and DNA methylation levels were compared, and mRNA expression levels between paired normal and GC tissues were analyzed. Upregulation of most “writers” was observed in tumor tissue than normal tissue (Figure 1G). Through CNV frequency analysis, CNV gains were more frequent than losses for many “writers”, including ADAR, CPSF1, CPSF4, KIAA1429, CSTF1, and CSTF3 (Figure 1C). CNV sites’ locations for all “writers” on the chromosomes can be found in Figure 1D. Some “writers” showed significant upregulation in tumor tissue without a high frequency of CNV gains or CNV losses. Eight “writers” were selected with obvious inconformity between CNV status and mRNA expression. These were divided into four groups based on their CNV status, including normal tissue, tumors with CNV gains, CNV losses, and non-CNV alterations. We compared the expression levels of the eight “writers” between the four groups (Supplementary Figure 2D–2K). Patients with CNV gains showed significantly higher expression of all “writers,” excluding ADARB2. Besides ADARB2, the expression of all “writers” were similarly upregulated in tumors with no CNV alterations and upregulated or not significantly changed in tumors with CNV losses. The expression of ADARB2 was downregulated in all three tumor groups compared with normal samples regardless of CNV status. Thus, expression levels may not be regulated by CNV alterations. DNA methylation analysis for these “writers” identified that the DNA methylation levels of CPSF3 were downregulated, whilst significant upregulation of CFI, ADAR, and ADARB2 in the tumors was observed (Figure 1E, 1F). Considering the mRNA expression of ADARB2 was markedly down regulated in tumors and unrelated to CNV alterations, we hypothesize that the expression of ADARB2 is regulated by hypermethylation in GC. These analyses demonstrate that 26 RNA modification “writers” expression levels are regulated by gene mutations and DNA methylation, which is highly heterogeneous between normal ones and GC tumors.

Figure 1.
The landscape of genetic, epigenetic, and transcriptional alterations of RNA modification “writers” in gastric cancer. (A) Mutational frequency of the 26 RNA modification “writers” in 433 GC patients from the TCGA-STAD cohort. Each column represents individual patients. Upper bar plots show TMB, and numbers on the right indicate the mutational frequency in each “writer”. Right bar plots show the proportion of each variant type. Stacked bar plots below show the fraction of conversions in each sample. (B) Kaplan–Meier curve showing the overall survival status of patients with mutations (red curve) in one or more of the 26 writers or without mutations (blue curve) in the TCGA-STAD cohort. The grouping status of the STAD patients is indicated at the bottom of the chart. P < 0.05 in the Log-rank test was considered statistically significant. (C) CNV variation frequency of the 26 writers in the TCGA-STAD cohort. The height of the column represents the alteration frequency. Red dots represent deletion frequency; blue dots represent amplification frequency. (D) Location of the CNV alterations of RNA modification writers on 23 chromosomes from the TCGA cohort. (E) The methylation level of some writers between normal and gastric cancer tissues. In the box plot, blue represents normal tissues, and red represents cancer tissue. The upper and lower ends of the boxes represent the interquartile ranges of values. Lines in the boxes represent median values. Blue or red dots show outliers. Asterisks above the boxes represent the p-value (*P < 0.05; **P < 0.01; ***P < 0.001). (F) Heatmap of the methylation level of some writers. Methylation levels of these writers increase gradually with the color changes from green to red. (G) Expression levels of the 26 writers were composed of four types of RNA modification between normal and gastric cancer tissues. The description of the box plots is as in (E).

RNA modification patterns by the 26 “writers”

RNA expression profiles of five GEO gastric cancer cohorts were merged with their clinical data into a single meta-GEO cohort containing 1051 GC patients. In the meta-GEO cohort, we have done Univariate Cox regression analysis to investigate the prognostic value of the 26 RNA modification “writers” for GC (Supplementary Figure 2L). A comprehensive landscape of mutual interactions, regulatory connections, and prognostic significance of the 26 “writers” in GC are presented in a regulatory network (Figure 2A). Within the network, a strong correlation between expressions was observed. TCGA-STAD cohort correlation analysis showed a stronger positive correlation of writer expression (Supplementary Figure 2M) amongst four types of RNA modification “writers”. GO BP enrichment analysis for these “writers” showed involvement in RNA processing and modification (Supplementary Figure 2N). Patients were grouped into distinct RNA modification patterns in accordance with the 26 “writers” in the meta-GEO cohort expressions (Supplementary Figure 3A–3F). Eventually, four distinct types of RNA modification patterns, including 329 patients in pattern A, 308 patients in pattern B, 204 patients in pattern C, and 210 patients in pattern D, were identified. These were named WM_Cluster A-D, respectively (Figure 2B, Supplementary Figure 3G, and Supplementary Table 4). The expression of these “writers” varied greatly between different clusters (Supplementary Figure 3H). Patients in the WM_Cluster D had significantly superior prognoses than patients in the other three patterns (Figure 2C).

Figure 2.
RNA modification patterns and different biological characteristics between these patterns. (A) Interactions between RNA modification writers in gastric cancer. Circle size represents the effect of each regulator on prognosis, and the range of values calculated by the Log-rank test was p <1e-04, p < 0.001, p < 0.01, p < 0.05, and p < 1, respectively. (B) Unsupervised clustering of the 26 writers in the meta-GEO cohort. WM_Cluster, GC cohort names, tumor stage, survival status, age, and gender were used as patient annotations. Red represents a high expression of the writers, and blue represents a low expression. (C) Kaplan–Meier curve showing the overall survival status of four types of RNA modification patterns based on 1051 patients of the meta-GEO cohort. WM_Cluster A, blue; WM_Cluster B, orange; WM_Cluster C, red; WM_Cluster D, green. The grouping status of the patients is indicated at the bottom of the chart. P < 0.05 in the Log-rank test was considered statistically significant. (D–F) GSVA enrichment analysis shows differentially activated biological pathways between the RNA modification patterns. Heatmaps were used to visualize these biological processes. Red represents activated pathways, whilst blue represents inhibited pathways. GC cohort names and WM_Cluster were used as patient annotations. (D) WM_Cluster A vs WM_Cluster D; (E) WM_Cluster B vs WM_Cluster D; (F) WM_Cluster C vs WM_Cluster D.

Biological pathways and TME cell infiltration in diverse RNA modification patterns

GSVA analysis was performed for the identification of distinctions in biological pathways between WM_Cluster D and WM_Cluster A-C (Figure 2D–2F and Supplementary Table 6). Pathways of stromal and carcinogenesis activation were inhibited, whilst cell senescence and apoptosis were activated in WM_Cluster D compared with WM_Cluster A-C. ECM receptor interactions were inhibited whilst cell cycle and mismatch repair pathways were activated in WM_Cluster D compared with WM_Cluster A (Figure 2D); TGF-beta signaling, MAPK signaling, and mTOR signaling were suppressed, but P53 signaling, cell cycle, and activation of mismatch repair pathways were observed in WM_Cluster D compared with WM_Cluster B (Figure 2E). Notch signaling and basal cell carcinoma were suppressed whilst P53 signaling, cell cycle, and nucleotide excision repair pathways were activated in the WM_Cluster D compared with the WM_Cluster C (Figure 2F). Surprisingly, samples in the WM_Cluster D were infiltrated with more abundantly immune-active cells, such as CD4 T cells, activated CD8 T cells, activated dendritic cells and macrophage M1 (Figure 3A). Subsequent analysis revealed that an array of immune-related functions was more active in WM_Cluster D (Figure 3B and Supplementary Tables 9, 10). Stromal-related and carcinogenic pathways were also inhibited, and tumor-suppressive and immune-active pathways were activated in the WM_Cluster D (Figure 3C). T cell enhancers were most highly expressed in the WM_Cluster D (Figure 3D). The immune-activation and immune-checkpoint gene expressions were also upregulated, whilst the stromal-activation gene expressions were downregulated in the WM_Cluster D (Supplementary Figure 3J–3L). Patients in the WM_Cluster D showed high immune but low stromal scores (Figure 3E, 3F). The WM_Cluster B was infiltrated by more stromal and eosinophils, MDSCs, mast cells, M2 macrophages, and regulatory T cells (Figure 3A). The WM_Cluster B was strongly associated with stromal activation based on the former GSVA analysis. Stromal-related pathways were significantly activated in the WM_Cluster B (Figure 3C). The WM_Cluster B was also the highest according to stromal score and stromal-related gene expression amongst the four groups (Figure 3F and Supplementary Figure 3L). We speculated that immunosuppressive cell infiltration and stromal activation inhibited antitumor immune responses, so patients in the WM_Cluster B showed no preferential prognosis. WM_Cluster C was the lowest for immune cell infiltration and immune-related functions (Figure 3A–3C). The WM_Cluster C mirrored immune suppression. The WM_Cluster A was intermediate with no significant characteristics of TME infiltration. Correlation analysis between TME infiltration cell type and immune activation showed a high correlation with WTAP, RBM15, NUDT21, CSTF3, and CPSF2 (Supplementary Figure 3I).

Figure 3.
TME characteristics and clinical features of distinct RNA modification patterns. (A) Relative abundance of 24 TME infiltrating cells in four RNA modification patterns. Names of the infiltrating cells are listed at the bottom of the chart. (B) Score of functions in immune regulation of four RNA modification patterns. Immune functions are listed. (C) Relative enrichment of stromal-related, immune-related and tumor-related pathways of the four RNA modification patterns. Pathway names are listed below. (D) Expression of T cell function enhancers. Gene names of these enhancers are listed. In both the box plots of (A–D), the color of WM_Cluster A is blue, WM_Cluster B is orange, WM_Cluster C is red, and WM_Cluster D is green. The upper and lower ends of the boxes represent an interquartile range of values. Lines in the boxes represent the median value. Four colors of the dots are outliers. Asterisks above the boxes represent the statistical p-value (*P < 0.05; **P < 0.01; ***P < 0.001). (E) Immune score and (F) stromal score of four types of RNA modification patterns. P < 0.05 in the Kruskal-Wallis test was regarded as statistically significant. (G–I) Bar plots showing RNA modification patterns in different clinical stages (G), ACRG molecular subtypes (H), and Lauren subtypes (I).

Association between clinical features and RNA modification patterns

We analyzed their association with clinical features. A proportion of tumor stages differed between the four clusters in the meta-GEO cohort, and tumor stages of the patients in the WM_Cluster D were less advanced (Figure 3G). Lauren’s subtypes were also related to the WM_Cluster, and the intestinal subtype of GC counts was predominant in the WM_Cluster D, whilst the diffuse subtype counts were lower. Assessment of the GSE62254 ACRG cohort classified GC into four types: EMT, MSI, MSS with TP53-active, and MSS with TP53-inactive. We found that the proportion of the MSI subtype of tumors was highest whilst the proportion of the EMT subtype was lowest in the WM_Cluster D. These results described how patients in the WM_Cluster D have superior overall prognosis over patients in the other three clusters.

Identifying RNA modification-associated genes and assessment of their functions

We identified 1801 DEGs common between WM_Cluster D and the other three clusters (Figure 4A). The biological functions of these DEGs were assessed. Genes were enriched in pathways related to cell cycle and cell senescence, PD-L1 expression, PD-1 checkpoints, carcinogenesis, and cancer (Figure 4B and Supplementary Table 14). Unsupervised clustering analysis was done based on the expression of the 1801 DEGs were performed to classify patients into different genomic subtypes (Supplementary Figure 4A–4F). We identified four distinct genomic subtypes and termed them gene clusters A-D, respectively (Figure 4C and Supplementary Table 4). Differences in the expression of the 26 “writers” between the four gene clusters were to the four RNA modification patterns (Figure 4D). Patients in gene cluster D had the highest overall survival, consistent with the RNA modification patterns’ prognosis analysis (Figure 4E). Mechanistically, stromal-related and carcinogenic pathways were significantly inactivated, whilst pathways for anti-tumor immune activation, cell senescence, and apoptosis were highly activated in gene cluster D (Figure 4F). T cell function enhancer and immune activation genes were also highly expressed, whilst genes regulating stromal activation were downregulated (Supplementary Figure 4I, 4J). Patients in gene cluster D also benefited from immune-checkpoint block therapy (Figure 4F and Supplementary Figure 4I).

Figure 4.
Construction of gene signatures and functional annotation. (A) Venn diagram showing the 1801 overlapped DEGs between WM_Cluster D and the other three clusters. (B) KEGG enrichment analysis of the overlapped 1801 RNA modification patterns. The length of the bars represents the number of genes enriched. (C) Unsupervised clustering of overlapping RNA modification-related genes in the meta-GEO cohort to classify patients into different genomic subtypes, defined as gene clusters A-D, respectively. Gene clusters, WM_Clusters, GC cohort names, tumor stage, survival status, age, and gender were used as patient annotations. (D) Kaplan–Meier curve showing the overall survival status of the four gene clusters based on the meta-GEO cohort. Gene cluster A, blue; gene cluster B, orange; gene cluster C, red; gene cluster D, green. The grouping status of the patients is indicated at the bottom of the chart. P < 0.05 in the Log-rank test was considered statistically significant. (E) Expression levels of the 26 writers in the four gene clusters. (F) Relative enrichment of some stromal-related, immune-related, and tumor-related pathways of the four gene clusters. For (E) and (F), the upper and lower ends of the boxes represent interquartile ranges of values. Lines in the boxes represent the median value. The four colors of the dots are outliers. Asterisks above the boxes represent the statistical p-value (*P < 0.05; **P < 0.01; ***P < 0.001).

Quantifying RNA modification patterns

We termed the scoring system as a WM_Score (Supplementary Table 4). Connections among RNA modification patterns, gene clusters, WM_Score groups, and ACRG subtypes in the GSE62254 cohort are shown through a Sankey diagram (Figure 5A). Consistent with previous data, we classified the majority of patients in WM_Cluster D into gene cluster D and WM_Score low groups. No patients in the WM_Score low group were classified into the EMT subtype (Figure 5A). The WM_Score was found to vary across different RNA modification patterns. WM_Cluster B had the highest WM_Score, whilst WM_Cluster D had the lowest (Figure 5B). Consistently, the WM_Score was highest in gene cluster B but lowest in gene cluster D (Figure 5C).

Figure 5.
Quantifying RNA modification patterns according to WM_Score. (A) Sankey plot showing the relationship between WM_Clusters, gene cluster, WM_Score and ACRG molecular subtypes. (B, C) WM_Scores according to different RNA modification patterns and gene clusters in the meta-GEO cohort. P < 0.05 in the Kruskal-Wallis test was considered statistically significant. (D, E) Kaplan–Meier curves showing the overall survival status of WM_Score-low (blue) and WM_Score-high (red) groups of patients in the meta-GEO cohort (D) and TCGA-STAD cohort (E), grouping status of the patients are indicated at the bottom of the chart. P < 0.05 in the Log-rank test was considered statistically significant. (F) Kaplan–Meier curves of WM_Score patients with or without adjuvant chemotherapy. H, high; L, low; ADJC, adjuvant chemotherapy. P < 0.001 by Log-rank test. (G, H) Multivariate Cox regression model of WM_Score, patient age, gender, TNM status, and patient outcomes in the meta-GEO cohort (G) and TCGA-STAD cohort (H) shown by the forest plots. (I) Predictive value of the WM_Score in patients amongst the meta-GEO cohort. The AUC of five, six, and seven-year survival was 0.612, 0.625, and 0.633, respectively.

WM_Score associated with genetic and epigenetic alterations

We compared the most frequently mutated genes and known effectors of targeted therapy and found that the mutation rates of all these genes were significantly higher (Figure 7A, 7B). According to TMB quantification analysis, a high WM_Score was significantly associated with a lower TMB (Figure 7C, 7D). Patients with low WM_Scores may therefore benefit more from targeted therapies for GC. Patients with low TMB and high WM_Scores had the worst prognosis than the other three groups (Figure 7E). The distribution of CNVs is shown in Figure 7F, 7G. The general CNV frequency was higher in the low WM_Score group (Figure 7H). The frequency of 9p21.3 deletions was higher in the high WM_Score group, whilst 19q12 amplifications occurred more frequently in the low WM_Score group (Figure 7F, 7G). Two important tumor suppression genes, CDKN2A, and CDKN2B, are transcribed from 9p21.3, the deletion of which contributes to tumorigenesis and tumor development, including GC. Given these data, we considered that TME immune suppression and inactivation of tumor suppressor genes caused by the higher frequency of 9p21.3 deletions might explain the poor prognosis of the low WM_Score group. However, the important oncogene CCNE1 is transcribed from 19q12, so its amplification frequency may be higher in the low WM_Score group. Factors affecting the biological and clinical outcomes were also variable. Expression level analyses for these genes identified CDKN2B and CCNE1 as more significantly highly expressed in the low WM_Score group, consistent with 9p21.3 locus deletion and 19q12 locus amplification (Figure 7I). Unsupervised clustering could divide patients in the TCGA-STAD cohort based on the methylation levels of prognosis-related gene sites (Supplementary Figure 7A–7F). Amongst the three DNA methylation clusters, the WM_Score was unfavorably related to patient prognosis (Figure 7J, 7K). DNA methylation levels of CDH1, DAPK1, and RASSF1 were positively correlated with WM_Score and adversely associated with patient survival, whilst DNA methylation levels of stromal-related TGFB2 negatively correlated with WM_Score but positively associated with patient survival (Figure 7L).

Figure 7.
Relationships between WM_Score and genetic and epigenetic alterations. (A, B) Waterfall plot of tumor somatic mutations established by those with high WM_Score (A) and low WM_Score (B). Each column represents individual patients. Upper bar plots show TMB. The numbers on the right indicate mutation frequency in each gene. The right bar plot shows the proportion of each variant type. (C) Difference of TMB between WM_Score high and low groups (Wilcoxon test). (D) Correlation between WM_Score and TMB (R = -0.46, p < 2.2e-16, Pearson’s Correlation test). (E) Kaplan–Meier curves showing the overall survival status of WM_Score high and low patients with low or high TMB. H, high; L, low. P = 0.02 by Log-rank test. (F, G) Distribution of copy number amplification and deletion sites in the WM_Score high group (F) and WM_Score low group (G). The Abscissa axis represents the location of CNV on the chromosome. The ordinate axis represents the G-score. Amplifications are marked with red, and deletions with blue. (H) CNV frequency in WM_Score high and low groups (p = 0.006, Wilcoxon test). (I) Box plot showing the expression of genes with significant differences in CNVs between the two groups. (J) WM_Scores according to DNA methylation subtypes (Kruskal-Wallis test). (K) Kaplan–Meier curve showing the overall survival status of three methylation subtypes of patients in the TCGA-STAD cohort (p = 0.011, Log-rank test). (L) Box plot showing the methylation level of GC survival-related genes (*P < 0.05; **P < 0.01; ***P < 0.001).

WM_Score-related transcriptional and post-transcriptional regulation

To investigate the effects of the WM_Score on transcriptional and post-transcriptional regulation, we studied the miRNA regulation to bridge the investigation of these two aspects in the TCGA-STAD cohort. We first screened 112 differentially expressed miRNAs (DEmiRNAs) and 8480 mRNAs (differentially expressed) between high and low WM_Score groups. We then predicted the target genes of the differentially expressed miRNAs using the miRTarBase database. KEGG enrichment analysis for the target genes was then performed amongst DEmRNAs. The enrichment of target genes was observed in cell cycle, cell senescence, PD-1 checkpoint, and cancer-related pathways (Figure 8A and Supplementary Table 16). Most genes related to tumor suppressor pathways, such as cell cycle status, cell senescence, P53 signaling, AMPK signaling, and Hippo signaling, were expressed to low levels in the high WM_Score group. Genes associated with the PD-1 checkpoint pathway were also downregulated in the high WM_Score group. These genes were mainly regulated by highly expressed miRNAs in the high WM_Score group (Figure 8A). Genes upregulated in the WM_Score high group were related to carcinogenic pathways, including gastric cancer, MAPK signaling, PI3K-Akt, proteoglycans in cancer, Rap1 signaling, and TGF-beta signaling. These genes are targeted by low-expressing miRNAs in the high WM_Score group (Figure 8A).

Figure 8.
Relationships between WM_Score and transcriptional and post-transcriptional regulation. (A) Differentially expressed miRNAs and their associated pathways between WM_Score high and low groups. Red lines and dots represent highly expressed genes in the WM_Score high group; blue lines and dots represent highly expressed genes in the WM_Score low group. Circles represent targeted mRNA-enriched pathways. (B) Differences in the percentage of Distal poly(A) site Usage Index (PDUI) of the genes between WM_Score high and low groups. Red, PDUI lengthened genes in WM_Score high group; blue, PDUI shortened genes in WM_Score high group; grey, genes with no significant difference in PDUI. (C, D) Kaplan–Meier curves showing the overall survival status of PDUI lengthening (red) and PDUI shortening (blue) of ARF6 (C) and C2CD4A (D). (E) Differences of A-I editing frequency of genes between WM_Score high and low groups. Red, high A-I editing genes in the WM_Score high group; blue, low A-I editing genes in WM_Score high group; grey, genes with no significant difference in A-I editing. (F, G) Kaplan–Meier curves showing the overall survival status of A-I editing high (red) and A-I editing low (blue) of gene ECH1 (F) and SAMHD1 (G).

In this study, we identified genes with significantly diverse 3’UTR lengths caused by APA in the high and low WM_Score groups and analyzed their association with patient prognosis (Figure 8B and Supplementary Table 17). We found that ARF6 and C2CD4A showed significant 3’UTR shortening in the WM_Score high group (Supplementary Table 17), which was associated with poor prognosis (Figure 8C, 8D). In our analyses, we speculated that the shortening of mRNA 3’UTRs of these genes in the high WM_Score group could decrease miRNA targeting efficacy and upregulate the expression of these oncogenes, contributing to tumorigenesis and GC development. We also screened out genes with significant differential A-to-I editing rates and analyzed their association with patient prognosis (Figure 8E and Supplementary Table 18). The A-to-I editing rate of ECH1 was higher, whilst that of KMT2C was lower in the high WM_Score group (Supplementary Table 18), which was significantly related to poorer survival (Figure 8F, 8G).

Predictive value of the WM_Score model in immunotherapy

In the Kim et al. cohort, patients in the response group had a lower WM_Score (Figure 9A, 9B). Consistent with his study, EBV-positive and MSI-H subtypes of GC were lower according to WM_Score in this study (Figure 9C). This GC cohort was the only sample group receiving ICB therapy from which gene expression profiles and the responses of immunotherapy could be obtained. We could not identify the survival information of the patients and were unable to make predictions for patient prognosis. We adopted three cohorts of ICB immunotherapy for other cancer types with accessible survival information. Patients in the WM_Score low group exhibited better overall prognosis in both the IMvigor210 cohort (Figure 9D), Liu et al. cohort (Figure 9K), and GSE78220 cohort (Figure 9O). Patients with significant therapeutic advantages and clinical responses to anti-PD-1/L1 immunotherapy also tended to have lower WM_Score in both cohorts (Figure 9E, 9F, 9P). To explore the mechanisms of the superior clinical outcomes in immunotherapy for the low WM_Score group, we analyzed TME cell infiltration, immune and stromal-related gene expression, and pathway enrichment in the two groups in the IMvigor210 cohort. Consistent with our previous findings, the low WM_Score group was infiltrated with a higher number of immune-active cells, such as activated CD4 and CD8 cells and M1 macrophages, but fewer stromal-activation cells such as eosinophils, MDSCs, mast cells and regulatory T cells (Supplementary Figure 8A). An array of T cell enhancer genes and immune activation and immune checkpoint genes were up regulated in the low WM_Score group (Supplementary Figure 8B, 8D, 8E), while stromal-activation-related genes were downregulated (Supplementary Figure 8F). Stromal-related and carcinogenic pathways were inhibited, whilst pathways related to immune activation and cell senescence were activated in the low WM_Score group (Supplementary Figure 8C). Moreover, the inflamed immune phenotype of the tumors had the lowest WM_Score compared with immune excluded and immune desert phenotypes (Supplementary Figure 8G). The expression of PD-L1 on tumor (TC, detected by SP142) and immune cells (IC, detected by SP142) correlated with the WM_Score. TC1 and TC2+ had lower WM_Scores than TC0, whilst IC2+ had significantly lower WM_Scores than IC0 and IC1 (Supplementary Figure 8H, 8I). Patients in the low WM_Score group showed higher expression levels of immune-checkpoint genes PD-L1 and CTLA4 (Figure 9M, 9N). This was consistent with the efficacy of immunotherapy, which was positively associated with PD-L1 expression. A low WM_Score was related to higher tumor neoantigen burden and tumor mutational burden (Figure 9G, 9H). Patients with low WM_Scores and high TMB or high neoantigen levels had a higher survival advantage over other groups (Figure 9I, 9J). The WM_Score is a robust model for clinical responses in immunotherapy (Figure 9Q). The mechanisms underlying its predictive ability were related to TME cell infiltration and PD-L1/PD-1-related pathway enrichment.

Figure 9.
Relationship between WM_Score and efficacy of immunotherapy. (A, B) Violin plots showing the WM_Score in patients with different responses to PD-L1 blockade in the Kim et al. immunotherapy cohort. Kruskal-Wallis tests were used to evaluate the significance. (C) Box plot showing the WM_Score in patients of different TCGA subtypes (Kruskal-Wallis test). (D) Kaplan–Meier curve showing the overall survival status of patients receiving PD-L1 blockade immunotherapy in the IMvigor210 cohort (p = 0.003, Log-rank test). (E) The proportion of patients responding to PD-L1 blockade immunotherapy in low and high WM_Score groups in the IMvigor210 cohort. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. (F) Violin plot showing the WM_Score in patients with different responses to PD-L1 blockade immunotherapy in the IMvigor210 cohort. (G, H) Box plots showing the neoantigen burden (G) and mutational burden (H) of high and low WM_Score groups in the IMvigor210 cohort. Wilcoxon tests were used to evaluate the significance. (I) Kaplan–Meier curve showing the overall survival of high or low WM_Score patients with high or low TMB (p = 0.007, Log-rank test). (J) Kaplan–Meier curves showing the overall survival of high or low WM_Score patients with high or low neoantigen burden in the IMvigor210 cohort (p = 0.005, Log-rank test). (K) Kaplan–Meier curve showing the overall survival status of patients receiving PD1 blockade immunotherapy in high and low WM_Score groups in the Liu et al. immunotherapy cohort (p = 0.018, Log-rank test). (L) The proportion of patients responding to PD1 blockade immunotherapy in low and high WM_Score groups. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. (M, N) Box plots showing the expression of PD-L1 (M) and CTLA4 (N) of high and low WM_Score groups. Wilcoxon tests were used to evaluate significance. (O) Kaplan–Meier curve showing the overall survival status of patients receiving anti-PD1 immunotherapy in high and low WM_Score groups in the GSE78220 cohort (p = 0.053, Log-rank test). (P) The proportion of patients responding to anti-PD1 immunotherapy in low and high WM_Score groups in the GSE78220 cohort. PD, progressive disease; CR, complete response; PR, partial response. (Q) Predictive value of the WM_Score in patients receiving anti-PD-1 immunotherapy in the GSE78220 cohort (AUC = 0.827).

Validating key genes of RNA modification and TME at the protein level

We analyzed GC immunohistochemistry images of three RNA modification “writers”, two immune activation-related genes, one stromal activation-related gene, and two immune checkpoints of three patients in the HPA database. Each gene among the three patients was stained using identical antibodies. CSTF3, TRMT6, and ZC3H13 were higher expressed in WM_Cluster D (Supplementary Figure 3H), which was associated with immune activation, stromal inhibition, and activation of PD-L1/PD-1 checkpoints (Figure 3C). In the images, CSTF3, TRMT6, and ZC3H13 expression were higher in the GC tissues of patient 2473 compared to patients 2626 and 3044 (Figure 10A–10C). Therefore, the RNA modification patterns of patient 2473 were more closely to WM_Cluster D. TNF, and GZMA were more highly expressed in patient 2473, but ACTA2 was expressed at lower levels. Surprisingly, the expression of PD-L1 and IDO1 was only detected in patient 2473 (Figure 10A–10C).

Figure 10.
Staining patterns of RNA modification writers, immune-related genes, stromal-related genes, and immune checkpoints. (A–C) Immunohistochemical images of patient ID 2473 (A), 2626 (B), and 3044 (C). Images were downloaded from The Human Protein Atlas (HPA). Genes, antibodies, and staining degree and intensity are listed on the right of each image.

Exploring WM_Score performance across pan-cancer types

The WM_Score correlated with the infiltrating fraction of 22 immune cell types in nearly all cancer types except for ACC, and the correlation trend for each infiltrating immune cell type varied amongst diverse cancer types (Figure 11A). The correlation trends between the WM_Score and stemness indices were also different across pan-cancer types but generally positively correlated with cancer stemness (Figure 11B). The WM_Score could also predict the overall survival and prognosis-free interval for most cancer types and represented a risk factor for ESCA, LAML, LIHC, LUSC, and UCEC (Figure 11B, 11C). The WM_Score in these cancer types negatively correlated with immune activation-related cell infiltration and positively correlated with stromal activation-related cell infiltration, such as Tregs, follicular helper T cells, and M2 macrophages (Figure 11A) and stemness indices (Figure 11B). TMB, MSI, and PD-L1 expression are biomarkers to evaluate the efficacy of ICB immunotherapy. Correlations between the WM_Score and TMB, MSI, and PD-L1 expression in pan-cancer types are shown using radar charts (Figure 11E–11G).

Figure 11.
Performance of the WM_Score across pan-cancer types. (A, B) Correlation between WM_Score and immune cell fractions (A) and stemness indices (B) for each cancer type. Upper regions of grids show p-values. The bottom regions show the correlation coefficient. (C, D) Overall survival (C) and prognosis-free interval (D) analyses for the WM_Score in TCGA pan-cancer types using a univariate Cox regression model. A hazard ratio < 1 represents protective factors for survival, while hazard ratio > 1 represents risk factors for survival. (E–G) Radar plots of the correlation between WM_Score and tumor mutation burden (E), microsatellite instability (F), and PD-L1 expression level (G). Dots in the radar plots represent the R-value of correlation: R > 0, positive correlation, and R < 0, negative correlation. Asterisks or the black dots represent the statistical p-value (Pearson’s Correlation test, *P < 0.05; **P < 0.01; ***P < 0.001).

Discussion

RNA modifications act in TME immune cell infiltration and the response to immunotherapy. These RNA modifications occur due to the activity of “writers,” which play a key role in tumor immunology and therapy. Previous reports have focussed on the activity of a single “writer” or one RNA modification type. In this study, we 4 major types of RNA modification “writers” and comprehensively analyzed their roles in TME cell infiltration, pathway enrichment, and clinical outcomes of patients with or without ICB immunotherapy in GC.

We identified frequent alterations in genetic, epigenetic, and RNA expression levels of these “writers” in GC. We further identified 4 RNA modification patterns according to the expression levels of the 26 “writers”. We found WM_Cluster D had superior overall survival over the other three RNA modification patterns. The WM_Cluster D was associated with higher immune-activation-related adaptive immune cell infiltration, characterized as the immune-inflamed phenotype of cancer. The WM_Cluster D was more active in CD8 T effectors, T cell receptors, APM, PD-L1 expression, and PD-1 checkpoint pathways. The WM_Cluster B was related to higher stromal and innate immune cell infiltration, characterized as an immune-excluded phenotype of cancer. EMT1/2/3, Pan-F-TBRS, angiogenesis pathways, and stromal-activation-related gene expression were more activated in WM_Cluster B, which validated the TME infiltration analysis, and explained why the WM_Cluster B conferred no survival advantage. WM_Cluster C was associated with low infiltration of nearly all immune cells, characterized as an immune-desert phenotype. No survival advantage was observed due to weak anti-tumor immunology. WM_Cluster A was intermediate with no significant survival characteristics. The expression levels of the 26 RNA modification “writers” were significantly distinct in the WM_Cluster D, and most of the “writers” were more highly expressed in the WM_Cluster D compared with the other three clusters. Some of the “writers” have been shown to facilitate T cell functions. For example, WTAP-mediated m6A modifications could prevent T cells from T cell antigen-receptor (TCR) signaling-induced cell death and promote the activation of T cells by destabilizing ORAI1 and RIPK1 mRNAs [24]. m1A “writer” TRMT61A and TRMT6 can strengthen MYC protein synthesis, which guides naive T cells from a quiescent state into a proliferative state, promoting rapid T cell expansion [25]. However, whether these “writers” were upregulated in the WM_Cluster D and could facilitate the anti-tumor immune functions of T cells and the mechanisms of TME immune activation by the upregulation of these “writers” in GC are rarely reported and require further exploration.

Further analysis showed that the DEGs in WM_Cluster D were also enriched in many cancer-related pathways and PD-L1 expression. Four gene clusters were identified based on the expression of DEGs. This demonstrated how RNA modifications were of great importance in shaping the distinct TME landscape and prognosis of GC. Given that RNA modifications are heterogeneous in individual patients, we constructed a WM_Score model to analyze the relationship between RNA modification patterns, TME landscape, and patient prognosis. High WM_Scores were associated with stromal activation, carcinogenic pathway activation, enhanced cancer stemness, and poor clinical outcomes, whilst low WM_Scores were associated with immune activation, tumor suppression, and improved clinical outcomes. These results were not observed in patients without ICB immunotherapy in the meta-GEO and TCGA-STAD cohort but were validated in four immunotherapeutic cohorts. TGFβ-activated stroma in TME inhibits T-cell responses against tumor cells and tumor susceptibility to anti-PD-1-PD-L1 therapy [26]. Many carcinogenic pathways, including TGF-β, promote the self-renewal and growth of cancer stem cells [27]. From these results. A close relationship between WM_Score, TME characteristics, cancer stemness, ICB immunotherapy, WM_Score, and prognosis of GC patients with and without immunotherapy was observed. Except for immunotherapy, we found that a low WM_Score and adjuvant chemotherapy had a superimposed effect in improving patient prognosis. Some RNA modification “writers” can regulate the sensitivity to chemotherapeutic drugs for certain types of cancer by mediating RNA modifications of some genes. The METTL3/YTHDF2 axis can lead to cisplatin resistance of ovarian cancer by regulating the mRNA stability of IFFO1 in an m6A-dependent manner and inhibiting IFFO1 expression [28]. METTL14 upregulates pri-miR-19a m6A and increases NSLCC resistance to cisplatin through the miR-19a-5p/RBM24/AXIN1 axis [29]. Similar mechanisms mediated by RNA modifications that impact tumor sensitivity to chemotherapy have also been observed in breast cancer [30], bladder cancer [31], and colorectal cancer [32].

The WM_Score is related to many genetic and epigenetic alterations, which are associated with TME features and immunotherapy efficacy. WM_Score negatively relates to TMB, and mutation frequencies of MUC16, KMT2D, and PIK3CA are higher in the low WM_Score group. These gene mutations are associated with high anti-tumor immune cell infiltration, improved overall survival, and enhanced efficacy of PD-L1/PD-1 blockers in GC [33–35]. In the CNV analysis, the frequency of 9p21.3 deletions was higher in the high WM_Score group, which is reported to be related to immune suppression in TME, a poor response to ICB immunotherapy, and poorer prognosis [36, 37]. DAPK1 can enhance anti-tumor immunology, and DNA is hypermethylated in high WM_Score patients, whilst stromal-activation-related genes such as TGFB2 are hypomethylated.

The activation of MAPK, PI3K-Akt, and TGF-beta pathways was done in the high WM_Score high group because of low miRNA expression that target these signaling components. Assessment of the WM_Score model across pan-cancer types identified different predictive values amongst cancer types. Thus, the application of the WM_Score model is cancer-specific and requires further exploration into individual cancer types.

Materials and Methods

Collection and pre-processing of attainable expression datasets

The workflow is shown in Supplementary Figure 1. Public gene expression results, and relative clinical information of the patients were downloaded from the Gene-Expression Omnibus and Cancer Genome Atlas databases. Patients lacking survival information were removed for further analysis. In total, 5 gastric cancer cohorts in GEO (GSE15459, GSE34942, GSE57303, GSE62254, and GSE84437) and TCGA-STAD, including 1422 patients, were included. Basic information on each dataset is shown in Supplementary Table 1. The “ComBat” algorithm of the sva package [38] was used for the correction of batch effects from non-biological technical bias by R software (version 4.1.3) and R Bioconductor packages.

Acquisition of the DNA methylation data

DNA methylation data were got from the Cancer Genome Atlas database. The “Methylation Beta Value” was selected. Data from Illumina human methylation 27 and 450 platforms were included.

Unsupervised clustering of 26 RNA modification writers

An unsupervised clustering algorithm was applied for the cluster analysis of 1051 gastric cancer patients. These RNA “writers” consisted of seven m6A modification enzymes (METTL3, METTL14, WTAP, RBM15, RBM15B, ZC3H13 and KIAA1429), four m1A modification enzymes (TRMT61A, TRMT61B, TRMT10C and TRMT6), twelve APA modification enzymes (CPSF1/2/3/4, CSTF1/2/3, PCF11, CFI, CLP1, NUDT21 and PABPN1) and 3 A-I modification enzymes (ADAR, ADARB1 and ADARB2). The consensus clustering algorithm was used to determine the number of clusters. Stability was assessed using the ConsensusClusterPlus package.

GSVA and functional annotation

The R package “GSVA” was used to investigate differences in biological signaling and RNA modification patterns [39]. The gene set “c2.cp.kegg.v7.4.symbols” for GSVA analysis was obtained from the MSigDB database. Adjusted P-values ≤ 0.05 were known as significant. We used the R package “clusterProfiler” for performing functional annotation for RNA modification-related genes with |log2FC| > 0.1 and adjusted P-values ≤ 0.05 [40].

Evaluation of TME cell infiltration

The ssGSEA (single-sample gene-set enrichment analysis) algorithm was used to quantify the relative abundance of each cell infiltrate in the GC TME. Gene sets for marking each TME infiltration cell type were obtained from Charoentong et al., which consisted of various human immune cell subtypes, including activated B cells, activated CD4 T cells, activated CD8 T cells, activated dendritic cells, macrophages, mast cells, natural killer T cells, and regulatory T cells (Supplementary Table 7) [41, 42]. Enrichment scores were calculated using ssGSEA analysis and represent the relative abundance of each TME infiltrating cell per sample (Supplementary Tables 8, 19).

Association between RNA modification and cell signaling

Correlation analysis was performed to reveal the association between RNA modification patterns and cancer cell signaling. Gene sets of related biological pathways involved in this analysis were obtained from Mariathasan et al., including EMT markers EMT1/2/3, pan-fibroblast TGFb response signature (Pan-F-TBRS), Angiogenesis, antigen processing machinery (APM), CD8 T effectors and immune checkpoints [43]. Other gene sets used were from significantly activated signaling pathways in the GSVA analysis and DEGs between different WM_Cluster enriched pathways following KEGG analysis (Supplementary Tables 11, 12, 20).

Immune and stromal scores

We used the ESTIMATE method to calculate stromal and immune scores and to analyze tumor purity for each patient using the R package “estimate”.

Identification of DEGs

We subdivided RNA modification patterns according to the expression of 26 RNA modification “writers”. We also used the R package “limma”‘ the empirical Bayesian approach [44]. P-values ≤ 0.005 were deemed significant criteria for determining DEGs (Supplementary Table 13).

KEGG enrichment analysis

We performed DEGs’ enrichment analysis and functional annotation by means of the R packages “clusterProfiler” and “org.hs.eg.db”. R packages “enrichplot” and “ggplot2” were used to construct bar plots following KEGG analysis.

WM_Score scoring system construction

PCA was used for the construction of the WM_Score scoring model. Establishment of a matrix for the expression of DEGs was done in all patients. We selected PC1 and PC2 as signature scores. WM_Scores were calculated for each patient using a comparable method to GGI [45]:

WM_Score=∑(PC1i + PC2i)

where i denotes the expression of DEGs between different RNA modification patterns. A cutoff value was determined using the R package “survminer” according to WM_Scores and the survival status of the patients. We divided patients into high and low WM_Score groups in accordance with the cutoff value.

Acquisition of stemness indexes

Profiles of the stemness index were downloaded from previous studies [46]. Four stemness indices were assessed through the analysis of transcriptomic and epigenetic feature sets.

Analysis of Copy number variant (CNV)

GISTIC 2.0 was employed to detect the common copy number alterations in all samples based on SNP6 CopyNumber segment data [48]. Masked Copy Number Segment data from the TCGA was first downloaded and disposed of using an R code. Disposed of data was uploaded onto the online analysis tool GenePattern (https://cloud.genepattern.org/gp/pages/login.jsf) for GISTIC analysis with confidence levels set at 0.95. Finally, the R package “maftools” was applied to visualize the GISTIC analysis.

Relationship between WM_Score and miRNA regulation

Data regarding miRNA expression was downloaded from the TCGA-STAD cohort. DEmiRNAs between high and low WM_Score groups were determined using the R package “edgeR”. MiRNAs with |log2FC| >1 and adjusted P-values ≤ 0.05 were considered significantly differentially expressed. We screened DEGs with |log2FC| >1 and adjusted P-values ≤ 0.05 between WM_Score high and low groups using “edge R”. The miRTarBase (https://miRTarBase.cuhk.edu.cn/) database was used to predict target genes amongst the DEGs of DEmiRNAs [49]. KEGG enrichment analysis was performed. Finally, R packages “magrittr”, “tidyverse”, “crosslink” and “ggplot2” were employed to construct a circLink plot to directly reveal differences in DEmiRNA-targeted signaling pathways between high and low WM_Score GC patients.

Association between WM_Score and APA events

APA data of STAD were obtained from the synapse database (https://www.synapse.org/#!Synapse:syn11511914). Changes in APA usage in each tumor were quantified as a change in the Percentage of Distal polyA site Usage Index (PDUI), which can identify 3’UTR lengthening (positive index) or shortening (negative index). The R package “edgeR” was used to distinguish transcripts with significantly differential PDUI values between high and low WM_Score patients.

Association between WM_Score and A-to-I editing

The profile of A-to-I RNA editing frequency in STAD was obtained from the synapse database (https://www.synapse.org/#!Synapse:syn4382531). We used the R package “edgeR” to screen sites with significantly different A-to-I RNA editing frequencies between high and low WM_Score patients.

Data from immunotherapeutic cohorts

Gene expression profiles and clinical information of immunotherapeutic cohorts that could be publicly obtained were analyzed. Four immunotherapeutic cohorts were included: (1) pembrolizumab treatment for metastatic gastric cancer [50]; (2) IMvigor210 cohort, mUC treated with atezolizumab [43]; (3) melanoma patients treated with anti-PD1 ICB Nivolumab or Pembrolizumab [51]; (4) GSE78220 cohort, melanomas undergoing pembrolizumab therapy [52]. Basic information about these cohorts is listed in Supplementary Table 1.

Validation of expression levels

Expression analysis was performed. Data were downloaded from reverse-phase protein arrays (RPPAs) of the 392 STAD patients in the TCGA. We performed correlation analysis between biological signaling pathways and the WM_Score [53]. Pathologic immunohistochemistry images of RNA modification “writers”, immune-activation-related genes, TGFβ-EMT pathway-related genes, and immune-checkpoint-related genes of GC patients were further performed using the Human Protein Atlas (HPA, https://www.proteinatlas.org/). Staining of each gene within the pathological tissue was obtained using the same antibodies [54].

Statistical analysis

All statistical analyses were generated by R version 4.1.3. The R package “maftools” was utilized to plot waterfall maps of the mutation landscape of 26 “writers” in the whole TCGA-STAD cohort and genes for targeted therapy in high and low WM_Score groups. R package “pheatmap” was used to generate heat maps of the methylation profile of the “writers”. R package “RCircos” was used to present the copy number variation landscape of 26 RNA modification “writers” on 23 pairs of chromosomes. Spearman correlation analysis was applied to calculate the correlation coefficient of expression of these “writers” and correlation tests to evaluate the P-values. Wilcoxon tests were used to perform comparisons between the two groups, whilst Kruskal-Wallis tests were performed for the analysis of three or more groups. The Univariate Cox regression model was adopted to calculate hazard ratios (HR) for single RNA modification “writers” in GC and for WM_Scores across pan-cancer types. The R package “forestplot” was used for generating forest maps for visualization of the data. We used the Kaplan-Meier method. R package “survminer” was used to determine the cut-off point of subgroups for each dataset. The function “surv-cutpoint” was used to dichotomize the WM_Score and divide patients into high and low WM_Score groups. WM_Score’s specificity and sensitivity were assessed through ROC curves generated using the R package “survivalROC”.

Supplementary Materials

Supplementary Figures
Supplementary Table 1
Supplementary Table 2
Supplementary Table 3
Supplementary Table 4
Supplementary Table 5
Supplementary Table 6
Supplementary Table 7
Supplementary Table 8
Supplementary Table 9
Supplementary Table 10
Supplementary Table 11
Supplementary Table 12
Supplementary Table 13
Supplementary Table 14
Supplementary Table 15
Supplementary Table 16
Supplementary Table 17
Supplementary Table 18
Supplementary Table 19
Supplementary Table 20

Author Contributions

ZQL designed this work, analyzed the data, and wrote the original draft. XHZ and WJW disposed of the data. GZ and QWR edited the figures. YT supervised this program and revised the draft. ZQL, XHZ, WJW, and YT validated the data. All authors read and approved the submitted version.

Acknowledgments

We are very thankful to Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) databases for providing the transcriptome and clinical information.

Conflicts of Interest

There are no conflicts of interest regarding the publication of this paper.

Funding

This work was supported by the Clinical Research Fund of Shandong Medical Association Qilu Special Project [YXH2022ZX02016, Yuan Tian].

References

  • 1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021; 71:209–49. https://doi.org/10.3322/caac.21660 [PubMed]
  • 2. Chia NY, Tan P. Molecular classification of gastric cancer. Ann Oncol. 2016; 27:763–9. https://doi.org/10.1093/annonc/mdw040 [PubMed]
  • 3. Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer. 2020; 20:303–22. https://doi.org/10.1038/s41568-020-0253-2 [PubMed]
  • 4. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019; 112:108613. https://doi.org/10.1016/j.biopha.2019.108613 [PubMed]
  • 5. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020; 19:53. https://doi.org/10.1186/s12943-020-01170-0 [PubMed]
  • 6. Yue B, Song C, Yang L, Cui R, Cheng X, Zhang Z, Zhao G. METTL3-mediated N6-methyladenosine modification is critical for epithelial-mesenchymal transition and metastasis of gastric cancer. Mol Cancer. 2019; 18:142. https://doi.org/10.1186/s12943-019-1065-4 [PubMed]
  • 7. Wang Q, Chen C, Ding Q, Zhao Y, Wang Z, Chen J, Jiang Z, Zhang Y, Xu G, Zhang J, Zhou J, Sun B, Zou X, Wang S. METTL3-mediated m6A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut. 2020; 69:1193–205. https://doi.org/10.1136/gutjnl-2019-319639 [PubMed]
  • 8. Safra M, Sas-Chen A, Nir R, Winkler R, Nachshon A, Bar-Yaacov D, Erlacher M, Rossmanith W, Stern-Ginossar N, Schwartz S. The m1A landscape on cytosolic and mitochondrial mRNA at single-base resolution. Nature. 2017; 551:251–5. https://doi.org/10.1038/nature24456 [PubMed]
  • 9. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017; 18:31–42. https://doi.org/10.1038/nrm.2016.132 [PubMed]
  • 10. Dominissini D, Nachtergaele S, Moshitch-Moshkovitz S, Peer E, Kol N, Ben-Haim MS, Dai Q, Di Segni A, Salmon-Divon M, Clark WC, Zheng G, Pan T, Solomon O, et al. The dynamic N(1)-methyladenosine methylome in eukaryotic messenger RNA. Nature. 2016; 530:441–6. https://doi.org/10.1038/nature16998 [PubMed]
  • 11. Wang Y, Wang J, Li X, Xiong X, Wang J, Zhou Z, Zhu X, Gu Y, Dominissini D, He L, Tian Y, Yi C, Fan Z. N1-methyladenosine methylation in tRNA drives liver tumourigenesis by regulating cholesterol metabolism. Nat Commun. 2021; 12:6314. https://doi.org/10.1038/s41467-021-26718-6 [PubMed]
  • 12. Su Z, Monshaugen I, Wilson B, Wang F, Klungland A, Ougland R, Dutta A. TRMT6/61A-dependent base methylation of tRNA-derived fragments regulates gene-silencing activity and the unfolded protein response in bladder cancer. Nat Commun. 2022; 13:2165. https://doi.org/10.1038/s41467-022-29790-8 [PubMed]
  • 13. Zhang Y, Liu L, Qiu Q, Zhou Q, Ding J, Lu Y, Liu P. Alternative polyadenylation: methods, mechanism, function, and role in cancer. J Exp Clin Cancer Res. 2021; 40:51. https://doi.org/10.1186/s13046-021-01852-7 [PubMed]
  • 14. Lai DP, Tan S, Kang YN, Wu J, Ooi HS, Chen J, Shen TT, Qi Y, Zhang X, Guo Y, Zhu T, Liu B, Shao Z, Zhao X. Genome-wide profiling of polyadenylation sites reveals a link between selective polyadenylation and cancer metastasis. Hum Mol Genet. 2015; 24:3410–7. https://doi.org/10.1093/hmg/ddv089 [PubMed]
  • 15. Zhu Y, Zhang R, Zhang Y, Cheng X, Li L, Wu Z, Ding K. NUDT21 Promotes Tumor Growth and Metastasis Through Modulating SGPP2 in Human Gastric Cancer. Front Oncol. 2021; 11:670353. https://doi.org/10.3389/fonc.2021.670353 [PubMed]
  • 16. Nishikura K. Functions and regulation of RNA editing by ADAR deaminases. Annu Rev Biochem. 2010; 79:321–49. https://doi.org/10.1146/annurev-biochem-060208-105251 [PubMed]
  • 17. Chan TH, Qamra A, Tan KT, Guo J, Yang H, Qi L, Lin JS, Ng VH, Song Y, Hong H, Tay ST, Liu Y, Lee J, et al. ADAR-Mediated RNA Editing Predicts Progression and Prognosis of Gastric Cancer. Gastroenterology. 2016; 151:637–50.e10. https://doi.org/10.1053/j.gastro.2016.06.043 [PubMed]
  • 18. Bartosovic M, Molares HC, Gregorova P, Hrossova D, Kudla G, Vanacova S. N6-methyladenosine demethylase FTO targets pre-mRNAs and regulates alternative splicing and 3'-end processing. Nucleic Acids Res. 2017; 45:11356–70. https://doi.org/10.1093/nar/gkx778 [PubMed]
  • 19. You XJ, Zhang S, Chen JJ, Tang F, He J, Wang J, Qi CB, Feng YQ, Yuan BF. Formation and removal of 1,N6-dimethyladenosine in mammalian transfer RNA. Nucleic Acids Res. 2022; 50:9858–72. https://doi.org/10.1093/nar/gkac770 [PubMed]
  • 20. Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, Cheng T, Gao M, Shu X, Ma H, Wang F, et al. VIRMA mediates preferential m6A mRNA methylation in 3'UTR and near stop codon and associates with alternative polyadenylation. Cell Discov. 2018; 4:10. https://doi.org/10.1038/s41421-018-0019-0 [PubMed]
  • 21. Xiang JF, Yang Q, Liu CX, Wu M, Chen LL, Yang L. N6-Methyladenosines Modulate A-to-I RNA Editing. Mol Cell. 2018; 69:126–35.e6. https://doi.org/10.1016/j.molcel.2017.12.006 [PubMed]
  • 22. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol. 2020; 17:807–21. https://doi.org/10.1038/s41423-020-0488-6 [PubMed]
  • 23. Wang G, Xie Z, Su J, Chen M, Du Y, Gao Q, Zhang G, Zhang H, Chen X, Liu H, Han L, Ye Y. Characterization of Immune-Related Alternative Polyadenylation Events in Cancer Immunotherapy. Cancer Res. 2022; 82:3474–85. https://doi.org/10.1158/0008-5472.CAN-22-1417 [PubMed]
  • 24. Ito-Kureha T, Leoni C, Borland K, Cantini G, Bataclan M, Metzger RN, Ammann G, Krug AB, Marsico A, Kaiser S, Canzar S, Feske S, Monticelli S, et al. The function of Wtap in N6-adenosine methylation of mRNAs controls T cell receptor signaling and survival of T cells. Nat Immunol. 2022; 23:1208–21. https://doi.org/10.1038/s41590-022-01268-1 [PubMed]
  • 25. Liu Y, Zhou J, Li X, Zhang X, Shi J, Wang X, Li H, Miao S, Chen H, He X, Dong L, Lee GR, Zheng J, et al. tRNA-m1A modification promotes T cell expansion via efficient MYC protein synthesis. Nat Immunol. 2022; 23:1433–44. https://doi.org/10.1038/s41590-022-01301-3 [PubMed]
  • 26. Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia-Ramentol J, Iglesias M, Sevillano M, Ibiza S, Cañellas A, Hernando-Momblona X, Byrom D, Matarin JA, Calon A, et al. TGFβ drives immune evasion in genetically reconstituted colon cancer metastasis. Nature. 2018; 554:538–43. https://doi.org/10.1038/nature25492 [PubMed]
  • 27. Yang L, Shi P, Zhao G, Xu J, Peng W, Zhang J, Zhang G, Wang X, Dong Z, Chen F, Cui H. Targeting cancer stem cell pathways for cancer therapy. Signal Transduct Target Ther. 2020; 5:8. https://doi.org/10.1038/s41392-020-0110-5 [PubMed]
  • 28. Zhang Y, Qiu JG, Jia XY, Ke Y, Zhang MK, Stieg D, Liu WJ, Liu LZ, Wang L, Jiang BH. METTL3-mediated N6-methyladenosine modification and HDAC5/YY1 promote IFFO1 downregulation in tumor development and chemo-resistance. Cancer Lett. 2023; 553:215971. https://doi.org/10.1016/j.canlet.2022.215971 [PubMed]
  • 29. Gong S, Wang S, Shao M. Mechanism of METTL14-mediated m6A modification in non-small cell lung cancer cell resistance to cisplatin. J Mol Med (Berl). 2022; 100:1771–85. https://doi.org/10.1007/s00109-022-02268-2 [PubMed]
  • 30. Jing L, Lan L, Mingxin Z, Zhaofeng Z. METTL3/LINC00662/miR-186-5p feedback loop regulates docetaxel resistance in triple negative breast cancer. Sci Rep. 2022; 12:16715. https://doi.org/10.1038/s41598-022-20477-0 [PubMed]
  • 31. Wei W, Sun J, Zhang H, Xiao X, Huang C, Wang L, Zhong H, Jiang Y, Zhang X, Jiang G. Circ0008399 Interaction with WTAP Promotes Assembly and Activity of the m6A Methyltransferase Complex and Promotes Cisplatin Resistance in Bladder Cancer. Cancer Res. 2021; 81:6142–56. https://doi.org/10.1158/0008-5472.CAN-21-1518 [PubMed]
  • 32. Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, Wang G, Ni D, Xun Z, Niu X, Ye Y, Li HB. Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer. 2021; 20:29. https://doi.org/10.1186/s12943-021-01322-w [PubMed]
  • 33. Li X, Pasche B, Zhang W, Chen K. Association of MUC16 Mutation With Tumor Mutation Load and Outcomes in Patients With Gastric Cancer. JAMA Oncol. 2018; 4:1691–8. https://doi.org/10.1001/jamaoncol.2018.2805 [PubMed]
  • 34. Li Z, Jia Y, Zhu H, Yuan H, Xing X, Xin Y, Ma T, Pang F, Zhang Y, Hu Y, Jia S, Ji J. Genomic landscape of microsatellite instability in Chinese tumors: A comparison of Chinese and TCGA cohorts. Int J Cancer. 2022; 151:1382–93. https://doi.org/10.1002/ijc.34119 [PubMed]
  • 35. Zhao R, Wan Q, Wang Y, Wu Y, Xiao S, Li Q, Shen X, Zhuang W, Zhou Y, Xia L, Song Y, Chen Y, Yang H, Wu X. M1-like TAMs are required for the efficacy of PD-L1/PD-1 blockades in gastric cancer. Oncoimmunology. 2020; 10:1862520. https://doi.org/10.1080/2162402X.2020.1862520 [PubMed]
  • 36. Spiliopoulou P, Yang SYC, Bruce JP, Wang BX, Berman HK, Pugh TJ, Siu LL. All is not lost: learning from 9p21 loss in cancer. Trends Immunol. 2022; 43:379–90. https://doi.org/10.1016/j.it.2022.03.003 [PubMed]
  • 37. Han G, Yang G, Hao D, Lu Y, Thein K, Simpson BS, Chen J, Sun R, Alhalabi O, Wang R, Dang M, Dai E, Zhang S, et al. 9p21 loss confers a cold tumor immune microenvironment and primary resistance to immune checkpoint therapy. Nat Commun. 2021; 12:5606. https://doi.org/10.1038/s41467-021-25894-9 [PubMed]
  • 38. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012; 28:882–3. https://doi.org/10.1093/bioinformatics/bts034 [PubMed]
  • 39. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
  • 40. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 41. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
  • 42. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016; 17:218. https://doi.org/10.1186/s13059-016-1070-5 [PubMed]
  • 43. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE II, Koeppen H, Astarita JL, Cubas R, Jhunjhunwala S, Banchereau R, Yang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018; 554:544–8. https://doi.org/10.1038/nature25501 [PubMed]
  • 44. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 45. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, Nordgren H, Farmer P, Praz V, Haibe-Kains B, Desmedt C, Larsimont D, Cardoso F, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006; 98:262–72. https://doi.org/10.1093/jnci/djj052 [PubMed]
  • 46. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, Kamińska B, Huelsken J, Omberg L, Gevaert O, Colaprico A, Czerwińska P, Mazurek S, et al, and Cancer Genome Atlas Research Network. Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell. 2018; 173:338–54.e15. https://doi.org/10.1016/j.cell.2018.03.034 [PubMed]
  • 47. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9:559. https://doi.org/10.1186/1471-2105-9-559 [PubMed]
  • 48. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011; 12:R41. https://doi.org/10.1186/gb-2011-12-4-r41 [PubMed]
  • 49. Huang HY, Lin YC, Cui S, Huang Y, Tang Y, Xu J, Bao J, Li Y, Wen J, Zuo H, Wang W, Li J, Ni J, et al. miRTarBase update 2022: an informative resource for experimentally validated miRNA-target interactions. Nucleic Acids Res. 2022; 50:D222–30. https://doi.org/10.1093/nar/gkab1079 [PubMed]
  • 50. Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, Liu XQ, Sher X, Jung H, Lee M, Lee S, Park SH, Park JO, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018; 24:1449–58. https://doi.org/10.1038/s41591-018-0101-z [PubMed]
  • 51. Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, Zimmer L, Gutzmer R, Satzger I, Loquai C, Grabbe S, Vokes N, Margolis CA, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019; 25:1916–27. https://doi.org/10.1038/s41591-019-0654-5 [PubMed]
  • 52. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, Seja E, Lomeli S, Kong X, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2016; 165:35–44. https://doi.org/10.1016/j.cell.2016.02.065 [PubMed]
  • 53. Cancer Genome Atlas Research Network. Electronic address: [email protected], and Cancer Genome Atlas Research Network. Integrated Genomic Characterization of Pancreatic Ductal Adenocarcinoma. Cancer Cell. 2017; 32:185–203.e13. https://doi.org/10.1016/j.ccell.2017.07.007 [PubMed]
  • 54. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson Å, Kampf C, Sjöstedt E, Asplund A, Olsson I, Edlund K, Lundberg E, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015; 347:1260419. https://doi.org/10.1126/science.1260419 [PubMed]