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) [27]. 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 [30]. 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.