The Gene Expression Omnibus (GEO) database was used to retrieve expression profiling datasets using the search term ‘gastric cancer' (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2012). The criteria for dataset selection are as follows:

i. Datasets from high-throughput sequencing and single-channel microarray experiments with proper gene annotations were included.

For validation of differential expression results and survival analyses, a large cancer genomics cohort from TCGA, as well as at least 2 external GEO datasets from geographically distinct populations, was obtained.

The integrative approach considered RNA sequencing experiments due to their high sensitivity and dynamic range. Additionally, it included extensive data obtained from microarray experiments to produce reliable results. Figure 1 depicts the workflow used in this study.

Selected datasets were preprocessed and normalized individually, and differentially expressed genes (DEGs) were assessed using methods such as limma, DeSeq2, and limma-voom (Law et al., 2014; Love et al., 2014; Ritchie et al., 2015). All analyses were conducted in R using Bioconductor packages. For each dataset, genes were given raw p-values, log2 fold changes, and adjusted p-values. All genes were sorted primarily by adjusted p-values and then by absolute log fold change to indicate statistical significance.

The ‘RobustRankAggreg’ package in R was used for the meta-analysis of DEG results, as it is useful for combining gene expression results from multiple datasets. Here, the RRA algorithm integrated and ranked the genes consistently across all datasets, ensuring that variation across different platforms was captured (Kolde et al., 2012). The parameter was set to N = 20,000 to include all the significant genes and reduce the noise. Furthermore, Bonferroni correction, as shown in Equation 2.1, was applied, and the genes with scores below this threshold were selected for further analysis. α c o r r e c t e d = α m (2.1) α corrected is the adjusted significance threshold

α is the original threshold

m is the number of tests (genes).

The WGCNA algorithm was applied to all datasets that met the minimum sample size of 15 to identify significant gene co-expression networks and their correlation with the tumor phenotype (Langfelder and Horvath, 2008). Each dataset was preprocessed and normalized individually using appropriate methods, and the top 15,000 most variable genes were selected for the analysis using the WGCNA package in R. For each dataset, sample clustering was performed to remove outliers, ensuring high sample quality. Appropriate soft-threshold powers were chosen from the scale-free topology fit plot where R2 crosses 0.8. This was followed by signed weighted gene co-expression network construction and the detection of modules through hierarchical clustering with the dynamic tree-cut method. Furthermore, highly similar modules were merged to refine the assignment of modules. After analyzing module-trait relationships, modules with significant correlations with tumor samples were selected for further analysis to capture robust biological signals from all the datasets. From the tumor-correlated modules with module trait p-value less than 0.05 of each dataset, genes with module membership (MM) greater than 0.8 and gene significance (GS) more than 0.2 were chosen as hub genes. Unique gene symbols of tumor-associated hub genes from each dataset were combined to create a master hub gene list and a presence-absence logical matrix to indicate the gene occurrences across the selected datasets. After calculating the frequency of occurrence, genes were divided into three confidence categories: high-confidence genes, which appear in 7 or more datasets; medium-confidence genes, which appear in 3–6 datasets; and exploratory genes, which appear in at least 2 datasets. Genes present only in one dataset were considered low-support genes and were excluded from further analysis.

A common background universe (U) was defined as the set of unique genes common to the ranked genes from the RRA analysis and to the genes considered in the WGCNA analysis. All gene lists, including the universe list, were mapped to HGNC symbols using limma, HGNChelper and org.Hs.eg.db and duplicates were removed (Oh et al., 2020; Carlson et al., 2019). Later, for each category (High, medium, and exploratory), 2 × 2 contingency tables were constructed to count.

i. Genes present in both the RRA list and the WGCNA category (a)

Then, a one-sided Fisher’s exact test was applied as in Equation 2.2 to determine if the observed overlap is more than what would be expected by chance (Fisher, 1922). Along with the p-value, an odds ratio representing the strength of association between WGCNA category membership and presence in the RRA list was reported. Furthermore, a representation factor, calculated as the observed overlap divided by the expected overlap under random sampling, was calculated; values greater than 1 indicated over-representation. P-values across the three categories were adjusted by Benjamini-Hochberg false discovery rate (BH-FDR), and only those categories with significant over-representation (FDR<0.05) were chosen for downstream analyses (Benjamini and Hochberg, 1995). p greater = ∑ i = a min T , R R A R R A i U − R R A T − i U T (2.2) U = Background universe

R R A = set of RRA hits

T = WGCNA tier

a = Observed overlap between the WGCNA tier and RRA hits

Overlapping genes with statistical significance were analyzed for protein-protein interactions using the STRING (https://string-db.org/) and STITCH (http://stitch.embl.de/) databases (Kuhn et al., 2007; Szklarczyk et al., 2019). Interaction networks with a minimum score of 0.400 or medium confidence were imported to Cytoscape software (version 3.10.3) (Shannon et al., 2003). All 12 metrics, such as Maximal Clique Centrality (MCC), Degree, Betweenness, and Maximum Neighborhood Component (MNC) and others from the CytoHubba plugin were used to select the top 40 genes for prioritizing the hub genes. The genes identified consistently by at least 50% of the metrics were chosen as hub genes from Cytohubba (Chin et al., 2014). Simultaneously, the PPI network was evaluated for overlapping protein complexes using ClusterONEweb (https://paccanarolab.org/clusteroneweb/) with a minimum module size of 3, a minimum density of 0.3, a maximum overlap of 0.6, and a penalty of 2 for larger clusters (Baez et al., 2025). The results with p-values were downloaded, and the top modules were selected based on size, density, and p < 0.05. Furthermore, the common genes from both analyses were visualized using the VennDiagram (1.7.3) package in R for finalizing targets for enrichment analyses, such as Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) (Oliveros, 2007).

Enrichment analyses, including GO molecular function (MF), cellular component (CC), and biological process (BP), KEGG, and Reactome pathway enrichment analyses, were conducted using the clusterProfiler, ReactomePA packages, org.Hs.eg.db, AnnotationDbi, and enrichplot in R (Yu et al., 2012; Yu and He, 2016; Pagès et al., 2017; Guangchuang et al., 2023). The overlapping genes from the cluster and topology PPI analyses were considered as target genes, and the initial set of genes from the PPI network was chosen as the background. Gene symbols were mapped to their Entrez IDs using the bitr function of the Annotation DBI package. Enrichments were calculated using the mapped background, BH correction, and thresholds such as p. Adjust ≤0.05, q-value ≤0.20, minimum gene set size of 3, and a maximum of 200. Gene ratios, background ratios, and fold enrichment values were reported; dot plots were generated for top terms, and the terms with p.adjust ≤0.05, q-value ≤0.2, fold enrichment ≥2, and a minimum gene count of 3 were considered significant. The candidates that contributed to significant terms in functional enrichment analysis were prioritized for validation and survival analyses.

To validate the differential expression of the chosen candidates, datasets from distinct clinical populations were chosen. It included TCGA-STAD (RNA seq) (Colaprico et al., 2016), GSE54129 (GPL570; China), GSE13911 (GPL570; Italy) (D’Errico et al., 2009) and GSE66229 (GPL570; Asian Cancer Research Group (ACRG cohort)) (Oh et al., 2018). These cohorts were not utilized in the discovery data and were preprocessed with platform-specific methods and mapped to HGNC gene symbols using org.Hs.eg.db and hgu133plus2.DB before validating the dysregulation of the selected genes (Carlson et al., 2016). For the TCGA cohort, STAR counts were imported, and samples with ambiguous labels from the clinical metadata and genes with zero total counts were removed, and duplicates were collapsed after mapping to HGNC symbols. Later, the RNA sequencing data were analyzed for differential expression with DESeq2, which called on unshrunken log2 fold-changes with BH-FDR. The raw CEL files of the GEO microarray cohorts were RMA normalized with affy, and for genes with multiple probes, the probe with the largest interquartile range (IQR) across samples was retained (Irizarry et al., 2003). They were tested with limma and log2FC, standard errors (SE), p, and FDR values were recorded. Regardless of the significance of each study, meta-analysis combined the log2FC ± SE for each gene from each cohort using a random-effects model (REML) using metafor (Viechtbauer, 2010). This meta-analysis reported heterogeneity (Q, inconsistency index I², and τ²), two-sided p-values, pooled log2FC, and 95% confidence intervals. Furthermore, BH-FDR was applied to the selected genes, and the REML model was refitted by excluding one cohort at a time (Leave-one-out; LOO) to assess sensitivity. Genes with no changes in the pooled effect’s sign were labeled as stable under LOO. Genes were validated as differentially expressed if they had FDR ≤0.05 and a CI excluding 0 in the REML model, and had a consistent sign in the LOO analysis and p ≤ 0.05.

Three independent gastric cancer cohorts were used to evaluate the effect of the selected genes on overall survival (OS): TCGA-STAD (RNA sequencing), GSE62254 (microarray), and GSE15459 (microarray) (Ooi et al., 2009; Cristescu et al., 2015). OS was measured as the time from diagnosis to death, with event 1 indicating death and 0 indicating censoring. For TCGA, STAR counts were collapsed to HGNC symbols, normalized, and transformed using the standard size factor normalization of DESeq2 (Variance Stabilizing Transformation, VST), and RMA normalization (affy) was used for microarray datasets. Gene expression was standardized to a mean of 0 and a standard deviation of 1 before modeling in each cohort. Kaplan-Meier (KM) plots were created for visualization using the top 25% vs. the bottom 25% of expression (Q4 vs. Q1) within the cohorts to avoid dilution around the median, and log-rank p-values, along with n-at-risk tables, were added to plots for all genes.

Cox proportional-hazards models were applied to all candidate genes in each cohort. Univariate models related the OS to the standard expression of each gene, followed by multivariate models that adjusted for clinical covariates. For TCGA, the multivariate models were adjusted for age and sex, whereas for the GEO cohorts, age, sex, and stage were adjusted. Across all the models, hazard ratio (HR), 95% CI, the Wald p-value, a global Schoenfeld test p-value that checks the proportional hazards, and model discrimination using Harrell’s c-index were recorded. Moreover, BH-FDR correction was applied to each cohort in both univariate and multivariate analyses. For cohorts with stage data, gene × stage interactions were evaluated with stage Ⅰ as the reference and stages Ⅱ, Ⅲ, and Ⅳ as contrasts, followed by BH-FDR correction across genes. Interactions were considered significant only if they had FDR ≤0.05. The univariate Cox HRs for each gene across all three cohorts were combined in a meta-analysis using metafor, with a random-effects REML model. This analysis reported pooled HRs with 95% CI, p-values, and BH-FDR testing across the genes, along with heterogeneity values (Q, I2, and τ2). Genes were considered significant in the meta-analysis if they had an HR > 1 in at least 2 of the 3 cohorts, indicating that their higher expression was associated with worse overall survival. The REML model was refitted by dropping one cohort at a time (LOO), and genes were considered stable under LOO if there was no loss of significance and sign changes. The final candidate genes with survival significance were selected if the pooled HR > 1 and the pooled FDR ≤0.05, if the HR > 1 in ≥2/3 cohorts with no sign changes in LOO.

For the final candidate genes, a targeted literature survey was conducted to exclude the already validated therapeutic targets. Only those genes without prior in vitro or in silico targeting were chosen as final therapeutic candidates. TCMBank was selected as the compound library for docking, and the phytocompounds for the search term ‘stomach neoplasm’ were retrieved for ADMET profiling (Lv et al., 2023). The drug-likeness, pharmacochemical, and pharmacokinetic properties were evaluated using SWISS ADME (Daina et al., 2017). Toxicity profiling was performed using ProTox-III, and compounds with acceptable ADMET profiles were shortlisted; their 3D structures were downloaded from PubChem (Banerjee et al., 2024; Kim et al., 2025). The shortlisted phytocompounds were optimized and converted to the required PDBQT format using OpenBabel (O’Boyle et al., 2011).

i. When experimentally validated structures were available, 3D coordinates with high resolution were downloaded from the RCSB PDB (Berman, 2000). The coordinates were cleaned and protonated using MGL Tools.

Target preparations were performed using MGL tools, and primary binding pockets of the finalized receptors were predicted using the PrankWeb, which provided the grid center coordinates for further docking analysis (Polák et al., 2025). To maintain consistency across all receptor sets, no bound metal ions or cofactors were retained in any of the receptor files. AutoDock Vina was used to perform the docking analysis, where the grid of each target was centered on the Prankweb-predicted primary binding pocket (Trott and Olson, 2010). A uniform cubical box of 24 × 24 × 24 Å, exhaustiveness of 32, energy range of 5 kcal mol^-1 and 20 modes were used consistently across all docking runs, and the highest-scoring Vina pose of all the ligands was retained. The targets and the ligands were docked in triplicate with different random seeds to confirm that the binding affinities did not vary much, and the top-scoring ligand pose was visualized with the corresponding target proteins in Biovia Discovery Studio Visualizer (v21.1.0.20298) (BIOVIA, 2021).

A total of 31 publicly available datasets that matched the inclusion criteria (discussed in Section 2.1) were retrieved from the GEO website for the comparative transcriptomic analysis. Of these, eighteen datasets were submitted from single-channel microarray experiments, and thirteen were generated from high-throughput RNA sequencing experiments. Across the chosen datasets, 927 tumor samples and 919 matched adjacent normal stomach tissues were included, providing a total of 1846 samples for the transcriptomic analysis. Cohorts chosen included different geographic regions and multiple platforms; for instance, microarray studies were frequently generated on Affymetrix platforms, such as GPL570 and GPL96, whereas RNA sequencing studies were mostly obtained from multiple Illumina-based platforms. Sample sizes varied widely, ranging from 3 pairs to 230 pairs of tumor-normal samples. This ensured that the study covered different technical and individual study settings. More details, including the accession IDs and individual sample counts for each dataset, are provided in Table 1.

To validate the differential expression of candidate genes, four independent external cohorts were utilized that covered distinct clinical populations: TCGA-STAD, GSE54129, GSE13911, and GSE66229. The prognostic relevance of selected genes was evaluated across three independent cohorts with survival data: TCGA-STAD, GSE62254, and GSE15459.

Differential expression testing was conducted using appropriate methods across all 31 datasets after proper preprocessing based on their data types. Data normalization in RNA sequencing studies generally included VST or log2 transformations of FPKM/CPM (Fragments per kilobase of transcript per million mapped reads/counts per million) values. Microarray data were frequently normalized using quantile normalization and robust multiarray average (RMA). Before differential expression (DE) analysis, outliers were searched for, and two datasets (GSE84787 and GSE79973) underwent outlier removal. Across datasets, Limma and DeSeq2 were primarily used for DE analysis, and the results exhibited considerable heterogeneity and size variations. Table 2 provides the preprocessing steps, testing methods, and resulting gene counts of each dataset.

The RobustRankAggreg algorithm (RRA) was used to integrate results across studies by ranking genes by adjusted p-value and then by log2 fold change (log2FC) in each dataset, and then combining the 31 ranked gene lists. Then, a score (RRA score) was assigned to all the genes, followed by Bonferroni correction to control multiple testing, where α was 0.05, m was 20,000 genes, producing a Bonferroni-adjusted threshold of 2.5 × 10⁻⁶. Following the RRA meta-analysis and the Bonferroni correction, 254 genes with consistent differential expression across both RNA-seq and microarray datasets were identified with significance. The uncorrected RRA scores of all the genes are provided in Supplementary Table 1A, and the list of Bonferroni-significant 254 genes and their scores is provided in Supplementary Table 1B.

A total of 19 datasets met the criterion of a minimum sample size of 15 samples; each of these datasets was preprocessed and analyzed using the WGCNA R package. The pickSoftThreshold function was used to determine the optimal soft threshold power for each dataset. Signed networks were then constructed, modules were identified, and highly similar modules were merged. To identify tumor-correlated modules, module eigengenes (MEs) were correlated with the tumor-versus-normal trait, and genes with GS > 0.2 and MM > 0.8 were selected as hub genes from these modules. Table 3 summarizes the preprocessing methods and soft-threshold powers used, along with the number of hub genes from tumor-correlated modules in each dataset. A presence-absence matrix was constructed using the gene symbols of these hub genes to calculate their frequency of occurrence across datasets, and the genes were then categorized into three tiers: high confidence, medium-confidence and exploratory. The sizes of these three categories were defined at three levels through pre-processing:

i. Original sizes based on the presence-absence matrix: High/Medium/Exploratory- 655/3390/3273 (Supplementary Tables2A-C provides the dataset details and the frequencies of the genes from all three categories)

To determine whether the overlaps were different from what would be expected by chance, one-sided Fisher tests were conducted to compare each WGCNA category with the RRA gene set, and the summarized effect sizes, along with OR and RF, were recorded. Of the 653 genes from the high-confidence WGCNA category, 142 overlapped with the RRA list with an OR of 66.80, RF of 23.75, p-value of 1.75 × 10 (⁻¹⁶⁵), and FDR of 5.25 × 10 (⁻¹⁶⁵), thus exhibiting significant enrichment. The 3355 genes from the medium-confidence WGCNA category exhibited a significantly enriched overlap of 94 genes with the RRA set with an OR of 4.36, RF of 3.06, p-value of 1.08 × 10 (⁻²⁴), and FDR of 1.6 × 10 (⁻²⁴). In contrast, among the 3,229 genes from the exploratory subset of WGCNA, only 10 genes overlapped with the RRA set, with an OR of 0.31, an RF of 0.34, a p-value, and an FDR approximately equal to 1. The contingency counts and other statistical values, including the a, b, c, and d values for each category and the lower bounds of the CI for the odds ratio, are provided in Supplementary Table 2D. The results of Fisher’s test demonstrated that the exploratory category did not exhibit significant enrichment, whereas the overlap between the high-confidence and medium categories was statistically significant. Thus, a total of 236 genes (142 from high-confidence overlap and 94 from medium-confidence overlap) were selected for further analysis, and Supplementary Table 2E lists these genes along with their WGCNA category and RRA scores.

The STRING database was employed to retrieve protein-protein interactions for the 236 overlapping genes, considering a minimum interaction score of 0.4. As STRING lacked entries for Aquaporin-4 (AQP4), its interactions were supplemented using the STITCH database, and the combined network was then imported into Cytoscape. Figure 2 shows the combined PPI network of the 236 hub genes.

All twelve metrics of the Cytohubba plugin were applied to the network to analyze and rank the topology of the nodes, and the top 40 genes identified by each measure are depicted in Supplementary Figures 1A-L. The 41 genes identified by at least six of the twelve metrics (≥50% consensus) were selected as Cytohubba-derived hub genes, and their details are provided in Supplementary Table 3A. Simultaneously, the same network was used to identify densely connected protein complexes using the ClusterONE web. Ten significant clusters were identified with a minimum density of 0.3 and p-values <0.05, as shown in Supplementary Table 3B, and were then prioritized based on a minimum size of 7 genes. Four clusters (Clusters 1,2, 5, and 6) with 65 component genes were chosen as key clusters due to their high quality and low p-values. For robust selection of hub genes, the 35 genes that overlapped between Cytohubba hub genes and genes from key ClusterONE modules were selected for further analysis. Figure 3 shows the Venn diagram used to visualize the overlapping hub genes from the analyses above.

To understand the biological roles of the 35 hub genes, GO (BP, CC, and MF), KEGG, and Reactome pathway enrichment analyses were conducted. In the GO-BP enrichment, genes involved in extracellular matrix (ECM) organization, cell adhesion, angiogenesis, blood vessel morphogenesis, cell migration and other related processes in tissue remodeling and tumor progression were significantly enriched. Some of the recurrently enriched genes include multiple collagen family members and other genes such as MMP14, TGFBI, SERPINH1, THY1, FAP, SULF1, PDGFRB, SPARC, NID2, LOXL2, and THBS2. Figure 4A shows the dot plot representation of the top enriched GO biological process terms.

The cellular localization of the enriched genes was primarily in the extracellular matrix category, along with collagen-containing ECM, extracellular space, basement membrane, endoplasmic reticulum lumen, and collagen trimers. The key genes enriched in this category include VCAN, MMP14, SERPINH1, BGN, SPARC, TGFB1, CTSB, LOXL2, NID2 and collagen-family genes such as COL1A1, COL3A1, COL4A2, COL5A2, COL4A2 and COL18A1. Figure 4B displays the dot plot representing the enriched GO cellular component terms.

The molecular functions (MF) category mainly consisted of structural constituents of the ECM and binding activities to collagen and growth factors, both of which are related to matrix interactions and cell signaling. The key genes recurring in the top GO-MF terms were the collagen and ECM proteins like COL6A3, COL1A1, COL1A2, COL4A1, COL3A1, SPARC, NID2, and VCAN, contributing to the structural component function and genes like FAP, SERPINH1, PDGFBR, MMP14, LUM and CTSK contributing to growth factor and collagen binding functions. Figure 4C shows the dot plot of the top GO Molecular Function terms that were enriched in the analysis.

The KEGG pathway analyses revealed key pathways in cell-matrix interactions and cancer signaling, including the PI3K-Akt pathway, focal adhesion, ECM-receptor interactions, and protein digestion and absorption. Collagen family members like COL6A3, COL4A1, COL4A2, COL1A2, and COL1A1, along with genes such as THBS2, PDGFRB, LUM, and SPP1, were among the key genes in KEGG-enriched pathways. Figure 4D depicts the dot plot of the key pathways that were enriched in the KEGG analysis.

The Reactome pathway enrichment analysis revealed major pathways in the organization and degradation of ECM, the formation, assembly, biosynthesis, and degradation of collagen, and interactions, including integrin/non-integrin interactions and PDGF/RTK signaling pathways. Apart from major collagen genes such as COL1A2, COL1A1, COL3A1, COL4A2, COL5A2, COL18A1, and COL4A1, which are part of multiple pathways, other genes, such as BMP1, PLOD1, PDGFRB, SPP1, CTSK, and MMP14, were also members of the enriched ECM/growth factor pathways. Figure 4E displays the dot plot of the top Reactome pathways that were enriched.

The genes that were not enriched in any of the above analyses (p.adj ≤0.05) and the ones without at least 10 occurrences across the ECM/cancer-related significant terms and pathways were excluded from being considered as hub genes. Supplementary Tables 4A-E provides the results of the GO, KEGG and Reactome enrichment analyses. The final hub genes that were nominated as the final 29 potential therapeutic candidates were BMP1, COL1A2, LUM, SPP1, THBS2, VCAN, COL1A1, COL3A1, CTSK, MMP14, PDGFRB, SERPINH1, THY1, BGN, COL4A2, COL5A2, FAP, SPARC, TGFB1, COL12A1, COL18A1, COL4A1, CTSB, LOXL2, NID2, COL6A3, PLAU, PLOD1, and SULF1.

For validation of the differential expression of the 29 hub genes, four independent cohorts from distinct populations were used: RNA sequencing data from TCGA-STAD (412 tumors and 36 normal samples) and three microarray datasets from GPL570. GSE54129 included 111 tumors and 21 normal samples from China; GSE13911 included 38 tumors and 31 normal samples from Italy; and GSE66229 included 300 tumors and 100 normal samples from the Asian Cancer Research Group (ACRG). All 29 hub genes were detected in the validation cohorts and quantified. Most of the candidates satisfied the |log2FC| ≥ 1 and FDR ≤0.05 in individual cohorts, and their expression results are provided in Supplementary Tables 5A-D. The volcano plots of these genes from all four cohorts are given in Figures 5A–D. Regardless of the individual results, a meta-analysis of the 29 genes was conducted using the REML model; 22 genes were significant (BH-FDR ≤0.05) with 95% CI and showed consistent dysregulation across cohorts. ECM and stromal genes, such as COL1A1, BGN, SPP1, MMP14, and PLOD1, were consistently higher in tumors with positive pooled log2FC values. Heterogeneity (I²) values were often high because the study combined RNA sequencing and microarray data from distinct populations, yet the pooled effects showed consistent direction. LOO studies were conducted for sensitivity, and most genes were retained as they did not exhibit any sign flips or p-values ≥0.05, regardless of which cohort was removed. CTSK, PLAU and TGFB1 were excluded from downstream analyses due to insignificant REML results, and COL3A1, COL4A2, COL6A3, and LUM were excluded due to insignificant p values in the LOO analysis. Supplementary Table 5E displays the REML and LOO results of the 29 genes. The 22 genes out of the 29 candidates with validated differential expression include BGN, BMP1, COL12A1, COL18A1, COL1A1, COL1A2, COL4A1, COL5A2, CTSB, FAP, LOXL2, MMP14, NID2, PDGFRB, PLOD1, SERPINH1, SPARC, SPP1, SULF1, THBS2, THY1, and VCAN. Forest plots of these genes were generated with log2FC per study, with pooled REML at the bottom, and are provided in Supplementary Figures 2A-Z, 2a-c.

For analyzing the prognostic effects of the selected 22 genes, three independent cohorts were utilized: TCGA-STAD had 445 tumor samples with OS and covariates mentioned such as age (437/445), sex (445/445) and stage (77/445); GSE62254 had 300 tumors, GSE15459 had 192 tumors and covariates such as age, sex and stage were available for all the samples from both datasets. Gene expression for the 22 genes was standardized within each cohort, and KM plots were generated for visualization, with Q4 vs. Q1 for each gene. Supplementary figures 3A-V, 4A-V, and 5A-V show the KM plots of the 22 genes with the n-at-risk tables from all three cohorts.

Across all three cohorts, higher expression of the 22 genes was often associated with worse OS. In the univariate and multivariate Cox models adjusted for age and sex, TCGA data showed that many genes, including BGN, FAP, NID2, PDGFRB, THY1, VCAN and some collagen genes were significant with adjusted-HR >1 and FDR ≤0.05. Similarly, univariate and multivariate models adjusted for age, sex, and stage across both GSE62254 and GSE15459 showed independent associations between gene expression and adverse survival outcomes for many genes, including COL18A1, COL4A1, BGN, NID2, VCAN, and THBS2, among others. Supplementary Tables 6A-F provides the univariate and multivariate results of all three cohorts. The results from these cohorts indicated consistent associations between ECM-related gene dysregulation and poorer survival outcomes. The power of the interaction analyses among genes and stages in TCGA was restricted by the small number of cases with stage details (n = 77), and by the lack of significant interactions between any genes and stages with FDR ≤0.05. Within the GSE62254 cohort, some genes, like CTSB and SPP1, showed gene-level associations with multiple stages, but after BH correction, no genes had significant interactions. For the GSE15459 cohort, no interactions were significant after BH correction within genes or across genes. Thus, there was no consistent evidence that the pathological stages could change the prognostic significance of any of the genes.

REML analysis of the univariate results across all cohorts identified 17 of the 22 validated genes as consistent markers (FDR ≤0.05) whose higher expression was associated with worse OS. The percentage of variation across cohorts (I²) was low to moderate for most of these genes, and in the LOO analysis, there was no loss of significance or sign changes. Even though they exhibited consistent tumor-normal dysregulation, five genes (COL12A1, CTSB, PLOD1, SPP1, and SULF1) were excluded from the final candidate list as they did not satisfy the meta-analysis criteria and exhibited loss of prognostic significance. This discordance suggests that differential expression does not necessarily mean the gene would contribute to the prognostic significance. Supplementary Tables 6G-H provides the results of the REML and LOO analyses. Forest plots of all the 22 genes displaying the cohort-specific HRs with 95% CIs, pooled HRs, and other statistical values like I², and c-indexes are provided as Supplementary Figures 6A-F. The significant genes in the survival analysis were shown to be key players in ECM remodeling with stromal significance, including collagen genes, matrix-processing genes, and cancer-associated fibroblasts (CAFs) or stromal markers. Thus, the 17 final potential therapeutic candidates with significant prognostic associations (the pooled HR > 1, pooled FDR ≤0.05, HR > 1 in ≥2/3 cohorts with no sign changes in LOO) are BGN, BMP1, COL18A1, COL1A1, COL1A2, COL4A1, COL5A2, FAP, LOXL2, MMP14, NID2, PDGFRB, SERPINH1, SPARC, THBS2, THY1, and VCAN.

For the 17 chosen targets, a systematic literature search was conducted to identify prior in silico, in vitro, or animal model-based experiments targeting them in gastric cancer. It was found that there is limited data on reports of natural compound or drug testing for 4 of the selected genes in gastric cancer. As a result, those genes were retained for docking analysis, and the other 13 genes, which had been previously targeted, were excluded. The source of phytochemicals was TCMBank, where 58 compounds were retrieved for the search term ‘stomach neoplasm’. Figure 6 shows the interaction diagram of the 58 compounds obtained from TCMBank. SMILES strings for the 58 compounds were utilized for ADMET evaluation, and SwissADME was used for ADME evaluation to prioritize compounds that obeyed Lipinski’s rule of five, exhibited high gastrointestinal (GI) permeability, and showed favorable bioavailability scores. ProTox-III was used for toxicity profiling to retain compounds in toxicity classes IV and V (low acute toxicity). Only 14 compounds exhibited acceptable ADMET profiles, and their 3D structures were retrieved in SDF format from PubChem for docking. Supplementary Table 7 presents the ADMET profiles of the selected 14 compounds. The reference compound chosen was 5-fluorouracil (5-FU), a common small-molecule chemotherapeutic agent in the treatment of gastric cancer. These compounds were energy-minimized and converted to formats suitable for docking (PDBQT format) using OpenBabel.

After the targeted literature search, it was shown that most of the final therapeutic candidates had already been studied as targets in in silico, in vitro, or in vivo settings, and four targets (THBS2, VCAN, THY1, and NID2) had limited information on their use as targets. For THBS2, PDB entry 1YO8 containing the C-terminal domain was used as it was the highest-quality entry available in the RCSB PDB. This PDB model was directly cleaned and protonated using MGL tools. As experimental structures were not available, AlphaFold models corresponding to the Uniprot entries were used for VCAN (D6RGZ6), THY1 (P04216), and NID2 (Q14112). For these models, residues with confidence scores less than 50 (pLDDT<50) were removed using PyMOL, and missing loops were remodeled using SWISS-MODEL. Additionally, steric clashes were repaired using PRAS for NID2, and the energy minimization was performed in SwissPDBViewer for the THY1 model (Nnyigide et al., 2022; Guex and Peitsch, 1997). The primary binding pocket information of all the finalized structures was obtained using PrankWeb. Finally, after protein preparation, all structures were saved in PDBQT format for docking.

AutoDock Vina was used for docking simulations against each protein’s primary binding pocket. Additionally, 5-FU was docked to these sites as a reference compound, and each ligand-target pair was run in triplicate with different random seeds. The best affinities were consistent across the runs (variation less than 1 kcal mol^-1), and lead phytocompounds were chosen based on the lowest binding affinities. All selected compounds exhibited better binding scores than the reference drug 5-FU, and the phytocompound showed the highest affinity for VCAN, NID2, and THBS2. For VCAN, along with evoden, ursolic acid scored the best affinity of −9.5 kcal mol^-1 in all triplicates. Evoden showed high affinities of −10.4, −10.4, −10.3 kcal mol^-1 with NID2 and a binding affinity of −9.5 kcal mol^-1 with THBS2 across all the runs. Ursolic acid scored best with the THY1 model with −6.9 kcal mol^-1 across the runs, even though it is higher than the reference drug (−4.4 kcal mol^-1); it is comparatively low, probably due to the shallow ectodomain of the protein. Table 4 summarizes the docking and interaction results with binding affinities (triplicate mean ± standard deviation, SD) of the four targets. The interactions between the top-scoring ligand model and targets were visualized using Discovery Studio Visualizer, and the interactions were noted. Figures 7A–E depicts the interaction diagrams of all the target-ligand complexes.

Across these targets, the docked complexes exhibited non-classical interactions, including alkyl, pi-alkyl, and carbon-hydrogen bonds in addition to the standard donor-acceptor hydrogen bonds. Unfavorable bumps may be due to shallow geometry, imperfections in the modeled regions, or the absence of cofactors and metal ions, suggesting areas for further model refinement and experimental validation. However, triplicate runs produced reproducible results across the targets, where ursolic acid and evoden showed good binding affinities to the chosen targets.