- Primary research
- Open Access
Cross talk between RNA modification writers and tumor development as a basis for guiding personalized therapy of gastric cancer
Human Genomics volume 16, Article number: 14 (2022)
Gastric cancer (GC) shows high metastasis and low survival. RNA modification writers play critical roles in tumor development. This study examined the clinical significance of RNA modification writers in GC prognosis based on four types of adenosine modifications (m1A, m6A, APA and A-to-I).
Writers demonstrated high mutation and expression in GC patients. Different expressions of 26 RNA modification writers were differentially associated with GC prognosis. High-WM score group appeared worse overall survival, higher immune infiltration and activation of EMT pathways than low-WM score group. WM score was correlated with both miRNAs-targeted signaling pathways and patients’ sensitivity to chemotherapeutic drugs and efficacy of immunotherapy.
This study further revealed the close association between adenosine-related RNA modifications and progression of GC. A cross talk between EMT and RNA modification was identified to be one of the mechanisms underlying GC development. Our WM scoring system could serve as a clinical indicator for predicting GC prognosis. Importantly, the WM score could guide personalized treatments such as chemotherapy and immunotherapy for GC patients.
Gastric cancer (GC) shows a high incidence rate each year worldwide. In 2020, 1,089,103 new GC cases were diagnosed, accounting for about 11.1% of all new cancer cases that year . The incidence and mortality demonstrate regional variations. GC has the highest diagnosis rate in Japanese male population (32.5/100,000) and Mongolian female population (13.2/100,000) . The top 5 regions with a high incidence are Eastern Asia, Eastern Europe, South America, Western Asia and South Europe. In China, the morbidity rate of GC reaches 26.54 per 100,000 in male and 11.09 per 100,000 in female . An estimated number of 478,508 new cases were diagnosed in Chinese population in 2020, with an age-standardized incidence rate (ASR) of 20.6/100,000 . Depending on different stages, the 5-year overall survival (OS) of GC is between 10 and 93% . Even if patients have been properly treated with chemotherapies, the median OS of metastatic GC patients is shorter than 1 year . Thus, an effective indicator predictive of GC prognosis or capable of guiding targeted therapy for GC patients is highly needed.
In the past decades, epigenetic modification has been illustrated to be closely associated with cancer development . A number of studies applied the epigenetic modification in clinical diagnosis or treatment of GC. For example, Meng et al. discovered a strong association between DNA methylation and tumor microenvironment (TME) in GC and explored a DNA methylation score (DMS) that may help guide personalized immunotherapy . Jing et al. explored the association between N6-methyladenosine (m6A) enzymes and GC prognosis and found that m6A writers (METTL3, RBM15 and WTAP), m6A erasers (FTO and ALKBH5) and m6A readers (YTHDF3 and YTHDC2) were related to GC development and prognosis . Based on m6A RNA methylation, a prognostic signature composed of FTO, RBM15 and ALKBH5 has been constructed to predict GC overall survival (OS) . A review on m6A modifications summarized that m6A writers, erasers and binding proteins play important roles in gastrointestinal tract cancers . These findings suggested that biomarkers identified based on epigenetic modification could assist clinical decision-making.
Epigenetic modification, especially the most frequent RNA modification, is an essential part of gene regulation. Dysregulation of epigenetic modification can result in aberrant gene expression and dysregulated functional pathways [5, 10, 11]; moreover, dysfunction of methylation signaling pathways will increase or reduce modifications, thereby affecting translation or post-translational processing. It was demonstrated that reduced m6A methylation could activate oncogenic Wnt/PI3K-Akt pathway, thereby promoting malignant transformation of GC cells . Zhang et al. identified three m6A clusters (cluster A, B and C) corresponding to three immune phenotypes in GC . Specifically, cluster A (immune-excluded phenotype), which is also known as T cell-suppressive type, shows high activity of angiogenesis, epithelial–mesenchymal transition (EMT) and TGF-β pathways .
Although m6A modification has been found to be closely involved in cancer progression through altering TME or oncogenic pathways [14,15,16,17,18], other less studied but common RNA modifications such as N1-methyladenosine (m1A), alternative polyadenylation (APA) and adenosine-to-inosine (A-to-I) editing are also critical in cancer development . Zhao et al. first demonstrated that m1A enzymes were dysregulated in gastrointestinal cancer and was involved in tumorigenesis possibly via ErbB and mTOR pathways . APA as an essential step of processing mRNA maturation is present in at least 70% of human genes . Evidence revealed that 95% APA genes have shorter 3′-untranslated regions (3′ UTRs) in cancers, which can facilitate tumor cell proliferation through inhibiting microRNAs (miRNAs)-mediated suppression [22, 23]. CPSF6 involved in APA process has been proven to induce tumor cell proliferation and inhibit apoptosis through APA-mediated 3′UTR shorting in gastric cells . A-to-I RNA editing is one of posttranscriptional modifications that is thought to be altered in thousands of miRNA binding sites in 32 cancer types . In GC, ADAR1 and ADAR2, enzymes of A-to-I editing, are dysregulated and its mis-editing is highly responsible for GC pathogenesis , suggesting that the activity of ADAR1 and ADAR2 may serve as potential biomarkers for GC.
Based on previous findings, this study elucidated the relation between four major RNA modifications (m1A, m6A, APA and A-to-I editing) related to adenosine editing and GC through applying bioinformatics analysis. We developed an RNA modification writers-based signature strongly correlated with GC prognosis, TME and EMT pathway. The signature demonstrated strong potentials in guiding chemotherapy and immunotherapy.
Different gene mutation patterns and expressions of RNA modification writers in cancer cells and normal cells
The workflow of this study was shown (Additional file 1: Fig. S1). The study included four types of RNA modifications writers (writers of m1A, m6A, APA and A-to-I editing) and analyzed a total of 26 genes of writers. In TCGA-STAD dataset, 353 samples were recruited to analyze gene mutations and mutation patterns of the 26 writers, and 118 of them were mutant (Fig. 1A). ZC2H13 was the most frequent mutated gene, showing a mutation frequency of 30%, and there were 10 genes that reached a mutation frequency up to 10% (Fig. 1A). The majority of mutations were missense mutation and frame-shift deletion. Interestingly, patients with mutations showed a longer OS than those without, suggesting a potential relationship between the writers and GC (p = 0.0245, Fig. 1B). GSEA was applied to enrich hallmark gene sets in the mutated samples. Here, we detected a series of tumor-related gene sets, such as E2F targets, Myc targets, G2/M checkpoint, pancreas beta cells and interferon-gamma response (Fig. 1C). In addition, by using PathScore tool, we observed that these significantly mutated genes were highly enriched in BARD1 pathway (56.5% patients, p < 1e−16) and interferon alpha–beta signaling pathway (33% patients, p = 0.048, Additional file 2: Fig. S2).
In addition to small-scale mutations, copy number variation (CNV) is a major factor contributing to aberrant gene expression. Our result showed that ADAR, CPSF1 and CPSF4 had high frequency of increased copy numbers, while RBM15 and RBM15B showed a high proportion of decreased copy numbers; interestingly, no CNV was detected in KIAA1429 (Fig. 1D). Compared with normal samples, GC cancer samples exhibited obviously higher expression of m6A, m1A and APA writers (Fig. 1E), but ADAR only showed significantly higher expression in 3 A-to-I writers (Fig. 1E). Notably, RBM15B with loss of copy numbers still expressed higher in cancer samples. To further understand the relation of CNV and gene expression, the samples were classified into CNV amplification and CNV deletion groups. In both amplification and deletion groups, we observed a higher expression level in most writer genes in GC samples than in normal samples (Fig. 1F). The above results indicated that the expression of writers was involved in a complicated regulation network and may be affected by small- or large-scale mutations. Additionally, these data also suggested that abnormal expression was a potential factor for GC tumorgenesis and development.
Expression patterns of RNA modification writers were associated with GC prognosis and tumor immune cell infiltration
To investigate the relation between each writer and GC prognosis, univariate Cox regression analysis was performed on 300 samples of GSE62254 dataset. Twelve out of 26 writers including ADARB1, PABPN1, PCF11 and METTL3 were positively related to hazard ratio (HR) and significantly related to the cancer prognosis, and NUDT21, CLP1, RBM15B, CSTF2, WTAP, CPSF3, CPSF2 and RBM15 were negatively related to HR (Fig. 2A). To validate their prognostic performance, we implemented immunohistochemistry (IHC) of two writers (ADARB1 and RBM15B) on GC tissues. Consistent with the above analysis, ADARB1 was higher expressed and RBM15B was lower expressed in advanced gastric cancer tissue (Additional file 3: Fig. S3), but both the two showed opposite expressions in early GC tissue, demonstrating the reliability of these writers of serving as prognostic indicators. Spearman rank analysis on the correlations among writers revealed interactions among different writers (Fig. 2B).
As the specific roles of writers are different in the development of GC prognosis and interactions among writers, we speculated whether different expression patterns of all writers would lead to different prognostic results. Thus, unsupervised consensus clustering was performed to classify 300 samples into two groups based on different expression patterns of 26 writers (Fig. 2C). Survival analysis showed a significant difference between the two groups, and cluster 1 group tended to develop a more favorable OS than cluster 2 (p = 0.002, Fig. 2D). Then, we used GSEA to analyze the correlation between the two clusters and functional pathways. Cluster 1 group was enriched to pathways related to DNA repair and cell cycle, such as base excision repair pathway, homologous recombination pathway, mismatch repair pathway, DNA replication pathway and cell cycle pathway (Fig. 2E). Cluster 2 group presented a lower GSVA score in the above pathways, but many metabolism pathways, for example, arachidonic acid metabolism pathway, and metabolism of xenobiotics were enriched by cytochrome P450 pathway and drug metabolism cytochrome pathway (Fig. 2E). On biological process (BP) and molecular function (MF) terms, the two clusters also exhibited apparently distinct enrichment (Additional file 4: Fig. S4). A total of 1046 (21.0%) and 1770 (35.5%) BP terms were significantly annotated in clusters 1 and 2, respectively (FDR < 0.05). Metabolism and RNA modification-related BP terms were highly enriched in Cluster 1, while Cluster 2 had a high correlation with cell proliferation, migration and autophagy (Additional file 4: Fig. S4A). Moreover, 328 and 251 MF terms were significantly enriched in clusters 1 and 2, respectively (FDR < 0.05); specifically, cluster 1 was more correlated with RNA synthesis, and cluster 2 was more associated with protein binding and extracellular matrix (Additional file 4: Fig. S4B). The results further supported the conclusion that Cluster 2 had worse prognosis, which may be resulted from its high activity of tumor cell proliferation and migration.
Previous studies showed that RNA modifications are correlated with tumor immune infiltration; therefore, we performed CIBERSORT to calculate the enrichment score of immune-related cells and analyzed its correlation with writers. Of the four types of writers, the expression of m6A and APA writers, especially gamma delta T cells, activated memory CD4 T cells, M0 macrophages, activated dendritic cells, resting and activated mast cells, were found to be closely associated with immune infiltration (Additional file 5: Fig. S5). Activated mast cells and M1 macrophages were greatly enriched in cluster 1, and gamma delta T cells, M2 macrophages and resting memory CD4 T cells were noticeably enriched in cluster 2 (Fig. 2F). Through using Timer method, we also observed similar results that cluster 2 had higher immune infiltration than cluster 1 (Additional file 6: Fig. S6).
Construction of a prognostic signature based on genes related to RNA modification writers
Variations in expression patterns of RNA modification writers normally result in different influence on GC prognosis; we therefore developed a signature for predicting the prognosis based on the genes related to RNA modification writers. A total of 194 differentially expressed genes related to RNA modification were identified through univariate Cox regression analysis. Gene ontology (GO) analysis revealed that these 194 genes were involved in multiple biological functions (Additional file 7: Fig. S7A–C). Next, according to unsupervised clustering analysis, 300 samples of GSE62254 dataset were stratified into two groups (cluster A and cluster B) (Additional file 9: Table S1). Survival analysis exhibited significantly differential OS of the two groups, with higher OS in cluster B (p = 0.0032, Additional file 7: Fig. S7D). Two types of systems for classifying GC patients were constructed based on 26 writers and 194 genes related to RNA modification writers. The performance of the two systems in clustering 300 samples were similar (Fig. 3A), indicating the reliability of developing a prognostic signature using RNA modification writers.
Based on 194 genes, we established a WM (Writers of RNA Modification) scoring system, which was shown as WM score = Σ(coefficient i) * (expression i) . The scoring system showed robust performance in different groups, in which cluster 2 had higher WM score than cluster 1 (p = 3.40e−29, Fig. 3B). Consistently, cluster B also had higher WM score than cluster A (p = 9.70e−49, Fig. 3C). Samples in GSE62254 dataset were classified into high-risk and low-risk groups with distinct OS using the WM score (p = 0.0015, Fig. 3D). According to receiver operator curve (ROC), WM score showed an area under ROC curve (AUC) of 0.64, 0.71 and 0.71, respectively, in predicting 1-year, 3-year and 5-year OS (Fig. 3E). Similarly, 353 samples in TCGA-STAD dataset were divided into high-risk and low-risk groups with differential prognosis (p = 0.045, Fig. 3F), and ROC analysis also demonstrated a robust performance of the WM score (Fig. 3G). Furthermore, multivariate Cox regression analysis in both datasets verified that WM score could serve as an effective indicator for GC prognosis prediction, with a hazard ratio of 1.5 (95% CI = 1.3–1.7, p < 0.0001) in TCGA-STAD dataset and 1.2 (95% CI = 1.0–1.4, p = 0.0098) in GSE62254 dataset (Fig. 3H, I).
WM score was correlated with tumor microenvironment
To examine whether WM score was related to tumor microenvironment (TME), CIBERSORT was applied to calculate the proportion and enrichment score of 22 immune cells. In GSE62254 dataset, the proportion of 22 immune cells in high-risk and low-risk groups was shown in the heatmap and histogram (Fig. 4A–C). Among the 22 immune cells, gamma delta T cells, M1 macrophages and activated mast cells consisted of a large proportion (Fig. 4B), and 11 immune cells were differentially enriched in the two groups (Fig. 4C). Particularly, activated memory T cells, M0 macrophages, activated dendritic cells and activated mast cells were more enriched in the low-risk group, while gamma delta T cells, M2 macrophages, resting dendritic cells and eosinophils were more enriched in the high-risk group (p < 0.05, Fig. 4C). After excluding seven immune cells without expression from more than half of the samples, 15 immune cells remained and their distributions in the two groups are shown in Fig. 4D. Immune infiltration analysis using CIBERSORT showed a higher enrichment score of stromal score, immune score and ESTIMATE score in high-risk group than low-risk group, indicating that WM score was positively related to immune infiltration (Fig. 4E, F).
The relation between WM score and clinical features, epithelial–mesenchymal transition
Next, we identified the relation between WM score and clinical features including T stage, N stage, M stage, stage I to IV, gender and age. In T stage, T3 had the highest WM score, and the score of T2 was significantly different from that of T3 and T4 (Fig. 5A). Though the difference of WM score was not obviously different in T stage or M stage, there was still a tendency of the higher WM score correlating with the progression of stage (Fig. 5B, C). Moreover, such a tendency was also observed in stage I to IV, which manifested significant differences between stage I and III, stage II and III, stage II and IV (Fig. 5D). Notably, the WM score of female group was higher than male group, and age was another important factor leading to differential WM score (Fig. 5E, F).
As epithelial–mesenchymal transition plays a crucial role in tumor progression, we also assessed the enrichment of EMT-related pathways in high-risk and low-risk groups using GSVA. In GSE62254 dataset, EMT1, EMT2, Pan-F-TBRS, EMT3 and angiogenesis pathways were highly enriched to high-risk group, while DNA replication, DNA repair and cell cycle pathways were more enriched to low-risk group (Fig. 5G, H). The majority of pathways showed significant differences between the two groups. In addition, similar result was also observed in TCGA-STAD dataset (Fig, 5I, J), indicating that the WM score was strongly associated with EMT. Due to the important role of EMT in cancer development, a number of studies have explored a series of EMT-related signatures [28,29,30]; for instance, Wang et al. built an EMT-CYT Index (ECI) model as an indicator for pan-cancer prognosis . We also compared our WM score with EMT score developed by Wang et al. according to hazard ratio, C-index and AUC. The EMT score was calculated for each sample in TCGA-STAD dataset. Although the EMT signature showed a more favorable result in hazard ratio, our WM score had higher C-index and AUC, which supported its robustness (Additional file 8: Fig. S8).
WM score was related to tumor mutation burden
Tumor mutation burden (TMB) has been considered to be correlated with tumor immune infiltration, and this can provide a guidance to immunotherapy to some extent. The group with low WM score had higher TMB than that with high WM score in TCGA-STAD dataset (Fig. 6A), and a negative correlation of R = − 0.424 (p = 9.33e−17, Fig. 6B) between WM score and TMB was observed. The top 20 significantly mutated genes were listed with their p value and proportion (Fig. 6C). Missense mutations accounted for the greatest proportion among these mutation types. Frame-shift deletions frequently occurred in TTN, TP53 and KMT2D, and nonsense mutations were found the most in TTN and TP53. Tumor-related genes, such as TP53, LRP1B and MUC16, showed great mutations in both groups but apparently more in low-risk group.
The relation between WM score and micro-RNA regulation
Previous studies demonstrated that 3′-UTR alterations caused by APA will affect posttranscriptional regulation and translation, which could result in loss of miRNA binding sites [31, 32]. The WM score system was developed based on 26 types of RNA modification writers, and we speculated that there was a correlation between WM score and miRNA regulation. Firstly, 46 differentially expressed miRNAs were screened between high-WM score and low-WM score groups. Then, miRNA-targeted genes of these miRNAs were enriched in KEGG pathways using GSVA. In the top 10 enriched pathways, the target genes of miRNAs were highly enriched in PI3K-Akt signaling pathway, focal adhesion and MAPK signaling pathway (Fig. 7A). Moreover, miRNA-targeted genes were differentially expressed between the two groups in these signaling pathways (Fig. 7B), suggesting that WM score was related to the regulation of miRNA-related signaling pathways.
The sensitivity of chemotherapeutic drugs was related to WM score
To examine whether WM score and patients’ sensitivity of chemotherapeutic drugs was correlated, we exported a series of drugs for treating GC from GDSC database. Spearman rank correlation analysis identified 12 chemotherapeutic drugs significantly correlated with WM score. Among these drugs, WM score showed close relations to the sensitivity of two drugs (p38 MAPK inhibitor doramapimod and IGF-1R inhibitor BMS-754807) and to the resistance of 10 drugs (Fig. 8A). Furthermore, pRRophetic R package was used to evaluate the drug response of TCGA-STAD samples to cisplatin, paclitaxel, docetaxel and 5-fluorouracil. Low score of estimated IC50 represented a high sensitivity to drugs. Low-WM score group showed high sensitivity to cisplatin and docetaxel, and patients with a high-WM score were more sensitive to paclitaxel and 5-fluorouracil (Fig. 8B–E). The targeted pathways of chemotherapeutic drugs were also investigated. Here, each drug targeted different pathways, such as IGF-1R inhibitor, and BMS-754807 and oxaliplatin were related to IR signaling pathway and DNA replication pathway, respectively (Fig. 8F).
The efficacy of anti-PD-1 therapy is associated with WM score
For metastatic cancer, immunotherapy, especially anti-PD-1/PD-L1 therapy, is a promising strategy in treating many cancer types. We therefore examined the WM score in predicting the efficacy of immunotherapy. In GSE78220 dataset, 27 patients with four outcome information of complete response (CR), partial response (PR), stable disease (SD) and progressive disease (PD) treated by anti-PD-1 therapy were included. The result presented that CR and PR patients had a lower WM score than SD and PD patients (p < 0.01, Fig. 9A), and a higher proportion of CR and PR patients was in low-WM score group (Fig. 9B). Kaplan–Meier survival plot also displayed a favorable OS in low-WM score group (p = 0.021, Fig. 9C), suggesting that WM score was closely associated with the outcomes of immunotherapy.
Writers, readers and erasers are the three different types of RNA modifications, and their effect on stimulating or suppressing tumor growth is regarded as hallmarks of cancers . Studies have proposed various signatures of RNA modification writers such as METTL3 , NSUN2  and a signature of 3 m6A regulators (FTO, RBM15, ALKBH5)  to predict GC prognosis. This study focused on a series of writers including m1A, m6A, APA and A-to-I writers capable of modifying adenosine. A total of 26 writers were included in the present study, and the assessment of multiple aspects on the writers confirmed their clinical significance in GC survival and treatment.
Compared with normal samples, high expression of 24 writers in GC samples was detected (Fig. 1E). Mutant group and WT group showed differential OS (Fig. 1B), suggesting that these mutations on RNA writers possibly had an important effect on their expression level or function, which also contributed to tumor development. Exploration on the GEO database showed that many writers differentially expressed in gastric cancer were also differentially expressed in other cancer types. For example, METTL3 is significantly up-regulated in esophageal squamous cell carcinoma, and this leads to a worse prognosis through activating Notch signaling pathway . Wang et al. discovered that the inhibition of METTL3 or METTL14 can help colorectal cancer and melanoma cell acquire an enhanced immune response to anti-PD-1 therapy through increasing CD8 T cells and IFN-γ secretion . Such an observation revealed that METTL3 or METTL14 could serve as potential targets for assisting immune checkpoint blockade therapy. WTAP is a nuclear protein that facilitates cell proliferation and apoptosis. WTAP expression is elevated in cholangiocarcinoma, and WTAP knockdown significantly suppresses tumorigenesis in a xenograft model . KIAA1429, the largest known component in the m6A methyltransferase complex, is overexpressed in hepatocellular carcinoma tissues and is associated with poor prognosis , which is consistent in GC.
Based on the discovery that different expression patterns of RNA modification writers were related to different OS, we constructed a WM scoring system with high-performance prognosis prediction. Previous researches have illustrated that RNA modifications can regulate immune functions such as immune recognition and immune responses [39, 40]. TME is a decisive factor of tumor development and could affect the efficacy of targeted treatment such as immunotherapies. According to the infiltration of cytotoxic lymphocytes (CTLs) and tumor-associated macrophages (TAMs) in tumor cells, TME can be divided into immune-inflammation, immune-exclusion and immune-desert . Immune-inflammation type shows abundant expressions of PD-L1 and activated CTLs; therefore, patients of this type can benefit the most from the anti PD-1/PD-L1 therapy . This study presented that high-WM score group had a higher enrichment score of immune cells, indicating that patients in high-WM score group could benefit more from the anti-PD-1/PD-L1 blockade. Moreover, as high- and low-WM groups demonstrated significantly different OS after receiving anti-PD-1/PD-L1 therapy, the WM scoring system may have a potential in predicting patients who may develop favorable outcome from the targeted therapy. Although no treatment data of GC patients with anti-PD-1/PD-L1 therapy could be obtained, the results still provided a possibility that WM score can guide immunotherapy.
EMT plays an essential role in embryonic development and allows epithelial cells to transform into mesenchymal cells . EMT also participates in tumor progression, metastasis, recurrence and chemoresistance . EMT process varies greatly in different cancer types, which is possibly related to different EMT-related pathways activated by cancers. Notably, cytokines, chemokines and growth factors secreted by stromal cells and immune cells can induce EMT process, showing a close relation between tumor microenvironment and EMT [44, 45]. A number of studies developed various signatures related to EMT in pan-cancer or different cancer types [29, 46, 47]. For instance, Dai et al. established an EMT-related signature for predicting GC prognosis, and five EMT-related prognostic genes were found to be highly expressed in GC tissues . Wang et al. proposed an EMT-CYT Index (ECI) model capable of predicting prognosis for pan-cancer patients . Our WM score manifested a better performance in prognosis prediction when compared with the EMT signature of Wang et al. Notably, EMT has a correlation with ECM remodeling , and this also supports our result that ECM–receptor interaction pathway was highly enriched in differentially expressed miRNAs between high-WM and low-WM score groups.
Although various EMT-related signatures were developed for many cancer types in previous studies, our WM scoring system further validated the importance of EMT in cancer progression and discovered the cross talk between RNA modification and EMT. Consistent with the previous findings, in this study, EMT pathways were highly activated in the high-WM score group, which also had a high immune enrichment score. The result indirectly proved the close relation between TME and EMT. In accordance with the previous observation that EMT process can promote tumor metastasis, the current research confirmed the relation between worse prognosis and highly activated EMT-related pathways in the high-WM score group.
RNA modification writers play a critical role in regulating biological process in cancer. Writers such as METTL3 and METTL14 show cancer-promoting or cancer-suppressive effects on the same or different cancer types [50,51,52,53]. To some extent, this may interfere the development of WM scoring system. Another limitation of the study was that the anti-PD-1/PD-L1 treatment data came from melanoma, as no available anti-PD-1/PD-L1 treatment data can be obtained from GC.
Nevertheless, this study still expanded the research on the effect of adenosine-related RNA modifications on GC development and pathogenesis. These RNA modifications were closely involved in altering TME and TME through oncogenic pathways such as proteoglycans, PI3K-Akt signaling, MAPK signaling, and ECM-receptor interaction pathways. High-WM and low-WM score groups manifested significant differences in GC patients’ sensitivity to chemotherapy and immunotherapy. Hence, the WM scoring system in this study showed a great potential in predicting GC prognosis and could assist the decision-making of targeted therapies.
In conclusion, this study constructed a WM scoring system based on four types of RNA modification writers. The WM scoring system exhibited robust performance in all the datasets explored. Overall, high-WM and low-WM score groups classified by the WM scoring system manifested distinct prognostic difference, immune infiltration, EMT-related pathways, gene mutation and sensitivity to chemotherapy and immunotherapy. We concluded that adenosine-related RNA modifications were closely associated with tumorigenesis and pathology in GC possibly through modulating tumor microenvironment and facilitating EMT. Notably, the WM scoring system had the ability to assist and guide chemotherapy and immunotherapy, especially for metastatic GC patients.
GC samples were exported from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas database (TCGA, https://portal.gdc.cancer.gov/). GSE62254 dataset containing gene expression data and clinical information were obtained from GEO (Additional file 9: Table S1). TCGA-STAD dataset with mRNA expression, miRNA expression, somatic cell mutations, somatic copy number amplifications (SCNAs) and clinical features including stages, sex, age and survival data was downloaded from TCGA (Additional file 10: Table S2). GSE78220 dataset acquired from GEO contained the treatment data of anti-PD-1 therapy for melanoma.
Unsupervised consensus clustering
Unsupervised consensus clustering is a widely used method for providing quantitative and visual stability evidence for estimating the number of unsupervised classes in a dataset . We used ConsensusClusterPlus R package to conduct unsupervised consensus clustering on 26 types of writers in GSE62254 dataset . Thousand repeats were performed to produce stable clusters. Writers included 7 m6A enzymes (KIAA1429, METTL3, METTL14, RBM15, RBM15B, WTAP and ZC3H13), 4 m1A enzymes (TRMT6, TRMT10C, TRMT61A and TRMT61B), 12 APA enzymes (CFI, CLP1, CPSF1-4, CSTF1/2/3, NUDT21, PABPN1 and PCF11) and 3 A-to-I enzymes (ADAR, ADARB1 and ADARB2).
Functional analysis on clusters 1 and 2
Gene set variation analysis (GSVA) was performed to assess the enrichment of mutant and wild-type RNA modification genes in “h.all.v7.2” and the enriched KEGG pathways of clusters 1 and 2 in “c2.cp.kegg.v7.1” using GSVA R package . PathScore  was employed to evaluate the enrichment of molecular function and biological process of RNA modification genes in “c5.bp.v7.1”, “c5.mf.v7.1.” Reference gene sets of “c2.cp.kegg.v7.1,” “c5.bp.v7.1,” “c5.mf.v7.1” and “h.all.v7.2” were downloaded from MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb). ClusterProfiler R package was used to annotate gene ontology terms on differentially expressed genes between clusters 1 and 2 . Benjamini and Hochberg (BH) method was used to correct false discovery rate (FDR) in function annotations, with FDR < 0.05 as a significance.
Immunohistochemistry (IHC) for evaluating RNA writer expression in GC tissues
IHC could identify cellular or subcellular distribution and expression of one objective protein based on monoclonal or polyclonal antibodies. This study performed IHC to verify the expression of ADARB1 (Affbiotech, Catalog number: DF3227, dilution: 1:100) and RBM15 (Proteintech, Catalog number: 10587-1-AP, dilution: 1:200) in early and advanced GC tissues following the general procedures. Twenty pairs of GC and paired normal tissue samples were obtained from The Second Affiliated Hospital of Guangzhou Medical University. Written informed consent was obtained from all patients. The experimental procedures were approved by the Institutional Ethics Committee of The Second Affiliated Hospital of Guangzhou Medical University.
Analysis of immune infiltration
CIBERSORT (http://cibersort.stanford.edu/) was used to analyze the distribution and enrichment score of the 22 immune cells . The CIBERSORT tool could estimate the fraction of 22 immune cells in complex tissues according to gene expression profiles, and it has been applied in a number of studies for characterizing tumor microenvironment in cancers.
Construction of a WM scoring system
Firstly, differentially expressed genes (DEGs) between cluster 1 and cluster 2 in GSE62254 dataset were identified using Limma R package . The expressions of 194 DEGs were used as input to conduct unsupervised consensus clustering in ConsensusClusterPlus R package . Cluster number k = 2 was selected to construct a consensus matrix. A WM scoring system based on Writers of RNA Modifications was presented as WM score = Σ(coefficient i) * (expression i); here, i represented differentially expressed genes significantly related to prognosis . High- and low-WM score groups were classified by the median of WM score. The effectiveness of the WM scoring system was validated by ROC through timeROC package and multivariate Cox regression analysis in survival R package . TCGA-STAD was used as a validation dataset to assess the robustness of the WM scoring system.
Gene mutation analysis
To characterize the difference of gene mutations between high-WM and low-WM groups, mutect2 software (https://software.broadinstitute.org/cancer/cga/mutect) was employed to analyze the gene mutations in TCGA-STAD dataset. Mutect2 could preprocess high-throughput sequencing data and detect mutations with a very low false positive rate. CancereffectsizeR package  was used to calculate mutation frequency, and the top 20 significantly mutated genes were screened by Fisher’s exact test.
Spearman rank correlation analysis
Spearman rank correlation analysis was applied to analyze the correlation among 26 writers, WM score and the sensitivity to chemotherapeutic drugs. The data of drug sensitivity were downloaded from Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database . |rs| > 0.2 was considered to show a correlation, where rs = Spearman rank correlation coefficient. FDR was corrected by BH correction (FDR < 0.05 was significant). Drug response was evaluated by pRRophetic R package .
Analysis of miRNA expression, miRNA-targeted genes and pathways
The miRNA expression data were exported from TCGA-STAD dataset. Differential expressed miRNAs between high- and low- WM score groups were screened by Limma R package. Then, GSVA R package was applied to analyze miRNA-targeted genes and KEGG functional pathways in “c2.cp.kegg.v7.1.” Wilcoxon test was performed to assess differences between the two groups. FDR was corrected by BH correction, and FDR < 0.05 was considered as a significance.
All the statistical analyses were performed in the R (3.6.2). Parameters of measurements were default if there was no special description. Statistical methods were shown in the corresponding figure legends. p < 0.05 was considered to be significant.
Availability of data and materials
The datasets generated during and/or analyzed during the current study are available in the [GSE62254] repository [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62254] and in [GSE78220] repository [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi].
Acute myeloid leukemia
An age-standardized incidence rate
Area under ROC curve
Benjamini and Hochberg
Copy number variation
Differentially expressed genes
A DNA methylation score
False discovery rate
Genomics of Drug Sensitivity in Cancer
Gene Expression Omnibus
Get set variation analysis
The biochemical half maximal inhibitory concentration
Kyoto Encyclopedia of Genes and Genomes
The Molecular Signatures Database
Programmed cell death protein 1
Programmed cell death ligand 1
Receiver operating characteristic curve
Somatic copy number amplifications
The Cancer Genome Atlas
Tumor mutation burden
Writers of RNA Modifications
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality Worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49.
Wang SM, Zheng RS, Zhang SW, Zeng HM, Chen R, Sun KX, et al. Epidemiological characteristics of gastric cancer in China, 2015. Zhonghua liu xing bing xue za zhi = Zhonghua liuxingbingxue zazhi. 2019;40(12):1517–21.
Padmanabhan N, Ushijima T, Tan P. How to stomach an epigenetic insult: the gastric cancer epigenome. Nat Rev Gastroenterol Hepatol. 2017;14(8):467–78.
Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, Lordick F. Gastric cancer. Lancet. 2020;396(10251):635–48.
Yang X, Liu M, Li M, Zhang S, Hiju H, Sun J, et al. Epigenetic modulations of noncoding RNA: a novel dimension of cancer biology. Mol Cancer. 2020;19(1):64.
Meng Q, Lu YX, Ruan DY, Yu K, Chen YX, Xiao M, et al. DNA methylation regulator-mediated modification patterns and tumor microenvironment characterization in gastric cancer. Mol Ther Nucleic Acids. 2021;24:695–710.
Jing JJ, Zhao X, Li H, Sun LP, Yuan Y. Expression profiles and prognostic roles of m6A writers, erasers and readers in gastric cancer. Future Oncol. 2021;17(20):2605–20.
Su Y, Huang J, Hu J. m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gastric cancer. Front Oncol. 2019;9:1038.
Hu BB, Wang XY, Gu XY, Zou C, Gao ZJ, Zhang H, et al. N(6)-methyladenosine (m(6)A) RNA modification in gastrointestinal tract cancers: roles, mechanisms, and applications. Mol Cancer. 2019;18(1):178.
Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell. 2017;169(7):1187–200.
Zhao W, Qi X, Liu L, Ma S, Liu J, Wu J. Epigenetic regulation of m(6)A modifications in human cancer. Mol Ther Nucleic Acids. 2020;19:405–12.
Zhang C, Zhang M, Ge S, Huang W, Lin X, Gao J, et al. Reduced m6A modification predicts malignant phenotypes and augmented Wnt/PI3K-Akt signaling in gastric cancer. Cancer Med. 2019;8(10):4766–81.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53.
Wang X, Huang J, Zou T, Yin P. Human m(6)A writers: two subunits, 2 roles. RNA Biol. 2017;14(3):300–4.
Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019;112:108613.
Chen XY, Zhang J, Zhu JS. The role of m(6)A RNA methylation in human cancer. Mol Cancer. 2019;18(1):103.
Zhou H, Yin K, Zhang Y, Tian J, Wang S. The RNA m6A writer METTL14 in cancers: roles, structures, and applications. Biochim Biophys Acta Rev Cancer. 2021;1876(2):188609.
Zeng C, Huang W, Li Y, Weng H. Roles of METTL3 in cancer: mechanisms and therapeutic targeting. J Hematol Oncol. 2020;13(1):117.
Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer. 2020;20(6):303–22.
Zhao Y, Zhao Q, Kaboli PJ, Shen J, Li M, Wu X, et al. m1A regulated genes modulate PI3K/AKT/mTOR and ErbB pathways in gastrointestinal cancer. Transl Oncol. 2019;12(10):1323–33.
Tian B, Manley JL. Alternative polyadenylation of mRNA precursors. Nat Rev Mol Cell Biol. 2017;18(1):18–30.
Mayr C, Bartel DP. Widespread shortening of 3′UTRs by alternative cleavage and polyadenylation activates oncogenes in cancer cells. Cell. 2009;138(4):673–84.
Xia Z, Donehower LA, Cooper TA, Neilson JR, Wheeler DA, Wagner EJ, et al. Dynamic analyses of alternative polyadenylation from RNA-seq reveal a 3′-UTR landscape across seven tumour types. Nat Commun. 2014;5:5274.
Shi X, Ding K, Zhao Q, Li P, Kang Y, Tan S, et al. Suppression of CPSF6 enhances apoptosis through alternative polyadenylation-mediated shortening of the VHL 3′UTR in gastric cancer cells. Front Genet. 2021;12:707644.
Pinto Y, Buchumenski I, Levanon EY, Eisenberg E. Human cancer tissues exhibit reduced A-to-I editing of miRNAs coupled with elevated editing of their targets. Nucleic Acids Res. 2018;46(1):71–82.
Chan TH, Qamra A, Tan KT, Guo J, Yang H, Qi L, et al. ADAR-mediated RNA editing predicts progression and prognosis of gastric cancer. Gastroenterology. 2016;151(4):637-50.e10.
Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, et al. Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer. 2021;20(1):29.
Chen T, Zhao Z, Chen B, Wang Y, Yang F, Wang C, et al. An individualized transcriptional signature to predict the epithelial–mesenchymal transition based on relative expression ordering. Aging. 2020;12(13):13172–86.
Mak MP, Tong P, Diao L, Cardnell RJ, Gibbons DL, William WN, et al. A patient-derived, pan-cancer EMT signature identifies global molecular alterations and immune target enrichment following epithelial-to-mesenchymal transition. Clin Cancer Res. 2016;22(3):609–20.
Wang G, Xu D, Zhang Z, Li X, Shi J, Sun J, et al. The pan-cancer landscape of crosstalk between epithelial–mesenchymal transition and immune evasion relevant to prognosis and immunotherapy response. NPJ Precis Oncol. 2021;5(1):56.
Turner RE, Pattison AD, Beilharz TH. Alternative polyadenylation in the regulation and dysregulation of gene expression. Semin Cell Dev Biol. 2018;75:61–9.
Fu Y, Chen L, Chen C, Ge Y, Kang M, Song Z, et al. Crosstalk between alternative polyadenylation and miRNAs in the regulation of protein translational efficiency. Genome Res. 2018;28(11):1656–63.
Wang Q, Chen C, Ding Q, Zhao Y, Wang Z, Chen J, et al. METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut. 2020;69(7):1193–205.
Mei L, Shen C, Miao R, Wang JZ, Cao MD, Zhang YS, et al. RNA methyltransferase NSUN2 promotes gastric cancer cell proliferation by repressing p57(Kip2) by an m(5)C-dependent manner. Cell Death Dis. 2020;11(4):270.
Han H, Yang C, Zhang S, Cheng M, Guo S, Zhu Y, et al. METTL3-mediated m(6)A mRNA modification promotes esophageal cancer initiation and progression via Notch signaling pathway. Mol Ther Nucleic Acids. 2021;26:333–46.
Wang L, Hui H, Agrawal K, Kang Y, Li N, Tang R, et al. m(6) A RNA methyltransferases METTL3/14 regulate immune responses to anti-PD-1 therapy. EMBO J. 2020;39(20):e104514.
Jo HJ, Shim HE, Han ME, Kim HJ, Kim KS, Baek S, et al. WTAP regulates migration and invasion of cholangiocarcinoma cells. J Gastroenterol. 2013;48(11):1271–82.
Lan T, Li H, Zhang D, Xu L, Liu H, Hao X, et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer. 2019;18(1):186.
Shulman Z, Stern-Ginossar N. The RNA modification N(6)-methyladenosine as a novel regulator of the immune system. Nat Immunol. 2020;21(5):501–12.
Eisenberg E, Levanon EY. A-to-I RNA editing—immune protector and transcriptome diversifier. Nat Rev Genet. 2018;19(8):473–90.
Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018;24(5):541–50.
Teng MW, Ngiow SF, Ribas A, Smyth MJ. Classifying cancers based on T-cell infiltration and PD-L1. Can Res. 2015;75(11):2139–45.
Yang J, Weinberg RA. Epithelial–mesenchymal transition: at the crossroads of development and tumor metastasis. Dev Cell. 2008;14(6):818–29.
Iwatsuki M, Mimori K, Yokobori T, Ishi H, Beppu T, Nakamori S, et al. Epithelial–mesenchymal transition in cancer development and its clinical significance. Cancer Sci. 2010;101(2):293–9.
Dongre A, Weinberg RA. New insights into the mechanisms of epithelial–mesenchymal transition and implications for cancer. Nat Rev Mol Cell Biol. 2019;20(2):69–84.
Gibbons DL, Creighton CJ. Pan-cancer survey of epithelial–mesenchymal transition markers across the Cancer Genome Atlas. Dev Dyn. 2018;247(3):555–64.
Rokavec M, Kaller M, Horst D, Hermeking H. Pan-cancer EMT-signature identifies RBM47 down-regulation during colorectal cancer progression. Sci Rep. 2017;7(1):4687.
Dai W, Xiao Y, Tang W, Li J, Hong L, Zhang J, et al. Identification of an EMT-related gene signature for predicting overall survival in gastric cancer. Front Genet. 2021;12:661306.
Peixoto P, Etcheverry A, Aubry M, Missey A, Lachat C, Perrard J, et al. EMT is associated with an epigenetic signature of ECM remodeling genes. Cell Death Dis. 2019;10(3):205.
Ma JZ, Yang F, Zhou CC, Liu F, Yuan JH, Wang F, et al. METTL14 suppresses the metastatic potential of hepatocellular carcinoma by modulating N(6) -methyladenosine-dependent primary MicroRNA processing. Hepatology. 2017;65(2):529–43.
Lin S, Choe J, Du P, Triboulet R, Gregory RI. The m(6)A methyltransferase METTL3 promotes translation in human cancer cells. Mol Cell. 2016;62(3):335–45.
Lin X, Chai G, Wu Y, Li J, Chen F, Liu J, et al. RNA m(6)A methylation regulates the epithelial mesenchymal transition of cancer cells and translation of Snail. Nat Commun. 2019;10(1):2065.
Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, et al. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology. 2018;67(6):2254–70.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7.
Gaffney SG, Townsend JP. PathScore: a web tool for identifying altered pathways in cancer data. Bioinformatics. 2016;32(23):3688–90.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32(30):5381–97.
Cannataro VL, Gaffney SG, Townsend JP. Effect sizes of somatic mutations in cancer. J Natl Cancer Inst. 2018;110(11):1171–7.
Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic acids Res. 2013;41(Database issue):D955–61.
Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9(9):e107468.
This study was supported by [Guangzhou Medical University Dr's Start-up Project] under Grant Number [2016C18].
Ethics approval and consent to participate
The experimental procedures were approved by the Institutional Ethics Committee of The Second Affiliated Hospital of Guangzhou Medical University.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Enriched pathways of significantly mutated genes evaluated by PathScore.
IHC results of ADARB1 and RBM15B in gastric cancer tissues. (A-B) IHC of ADARB1 in early gastric cancer tissue. (C–D) IHC of ADARB1 in advanced gastric cancer tissue. (E–F) IHC of RBM15B in early gastric cancer tissue. (G–H) IHC of RBM15B in advanced gastric cancer tissue. (I–J) Sections were semi-quantitatively scored for ADARB1 or RBM15B staining patterns as follows: the staining extent in each core was scored as 1+ (< 25% staining of tumor cells), 2+ (25–50% staining of tumor cells), 3+ (50% to 75% staining of tumor cells), or 4+ (> 75% staining of tumor cells). Additionally, the staining intensity was quantified as 0 (negative), 1+ (weak), 2+ (intermediate), or 3+ (strong). The final immunoreaction score was obtained by multiplying the intensity and extension values (range 0–12) and the samples were grouped as 1+ (score 0–2), 2+ (score 3–4), 3+ (score 6–8) and 4+ (score 9–12) staining. Meanwhile, for statistical purposes, scores of 3+ and 4+ was defined as high expression and the other final scores were considered as low expression, and then chi-squared test was used to compare the differences between high-expression and low-expression groups.
The top 10 significantly enriched BP terms (A) and MF (B) terms of clusters 1 and 2. FDR < 0.05.
The correlation between 26 RNA modification writers and immune infiltration. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Evaluation of immune infiltration of two clusters by Timer. ***p < 0.001.
Functional analysis of 194 RNA modification related genes and survival analysis of 300 samples in GSE62254 dataset. (A-C) GO analysis of 194 genes on molecular function (A), cellular component (B) and biological process (C). (D) Kaplan-Meier survival curve of cluster A and cluster B in GSE62254 dataset. Log-rank test was performed.
Comparison of WM score and EMT score on hazard ratio, C-index, and 1-year, 3-year and 5-year AUC. CI, confidence interval. AUC, area under ROC curve.
Comparison of WM score and EMT score on hazard ratio, C-index, and 1-year, 3-year and 5-year AUC. CI, confidence interval. AUC, area under ROC curve.
Clinical information of TCGA-STAD dataset.
Clinical information of GSE62254 dataset.
About this article
Cite this article
Zhang, S., Kuang, G., Huang, Y. et al. Cross talk between RNA modification writers and tumor development as a basis for guiding personalized therapy of gastric cancer. Hum Genomics 16, 14 (2022). https://doi.org/10.1186/s40246-022-00386-z
- Gastric cancer
- RNA modifications
- Epithelial–mesenchymal transition (EMT)
- Tumor microenvironment (TME)
- Bioinformatics analysis