Integrated machine learning and bioinformatic analyses constructed a sulfur metabolism-related breast cancer risk model and identified heat-shock protein A9 as a potential therapeutic target for human breast cancer.

Introduction
Breast cancer is one of the leading causes of mortality among women globally [1], presently surpassing lung cancer; the disease became the most commonly diagnosed malignancy in 2020 [2]. The disparities in breast cancer mortality rates are clear [3–5]. Although mortality rates have declined in North America and the European Union, breast cancer remains a predominant cause of cancer-related mortality in a few developed nations [6, 7]. Breast cancer is characterized by its heterogeneity, including molecular diversity resulting from genetic alterations that affect cell metabolism, growth, and tumorigenesis [8, 9]. Therefore, it is crucial to develop a new model that can distinguish tumor subtypes based on comparable molecular characteristics to improve clinical treatments of breast cancer [10].
Sulfur is an essential non-metallic element required by living organisms [11]. Iron-sulfur clusters formed by iron and sulfur fundamentally influence electron transfer and catalysis [12]. Furthermore, sulfur influences the synthesis of certain amino acids (cysteine and methionine) in the human body. Metabolic reprogramming is a hallmark of cancer cells that is often highly dependent on specific substrates or metabolic pathways [13, 14]. As metabolic reprogramming has become increasingly recognized, the effects of some essential metabolic pathways on cancer have attracted attention [15]. Ferroptosis is a form of programmed cell death induced by disruptions in iron metabolism associated with malignant features, immunosuppression, and therapeutic responses in various tumors [16]. Dysregulated sulfur metabolism has been implicated in tumorigenesis [17]. Recent studies have identified a novel mechanism of controlled cell death based on disulfide precipitation, which involves nicotinamide adenine dinucleotide phosphate depletion and proline accumulation within the cell [18]. Foods rich in organosulfur compounds, including garlic and onion, exert anti-tumor effects due to the generation of reactive sulfur species [19].
Increasing evidence has revealed the therapeutic potential of modulating the sulfur metabolic pathway [20, 21]. Although previous studies have focused on excess sulfur, emerging evidence suggests that genes regulating sulfur metabolism are also involved in tumorigenesis and progression [17, 22]. Currently, there are no breast cancer prediction models targeting sulfur metabolism that can help physicians predict patient outcomes. Consequently, it is essential to develop an innovative prognosis framework based on the molecular features associated with sulfur metabolism among patients with breast cancer. This study aimed to develop a clinical prognostic model based on genetic and clinical data from the Cancer Genome Atlas (TCGA) for cases of breast invasive carcinoma (BRCA). Furthermore, we aim to identify the most pertinent risk gene for breast cancer and validate its functionality within MDA-MB-231 and MDA-MB-468 breast cancer cells.

Results

Identification of differentially expressed sulfur metabolism-related genes in BRCA
Figure 1 depicts the comprehensive workflow of our study. We examined data from TCGA, which included 1099 cases of BRCA, and GTEx, which included 292 normal tissues. After rigorous analysis, we identified 153 distinct sulfur metabolism-related genes. These were selected based on stringent criteria (|log2 (fold change)| > 2.0 and adjusted p < 0.05) from a pool of 332 sulfur metabolism-associated genes in comparisons between tumor tissues and normal BRCA patients (Fig. 2A). Supplementary Tables 1 and 2 present detailed information about these genes. Figure 2B depicts the possible interaction of these differentially expressed genes (DEGs). Subsequently, we utilized univariate Cox regression analysis to identify 12 sulfur metabolism-related genes with a significance level of p < 0.05. Among these, six were potential risky genes (CHPF, ENOPH1, HLCS, HSPA9, SLC38A7, and SLC7A5), whereas the other six were deemed potential protective genes (ACOT2, ACOT4, ELOVL2, MICAL1, SPOCK2, and TCF7L2) (Fig. 2C). Furthermore, a multivariate Cox regression analysis revealed two independent risk genes (CHPF, and HSPA9) and three independent protective genes (ACOT4, ELOVL2, and SPOCK2). As illustrated in Fig. 2D and E, boxplots and heat maps indicate that ACOT4, CHPF, ELOVL2, ENOPH1, HLCS, HSPA9, SLC38A7, SLC7A5, and SPOCK2 were upregulated in tumor tissues, whereas ACOT2 and TCF7L2 were downregulated.

Development of a prognostic model for breast cancer utilizing genes associated with sulfur metabolism
We identified specific genes strongly associated with BRCA (refer to Fig. 3A, B). Subsequently, we developed a prognostic model based on sulfur metabolism-related genes. This model included ACOT2, ACOT4, CHPF, ELOVL2, HLCS, HSPA9, MICAL1, SPOCK2, and TCF7L2. We carefully eliminated weak predictors whose coefficients gravitated toward zero during the least absolute shrinkage and selection operator (LASSO) regression process. The prognostic index for BRCA samples was determined using the formula: Risk score = (ACOT2 expression) × (-0.192104468) + (ACOT4 expression) × (-0.527345667) + (CHPF expression) × (0.517461254) + (ELOVL2 expression) × (-0.488079095) + (HLCS expression) × (0.262180346) + (HSPA9 expression) × (0.452787841) + (MICAL1 expression) × (-0.177572806) + (SPOCK2 expression) × (-0.487639255) + (TCF7L2 expression) × (-0.209252575). To validate the sulfur metabolism-related model’s prognostic efficacy for BRCA patients, the patients were stratified into high- and low-risk groups based on the median risk score. High scores were associated with worse prognoses for BRCA patients. Notably, in the high-risk group, HLCS, HSPA9, and CHPF exhibited higher expression, and ELOVL2, ACOT2, ACOT4, TCF7L2, MICAL1, and SPOCK2 exhibited lower expression levels (Figs. 3C, D). Kaplan–Meier curves further validated this observation, demonstrating reduced survival durations for patients in the high-risk category (Fig. 3E). Furthermore, time-dependent receiver operating characteristic (ROC) curves highlighted the risk score’s prognostic precision for overall survival (OS) at 1, 2, and 3 years, yielding area under the curve (AUC) values of 0.736, 0.755, and 0.711, respectively (Fig. 3F). These findings highlight the strong predictive potential of the sulfur metabolism-associated gene signature in our model for clinical outcomes in BRCA.

Construction of nomogram and calibration curves
To enhance prognostic assessment of BRCA patients, we developed a nomogram for OS that included the risk score and additional prognostic indicators, including T, N, M stages, and age [23]. Figure 4A depicts that the risk score emerged as a pivotal determinant amid various clinical parameters. To evaluate the nomogram’s effectiveness, we examined the calibration plot, which indicated a strong correlation between our nomogram and the actual BRCA patients’ survival (Fig. 4B). Subsequently, ROC analysis revealed higher accuracy (AUC = 0.659) compared to conventional prognostic factors (Fig. 4C). These results demonstrate that the nomogram incorporating our risk index system is effective in accurately estimating the OS of BRCA patients.

Potential therapeutic response in the different risk groups
The R software and the pRRophetic package were utilized to predict chemotherapeutic responses based on the genomics of drug sensitivity in the cancer database for BRCA patients. Five small-molecule compounds with significantly different responses between the high- and low-risk groups were identified. The results are: NU7441 (p = 0.0082), cyclopamine (p = 0.0015), embelin (p = 0.0042), A-443,654 (p = 0.013), and docetaxel (p = 0.047) (Supplementary Fig. 1 A, C, E, G, and I); 3D data for these five compounds was obtained from PubChem (Supplementary Fig. 1B, D, F, H, and J). Our data indicate that patients in the high-risk group exhibited increased sensitivity to all five small molecular compounds, suggesting that some specific chemotherapies may confer therapeutic benefits for those with BRCA mutations, warranting additional investigation into this matter. These results indicate novel potential candidate compounds for BRCA treatment.

Analyzing the enrichment of DEGs between the high and low-risk groups
We identified 3,255 DEGs, including 1,383 upregulated genes and 1,417 downregulated genes (|log2(FC)| > 2, p.adj < 0.05) (Fig. 5A), between high- and low-risk groups of patients. Supplementary Table 3 presents detailed information. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses of DEGs combined with log2(FC) values were also performed (Figs. 5B, C). The radar chart indicates that melanosomes, pigment granules, cadherin binding, protein folding, endoplasmic reticulum-protein-containing complex, unfolded protein binding, and ATP metabolic process were activated in the high-risk group, whereas terms including cytosolic ribosome, ribosomal subunit, and structural constituent of ribosomes were activated in the low-risk group. KEGG analysis revealed that protein processing in the endoplasmic reticulum, Parkinson’s disease, prostate cancer, and HIF-1 signaling pathways were significantly enriched in the high-risk group. Gene set enrichment analysis revealed that the specific GO terms, including unfolded protein binding, vesicle lumen, cadherin binding, endoplasmic reticulum-protein-containing complex, and protein stabilization, were positively enriched, whereas immune-related terms were negatively enriched in upregulated DEGs (Figs. 5D, E). These results revealed that the cell-protein folding system may be dysfunctional in high-risk patients.

Evaluation of the infiltration of immune cells within the tumor
The tumor microenvironment exhibits a complex interplay among malignant cells, immunological components, and the extracellular matrix. This dynamic environment significantly influences tumor progression, metastasis, and invasion. We conducted a thorough analysis of key factors, including tumor purity, stromal score, immune score, and ESTIMATE score, to gain deeper insights. The results, depicted in Fig. 6A, reveal that the high-risk group exhibited higher tumor purity, while the low-risk group exhibited significantly higher stromal, immune, and ESTIMATE scores compared to their high-risk counterparts (Figs. 6B–D). Furthermore, we investigated the infiltration levels for 28 distinct types of immune cells in high- and low-risk groups. Notably, almost all types of immune cells demonstrated increased infiltration levels in the low-risk group compared to the high-risk group. These include activated B cells, activated CD8 T cells, natural killer cells, central memory CD4 T cells, effector memory CD4 T cells, effector memory CD8 T cells, eosinophils, immature B cells, immature dendritic cells, macrophages, mast cells, myeloid-derived suppressor cells, memory B cells, monocytes, natural killer T cells, plasmacytoid dendritic cells, T follicular helper cells, and type 1 T helper cells (p < 0.05) (Fig. 6E). These findings imply that high-risk BRCA patients may exhibit lower immune cell infiltration levels, potentially indicating a higher resistance to immunotherapeutic interventions compared to low-risk patients.

Analysis of the association between sulfur metabolism-related genes and prognosis in BRCA
In our Kaplan–Meier survival plots, ACOT2, ACOT4, SPOCK2, ELOVL2, TCF7L2, HLCS, HSPA9, and CHPF were significantly correlated with OS in BRCA (p < 0.05) (Fig. 7A), and only ACOT2, ACOT4, ELOVL2, HSPA9, and CHPF expression exhibited significant association with patients’ disease-specific survival (DSS) (p < 0.05) (Fig. 7B). The risk genes, HSPA9 and CHPF, were associated with prognosis, regardless of OS or DSS.

HSPA9 knockdown inhibited BRCA cell malignant phenotypes
Conforming to our preceding univariate Cox regression assessment, HSPA9 emerged as the most significant risk gene associated with BRCA (Fig. 2C). To investigate the role of HSPA9 in BRCA, we initiated our analysis by assessing mRNA and protein expression levels. This was achieved using qRT-PCR and Western blotting techniques. HSPA9 was significantly upregulated compared to normal controls (Fig. 8A, B). Subsequently, we developed two specific siRNAs designed to inhibit HSPA9 expression in MDA-MB-231 and MDA-MB-468 cells. The efficacy of this knockdown approach was confirmed using Western blotting (Fig. 8C). A cell counting kit-8 assay was subsequently performed on MDA-MB-231 and MDA-MB-468 cells post-siRNA transfection. The findings highlighted a significant impairment in the proliferative ability of cancer cells in si-HSPA9 1 and 2 groups, compared to the si-NC group, at 24, 48, and 72 h post-transfection (Fig. 8D). Furthermore, the impact of HSPA9 knockdown on cell migration was evaluated through a wound healing assay. Figure 8E demonstrates that the migration rate of breast cancer cells was significantly inhibited at 24 and 48 h, a phenomenon concurrent with the silencing of HSPA9. In addition, we also tested the (epithelial-mesenchymal transition) EMT marker in HSPA9 knockdown breast cancer cells, and not surprisingly, the results showed us that epithelial cell marker E-cadherin was upregulation accompanied by decreasement of mesenchymal markers like N-cadherin and Vimentin after HSPA9 silence, which means EMT progress was reversed (Fig. 8F). The targeted inhibition of HSPA9 expression, accomplished using siRNA, markedly attenuated the malignant behavior exhibited by triple-negative breast cancer cells.

Discussion
Sulfur, a crucial element in the human body, influences energy metabolism and protein synthesis [24, 25]. Liu et al. recently identified a new form of cell death, termed disulfidptosis, which is associated with sulfur-based disulfide stress caused by the excessive intracellular accumulation of disulfides during glucose starvation conditions. They reported the latent effect of sulfur on cancer cell death [26, 27]. This novel cell death modality directly depends on sulfur-containing amino acids and results in cell death through the formation of aberrant disulfide bonds in the cytoskeleton, which is expected to draw increasing attention to sulfur metabolism [28–30].
Our risk model, developed using univariate Cox and LASSO Cox regression analyses, harnessed the potential of nine genes (TCF7L2, SPOCK2, MICAL1, HSPA9, HLCS, ELOVL2, CHPF, ACOT4, and ACOT2) derived from the TCGA BRCA cohort. We stratified patients into several groups based on the median risk score, as validated using Kaplan–Meier survival analysis and time-dependent ROC curve. This model emerges as a powerful prognostic tool for BRCA patients. By analyzing the detailed molecular alterations between these risk cohorts, we identified 1,417 downregulated genes and 1,383 upregulated genes in the high-risk group. This complex array of genetic alterations elucidates the underlying mechanisms. Furthermore, we identified five potential small-molecule compounds—docetaxel (p = 0.047), A-443,654 (p = 0.013), embelin (p = 0.0042), cyclopamine (p = 0.0015), and NU7441 (p = 0.0082)—demonstrating increased sensitivity in high-risk groups. These compounds exhibit potential for focused therapeutic approaches. Furthermore, our nomogram, which combines the risk score with clinical features (age and T, N, and M stages), enables clinicians to make accurate prognostic assessments for patient outcomes. This tool may aid clinicians in delivering more tailored medical interventions. Enrichment analysis revealed that the major functional differences between low- and high-risk groups included cellular biological processes, signal transduction pathways, protein synthesis, and immune-related terms. Immunotherapy represents a promising treatment for patients with breast cancer [31–33]. Future research may explore a combination of immunotherapy and other modalities as a novel treatment approach for breast cancer [33, 34], with the efficacy of immunotherapy dependent on the assessment of the patient’s immune infiltration level [35, 36]. The two different groups divided by our prognostic model demonstrated varying degrees of immune infiltration, indicating that our risk model may serve as a useful indicator of the patients’ immune status.
Univariate Cox analysis identified HSPA9 as the most relevant gene for patients with BRCA. HSPA9, or GRP-75 [37], is a chaperone protein found in the mitochondria, cell membrane, and endoplasmic reticulum [38, 39]. Previous studies have demonstrated that HSPA9 is crucial for tumor metastasis, epithelial-mesenchymal transition, and drug resistance [40–43]. Consequently, we assessed the role of HSPA9 in enhancing the malignant phenotype of breast cancer cells by modulating its expression. However, this study still has some shortcomings. The molecular mechanism by which HSPA9 promotes breast cancer cell proliferation and migration is unclear, and a specific difference in the immune response between the low- and high-risk groups requires further investigation. In follow-up investigation, the research point would be paied attention to validate the realiability and accuracy of the sulfur metabolism-related risk genes model by collecting multicenter clinical samples, and reveal the regulatory network of HSPA9 which impacted breast cancer cell malignant biological phenotypes and cancer cells development.
Recently, patients with breast cancer have increasingly benefited from early individualized treatments that depend on the investigation of cancer biomarkers, prognostic indicators, and prognostic model development [44–46]. The essence of the gene prognosis model is its ability to predict patients’ survival and treatment efficacy by analyzing gene expression in tumor specimens. For instance, certain genes related to ferroptosis are crucial in breast cancer prognosis, and the developed ferroptosis-based prognosis model can effectively distinguish between high-risk and low-risk patients, providing guidance for personalized therapy [47]. Moreover, using transcriptome data from the TCGA database, researchers identified a series of inflammation-related genes associated with breast cancer prognosis, and the expression patterns of these genes can be utilized to construct effective prognosis models [48]. Advancements in bioinformatics technology have rendered machine learning algorithms increasingly vital for identifying disease-related risk genes and developing disease prognosis models owing to their accuracy and reliability in prediction. Herein, we identified sulfur metabolism-related risk genes among patients with breast cancer for the first time through Cox regression analysis combined with the LASSO algorithm and developed a breast cancer prognosis model using the relevant risk genes. This enhances the biological markers of breast cancer for clinical practice, paving the way for targeted therapy of breast cancer, while also establishing a basis for future personalized treatment. However, our research still has some unresolved issues. Although our model performed excellently on specific datasets, its predictive effectiveness may be inconsistent across different populations or clinical environments. Previous studies have demonstrated that prognosis models built on specific genomic data may have good predictive capabilities on training datasets but fail to replicate these results on external validation sets [49]. Overfitting limits the clinical applications of the model, especially in diverse patient populations. As a result, enhancing the model’s generalization ability to maintain stable predictive performance in different clinical environments is essential for future research. Furthermore, many models lack sufficient clinical trial data to validate their effectiveness and reliability. Therefore, future research should focus on developing a localized gene bank for patients with breast cancer and on effective clinical validation of our sulfur metabolism-related risk gene model for breast cancer.

Conclusion
This study identified a reliable sulfur metabolism-related prognostic model for predicting the clinical outcomes and immune infiltration status of patients with breast cancer based on the genes ACOT2, ACOT4, CHPF, ELOVL2, HLCS, HSPA9, MICAL1, SPOCK2, and TCF7L2. In vitro analyses revealed that HSPA9 is a potential pro-oncogene and may serve as a novel target for breast cancer therapy.

Materials and methods

Patient’s expression profile data acquisition and preprocessing
We utilized the University of California, Santa Cruz (UCSC) XENA browser accessible via Genotype-Tissue Expression and TCGA databases (https://www.gtexportal.org/home/ and https://portal.gdc.cancer.gov/) to extract both breast invasive carcinoma (BRCA) samples and normal samples of RNA sequencing data for in-depth analysis. We applied the ‘rma’ function within the R package (R version: 3.6.3) to execute a log2 (TPM + 1) conversion on the sequencing data. Finally, we exclude missing values and duplicate results to the data.

Recognition of differentially expressed sulfur metabolism-related genes
Through the website (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), we identified a total of 332 genes associated with sulfur metabolism. These genes were sorted into gene sets labeled as GO:0051176, GO:0042762, GO:0000101, and GO:0006790. Subsequently, we utilized the R package DESeq2 to pinpoint differentially expressed genes related to sulfur metabolism between BRCA samples and normal samples. Genes meeting the criteria of P < 0.05 and |log2 fold change (FC)|> 2 were considered significantly different. Visual representations, including volcano plots, box plots, and heatmaps for the differentially expressed genes linked to sulfur metabolism, we first employed the R packages ggplot2 and ComplexHeatmap to generate visual representations. Following this, we utilized the String database (https://cn.string-db.org/) to establish a Protein-Protein Interaction Network (PPI). Key parameters of the network included assessing the significance of network edges (“evidence”), setting a minimum interaction score [“medium confidence (0.400)”], and determining the maximum number of reactors displayed (“no more than 50 reactors in the first layer”), and the active sources of interactions (“experiment”).

Identification of prognostic genes
In pursuit of genes correlated with Overall Survival (OS) in BRCA patients, we leveraged the R package ‘survival’, performed an analysis based on univariate Cox regression. Subsequently, we applied LASSO - Cox regression to address the potential impact of gene collinearity, resulting in the selection of more reliable genes. This was succeeded by a multivariate Cox regression analysis, drawing on the findings from the initial univariate Cox regression.

Construction and validation of prognostic model based on sulfur metabolism-related genes signature
Using LASSO Cox regression analysis, we developed a risk assessment model involving eight genes linked to sulfur metabolism, which holds relevance for predicting the prognosis of BRCA patients. The result was derived from the standardized mRNA expression profile specific to BRCA [50].
Within the scope of LASSO Cox regression analysis, X signifies the coefficients associated with pertinent metabolism-related genes, whereas Y represents the corresponding levels of gene expression. Breast cancer (BRCA) patients were then stratified into high and low-risk cohorts based on their average risk scores. Furthermore, the study investigated the disparity in OS between these groups. The model’s predictive performance can be assessed via the ggplot2 package.

Construction of nomogram and independent prognostic analysis
RMS package is used to establish the Nomogram for predicting individual survival results, the Nomogram model’s effectiveness was assessed through calibration plots, the ROC curve was utilized to assess the accuracy of both the risk score and various clinical features in predicting the outcomes of BRCA patients. This comprehensive evaluation provided insights into the model’s predictive performance.

Drug sensitivity analysis
We utilized the prophetic tool within the R software to calculate sensitivity scores of diverse small molecular compounds for individuals. This approach provides valuable insights into potential therapeutic options tailored to each patient group and identified 5 compounds in which IC50 in two groups has a significant difference. Then, we obtained the structure conformations of these five drugs from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) [51].

Enrichment analysis of differentially expressed genes based on risk score
Utilizing the DESeq2 package within R, we identified genes demonstrating significant disparities in expression levels between the high-risk and low-risk groups. These genes had to meet stringent criteria: P < 0.05 and |log2 fold change (FC)|> 2. Subsequently, we employed the ClusterProfiler and GOlot packages to scrutinize the genetic data, drawing from both Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Moreover, genome enrichment analysis (GSEA) was conducted to unveil distinct biological processes distinguishing the high-risk and low-risk groups. The reference gene sets were meticulously curated from the C2 subcollection (c5.go.all.v2022.1.Hs.symbols.gmt), providing a robust foundation for our analyses. All results were visualized by ggplot2 package.

Analysis of the association between sulfur metabolism-related genes and prognosis
We utilized Kaplan–Meier (KM) plots to evaluate the impact of genes associated with the sulfur metabolism signature on the prognosis of BRCA patients. These plots were generated through the application of the ggplot2 and survminer packages, providing clear visual representations of our findings.

Assessment of immune cell infiltration and the tumor microenvironment
The infiltration status of 28 immune cells with HSPA9 expression level in BRCA was estimated utilizing ssGSEA through the GSVA. It was calculated by Spearman correlation analysis.

Cell culture
The MDA-MB-231 cell line, representing human triple-negative breast cancer (TNBC), and the MCF-10 A human mammary epithelial cell line, were procured from the Peking Union Medical College Cell Bank in China. Additionally, the TNBC MDA-MB-468 cell line was sourced from Procell Life Science&Technology Co., Ltd (CL-0290B, China). Cultivation was conducted using DMEM medium (Oshima, NY, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco, Grand Island, NY, USA) for MDA-MB-231 and MDA-MB-468. For MCF-10 A, a specialized medium (CM-0525, China) was employed for maintenance.

Reverse transcription quantitative polymerase chain reaction (RT-qPCR)
Utilizing Trizol reagent, we proceeded to extract total cellular RNA. Post-extraction, NanoDrop 2000 was employed to evaluate both its quality and concentration. Subsequently, we employed the Hifair II 1st Strand cDNA Synthesis Kit for the synthesis of complementary DNA (cDNA). The next step involved conducting qRT-PCR utilizing SYBR green qPCR mix. The qRT-PCR procedure included cycles of heating and cooling, with specific primers used for target genes (HSPA9 and GAPDH). PCR primers are 5′- CTTGTTTCAAGGCGGGATTATGC − 3′ (forward) and 5′- GCAGGAGTTGGTAGTACCCAAA − 3′ (reverse) for HSPA9, 5′-GGACTCATGACCACAGTCCATG-3′ (forward) and 5′-CAGGGATGATGTTCTGGAGAGC-3′ (reverse) for GAPDH.

SiRNA transfection
Cells from MDA-MB-231 and MDA-MB-468 lines were placed in 6-well plates at a concentration of 1.5 × 10^5 cells per well and left to incubate overnight. Utilizing Lipofectamine 2000, small RNA interfering (siRNA) transfections were conducted. The siRNA specific to HSPA9 was acquired from RiboBio (Guangzhou, China), while a random siRNA was utilized as a control with no specific target. The siRNA sequence was as follows: siHSPA9-1: 5′-UUCAGGAUCAUCUUAUCGC − 3′, siHSPA9-2: 5′- UGACUAAGGCCAUUCCAGC − 3′.

CCK-8 assay
After the introduction of siRNA, we evaluated the proliferative capacity of MDA-MB-231 and MDA-MB-468 cells using a CCK-8 kit sourced from Lablead Biotech (China). To elaborate, siRNA-transfected MDA-MB-231 and MDA-MB-468 cells were seeded into 96-well plates (7000 cells per well). Then, 100µL of a 10% CCK-8 solution was introduced, and the absorbance was measured at 450 nm at time intervals of 0, 24, 48, and 72 h employing Tecan Infinite 200 pro. Subsequently, the acquired data underwent thorough analysis utilizing GraphPad Prism 9.0.

Western blotting
Cell proteins were confirmed through western blot analysis. To initiate cell lysis, we employed RIPA buffer and gauged the total protein concentrations utilizing a BCA Kit (Thermo, USA). Following this, proteins were subjected to separation. Subsequent to a one-hour blockade with 5% skim milk at room temperature, the membrane underwent an overnight incubation at 4 °C with specific primary antibodies. Following three TBST washes, the membrane was treated for an hour. The ultimate step involved visualizing the signals through an enhanced ECL chemiluminescent substrate kit (Yeasen Biotech). The quantification of target protein expression was meticulously executed utilizing the Image J tool, and the values were normalized to β-actin. Here are the details of the antibodies employed: HSPA9 (A11256, Abclonal, China), β-actin (#4967, CST), E-cadherin (A24874, Abclonal, China), N-cadherin (A0433, Abclonal, China), Vimentin (A2584, Abclonal, China) and HRP-conjugated goat anti-rabbit IgG.

Cell migration assessment
A wound-healing experiment was conducted to assess cellular motility. MDA-MB-231 and MDA-MB-468 cells were seeded in six-well plates after siRNA transfection and cultured until they reached confluence. Using a sterile 200 µL pipette tip, a controlled incision was made, followed by the removal of cellular debris using 1 × PBS. The incised area was then documented at intervals of 0 h, 24 h, and 48 h using a Leica DMi8 inverted microscope, and subsequently analyzed using Image J The rate of cell migration was computed as [(initial incision area - incision area at 24–48 h)/initial incision area] × 100%. Migration rates at distinct time points were visualized and assessed using GraphPad Prism 9.0.

Statistical analysis
Statistical differences were evaluated utilizing both the two-sided t-test and one-way ANOVA functions available in the GraphPad Prism software. Each experiment underwent a minimum of three repetitions.

Supplementary Information