Construction of an RNA modification-related gene predictive model associated with prognosis and immunity in gastric cancer
BMC Bioinformatics volume 24, Article number: 147 (2023)
Gastric cancer (GC) is one of the most common causes of cancer-related fatalities worldwide, and its progression is associated with RNA modifications. Here, using RNA modification-related genes (RNAMRGs), we aimed to construct a prognostic model for patients with GC.
Based on RNAMRGs, RNA modification scores (RNAMSs) were obtained for GC samples from The Cancer Genome Atlas and were divided into high- and low-RNAMS groups. Differential analysis and weighted correlation network analysis were performed for the differential expressed genes (DEGs) to obtain the key genes. Next, univariate Cox regression, least absolute shrinkage and selection operator, and multivariate Cox regression analyses were performed to obtain the model. According to the model risk score, samples were divided into high- and low-risk groups. Enrichment analysis and immunoassays were performed for the DEGs in these groups. Four external datasets from Gene Expression Omnibus data base were used to test the accuracy of the predictive model.
We identified SELP and CST2 as key DEGs, which were used to generate the predictive model. The high-risk group had a worse prognosis compared to the low-risk group (p < 0.05). Enrichment analysis and immunoassays revealed that 144 DEGs related to immune cell infiltration were associated with the Wnt signaling pathway and included hub genes such as ELN. Overall mutation levels, tumor mutation burden, and microsatellite instability were lower, but tumor immune dysfunction and exclusion scores were greater (p < 0.05) in the high-risk group than in the low-risk group. The validation results showed that the prediction model score can accurately predict the prognosis of GC patients. Finally, a nomogram was constructed using the risk score combined with the clinicopathological characteristics of patients with GC.
This risk score from the prediction model related to the tumor microenvironment and immunotherapy could accurately predict the overall survival of GC patients.
Gastric cancer (GC) is the fifth most diagnosed malignancy worldwide . The most common pathological type of GC is stomach adenocarcinoma (STAD), which originates from the epithelial cells of the mucosa . Surgery is the cornerstone treatment for GC , and survival rates have been improved using chemotherapy . Although many individuals with GC have benefited from the recent development of tailored treatment options, such as immunological, targeted, and combination therapies , the mortality rate of GC remains high, making it the third most common cause of cancer-related death . It is therefore important to identify novel prognostic markers and therapeutic targets for patients with GC.
Over 170 different types of modifications for coding and non-coding RNA (ncRNA) have been identified, including N6-methyladenosine (m6A), 5-methylcytosine (m5C), and N1-methyladenosine (m1A) . The methyltransferases, demethylases, and binding proteins, together known as “writers,” “erasers,” and “readers,” can control the dynamics and reversibility of RNA epigenetic modifications [8, 9]. Accumulating evidence suggests that the aberrant expression of RNA modifications is associated with cell survival, proliferation, self-renewal, differentiation, stress adaptation, invasion, and resistance to therapy, all of which are hallmarks of cancer [10,11,12,13,14,15,16,17]. Abnormal post-transcriptional alterations contribute toward cancer cell migration, self-renewal, proliferation, and survival  and thus are promising therapeutic targets for cancer. The most prevalently distributed RNA post-transcriptional modification is m6A . It has been reported that the m6A writers METTL3 and METTL14 have an impact on GC progression [20,21,22], the expression of METTL3 is associated with GC prognosis . Therefore, studying the relationship between RNA modification-related genes (RNAMRGs) and GC would facilitate exploration of potential prognostic markers and therapeutic targets.
In this study, we identified 48 RNAMRGs associated with m6A, m5C, and m1A. Similar to bioinformatics methods of some studies [23,24,25,26], we used The Cancer Genome Atlas (TCGA) data to screen prognosis-related genes and built a predictive model according to the study flowchart shown in Fig. 1 to explore whether RNAMRGs were associated with GC prognosis.
Expression profiles and clinical data of patients with STAD were downloaded from TCGA. GC samples with an overall survival (OS) > 0 were retained, leaving a total of 328 samples (Table 1). "Masked Somatic Mutation" data were downloaded and analyzed. The tumor mutational burden (TMB) and microsatellite instability (MSI) were obtained. “Masked Copy Number Segment” data were used as copy number variation (CNV) data for 440 patients with STAD, who were divided into high- and low-risk groups. Segments were analyzed using GenePattern (https://cloud.genepattern.org) for GISTIC 2.0  with default settings and a confidence level of 0.99. The STAD datasets GSE26899 (93 samples) , GSE26901 (109 samples) , GSE84437 (357 samples) , and GSE62254 (300 samples)  were downloaded from Gene Expression Omnibus (GEO) , and 48 RNAMRGs (Additional file 1: Table S1) were extracted from an article .
RNA modification score (RNAMS) analysis
The RNAMS was calculated for TCGA samples based on the RNAMRGs using the “ssGSEA” algorithm. Samples were divided into high- and low-RNAMS groups, and differentially expressed genes (DEGs) were screened using adj. p < 0.01 and |log fold change (FC)|> 1.
Weighted gene correlation network analysis (WGCNA) was performed using TCGA data as the input, and 0.85 was calculated as the optimal soft threshold by the “pickSoftTreshold” function. With 200 as the minimum number of genes per module, gene modules were dynamically assigned to identify genes. We screened four modules with the strongest positive and negative correlations with the high-RNAMS group. Genes with module membership > 0.5 and significance > 0.2 were screened as key gene modules. DEGs in the high- and low-RNAMS groups were intersected with key gene modules to obtain key genes.
Consensus clustering was performed using TCGA data and the key genes to improve the differentiation between different GC subtypes. The number of clusters was set between 2 and 8, and 80% of all samples were drawn in 1000 repetitions with clusterAlg = “pam” and distance = “spearman.”
Risk prediction model construction
Univariate Cox regression analysis was used to calculate the association between the expression of each key gene and OS. The least absolute shrinkage and selection operator (LASSO) algorithm was then used to eliminate multicollinearity and screen for significant variables. Following multivariate Cox regression analysis, final screening was performed using stepwise regression. The risk score formula was calculated by considering the expression of optimized genes and the multivariate Cox regression coefficients, as follows:
TCGA patients were divided into high- and low-risk groups based on the calculated optimal risk score cutoff. Kaplan–Meier analysis and log-rank tests were performed, and time-dependent receiver operating characteristic (ROC) curves were used to assess survival prediction by calculating the area under curve (AUC).
GSE26901, GSE26899, GSE84437 and GSE62254 were used as validation datasets. After calculating risk scores using the formula above, data were grouped according to the optimal risk score cutoff and subjected to Kaplan–Meier analysis and log-rank tests.
Risk score analysis
Differential analysis was performed based on the high- and low-risk groups from TCGA data, with adj. p < 0.01 and |logFC|> 1.5 as screening thresholds. Gene Ontology (GO) analysis  and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis  were performed on the intersecting genes, the critical value of FDR < 0.05 was considered statistically significant, and the entry screening criteria were adj. p < 0.05 and q-value < 0.05, and the p-value correction method was the Benjamini–Hochberg method. Based on the gene expression profile dataset for patients with STAD in TCGA, we performed gene set enrichment analysis (GSEA) , and adj. p < 0.05 was considered statistically significant.
The STRING (https://string-db.org/) database  was used to construct a protein–protein interaction (PPI) network for DEGs in the high- and low-risk groups, with a coefficient of 0.4. PPI results were exported from STRING and visualized using Cytoscape , and the “CytoHubba” plugin  was used to analyze the bub genes in the PPI network. miRNA–mRNA interaction information was downloaded from the miRTarBas (https://mirtarbase.cuhk.edu.cn/) database . Based on the hub genes identified using the PPI network, an miRNA–mRNA regulatory network was constructed by predicting possible regulated miRNAs.
Immune-related genes were downloaded from a pan-cancer immunogenomic analysis article , which included 782 genes and 28 cell types. The degree of immune cell infiltration was analyzed using TCGA-STAD data, and an immune score is obtained for each tumor sample. Tumor Immune Dysfunction and Exclusion (TIDE)  (http://tide.dfci.harvard.edu) was used to analyze the treatment response in patients with high- and low-risk scores.
Finally, a nomogram was constructed using the risk score and clinicopathological characteristics significantly associated with OS. Calibration curves were generated to assess nomogram performance.
All data processing and analyses were conducted using R (v4.1.1). To compare two groups of continuous variables, the statistical significance of normally distributed variables was estimated using independent Student’s t-tests, and differences between non-normally distributed variables were analyzed using Mann–Whitney U-tests. Chi-square or Fisher’s exact tests were used to compare and analyze statistical significance between two groups of categorical variables. All statistical p values were bilateral and statistically significant at p < 0.05. Significance labeled as NS, p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.
The following packages were used in this study: TCGAbiolinks , maftools , TCGAmutations, GEOquery , GSVA , survival, survminer, limma , WGCNA , Venn, ConsensusClusterPlus , survivor, timeROC , clusterProfiler , ESTIMATE , rms, glmnet , Pheatmap, and ggplot2.
First, to analyze the effect of RNA modification on the process of gastric carcinogenesis, we downloaded the gene expression profiles of STAD patients from the TCGA database as the training set and the GSE26899, GSE26901, GSE84437 and GSE62254 gene expression profiles associated with STAD from the GEO database as the validation sets. The box line plots of the gene expression matrices of GSE26899 (Fig. 2a), GSE26901 (Fig. 2b), GSE84437 (Fig. S1a) and GSE62254 (Additional file 2: Fig. S1b) were plotted. The results showed the same expression trends between samples for these datasets, with no intra-group differences, and can be used for subsequent analysis.
Prognostic analysis of RNAMS in GC and identification of key genes
We calculated the RNAMS for patients with STAD in TCGA to indicate their RNA modification levels. Based on the optimal RNAMS cutoff value (1.207673), TCGA-STAD patients were then divided into high- and low-RNAMS groups. Survival analysis revealed that patients with a high RNAMS had a better prognosis than those with a low RNAMS (Fig. 3a), while clinical analysis showed that patients in the low-RNAMS group had a low age bias and those in the high-RNAMS group had a high age bias (Fig. 3b).
To identify key RNAMRGs in GC, we analyzed DEGs in the high- and low-RNAMS groups (Fig. 3c, d). Of the 275 DEGs identified, 77 were significantly upregulated and 198 were significantly downregulated (Additional file 3: Table S2). Next, we constructed a gene co-expression network to identify biologically significant gene modules as well as genes closely related to RNAMS. Eleven modules were obtained (Fig. 3e) and screened to identify the four modules with the strongest positive and negative correlations with the high-RNAMS group (blue; purple; brown; pink). Intersecting the strongly related genes in the four modules yielded 2,675 module key genes (Additional file 4: Table S3), of which 236 overlapped with the DEGs (Additional file 5: Table S4) and were considered key RNA modification-associated genes in GC (Fig. 3f).
Analysis of key genes and molecular typing of RNAMS-related genes
To investigate the overall expression of the 236 key genes in patients with STAD, we plotted a heatmap (Fig. 4a) and analyzed their correlation (Additional file 6: Table S5). Fifty genes were randomly selected for visualization (Fig. 4b). Most of the key genes were highly expressed in the low-RNAMS group, and the correlations between genes were generally positive. When we further analyzed the DNA level variation of key genes, we found that single nucleotide polymorphisms (SNPs) were present in key genes in 193 (60.12%) of TCGA-STAD samples (Fig. 4c).
Next, we used the expression profiles of the 236 key genes in all samples to perform unsupervised clustering using a consensus clustering algorithm. The optimal number of clusters was calculated and determined to be three (Fig. 4d–g); therefore, k = 3 was used to cluster all samples into three subtypes that were subjected to principal components analysis (Fig. 4h). The three subtypes were well differentiated, and prognostic analysis revealed that survival differed significantly among the three subtypes (Fig. 4i), indicating the accuracy of the clustering results.
RNAMS risk model construction and evaluation
To quantify the effect of the 236 key genes on prognosis, we combined their expression to construct a risk score model. First, the key genes were screened using univariate Cox regression, and 66 genes were retained for LASSO regression to eliminate covariance (Fig. 5a, Additional file 7: Table S6). After the best lambda values had been determined and crossed by ten-fold validation, six genes (SELP, APOD, MXRA8, CST2, RRAD, and GPX3; Fig. 5b) were subjected to multivariate Cox regression analysis, and the optimal combination was screened using stepwise regression to obtain two genes (SELP and CST2) for model construction (Fig. 5c). The model scoring formula was as follows:
To verify the validity of the model, the correlation between risk score and RNAMS for TCGA patients was plotted (Fig. 5d), and a significant negative correlation was found (cor = -0.49, p < 2.2e−16). A risk score cutoff value of 0.4448089 was used to classify the training set TCGA-STAD patients into high- and low-risk groups. Survival analysis showed that the low-risk group had significantly better survival than the high-risk group (Fig. 5e). Consequently, patients in the validation sets were classified into high- and low-risk groups based on their respective risk score cutoff values (GSE26899: 2.059207; GSE26901: 2.140648; GSE84437: 2.156697; GSE62254: 1.788483), and survival analysis showed that the low-risk groups had significantly better survival than the high-risk groups (Fig. 5f, g; Additional file 8: Fig. S2a, b). Analyses for risk score distribution, survival status, and characteristic gene expression patterns for TCGA (Fig. 5h), GSE26899 (Fig. 5i), GSE26901 (Fig. 5j), GSE84437 (Additional file 8: Fig. S2c), and GSE62254 (Additional file 8: Fig. S2d) revealed that patients in high-risk groups had worse survival and similar gene expression patterns. Time-dependent ROC analysis of risk scores in the five datasets revealed the AUC values for 1-, 3-, and 5-year OS in TCGA (Fig. 5k), GSE26899 (Fig. 5l), GSE26901 (Fig. 5m), GSE84437 (Additional file 8: Fig. S2e) and GSE62254 (Additional file 8: Fig. S2f) indicating that the risk score could accurately predict the OS of patients with GC.
To determine the ability of the model to predict cancer development in patients with STAD, we performed differential expression analysis in high- and low-risk patient groups. Of the 144 DEGs identified, 2 were upregulated and 142 were significantly downregulated (Fig. 6a, b). Next, we performed GO and KEGG enrichment analyses to identify biological processes, molecular functions, cellular components, and biological pathways related to the 144 DEGs (Fig. 6c–f, Additional file 9: Tables S7, Additional file 10: Table S8), which included the Wnt signaling pathway and vascular smooth muscle contraction.
To elucidate the effect of gene expression levels on GC, we analyzed the associations between biological processes involved in gene expression in TCGA data using GSEA. The MYC targets V2 pathway (Fig. 7a) was significantly enriched in high-risk patients, whereas KRAS signaling up (Fig. 7b), inflammatory response (Fig. 7c), interferon gamma response (Fig. 7d), IL2–STAT5 signaling (Fig. 7e), KRAS signaling dn (Fig. 7f), TGFβ signaling (Fig. 7g), apoptosis (Fig. 7h), and IL6–JAK–STAT3 signaling (Fig. 7i) were significantly enriched in low-risk patients (Additional file 11: Table S9).
Construction of PPI and related regulatory networks
To analyze the protein interactions among the 144 DEGs, we constructed a PPI network (Fig. 8a). The top 10 hub genes (ELN, DCN, MYH11, ACTA2, FBLN5, TAGLN, CLU, MFAP5, MFAP4, and FBLN1) were obtained based on local node density, with a darker color indicating node importance (Fig. 8b). Functional similarity (Friends) analysis revealed that ELN was an important hub gene (Fig. 8c), while correlation analysis indicated significant co-expression patterns between hub genes and risk scores, all of which correlated positively (Fig. 8d). Based on miRNA–mRNA interaction information downloaded from the miRTarBase database, an miRNA–mRNA regulatory network was constructed using the hub genes obtained from the PPI network (Fig. 8e), which contained a total of 202 miRNAs and 9 mRNAs worthy of further study.
Immune cell infiltration analysis
To evaluate the effect of the risk score on the overall immune characteristics and immune cell infiltration levels of patients with STAD, we calculated immune and stromal scores for each patient. A significant positive correlation was observed between risk score and both the immune score (cor = 0.37, p = 4.9e−12, Fig. 9a) and stromal score (cor = 0.73, p < 2.2e−16, Fig. 9b). Immune cell infiltration scores in TCGA were displayed using a heatmap to further examine immune cell expression (Fig. 9c). Most immune cells correlated positively with one another (Fig. 9d), with > 90% of immune cells in both high- and low-risk groups showing a significant difference and most immune cells displaying significantly higher infiltration in the high-risk group than in the low-risk group (Fig. 9e).
Next, we plotted the correlation heatmap between hub genes and immune cells (Fig. 9f) and between hub genes and immune genes (Fig. 9g) to determine the effect of hub genes on immune function in patients with STAD. Differential boxplots revealed that the expression of all immune genes differed significantly between the high- and low-risk groups (Fig. 9h).
Effect of risk score on genomic changes in patients with STAD
To determine the effect of the risk score on genetic variation in patients with STAD, we analyzed changes in SNPs and CNVs. The overall level of single nucleotide mutations in common tumorigenesis driver genes was lower in the high-risk group than in the low-risk group (Fig. 10a, b), whereas the frequency of CNVs was significantly increased in all high-risk patients, particularly for copy number deletions (Fig. 10c, d). By calculating and plotting the TIDE, TMB, and MSI, we found that TIDE scores were significantly higher in the high-risk group than in the low-risk group (p = 2.6e−06, Fig. 10e), whereas TMB (p = 1.8e−05, Fig. 10f) and MSI (p = 6.7e−05, Fig. 10g) were significantly lower in the high-risk group than in the low-risk group. Together, these results indicate that the model risk score can be a potential indicator for evaluating the efficacy of immunotherapy in patients with GC.
Construction of a clinical predictive nomogram
Finally, we generated boxplots to visualize the differences between risk score and clinical characteristics. Risk score differed significantly with pathological stage, T stage, and N stage (Fig. 11a), and the RNAMS was significantly lower in the high-risk group than in the low-risk group (p = 1.5e−08, Fig. 11b). To determine whether the risk score and clinicopathological characteristics were independent prognostic factors, we performed regression analyses. Univariate Cox regression analysis showed that risk score (p < 0.001), pathological stage (p = 0.00159), M stage (p = 0.00243), T stage (p = 0.00565), and gender (p = 0.0473) were all associated with OS (Fig. 11c, Table 2), and multivariate Cox regression analysis of these significant factors revealed that risk score (p = 0.001), M stage (p < 0.001), T stage (p = 0.014), and gender (p = 0.005) were significantly associated with OS (Fig. 11d, Table 3). Therefore, we combined the risk score with these three clinicopathological characteristics to construct a predictive nomogram to assess the OS of patients with STAD (Fig. 11e). Calibration curves revealed good uniformity between the OS estimates of the nomograph for 1-, 2- and 3-year OS of patients and the actual observed values (Fig. 11f), suggesting that the nomogram predictions were accurate.
The ROC curves of the prognostic discrimination of patients at 1 year, 2 years and 3 years of the nomogram were plotted (Fig. 11g), and the results showed that the AUC values of the nomogram for predicting the prognosis of patients at 1 year, 2 years and 3 years were 0.696, 0.683 and 0.655, respectively, indicating that the nomogram has good predictive ability.
Although the incidence and mortality rates of GC have declined in the past 5 decades [53,54,55], GC remains the third leading cause of cancer-related deaths . Cancer cells are characterized by genomic instability, which favors mutation accumulation and increased tumor heterogeneity [57,58,59]. RNA modifications have recently emerged as key post-transcriptional regulators of gene expression and are thought to be associated with various diseases, including cancer . However, the intrinsic relationship between RNA modifications and GC progression remains unknown. In this study, we analyzed the relationship between RNAMRGs and GC using bioinformatics analyses and statistical methods and constructed a prognostic model.
First, we scored each TCGA sample and grouped the samples according to their RNAMS and found that patients with high RNA modification correlation had a worse prognosis. Dynamic RNA modification events can promote tumor progression by promoting tumor cell proliferation or regulating invasive and metastatic potential, resulting in a poor prognosis, consistent with our findings. Therefore, we explored hub RNAMRGs and established a prognostic prediction model for patients with GC using selectin P (SELP) and cystatin SA (CST2). SELP encodes a protein stored on platelets and endothelial cells that is redistributed to the plasma membrane during platelet activation and degranulation and mediates the interactions between activated endothelial cells or platelets and leukocytes . Coxsackievirus B 1 and 3 interact with human platelets to trigger SELP expression , and SELP has predictive value in various cancers [63, 64]; however, the relationship between SELP and RNA modification in GC carcinogenesis and development remains unclear. CST2 overexpression enhances the growth, migration, and invasion of GC cells by regulating the epithelial–mesenchymal transition and TGF-β1 signaling pathways . According to the formula of our predictive model, higher SELP and CST2 expression were associated with a higher risk score, and survival analysis further indicated that the high-risk group had a worse prognosis. The high predictive efficacy of the model was further validated using GEO datasets. Taken together, these findings are consistent with the hypothesis that SELP and CST2 may act as cancer progression-promoting genes in various tumors.
Furthermore, enrichment analysis revealed the importance of the Wnt signaling pathway, which plays a key role in tumorigenesis [66, 67]. While constructing the PPI network, we identified ELN as an important hub gene. ELN encodes one of the two components of elastic fibers that form part of the extracellular matrix and confer elasticity to organs and tissues. Abnormal ELN levels have been observed in many fibrotic diseases, including kidney , lung , and liver fibrosis , and ELN accumulation is associated with the development of hepatocellular carcinoma . ELN has also been shown to regulate tumor development and the tumor microenvironment (TME) in colorectal cancer . Therefore, studying the mechanism of ELN and RNA modification in GC carcinogenesis is essential. miRNAs are epigenetic regulators that affect the levels of proteins encoded by target mRNAs without modifying gene sequences and are themselves regulated by epigenetic mechanisms . Reciprocal relationships between miRNAs and epigenetic regulation form miRNA-epigenetic feedback loops that can regulate cellular processes such as cell proliferation , apoptosis , and differentiation . To understand the regulatory relationship between the hub genes (e.g., ELN) and miRNAs, we constructed an miRNA–mRNA regulatory network to explore whether miRNAs and mRNAs are associated with GC development.
Recently, breakthroughs have been made in immunotherapy for cancers , including GC [78, 79]; however, not all patients with GC benefit from immunotherapy, possibly owing to the TME . Here, we found that immune cell infiltration was increased in the high-risk group. FBLN1 was highly expressed in most immune cells and correlated positively with the majority of immune genes, whereas FBLN5 was expressed at low levels in most immune cells. FBLN1 and FBLN5 are components of extracellular matrix fibronectin, and ectopic FBLN1 expression inhibits the growth of GC cells by inducing apoptosis . Additionally, FBLN5 is a potential indicator of therapeutic efficacy in patients with hepatocellular carcinoma . Thus, FBLN1 and FBLN5 may be potential immunotherapeutic targets for treating GC. A comparison of the high- and low-risk groups revealed that all immune genes exhibited significantly different expression patterns, suggesting that our predictive model is a potential indicator for immunotherapy in patients with GC.
We also found that overall mutation levels for common tumor driver genes were lower in the high-risk group than in the low-risk group, while CNVs, particularly due to deletion, were significantly increased in all patients in the high-risk group. TIDE, TMB, and MSI results further suggested that our model could predict the effect of immune checkpoint inhibition therapy in patients with GC. We therefore combined the model with clinicopathological features to construct a nomogram to assess the clinical prognosis of GC. Despite these promising findings, our study had some limitations. First, the model has not yet been applied in a clinical setting, and its actual predictive accuracy remains unknown. Second, the value and mechanism of action of the hub RNAMRGs identified in this study have not been experimentally verified in patients with GC. To overcome these limitations, we will apply the nomogram in future clinical work and conduct further experiments to explore the mechanism of key RNA modifications in the development of GC.
In conclusion, the RNAMRG risk prediction model developed in this study could effectively predict the overall survival of patients with GC. Thus, the model risk score can be used to study GC carcinogenesis and developing targeted therapies for GC.
Availability of data and materials
The datasets analyzed during the current study are available in The Cancer Genome Atlas (https://portal.gdc.cancer.gov) and the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26899, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26901, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84437 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62254).All data used in this study are available from public databases and can be provided by the authors upon reasonable request.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424.
Martínez-Barriocanal Á, Arango D, Dopeso H. PVT1 long non-coding RNA in gastrointestinal cancer. Front Oncol. 2020;10:38.
Kwon IG, Son T, Kim HI, Hyung WJ. Fluorescent lymphography-guided lymphadenectomy during robotic radical gastrectomy for gastric cancer. JAMA Surg. 2019;154(2):150–8.
Sexton RE, Al Hallak MN, Diab M, Azmi AS. Gastric cancer: a comprehensive review of current and future treatment strategies. Cancer Metastasis Rev. 2020;39(4):1179–203.
Li K, Zhang A, Li X, Zhang H, Zhao L. Advances in clinical immunotherapy for gastric cancer. Biochim Biophys Acta Rev Cancer. 2021;1876(2):188615.
Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, Lordick F. Gastric cancer. Lancet. 2020;396(10251):635–48.
Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer. 2020;20(6):303–22.
Niu Y, Zhao X, Wu YS, Li MM, Wang XJ, Yang YG. N6-methyl-adenosine (m6A) in RNA: an old modification with a novel epigenetic function. Genomics Proteom Bioinform. 2013;11(1):8–17.
Xue C, Zhao Y, Li L. Advances in RNA cytosine-5 methylation: detection, regulatory mechanisms, biological functions and links to cancer. Biomark Res. 2020;8:43.
Blanco S, Bandiera R, Popis M, Hussain S, Lombard P, Aleksic J, et al. Stem cell function and stress response are controlled by protein synthesis. Nature. 2016;534(7607):335–40.
Blanco S, Dietmann S, Flores JV, Hussain S, Kutter C, Humphreys P, et al. Aberrant methylation of tRNAs links cellular stress to neuro-developmental disorders. EMBO J. 2014;33(18):2020–39.
Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, et al. m6A RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells. Cell Rep. 2017;18(11):2622–34.
Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M, et al. Stem cells. m6A mRNA methylation facilitates resolution of naïve pluripotency toward differentiation. Science. 2015;347(6225):1002–6.
Janin M, Ortiz-Barahona V, de Moura MC, Martínez-Cardús A, Llinàs-Arias P, Soler M, et al. Epigenetic loss of RNA-methyltransferase NSUN5 in glioma targets ribosomes to drive a stress adaptive translational program. Acta Neuropathol. 2019;138(6):1053–74.
Su R, Dong L, Li C, Nachtergaele S, Wunderlich M, Qing Y, et al. R-2HG exhibits anti-tumor activity by targeting FTO/m6A/MYC/CEBPA signaling. Cell. 2018;172(1–2):90-105.e23.
Weng H, Huang H, Wu H, Qin X, Zhao BS, Dong L, et al. METTL14 inhibits hematopoietic stem/progenitor differentiation and promotes leukemogenesis via mRNA m6A modification. Cell Stem Cell. 2018;22(2):191-205.e9.
Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, et al. FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014;24(12):1403–19.
Boriack-Sjodin PA, Ribich S, Copeland RA. RNA-modifying proteins as anticancer drug targets. Nat Rev Drug Discov. 2018;17(6):435–53.
Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m6A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 2018;28(6):616–24.
Fan HN, Chen ZY, Chen XY, Chen M, Yi YC, Zhu JS, et al. METTL14-mediated m6A modification of circORC5 suppresses gastric cancer progression by regulating miR-30c-2-3p/AKT1S1 axis. Mol Cancer. 2022;21(1):51.
Wang Q, Chen C, Ding Q, Zhao Y, Wang Z, Chen J, et al. METTL3-mediated m6A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut. 2020;69(7):1193–205.
Yue B, Song C, Yang L, Cui R, Cheng X, Zhang Z, et al. METTL3-mediated N6-methyladenosine modification is critical for epithelial-mesenchymal transition and metastasis of gastric cancer. Mol Cancer. 2019;18(1):142.
Liu Y, Wang J, Li L, Qin H, Wei Y, Zhang X, et al. AC010973.2 promotes cell proliferation and is one of six stemness-related genes that predict overall survival of renal clear cell carcinoma. Sci Rep. 2022;12(1):4272.
Wei Y, Chen X, Ren X, Wang B, Zhang Q, Bu H, et al. Identification of MX2 as a novel prognostic biomarker for sunitinib resistance in clear cell renal cell carcinoma. Front Genet. 2021;12:680369.
Wu D, Yin Z, Ji Y, Li L, Li Y, Meng F, et al. Identification of novel autophagy-related lncRNAs associated with a poor prognosis of colon adenocarcinoma through bioinformatics analysis. Sci Rep. 2021;11(1):8069.
Yu L, Shen H, Ren X, Wang A, Zhu S, Zheng Y, et al. Multi-omics analysis reveals the interaction between the complement system and the coagulation cascade in the development of endometriosis. Sci Rep. 2021;11(1):11926.
Reich M, Tabor T, Liefeld T, Thorvaldsdóttir H, Hill B, Tamayo P, et al. The GenePattern notebook environment. Cell Syst. 2017;5(2):149-51.e1.
Oh SC, Sohn BH, Cheong JH, Kim SB, Lee JE, Park KC, et al. Clinical and genomic landscape of gastric cancer with a mesenchymal phenotype. Nat Commun. 2018;9(1):1777.
Yoon SJ, Park J, Shin Y, Choi Y, Park SW, Kang SG, et al. Deconvolution of diffuse gastric cancer and the suppression of CD34 on the BALB/c nude mice model. BMC Cancer. 2020;20(1):314.
Cristescu R, Lee J, Nebozhyn M, Kim KM, Ting JC, Wong SS, et al. Molecular analysis of gastric cancer identifies subtypes associated with distinct clinical outcomes. Nat Med. 2015;21(5):449–56.
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets–update. Nucl Acids Res. 2013;41(Database issue):D991–5.
Wang E, Li Y, Ming R, Wei J, Du P, Zhou P, et al. The prognostic value and immune landscapes of a m6A/m5C/m1A-related LncRNAs signature in head and neck squamous cell carcinoma. Front Cell Dev Biol. 2021;9:718974.
Gene Ontology Consortium. going forward. Nucl Acids Res. 2015;43(Database issue):1049–56.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4(Suppl 4)):S11.
Huang HY, Lin YC, Li J, Huang KY, Shrestha S, Hong HC, et al. miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res. 2020;48(D1):D148–54.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18(1):248–62.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8.
Silva TC, Colaprico A, Olsen C, D’Angelo F, Bontempi G, Ceccarelli M, et al. TCGA Workflow: analyze cancer genomics and epigenomics data using bioconductor packages. F1000Res. 2016;5:1542.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.
Davis S, Meltzer PS. GEOquery: a bridge between the gene expression omnibus (GEO) and bioconductor. Bioinformatics. 2007;23(14):1846–7.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7.
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.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
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.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22.
Arnold M, Karim-Kos HE, Coebergh JW, Byrnes G, Antilla A, Ferlay J, et al. Recent trends in incidence of five common cancers in 26 European countries since 1988: analysis of the European cancer observatory. Eur J Cancer. 2015;51(9):1164–87.
Ferro A, Peleteiro B, Malvezzi M, Bosetti C, Bertuccio P, Levi F, et al. Worldwide trends in gastric cancer mortality (1980–2011), with predictions to 2015, and incidence by subtype. Eur J Cancer. 2014;50(7):1330–44.
Luo G, Zhang Y, Guo P, Wang L, Huang Y, Li K. Global patterns and trends in stomach cancer incidence: Age, period and birth cohort analysis. Int J Cancer. 2017;141(7):1333–44.
Thrift AP, El-Serag HB. Burden of gastric cancer. Clin Gastroenterol Hepatol. 2020;18(3):534–42.
Böttcher R, Kweldam CF, Livingstone J, Lalonde E, Yamaguchi TN, Huang V, et al. Cribriform and intraductal prostate cancer are associated with increased genomic instability and distinct genomic alterations. BMC Cancer. 2018;18(1):8.
Perkhofer L, Schmitt A, Romero Carrasco MC, Ihle M, Hampp S, Ruess DA, et al. ATM deficiency generating genomic instability sensitizes pancreatic ductal adenocarcinoma cells to therapy-induced DNA damage. Cancer Res. 2017;77(20):5576–90.
Tubbs A, Nussenzweig A. Endogenous DNA damage as a source of genomic instability in cancer. Cell. 2017;168(4):644–56.
Nombela P, Miguel-López B, Blanco S. The role of m6A, m5C and Ψ RNA modifications in cancer: novel therapeutic opportunities. Mol Cancer. 2021;20(1):18.
de Stoppelaar SF, Van’t Veer C, Roelofs JJ, Claushuis TA, de Boer OJ, Tanck MW, et al. Platelet and endothelial cell P-selectin are required for host defense against Klebsiella pneumoniae-induced pneumosepsis. J Thromb Haemost. 2015;13(6):1128–38.
Negrotto S, Jaquenod de Giusti C, Rivadeneyra L, Ure AE, Mena HA, Schattner M, et al. Platelets interact with Coxsackieviruses B and have a critical role in the pathogenesis of virus-induced myocarditis. J Thromb Haemost. 2015;13(2):271–82.
Avan A, Avan A, Le Large TY, Mambrini A, Funel N, Maftouh M, et al. AKT1 and SELP polymorphisms predict the risk of developing cachexia in pancreatic cancer patients. PLoS ONE. 2014;9(9):e108057.
Powrózek T, Mlak R, Brzozowska A, Mazurek M, Gołębiowski P, Małecka-Massalska T. Relationship between -2028 C/T SELP gene polymorphism, concentration of plasma P-selectin and risk of malnutrition in head and neck cancer patients. Pathol Oncol Res. 2019;25(2):741–9.
Zhang WP, Wang Y, Tan D, Xing CG. Cystatin 2 leads to a worse prognosis in patients with gastric cancer. J Biol Regul Homeost Agents. 2020;34(6):2059–67.
Togasaki K, Sugimoto S, Ohta Y, Nanki K, Matano M, Takahashi S, et al. Wnt signaling shapes the histologic variation in diffuse gastric cancer. Gastroenterology. 2021;160(3):823–30.
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.
Sun Q, Baues M, Klinkhammer BM, Ehling J, Djudjaj S, Drude NI, et al. Elastin imaging enables noninvasive staging and treatment monitoring of kidney fibrosis. Sci Transl Med. 2019;11(486):eaat4865.
Mecham RP. Elastin in lung development and disease pathogenesis. Matrix Biol. 2018;73:6–20.
Masugi Y, Abe T, Tsujikawa H, Effendi K, Hashiguchi A, Abe M, et al. Quantitative assessment of liver fibrosis reveals a nonlinear association with fibrosis stage in nonalcoholic fatty liver disease. Hepatol Commun. 2017;2(1):58–68.
Yasui Y, Abe T, Kurosaki M, Higuchi M, Komiyama Y, Yoshida T, et al. Elastin fiber accumulation in liver correlates with the development of hepatocellular carcinoma. PLoS ONE. 2016;11(4):e0154558.
Li J, Xu X, Jiang Y, Hansbro NG, Hansbro PM, Xu J, et al. Elastin is a key factor of tumor development in colorectal cancer. BMC Cancer. 2020;20(1):217.
Yao Q, Chen Y, Zhou X. The roles of microRNAs in epigenetic regulation. Curr Opin Chem Biol. 2019;51:11–7.
Kiga K, Mimuro H, Suzuki M, Shinozaki-Ushiku A, Kobayashi T, Sanada T, et al. Epigenetic silencing of miR-210 increases the proliferation of gastric epithelium during chronic Helicobacter pylori infection. Nat Commun. 2014;5:4497.
Matsushima K, Isomoto H, Yamaguchi N, Inoue N, Machida H, Nakayama T, et al. MiRNA-205 modulates cellular invasion and migration via regulating zinc finger E-box binding homeobox 2 expression in esophageal squamous cell carcinoma cells. J Transl Med. 2011;9:30.
Chen JF, Mandel EM, Thomson JM, Wu Q, Callis TE, Hammond SM, et al. The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation. Nat Genet. 2006;38(2):228–33.
Rahman M, Sawyer WG, Lindhorst S, Deleyrolle LP, Harrison JK, Karachi A, et al. Adult immuno-oncology: using past failures to inform the future. Neuro Oncol. 2020;22(9):1249–61.
Lote H, Cafferkey C, Chau I. PD-1 and PD-L1 blockade in gastrointestinal malignancies. Cancer Treat Rev. 2015;41(10):893–903.
Procaccio L, Schirripa M, Fassan M, Vecchione L, Bergamo F, Prete AA, et al. Immunotherapy in gastrointestinal cancers. Biomed Res Int. 2017;2017:4346576.
Davidson M, Okines AF, Starling N. Current and future therapies for advanced gastric cancer. Clin Colorectal Cancer. 2015;14(4):239–50.
Cheng YY, Jin H, Liu X, Siu JM, Wong YP, Ng EK, et al. Fibulin 1 is downregulated through promoter hypermethylation in gastric cancer. Br J Cancer. 2008;99(12):2083–7.
Chen X, Wang P, Yang M, Zhou W, Chen J, Meng Z, et al. Therapeutic effect of Jianpi Liqi Fang combined with transcatheter arterial chemoembolization in patients with hepatocellular carcinoma and spleen deficiency syndrome. J Tradit Chin Med. 2021;41(1):157–66.
This work was supported by People’s Hospital of Xinjiang Uygur Autonomous Region.
Ethics approval and consent to participate
Consent for publication
None of the authors have any conflicts of interest to disclose.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
RNA modification-related genes.
Boxplots of GSE84437 (a) and GSE62254 (b) expression profile data.
Differentially expressed genes in the high- and low-RNA modification score groups.
Key gene modules.
Correlation among all key genes.
Univariate Cox regression analysis results for all key genes.
Evaluation of RNA modification score (RNAMS) risk prediction model. a–b Survival curves for GSE84437 (a) and GSE62254 (b); c–d Risk score distribution for GSE84437 (i) and GSE62254 (j) expression profiles; e–f ROC curves for GSE84437 (e) and GSE26901 (f).
Gene Ontology enrichment analysis results.
Kyoto Encyclopedia of Genes and Genomes enrichment analysis results.
Gene set enrichment analysis results.
About this article
Cite this article
Tuhongjiang, A., Wang, F., Zhang, C. et al. Construction of an RNA modification-related gene predictive model associated with prognosis and immunity in gastric cancer. BMC Bioinformatics 24, 147 (2023). https://doi.org/10.1186/s12859-023-05283-3