Identification of disulfidptosis-related subtypes, characterization of tumor microenvironment infiltration, and development of a prognosis model in breast cancer

Introduction: Breast cancer (BC) is now the most common type of cancer in women. Disulfidptosis is a new regulation of cell death (RCD). RCD dysregulation is causally linked to cancer. However, the comprehensive relationship between disulfidptosis and BC remains unknown. This study aimed to explore the predictive value of disulfidptosis-related genes (DRGs) in BC and their relationship with the TME.

Methods: This study obtained 11 disulfidptosis genes (DGs) from previous research by Gan et al. RNA sequencing data of BC were downloaded from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus database (GEO) databases. First, we examined the effect of DG gene mutations and copy number changes on the overall survival of breast cancer samples. We then used the expression profile data of 11 DGs and survival data for consensus clustering, and BC patients were divided into two clusters. Survival analysis, gene set variation analysis (GSVA) and ss GSEA were used to compare the differences between them. Subsequently, DRGs were identified between the clusters used to perform Cox regression and least absolute shrinkage and selection operator regression (LASSO) analyses to construct a prognosis model. Finally, the immune cell infiltration pattern, immunotherapy response, and drug sensitivity of the two subtypes were analyzed. CCK-8 and a colony assay obtained by knocking down genes and gene sequencing were used to validate the model.

Result: Two DG clusters were identified based on the expression of 11DGs. Then, 225 DRGs were identified between them. RS, composed of six genes, showed a significant relationship with survival, immune cell infiltration, clinical characteristics, immune checkpoints, immunotherapy response, and drug sensitivity. Low-RS shows a better prognosis and higher immunotherapy response than high-RS. A nomogram with perfect stability constructed using signature and clinical characteristics can predict the survival of each patient. CCK-8 and colony assay obtained by knocking down genes have demonstrated that the knockdown of high-risk genes in the RS model significantly inhibited cell proliferation.

Discussion: This study elucidates the potential relationship between disulfidptosis-related genes and breast cancer and provides new guidance for treating breast cancer.

1 Introduction

In recent years, the incidence of breast cancer in women has continued to grow at a rate of 0.5% per year. By 2022, breast cancer has surpassed lung cancer and has become the world’s highest incidence of tumors. According to the latest estimates, there will be 297,790 new breast cancer cases and 43,170 deaths among women in the United States in 2023, making it the leading cause of cancer death among women aged 20-49 (1). Many treatment options for BC have been developed, including surgery, chemotherapy, endocrine therapy, immune antibody (trastuzumab) therapy, and radiotherapy, based on disease stage and pathological characteristics (2, 3). Although the prognosis of breast cancer has improved significantly in the past few decades, even if patients pass standard diagnosis and treatment, 20–30% of BC patients still have distant metastasis, accounting for about 90% of all breast cancer deaths (2).

As a highly complex and heterogeneous disease with different molecular spectrums, breast cancer limits the broad application of classification and standard treatment to some extent. Furthermore, it is difficult to predict the prognosis of BC (2, 4). Therefore, we urgently need to explore the characteristics of the high-risk population of breast cancer patients to obtain the key markers of prognosis and to find potential therapeutic targets to designate individualized treatment to improve the prognosis of patients.

SLC7A11-mediated cystine reduction to cysteine highly depends on the reduced nicotinamide adenine dinucleotide phosphate (NADPH) generated by the glucose–pentose phosphate pathway. Under glucose starvation, the NADPH in SLC7A11-high-expression cells is consumed in large quantities, and the abnormal accumulation of disulfides, such as those in cystine, induces disulfide stress. This causes actin filaments to aggregate and contract rapidly, peeling off the plasma membrane before apparent cell death, called disulfidptosis. During this process, SLC7A11 and SLC3A2, which encode the SLC7A11 chaperone protein, mediate the reduction of ingested cystine to cysteine. Next, the WAVE regulatory complex (WRC) can activate seven subunits of actin-related protein 2 and 3 (Arp2/3) complexes to promote actin polymerization and plate pseudopod formation. This produces a branched cortical actin network under the plasma membrane, thereby stripping actin filaments from the plasma membrane. Nck-related protein 1 encoded by NCKAP1 is a WRC component. Its deletion will reduce the protein levels of other components in the WRC, including WAVE-2, CYFIP 1, Abi 2, and HSPC 300, and inhibit disulfidptosis. Similarly, knocking out the other four components of the WRC will also repress disulfidptosis. Finally, Rac can activate WRC to promote plate pseudopod formation and disulfidptosis. Knockout of the RPN1 gene, which encodes an N-oligosaccharide transferase in the endoplasmic reticulum, will offer UMRC6 cells stronger resistance to disulfidptosis.

However, DRGs have not been demonstrated in the survival prognosis, tumor immune microenvironment, TMB, immunity, and clinical treatment of BC patients. We still lack direct evidence of the predictive ability of DRGs for BC prognosis and immunotherapy. This study aimed to construct a risk score (RS) model based on DRGs to predict the prognosis of BC patients. In this study, TCGA and GEO databases were used to obtain DRGs by analyzing the difference in consensus clustering of DGs. The RS model of DRGs was built using the LASSO-Cox method. Further, we verified that this feature could be used as a reliable, independent predictor of prognosis and immune-sensitive response and could predict the prognosis of BC patients. Through this model, patients were divided into high-risk and low-risk groups. The differences in survival outcomes, tumor immune microenvironment, and immunotherapy response of BC patients were analyzed. According to the differences between groups, drug sensitivity studies were conducted to find sensitive drugs in different populations, and individualized intervention was implemented to improve the prognosis of BC patients.

2 Methods

2.1 Data collection and processing

We downloaded and organized 1,180 samples (1,081 tumor and 99 normal samples) from the TCGA database using the TCGAbiolinks and SummaryExperiment package in R v.4.2.2 (Supplementary Table 1). We searched the GEO database with “breast cancer” and “survival information” as keywords and screened according to the following criteria: first, all samples were from humans; second, all data sets included matched cancer tissue samples and clinical information; and third, the data set contained at least 200 samples. We downloaded the GSE12276 and GSE20685 gene expression profile, GPL570 platform annotation information, and clinical information using GEOquery, the data transmission tool of the GDC application. Gene expression profiles were standardized using the scaling method provided in the limma R package (5). TCGA and GEO gene expression profile information and clinical data are publicly available and open access. Therefore, no ethical issues were involved. Clinicopathologic features include gender, age, T (tumor) stage, N (lymph node metastasis) status, M (distant metastasis) status, tumor grade, survival status, and survival time. The tumor mutation data were derived from TCGA [GDC (cancer.gov)], and the gene copy number was downloaded from the Xena database (UCSC Xena). We obtained 11 disulfidptosis genes (DGs) from a previous study (6). The workflow of the current study is shown in Figure 1.

FIGURE 1

Figure 1 Workflow diagram.

2.2 Survival, TMB, and CNV difference analysis between DGs

We used the Limma package using wilcox.test to test the difference in DG expression between tumor and normal samples (p < 0.05). R package survival and survminer were used for the survival analysis of DGs. DGs were divided into high- and low-expression groups according to gene expression. The survival differences between the two groups were compared (p < 0.05). Tumor mutation burden (TMB) refers to the number of somatic nonsynonymous mutations or all mutations per megabase in the gene region detected by whole exome sequencing or targeted sequencing in a tumor sample (7). Mutation data were downloaded from the Cancer Genome Atlas Breast Cancer (TCGA-BRCA) collection using the GDCquery package with 988 data downloads. According to the data.category = Simple Nucleotide Variation, data.type = Masked Somatic Mutation. The mutation data were analyzed with the R package maftools (8). The mutation probability of DGs in each tumor sample was calculated and analyzed. CNV (Copy Number Variation) is the increase or decrease in the copy number of genomic fragments caused by genome rearrangement. Xena downloaded and collated copy number variation data and used R-package RCircos to visualize the CNV of DGs.

2.3 Establishment of the disulfidptosis gene cluster

We used the R software package “ConsensusCluster Plus” for consensus clustering analysis and identified clusters of BC patients based on DGs. We set the cluster count (k) between two and nine and selected the optimal k value based on the sum of squared errors (SSE) inflection point. The stability of the DG group was verified by the PCA algorithm. Additionally, Kaplan-Meier survival analysis evaluated the OS of different DG clusters, and the heat map showed the degree of difference in DGs between the two groups.

2.4 TME infiltration and functional enrichment analysis of different clusters

Gene set variation analysis (GSVA) is a particular gene set enrichment method. This method works on single samples and enables pathway-centric analyses of molecular data by performing a conceptually simple but powerful change in the functional unit of analysis from genes to gene sets. We used the GSVA and GSEABase packages to evaluate the pathways enriched in groups A and B. We attempted to explain the reasons for the difference in survival between the groups from the bioinformatics perspective (9). Single sample gene set enrichment analysis (ssGSEA) was used to quantitatively analyze the immune infiltration of the overall sample and to observe the difference in component immune infiltration (10). We used the “limma” R package to analyze the DRGs between the two disulfidptosis clusters (p < 0.05, |logFC| = 0.585). Finally, we used the “ggplot2” package. The gene ontology (GO) analysis of DRGs was performed, and the histogram, bubble diagram, and circle diagram were drawn to explore its function and biological process. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, histogram, and bubble diagram were drawn to explore its function and biological processes.

2.5 Identification of DRGs and construction of the prognostic signature

Before constructing the risk prediction model, we screened the features. First, a univariate Cox model was used to explore the relationship between 225 DRGs and OS of patients. A total of 114 single-factor DRGs related to BC prognosis were obtained. The least absolute shrinkage and selection operator (LASSO) was used to avoid overfitting in the TCGA training cohort (11). The prognostic DRG model was established using multivariate Cox regression analysis and the step Akaike information standard (stepAIC) value. The RS for each patient was calculated by combining the expression of each gene (Ei), LASSO coefficients (Li), and RS = Σⁿ_i Ei * Li. The patients were divided into high-risk and low-risk groups according to the median RS. Finally, the Kaplan–Meier survival curve was used to analyze the difference in overall survival between the two groups.

The sensitivity and specificity of the prognostic indicators were evaluated using the receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC). A bootstrap method based on 1,000 resamplings was used to obtain the test set (12, 13) to verify the effectiveness of the prediction model. The training set was a combined dataset established using TCGA and GEO based on common genes. In the test set, the AUC of the prognostic model was calculated using the R package “riskRegression.” Subsequently, the stability of the ROC model was tested in the separate TCGA and GEO datasets, and the TCGA and GEO joint datasets merged based on common genes. Analyze the relationship between RSs and clinical factors, verify the effectiveness of risk markers, and further compare the survival prediction ability of prognostic factors. The independence of the prognostic model was verified using univariate and multivariate Cox regression analyses by comparing the clinical characteristics of the patients. At the same time, the nomogram was constructed with the Cox regression coefficient of the package “rms,” and the calibration curve was drawn.

2.6 Establishment of DRG cluster

We used the “ConsensusCluster Plus” package for consistency clustering analysis and identified clusters of BC patients based on DRGs. We set the cluster count (k) between two and nine and selected the optimal k value based on the sum of the squared errors (SSE) inflection point. The stability of the DG group was verified using a PCA algorithm. Additionally, Kaplan–Meier survival analysis evaluated the OS of different DG clusters, and the heat map showed the difference in DGs between groups.

2.7 Different tumor immune microenvironment patterns with RS

We used the CIBERSORT algorithm to calculate the proportion of tumor-infiltrating immune cells (14). The difference in the proportion of tumor-infiltrating immune cells between the high- and low-risk groups was compared. The ESTIMATE algorithm was used to evaluate the differences in immune, stromal, and tumor purity scores between the high- and low-risk groups. The tumor mutation burden was assessed based on whether the sample was considered high or low risk. Additionally, we used the Maftools package to perform somatic mutation analysis on breast cancer patients to view and analyze somatic mutation data. We also studied the relationship between BC RSs and cancer stem cells.

2.8 The role of RS based on DRGs in predicting drug sensitivity and clinical immune efficacy

For a long time, the research and development of new drugs have been a hot spot in breast cancer treatment. We used the pRRophetic package to calculate the half inhibitory concentration (IC50) of commonly used drugs in breast cancer patients. We screened out potential drugs with sensitivity differences between the two groups according to the risk level of breast cancer patients (p < 0.001) and drew a box plot. We used ggplot2, ggpubr, limma, and reforme2 R software packages to analyze the statistical differences in the expression levels of 79 common ICI-related immunosuppressive molecules (15). Tumor Immune Dysfunction and Exclusion (TIDE) [Tumor Immune Dysfunction and Exclusion (TIDE) (harvard.edu)] is a simple method to predict the immune escape of patients based on the evaluation of the tumor microenvironment using gene expression profiles. Patients with high TIDE scores have a high chance of antitumor immune escape (16). We obtained information on the immune escape ability after submitting the transcriptome data of TCGA-BRCA patients to the website. TCGA is a quantitative scoring scheme developed by developers using machine learning: it is a better predictor of anti-cytotoxic T lymphocyte antigen 4 (CTLA-4) and anti-programmed cell death protein 1 (anti-PD-1) antibody responses. We downloaded BC immunophenotype score (IPS) from TCIA database (https://tcia.at/). In order to predict the sensitivity of immunotherapy, we compared the IPS of high and low risk groups in different immunotherapy decisions (17).

2.9 Cell lines, cell culture, cell transfection, and real-time quantitative PCR

The Tumor Cell Line Comprehensive Analysis Database (DepMap Portal) was utilized to screen for cell lines for further experimental validation (18). Breast cancer cell lines (MDA-MB-468) were obtained from Sichuan Huijixin Biotechnology Co., Ltd. The MDA-MB-468 cells were grown in Dulbecco’s Modified Eagle Medium (DMEM) culture medium, supplemented with 10% Fetal Bovine Serum (FBS) in a standard humidified incubator with 5% CO2 at 37°C. The TMEM45A and SHCBP1 specific short hairpin RNAs (shRNAs) were synthesized from Chengdu Youkangjianxing Biotechnology Co., Ltd.

The sequences of shRNAs are as the following: sh- TMEM45A -1: 5′- TGCTGTTGACAGTGAGCGCGGTTAAAGTATTTGAATTTAATAGTGAAGCCACAGATGTATTAAATTCAAATACTTTAACCATGCCTACTGCCTCGGA-3′; sh- TMEM45A -2: 5′- TGCTGTTGACAGTGAGCGCGGTGTACAAAGAGTATTCTGATAGTGAAGCCACAGATGTATCAGAATACTCTTTGTACACCATGCCTACTGCCTCGGA′. sh- SHCBP1 -1: 5′-TGCTGTTGACAGTGAGCGCCACATTGATTTTTCAATTGAATAGTGAAGCCACAGATGTATTCAATTGAAAAATCAATGTGATGCCTACTGCCTCGGA-3′; sh- SHCBP1 -2: 5′- TGCTGTTGACAGTGAGCGCCAGCCAAATGTTGATATTAAATAGTGAAGCCACAGATGTATTTAATATCAACATTTGGCTGATGCCTACTGCCTCGGA -3′.

The knockdown efficiency was evaluated using real-time quantitative PCR (RT-qPCR) after 48 h transfection. The primer sequences used in the experiment are as follows. For the TMEM45A gene, the primers include qpcr-TMEM45A-F (forward primer): TTGGATGCCCACACTATGA and qpcr-TMEM45A-R (reverse primer): TCCATGGTCAAGGAGTTACA. For the SHCBP1 gene, the primers consist of qpcr-SHCBP1-F (forward primer): CTGGAGTTACAGAAGGATGGTG and qpcr-SHCBP1-R (reverse primer): CCATAGAAGCCTGTGGAATGT. After the knockdown of TMEM45A and SHCBP1 in cell line MDA-MB-468, total mRNA from cells was extracted with TRIzol reagent (TaKaRa, Japan). Then, concentration and purity were evaluated by Nanodrop 2000 (Thermo Fisher, USA). After the RNA was reversely transcribed into cDNA with PrimeScript RT kit (TaKaRa, Japan) according to the instructions, SYBR Premix Ex Taq TM kit (TaKaRa, Japan) was applied for RT-qPCR, with β-actin as the endogenous control gene. The RT-qPCR amplification instrument (ABI StepOne Plus) was used to detect the SYBR Green fluorescence signal level after each amplification cycle. Data processing was performed using GraphPad Prism 10.0.0, and T-test was conducted to compare the experimental group to the control group. This supplementary material has been reflected in the preceding text.

2.10 Cell proliferation assay

For three days, proliferation assays were conducted daily on BC cells in 96-well plates using Cell Counting Kit 8 (CCK8) reagent (Beyotime, China). Incubation occurred at 37°C for 2 h, and the plate was analyzed with a microplate reader at 450 nm to measure absorbance.

2.11 Colony formation assays

Approximately 500 cells per well were seeded into a 6-well culture plate and incubated at 37°C for two weeks. After washing with PBS twice, cells were fixed with 4% paraformaldehyde for 15 min and then dyed with crystal violet. Each experiment was repeated three times. ImageJ was used for image analysis to convert images into cellular count data (19). The acquired counts were normalized by dividing them by the corresponding cell count in the control group, yielding percentage data. Data and image processing were performed using GraphPad Prism 10.0.0. The statistical analysis consisted of a t-test conducted on three replicate datasets, comparing the experimental and control groups.

2.12 RNA-seq

Cells were used for RNA sequencing after the knockdown of TMEM45A and SHCBP1 in cell line MDA-MB-468. Approximately 2 μg of total RNA was extracted from each specimen and pretreated with Epicentre Ribo-zeroTM rRNA Removal Kit. Then, the RNA expression profile library was constructed in line with the manufacturer’s protocol of NEBNext R Ultratdirectional RNA Library Prep Kit(NEB, USA). The steps are as follows: First, RNA was lysed into small fragments after being treated with NEBNext first strand synthesis reaction buffer at high-temperature treatment, and the first strand cDNA was synthesized using random hexmer primers and M-MULv reverse transcriptase. The second strand dsDNA was then obtained, and the fragment residues were converted into blunt ends by exonuclease or polymerase. Subsequently, the 3’end of each dsDNA fragment was adenylated and connected to the NEBNext adapter with a hairpin structure. After purification using the AMPure XP system (Beckman Coulter, Beverly, USA), 150–200 bp DNA fragments were obtained and sequenced using HiSeq 2500 (Illumina, CA, USA).

2.13 RNA-seq data processing and analysis

FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was utilized to check the sequencing quality of all the sample data trimmed using the FASTX Toolkit. The sequencing reads against the human assembly GRCh37 were mapped using TopHat (v 2.0.9). We perform differential analysis using the gene expression matrix in counts format. We employed the R package edgeR to perform differential analysis between snSHCBP1 and NC, as well as snTMEM45A and NC, using predetermined criteria (fold change > 1, padj < 0.05). Following this analysis, we created volcano plots to visualize the differential expression of genes. We conducted an intersection analysis to identify genes that exhibited consistent differential expression in both replicates. Subsequently, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (p-value < 0.05) on the genes found in the intersection.

3 Results

3.1 11DGs in breast cancer: expression, genetic variants, and prognostic values

We analyzed the differences between 11 DGs in breast cancer and normal breast tissue samples. Except for NCKAP1, other genes significantly differed between tumor and normal groups (p < 0.05). The expression levels of SLC3A2, RPN1, BRK1, ACTR2, ACTR3, RAC1, SLC7A11, and WASF2 in the tumor samples were higher than those in the normal samples. On the other hand, the expression levels of WASF2, CYFIP1, and ABI2 in the normal samples were higher (Figure 2A). Of the 911 breast cancer samples of TCGA, 33 samples had DG mutations. The genes with the highest mutation frequencies were NCKAP1 and CYFIP1. Still, no base mutation was found in WASF2 (Figure 2B). The location and copy number changes of DGs on chromosomes are shown in Figure 2C. Among them, SLC3A2, BRK1, ABI2, ACTR2, NCKAP1, RAC1, RPN1, SLC7A11, and ACTR3 are upregulated in breast cancer, while CYFIP1 and WASF2 were downregulated. COX regression analysis and K-M survival analysis were performed on the expression of 11 DGs and the survival time of patients, and the tumor samples were divided into high- and low-risk groups according to the median value of a single gene. It can be concluded that SLC7A11, SLC3A2, RPN1, NCKAP1, BRK1, ACTR2, ACTR3, and RAC1 are single-factor genes that can predict the survival of patients (p < 0.05, Supplementary Table 2). The survival analysis was statistically significant, and the gene expression was negatively correlated with the survival time (Figure 2D). The 11 DGs positively regulated each other (Figure 2E).

FIGURE 2

Figure 2 The analysis of 11 DGs' expression and association in the TCGA cohort. (A) The expression of the 11 DGs in BC tissues and healthy breast tissues (*p<0.05; ***p< 0.001). (B) Data on the frequency of DGs' mutations for 991 BC patients. (C) The sites of CNV variation in DGs on the 23 chromosomes. (D) The relationship between 8 DGs and overall s survival. (E) The interactions between DGs in BC (the red and blue strings denote positive and negative correlation, respectively; the intensity of the correlation is indicated by the color shades).

3.2 Clusters of DGs identified in breast cancer

We performed a consistent cluster analysis of the expression levels of 11 DGs to explore their roles in the occurrence and development of breast cancer. Among the clustering variables, k = 2 had better clustering stability, the highest intra-group correlation, and the lowest inter-group correlation (Figures 3A–C). Therefore, we divided all tumor samples into clusters A (n = 339) and B (n = 850). PCA showed a significant interval between the two groups (Figure 3D). Survival analysis between groups A and B showed that the prognosis of group A was significantly better than that of group B, and the survival difference between the two groups was statistically significant (p < 0.05, Figure 3E). Additionally, the heat map between the two AB groups showed significant differences in the expression of 11 DGs between the two groups, especially between the SLC7A11 groups (Figure 3F).

FIGURE 3

Figure 3 Biological and clinicopathological characteristics of DG subtypes. (A) The consensus matrix's heatmap of two clusters (κ = 2). (B) tracking plot. (C) Consensus CDF. (D) A considerable transcriptome divergence between the two subtypes is seen by PCA analysis. (E) Subtype-specific Kaplan-Meier OS curves. (F) DGs expression levels and clinicopathological traits vary across subtypes.

3.3 TME infiltration and functional enrichment analysis of different clusters

The samples in the A and B groups were analyzed using GSVA. From Figure 4A, it can be seen that group B was more active in arachidonic acid metabolism and taurine and hypotaurine metabolism pathways than group A. At the same time, group A was active in nonhomologous end connections, ubiquitin-mediated proteolysis, and other pathways. We performed a single sample gene set enrichment analysis (ssGSEA) on groups A and B. From Figure 4B, we can see Activated.B.cell, Activated.CD8.T.cell, CD56dim.natural.killer.cell, Eosinophil, MDSC, Macrophage, Mast.cell, Monocyte. Natural.killer.cell, Neutrophil, Plasmacytoid.dendritic.cell, T.follicular.helper.cell, Type.1.T.helper.cell, and Type.17.T.helper.cell were highly expressed in Cluster A and Activated.CD4.T.cell, Activated.dendritic.cell, Gamma.delta.T.cell, Immature.B.cell, Immature.dendritic.cell, Regulatory.T.cell, and Type.2.T.helper.cell were highly expressed in cluster B. We found 225 differential genes between groups A and B to explore the differences in biological processes between groups A and B according to the DG grouping.

FIGURE 4

Figure 4 Disulfidptosis subtypes linked to TME invasion. (A) GSVA of two disulfidptosis subtype-related cellular pathways (Red means activated and blue means inhibited). (B) Correlations between immune cell infiltration levels in the two subtypes associated with disulfidptosis. (C, D) The GO function enrichment analyses. (E, F) The KEGG function enrichment analyses.

We performed GO and KEGG enrichment analyses on these differential genes. As shown in Figures 4C, D, DRGs are mainly enriched in the nuclear division, chromosome segregation, and nuclear chromosome segregation in biological processes (BP). In addition, DRGs were primarily enriched in the spindle, chromosomal region, and condensed chromosome. Regarding molecular function, DRGs were mainly enriched in microtubule binding, tubulin binding, and ATP hydrolysis activity. We also performed a KEGG enrichment analysis to explore the differential genes in the pathway. The most enriched pathways were in the cell cycle. Cytokine−cytokine receptor interaction and human papillomavirus infection were the most significant enrichments in the cell cycle, PPAR signaling pathway, and ECM−receptor interaction (Figures 4E, F).

3.4 Clusters of DRGs identified in breast cancer

Based on a consensus cluster analysis of 225 DRGs, the relationship between disulfidptosis and BRCA subtypes was explored. K = 3 is the appropriate choice for the most stable aggregation (Figures 5A–D). Therefore, we divided all tumor samples into three subgroups: cluster A (n = 525), cluster B (n = 379), and cluster C (n = 285). 11DGs exhibit significant differences among clusters A, B, and C (Figure 5E). Figure 5F shows that the prognosis of group A is the best, followed by group B, and the prognosis of group C is relatively poor. From Figure 5G, it can be seen that there were significant differences in the expression levels of DRGs between different clusters. The samples of DRG cluster A were mostly in cluster A of 11 DGs, and the survival prognosis of these two clusters was significantly better than that of other groups.

FIGURE 5

Figure 5 Biological and clinicopathological characteristics of DRG subtypes. (A) The consensus matrix's heatmap of two clusters (κ = 3). (B) Consensus CDF (C) Delta area. (D) tracking plot Consensus CDF. (E) DRGs expression levels and clinicopathological traits vary across subtypes. (F) Subtype-specific Kaplan-Meier OS curves. (G) DRGs expression levels and clinicopathological traits vary across subtypes.

3.5 Creating and confirming the predictive RS

We used univariate Cox regression analysis to extract 114 DRGs associated with BC prognosis in the preliminary screening (Supplementary Table 3). Using the LASSO regression algorithm, 14 BC-related genes were identified based on the minimum partial likelihood of the best λ value and deviation (Figures 6A, B, Supplementary Table 4). Multivariate Cox regression analysis was performed on these 14 genes, and a risk model consisting of six genes was obtained (Figure 6C). Among them, SHCBP and TMEM45A were molecules that improved prognosis. PIGR, IGLV6-57, TCN1, and GFRA1 were risk factors. The molecular formula of the model was as follows: RS = (0.1917 * SHCBP1 + 0.0836 * TMEM45A–0.0721 * PIGR–0.1633 * IGLV6–57–0.0489 * TCN1-0.0676 * GFRA1). In the bootstrap set, AUC at one year was 0.767, AUC at three years was 0.717, and AUC at five years was 0.694 (Figure 6D, Supplementary Table 5). Besides, we used the GEO database GSE20685 for external data verification. The AUC values of breast cancer patients predicted by our model were 0.743, 0.650, and 0.615 at one, three, and five years, respectively (Figure 6E). We validated it separately in the GEO dataset: AUC at one year was 0.762, AUC at three years was 0.734, AUC at five years was 0.758, AUC at one year was 0.771, AUC at three years was 0.712, and AUC at five years was 0.652 in TCGA. In the combined dataset of TCGA and GEO, the AUC at one year was 0.766. The AUC at three years was 0.716, and the AUC at five years was 0.686 (Figures 6F–H). BC patients were randomly selected for scoring, total point = point (T) + point (N) + point (RS) + point (stage) + point (age) below by combining RS and clinicopathological features, using the nomogram (a quantitative method), as shown in Figure 6I. The total score corresponds to the scale in the figure. It can predict 1-year, 3-year, and 5-year OS of BC patients. The calibration curve showed adequate consistency between the predicted values of the 1-year, 3-year, and 5-year OS nomograms and the actual observed values (Figure 6J).

FIGURE 6

Figure 6 (A, B) LASSO variable trajectory plot for 1,000 cross validations (A) and LASSO coefficient profile (B). (C) Forest Plot for Multifactorial Cox Regression Analysis. (D) ROC curves and AUCs for 1-, 3-, and 5-year survival rates. (E) ROC curves and AUCs for 1-, 3-, and 5-year survival rates in another independent GEO dataset. (F) ROC curves and AUCs for 1-, 3-, and 5-year survival rates in the TCGA and GEO merged datasets. (G) ROC curves and AUCs for 1-, 3-, and 5-year survival rates in TCGA datasets. (H) ROC curves and AUCs for 1-, 3-, and 5-year survival rates in GEO datasets. (I) The nomogram used to calculate the survival rates of 1-, 3-, and 5-years for patients with BC. (J) Calibration curve for nomograms.

3.6 Dividing the high- and low-risk groups and observing their distribution in clusters

The patients in the TCGA, GEO, and total datasets were divided into high-risk and low-risk groups according to the median RS value. The differences between the groups were compared. Figure 7A shows that the genes (SLC7A11, SLC3A2, BRK1, ACTR2, ACTR3, RPN1, and NCKAP1) have the K-M survival analysis between the high- and low-risk groups showed that the survival differences between the high- and low-risk groups were statistically different (p < 0.001) in the TCGA dataset, GEO dataset, combined dataset and external validation of the GEO database. The prognosis of the significantly high-risk group was poor (Figure 7B). Seven of the 11 disulfide death genes showed significant differences in gene expression between high-risk and low-risk individuals (Figure 7C). Figure 7D illustrates the proportion of surviving and dying patients in two RS high and low groups, two DGs, and three DRG subtypes. Figures 7E, F show the RS distribution of two DG subtypes and three DRG subtypes, respectively. The PCA diagram shows that RS has an adequate grouping function (Figure 7G).

FIGURE 7

Figure 7 (A) Ranked dot, scatter plots and heat map of the model gene expressions in the TCGA and GEO merged datasets, TCGA datasets and GEO datasets. (B) Kaplan–Meier analyses of the OS between the TCGA and GEO merged datasets, TCGA datasets, GEO datasets and another independent GEO dataset. (C) RS score differences in eleven DGs. *p < 0.05, ***p < 0.001. (D) The subtype distributions among groups, risk scores and survival outcomes. (E) Variations in risk scores among DGs subtypes. (F) Variations in risk scores among CRGs subtypes. (G) Through PCA analysis, it can be seen that there are large transcriptome differences between high and low groups.

3.7 Different immune landscapes in the two risk groups

The CIBERSORT algorithm evaluated the relationship between RS and the relative number of immune cells. RS was positively correlated with the number of immune cells, such as Macrophages M0 and M2, Dendritic cells resting, Mast cells activated, B cells memory, NK cells resting, T cells CD4 memory resting, B cells naïve, Mast cells resting, and Monocytes. The expression of immune cells, such as dendritic cells activated, Plasma cells, Macrophages M1, and T cells CD8, was negatively correlated (p < 0.05). The stromal score, immune score, and ESTIMATE score of the high-risk group were higher. The difference between the groups was statistically significant (Figure 8A). Our study also examined the association between six genes and the number of immune cells (Figure 8B). According to our research, these six genes affect most immune cells. By observing the mutation frequency of each tumor sample and the mutation frequency of the gene, the TMB difference between the high- and low-risk groups was compared. From Figures 8C, D, it can be seen that the TMB of patients in the high-risk group is higher than that in the low-risk group, RS is positively correlated with TMB, and CRGsClusteA is significantly enriched in the low-risk area.

FIGURE 8

Figure 8 Comprehensive analysis of the risk scores in BC. (A) Correlations between immune cell types and risk score. (B) The six genes from the proposed model are correlated with the number of immune cells. (C) risk score and TMB spearman correlation analysis. (D) The differences in TMB between high- and low-risk groups. (E) Kaplan–Meier survival curves of BC patients between the H-TMB and L-TMB groups. (F) Kaplan–Meier survival curves of BC patients across H-TMB + high risk, H-TMB + low risk, L-TMB + high risk, and L-TMB + low risk. TMB, tumor mutational burden; H, high; L, low. (G, H) The somatic mutation features waterfall plot determined by high and low risk scores. One patient was represented by each column. The correct number represented each gene's frequency of mutation, and the upper barplot displayed TMB. The proportion of each variant type was displayed in the right barplot.

The TMB survival curve showed that patients with low TMB had a better prognosis (Figure 8E). Compared with other groups, BC patients with low risk and low TMB had the best prognoses (Figure 8F). The mutation rate of the low-risk group was 85.81% (Figure 8G). The mutation rate of the high-risk group was 85.78% (Figure 8H). The waterfall diagram shows that the mutation genes in the high- and low-risk groups are mainly PIK3 CA, TP53, TTN, CDH1, GATA3, MUC16, and MAP3K1, and the mutation rates of these genes are different in the high- and low-risk groups. PIK3CA had the highest mutation frequency in the low-risk samples, while TP53 had the highest in the high-risk groups. The mutation probability of TP53 in the high-risk group was 43%, while the mutation frequency of TP53 in the low-risk group was only 29%.

3.8 Drug sensitivity and different immunotherapy responses in the two risk groups

There was a positive correlation between RS and stem cells. The higher the RS, the higher the content of stem cells (p < 0.05) (Figure 9A). The TIDE score of the high-risk group was lower than that of the low-risk group, and the tumor immune escape ability was weak (Figure 9B). The ESTIMATE results showed that the stromal immunity and estimated scores of the high-risk group were low (Figure 9C). The drug treatment of breast cancer has broad research prospects and has attracted much attention. Therefore, the IC50 value of chemotherapeutic drugs for BC was calculated, and the relationship between RS and drug resistance was analyzed. We noted that in addition to docetaxel (microtubule depolymerization inhibitors) and parthenolide (NF-κB inhibitors) in high-risk patients with lower IC50. In contrast, other drugs [ABT.888 (Veliparib, PARP inhibitors)], AG.014699 (Rucaparib, PARP inhibitors), AMG.706 (Motesanib, VEGFR inhibitors), ATRA, AUY922 (Luminespib, HSP90 inhibitors), GDC0941 (Pictilisib, PI3K inhibitors), Metformin, Methotrexate, Nilotinib, Nutlin.3a (MDM2 inhibitors), Roscovitine, Temsirolimus, and Tipifarnib) had lower IC50 in low-risk patients (Figure 9D). We performed immune checkpoint analysis on high-risk and low-risk groups to explore the precise use of immune checkpoint inhibitors (ICI) in breast cancer patients. Except for the high expression of the PVR gene in the high-risk group, the other immune checkpoint genes were increased in the low-risk group. The IPS of the low-risk group was significantly higher than that of the high-risk group, and the immunotherapy efficacy of the low-risk group was better (Figure 9E).

FIGURE 9

Figure 9 (A) Associates between the CSC index and the risk score. (B) TIDE score between high- and low-risk group. (C) The boxplots for the drug sensitivity analysis. (D) Differential expression analysis of the immune checkpoint genes between the high-risk and low-risk groups. IC50, the half-maximal inhibitory concentration; *, p < 0.05; **, p < 0.01; ***, p < 0.001. (E) The difference between anti-PD1 treatment and anti-CTLA-4 treatment between high and low risk groups.

3.9 Knockdown of TMEM45A and SHCBP1 inhibited BC cell proliferation

We conducted a search on the DepMap Portal to assess the expression levels of TMEM45A and SHCBP1 in various breast cancer cell lines. Ultimately, we observed that these two genes exhibited significantly high expression in the MDA-MB-468 cell line (Supplementary Figure 1). In our study, we constructed a predictive model using six dual sulfur death-related genes. We selected the high-risk genes TMEM45A and SHCBP1 in our model to validate the potential for targeted gene therapy. TCGA and GEO data analyses revealed that TMEM45A and SHCBP1 were highly expressed in breast cancer. Therefore, we knocked down TMEM45A and SHCBP1 and further investigated the role of TMEM45A and SHCBP1 in MDA-MB-468 breast cancer cells in vitro.

First, we established MDA-MB-468 shTMEM45A and shSHCBP1 cell lines through lentiviral transduction. As shown in Figure 10A, the knockdown of TMEM45A and SHCBP1 in MDA-MB-468 was satisfactory. We further investigated the impact of TMEM45A and SHCBP1 knockdown on the functionality of breast cancer cells through CCK-8 and colony formation assays to determine their effects on BRCA cell proliferation.

FIGURE 10

Figure 10 (A) The efficiency of silencing TMEM45A and SHCBP1 was indicated by RT-qPCR in MDA-MB-468 cell lines. (B) The MTT assay revealed that, in comparison to the control group, the proliferation capability of the MDA-MB-468 cell line with TMEM45A and SHCBP1 knockdown significantly diminished. (C, D) The clonogenic assay revealed that the depletion of TMEM45A and SHCBP1 attenuated the proliferative capacity of MDA-MB-468 cells. The data is presented as the mean from at least three independent experiments. (*p<0.05; **p<0.01; ***p< 0.001; ****p< 0.0001).

As expected, the CCK-8 assay (Figure 10B) showed that the knockdown of TMEM45A and SHCBP1 in MDAMB468 similarly inhibited the proliferation of breast cancer cells. The colony formation assays (Figure 10C, D) demonstrated that the knockdown of TMEM45A and SHCBP1 significantly and independently inhibited cell proliferation (p < 0.05). These results collectively indicate that TMEM45A and SHCBP1 influence the proliferation of BC cells.

3.10 Sequencing and functional enrichment analysis

The knockdown of SHCBP1 and TMEM45A was achieved by generating shSHCBP1_1, shSHCBP1_2, shTMEM45A_1, and shTMEM45A_2. Differential analysis was conducted using the average counts from three replicate sequencing experiments (Supplementary Table 6). shSHCBP1_1 exhibited differential expression in 1171 genes (687 upregulated, 484 downregulated) compared to the NC group, while shSHCBP1_2 showed differential expression in 593 genes (421 upregulated, 172 downregulated) compared to the NC group (Figures 11A–C). shTMEM45A_1 exhibited differential expression in 7710 genes (3896 upregulated, 3814 downregulated) compared to the NC group, while shSHCBP1_2 showed differential expression in 7272 genes (3683 upregulated, 3589 downregulated) compared to the NC group (Figures 11F–H). We identified the intersecting genes in the upregulated and downregulated gene sets from the two independent replicate experiments and subjected these genes to GO and KEGG enrichment analyses (Figures 11D, E, I, J).

FIGURE 11

Figure 11 The volcano plot displays the differentially expressed genes between shSHCBP1_1 (A) and shSHCBP1_2 (B) compared to the NC group. (C) The Venn diagram shows the overlap in differentially expressed genes between shSHCBP1_1 and shSHCBP1_2 compared to the NC group. (D, E) GO and KEGG Enrichment Analysis of Overlapping Differentially Expressed Genes in shSHCBP1_1 and shSHCBP1_2 vs. NC. The volcano plot displays the differentially expressed genes between shTMEM45A_1 (F) and shTMEM45A_2 (G) compared to the NC group. (H) The Venn diagram shows the overlap in differentially expressed genes between shTMEM45A_1 and shTMEM45A_2 compared to the NC group. (I, J) GO and KEGG Enrichment Analysis of Overlapping Differentially Expressed Genes in shTMEM45A_1 and shTMEM45A_2 vs. NC.

4 Discussion

Regulated cell death (RCD) is a type of cell death controlled by specific molecular pathways and regulated by genetic or pharmacological manipulation (20). Recently, disulfidptosis has been defined as a new RCD. Previous studies have suggested that SLC7A11-mediated cystine intake is critical in promoting glutathione biosynthesis and inhibiting oxidative stress and ferroptosis. Subsequently, SLC7A11 was found to significantly promote cell death under glucose starvation (21–23). Subsequently, it was found that SLC7A11-mediated reduction of cystine to cysteine was highly dependent on the reduced nicotinamide adenine dinucleotide phosphate (NADPH) produced by the glucose–pentose phosphate pathway (24). Recently, Liu et al. proposed that disulfidptosis is an abnormal accumulation of disulfides, such as cystine, which induces disulfide stress, causes disulfide bond cross-linking and cytoskeleton contraction of the actin cytoskeleton, and ultimately induces cell death6. Not only does RCD play a key role in body development and cell homeostasis, but its dysregulation is also causally linked to many diseases, including cancer. Escape from cell death is considered to be one of the core markers of cancer. The importance of other RCDs in BC has been revealed, but the role of disulfidptosis in BC remains unclear (25–29). Our study explored the importance of disulfidptosis in predicting the prognosis, survival time, immunotherapy response, and chemosensitivity of BC patients. This result may lay the foundation for precisely treating BC breast protrusion-related diseases.

We first performed a differential gene expression analysis, gene copy number mutation between tumor and normal tissues, and tumor mutation load analysis on 11 DGs. We found that, except for NCKAP1, the remaining DGs had significant expression differences between the two groups. Most of them had increased gene expression in breast cancer, meaning these genes may be involved in some BC generation and development processes. In addition to the decrease in CYFIP1 and WASF2 gene copy numbers in tumor samples, the copy number of other DGs increased to varying degrees, consistent with the decline in CYFIP1 and WASF2 gene expression in BC patients. It has been reported that ABI2 can promote the growth and metastasis of HCC. In BC patients, the ABI2 gene copy number increased. Still, gene expression decreased, indicating that there may be other mechanisms in the body, such as methylation, acetylation, ubiquitination, and so on, to regulate ABI2 expression. Among the 911 breast cancer mutation data samples, only 33 had DG mutations. The highest mutation probability was 1% for NCKAP1 and CYFIP1. DGs showed more specific genetic stability than the 53% mutation rate of high-mutation genes, such as TP53 and PIK3CA (30). K-M survival analysis showed that the expression levels of SLC7A11, SLC3A2, RPN1, NCKAP1, BRK1, ACTR2, ACTR3, and RAC1 could independently predict the prognosis of patients (p < 0.05). Among these 11 DGs, SLC7A11 can affect non-small cell lung cancer (31), renal cell carcinoma (32), prostate cancer progression through ferroptosis, and WASF2 is associated with poor ovarian cancer prognosis (33). These results further indicate that DGs play an important role in the development and progression of tumors.

Consensus clustering is a standard unsupervised clustering method for cancer subtype classification research. It can distinguish samples into several subtypes according to different omics datasets to find new disease subtypes or compare different subtypes (34–36). We consistently clustered 11 prognostic DGs in the TCGA and GEO breast cancer data and identified two DG clusters. The OS of BC patients was statistically different between the two clusters, and the OS of group A was longer than that of group B. GSVA analysis showed that some pathways were enriched differently between the two DG groups, such as arachidonic acid metabolism, ubiquitin-mediated proteolysis, and taurine and hypotaurine metabolism. The ssGSEA analysis showed that most of the immune cells in group A were widely enriched and infiltrated, which may inhibit tumor cells through an immune response. In contrast, the lack of immune cells and immunosuppression in group B may be related to the poor prognosis of patients. PCA analysis showed a considerable difference between groups A and B.

A total of 225 DRGs were obtained, and GO and KEGG enrichment analyses were performed on these genes. GO analysis showed that these DRGs were mainly enriched in microtubule binding, tubulin binding, and other biological behaviors closely related to the construction of the cytoskeleton. Therefore, it can be speculated that these DRGs may mediate changes in the cytoskeleton and cause cell death, which is consistent with previous studies 6. Furthermore, through KEGG enrichment analysis, these differential genes were mainly enriched in signal transduction pathways, such as the cell cycle, PPAR signaling pathway, ECM-receptor interaction, and other pathways, indicating that intercellular interaction may be a critical link in disulfidptosis-induced cell death.

One-hundred-and-fourteen of the 225 differential genes were associated with BC prognosis. According to the expression of these 114 DRGs, breast cancer was divided into three subtypes by consensus clustering. The prognosis of group A was the best, followed by group B, and the prognosis of group C was relatively poor. In this study, we further developed six DRGs (SHCBP, TMEM45A, PIGR, IGLV6-57, TCN1, GFRA1) to construct RS to predict the prognosis of BC patients and used TCGA, GEO separated and TCGA and GEO joint databases to evaluate the prognostic value of RS through a survival curve, RS map, survival state map, and heatmap. RS = (0.1917 * SHCBP1 + 0.0836 * TMEM45A-0.0721 * PIGR-0.1633 * IGLV6-57-0.0489 * TCN1-0.0676 * GFRA1). SHCBP1 has been previously reported as an immune-related biomarker for cancer diagnosis and prognosis and a potential therapeutic target for tumor immunotherapy (37). Jing et al. proposed that TMEM45A can be used as an oncogene in ovarian cancer and that inhibition of TMEM45A may be a therapeutic strategy for ovarian cancer. Wichitra et al. proposed that M1 macrophages can cause high expression of PIGR in breast cancer cells, and high expression of polymeric immunoglobulin receptor (PIGR) in breast cancer is associated with an increased 5-year survival rate (38), indicating that PIGR may be a protective factor for breast cancer. IGLV6-57 is also widely used in cancer diagnosis (39). TCN1 may play a carcinogenic role in colorectal cancer by regulating the ITGB4 signaling pathway leading to cytoskeleton damage and promoting cell death (40). This type of cell death may be disulfidptosis. Sunil Bhakta et al. proposed that GFRA1 is associated with targeted therapy for hormone receptor–positive breast cancer (41). The six genes constructing the model are inseparable from regulating tumor life activities.

The area under the ROC curve (AUC) was used to evaluate the predictive ability of RS for patient prognosis (42). A considerable AUC indicates that the model has good classification ability and can compare features (43). The AUC of our model was 0.762, 0.734, and 0.758 at one year, three years, and five years, respectively, which was significantly higher than most prediction models. We integrated all the essential information of a series of prognostic models, including author, year, and genetic characteristics, and constructed AUC to verify the diagnostic performance of the model. After comparison, we found that our model had the best diagnostic performance (Table 1). Table 1 shows a metabolic-related 4-gene prognostic model with AUCs of 0.764, 0.689, and 0.612 for 1-year, 3-year, and 5-year, respectively (44). Another pyroptosis model consisted of 16 genes with AUC of 0.756, 0.752, and 0.723 for 1, 3, and 5 years (45). Its diagnostic value is almost the same as our model, but it is composed of 16 genes, while our model has only six genes that can be better applied in clinical practice. In two cuproptosis-related prognostic models, 1-year, 3-year, and 5-year AUC values were 0.685, 0.678, and 0.678 (46), and for another 1-, 3-, and 5-year model, the AUCs were 0.554, 0.527, and 0.649 (47). The AUC of the ferroptosis model constructed by Zhu et al. was 0.821, 0.678, and 0.651 in 1-, 3-, and 5-years (48). The AUCs of the ferroptosis model produced by Zhu et al. were 0.821, 0.678, and 0.651 in 1-, 3-, and 5-year models. Although it had a good prediction effect in the first year, our model AUC was above 0.734, which was more stable. We also list other models more concerned with short-term survival rates, such as the AUC of the necroptosis-related model at 1-, 2-, and 3-year patients, which were 0.701, 0.716, and 0.708 (49). The AUC of the necroptosis-related prognostic model constructed by Chen et al. in 1-, 3-, and 5-year patients were 0.731, 0.643, and 0.641, respectively (50). Our prognostic model has adequate predictive value. Our model involves only six genes, while other models tend to have more genes. To a certain extent, our RS model is more convenient to use. The C-index of the nomogram was 0.762 (95% CI: 0.711–0.813), indicating that the predicted results of 1-, 3-, and 5-year patients were consistent with the actual results. The expression of DGs in the high-risk group was significantly increased, suggesting that genes have an adverse effect on the prognoses of patients. Therefore, the prognostic model constructed by DRGs is reliable and accurate. On the other hand, it also shows the critical role of DGs in the occurrence, development, and life processes of tumors. In summary, the model we constructed can be considered a suitable prognostic signal, and its mechanism of action in BC deserves further exploration and verification.

TABLE 1

Table 1 The area under the ROC curve (AUC) showed the sensitivity and specificity of the known gene signatures in predicting the prognosis of BC patients.

Subsequently, we found the immune factors that determined the prognoses of the high- and low-risk groups through the immune analysis of the patients. The high-risk group had a large number of macrophage M2 cell infiltrations. Macrophages M2 increase is associated with tumor growth and poor prognosis of cancer 48, and are considered essential biomarkers in cancer diagnosis and a potential target for cancer treatment (51, 52). T cells CD8 is a significant member of the low-risk group. In previous reports, T cell CD8 can prevent tumor growth and promote immune response and immunotherapy (53). This may be one of the reasons for the survival difference between high- and low-risk groups. The TMB study showed that PIK3CA had the highest mutation frequency in the low-risk samples, while TP53 had the highest mutation frequency in the high-risk groups. The mutation rate of TP53 in patients with an increased risk of BC is 43%, while the mutation rate of TP53 in patients with low risk is only 29%. Previous studies have reported that the mutation of the TP53 gene is related to the poor therapeutic effect and prognosis of BC (54, 55), confirming the poor prognosis of patients with a high TP53 mutation rate in our high-risk group.

Tumor treatment is an essential field of concern. Through IC 50 screening analysis of potential chemical drugs, we realized that low-risk patients might be more sensitive to chemotherapy, ABT.888 (Veliparib, PARP inhibitors), AG.014699 (Rucaparib, PARP inhibitors), AMG.706 (Motesanib, VEGFR inhibitors), ATRA, and AUY922 (Luminespib, HSP90 inhibitors). GDC0941 (Pictilisib, PI3K inhibitors), Metformin, Methotrexate, Nilotinib, Nutlin.3a (MDM2 inhibitors), Roscovitine, Temsirolimus, Tipifarnib, and other drugs had a lower IC50 in the low-risk group. High-risk patients resisted PARP and VEGFR inhibitors and other drugs and were sensitive to docetaxel (a microtubule depolymerization inhibitor) and parthenolide (NF-κB inhibitors). Immunotherapy also occupies a crucial position in the treatment of BC patients. The ESTIMATE results showed that the stromal immunity and estimated scores of the high-risk group were low. The TIDE analysis showed that the TIDE immune escape score of the high-risk group was low. The combination of PD-1 and PDL-1 caused T cells to lose the ability to attack cancer cells, resulting in the immune escape of tumor cells. The expression of immune checkpoint genes, such as PD-1 and PDL-1, in high-risk patients, was low, which may be the reason for the low TIDE immune escape score in the high-risk group. However, based on the difference in immune checkpoints between the high- and low-risk groups, distinguished by the RS model constructed by DRGs, we found that the expression of immune checkpoint genes in the high-risk group was significantly lower than that in the low-risk group. It is speculated that patients in the low-risk group are more likely to induce an antitumor immune response through immune genes, thus benefiting from immunotherapy, which is consistent with the high Stromal Score, Immune Score, and ESTIMATE Score results of the low-risk population in the previous article. TMB is a reliable biomarker for predicting treatment outcomes in cancer patients treated with ICI (56, 57). This is consistent with our earlier results that ‘RS is significantly associated with TMB. Patients in the high-risk group had higher TMB, antitumor immune dysfunction, poor ICI treatment, and poor prognosis. In summary, combined with drug sensitivity and immune efficacy analysis, we predict that low-risk patients will benefit more from combining chemotherapy and immunotherapy, providing a basis for individualized treatment of BC patients. More drugs are worthy of selection and development for low-risk patients. The effects of chemotherapy and immunotherapy in the high-risk group were poor. The high expression of disulfidptosis genes (e.g., SLC3A2, RPN1, BRK1, ACTR2, ACTR3, SLC7A11, and NCKAP1) in the KM analysis of breast cancer indicated the poor prognosis of patients. The difference analysis between the high- and low-risk groups differed. These genes were highly expressed in the high-risk group. Through our study, targeted therapy of the disulfidptosis gene will benefit the high-risk group.

In addition, we conducted in vitro validation through cell formation assays and CCK-8 on the high-risk genes (TMEM45A and SHCBP1) within the RS model. This initial validation provided preliminary evidence of the feasibility of utilizing these two genes for targeted therapy in breast cancer (BC). Subsequently, we conducted cellular sequencing following the gene knockout of TMEM45A and SHCBP1 to infer their roles in the biological mechanisms.

This study systematically analyzed the role of disulfidptosis-related genes in the prognosis of breast cancer and their correlation with the tumor microenvironment and clinical characteristics and constructed a better prognosis prediction model. The model performs well in predicting the survival outcome, tumor immune microenvironment, and immunotherapy response of BC patients. Furthermore, it expounds on its correlation with TMB, immunity, and clinical treatment (chemotherapy and immunotherapy), providing a reference for individualized and precise breast cancer treatment. Finally, Experimental validation and sequencing were also conducted to substantiate these findings.

5 Conclusion

This study has the following contributions. First, this study is the first to identify the subtypes associated with disulfidptosis and to create a predictive model based on breast cancer DRGs. Although disulfidptosis differs from other recognized cell death methods, it may provide new therapeutic possibilities for cancer treatment. Second, a variety of different technologies and databases are used to improve the reliability of the results. We also defined subtypes associated with disulfidptosis and created a predictive model for the screening and testing processes. Finally, we first proposed a disulfidptosis gene-targeted therapy for high-risk BC groups. We intend to gather further patient samples during subsequent clinical investigations to substantiate the dependability of our model.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/Supplementary Materials. The sequencing raw data is stored in the GEO database under GSE246486.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

JHL, ZR, and MR designed and conceived this project; XW and JY completed these experiments; JS, JHL, and XW interpreted and analyzed the data; JHL and JCL revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was supported by Wu Jieping Medical Foundation (No. 320.6750.2020-7-1), Anhui Provincial Natural Science Foundation (No. 2308085MH276), Anhui Medical University Clinical Science Foundation (No. 2021xkj140). Research Fund of Anhui Institute of translational medicine (No. 2022zhyx-C29).

Acknowledgments

We acknowledge TCGA and GEO database for providing their platforms and contributors for uploading their meaningful datasets.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1198826/full#supplementary-material

References

1. Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. CA Cancer J Clin (2023) 73(1):17–48. doi: 10.3322/caac.21763

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Waks AG, Winer EP. Breast cancer treatment: A review. Jama (2019) 321(3):288–300. doi: 10.1001/jama.2018.19323

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Britt KL, Cuzick J, Phillips KA. Key steps for effective breast cancer prevention. Nat Rev Cancer (2020) 20(8):417–36. doi: 10.1038/s41568-020-0266-x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Perou CM, Sørlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, et al. Molecular portraits of human breast tumours. Nature (2000) 406(6797):747–52. doi: 10.1038/35021093

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lohr M, Hellwig B, Edlund K, Mattsson JS, Botling J, Schmidt M, et al. Identification of sample annotation errors in gene expression datasets. Arch Toxicol (2015) 89(12):2265–72. doi: 10.1007/s00204-015-1632-4

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Liu X, Nie L, Zhang Y, Yan Y, Wang C, Colic M, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol (2023) 25(3):404–14. doi: 10.1038/s41556-023-01091-2

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zhang C, Li Z, Qi F, Hu X, Luo J. Exploration of the relationships between tumor mutation burden with immune infiltrates in clear cell renal cell carcinoma. Ann Transl Med (2019) 7(22):648. doi: 10.21037/atm.2019.10.84

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

10. Jin Y, Wang Z, He D, Zhu Y, Chen X, Cao K. Identification of novel subtypes based on ssGSEA in immune-related prognostic signature for tongue squamous cell carcinoma. Cancer Med (2021) 10(23):8693–707. doi: 10.1002/cam4.4341

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Friedman JH, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw (2010) 33(1):122. doi: 10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shi X, Li Y, Pan S, Liu X, Ke Y, Guo W, et al. Identification and validation of an autophagy-related gene signature for predicting prognosis in patients with esophageal squamous cell carcinoma. Sci Rep (2022) 12(1):1960. doi: 10.1038/s41598-022-05922-4

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Liu XS, Pompey KT. Bootstrap estimate of bias for intraclass correlation. J Appl Meas (2020) 21(1):101–8.

PubMed Abstract | Google Scholar

14. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hu FF, Liu CJ, Liu LL, Zhang Q, Guo AY. Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Brief Bioinform (2021) 22(3). doi: 10.1093/bib/bbaa176

CrossRef Full Text | Google Scholar

16. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Tsherniak A, Vazquez F, Montgomery PG, Weir BA, Kryukov G, Cowley GS, et al. Defining a cancer dependency map. Cell (2017) 170(3):564–76. doi: 10.1016/j.cell.2017.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Grishagin IV. Automatic cell counting with ImageJ. Anal Biochem (2015) 473:63–5. doi: 10.1016/j.ab.2014.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Galluzzi L, Vitale I, Aaronson SA, Abrams JM, Adam D, Agostinis P, et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ (2018) 25(3):486–541. doi: 10.1038/s41418-017-0012-4

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Shin CS, Mishra P, Watrous JD, Carelli V, D'Aurelio M, Jain M, et al. The glutamate/cystine xCT antiporter antagonizes glutamine metabolism and reduces nutrient flexibility. Nat Commun (2017) 8:15074. doi: 10.1038/ncomms15074

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Koppula P, Zhang Y, Shi J, Li W, Gan B. The glutamate/cystine antiporter SLC7A11/xCT enhances cancer cell dependency on glucose by exporting glutamate. J Biol Chem (2017) 292(34):14240–9. doi: 10.1074/jbc.M117.798405

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Goji T, Takahara K, Negishi M, Katoh H. Cystine uptake through the cystine/glutamate antiporter xCT triggers glioblastoma cell death under glucose deprivation. J Biol Chem (2017) 292(48):19721–32. doi: 10.1074/jbc.M117.814392

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Liu X, Olszewski K, Zhang Y, Lim EW, Shi J, Zhang X, et al. Cystine transporter regulation of pentose phosphate pathway dependency and disulfide stress exposes a targetable metabolic vulnerability in cancer. Nat Cell Biol (2020) 22(4):476–86. doi: 10.1038/s41556-020-0496-x

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhu L, Tian Q, Jiang S, Gao H, Yu S, Zhou Y, et al. A novel ferroptosis-related gene signature for overall survival prediction in patients with breast cancer. Front Cell Dev Biol (2021) 9:670184. doi: 10.3389/fcell.2021.670184

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Hu T, Zhao X, Zhao Y, Cheng J, Xiong J, Lu C. Identification and verification of necroptosis-related gene signature and associated regulatory axis in breast cancer. Front Genet (2022) 13:842218. doi: 10.3389/fgene.2022.842218

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Xu L, Hu Y, Liu W. Pyroptosis-mediated molecular subtypes are characterized by distinct tumor microenvironment infiltration characteristics in breast cancer. J Inflammation Res (2022) 15:345–62. doi: 10.2147/jir.S349186

CrossRef Full Text | Google Scholar

29. Wang D, Wei G, Ma J, Cheng S, Jia L, Song X, et al. Identification of the prognostic value of ferroptosis-related gene signature in breast cancer patients. BMC Cancer (2021) 21(1):645. doi: 10.1186/s12885-021-08341-2

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lang GT, Jiang YZ, Shi JX, Yang F, Li XG, Pei YC, et al. Characterization of the genomic landscape and actionable mutations in Chinese breast cancers by clinical sequencing. Nat Commun (2020) 11(1):5679. doi: 10.1038/s41467-020-19342-3

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Lu X, Kang N, Ling X, Pan M, Du W, Gao S. MiR-27a-3p promotes non-small cell lung cancer through SLC7A11-mediated-ferroptosis. Front Oncol (2021) 11:759346. doi: 10.3389/fonc.2021.759346

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Li YZ, Zhu HC, Du Y, Zhao HC, Wang L. Silencing lncRNA SLC16A1-AS1 Induced Ferroptosis in Renal Cell Carcinoma Through miR-143-3p/SLC7A11 Signaling. Technol Cancer Res Treat (2022) 21:15330338221077803. doi: 10.1177/15330338221077803

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yang X, Ding Y, Sun L, Shi M, Zhang P, He A, et al. WASF2 serves as a potential biomarker and therapeutic target in ovarian cancer: A pan-cancer analysis. Front Oncol (2022) 12:840038. doi: 10.3389/fonc.2022.840038

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Shao W, Yang Z, Fu Y, Zheng L, Liu F, Chai L, et al. The pyroptosis-related signature predicts prognosis and indicates immune microenvironment infiltration in gastric cancer. Front Cell Dev Biol (2021) 9:676485. doi: 10.3389/fcell.2021.676485

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Chen D, Huang H, Zang L, Gao W, Zhu H, Yu X. Development and verification of the hypoxia- and immune-associated prognostic signature for pancreatic ductal adenocarcinoma. Front Immunol (2021) 12:728062. doi: 10.3389/fimmu.2021.728062

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Coleman S, Kirk PDW, Wallace C. Consensus clustering for Bayesian mixture models. BMC Bioinf (2022) 23(1):290. doi: 10.1186/s12859-022-04830-8

CrossRef Full Text | Google Scholar

37. Wang N, Zhu L, Wang L, Shen Z, Huang X. Identification of SHCBP1 as a potential biomarker involving diagnosis, prognosis, and tumor immune microenvironment across multiple cancers. Comput Struct Biotechnol J (2022) 20:3106–19. doi: 10.1016/j.csbj.2022.06.039

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Asanprakit W, Lobo DN, Eremin O, Bennett AJ. M1 macrophages evoke an increase in polymeric immunoglobulin receptor (PIGR) expression in MDA-MB468 breast cancer cells through secretion of interleukin-1β. Sci Rep (2022) 12(1):16842. doi: 10.1038/s41598-022-20811-6

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Miao Y, Wang J, Ma X, Yang Y, Mi D. Identification prognosis-associated immune genes in colon adenocarcinoma. Biosci Rep (2020) 40(11). doi: 10.1042/bsr20201734

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhu X, Jiang X, Zhang Q, Huang H, Shi X, Hou D, et al. TCN1 deficiency inhibits the Malignancy of colorectal cancer cells by regulating the ITGB4 pathway. Gut Liver (2022) 17(3):412–29. doi: 10.5009/gnl210494

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Bhakta S, Crocker LM, Chen Y, Hazen M, Schutten MM, Li D, et al. An anti-GDNF family receptor alpha 1 (GFRA1) antibody-drug conjugate for the treatment of hormone receptor-positive breast cancer. Mol Cancer Ther (2018) 17(3):638–49. doi: 10.1158/1535-7163.Mct-17-0813

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zhang D, Li Y, Yang S, Wang M, Yao J, Zheng Y, et al. Identification of a glycolysis-related gene signature for survival prediction of ovarian cancer patients. Cancer Med (2021) 10(22):8222–37. doi: 10.1002/cam4.4317

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Lv W, Yu H, Han M, Tan Y, Wu M, Zhang J, et al. Analysis of tumor glycosylation characteristics and implications for immune checkpoint inhibitor's efficacy for breast cancer. Front Immunol (2022) 13:830158. doi: 10.3389/fimmu.2022.830158

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Lu J, Liu P, Zhang R. A metabolic gene signature to predict breast cancer prognosis. Front Mol Biosci (2022) 9:900433. doi: 10.3389/fmolb.2022.900433

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Chen H, Luo H, Wang J, Li J, Jiang Y. Identification of a pyroptosis-related prognostic signature in breast cancer. BMC Cancer (2022) 22(1):429. doi: 10.1186/s12885-022-09526-z

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Li Z, Zhang H, Wang X, Wang Q, Xue J, Shi Y, et al. Identification of cuproptosis-related subtypes, characterization of tumor microenvironment infiltration, and development of a prognosis model in breast cancer. Front Immunol (2022) 13:996836. doi: 10.3389/fimmu.2022.996836

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sha S, Si L, Wu X, Chen Y, Xiong H, Xu Y, et al. Prognostic analysis of cuproptosis-related gene in triple-negative breast cancer. Front Immunol (2022) 13:922780. doi: 10.3389/fimmu.2022.922780

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhu J, Chen Q, Gu K, Meng Y, Ji S, Zhao Y, et al. Correlation between ferroptosis-related gene signature and immune landscape, prognosis in breast cancer. J Immunol Res (2022) 2022:6871518. doi: 10.1155/2022/6871518

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Yu H, Lv W, Tan Y, He X, Wu Y, Wu M, et al. Immunotherapy landscape analyses of necroptosis characteristics for breast cancer patients. J Transl Med (2022) 20(1):328. doi: 10.1186/s12967-022-03535-z

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Chen F, Yang J, Fang M, Wu Y, Su D, Sheng Y. Necroptosis-related lncRNA to establish novel prognostic signature and predict the immunotherapy response in breast cancer. J Clin Lab Anal (2022) 36(4):e24302. doi: 10.1002/jcla.24302

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Cha YJ, Koo JS. Role of tumor-associated myeloid cells in breast cancer. Cells (2020) 9(8). doi: 10.3390/cells9081785

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Komohara Y, Jinushi M, Takeya M. Clinical significance of macrophage heterogeneity in human Malignant tumors. Cancer Sci (2014) 105(1):1–8. doi: 10.1111/cas.12314

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Farhood B, Najafi M, Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunotherapy: A review. J Cell Physiol (2019) 234(6):8509–21. doi: 10.1002/jcp.27782

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Olivier M, Hollstein M, Hainaut P. TP53 mutations in human cancers: origins, consequences, and clinical use. Cold Spring Harb Perspect Biol (2010) 2(1):a001008. doi: 10.1101/cshperspect.a001008

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Takahashi S, Fukui T, Nomizu T, Kakugawa Y, Fujishima F, Ishida T, et al. TP53 signature diagnostic system using multiplex reverse transcription-polymerase chain reaction system enables prediction of prognosis of breast cancer patients. Breast Cancer (2021) 28(6):1225–34. doi: 10.1007/s12282-021-01250-z

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Addeo A, Friedlaender A, Banna GL, Weiss GJ. TMB or not TMB as a biomarker: That is the question. Crit Rev Oncol Hematol (2021) 163:103374. doi: 10.1016/j.critrevonc.2021.103374

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Cao D, Xu H, Xu X, Guo T, Ge W. High tumor mutation burden predicts better efficacy of immunotherapy: a pooled analysis of 103078 cancer patients. Oncoimmunology (2019) 8(9):e1629258. doi: 10.1080/2162402x.2019.1629258

PubMed Abstract | CrossRef Full Text | Google Scholar