Abstract Background Gastric cancer (GC) is a malignancy known for its high fatality rate. Disulfidptosis, a potentially innovative therapeutic strategy for cancer treatment, has been proposed. Nevertheless, the specific involvement of disulfidptosis in the context of GC remains uncertain. Methods The mRNA expression profiles were obtained from the TCGA and GEO databases. Univariate and LASSO Cox regression analyses were employed to identify differentially expressed genes and develop a risk model for disulfidptosis-related genes. The performance of the model was evaluated using Kaplan-Meier curve, ROC curve, and nomogram. Univariate and multivariate Cox regression analyses were conducted to determine if the risk model could serve as an independent prognostic factor. The biological function of the identified genes was assessed through GO, KEGG, and GSEA analyses. The prediction of drug response was conducted employing the package “pRRophetic”. Furthermore, gene expression was determined using qRT-PCR. Results An eight-gene signature were identified and utilized to categorize patients into low- and high-risk groups. Survival, receiver operating characteristic (ROC) curve, and Cox analyses provided clarification that these eight hub genes served as a favorable independent prognostic factor for patients with GC. A nomogram was constructed by integrating clinical parameters with the risk signatures, demonstrating high precision in predicting 1-, 3-, and 5-year survival rates. Additionally, drug sensitivity was different in the high-risk and low-risk groups, and the expression of three genes was verified by qRT-PCR. Conclusion The prognostic risk model developed in this study demonstrates the potential to accurately forecast the prognosis of patients with GC. Keywords: Gastric cancer, Disulfidptosis, Prognostic signature, Immunotherapy response, Drug sensitivity 1. Introduction Gastric Cancer (GC) ranks as the fifth most prevalent malignant tumor globally and the fourth leading cause of cancer-related mortality [[25]1]. Epidemiological data reveal higher age-standardized cancer mortality for gastric and liver cancer in China compared to the United States and the United Kingdom for both genders [[26]2]. Notably, the occurrence of GC among individuals under 50 is progressively increasing, this rise varying with the geographical region [[27]3]. A significant proportion of patients with GC receive a diagnosis at an advanced stage due to insufficient early detection measures [[28]4]. Several risk factors associated with GC have been identified, including Helicobacter pylori infection, age, sex, smoking, obesity, metabolic dysfunction, dietary factors, alcohol use, medication, and host genetics [[29]5]. Immunotherapy has emerged as an effective therapeutic approach for combating tumors. Specifically, immune checkpoint inhibitors (ICI) enhance T cell anti-tumor effects by inhibiting their corresponding negative regulatory mechanisms [[30]6]. Despite the transformative impact of ICI on the treatment of advanced GC, its efficacy is limited in most patients [[31]7]. As stated earlier, immunotherapy is currently not good for GC, so our research needs to improve the efficiency of immunotherapy for GC. Cell death is highly intricate, encompassing various mechanisms by which cells perish. These include abrupt physiological impairments due to intense and rapid external influences, commonly referred to as accidental cell death (ACD). Conversely, cells can undergo regulated cell death (RCD) through specific molecular pathways that can be manipulated pharmacologically or genetically. Under physiological circumstances, RCD encompasses various types of programmed cell death mechanisms, including programmed cell death [[32]8]. Recent investigations unveiled a novel form of cell death termed disulfidptosis. Disulfidptosis primarily occurs in cancer cells with high expression of SLC7A11. This is mainly due to insufficient NADPH supply, which is needed to reduce cystine to cysteine, causing disulfide stress. Furthermore, introducing additional cystine to a glucose-free medium also leads to the excessive consumption of NADPH in low SLC7A11 expression cells. This, in turn, induces actin protein crosslinking and cytoskeletal contraction, ultimately leading to disulfidptosis [[33]9]. The mechanism of disulfidptosis in gastric cancer is still in its nascent stage in the existing literature, necessitating further investigations to elucidate whether activating disulfidptosis could serve as a viable therapeutic approach for the treatment of gastric cancer. However, the precise mechanism of disulfidptosis in GC remains unknown. In this study, we employed lasso and Cox regression techniques to identify prognostic disulfidptosis signature genes. The risk scores derived from differentially expressed genes (DEGs) demonstrated high accuracy in predicting various clinical outcomes, including prognosis, immune invasion, microsatellite instability (MSI), tumor mutational burden (TMB), stem cell analysis, immunotherapy response, and drug sensitivity. Additionally, we performed molecular docking and molecular dynamics simulations on shikonin, predicted to possess anticancer properties. Our study effectively illustrated the significance of disulfidptosis patterns in GC for prognostic prediction, immune infiltration, immunotherapy response, and drug sensitivity. In summary, our findings offer a novel discovery, as they demonstrate for the first time that assessing the risk associated with disulfidptosis can unveil its significance in the context of GC. This discovery provides a fresh perspective for clinical immunotherapy and medication strategies in patients with GC. 2. Materials and methods 2.1. Sample collection and processing Data of patients with GC were collected from The Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) and Gene Expression Omnibus (GEO) [34]GSE84437 and [35]GSE84433 databases. Complete RNA-seq data were downloaded from both databases, along with mutation information for TCGA samples. The “SVA” R software package merged standardized TCGA-STAD (407 samples) and [36]GSE84437 (433 samples) patient data, retaining only one copy of each transcript. The data is batch corrected and the corrected results are outputted using code commands. [37]GSE84433 from the GEO database as an external data validation model. This approach aimed to maximize the availability of transcription data for future studies while maintaining analytical accuracy. 2.2. Consistent gene cluster analysis for disulfidptosis In Liu's study [[38]9], ten disulfidptosis-associated genes were identified. The ConsensusClusterPlus R software package was utilized to cluster and evaluate the data, considering the expression profiles of these genes for consistency analysis. The “survival” R software package was employed to examine survival disparities among the identified clusters. Additionally, the R software packages “igraph”, “psych”, “reshape2″, and “RColorBrewer” were utilized to construct a network representation of disulfidptosis-related genes. A heatmap was generated to visualize clinical characteristic dissimilarities between the clusters, using the “pheatmap” function. 2.3. Evaluation of tumor microenvironment (TME) of molecular subtypes To investigate the relationship between identified subtypes through clustering and the TME, we utilized an R software package to assess the scores of all samples. The R “ESTIMATE” package was employed to export TME and immune scores for all patients with GC. To examine the differential immune signatures within the clusters, we employed a single-sample gene set enrichment analysis (ssGSEA). The “GSEABase” and “GSVA” packages were utilized to assess immune characteristics. Cluster analysis was conducted using the “limma” R software package, and the resulting outcomes were visualized using the “ggpubr” box plot. 2.4. Functional enrichment analysis between molecular subtypes To investigate potential biological functions between clusters, gene set variation analysis was conducted. Function expansion was performed using MSigDB data (“c2. Cp. Kegg. V7.4. Symbols. GMT” and “c5. Go. V2023.1. Hs. Symbols. GMT”). The genomic enrichment pathway was determined using the “GSVA” package. Statistical significance between clusters was defined as an adjusted P-value <0.05. 2.5. Identification of DEGs in disulfidptosis subtype The DEGs within the disulfidptosis subtypes were identified using the “limma” software package. In this study, |log2FC| > 0.585 was selected, that is, the genes up-regulated or down-regulated by 1.5 times [[39]10,[40]11]. DEGs between clusters were determined based on an adjusted |log2FC| > 0.585 and an adjusted P-value <0.05. “ClusterProfiler” was used to analyze the functional enrichment of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for DEGs in each disulfidptosis subtype. 2.6. Establishment and verification of disulfidptosis-related prognostic model To quantitatively assess the correlation between disulfidptosis pattern and GC, we implemented a disulfidptosis-related prognosis model. Survival-related genes within DEGs were identified through clustering DEGs and subsequent univariate Cox analysis. Additionally, a lasso-Cox analysis was used to develop a prognostic risk model for disulfidptosis, employing the “caret” and “glmnet” R software packages. GC samples were evenly divided into training and test sets. The disulfidptosis risk score was calculated by applying the model equations based on the expression profiles of pivotal genes. The model formula utilized in this study was the following: Risk score = ∑n i = 1 expi * coefficient (where expi represents gene expression, and coef represents gene risk factor). Risk score classification followed common rules, specifically, those of the median method, categorizing patients with GC into high- or low-risk groups. The clinical prognosis significance of the risk score was assessed through Kaplan-Meier survival analysis. Model predictive accuracy was validated using a receiver operating characteristic curve (ROC). 2.7. TME and immune status in different risk groups To estimate tumor purity and infiltration of immune cells across different cancer types, we computed the immune score for each sample within the TCGA-STAD cohort, using the “estimate” package in R. The CIBERSORT algorithm in R software was utilized to evaluate 22 immune cell types within the TME of each sample. 2.8. Analysis of TMB and drug sensitivity The “maftools” R software package was employed to perform mutation spectrum mapping for the two risk groups, facilitating the visualization of the frequency and types of mutated genes [[41]12]. Box diagrams were utilized to visually depict the disparities between TMB risk groups. The R software package “pRRophetic” was employed to compute the maximum half inhibitory concentration (IC50) values of chemotherapy agents and targeted agents for high- and low-risk STAD groups, forecasting therapeutic efficacy. The cut-off criteria for drug sensitivity were set at P < 0.001 [[42]13]. 2.9. Molecular docking and molecular dynamics simulation Based on network pharmacology results, the 3D structures of the target proteins were obtained from the Protein Data Bank (PDB) ([43]http://www.rcsb.org/) [[44]14]. They were conFig.d in AutoDock4 as follows: water was removed, hydrogen was replaced, and a receptor designation was made. The drug was similarly prepared and set as a ligand. The software automatically conFig.d the torsion tree. The structure was saved as a PDBQT protein-ligand file. Receptor and ligand PDBQT structures were imported into AutoDock4, and docking parameters were defined by setting their docking ranges. Molecular docking in AutoDock4 involved inserting small drug molecules into macromolecular protein constructs, along with configuring operating methods and docking parameters. Subsequently, docking is performed. Molecular dynamics simulations were performed as previously described [[45]15]. 2.10. Correlation analysis of satellite instability and stem cells The prognosis of GC can be influenced by the status of MSI [[46]16]. The correlation between risk score and MSI status was examined using the R software packages “plyr”, “ggpubr”, and “ggplot2” [[47]17]. Stem cell correlation analysis was performed using R software packages “limma”, “ggpubr”, “ggplot2″, and “ggExtra”, employing Spearman's test [[48][18], [49][19], [50][20]]. 2.11. Construction and verification of nomogram A nomogram was developed to assess the risk of disulfidptosis and predict the prognosis of GC. The nomogram was constructed using various clinical features and risk scores of patients with GC, employing the “rms” R software package. Subsequently, the survival probabilities at 1, 3, and 5 years were determined based on individual patient scores. Calibration and ROC curves were employed to evaluate the predictive accuracy of the nomogram. 2.12. Cell culture and validation of gene expression The GES-1 immortalized gastric cell line and two human gastric cancer cell lines, MGC-803 and SGC‐7901, were procured from the Chinese Academy of Science in Shanghai, China. These cell lines were cultured in Roswell Park Memorial Institute (RPMI) 1640 medium (Hyclone, Logan, Utah, USA) and maintained in a humidified chamber at 37 °C with 5% CO[2], adhering to the guidelines provided by the manufacturer. Cell samples were subjected to RNA extraction using TRIzol (Invitrogen) followed by reverse transcription into cDNA utilizing the Reverse Transcription Kit (Vazyme, Nanjing, China) as per the manufacturer's guidelines. Subsequently, qRT-PCR was performed using the HiScript® Q RT SuperMix for the qPCR kit (Vazyme) along with the appropriate primers. The relative expression levels were determined using the 2^−ΔΔCt method. The primer sequences were as follows: OXSM: Forward: 5′-TGGTGTTGGAACTCACCTGG -3′, Reverse 5′- CCGGGCTCTTCATCACTACC -3’; SLC3A2: Forward: 5′- AGCTGGAGTTTGTCTCAGGC-3′, Reverse 5′-GGCCAATCTCATCCCCGTAG -3’; SLC7A11: Forward: 5′- AGATGATGATACCTGCCTGTCTG-3′, Reverse 5′-CTGAATTGAGCAATACAAGGAAGC-3’; GAPDH: Forward: 5′- GACAGTCAGCCGCATCTTCT-3′, Reverse 5′-GCGCCCAATACGACCAAATC-3’. 2.13. Statistical analysis All statistical analysis and visualization were conducted using R software (version 4.3.0). A P-value <0.05 was considered statistically significant. 3. Results 3.1. Genetic study of genes related to disulfidptosis in GC To investigate the involvement of disulfidptosis in GC, this study incorporated 10 disulfidptosis-related genes. The expression levels of these genes exhibited substantial differences between the normal and GC cohorts. Notably, NDUFS1 exhibited no statistically significant difference in expression between unpaired and paired samples. Similarly, the levels of expression of NUBPL, NDUFA11, NCKAP1, and GYS1 showed no significant differences in paired samples. Conversely, the expression of other genes exhibited a significant increase in GC ([51]Fig. 1A and B). Correlation analysis of disulfidptosis-related genes revealed positive correlations among all genes except NDUFS1, OXSM, LRPPRC, NUBPL, NCKAP1, RPN1, and NUDFA11, which exhibited a negative correlation ([52]Fig. 1C). Notably, the highest positive correlation (0.72) was observed between NCKAP1 and NDUFS1. Additionally, the regulatory analysis network provided evidence of gene interactions and highlighted the potential prognostic significance of this relationship for GC ([53]Fig. 1D). Furthermore, the genetic map ([54]Fig. 1E) displayed the chromosomal locations of the copy number variation (CNV) genes associated with disulfidptosis, and this CNV occurrence was observed across all 10 genes. Notably, NUBPL and SLC3A2 exhibited high amplification, while NCLAP1, LRPPRC, OXSM, and NDUFA11 exhibited a higher frequency of deletions ([55]Fig. 1F). Subsequent analysis revealed that 13.16% of the 433 samples exhibited 57 mutations in total. SLC3A2 displayed the highest number of mutations, mainly frameshift deletions, followed by GYS1 and NCKAP1, which exhibited missense mutations. NDUFA11 exhibited no mutations ([56]Fig. 1G). Fig. 1. [57]Fig. 1 [58]Open in a new tab Expression of disulfidptosis-regulatory genes in GC. A Comparison of expression between normal and tumor tissues of unpaired samples. B Comparison of expression between normal and tumor tissues of paired samples. C Spearman correlations between genes. The color blue represents negative regulation; the color red represents positive regulation. D Comprehensive network map combining disulfidptosis regulator interactions and prognosis. E Distribution of disulfidptosis-related genes on chromosomes. F CNV values for disulfidptosis regulators in GC specimens. G Mutation frequencies of the 10 regulatory genes in 433 GC specimens. *P < 0.05; **P < 0.01; ***P < 0.001. 3.2. GC classification pattern based on disulfidptosis To investigate disulfidptosis gene expression characteristics in GC, we analyzed 407 study samples from the TCCA-STAD dataset using clustering analysis. Employing a commonly used algorithm, we based this analysis on the expression data of the 10 disulfidptosis genes. The results of the consensus clustering, visualized through a cumulative distribution function (CDF) diagram, indicated an optimal cluster number of k = 2. Notably, each sample within the cluster exhibited a strong correlation ([59]Fig. 2A), leading to the categorization of patients with GC into two distinct groups: geneCluster A and B. Principal component analysis highlighted distinct directions within the cluster ([60]Fig. 2B). Furthermore, a notable disparity in patient outcomes was observed, with geneCluster B patients exhibiting superior overall survival ([61]Fig. 2C). The heat map results revealed higher expression levels of the disulfidptosis-related genes in geneCluster B compared to geneCluster A ([62]Fig. 2D). Additionally, GSEA enrichment analysis unveiled significant differences in GO between the clusters, particularly in metabolic, tumor signaling pathways, and immune-related pathways ([63]Fig. 2E). Cluster A exhibited a higher degree of enrichment in metabolic processes such as monoatomic ion transmembrane transport, chemical synaptic transmission, and regulation of postsynaptic membrane potential. On the other hand, cluster B demonstrated a greater enrichment in cysteine-type peptidase activity, positive regulation of telomere maintenance, and ER-nucleus signaling pathways. Pathway enrichment analysis indicated that cluster A displayed increased activity in neuroactive ligand-receptor interaction and calcium signaling pathway, while cluster B was primarily associated with amino sugar and nucleotide sugar metabolism, protein export, basal transcription factors, glyoxylate- and dicarboxylate-metabolism, and the pentose phosphate pathway ([64]Fig. 2F). Subsequently, an examination of the infiltrating immune cells within the cluster ([65]Fig. 2G) revealed significant variations in B cells, different types of T cells, mast cells, monocytes, and neutrophils. Fig. 2. [66]Fig. 2 [67]Open in a new tab Identification of disulfidptosis clusters. A Consensus matrix based on disulfidptosis regulator expression in the GC cohort at k = 2. B PCA of consensus matrix when k = 2. C K-M survival analysis between the two disulfidptosis clusters. D Heatmap of clinical characteristics between the two subtypes. E Biological pathways analysis between the two disulfidptosis clusters. F Multiple pathways analysis between the two disulfidptosis clusters. G Infiltration of immune cells between the two disulfidptosis clusters. *P < 0.05; **P < 0.01; ***P < 0.001. 3.3. Characterization of phenotypes associated with disulfidptosis Differential genes were identified and classified based on the type of disulfidptosis, with optimal selection indicated at k = 2 ([68]Fig. 3A). Expression analysis of disulfidptosis genes in different subtypes revealed significant differences in the expression levels of GYS1, NDUFS1, OXSM, LRPPRC, NUBPL, RPN1, SLC3A2, and SLC7A11 ([69]Fig. 3B). Furthermore, gene cluster B exhibited significantly increased survival time compared to gene cluster A ([70]Fig. 3C). To identify prognostic-related genes, we utilized R-packages “limma” and “survival” to analyze patient clinical data, resulting in the generation of 194 prognostic-related DEGs ([71]Fig. S1). GSEA analysis indicated that DEG enrichment was associated with DNA repair, regulation of the cell cycle, and immune cell regulation. This finding confirms the close relationship between disulfidptosis and processes such as the cell cycle, DNA replication, and the regulation of the TME ([72]Fig. 3D and E). These findings suggest the existence of two distinct regulatory patterns associated with disulfidptosis in GC, with functional enrichment differences between subtypes possibly influencing observed patient outcome variations. Fig. 3. [73]Fig. 3 [74]Open in a new tab Gene cluster determination. A Consensus matrix based on disulfidptosis-related gene expression in the GC cohort at k = 2. B Differential expression of disulfidptosis regulators in the two clusters. C K-M survival analysis of disulfidptosis regulators in the two clusters. D Multiple GSVA analysis between the two clusters of biological pathways. E Multiple GSVA analysis between the two clusters of pathways. **P < 0.01; ***P < 0.001. 3.4. Construction and verification of disulfidptosis prognosis model Classifying disulfidptosis holds significant potential for clinical prognosis in GC. To further investigate the characteristics of disulfidptosis, we developed a prognostic model. Lasso-Cox regression analysis identified 8 key genes from 194 prognostic genes ([75]Fig. 4A and B), forming the basis of a disulfidptosis risk score. These 194 genes were p < 0.05 genes in cross-validation, namely lasso-cox regression significant genes, which were obtained by disulfidptosis typing, equivalent to new disulfidptosis genes. The risk score formula was established as: (−0.1866 × SLC27A2) + (0.1240 × PPP1R14A) + (0.1623 × PHLDA1) + (0.1156 × LAMC2) + (0.1157 × KRT80) + (0.1571 × NMU) + (0.1111 × GPC3) + (0.05 53 × SFRP2). Subsequently, GC samples were categorized into high-risk and low-risk groups based on these risk scores, validated by Kaplan-Meier survival curves ([76]Fig. 4C–E). The analysis of the risk curve indicated that an increase in the risk score was associated with an elevated risk of death ([77]Fig. 4F–H). The disulfidptosis risk score demonstrated notable accuracy in predicting 1-, 3-, and 5-year outcomes, with values of the area under the curve (AUC) for the overall samples of 0.621, 0.644, and 0.668, respectively ([78]Fig. 4I). Internal training set AUC values were 0.653 (1 year), 0.652 (3 years), and 0.703 (5 years) ([79]Fig. 4J). Similarly, test set AUC values were 0.586 (1 year), 0.630 (3 years), and 0.625 (5 years) ([80]Fig. 4K). The AUC values obtained from the external dataset [81]GSE84433 were 0.704 (1 year), 0.663 (3 years), and 0.688 (5 years) ([82]Fig. S2). These results demonstrate a higher level of accuracy in model construction. Enhancing prognostic prediction, we constructed a histogram of the disulfidptosis risk score. The characteristics of the patients with GC included in this study can be found in [83]Fig. S3. Fig. 4. [84]Fig. 4 [85]Open in a new tab Construction of the risk signature. A Cross-validation for selection of tuning parameters in LASSO regression. B LASSO coefficient profiles of candidate DEGs. C-E Survival curves for risk score groups in total sample, train, and test data. F–H Survival status curves for risk scores in total sample, train, and test data. I–K ROC curves for total sample, train, and test data. L Columnar plots predicting one, three, and five years overall survival in patients with GC. J Calibration curves for nomogram-predicted survival outcomes. ***P < 0.001. Notably, survival data from the GC cohort in this study did not include patients with a survival time of 0. We considered age, sex, T stage, and N stage as clinical characteristics, calculating 1-, 3-, and 5-year survival rates based on the corresponding scores ([86]Fig. 4L). The obtained survival rates were 0.73, 0.401, and 0.287, respectively. A calibration curve demonstrated alignment between the predicted probabilities from the column graph and the actual probabilities ([87]Fig. 4M). In conclusion, our findings highlight the promising potential of the disulfidptosis risk score for prognostic assessment. 3.5. Independent predictive value of disulfidptosis risk prognostic model We developed a prognostic model for disulfidptosis risk, revealing diverse survival outcomes. Expanding our investigation, we integrated comprehensive clinical data across all samples, identifying four common clinical features: age, sex, grade, and stage. Our analysis confirms risk scores as potent prognostic indicators for GC ([88]Fig. 5A and B). Additionally, genecluster B exhibited significantly lower risk than genecluster A ([89]Fig. 5C). Furthermore, significant differences in gene expression used in constructing the disulfidptosis model were observed between the two risk groups ([90]Fig. 5D). Sankey's map demonstrated that the low-risk cohort primarily consisted of patients with cluster B, including a significant proportion of survivors. Conversely, higher risk scores and association with cluster A were prominent among patients who succumbed to the disease ([91]Fig. 5E), aligning with previous analyses. Fig. 5. [92]Fig. 5 [93]Open in a new tab Survival predictive value of disulfidptosis risk model. A Forest plot of risk scores and clinical characteristics for univariate Cox regression analysis. B Forest plot of risk scores and clinical characteristics for multivariate Cox regression analysis. C Differences in risk scores between the two geneclusters. D Differential expression of disulfidptosis regulators in the two risk groups. E Alluvial diagram showing disulfidptosis clusters, gene clusters, risk grouping, and survival status. ***P < 0.001. 3.6. The disulfidptosis risk score has considerable potential in predicting tumor treatment outcomes Immune infiltration analysis was performed across GC risk groups ([94]Fig. 6A), revealing significant differences in factors such as cytolytic activity, human leukocyte antigen type, mast cells, natural killer cells, inflammation-promoting factors, major histocompatibility complex class I, and T cell co-inhibition scores between the risk groups. Additionally, T cell co-stimulation, Th1 cells, Th2 cells, and type II interferon response also varied notably. The correlation between immune cells and risk genes is depicted in [95]Fig. 6B, demonstrating their significant association. Furthermore, the tumor immune dysfunction and exclusion (TIDE) score serves as an evaluative tool for the clinical efficacy of immune checkpoint inhibitor (ICI) therapy, being employed to appraise the response of both high-risk and low-risk populations to immunotherapy. A greater TIDE predictive score corresponds to an increased probability of immune evasion, indicating that patients are less inclined to derive benefits from immunotherapy [[96]21]. Our analysis revealed a substantial disparity in TIDE scores between the low-risk and high-risk cohorts ([97]Fig. 6C). Risk scores correlated significantly with TIDE scores. Furthermore, the analysis of tumor stemness revealed a significant inverse correlation between risk score and replicative senescence-associated secretory phenotype ([98]Fig. 6D). Fig. 6. [99]Fig. 6 [100]Open in a new tab Immune landscape analysis. A Immune cell scores between high-risk and low-risk groups. B Correlation between immune immune cells and disulfidptosis regulators. C TIDE scores in different risk groups. D Relationships between risk scores and RNAss values. E The proportion of MSI-L + MSI-H in the low-risk group was higher than that in the high-risk group. F Comparison of the risk score among MSS, MSI-L samples, and MSI-H samples. *P < 0.05; **P < 0.01; ***P < 0.001. Based on the microsatellite instability analysis, patients were categorized into three groups: microsatellite stabilization group (MSS), low-frequency microsatellite instability group (MSI-L), and high-frequency microsatellite instability group (MSI-H). In the low-risk group, a higher proportion of patients belonged to the MSI-H subgroup, while the high-risk group had more patients belonging to the MSS and MSI-L subgroups ([101]Fig. 6E). There was no statistically significant difference in risk values between MSS and MSI-L groups. Risk values differed significantly between MSS and MSI-H subgroups, with the MSI-H subgroup having a lower risk ([102]Fig. 6F). TMB serves as an indicator of mutation counts in cancer [[103]22]. A higher TMB increases self-neoantigens and immunogenic recognition, enhancing the likelihood of T cell recognition, facilitating T cell responses against tumor cells, and indicating improved outcomes with ICI [[104]23]. “Maftools” software package was used to analyze somatic mutation distribution in TCGA and [105]GSE84437 cohorts. [106]Fig. 7A and B revealed a higher TMB in the low-risk group (mutation rate 92.77% with 154 alterations in 166 samples) compared to the high-risk group (mutation rate 90.26% with 176 alterations in 195 samples). Fig. 7. [107]Fig. 7 [108]Open in a new tab Mutation profile and drug sensitivity analysis of disulfidptosis risk score in GC. A Top 20 genes in the low risk group in terms of mutation frequency. B Top 20 genes in the high risk group in terms of mutation frequency. C Differences in stromal score, immune score and ESTIMATE score between subtypes. D TMB values in the groups with high and low risk. E Association between risk scores and TMB. F-G K-M curves for OS in the groups with high and low TMB. H–I The relationship between risk groups and drug sensitivity. *P < 0.05; **P < 0.01; ***P < 0.001. Furthermore, our analysis revealed notable differences in TME between the high-risk and low-risk groups ([109]Fig. 7C). Correlation analysis revealed variations in tumor mutation load between the two risk groups, with the low-risk group exhibiting a higher TMB ([110]Fig. 7D). Furthermore, a negative correlation was observed between the two genotypes ([111]Fig. 7E). Kaplan-Meier survival curves indicated superior survival for patients in the high TMB group than those in the low TMB group ([112]Fig. 7F and G), indirectly validating the predictive efficacy of the risk score in immunotherapy. For drug sensitivity, analysis using the “pRRophetic package” revealed that patients in the low-risk group were more sensitive to bexarotene, bicalutamide, bortezomib, dasatinib, and imatinib ([113]Fig. 7H). These drugs have the potential to enhance outcomes in patients with low-risk GC. Conversely, patients in the high-risk group demonstrated greater sensitivity to gemcitabine, gefitinib, bosutinib, sorafenib, and vorinostat, as illustrated in [114]Fig. 7I. Treatment with these agents may be beneficial for patients in the high-risk group. 3.7. Molecular docking analysis Through the use of molecular docking analysis, it was determined that the binding free energies of SLC27A2, PPP1R14A, PHLDA1, LAMC2, KRT80, NMU, GPC3, and SFRP2 proteins with shikonin were calculated to be −7.6, −8.3, −5.9, −6.1, −5.2, −4.6, −7.3, and −5.6 kcal/mol, respectively ([115]Fig. 8). Specifically, THR242 and LYS465 of the GPC3 protein were observed to interact with shikonin, resulting in the formation of a conventional hydrogen bond. Furthermore, LYS247 and VAL458 of the proteins were found to interact with shikonin, leading to the generation of an alkyl interaction. Additionally, THR242 exhibited a pi-sigma interaction with shikonin. Moreover, GLU238, HIS245, LEU246, and GLN461 were observed to interact with shikonin through van der Waals forces. Fig. 8. [116]Fig. 8 [117]Open in a new tab Representation of the interaction between the shikonin and the different proteins. left side, illustrate the interactions in a three-dimensional format, while on the right side, depict the same interaction in a two-dimensional format. A GPC. B KRT8. C LAMC2. D NMU. E PHLDA1. F PPP1R14. G SFRP. H SLC27A2. 3.8. Molecular dynamics simulation analysis The root-mean-square deviation (RMSD) is a metric utilized in assessing the stability of a simulation system by quantifying the deviation of a particular atom's coordinates from a reference structure [[118]24]. A stable RMSD indicates the stability of the corresponding atom, whereas a fluctuating RMSD signifies variability. [119]Fig. 9A illustrates the attainment of equilibrium by the protein ligand after a duration of 35 ns. Fig. 9. [120]Fig. 9 [121]Open in a new tab Comparing post-dynamic data of key shikonin compounds-GPC3 systems, including A (RMSD), B (SASA), C (Rg), D (RMSF), and E (number of intramolecular H-bond plots), over varying ns simulation periods. The solvent-accessible surface area (SASA) is determined through the interplay between van der Waals forces and solvent molecules [[122]25]. SASA of a protein diminishes as its compactness increases, thereby enabling the prediction of structural alterations in the protein. [123]Fig. 9B illustrates a consistent decline in SASA values throughout all ligand recombination simulations, signifying an augmented protein tightness. Radius of gyration (Rg) is employed as a measure of the compactness of the protein structure in the simulation procedure, representing the distance between the center of mass of all atoms and their respective termini within a designated time interval [[124]26]. As depicted in [125]Fig. 9C, the protein Rg values exhibited a consistent decline throughout the entirety of the complex molecular dynamics (MD) simulation, aligning with the observed trend in SASA values. This concurrence suggests an augmentation in the tightness of the protein structure. Root mean square fluctuation (RMSF) is employed to assess the variation of individual atoms in relation to their mean positions, providing insight into the average structural alterations over time and offering an understanding of the flexibility exhibited by different protein regions [[126]27]. Analysis of [127]Fig. 9D reveals substantial fluctuations in the amino acid residues within the protein ligand complex group at the 100–200 stage, suggesting a higher degree of flexibility in the protein during this particular stage. To investigate the interaction between protein and ligand, an initial analysis of hydrogen bonding between the protein and ligand was conducted. As depicted in [128]Fig. 9E, the average count of hydrogen bonds formed between proteins and small molecules was determined to be 2.15, suggesting the presence of hydrogen bond interactions between proteins and ligands. To elucidate the interaction between protein and ligand, the gmx_mmpbsa script ([129]https://jerkwin.github.io/gmxtool/) is employed to determine the binding energy of all protein and ligand complexes in the equilibrium phase. The MMPBSA method decomposes the total binding energy into four distinct components, namely electrostatic interaction, van der Waals interaction, polar solvation, and non-polar solvation interaction, with the latter commonly referred to as SASA. The binding energies of the protein and ligand are presented in [130]Table 1. In the intricate protein-ligand system, the protein and small molecule exhibit a negative binding free energy of −47.493 kJ/mol, respectively, suggesting a stable binding capability. The primary interaction is likely attributed to van der Waals forces. Table 1. Protein ligand MMPBSA analysis. Energy Protein-ligand Van der Waals Energy (KJ/mol) −130.484 Electrostatic energy (kJ/mol) −46.699 Polar solvation energy (KJ/mol) 121.904 Nonpolar solvation Energy (KJ/mol) −19.422 Total Binding Energy (KJ/mol) −74.701 TΔS(KJ/mol) 27.207 Total Binding Free Energy (KJ/mol) −47.493 [131]Open in a new tab 3.9. Validation of the expression levels of three genes in GC We used qRT-PCR to validate the mRNA levels of three genes between GC cell line (MGC-803 and SGC‐7901) and normal gastric cell line (GES-1). The results revealed that OXSM, SLC3A2, and SLC7A11 mRNA levels were dramatically higher in GC cells than in normal gastric cells ([132]Fig. 10). Fig. 10. [133]Fig. 10 [134]Open in a new tab The mRNA levels of OXSM, SLC3A2, and SLC7A11 were quantified using qPCR in the GES-1 immortalized gastric cell line, as well as in two human gastric cancer cell lines, MGC-803 and SGC‐7901. *P < 0.05, **P < 0.01, and ***P < 0.001. 4. Discussion GC is a malignancy with high global prevalence, ranking fifth worldwide in frequency and third in cancer-related mortality, leading to approximately 800,000 deaths annually [[135]28]. Incidence among young individuals is gradually increasing [[136]29]. High-incidence regions include East Asia, Eastern Europe, and South America, while North America and Africa have the lowest rates [[137]30]. Furthermore, males have double the prevalence of GC compared to females [[138]31]. Patients with GC present “three high and three low” attributes: the prevalence, metastasis rate, and mortality being high, while the rates of early diagnosis, radical resection, and 5-year survival are notably low [[139]32]. Surgical intervention is the preferred approach for advanced cases of the disease [[140]33]. However, for patients who are unsuitable for surgery, combining neoadjuvant chemotherapy, targeted therapy, immunotherapy, radiotherapy, and other treatments can potentially extend survival and enhance the overall quality of life [[141]34]. Immunotherapy, including tumor vaccines, monoclonal antibodies, and ICI, holds significant promise in the treatment of GC. Disulfidptosis, a newly identified form of programmed cell death, distinguishes itself from previously investigated forms of cell death such as apoptosis, necrosis, pyroptosis, autophagy, ferroptosis, and cuproptosis [[142]35]. Although disulfidptosis shows promise as a treatment strategy for tumors, including GC, its mechanism of action remains unclear. The limitation of glucose supply to cancer cells with high expression of SLC7A11 triggers substantial intracellular glycoside accumulation, disrupting the oxidoreductase system and inducing cell death [[143]36]. This presents a novel therapeutic strategy for cancer treatment. While the exact molecular mechanisms of disulfidptosis in tumors require further investigation, a tumor prognosis model based on it has been outlined [[144][37], [145][38], [146][39]]. Our study classified patients as belonging to two subtypes of disulfidptosis through unsupervised clustering, revealing distinct patterns of tumor progression and treatment in GC. Significantly increased expression of the disulfidptosis-related genes except NDUFA11, was observed in GC. This elevated expression is associated with notable differences in patient survival, except for RPN1 and SLC3A2. Additionally, we observed that the upregulation of DRG, which is associated with copy number amplification, was prevalent in GC, suggesting that copy number amplification may contribute to the overexpression of this gene. We included 871 patients in this study. Through an analysis of prognosis, we observed that high disulfidptosis subtypes indicated poor prognosis, while low disulfidptosis subtypes exhibited a more favorable prognosis. The GSEA enrichment analysis revealed significant differences between clusters A and B. Notably, cluster B exhibited enrichment in processes positively regulating telomere maintenance, a key aspect in cancer cell perpetuation. This process enables the unlimited proliferation of human malignant cells and the maintenance of other cancer-related characteristics [[147]40]. Cancer cells achieve perpetual replication by activating the telomere maintenance mechanism through the activation of telomerase or telomere lengthening pathway [[148]41]. Consequently, enhancing this process could potentially improve survival rates in cluster B patients. The analysis of signaling pathways revealed that calcium signaling pathways were notably active in cluster A. Calcium ions (Ca^2+) serve as crucial second messengers in numerous physiological and pathological processes. Within tumor cells, Ca^2+ signaling influences proliferation, invasion, and metastasis. Several existing chemotherapeutic agents such as cisplatin, can modify intracellular Ca^2+, prompting tumor cell apoptosis [[149]42]. Furthermore, the calcium signaling pathway significantly contributes to drug resistance in tumors [[150]43]. Consequently, proteins that regulate Ca^2+ such as transmembrane cation channels, hold promise as potential targets for cancer treatment. The investigation revealed that the pentose phosphate pathway (PPP), either directly or indirectly, supplies the necessary reducing capacity for lipid, nucleotide biosynthesis, and antioxidant reactions, supporting cellular survival and proliferation. Cancer cells undergo metabolic reprogramming to redirect glucose toward the PPP, enhancing tumor growth and survival, particularly under stress. Notably, cancer cells exhibit downregulation of the glycolytic pathway during oxidative stress, increasing glucose flux through the PPP. This metabolic shift results in the production of higher levels of NADPH, which serves as an essential antioxidant defense element [[151]44]. The enzyme glucose-6-phosphate dehydrogenase (G6PD) is pivotal in the PPP, and recent studies reveal that tumor suppressor p53 acts as a negative regulator of G6PD. Specifically, p53 directly binds to G6PD, preventing its activation. Silencing the G6PD-encoding gene in the PPP results in the accumulation of p53 protein in cells, leading to cellular senescence in lung cancer and slowing tumor growth in mouse xenotransplantation models [[152]45]. Recent research has demonstrated that cancer cells with high expression of SLC7A11 heavily rely on glucose and the PPP. This reliance results from the rapid conversion of cystine to cysteine, demanding a significant NADPH supply [[153]46]. Additionally, disulfidptosis has been observed in cancer cells with elevated SLC7A11 expression. Further investigation is needed to clarify the link between disulfidptosis and the PPP pathway. Notably, our findings also indicate enrichment in the ubiquitin-mediated proteolytic pathway and cell cycle processes. Ubiquitin signaling plays a vital role in regulating the cell cycle to safeguard genomic integrity and maintain the identity of progeny cells. Disrupting ubiquitination perturbs the cell cycle, leading to tumorigenesis [[154]47]. Thus, targeting the ubiquitin system may represent a promising therapeutic approach for the treatment of cancer. Currently, the understanding of the disulfidptosis mechanism is still in a preliminary stage, requiring further investigation to elucidate the intricate regulatory mechanism of protein ubiquitination in disulfidptosis. To examine the significance of the disulfidptosis in GC, we screened DEGs between subtypes, identifying eight key genes: SLC27A2, PPP1R14A, PHLDA1, LAMC2, KRT80, NMU, GPC3, and SFRP2. Subsequently, a prognostic risk model was developed based on these eight genetic traits, and its accuracy was validated both internally and externally. The disulfidptosis model independently predicts prognosis, tumor immune microenvironment, TMB, tumor stemness, and microsatellite instability. Compared to conventional staging, the disulfidptosis risk score better predicts prognosis. Furthermore, we explored how belonging to high or low-risk groups impacts the response to immunotherapy. TIDE algorithm aids oncologists to identify patients who are more likely to benefit from ICI and has been validated as a tool to predict immune response to glioma treatment through pyroptosis-related genes [[155]48]. TIDE score findings confirm our hypothesis that high-risk groups are associated with a higher prevalence of immunosuppression scores. Thus, the low-risk group may derive therapeutic advantages from immunotherapy. TMB indicates the mutational burden in cancer [[156]23]. Patients with untreated advanced cancer patients present a significant correlation between TMB-H and extended survival [[157]49]. Similarly, our study associates high-risk scores with elevated TMB levels and prolonged survival. Furthermore, we analyzed chemotherapy response in patients with GC based on risk. Low-risk patients were more sensitive to bexarotene, bicalutamide, bortezomib, dasatinib, docetaxel, imatinib, lapatinib, and pazopanib. High-risk patients were sensitive to all-trans retinoic acid, bosutinib, gefitinib, gemcitabine, lenalidomide, methotrexate, sorafenib, and roscovitine. Sorafenib, a tyrosine kinase inhibitor, exhibits significant antitumor properties across various cancer types, including gastric cancer. Sorafenib effectively impedes the progression of GC by inducing ferroptosis [[158]50]. Furthermore, identifying disulfidptosis-associated genes could lead to the development of disulfidptosis inducers, thereby enhancing the predictive capabilities of sorafenib and other chemotherapy drugs. We identified shikonin in Arnebia euchroma (Royle) Johnst as a potent agent with anti-inflammatory, antioxidant, anticancer, and wound healing properties. Additionally, it has beneficial effects on autoimmune disorders [[159]51]. A recent investigation has revealed that shikonin effectively impedes the progression of counteracts ovarian cancer resistance through the induction of ATF3-dependent ferroptosis [[160]52]. Thus, shikonin has demonstrated significant promise in the field of anticancer therapy. Its limited bioavailability led to the development of diverse nanosystems of drug delivery to enhance its efficient delivery and augment its anti-tumor effect [[161]53]. Molecular docking, a computer-based approach extensively employed in drug discovery, plays a crucial role in identifying novel compounds of therapeutic relevance. Compounds exhibiting binding energies below −5.6 kcal/mol to receptor proteins bind robustly to their targets. Most of the eight proteins used to construct the risk model exhibited a binding energy to shikonin below −5.6 kcal/mol [[162]54]. Notably, the docking score between shikonin and PPP1R14A was the highest. As the expression of GPC3 is significantly different in the survival of gastric cancer patients, we further simulated the binding of shiksin-GPC3 by molecular dynamics, and the results revealed that the protein ligand binds stably, and the main interaction is van der Waals interaction. These drugs effectively inhibit tumors, benefiting high/low disulfidptosis risk patients, and potentially acting as disulfidptosis inducers. 5. Conclusion In summary, this study represents a pioneering effort to integrate the disulfidptosis subtype pattern with immune infiltration in individuals diagnosed with GC. The identification of distinct clusters aids in understanding the varying degrees of immune infiltration among patients. Moreover, we have successfully developed a disulfidptosis risk score for prognosis and response to immunotherapy prediction, serving as a novel biomarker. Significant prognostic genes associated with the development of GC were identified, with validated expressions in cancer cells, offering novel targets and avenues for the prognosis and treatment of gastric cancer. To rationally optimize drug use in patients with varying risks, a drug sensitivity analysis was conducted, revealing the potential of shikonin to induce disulfidptosis. Overall, this research offers a robust foundation for understanding the prognostic significance of disulfidptosis-related genes and uncovering novel therapeutic approaches. Nevertheless, our investigation is subject to certain limitations, necessitating analysis of clinical samples, pharmacological experiments, or animal studies to authenticate its results. Availability of data and material The datasets used and analyzed during the current study are available from the corresponding author upon reasonable request. Funding support This work received grants from the National Natural Science Foundation of China (82304796). CRediT authorship contribution statement Xing Liu: Writing – original draft, Methodology, Funding acquisition, Formal analysis, Conceptualization. Jianghong Ou: Writing – review & editing, Visualization, Methodology, Formal analysis, Conceptualization. Declaration of competing interest None Declared. Acknowledgment