The whole steps performed in the analysis were displayed in the Supplementary Figure S1 The transcription profile and clinical characteristics of gastric cancer were extracted from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and the National Center for Biotechnology Information Gene Expression Omnibus (GEO) databases (https://www.ncbi.nlm.nih.gov/geo/). The GEO datasets that contain more than 200 human samples and have complete expression data and clinical information were selected as independent test cohorts. Finally, the GSE84433 as well as the GSE84437 which is a SuperSeries that composed of the GSE84426 and GSE84433 Subseries were chosen. The TCGA-STAD project containing 375 tumor and 32 normal tissues, GSE84437 containing 433 human samples, and GSE84433 containing 357 human tissues were selected for further analysis (Yoon et al., 2020). The transcription profile was converted into fragments per kilobase million (FPKM), and the batch effect of the data in GSE84437 was eliminated. The clinical information of TCGA and GEO datasets were presented in Supplementary File S2. The transcriptomic data of TCGA-STAD was downloaded using the “TCGAbiolinks” R package. The expression difference between tumor and normal tissues was identified using the “limma” algorithm and presented with a heatmap and a volcano plot (Ritchie et al., 2015). The log2 fold-change (log2FC) > 1 and false discovery rate (FDR) < 0.05 were set as the cut-off to screen the differently expressed genes (DEGs). ROS-related genes (Supplementary File S1) were extracted from the Molecular Signatures Database (MSigDB) database (http://www.gsea-msigdb.org/gsea/index.jsp).

A protein-protein interaction (PPI) network of ROS-related genes was constructed using the STRING database (http://string-db.org/), revealing their putative connections. Then, the correlation between these genes and their mRNA expression was further explored by R packages “igraph” and “reshape2 (Csardi and Nepusz, 2006)”. Red lines show positive correlations, whereas blue lines indicate negative correlations.

First, the consensus clustering was conducted with the expression of prognostic ROS-related genes from the “limma” as well as “ConsensusClusterPlus” R packages (Wilkerson and Hayes, 2010; Ritchie et al., 2015). Through consistent cluster analysis, we obtained a cluster Max K value = 9 and cluster Num = 2. Then, K-M analysis of ROS-related clusters was conducted using the “survival” and “survminer” R packages (Rich et al., 2010). Additionally, a heatmap was applied to show the relationship between clusters and age, sex, tumor grade, clinical and T, N, and M stages, as well as the expression profiles of these genes in different clusters, demonstrating the success of this clustering. Sequentially, the immune cell infiltrations between different clusters were evaluated utilizing the “MCPcounter” R package, including B linage, CD8⁺ T cell, cytotoxic lymphocyte, fibroblasts, T cell, NK cell, monocytic linage, myeloid dendritic cell, and neutrophils (Becht et al., 2016).

TTCGA-STAD was used as the training set. The GSE84437 and GSE84433 were used as the validation sets. DEGs were identified in ROS-related clusters. Then, least absolute shrinkage and selection operator (LASSO) Cox regression was applied with the DEGs in the TCGA-STAD project to obtain a ROS-related signature (Tibshirani, 1997). The risk score was derived using the formula = ∑ i = 1 n C o e f i * E x p r i . Univariable Cox regression analyses of genes were conducted to assess the prognostic value of each single gene, and the mutation profiles of these genes were presented using a waterfall plot. The “maftools” R package was used to analyze and visualize the mutation including the missense mutation, non-sense mutation, frame-shift mutation, and multi-hit mutation. Patients were grouped into low- or high-risk subgroups according to the median score of the entire GC. Scatter plots were visualized to establish the correlation between survival status and risk score. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis were applied to assess the capability of the risk score to distinguish low- and high-risk patients (Ringnér, 2008; Cieslak et al., 2020). The receiver operating characteristic (ROC) curves were conducted to evaluate the predictive ability of the signature by “SurvivalROC” R package (Mandrekar, 2010). K-M survival analyses were conducted between different risk GC groups to verify differences in survival. Besides, three other signatures were used to conduct a comparation analysis with the ROS-related signature, comparing the C-index (Liu et al., 2022; Mak et al., 2022; Zhou et al., 2022).

Uni- and multivariable Cox regression analyses were performed on TCGA-STAD and GSE84437 datasets to identify independent prognostic factors. Sequentially, the relationship between the risk score and clinical characteristics was presented in the heatmap. Then we constructed a prognostic nomogram using the “rms” R package (Iasonos et al., 2008). Calibration plots were used to assess the discriminative power of the nomogram (Van Calster et al., 2016). The accuracy of this nomogram was verified using ROC and decision curves (Mandrekar, 2010; Kerr et al., 2016).

The limma algorithm was used to reveal the distinct risk of GC in different subgroups, including age, sex, tumor grade, clinical and T, N, and M stages, as well as immune subtypes (Ritchie et al., 2015). Moreover, Survival differences in the low- and high-risk subgroups were established using the K-M analysis. The ImmuneScore, StromalScore, and EstimateScore were assessed by “estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE)” algorithm, as well as immune cell infiltrations were evaluated using the “tumor immune estimation resource (TIMER),” “cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT),” “CIBERSORT-absolute mode (ABS),” “QUANTISEQ,” “microenvironment cell populations-counter (MCPCOUNTER),” “XCELL,” and “estimating the proportion of immune and cancer cells (EPIC)” algorithms (Becht et al., 2016; Aran et al., 2017; Li et al., 2017; Chen et al., 2018; Plattner et al., 2020). In addition, single-sample gene set enrichment analysis (ssGSEA) and immune function analyses have been conducted for different risk patients (Zhu L. et al., 2021).

The expression of immune-related genes was assessed using the “limma” algorithm. The tumor mutation burden (TMB), microsatellite instability (MSI), and tumor immune dysfunction and exclusion (TIDE) scores were then calculated. TMB score was conducted by R package “maftools” (Mayakonda et al., 2018), MSI score was obtained from previous research (Bonneville et al., 2017), and TIDE score was conducted by online database (http://tide.dfci.harvard.edu/). Finally, the immunophenoscore (IPS) of the different risks was evaluated in different subgroups. IPS refers to four main parts (effector cells, immunosuppressive cells, MHC molecules, and immunomodulators) determining the immunogenicity, and is calculated without bias using machine learning methods. The IPS of STAD patients were downloaded from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home).

To further explore the relevant factors of the ROS-related signature, we established the risk score in distinct wild and mutated GC types, including ARID1A, CSMD3, FAT4, FLG, LRP1B, MUC16, SYNE1, TP53, and TTN. Afterward, the relationship between signature and 5-fuorouracil, gemcitabine, cytarabine, dasatinib, etoposide, GSK690693, masitinib, and tipifarnib were evaluated according to inhibitory concentration (IC50) by the “pRRophetic” R package (Geeleher et al., 2014) based on the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to evaluate the functions of DEGs in low- and high-risk GC (Kanehisa and Goto, 2000; The Gene Ontology Consortium, 2019). Gene set enrichment analysis (GSEA) was applied between different risk groups. The threshold of GSEA analysis was log2FC > 1 and p < 0.05, and the top five enriched pathways were identified (Subramanian et al., 2005). The “c2. cp.kegg.v7.4. symbols. gmt” file downloaded from the GSEA database (https://www.gseamsigdb.org/gsea/index.jsp). Moreover, gene set variation analysis (GSVA) was performed with reference to the methods of previous research, and the enriched functions of different risk groups were visualized using a heat map (Hänzelmann et al., 2013). The enrichment analyses were performed by the “limma” and “GSVA” R packages (Yu et al., 2012; Ritchie et al., 2015).

HGC-27 and GES-1 cells were obtained from the Cell Bank of Shanghai Institute of Biochemistry and Cell Biology (Shanghai, China). The cell lines were cultured in DMEM or RPMI-1640 supplemented with 10% fetal bovine serum (FBS), 100 μg/mL streptomycin, and 100 U/mL penicillin (Gibco).

All samples were obtained with the approval of the Ethics Committee of Ningbo First Hospital. A total of 4 GC patients who signed the informed consents were recruited from the Ningbo First Hospital. The tumor samples and normal tissues were collected from the surgical specimen of the patients.

Total RNA was extracted using TRIzol reagent (Invitrogen) and reverse-transcribed into cDNA templates. β-Actin was used as an endogenous reference. The primer sequences used for amplification are listed in Supplementary File S3.

All statistical analyses were performed using the R software (version 4.1.3) (Supplementary File S4). Group comparisons were performed using the Wilcoxon test or Student’s t-test. Correlations between two continuous variables were assessed using Spearman’s correlation analysis. The Kruskal-Wallis test was used to compare the three groups. Adjusted p-value was calculated using Benjamini-Hochberg FDR. Statistical significance was set at p < 0.05.

First, 87 ROS-related genes were obtained from the GSEA database. Using the “limma” package, 31 differentially expressed ROS-related genes were identified (Figures 1A, B). The PPI network was constructed to identify ROS-related DEG interactions (Figure 1C). Meanwhile, Figure 1D shows the correlation network of DEGs, in which different colors indicate different correlation coefficients.

The effect of ROS-related genes on GC samples was assessed using consensus clustering analysis. When k = 2, the cumulative distribution curve was most horizontal in the middle section and the heatmap of clustering showed a clear edge. Therefore, GC samples in TCGA-STAD were grouped into two clusters based on ROS gene expression (Figure 2A). The overall survival (OS) rates of the C1 and C2 clusters were significantly different, and the prognosis of the C1 cluster was worse than that of the C2 cluster (Figure 2B). Besides, the K-M survival curve of progression free survival (PFS) also demonstrated the poorer outcome of C1 cluster (Supplementary Figure S2). Figure 2C shows the gene expression profiles between the C1 and C2 clusters. We also analyzed differences in immune cell infiltration between the two clusters and found more immune cell infiltration in the C1 cluster, including B lineage, CD8 T cells, cytotoxic lymphocytes, T cells, NK cells, neutrophils, fibroblasts, monocytic lineage, and myeloid dendritic cells, than in the C2 cluster (Figure 3).

A 4-gene prognosis model including GPX3, DUSP1, NOS3, and TCIRG1 was constructed using LASSO Cox regression (Figures 4A, B). All of the signature genes were identified prognostic genes by univariable Cox analysis, including GPX3, DUSP1, NOS3, and TCIRG1 (Figure 4C). We also analyzed mutations in these genes in GC (Figure 4D). GC patients were divided into high- and low-risk groups according to the median risk score in the training cohort (Figures 4E, F) and test cohort (Figures 4G, H). Patients in the high-risk group had a higher risk of death and poorer prognosis than those in the low-risk group. As shown in Figures 5A, B, PCA demonstrated clear distinctions between the two risk groups.

Low-risk patients were linked to better prognoses based on survival curves (Figure 5C) in the TCGA cohort. Similarly, in the test cohort, patients with lower-risk scores showed better survival (Figure 5D and Supplementary Figure S3). Additionally, the area under the curve (AUC) value indicated that our model had a good predictive power (Figures 5E, F). According to Cox regression analyses, the risk score was regarded as an independent prognostic factor for the training and testing cohorts (Figures 6A–D). The landscape of the four signature genes in the training group is shown in the heatmap (Figure 6E). The comparation analysis indicated that the stability of our signature is better than signatures of Mak et al., Zhou et al., and Liu et al. (C-index, 0.614 to 0.568, 0.597, and 0.607) (Supplementary Figure S4). To explore the genomic alterations in low- and high-risk GC, gene mutation analysis was conducted. Result showed that the TMB of low-risk GC is significantly higher than high-risk GC. Besides, some vital cancer-related genes such as TTN, TP53, MUC16, and ARID1A more frequently mutated in low-risk GC (Supplementary Figure S5).

A nomogram was developed to predict survival rates in GC (Figure 7A). The predicted results were in good agreement with the actual results according to the calibration plots (Figure 7B). Moreover, the nomogram was regarded as an independent prognostic factor (Figures 7C,D). The AUC of the nomogram was 0.747 (Figure 7E).

Next, the correlation between the risk scores and clinical features was studied. Different subgroups, including age, grade, T and TNM stage, as well as immune subtype, had significantly distinct risk scores (Figure 8). In addition, subgroup analysis showed differences in survival between different risk groups in different GC subgroups, including age >65 years, female, male, M0, age ≤65 years, Grade III, stage III-IV, N1-3 and T3+4, further verifying the reliability of this model (Figure 9).

The immune microenvironment plays a critical role in the development of GC. In our study, we found that patients in the high-risk group had higher immune, stromal, and ESTIMATE scores (Figure 10A). Next, the differences in immune cell infiltration between the different risk groups were assessed (Figure 10B). Infiltration of Monocytes, Macrophages M2, Mast cells resting, naïve B cells, and eosinophils was more abundant in the high-risk group, while infiltration of CD4 memory activated T cells, resting NK cells resting and Macrophages M1 was more abundant in the low-risk group (Figure 10C). In addition, APC co-stimulation, CCR, and Type II IFN responses were usually more significant in the high-risk group (Figure 10D).

Immune checkpoint molecules are also involved in the occurrence and development of GC. We found that the expression of most checkpoint genes was higher in the high-risk group (Figures 11A, B). TIDE, TMB, and MSI were used to predict the response to tumor immunotherapy. As shown in Figure 11C, the TIDE, exclusion, and dysfunction scores were all higher in the high-risk group than in the low-risk group. Moreover, the TMB and MSI scores were higher in the low-risk group (Figures 11D, E). Furthermore, we found that patients in the low-risk group responded better to immunotherapy (Figure 11F).

Gene mutations are crucial factors in tumor development. In the current study, the correlation between prognostic models and genetic mutations was analyzed. We found that the risk score was lower in the mutation group than in the wild-type group (Figure 12). Next, we predicted the chemoresponses of the subgroups to common chemotherapeutics (Figure 13). The results revealed that high-risk patients were more sensitive to 5-Fluorouracil, Gemcitabine, cytarabine, dasatinib, etoposide, GSK690693, masitinib, and tipipifarnib, suggesting that these patients could benefit from these chemotherapeutic agents.

To further explore the functional enrichment of the signature, we performed an enrichment analysis of DEGs between the low- and high-risk groups. GO enrichment indicated that DEGs were enriched in “muscle system process,” “collagen−containing extracellular matrix,” and “extracellular matrix structural constituent” (Figures 14A, C). KEGG enrichment analysis showed significant enrichment of the “PI3K-AKT” and “cGMP-PKG” signaling pathways (Figures 14B, D). Additionally, from the heatmap of GSVA, significant difference of enriched functions between different risk groups were presented (Figure 14E).

As shown in Figure 15A, the expression of GPX3, DUSP1, NOS3, and TCIRG1 was explored using the GEPIA2 site (http://gepia2.cancer-pku.cn/#index). Low expression of GPX3 and DUSP1 was observed in GC tissues as compared to that in normal tissues, whereas NOS3 and TCIRG1 were highly expressed in GC tissues. RT-qPCR was performed to verify the expression levels of these genes. The results demonstrated that NOS3 and TCIRG1 were highly expressed in tumor cells and samples compared to normal cells and tissues, while GPX3 and DUSP1 were expressed at low levels in tumor cells and samples (Figures 15B, C).