Initially, we identified differentially expressed proteins between groups with different pathological metastasis statuses. Simultaneously, we conducted a prognostic analysis across all proteins to identify those associated with prognosis. The intersection of these two protein sets yielded the prognosis-related differentially expressed proteins, from which we selected the top 50 for further analysis. Next, we determined the intersection between prognosis-related proteins and extracellular matrix-related proteins, thereby identifying the prognosis-related extracellular matrix proteins. Using the top 50 prognosis-associated differentially expressed proteins and the prognosis-associated extracellular matrix proteins, we constructed a protein interaction network. This network was visualized using Cytoscape (28) (v3.9.1), and its parameters were subsequently analyzed with the cytoHubba plugin. From the network, we extracted the top ten genes based on Degree, Closeness, and Betweenness centrality measures. The intersection of these genes with the top 50 prognosis-related differentially expressed proteins revealed the core hub proteins. Finally, we generated box plots and prognostic curves to illustrate our findings.

Constructing a gene regulatory network involved several key steps. First, potential target genes for each transcription factor were identified through co-expression analysis, using a gene expression matrix and algorithms like GENIE3/GRNBOOST to generate co-expression modules. Subsequently, DNA motif analysis was performed to pinpoint potentially direct targets, which requires defining cellular states and their respective regulators. Network activity was then analyzed by scoring regulators within each cell, calculating the area under the curve (AUC), and transforming this activity into a two-dimensional matrix. Finally, the Python package pysceinc was employed to identify metastatic intergroup-specific transcription factors and construct the gene regulatory network.

Following the arithmetic process, the cell’s transcription factor activity data were stored in the pyscenic_output.rds file, from which the area under the curve (AUC) value for each transcription factor was derived. These AUC values were then used to calculate the Connection Specificity Index (CSI) for each transcription factor. Based on the CSI, transcription factors were grouped into modules and visualized with a heatmap displaying CSI values across modules. Furthermore, UMAP and time series plots illustrated the mean CSI value for each module’s transcription factors. Subsequently, the Regulon Specificity Score (RSS) was calculated based on metastasis grouping, visualized with an RSS ranking diagram. UMAP plots then displayed the AUC of the top 10 regulons for both metastasis and non-metastasis groups. Finally, functional enrichment analysis of these top 10 regulons was performed using the clusterProfiler R package, with enriched pathways (p-value < 0.05 and q-value < 0.2) presented in a bubble plot.

The R package CellChat (34) (v 1.6.1) was further utilized for communication analysis of all cells. Single-cell data were first stratified into metastasis and non-metastasis groups, and CellChat objects were constructed for each group using the createCellChat() function with CellChatDB.human as the ligand-receptor database. Differentially expressed genes and ligand-receptor interactions were identified using the identifyOverExpressedGenes() and identifyOverExpressedInteractions() functions, respectively. Gene expression data were then projected onto a protein-protein interaction network using projectData(). Communication probabilities were calculated, and cell-cell communication networks were extrapolated with computeCommunProb(). A minimum cell threshold of 10 (min.cells = 10) was applied using filterCommunication(). Signaling pathway-level communication was inferred with computeCommunProbPathway(), and finally, the resulting cell-communication networks were integrated using aggregateNet().

After obtaining the communication information, enriched pathways between groups were compared using Fisher’s exact test to identify significant differences (p < 0.05). These differences were visualized using box plots. Subsequently, the top four significantly enriched pathways for each group were selected for circular plot representation. Moreover, network centrality indices were calculated for each cell population to determine their respective roles, and the corresponding importance was illustrated using a heatmap.

Enrichment scores for metastasis-related tumor cell subtypes were calculated using TCGA-STAD data, and samples were categorized into high and low-score groups based on the median enrichment score. Subsequently, samples were similarly divided into high and low-expression groups based on median hub protein expression. Drug response predictions were generated using the oncoPredict (35)(v0.2) R package, leveraging the CTRP2 and PRISM drug response databases. CTRP2 data were included with the oncoPredict package, while PRISM data were downloaded from the DepMap Portal (https://depmap.org/portal/). Gastric cancer cell data from both databases were processed, and drugs common to both were used for IC₅₀ prediction. Wilcoxon tests were performed to compare IC₅₀ values between high and low enrichment scores and expression groups. Drugs exhibiting consistent response trends across both databases and groupings, with statistically significant differences (p < 0.05), were visualized using box plots. These plots illustrated IC₅₀ differences between high and low groups, as well as correlations with both enrichment score and hub protein expression. Finally, Spearman correlations between hub protein expression levels were visualized using scatter plots.

In addition, the Drug Gene Interaction Database (DGIdb) (36)was queried to identify gene pairs interacting with differentially expressed drugs. After selecting approved interactions, expression matrices for the interacting genes were extracted. Spearman correlations were calculated between these genes and drug IC₅₀ values, enrichment scores, and hub protein expression levels. Bubble plots were then generated to visualize gene interactions within enrichment score subgroups, highlighting expression differences between hub protein expression subgroups.

Database analysis identified 1303 extracellular matrix proteins. Using histology data (PMC9522856), differential expression analysis across pM groups yielded 74 differentially expressed extracellular matrix proteins. Upregulated proteins (37 total) were defined as those significantly more abundant in the non-metastatic group(M0), while downregulated proteins were significantly less abundant in M0 (Figure 1A; Supplementary Table 1). GO and functional enrichment analyses of corresponding genes revealed distinct pathways for up- and downregulated proteins. Upregulated proteins were enriched in pathways such as negative regulation of protein localization to the plasma membrane, integrin-mediated signaling, Rap1 signaling, and natural killer cell-mediated cytotoxicity. Downregulated proteins were enriched in external encapsulating structural organization, cell-cell junction assembly, cell adhesion assembly, cell adhesion molecules, and tight junctions (Figure 1B). Hallmark enrichment scores, calculated from the protein abundance matrix (Supplementary Figure 1A), identified seven significantly enriched pathways differing across pM groups (Figures 1C, D). Prognostic analysis revealed that higher scores for HALLMARK_CHOLESTEROL_HOMEOSTASIS, HALLMARK_ESTROGEN_RESPONSE_EARLY, and HALLMARK_ANDROGEN_RESPONSE were associated with worse prognosis, while higher scores for HALLMARK_SPERMATOGENESIS and HALLMARK_G2M_CHECKPOINT correlated with better prognosis. HALLMARK_TGF_BETA_SIGNALING and HALLMARK_MITOTIC_SPINDLE showed no significant prognostic impact (Figure 1E). Differential expression analysis across pM groups identified 1535 proteins, and prognostic analysis across all proteins identified 734 prognostic-related proteins. The intersection of these sets yielded 282 prognostic-related differentially expressed proteins, from which the top 50 were selected for further analysis (Supplementary Table 2). Intersecting these with extracellular matrix-related proteins yielded 77 prognosis-related extracellular matrix proteins. A protein-protein interaction network, constructed using the top 50 prognostic-related differentially expressed proteins and the 77 prognosis-related extracellular matrix proteins, comprised 150 interactions (Supplementary Table 3). Analysis of degree, closeness, and betweenness centrality within this network identified four core hub proteins (Figures 1F–I). Further analysis demonstrated that RPS29 was significantly more abundant in M1 than M0, while IQGAP1, SH3GL1, and SCGN were more abundant in M0. Furthermore, high RPS29 abundance correlated with a worse prognosis, while low IQGAP1, SH3GL1, and SCGN abundance correlated with a poorer prognosis (Figures 1J, K).

Integrating two datasets, we analyzed four primary foci samples, two with metastases and two without. Single-cell sequencing identified eight initial cellular isoforms (Supplementary Table 4). Comparing cellular proportions between metastasis and non-metastasis samples revealed minimal differences in immune-like cell proportions, while epithelial cells, endothelial cells, and fibroblasts constituted larger proportions in non-metastasis samples (Figures 2A–C; Supplementary Figures 1B–E). Hub protein enrichment scores, calculated for each cell type and categorized as high or low, were significantly higher in T/NK cells and slightly lower in Myeloid cells, with a more even distribution across other cell types (Figure 2D). UMAP plots based on these score subgroups demonstrated higher T/NK and B cell scores in metastasis-free samples (Figure 2E), a finding confirmed by differential enrichment score analysis (Figure 2F). Given the observed heterogeneity in T/NK cells, myeloid cells, epithelial cells, and fibroblasts, we further identified subtypes within these populations (Supplementary Figures 2A–E; Supplementary Table 5). Examining the expression of four hub proteins (Supplementary Figures 2F–I), we found RPS29 highly expressed primarily in T/NK, B, and epithelial cells, and IQGAP1 highly expressed mainly in T/NK cells and some myeloid cells, with comparatively low expression of the remaining two genes across all cell types (Figure 2G).

T/NK cells were categorized into four subtypes, revealing a significant increase in naive CD4+ T cells within the metastatic group compared to the non-metastatic group, while the remaining subtypes were more prevalent in the non-metastatic group. Myeloid cells were classified into five subtypes, with Neutrophils1 and Monocytes primarily concentrated in the metastatic group, and Neutrophils2 mainly distributed in the non-metastatic group. Epithelial cells were identified as three subtypes, mainly distributed in the non-metastatic subgroup, particularly Proliferative cells1, which expressed both the tumor cell marker EPCAM and the proliferative cell marker MKI67. Fibroblasts, classified into two subtypes, were primarily located in the non-metastatic group (Figure 2H). Analysis of hub protein enrichment scores across cell types in metastatic and non-metastatic groups indicated higher scores in the non-metastatic groups for all cell types except B cells and Naive CD4+ T cells, which exhibited higher scores in the metastatic groups (Figure 2I).

KEGG enrichment scores were calculated for each cell, and the differences in scores between high and low hub protein enrichment score subgroups were analyzed. Among the top 20 enriched pathways, the high hub protein enrichment score subgroup was primarily enriched in KEGG-RIBOSOME, KEGG-SPLICEOSOME, and KEGG-ADHERENS-JUNCTION, while the low hub protein enrichment score group was mainly enriched in KEGG-COMPLEMENT-AND-COAGULATION-CASCADES (Figure 2J). At the cellular level, the top 20 pathways, displayed in a bubble plot, revealed that Proliferative cells1 were predominantly enriched in KEGG-PROTEASOME, KEGG-OXIDATIVE-PHOSPHORYLATION, and other pathways (Figure 2K; Supplementary Figure 3A). Expression of 60 T-cell co-stimulation-related genes (PMID: 37068318) and 15 HLA genes (PMID: 35533264) was visualized in heatmaps for Proliferative cells1 and Tumor_cells, showing lower HLA gene expression in some Tumor_cells (Figure 2L; Supplementary Figure 3B). Analysis of cell subtype distribution across metastatic subgroups and samples indicated that Neutrophils1 were mainly found in metastatic samples, particularly M1 samples, whereas Neutrophils2 were predominantly in non-metastatic samples, especially P1 samples (Figure 2M).

InferCNV analysis, initially performed on epidermal cells, endothelial cells, and fibroblasts, revealed substantial differential expression between Proliferative cells1 and Tumor_cells (Supplementary Figure 3C). Given the predominant origin of these cell types from M1 and P1 samples, respectively, CNV prediction was conducted separately for each sample. Expression heatmaps across chromosomal regions confirmed marked expression differences, suggesting a higher prevalence of CNV variants in these two cell types (Figure 3A). Subsequent CNV variant prediction and subtype identification facilitated the construction of clonal trees, each exhibiting seven branches representing distinct subtypes. Branches constituting less than 5% of the cell population in either sample were excluded from visualization. M1 samples exhibited predominant variants in chr11q_gain, chr5q_loss, chr6q_gain, chr1p_gain, and chr20q_gain, with branch B significantly enriched for chr6p_gain variants and branch M significantly enriched for chr22q_gain variants. CNV variant frequencies in P1 samples were significantly lower than those in M1 samples, indicating potential sample shifts with increased CNV variant frequencies (Supplementary Table 6). P1 samples primarily displayed variants in chr6p_gain, chr10q_gain, chr13q_loss, and chr15q_loss, with branch B significantly enriched for chr8q_gain variants (Figure 3B). Chi-square tests comparing branch frequencies between P1 and M1 samples revealed significantly higher frequencies of branches E and P in M1 samples, while branch G was significantly more frequent in P1 samples (Supplementary Table 7). UMAP plots illustrating cellular branch distributions demonstrated a significant difference in the distribution of branch P in M1 samples compared to RPS. Furthermore, the distribution of branch P in M1 samples exhibited greater overlap with high expression areas of RPS29 and IQGAP1, suggesting a potential association between these genes and metastatic grouping (Figure 3C).

To investigate the regulatory patterns and mechanisms by which tumor hub proteins drive metastasis, we focused on Proliferative cells1 and Tumor_cells, as prior analysis indicated their heightened malignancy and potential role in metastasis. Time-series analysis of these cell types’ expression matrices identified three distinct states. STATE2 exhibited the highest degree of stemness, suggesting its position at the beginning of the developmental trajectory, thus serving as our starting point (Figure 3D). Trajectory analysis of cells, grouped by metastatic potential, revealed that non-metastatic cells were primarily located in STATE2 and the other two states, while metastatic cells were predominantly found at the end of the latter two states, with a small presence in STATE2. This suggested that metastatic trajectories originated in STATE2 and progressed towards STATE1 and STATE3, mirroring our temporal analysis. Proliferative cells1 were mainly concentrated in STATE2, with a small fraction in late STATE1, while tumor cells were distributed across all states. The significantly higher stemness of Proliferative cells1 compared to Tumor_cells suggests that Proliferative cells1 may be associated with metastasis initiation (Figure 3E-G). Analysis of the four hub proteins’ expression across the trajectories, focusing on RPS29 and IQGAP1 due to low expression of SCGN and SH3GL1, showed that IQGAP1 expression gradually increased and stabilized with progression from low levels in early STATE2, with a slight increase at the end of STATE1 and a significant decrease at the end of STATE3. RPS29 was highly expressed throughout, with a slight decrease in the late stages and a significant decrease at the end of STATE3 (Figures 3H–I; Supplementary Figure 3D). Differential gene expression analysis across time, using the top 500 genes for subtype identification and heatmap visualization, yielded four gene expression subtypes. KEGG enrichment analysis revealed that cluster 1 genes, highly expressed in STATE1, were enriched in pathways like antigen processing and presentation and allograft rejection. Cluster 2 genes, highly expressed in STATE3, were enriched in protein processing in the endoplasmic reticulum and MAPK signaling pathways. Cluster 3 genes, also highly expressed in STATE3, were enriched in MAPK signaling pathways. Cluster 4 genes, highly expressed in STATE1 and STATE3, were enriched in ferroptosis and autophagy pathways. (Figure 3J; Supplementary Figure 3E; Supplementary Table 8).

Correlation analysis between isoform genes and hub proteins revealed a strong association of RPS29 with the cluster 4 gene group, while IQGAP1 exhibited a high correlation with all other gene groups (Figure 3K). Integrating cell subtype and temporal groupings, the low cell count in Proliferative cells1+State3 necessitated its merger with Proliferative cells1+State2, creating a combined Proliferative cells1+State23 group. Subsequent Hallmark functional enrichment analysis demonstrated elevated EPITHELIAL-MESENCHYMAL-TRANSITION pathway scores in state 2 and 3 cells, particularly Proliferative cells1+State1 and Tumor_cells+State1, suggesting an association between state 1 and metastasis initiation. Additionally, Proliferative cells1+State2 showed enrichment in NOTCH-SIGNALING, OXIDATIVE-PHOSPHORYLATION, and WNT-BETA-CATENIN-SIGNALING pathways (Figure 3L). Cell cycle analysis indicated a predominance of G1 phase cells across all three states, though state 2 exhibited a higher proportion of G2M/S phase cells, indicative of increased cell division activity. This observation was reinforced by cell subtype analysis, which revealed a majority of Proliferative cells1 in the G2M/S phase, consistent with rapidly dividing malignant tumor cells (Figures 3M, N). Finally, analysis of state distribution across cell subtypes and metastasis groupings showed a concentration of Proliferative cells1 in state 2, while metastatic original samples were predominantly represented in state 1, and non-metastatic samples in state 2, further supporting the association of state 1 cells with metastasis (Figure 3O).

To illustrate the distinct expression patterns across the four molecular subtypes, scatter plots were generated using the top 500 differentially expressed genes identified from the top 500 list related to timing (Figure 4A). Subsequently, utilizing signature genes characteristic of Proliferative cells1 and Tumor_cells at various time points, enrichment scores were calculated within the Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) cohort. Subsequent survival analyses revealed significant associations between Proliferative cells1 in State1 and poorer OS, DSS, and PFI; conversely, Proliferative cells1 in State2/3 correlated with worse OS, DFI, and PFI, but in this instance, lower enrichment scores were associated with worse outcomes. Finally, tumor cells in State1 demonstrated a significant association with worse OS, DSS, and PFI with higher enrichment scores indicative of poorer prognosis (Figures 4B–E).

Transcription factor activities were further calculated for all cells, and CSI values and RSS values were calculated from AUC values. Eight regulatory modules were identified based on CSI values (Figure 5A), and a UMAP plot, based on the mean CSI scores within each module, revealed significant differences in differentiation degree between modules, with distinct mean CSI values observed across the UMAP plot (Figure 5B). Intergroup-specific transcription factors were identified for metastatic and non-metastatic samples, and RSS ranking maps were generated. BACH2, EGR1, KLF2, DDIT3, and JUN were the top five RSS-ranked transcription factors in the metastatic group, while TCF7L2, CREB3L1, NR2F6, EHF, and TFAP2A were the top five in the non-metastatic group (Figure 5C). The top ten RSS-ranked transcription factors for each group were designated as group-specific, and their AUC activities were visualized via UMAP. Non-metastatic group-specific transcription factors exhibited higher activity in most B cells, T/NK cells, some fibroblasts, and myeloid cells, whereas metastatic group-specific transcription factors showed higher activity in most epithelial cells (Figure 5D). Functional enrichment analyses were performed on these group-specific transcription factors. In the non-metastatic group, enriched Gene Ontology Biological Processes (BP) included regulation of transcription from the RNA polymerase II promoter and DNA-templated transcription in response to stress; enriched Cellular Components (CC) included RNA polymerase II transcription regulator complex; and enriched Molecular Functions (MF) included DNA-binding transcription activator activity. Enriched KEGG pathways included viral carcinogenesis and the C-type lectin receptor signaling pathway. In the metastatic group, enriched BP included positive regulation of miRNA transcription; enriched CC included transcription regulator complex; and enriched MF included DNA-binding transcription activator activity. Enriched KEGG pathways included melanogenesis and the AMPK signaling pathway (Figures 5E, F).

Cell communication analysis, performed across all cells, revealed specific interactions detailed in Supplementary Table 9. These interactions were subsequently aggregated at the pathway level. In metastatic subgroups, intercellular interactions mapped 89 signaling pathways, while 92 pathways were identified in non-metastatic subgroups. Differences in pathway enrichment between these groups were assessed (Figure 6A; Supplementary Table 10), and the four pathways exhibiting the greatest intergroup variation were selected for network visualization and importance heatmap display. Metastatic subgroups showed significant enrichment in TNF, CD22, SPP1, and CLEC pathways. Notably, the TNF pathway, primarily involving Proliferative cells1+State1, Monocytes, Tregs, and Neutrophils1, was predominantly mediated by Proliferative cells2. Conversely, non-metastatic subgroups were significantly enriched in RESISTIN, VISFATIN, COLLAGEN, and NOTCH pathways. RESISTIN interactions primarily occurred between Monocytes and Macrophages, while the NOTCH pathway mainly involved epidermal cells, endothelial cells, and fibroblasts (Figure 6B). These findings highlight distinct cellular communication preferences between metastatic and non-metastatic subgroups at the pathway level.

Based on prior analysis, we identified Proliferative cells1+State1 and Tumor_cells+State1 as metastasis-associated cell types. We then combined the gene sets characterizing these two groups and calculated an enrichment score, stratifying samples into high- and low-enrichment groups based on the median score. Similarly, given the inference that RPS29 is a metastasis-associated hub protein, we also divided samples into high- and low-expression groups according to the median RPS29 expression. Predicting TCGA-STAD sample drug response IC₅₀ values using CTRP2 and PRISM database information and integrating the results revealed a significant negative correlation, specifically between TRIFLUOPERAZINE IC₅₀ values and both RPS29 expression and the aforementioned enrichment score, suggesting greater TRIFLUOPERAZINE sensitivity in patients within the high-expression/high-enrichment subgroups (Figures 6C, D). Subsequent database analysis identified 28 genes known to interact with TRIFLUOPERAZINE. Further correlation analysis of these interacting genes’ expression with TRIFLUOPERAZINE IC₅₀, enrichment score, and RPS29 expression revealed significant correlations for ABCD1, APOB, EPM2A, RGS4, and USP1 across all three metrics (Figures 6E, F).

Database analysis revealed that RPS29 may play a key role in the metastatic process of gastric cancer, while current reported studies are lacking. The higher expression of RPS29 protein in gastric cancer cell lines compared with normal human gastric mucosal epithelial cells was first identified by in vitro cellular experiments (Figure 7A). To explore the role of RPS29 in the malignancy of gastric cancer, transient transfection knockdown was performed in AGS and HGC27, gastric cancer cell lines with high expression of RPS29 (Figure 7B). Subsequent wound healing assays showed that knockdown of RPS29 significantly inhibited cell migration in AGS and HGC27 (Figure 7C). In further transwell migration and invasion assays, it was found that the number of migrating and invading cells in AGS and HGC27 was lower in the RPS29 knockdown group than in the control group (Figure 7D). Epithelial-mesenchymal transition (EMT) is a biological process that enables polarized epithelial cells to acquire mesenchymal characteristics, thereby promoting tumor migration and invasion. Therefore, it was investigated whether RPS29 expression affects EMT. Protein blotting results showed that the knockdown of RPS29 significantly decreased the expression of mesenchymal markers waveform protein and N-calmodulin but increased the expression of epithelial marker E-calmodulin in AGS and HGC27 (Figure 7E). The above experimental results indicated that RPS29 might promote migration, invasion, and EMT of gastric cancer cells (Figure 7F).

The technology roadmap (Figure 8) clearly illustrates the complete logical chain and key decision points. Data acquisition (proteomics, single-cell transcriptomics), Differential analysis and prognostic screening,PPI network construction and hub protein identification, Single-cell expression and functional association analysis, In vitro functional validation focusing on RPS29 (migration, invasion, EMT).