This study analyzed publicly available clinical and transcriptomic data of five GC cohorts, as discovery datasets. First, The Cancer Genome Atlas (TCGA) gastric cancer cohort (n = 345) (Bass et al., 2014) with RNA-seq data was included (Table 1). Transcriptome data of the TCGA cohort were generated using HiSeq 2000 (Illumina) and included information on 20,508 genes. Clinical and genetic information were retrieved from the Firebrowse data portal site (http://firebrowse.org/). Next, four transcriptome datasets based on the microarray were included: Asian Cancer Research Group (ACRG, GSE66229 dataset, n = 297) (Cristescu et al., 2015), Yonsei University Severance Hospital (YUSH, GSE84437 dataset, n = 433) (Kim et al., 2019), Korea University Guro Hospital (KUGH, GSE26899 data set, n = 93) (Oh et al., 2018) and the National Cancer Centre of Singapore (NCCS, GSE15459 data set, n = 192) (Ooi et al., 2009) (Table 1). Microarray data were retrieved from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). Data from ACRG and NCCS cohorts were generated using Human Genome U133 Plus 2.0 (Affymetrix) which contained probes for 23,520 genes. In comparison, YUSH and KUGH cohorts were generated using HumanHT-12 V3.0 (Illumina) which contained probes for 25,124 genes. As external and independent validation datasets, we used additional tumor specimens collected from the Kosin University College of Medicine (KUCM, GSE26901 dataset, n = 109) and The University of Texas MD Anderson Cancer Center (MDACC, GSE28541 data set, n = 40), described in a previous study (Oh et al., 2018).

After collecting gene expression data, we used survival receiver operating characteristic (ROC) curves to calculate the optimal gene expression cutoff value for each gene stratifying patients into high and low-expression level groups. About the sample size of high/low expression groups, we considered a cutoff value to preserve at least 20% of the total sample size in these groups avoiding possible overfitting, groups with small sample sizes, and biased selection of cutoff value. If multiple cutoffs were retrieved for each gene, the cutoff value closer to the median was selected. Overall survival (OS) events were used in this calculation. Then, these groups were compared by univariate and multivariate Cox proportional hazard regression for each gene in each cohort. Possible confounding variables were included in the multivariate analysis depending on data availability. Thus, we adjusted prognosis prediction with the following potential confounders: age (as continuous), sex, Lauren classification, T stage, N and M status, primary tumor site, and molecular subtype. We considered that a gene is associated with prognosis if the area under the survival ROC curve (AUC) was >0.5, and a p-value < 0.05 for Cox regressions, in the five cohorts in an independent analysis for each cohort.

To improve the prognosis prediction of the gene signatures, we modeled gene expression-based scores using the genes identified as described in the previous section. The scores were calculated by the weighted sum of gene expression values, using the β coefficient of Cox regression, according to the next formula: S c o r e = G e n e 1 β 1 + G e n e 2 β 2 + G e n e 3 β 3 + . . . + G e n e i β i

The β value corresponds to pooled β coefficient calculated from the five cohorts by gene and estimated by a fixed- and random-effects model according to the heterogeneity between studies. Heterogeneity was assessed using Higgin’s I2 statistic and Cochran’s Q-test (Higgins and Thompson, 2002). A score was considered highly predictive when univariate HR > 2 (p < 0.05), multivariate HR > 1.5 (p < 0.05), AUC > 0.5, sensitivity > 0.5, and specificity > 0.5, in the five cohorts in an independent analysis for each cohort. Then, we selected the best gene expression-based score after considering the HR values and the number of genes constituting the panel. Finally, to homogeneously apply the score to any cohort, we normalized the score by dividing the score values by the median of the score calculated for each discovery dataset. Then, we established three groups according to the score values: low (score < 0.9), intermediate (0.9 ≤ score < 1.1), and high (score ≥ 1.1) groups.

We performed the Gene Set Enrichment Analysis (GSEA) using the Broad Institute software (http://software.broadinstitute.org/gsea/index.jsp) to find biological processes associated with the score. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome databases were included in our analysis. Enrichment scores were calculated based on Kolmogorov-Smirnov statistics tested for significance using 1,000 permutations. Additionally, Pearson correlation was used as a metric for ranking genes. A pathway was considered enriched when a nominal p-value and FDR q-value were <0.05.

The immune stromal cell signature of tumor samples was evaluated by using the following algorithms: CIBERSORT (Newman et al., 2015), EPIC (Racle and Gfeller, 2020), xCell (Aran et al., 2017), MCP-counter (Becht et al., 2016), and quanTIseq (Finotello et al., 2019); from the TIMER web tool (http://timer.cistrome.org/). Next, we explored the drug response in gastric cancer cells with different score values. The gene expression score was calculated in the following gastric cancer cell lines: AGS, FU97, GCIY, HGC-27, Hs-746T, IM-95, MKN45, MKN7, NCI-N87, NUGC-3, RERF-GC-1B, SNU-1, SNU-16, SNU-5, TGBC11TKB, and 23132-87. Gene expression data were retrieved from the Expression Atlas database (www.ebi.ac.UK/gxa/experiments/E-MTAB-2770/). Then, we correlate the scores with the inhibitory concentration 50 (IC50) for 448 drugs (https://www.cancerrxgene.org/). Spearman rho coefficients >0.25 or < −0.25 with p < 0.05 were considered to identify drug responses.

The OS was defined as the time from diagnosis to death for any cause; those alive or lost to follow-up were censored. Relapse-free survival (RFS) or disease-free survival (DFS) was defined as the time from the achievement of response to the first adverse event. That is relapse or death from any cause, whichever occurs first. Kaplan-Meier (KM) plots were performed to demonstrate the power of stratification. The Mann-Whitney test and Fisher’s exact test were used to compare clinical and biological characteristics between patients. All p-values were two-sided with a significance level of 0.05. All calculations were performed using R 4.1.1 (The CRAN project, www.r-project.org) software.

To identify survival-related genes in gastric cancer (GC), we analyzed the transcriptome data of five GC cohorts (Table 1). We dichotomized cohorts according to gene expression using the survival ROC curve. Then, for prognosis evaluation, high and low expression categories were established for each gene in each cohort (Figure 1). We identified 48 genes whose high expression was associated with short overall survival rates. Furthermore, 40 genes were consistently associated with disease-free and relapse-free survival in TCGA, ACRG, and KUGH cohorts (Figure 2A). On the other hand, low expression of three genes was associated with short OS rates, although none of these could predict DFS or RFS (Figure 2B); consequently, these genes were not considered in further analysis. By merging these results, we identified 11 genes associated with prognosis (Figure 2C). Supplementary Tables S1–S5 show the association between gene expression and overall survival for all genes included in the platforms for each cohort. Next, based on the 11 genes we found, we modeled 2,047 gene expression scores, which refer to all unique combinations among these 11 genes. To uniformize our score calculations, we pooled the β coefficients of the five cohorts for each gene (Supplementary Figures S1–S3). After applying Cox regression and AUC discrimination, 9 scores showed an improved power of survival prediction in the five cohorts (Figure 2D; Supplementary Figure S4). Next, we selected the most parsimonious score, based on seven genes, here referred to as GES7. The GES7 is composed of CCDC91, DYNC1I1, FAM83D, LBH, SLITRK5, WTIP, and NAP1L3. Expression levels of these genes and the comparisons with their normal counterparts are shown in Figure 2E. Also, Kaplan-Meier plots showing the stratification power of these genes are shown in Supplementary Figures S5–S8 (Log-rank, p < 0.05).

The GES7 was highly efficient in stratifying overall survival in GC patients (Figure 3A). Thus, GES7^high patients had inferior survival rates, with median OS of 19.9 vs. 69 months for TCGA, 27.5 months vs. not reached for ACRG, 44 months vs. not reached for YUSH, 22.1 months vs. not reached for KUGH, and 20.3 months vs. not reached for NCCS. In all cases, the AUC values of the GES7^high patients were higher than 0.5 (Figure 3B). Consistent with these results, we found that the GES7 score covariate with pointwise estimates of the overall survival hazard ratios (HRs) (Figure 3C). Also, GES7^high patients were more likely to relapse according to the disease- and relapse-free survival analysis (Figure 3D). Furthermore, after including potential confounding variables in our Cox regression analysis, GES7 remained as a relevant prognostic marker with HRs of 1.75 (p = 0.0061), 1.88 (p = 0.0014), 2.01 (p < 0.0001), 3.8 (p = 0.001), and 3.59 (p < 0.0001), for TCGA, ACRG, YUSH, KUGH, and NCCS, respectively (Figures 4A–E). Although GES7^low patients included more MSI cases and GES7^high patients showed more CIN, GS, and EMT cases (Supplementary Tables S6, S7), our multivariate regression demonstrates that the molecular subtypes are not a confounder variable. Interestingly, we found that GES7^high was associated with an invasive phenotype as described for the NCCS cohort (Supplementary Table 10). No other clinical variable was consistently associated with the GES7 (Supplementary Tables S6–S10). Finally, we validated our GES7 in two independent cohorts, KUCM and MDACC. We found that GES7^high patients had inferior survival rates, with median OS of 17.5 months vs. not reached (log-rank p < 0.0001), and 9.3 vs. 24.6 months (log-rank p = 0.0013), respectively (Figure 4F). Also, ROC curve analysis demonstrated that the area under the curve values remained consistently high, exceeding 0.6, for up to 5 years of follow-up in both validation cohorts (Figure 4G). Therefore, the GES7 gene signature proved to be highly efficient in stratifying overall survival in gastric cancer patients across multiple cohorts.

Next, we combined the traditional anatomical staging system with the molecular signature of GES7. For this purpose, we normalized the GES7 values across the five cohorts to find widely applicable cutoffs despite the techniques used to measure gene expression (RNA-seq or microarray) (Supplementary Figure S9). Next, we established three categories, GES7 low, intermediate, and high, applying cut-off points of 0.9 and 1.1 for all cohorts. By considering these categories, we proposed a risk classification system that integrated the AJCC TNM staging and the GES7 (Supplementary Table 11). As seen in Supplementary Figure 10A, the three categories of GES7 identified patients with different survival rates in the ACRG cohort. As expected, the AJCC TNM staging (7th edition) also efficiently distributed patients into other prognostic groups (Supplementary Figure 10B). After integrating both systems, we observed a reclassification of patients in a best-fitted risk group (Supplementary Figures 10C–D). Considering the IIIA category as intermediate for prognosis (median OS_IIIA = 68.1 months) and IIIB and IV categories as adverse (median OS_IIIB = 28.2 months and median OS_IV = 17.4 months), we had that the GES7 restructured these groups by the formation of an adverse-risk group constituting of patients from IIIA, IIIB and IV (median OS_adverse = 25.2 months), leaving the remaining patients of IV group forming the very-adverse risk group (median OS_adverse = 12.7 months) (Supplementary Figures 10D, E). These results were validated in the TCGA cohort, where the adverse and very-adverse risk groups differed from the intermediate and favorable groups compared to the AJCC TNM system alone (Supplementary Figures 11A–E). Therefore, the integration of our gene signature with the classic anatomical staging system improves prognostication accuracy and provides clinicians with a more comprehensive assessment of the patient’s disease.

To explore GES7-related biological pathways, we performed a GSEA with the transcriptome from the five discovery datasets. Figures 5A, B show that 579 pathways were consistently enriched in GES7 high samples, while 12 were for GES7^low. The top enriched pathway was the transforming growth factor beta (TGF-β) activation (Figure 5C). In line with this result, pathways related to elastic fiber were enriched, indicating a motility phenotype. On the other hand, antigen processing and presentation-related pathways were negatively enriched, indicating low activity of the innate immune system in these tumors (Figure 5D). Then, we tested the hypothesis of non-tumor cell infiltrations being associated with GES7. We found that GES7^high tumors are enriched with stromal cells like cancer-associated fibroblasts and have a reduced presence of innate immune cells (Figure 6A). Finally, we found that GC cell lines with high GES7 respond well to in vitro treatment with PAK_5339. This drug inhibits the p21-activated kinases 1 (PAK1) and 2 (PAK2), Elesclomol (anti-HSP90), Pilaralisib (anti-PI3K), UNC0628 (anti-G9a and GLP methyltransferases), and NSC-207895 (anti-MDM4) (Figure 6B), that could be valuable treatment strategies for GES7^high patients.

Previous gene expression-based prognostic indexes have been proposed as predictors of survival in GC. Here, we compared the GES7 with 30 prognostic scores published between 2013 and 2023 (Supplementary Table 12). After retrieving the gene list of each score, we calculated them using the weighted sum of gene expression values, the base strategy of all the studies. To homogenize the comparison between scores, we divided the cohorts into tertiles, and considered tertile 3 as the high score group, while tertile 1 as the low score group (Figure 7A). After applying the scores to our discovery datasets by univariate Cox regression analysis (Figure 7B), we found that 9 of 31 indexes were able to discriminate between patients of different risk groups. Finally, when we added the evaluation of our validation datasets, only two scores, within which is the GES7, were shown as the robust indexes (Figure 7C; Supplementary Figure 12).