RNA sequencing data, along with corresponding clinical data, were acquired from the TCGA and GEO (GSE84433 and GSE62254) databases. The GSE62254 cohort (n = 300), also known as the ACRG cohort, was utilized as an independent external validation set to verify the generalizability of the prognostic model. For the TCGA-STAD cohort, FPKM values were transformed and normalized to TPM. We aggregated the two datasets, excluding patients lacking necessary clinical information. Consequently, 745 STAD patients were enrolled in our study. A set of 10 cuproptosis-associated genes was chosen for this study, aligning with previous publications (Supplementary Table S1).

Somatic mutation and copy-number datasets processed through the VarScan2 pipeline were obtained from TCGA. Mutation frequencies across the ten cuproptosis-related genes were summarized to highlight those most frequently altered, and CNV patterns were quantified in parallel. We then examined how copy-number status related to transcript abundance and mapped the chromosomal distribution of these genes.

ConsensusClusterPlus was applied to perform unsupervised subtype classification based on CRG expression patterns, allowing STAD cases to be partitioned into molecularly distinct groups. Subtype selection was guided by stability criteria, including a smoothly rising cumulative distribution function and avoidance of highly uneven sample allocation across clusters.

To determine the clinical significance of the molecular subgroups, we compared them against pathological features, CRG expression patterns, and patient outcomes. The variables examined included age, sex, tumor stage, and histological grade. Differences in overall survival among the subtypes were then evaluated using K-M analysis.

Immune and stromal characteristics for each STAD case were first inferred using the ESTIMATE scoring system. Next, we estimated relative infiltration levels across multiple immune cell populations via ssGSEA. Enrichment profiles were then used to evaluate multiple immune-related functional programs, including antigen-presentation mechanisms, checkpoint regulation, HLA activity, T-cell stimulatory or suppressive signaling, and IFN-γ–associated pathways.

Transcriptomic differences between the cuproptosis-related groups were first screened using limma in R, applying p < 0.05 and a 1.5-fold expression threshold. To explore pathway-level heterogeneity, GSVA was performed with KEGG-based gene sets (c2.cp.kegg.v7.4.symbols.gmt) sourced from MSigDB. The resulting gene list was then subjected to enrichment analysis using clusterProfiler and visualized with enrichplot to determine biologically relevant signaling patterns.

The cuproptosis-related prognostic risk score (CRG score) was formulated to numerically represent the cuproptosis pattern in each STAD specimen. Initial selection of DEGs linked to OS was performed via univariate Cox regression analysis. These prognostic DEGs were then used for unsupervised clustering, partitioning all patients into distinct cuproptosis gene clusters (A, B, C). All STAD cases (n = 745) from the integrated TCGA and GEO datasets were randomly split (1:1 ratio) into a training cohort (n = 373) and a validation cohort (n = 372). The training dataset served as the foundation for constructing the CRG score. Next, Lasso regression (R package “glmnet”) was employed, which was ultimately followed by multivariate Cox regression analysis to pinpoint the final set of candidates prognostic cuproptosis genes. Patients were subsequently divided into high- and low-risk strata using the median CRG score as the cutoff threshold. A correlation plot illustrated the relationship between OS status and the CRG score. To confirm the robustness and external applicability of the signature, the same risk scoring formula was applied to the independent validation cohort (GSE62254, n = 300). Patients in this external cohort were similarly stratified into high- and low-risk groups based on the median risk score calculated within that specific cohort.

Using the CRG-based risk measure, we assessed its linkage with broad immune cell composition via CIBERSORT and compared major stromal–immune features across risk categories with ESTIMATE. We next examined how individual prognostic CRGs aligned with immune infiltration patterns and evaluated the score in the context of immunotherapy-related indicators, including checkpoint expression, genomic instability, and tumor mutation characteristics. Finally, mutational differences between high- and low-risk tumors in the TCGA cohort were profiled using maftools.

To explore the potential of the CRG score in guiding treatment selection, we applied the pRRophetic framework to model drug sensitivity. Predicted IC50 values for commonly used chemotherapeutic and targeted agents were then contrasted between the two stratified risk categories.

Independent survival determinants were first screened using Cox regression, after which a nomogram combining the CRG score with key clinical features was generated to estimate long-term survival. Model performance was examined using calibration assessments and its potential clinical usefulness was evaluated through decision curve analysis (DCA).

Tissue samples consisted of 14 primary gastric carcinomas and their paired adjacent normal mucosa, obtained from patients undergoing resection at the First Affiliated Hospital of Zhengzhou University. Ethical clearance was granted by the institutional review committee (No. 2023-KY-0951-002). RNA was extracted using Trizol (Takara, Beijing, China), after which complementary DNA was produced with the PrimeScript RT kit equipped with a gDNA removal step. qRT-PCR amplification was then performed using SYBR Green reagents (Cowin Biosciences, Jiangsu, China) along with primers sourced from Sangon Biotech (Shanghai, China). GAPDH expression served as the internal reference for normalization.

The single-cell RNA-seq dataset (GSE163558) was acquired from the Gene Expression Omnibus (GEO) database. Raw data processing, quality control, and downstream analyses were performed using the Seurat package (version 5.0.1) within the R statistical environment (version 4.3.2). Strict quality control measures were implemented to filter out low-quality cells and potential doublets. Following data normalization and the identification of highly variable features, Principal Component Analysis (PCA) was applied for dimensionality reduction. Cluster visualization was subsequently achieved using the Uniform Manifold Approximation and Projection (UMAP) algorithm. Cell types were annotated based on the expression profiles of canonical marker genes. Finally, feature plots and violin plots were generated to visualize the specific expression patterns of risk genes across the identified cell clusters.

To investigate the potential protein-protein interaction between the core risk gene KRT17 and the cuproptosis regulator FDX1, molecular docking simulations were conducted using the HDOCK server (http://hdock.phys.hust.edu.cn/). This platform employs a hybrid algorithm combining template-based modeling and ab initio free docking. The protein structures were docked to calculate binding energies and docking scores. The model exhibiting the highest confidence score and the lowest binding energy was selected as the optimal binding conformation. The 3D visualization and structural analysis of the docking results were performed using PyMOL software (version 2.5.4; Schrödinger, Inc.).

All analyses were carried out in R (version 4.1.1). Group comparisons followed standard statistical procedures, with Student’s t-tests or Wilcoxon tests applied to two-group contrasts and ANOVA used when evaluating multiple groups. A p-value threshold of 0.05 was considered statistically significant.

Ten cuproptosis-related genes were examined in the STAD cohort. Somatic alterations occurred in roughly 13% of cases, with CDKN2A showing the highest mutation burden, followed by DLAT and several others at lower frequencies, while no mutations were detected in FDX1 or PDHA1 (Figure 2A). Representative mutation sites retrieved from cBioPortal included substitutions in PDHB and FDX1 (Figures 2B,C). Supplementary Figure S1 presents the common mutation site of other cuproptosis genes. Copy-number profiling revealed distinct alteration patterns, characterized by gains in GLS, MTF1, and LIPT1, whereas losses were observed in CDKN2A, DLAT, FDX1, PDHB, and DLD (Figure 2D). Most genes demonstrated significantly increased expression in tumor tissue compared with normal controls (Figure 2E). An analysis of copy number variations (CNVs) revealed that cuproptosis-associated genes exhibiting CNV gain, specifically GLS, MTF1, and LIPT1, showed a marked upregulation in their expression levels within STAD specimens. This finding suggests a potential mechanism where CNV influences the transcriptional activity of cuproptosis genes. Conversely, genes characterized by CNV deletion nonetheless presented with heightened messenger RNA (mRNA) abundance. Moreover, certain genes, such as LIAS, demonstrated no change in expression when contrasting tumor versus normal samples, despite frequent CNV gains or losses. This collective evidence implies that the control of cuproptosis gene expression is not solely dictated by CNV status.

Initial analysis employed univariate Cox regression to establish the prognostic relevance of the entire panel of ten cuproptosis genes in STAD patients (Supplementary Table S2). Drawing upon these regression findings, a prognostic network for cuproptosis was created, which mapped the overall pattern, reciprocal connections, and predictive capability of each cuproptosis gene within STAD (Figure 2F). Among the analyzed genes, five cuproptosis regulators-FDX1, LIAS, LIPT1, DLAT, and PDHA1-showed significant associations with overall survival. K-M analysis demonstrated that higher expression of these genes corresponded with improved survival outcomes in STAD patients (Figures 2G,H; Supplementary Figure S2). In addition, most cuproptosis genes displayed notable expression differences between tumor and adjacent normal tissue (Figure 2I).

Using consensus clustering of CRG expression profiles, the cohort was partitioned into two stable molecular groups (k = 2), designated Cluster A and Cluster B (Figure 3A; Supplementary Figure S3; Supplementary Table S3). This selected partition parameter ensured a comparative balance in the patient counts across the two resultant subgroups. Principal Component Analysis (PCA) demonstrated that the cuproptosis transcriptional landscape displayed a statistically significant distinction between Cluster A and Cluster B (Figure 3B). Despite this genetic separation, no significant differences in OS were noticed between the two groups (p = 0.487; Supplementary Figure S4). However, the subtypes demonstrated clear distinctions in CRG expression patterns and baseline clinical features (Figure 3C).

GSVA was used to characterize the biological divergence between the two cuproptosis-associated subgroups, and the resulting profiles indicated that Subtype B was predominantly linked to pathways governing cell-cycle progression and steroid biosynthesis, whereas Cluster A demonstrated a markedly different pattern defined by activation of apoptosis-related signalling and multiple immune pathways, including NK-cell cytotoxicity, adaptive receptor signalling, leukocyte migration, and chemokine networks (Figure 3D; Supplementary Table S4).

When ssGSEA was applied to quantify immune infiltration, a clear stratification emerged: most immune cell types showed substantially higher infiltration in Cluster A, matching the immune-enriched functional signature detected in GSVA (Figure 3E). This pattern was further reflected in immune activity scoring, where Cluster A consistently exhibited stronger enrichment across diverse immune functional categories, such as antigen presentation, immune checkpoint activity, cytolytic capacity, CCR signalling, HLA processes, inflammatory signalling, and type II IFN-γ responses (Figure 3F). To further contextualize these findings within the tumor microenvironment, ESTIMATE scoring was performed, and the results showed that Cluster A carried both elevated stromal and immune scores, suggesting a microenvironment characterized by higher stromal presence and heightened immune engagement compared with Cluster B (Figure 3G).

Differential expression analysis identified 103 cuproptosis-related differentially expressed genes (DEGs) distinguishing Clusters A and B, a process utilizing the R package “limma” (Supplementary Table S5). Functional annotation followed, beginning with GO analysis. This revealed significant enrichment of the DEGs in biological processes linked to the metabolism of intracellular components (Figure 3H).

Subsequently, a univariate Cox regression analysis was performed, yielding 26 prognostic DEGs that correlated with OS in STAD (p < 0.05, Supplementary Table S6). These prognostic DEGs then served as the basis for a further consensus clustering step aimed at uncovering underlying regulatory patterns. This refined analysis categorized STAD patients into three cuproptosis gene subtypes (k = 3, Figure 3I; Supplementary Figure S5; Supplementary Table S7). Expression levels of 6 key CRGs showed significant variations across these three resulting gene subtypes (Figure 3J). Across the three gene-defined groups, clinical characteristics were not evenly distributed, as illustrated in Figure 3K. When evaluating the tumor microenvironment, ESTIMATE scores indicated that Subtype C stood out with the highest stromal signal and a comparatively strong immune component (Figure 3L). Survival patterns further separated the subtypes, with Kaplan-Meier analysis showing that patients in Subtype A experienced the most favorable outcomes, while those in Subtype C had the shortest survival time (p < 0.001, Figure 3M).

The complete collection of STAD patient data from the TCGA and GEO cohorts was combined and subsequently partitioned into a training set and a validation set at an equal ratio (Figures 4A,B). This process concluded with the identification of five prognostic cuproptosis-related genes (CRGs) identified as contributing to high risk: ribosomal protein L39 like (RPL39L), paternally expressed gene 10 (PEG10), Synaptopodin 2 (SYNPO2), matrix metallopeptidase 11 (MMP11), and Keratin 17 (KRT17) (Supplementary Table S8). The CRG risk score was formulated using the coefficients derived from the multivariate Cox regression as: CRG risk score = (0.1098 * expression of RPL39L) + (0.1460 * expression of PEG10) + (0.1657* expression of SYNPO2) + (0.1701* expression of MMP11) + (0.0598* expression of KRT17). All STAD patients were then stratified into two cohorts, designated high-risk and low-risk, based on the median value of this calculated CRG risk score.

Analysis showed that the CRG risk score varied significantly across the two initial cuproptosis clusters and the three refined gene subtypes (Figures 4C,D). Gene subtype A exhibited the lowest mean CRG score. Notably, there was no meaningful statistical difference detected between gene subtypes B and C, a finding consistent with parallel K-M results. A Sankey plot was constructed to display how the two cuproptosis clusters, the three gene-defined subtypes, the calculated risk score, and patient survival status intersect across the cohort (Figure 4E).

Distribution analysis plots illustrated an inverse relationship, showing that OS survival time in STAD decreased corresponding to an elevation in the CRG risk score (Figures 4F,G). Similar correlation patterns were established when examining the entire patient cohort and the designated validation set (Supplementary Figure S6). Survival analysis demonstrated a clear separation between groups, with the low-risk population showing markedly longer OS compared with the high-risk group (p < 0.001, Figure 4H). ROC analysis further supported the model’s prognostic utility, with AUC values of 0.645, 0.698, and 0.704 for predicting 1-, 3-, and 5-year survival, respectively (Figure 4I), and similar performance was observed in the full and external validation sets (Supplementary Figure S7).

To further substantiate the robustness and clinical generalizability of the 5-gene signature, we employed the GSE62254 dataset ((also known as the ACRG cohort, n = 300) as a large-scale independent external validation cohort. We calculated the risk score for each patient in the ACRG cohort using the same formula constructed in the training set and stratified patients into high- and low-risk groups based on the median risk score. Consistent with our internal findings, Kaplan-Meier analysis revealed that patients in the high-risk group of the ACRG cohort exhibited significantly worse overall survival compared to those in the low-risk group (p < 0.0001, Supplementary Figure S8). These results strongly confirm that our cuproptosis-related prognostic model maintains high stability and predictive power across different populations and sequencing platforms.

Several cuproptosis-related genes showed distinct expression differences between high- and low-risk groups (Figure 4J). A heatmap of the training set demonstrated that all five prognostic markers were expressed at higher levels in the high-risk category (Figure 4K). Clinical comparison based on risk stratification showed that adverse features such as grade 3 tumors and T3/T4 staging occurred more frequently in the high-risk cohort (Table 1).

Protein expression data from HPA repository indicated differential protein abundance in STAD versus normal tissue for some of the signature genes. KRT17 protein levels were higher in STAD, implying its possible function as a risk factor. In contrast, RPL39L and SYNPO2 displayed reduced expression in tumor samples, suggesting they may serve protective roles. Notably, no significant differences in the protein expression of MMP11 and PEG10 were observed (Figures 5A–E).

To verify these observations at the transcriptional level, RT-qPCR was conducted using fourteen matched pairs of STAD and adjacent normal tissues. The transcript-level validation (Figures 5F–J) showed that KRT17, MMP11, and RPL39L were substantially overexpressed in the tumor samples, while SYNPO2 expression was diminished. Transcript levels of PEG10 showed no significant difference. Comprehensive RT-qPCR data and associated primer sequences are provided in Supplementary Tables S9, S10.

To further elucidate the cellular source and potential molecular mechanism of the risk signature, we performed single-cell RNA-seq and molecular docking analyses. Among the five genes included in the prognostic model, KRT17 was identified as the most critical risk component with the strongest correlation with poor survival, thus serving as the representative target for mechanistic exploration.

First, we explored the potential interaction between KRT17 and the cuproptosis machinery. Molecular docking simulation visualized a stable complex formed between the KRT17 protein and the core cuproptosis regulator FDX1, with a remarkably low binding energy of −204.18 kcal/mol (Figure 6A). To verify the robustness of this interaction, we screened the top 100 docking models, which consistently demonstrated high thermodynamic stability (Supplementary Table S11). This strong thermodynamic affinity suggests that KRT17 may physically sequestrate or modulate FDX1, thereby interfering with the copper-induced cell death process.

Next, to validate the cellular origin of KRT17, we analyzed the single-cell transcriptomic dataset GSE163558. Prior to downstream analysis, rigorous quality control was performed to filter out low-quality cells and doublets, ensuring the reliability of the dataset (Supplementary Figure S9). Using the UMAP algorithm, we successfully deconvoluted the tumor microenvironment into nine distinct cell clusters, including B cells, T cells, macrophages, endothelial cells, fibroblasts, and malignant epithelial cells (Figure 6B). Identity validation using canonical markers confirmed the accuracy of cell annotation (Figure 6F).

Visualizing the expression distribution of KRT17 revealed a striking pattern: unlike immune markers, KRT17 was predominantly enriched in the malignant epithelial cell cluster, with negligible expression in immune lineages such as T cells or macrophages (Figures 6C,D). Furthermore, pathway activity analysis (GSVA) highlighted significant heterogeneity in signaling pathway activation across these cell types, suggesting that KRT17-high tumor cells may harbor distinct oncogenic features (Figure 6E).

Collectively, these results provide robust evidence that the high-risk signature is driven by tumor-intrinsic alterations (specifically KRT17 overexpression) rather than primary immune cell dysfunction, which may subsequently influence the immune microenvironment.

Using CIBERSORT to quantify immune infiltration, a clear relationship emerged between the CRG-based risk signature and multiple immune cell populations (Figure 7A; Supplementary Table S12). Using the ESTIMATE algorithm, we observed that stromal scores increased in parallel with higher CRG risk classifications in STAD, while immune scores showed no meaningful separation between the high- and low-risk groups (Figure 7B; Supplementary Table S13).

To extend this analysis to the immune landscape, expression of the five prognostic cuproptosis genes was systematically correlated with inferred abundances of 22 distinct immune cell populations (Figure 7C). Critically, the SYNPO2 gene correlated positively with resting memory CD4⁺ T cells and M2 macrophages, yet negatively correlated with M1 macrophages. This observation suggests that SYNPO2 may be closely involved in the interaction between cuproptosis and the TIME. Further assessment of the CRG prognostic risk score’s relationship with immune checkpoint molecules was performed (Figure 7D). This analysis identified nine immune checkpoints exhibiting differential expression between the distinct cuproptosis risk groups: tumor necrosis factor superfamily member 4, inducible T-cell CO-Stimulator, CD44, Neuropilin-1, CD276, KIR3DL1, LGALS9, and CD160.

To explore how the CRG-based risk model relates to immunotherapy-relevant features, we compared somatic alterations and biomarker patterns between the two risk groups. The low-risk subgroup demonstrated a higher mutation rate (90.57%) than the high-risk subgroup (83.67%) (Figures 8A,B), with TTN and TP53 being the most common mutations. TMB was substantially reduced in the high-risk group and showed a negative association with the cuproptosis signature (Figure 8C; Supplementary Figure S10). The CRG score also correlated inversely with the CSC index, indicating that higher-risk tumors displayed fewer stem-like features and a more differentiated phenotype (Figure 8D). MSI patterns further supported this trend: MSI-H occurred more frequently in the low-risk group, and MSI-H cases carried lower risk scores than MSI-L or MSS tumors (Figures 8E,F). Analysis of antigen-presentation–related genes showed reduced expression of HLA-DMA and HLA-F in high-risk tumors compared with low-risk samples (Figure 8G).

Predicted drug response profiling revealed treatment-related divergence between subgroups. High-risk patients appeared more responsive to docetaxel, lapatinib, pazopanib, and imatinib, as reflected by lower estimated IC50 values, whereas sorafenib sensitivity was greater in the low-risk group (Figures 8H–L). Collectively, these data suggest that the CRG-based score may provide value for predicting immune features and informing therapy selection in STAD.

Survival determinants were explored by first fitting univariate and then multivariate Cox proportional hazards models, through which age, TNM stage, and the cuproptosis-based risk index emerged as variables with independent prognostic weight in STAD (Figures 8M,N). These factors were subsequently incorporated into a graphical prediction tool, allowing clinicians to read off overall survival probabilities at early, intermediate, and later follow-up intervals corresponding to approximately 1, 3, and 5 years (Figure 8O). The resulting nomogram yielded a C-index of 0.746, consistent with robust discrimination. Agreement between estimated and actual outcomes was further supported by well-aligned calibration curves across all evaluated time windows (Figure 8P). In decision curve analysis, this integrated model provided a higher net clinical benefit than use of the TNM staging system alone (Figure 8Q).