Expression and mutation profile of disulfidptosis-related genes (DRGs) in GC

The flow chart listed the detailed procedures to excavate the values of DRIGs in this study, based on quality control and identification of cell subtypes using single-cell data from GC, analyzing the characterization of immune cell-specific genes within the GC microenvironment, correlating them with DRGs, filtering for DRIGs, further selecting genes associated with GC prognosis, constructing a predictive model and validating it using external datasets ACRG and PUMCH (Fig. 1). First, the four DRGs were identified and obtained from the previous literature [18, 30]. The expression of SLC7A11, SLC3A2, RPN1 and NCKAP1 were found to be upregulated in the tumor tissues compared with that in the corresponding normal tissues from the TCGA data (Fig. 2A). The mutation profile of these four DRGs was detected in the GC samples. Among the four genes, SLC7A11, RPN1 and NCKAP1 exhibited minimal alterations, whereas SLC3A2 showed a higher mutation rate in the 431 GC samples (Fig. 2B). Further correlation analysis of expression levels in GC samples confirmed the synergistic effects of the four DRGs in disulfidptosis (Fig. 2C). Additionally, the correlation between the DRGs and immune cells was examined. NCKAP1 showed predominantly positive correlations with M2 macrophages, while displaying negative correlations to regulatory T cells. RPN1, SLC3A2 and SLC7A11 were found to positively activate mast cells, whereas RPN1 and SLC3A2 were negative correlations to memory B cells. And SLC7A11 was negative correlations to regulatory T cells (Fig. 2D).

Fig. 1

The analysis flow of this study

Fig. 2

The expression of disulfidptosis-related genes in gastric cancer (GC). A Violin plots of SLC7A11, SLC3A2, RPN1 and NCKAP1 expression in GC and adjacent normal tissues from the TCGA. B Mutation profile of SLC7A11, SLC3A2, RPN1 and NCKAP1 in 431 GC patients from the TCGA. C Correlation of SLC7A11, SLC3A2, RPN1 and NCKAP1 expression. D Correlation between SLC7A11, SLC3A2, RPN1 and NCKAP1 and immune cell infiltration. E The expression of SLC7A11 in 63 paired GC and adjacent normal tissues. F The IHC of SLC7A11 in GC and adjacent normal tissues (scale bar: 50 μm and 20 μm). Student’s t-test was used to determine statistical significance: ***p < 0.001

To validate the expression of SLC7A11 in GC tissues, a tissue array containing 63 paired samples was utilized. SLC7A11, known to play a role in the ferroptosis process and chemoresistance in various tumors [31, 32]. As shown in Fig. 2E, the IOD of IHC revealed that SLC7A11 was significantly increased in the GC tissues. The representative images of IHC showed that SLC7A11 was elevated in the GC tissues and mainly located in the cytoplasm (Fig. 2F). The correlation between SLC7A11 expression and clinicopathological characteristics revealed that higher expression of SLC7A11 was found to be associated with worse T stage, N stage, and AJCC stage in GC patients from the PUMCH cohort (Table 1).

Table 1 Relationship between clinicopathological characteristics and expression of VPS35, HIF-1α, GLA, CDC37 and SLC7A11 in gastric cancer (n = 63)The distribution and expression of four DRGs in GC at single-cell level

To explore the DRGs in GC at the single-cell level, scRNA-seq data was obtained from the GSE183904 dataset. First, the data was screened with quality control through detected gene numbers, the depth of sequencing, and the gene ratio of mitochondrial and hemoglobin in each specimen (Figures S1A and S1B). Following data normalization, 3000 variable genes were chosen to further analysis (Figure S1C). PCA analysis was used to reduce dimensionality and visualize the data, and the top 30 principal components (PCs) were selected as the input for UMAP analysis (Figure S1D).

Subsequently, UMAP dimensionality reduction resulted in the classification of cells into 37 clusters (Fig. 3A). These 37 clusters were classified into 20 cell types by using SingleR (basophils, CD8+ NKT − like cells, endothelial, ENS glia, goblet cells, ISG expressing immune cells, macrophages, memory B cells, memory CD8+ T cells, mesothelial cells, MUC13_DMBT1 positive cells, myeloid dendritic cells, naive B cells, naive CD4+ T cells, naive CD8+ T cells, neuroendocrine cells, non − classical monocytes, plasma B cells, stromal cells and vascular endothelial cells) (Fig. 3A). Additionally, the specific number and the proportion of each cell type in each GC specimen were calculated (Fig. 3B). The heatmap presented the top 5 marker genes of each cell type (Fig. 3C), while the bubble diagram presented the top 2 marker genes in each cell type (Fig. 3D). The UMAP distribution presented the most representative marker gene for each cell type (Figure S2).

Fig. 3

The analysis of single-cell RNA sequencing. A UMAP of 37 independent cell clusters and 20 cell types were identified by marker genes. B The total number of each cell type and distribution of cell types in different specimens. C The heatmap shows the top 5 marker genes in each annotated cell type. D The bubble diagram shows the top 2 marker genes in each annotated cell type

Examining the single-cell expression profiles of SLC7A11, SLC3A2, RPN1 and NCKAP1 in GC, as depicted in UMAP distribution (Fig. 4A) and violin plots (Fig. 4B), revealed a prevalent disulfidptosis signature. Each cell’s disulfidptosis score, derived from the expression of the four DRGs using the ‘AUCell’ function, was effectively illustrated in the UMAP distribution (Fig. 4C) and visually represented for each cell type through violin plots (Fig. 4D). Of particular note is the heightened disulfidptosis score observed in plasma B cells at the single-cell level. This observation suggests the potential significance of plasma B cells within the gastric cancer tumor microenvironment. These results contribute valuable insights into the distribution and relevance of disulfidptosis across diverse cell types, enriching our understanding of its role in GC.

Fig. 4

The expression of disulfidptosis-related genes in GC at single cell level. A The UMAP distribution of SLC7A11, SLC3A2, RPN1 and NCKAP1 abundance in GC at single cell level. The color variation represents the relative level of gene expression within different cellular clusters, with red and yellow areas indicating higher gene expression, and grey areas indicating lower gene expression or non-detection. B Violin plots of SLC7A11, SLC3A2, RPN1 and NCKAP1 abundance in GC at single cell level. C UMAP distribution of disulfidptosis score of cells by using AUCell function. D Violin plots of disulfidptosis score of each cell type

Screening of differentially expressed DRIGs in GC

To reveal the intrinsic correlations between IRGs and DRGs, a total of 474 differentially expressed DRIGs (459 upregulated DRIGs and 15 downregulated DRIGs) were screened out in GC via the criteria of p < 0.05 and Log2 (fold change) > 1 (Fig. 5A). These DRIGs exhibited enrichment in the establishment of protein localization to organelle in biological process, pigment granule in cellular component, cadherin binding in molecular function, viral carcinogenesis and cGMP-PKG signaling pathway to aggravate the progression of GC through using GO and KEGG analysis (Fig. 5B and C; Table 2). Then, univariate Cox and LASSO analyses were utilized to filter the prominent DRIGs from 474 differentially expressed DRIGs for the prognosis of GC patients (Fig. 5D and E). GLA, HIF-1α, VPS35 and CDC37 were successfully identified from 474 differentially expressed DRIGs to forecast the survival time of GC patients by using the abovementioned methods.

Fig. 5

The screen of disulfidptosis-related immune genes. A The volcanic plot of differentially expressed disulfidptosis-related immune genes in GC. B The GO analysis of differentially expressed disulfidptosis-related immune genes in GC. C The KEGG analysis of differentially expressed disulfidptosis-related immune genes in GC. D LASSO Cox regression analysis of the association between deviance and log(λ). E LASSO Cox regression analysis of the association between coefficients of genes and log(λ)

Table 2 The GO and KEGG enrichment analysis of disulfidptosis-related immune genesConstruction and validation of the DRIG signature

The formulation of the risk score formula and calculation of gene coefficients were performed using multivariate Cox analysis. The formula to score GC patients was obtained as follows: Risk score = (-0.25455 × GLA) + (0.29278 × HIF-1α) + (0.62703 × VPS35) + (-0.52510 × CDC37). Then, patients were classified into high-risk and low-risk groups based on the median value of all patients’ risk scores.

As shown in Fig. 6A and B, the proportion of death was higher in the GC patients with high risk scores. The heatmap displayed the expression of four DRIGs in two groups, in which HIF-1α and VPS35 were increased in the high-risk group while GLA and CDC37 were elevated in the low-risk group (Fig. 6C). In addition, KM analysis showed that patients with high scores had shorter survival times compared with patients with low scores (Fig. 6D). The risk signature was then utilized to predict the survival probability of GC patients. The accuracy of signature in predicting the 1-, 3-, and 5-year probability was 0.686, 0.647 and 0.632 in GC patients, respectively (Fig. 6E). A nomogram consisting of a risk model and various clinicopathological features was successfully established to accurately forecast the survival time of GC patients (Fig. 6F). The calibration curve was used to verify the validity and accuracy of the nomogram for the predictive probability of 1-, 3- and 5-year in GC patients (Fig. 6G). Further, a webserver (https://pumc.shinyapps.io/GastricCancer/) was constructed to make full use of the nomogram in prognostic prediction for GC patients. The quick response code provided a convenient entrance using the online tool in clinical practice for physicians (Fig. 6H). For the stringency of the study, the Asian Cancer Research Group cohort (GSE62254) was utilized to validate the clinical values of the signature. The signature not only well distinguished the overall survival (OS) and disease-free survival (DFS) between the two groups, but also accurately predicted the 1-, 3- and 5-year OS and DFS of GC patients (Fig. 6I and J). In summary, we preliminarily investigated and validated the potential use of DRIGs in predicting prognosis of GC patients in clinical translation.

Fig. 6

Construction and validation of prognostic models based on disulfidptosis-related immune genes. A Distribution of risk score. B The survival status and survival time of GC patients ranked by risk score. C The heatmap of GLA, HIF-1 A, VPS35 and CDC37 in two groups. D Kaplan-Meier analysis between high-risk group and low-risk group. E Time-dependent ROC curve of risk score predicting the 1-, 3-, and 5-year overall survival. F Details of the nomogram. G The calibration curve for predicting 1-, 3-, and 5-year overall survival. H The quick response code of online dynamic nomogram I The OS discrepancy between the high-risk and low-risk groups and the ROC curve of risk score predicting the 1-, 3-, and 5-year OS using the GSE62254 dataset. J The DFS discrepancy between the high-risk and low-risk groups and the ROC curve of risk score predicting the 1-, 3-, and 5-year DFS using the GSE62254 dataset

Gene Set Enrichment Analysis (GSEA) and immune cell infiltration in two groups

The discrepancy in survival time between high-risk and the low-risk groups suggests a significant heterogeneity of the genome in the two groups. To further investigate the underlying mechanisms contributing to this heterogeneity, GSEA analysis was used to further explore the gaps in underlying mechanisms between the two groups. As shown in Fig. 7A, multiple classical signaling pathways of tumors were involved in the high-risk group, including the activation of oxidative phosphorylation, antigen processing and presentation, DNA replication and cell cycle and inhibition of focal adhesion, calcium signaling pathway, adherens junction, Wnt signaling pathway, pathways in cancer, MAPK pathway, PPAR pathway, TGF beta pathway, mTOR pathway, Toll-like receptor pathway, JAK-STAT pathway and P53 pathway. The involvement of these signaling pathways resulted in the exacerbation of GC progression. In addition, the correlation between risk score and immune cell infiltration was also evaluated. The result revealed that risk score was positively correlated to the activation of dendritic cells and resting of mast cells, while negatively related to the plasma cells, CD8+ T cells and regulatory T cells (Fig. 7B). Furthermore, the correlation between 4 DRIGs and immune cell infiltration was also investigated. As shown in Fig. 7C, CDC37 was positively correlated to the activation of mast cells, M0 and M1 macrophages. GLA was positively correlated to the activation of memory CD4+ T cells while negatively related to memory B cells, HIF-1α was positively related to the activation of memory CD4+ T cells while negatively correlated to the activation of NK cells, VPS35 was negatively related to regulatory T cells while positively correlated to the activation of memory CD4+ T cells and M2 macrophages. The data revealed that the 4 DIRGs (CDC37, GLA, HIF-1α and VPS35) that constructed the risk signature profoundly regulated the immune microenvironment of GC.

Fig. 7

The GSEA results and correlation between the risk signature and immune cell infiltrations. A The results of GSEA analysis between high-risk group and low-risk group. B The correlation between the risk score and infiltrated immune cells. C The correlation between the GLA, HIF-1 A, VPS35 and CDC37 and infiltrated immune cells

Somatic mutation profile of 4 DIRGs and two groups

First, the mutation of CDC37, GLA, HIF-1α and VPS35 was investigated in the 431 GC specimens. The results indicated that 4 DIRGs were generally stable in GC and only 5.1% of samples exhibiting mutations in CDC37, GLA, HIF-1α and VPS35 (Fig. 8A). The results of somatic mutation in 431 GC specimens revealed that TTN, TP53, MUC16, ARID1A, LRP1B, CSMD3, SYNE1, FAT4, FLG and PCLO were among the top 10 genes with the highest mutation frequencies (Fig. 8B). Then, the mutation profile in the two groups was also detected. TTN, TP53, MUC16, ARID1A, LRP1B, CSMD3 and SYNE1 were the most frequently mutated genes in both groups. However, PIK3CA, OBSCN and PCLO were the three most frequently mutated genes in the low-risk group (Fig. 8C), while SPTA1, FAT4 and FLG were the frequently mutated genes in the high-risk group (Fig. 8D).

Fig. 8

Somatic mutations in the entire GC samples and two groups. A Mutation profile of GLA, HIF-1 A, VPS35 and CDC37 in 431 GC patients. B The mutation profile in low-risk group. C The mutation profile in high-risk group. D The mutation profile in entire GC samples

The clinical application of risk signature in GC

A cruel situation is that about 80% of hospitalized patients were initially diagnosed with locally advanced or metastatic GC in China [8, 33]. This phenomenon required postoperative chemotherapy and immunotherapy in treatment for GC patients. Therefore, we explored the relationship between 4 DIRGs and CD274 expression, as well as the response to classical chemotherapy drugs of GC in the two groups. As shown in Fig. 9A, the expression of CD274 was usually consistently elevated in parallel with the increased expression of CDC37, GLA, HIF-1α and VPS35 in GC. This phenomenon indirectly reflected that 4 DRIGs may be associated to the immune evasion of GC cells. In addition, patients with higher risk scores were generally more resistant to 5 − Fluorouracil, docetaxel, erlotinib, methotrexate and paclitaxel treatments in GC (Fig. 9B). These findings highlight a potential risk signature in guiding therapy selection for GC patients in the future.

Fig. 9

Immune infiltration level and drug sensitivity analysis based on the risk model. A The correlation between the GLA, HIF-1 A, VPS35 and CDC37 and CD274. B The results of drug sensitivity analysis between high-risk and low-risk groups

Validation of 4 DIRGs in GC using PUMCH cohort

The PUMCH cohort were then used to validate the expression of CDC37, GLA, HIF-1α and VPS35, as well as their relationships to SLC7A11 and clinicopathological characteristics in GC for the abovementioned multiple explorations in bioinformatics results. The expression of CDC37, GLA, HIF-1α and VPS35 was detected in the tissue array by IHC staining. The results of IHC demonstrated that the expression of CDC37, GLA, HIF-1α and VPS35 was upregulated in the tumor tissues compared with corresponding normal tissues (Fig. 10A and B). Additionally, separate IHC staining of the same GC tissue demonstrated that the expression of SLC7A11 increased with elevated level of CDC37, GLA, HIF-1α or VPS35 (Fig. 10C). Analysis of data from 63 GC tissues showed a positive correlation between CDC37, GLA, HIF-1α, VPS35 and SLC7A11 (Fig. 10D). Furthermore, patients with higher expressions of VPS35 and HIF-1α frequently associated with larger tumor size; while higher expression of GLA was indicative of worse T stage, N stage and AJCC stage. Increased expressions of CDC37 and HIF-1α often predicted worse N stage and advanced clinical stage of GC (Table 1). The close relationships between the 4 DIRGs and clinicopathological characteristics further support the validity and accuracy of the disulfidptosis-related immune genes signature in predicting the survival time of GC patients.

Fig. 10

The expression of GLA, HIF-1 A, VPS35, CDC37 and SLC7A11 in GC. A The expression of GLA, HIF-1 A, VPS35 and CDC37 in 63 paired GC and adjacent normal tissues. B The IHC of GLA, HIF-1 A, VPS35 and CDC37 in GC and adjacent normal tissues (scale bar: 50 μm and 20 μm). C The representative IHC images of GLA, HIF-1 A, VPS35, CDC37 and SLC7A11 in the same GC tissues (scale bar: 50 μm). D The correlation between GLA, HIF-1 A, VPS35, CDC37 and SLC7A11 in the GC tissues. Student’s t-test was used to determine statistical significance: ***p < 0.001