From previous systematic reviews and the Molecular Signatures Database, a total of 29 m7G-related genes were extracted (11–14). RNA sequencing datasets of 32 normal gastric tissues and 375 gastric cancer specimens were downloaded from The Cancer Genome Atlas (https://portal.gdc.cancer.gov/), together with the corresponding clinical data. All raw data were preprocessed with the limma R package and screened for differentially expressed genes (DEGs) associated with m7G. A Wilcox test was performed using |logFC(foldchange)| ≥ 1 and false discovery rate (FDR) < 0.05 as the cutoff criteria. The training group comprised data from the TCGA cohort, and the validation group comprised gene expression data and corresponding post-operative data collected from January 2015 to December 2020 in 70 human GC patients admitted at Taizhou People’s Hospital. All participants signed written informed consent forms prior to their participation. The study protocol was approved by the hospital ethics committee. Clinical information including age, gender, and TNM stage was collected ( Tables 1 , 2 ).

R software (R Foundation for Statistical Computing, Vienna, Austria) with the Survival package was used to conduct univariate and multivariate Cox analyses. The optimal prediction model was determined as per the Akaike Information Criterion (AIC). Based on the results of the univariate Cox regression analysis, we used a support vector machine (SVM) approach to further screen for prognostic genes with the help of the R package e1071(http://cran.r-project.org/package=e1071) (15). The risk score for each patient was accorded based on the expression level of m7G-related genes and the regression coefficient β of the weighted linear combination in the multifactorial analysis, calculated as: risk score = βgene1 × expr(gene1) + βgene2 × expr(gene2) +… + βgeneN × expr(geneN). Based on the median risk score, GC patients were divided into low-risk and high-risk groups. The Kaplan-Meier (K-M) method was used to clarify the relationship between risk values and survival time. Furthermore, ROC curves were drawn to assess the predictive effect of the model; and the risk and survival status were plotted for high- and low-risk groups. In addition, principal component analysis (PCA) was performed using the “prcomp” package to visualise the grouping of data and the distribution of different groups. Based on the above risk score formula, risk scores were calculated for patients in the validation group, who were subsequently divided into high- and low-risk groups and analysed using the K-M method.

Risk scores were combined with clinical information including the survival time, survival status, sex, age, tumour grade, pathological grade, TNM stage, and risk score. Clinical factors affecting GC prognosis were initially screened using univariate Cox regression analysis, and factors significantly associated with GC prognosis were subsequently subjected to multivariate Cox regression analysis.

The biological functions of these differentially expressed m7G-genes were examined integrally by Gene Ontology term enrichment analysis (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (16). All enrichment analyses were done using the org.Hs.eg.db, DOSE, clusterProfiler, and enrichplot packages. P<0.05 while FDR<0.05 was considered a statistically significant difference. The Drug Signatures Database (DSigDB) currently holds 17,389 unique compounds and 19,531 drug target genes, and is freely available at http://tanlab.ucdenver.edu/DSigDB (17). We used the DSigDB database to identify potential drugs that significantly interacted with genes (18). P<0.05 was considered a statistically significant difference.

Infiltration of immune cells into the TME is known to influence tumour progression (19). A systematic analysis of the relationship between tumour immune infiltration and hub genes can throw light on the mechanism of tumour immune escape and facilitate the development of new biomarkers and therapeutic strategies. Thus, we checked for correlation between the expression of hub genes and immune infiltration in GC. By single-sample gene set enrichment analysis (ssGSEA), the infiltration of 28 immune cells was analysed, with the ssGSEA score set as the standard (20).

The cancerous tissues and normal tissues adjacent to the tumour were surgically excised, added to RNAlater protective solution, and frozen at -80°C for further examination. The total RNA was extracted by the TRIzol method. After centrifugation, 1 μg of total RNA was extracted from the supernatant, and cDNA was synthesised using a cDNA reverse transcription kit. The cDNA was used as the template and the GAPDH gene as the internal reference for real-time PCR. A two-step standard PCR amplification procedure was used: the first step was pre-denaturation at 95°C for 30 s for 1 cycle; the second step was PCR reaction at 95°C for 5 s, 60°C for 15 s, 72°C for 30 s for 40 cycles. Three replicate wells were set up for each sample and the relative expression of the target gene was calculated using the 2-ΔΔCT method. The primers used are as follows: NUDT10 (forward 5′- CGGTCCGAGAGGTGTACGA -3′, reverse 5′- AATCTTCCCAATCCTCCAGCA -3′).

METTL1 (forward 5′-GGCAACGTGCTCACTCCAA-3′, reverse 5′- CACAGCCTATGTCTGCAAACT-3′); NUDT4 (forward 5′-ACCAGTGGATTGTCCCAGGA-3′, reverse 5′-CCCAGAAGTCTGCCTAGTTTTC-3′); GEMIN5 (forward 5′-CCTCCGTCTTCCTTGTCCG-3′, reverse 5′- CAGAGACCCTTTCGGTGTGTC-3′); EIF4E1B (forward 5′-GGACTTCTGGGCGCTATACAG-3′, reverse 5′-CTTGAAGAGGGCGTAGTCACA-3′); and DCPS (forward 5′- GCAGCTCCTCAACTAGGCAAG-3′, reverse 5′-GAAGCCGGAGAACGGTAAGC-3′).

The flow chart was shown in Figure 1 . Thirty-two normal tissue samples and 375 GC tissue samples were obtained from the TCGA database. As per the screening criteria, a total of 20 DEGs were identified out of 29 m7G-related genes ( Figure 2A , green: low expression level; red: high expression level). Figure 2B shows the relationship between different m7G-related DEGs, with a significant positive correlation observed between most m7G-related genes. Nine DEGs were associated with overall survival in the univariate regression analysis( Figure 3A ). After screening these genes by SVM, six genes were obtained( Figure 3B ). They were then subjected to a multivariate Cox analysis, weighted according to the relative coefficients in the multivariate Cox regression. The scoring equation was:

Risk score = NUDT10×0.021+METTL1×0.0223-NUDT4×0.00570+GEMIN5×1.44+EIF4E1B×0.175-DCPS×0.0254. GC patients were divided into high-risk and low-risk groups based on the median risk score ( Figure 4A , median risk score=1.09). In addition, the ROC curve analysis showed that the area under the curve (AUC) was 0.814 at 1 year, 0.798 at 3 years, and 0.714 at 5 years( Figure 4B ), all of which were greater than 0.7. The calibration curve ( Figure S1 ) showed that the prediction results of the prognostic risk score model were in good agreement with the actual observations. The risk model described in this study demonstrates a good degree of discrimination and calibration. The K-M survival analysis showed a significant difference in the survival of the high-risk and low-risk groups ( Figure 4C , p<0.05), with the high-risk group having a significantly lower survival rate than the low-risk group. The survival rate also decreased over time. PCA analysis showed that the different risk groups were distributed in two directions, suggesting that a six-gene prognostic risk score model could distinguish well between high- and low-risk groups ( Figure 4D ).

The β’s calculated from the TCGA cohort and the expression levels of m7G-related genes in the validation cohort were used to calculate risk scores for patients in the validation cohort. PCR results demonstrate that six genes are still differentially expressed in the validation cohort( Figure 5A ). Patients with GC were divided into high-risk and low-risk groups based on the median risk score( Figure 5B , median risk score=1.06). The results showed that compared with the low-risk group, the high-risk group had a worse prognosis ( Figure 5C ). Furthermore, the AUC of the 6-gene signature was 0.778, 0.761, and 0.714 at 1, 3, and 5 years, respectively( Figure 5D ). On a calibration graph, the validation group revealed a moderate calibration ( Figure S2 ).

Univariate and multivariate Cox regression analyses were performed to evaluate whether the prognostic risk score model could be used as an independent prognostic predictor. The risk score was found to be significantly associated with overall survival (OS) in patients with GC, and can thus be used as an independent prognostic factor for evaluating patients with GC ( Table 3 ). Considering the clinical utility of the risk model, we incorporated clinical parameters to establish a nomogram to predict the 1-, 3-, 5- OS ( Figure S3 ). The discrimination of the nomogram was assessed using the ROC curve. Figures S4 , S5 show the ROC curves of the nomogram in modeling group and validation group respectively, and the results prove that the nomogram have good discrimination.

The differentially expressed m7G-related genes were enriched in 144 biological processes (BP), 16 cellular components (CC), and 30 molecular functions (MF). Among GO-BP, there was significant enrichment of differential genes in the regulation of ‘cellular amide metabolic process,’ ‘regulation of translation,’ and ‘translational initiation.’ Significantly enriched GO-CC included ‘eukaryotic translation initiation factor 4F complex,’ ‘RNA cap binding complex,’ and ‘mRNA cap binding complex.’ The significantly enriched GO-MF included ‘translation initiation factor activity,’ ‘RNA 7−methylguanosine cap binding,’ and ‘RNA cap binding.’ Enrichment analysis by KEGG showed that the differentially expressed genes were mainly enriched in the ‘RNA degradation pathway,’ ‘EGFR tyrosine kinase inhibitor resistance pathway’, and ‘RNA transport pathway’ ( Figure S6 ).

To further explore the correlation between the prognostic risk score model and immune status, we quantified different immune cell subsets, cell functions, and pathways using ssGSEA. The scores of most immune cells, including CD8+_T_cells, iDCs, Mast_cells, pDCs, Th1_cells, Th2_cells, and macrophages, were significantly different between the high-risk and low-risk groups. Th1_cells, Th2_cells, CD8+ T cells, and APC co-stimulation scores, inflammation−promoting, T_cell_co−stimulation were higher in the low-risk group, indicating that the antigen presentation process and cellular and humoral immunity were better in the low-risk group ( Figures 6A, B ). Therefore, the weakening of anti-tumour immunity in high-risk patients may be a reason for their poor prognosis. We further analysed the relationship of the six genes with immune cells in gastric cancer, and found that DCPS, GEMIN5, and METTL1 had a significant positive correlation with Th2 cells, suggesting that they are synergistic in immune cell activation. In contrast, EIF4E1B and NUDT10 had a negative correlation with Th2 ( Figure 7 ).

Table S1 reveal the candidate drugs from the DSigDB database for the DEGs. Our results suggest that Trichostatin A and scriptaid may serve as a potential agent for the treatment of GC. Trichostatin A and scriptaid are histone deacetylase inhibitors (HDACi) (21). Studies have shown that some HDACi have the ability to resist tumour development through a variety of mechanisms, such as the induction of apoptosis, oxidative stress damage, and autophagy (22–24). Scott showed that HDACi blocked tumour cells in the G1/S phase and inhibited tumour growth (25). Vorinostat inhibits the growth of cancer cells by TGF‐β1 (26).