Identification and validation of a prognostic signature of drug resistance and mitochondrial energy metabolism-related differentially expressed genes for breast cancer

Abstract

Background

Drug resistance constitutes one of the principal causes of poor prognosis in breast cancer patients. Although cancer cells can maintain viability independently of mitochondrial energy metabolism, they remain reliant on mitochondrial functions for the synthesis of new DNA strands. This dependency underscores a potential link between mitochondrial energy metabolism and drug resistance. Hence, drug resistance and mitochondrial energy metabolism-related differentially expressed genes (DMRDEGs) may emerge as candidates for novel cancer biomarkers. This study endeavors to assess the viability of DMRDEGs as biomarkers or therapeutic targets for breast cancer.

Methods

We utilized the DRESIS database and MSigDB to identify genes related to drug resistance. Additionally, we sourced genes associated with mitochondrial energy metabolism from GeneCards and extant literature. By merging these genes with differentially expressed genes observed in normal and tumor tissues from the TCGA-BRCA and GEO databases, we successfully identified the DMRDEGs. Employing unsupervised consensus clustering, we divided breast cancer patients into two distinct groups based on the DMRDEGs. Consequently, we identified four hub genes to formulate a prognostic model, applying Cox regression, LASSO regression, and Random Forest methods. Furthermore, we examined immune infiltration and tumor mutation burden of the genes within our model and scrutinized divergences in the immune microenvironment between high- and low-risk groups. Small hairpin RNA and lentiviral plasmids were designed for stable transfection of breast cancer cell lines MDA-MB-231 and HCC1806. By conducting clone formation, scratch test, transwell assays, cell viability assay and measurement of oxygen consumption we initiated a preliminary investigation into mechanistic roles of AIFM1.

Results

We utilized DMRDEGs to develop a prognostic model that includes four mRNAs for breast cancer. This model combined with various clinical features and critical breast cancer facets, demonstrated remarkable efficacy in predicting patient outcomes. AIFM1 appeared to enhance the proliferation, migration, and invasiveness of breast cancer cell lines MDA-MB-231 and HCC1806. Moreover, by reducing oxygen consumption, it aids in the cancer cells' acquisition of drug resistance.

Conclusions

DMRDEGs hold promise as diagnostic markers and therapeutic targets for breast cancer. Among the associated mutated genes, ATP7B, FUS, AIFM1, and PPARG could serve as early diagnostic indicators, and notably, AIFM1 may present itself as a promising therapeutic target.

Background

Breast cancer has emerged as one of the most common malignant tumors and stands as the primary contributor to cancer-related mortality among women [1]. Despite continuous advancements in breast cancer treatment modalities and significant enhancements in therapeutic efficacy, a subset of patients still experiences adverse outcomes—regardless of the comprehensive treatment strategy applied—due to recurrence, metastasis, and resistance to treatment [2,3,4]. This outcome is largely attributed to the heterogeneity exhibited by breast cancer. Therefore, an extensive body of research is imperative to augment our understanding of the disease's etiology, pinpoint risk factors, advance methods of early detection, and cultivate efficacious molecular biomarkers.

The metabolic processes in cancerous tissues are complex, affecting both the genesis and progression of tumors, thus making metabolic intervention a crucial strategy in cancer management [5]. Within the intricate network of tumor metabolic pathways, numerous potential therapeutic targets exist, each integrated into a nexus of connections with pathways involved in immunity, metastasis, drug resistance, and more [6]. The mitochondrion is pivotal in energy production within all cells, encompassing cell signaling, metabolism, apoptosis, and calcium homeostasis. The important role of mitochondrial metabolism lies in its ability to transmute the chemical energy stored in the three main nutrients into energy that is utilizable by body, such as heat and adenosine triphosphate(ATP) [7]. Aberrations in the mitochondrial energy metabolism pathway are associated with oxidative stress, characterized by the overproduction of free radicals,the release of copious amounts of reactive oxygen species(ROS), and the body's intrinsic antioxidative capacity,which leads to an overaccumulation of ROS and disrupts the balance between oxidative and antioxidative abilities. Typically, cancer cells manifest higher levels of ROS in comparison to their normal cellular counterparts [8, 9]. Mitochondrial respiration serves as a principal source of ROS, and elevated levels of ROS can inflict damage upon organelles [10]. Owing to their robust metabolic activity, cancer cells thrive even in hypoxic conditions. Owing to their incessant growth and insufficient vascular perfusion, cancer cells are susceptible to elevated levels of ROS, which can permeate the mitochondrial membrane and ultimately cause damage to DNA [11]. Therefore, mitochondria assume a pivotal role in the occurrence, development, and metastasis of cancers. Mitochondrial energy metabolism-related genes(MRGs) may be implicated in aberrant energy production, potentially contributing to the resistance observed in tumor treatments [12, 13].

The precise molecular mechanisms that confer drug resistance in breast cancer have yet to be fully determined. This investigation examined genes related to mitochondrial energy metabolism in breast cancer that may be associated with drug resistance, as well as the relationship between hub genes and tumor immune invasion and its potential mechanism, with the goal of screening specific biomarkers to predict therapeutic resistance. We have established a prognostic model of breast cancer based on these genes and further verified the model's accuracy in predicting outcomes for breast cancer. Complementary cytological studies have also corroborated the role of AIFM1 in breast cancer. These findings will provide promising data for breast cancer patient stratification, accurate prognosis evaluation, and personalized treatment.

Methods and materials

Data download

The R package TCGAbiolinks [14] was used to download RNA-seq data from the The Cancer Genome Atlas(TCGA) (https://portal.gdc.cancer.gov/). The TCGA-BRCA expression matrix contained 1118 breast cancer samples and 113 Normal samples. A total of 1231 samples were included in the differential analysis study and were standardized into TPM (Transcripts Per Million) format. Among 1118 breast cancer samples, 1003 samples with complete clinical information were included in the prognostic analysis study.

The breast cancer related datasets GSE42568 [15], GSE86374 [16] and GSE10886 [17] were downloaded from the GEO database through the R package GEOquery [18]. The samples were all from Homo sapiens. GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, respectively. GPL6244 [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version] and GPL1390 Agilent Human 1A Oligo UNC custom Microarrays. Among them, dataset GSE42568 contained 104 Tumor samples and 17 Normal samples. GSE86374 dataset contained 124 Tumor samples and 35 Normal samples. Dataset GSE10886 contains 190 Tumor samples and 7 Normal samples. All breast cancer samples and normal samples from the above GEO datasets were included in this study. The R package sva [19] was used to debatch the Datasets GSE42568, GSE86374 and GSE10886 to obtain Combined Datasets. The Combined Datasets included a total of 418 Tumor samples and 59 Normal samples. The Combined Dataset was standardized by R package limma [20], and the annotation probes were standardized and normalized. The Tumor and Normal samples of the Combined Dataset were included in this study as the validation set.

By DRESIS database [21] (http://dresis.idrblab.net/), download The drug-resistant genes associated set of "The general information of molecule associated with hold", 1784 Drug resistant-related Genes (DRGs) were obtained after deduplication. In the MSigDB database [22], using "drug AND resistan*" as the search keyword, 234 DRGs including 18 gene sets were collected. The above genes were merged to remove the duplication, and 2002 DRGs were obtained in total. GeneCards [23] database provides comprehensive information on human genes, using "Mitochondrial Energy Metabolism" as the search keyword, Relevance score > 0 as the threshold, A total of 282 MRGs were obtained. In addition, the set of MRGs published in the PubMed website was obtained by using the term "Mitochondrial Energy Metabolism" as the search keyword [24, 25], which contained a total of 239 MRGs. A total of 495 MRGs were obtained by combining the above genes to remove the weight. Detailed information is provided in Tables S1 and S2.

Differentially expressed genes related to drug resistance and mitochondrial energy metabolism

The R package DESeq2 [26] was used for differential analysis of genes in breast cancer samples (Tumor) and Normal samples (Normal) in TCGA-BRCA. Set | logFC | padj < > 0.5 and 0.05 for the Differentially Expressed Genes (DEGs) threshold., logFC padj < > 0.5 and 0.05 genes differentially expressed genes to increase (up—regulated genes), logFC < 0.5 and padj < 0.05 genes differentially expressed genes to cut (the down regulated genes).Intersect the DEGs with DRGs and MRGs to generate a Venn diagram, thereby identifying the DMRDEGs. The R package ggplot2 was used to draw the volcano map, the R package pheatmap was used to draw the expression heatmap, and the R package corrplot was used to draw the correlation heatmap. In addition, the R package RCircos [27] was used to annotalize the location of these genes in the human chromosome.

Mutation analysis and copy number variation analysis

To analyze Somatic Mutation of TCGA-BRCA, "Masked Somatic Mutation" data was selected as the somatic mutation data of TCGA-BRCA through TCGA website. Finally, the R package maftools [28] was used to visualize the somatic mutation situation.

To analyze the Copy Number Variation (CNV) in TCGA-BRCA, The "Masked Copy Number Segment" data was selected as the CNV data of TCGA-BRCA through TCGA website. Then, the downloaded and processed CNV fragments were analyzed by GISTIC2.0 [29], and all the default parameters were used in the GISTIC2.0 analysis process.

Functional and pathway enrichment analysis

Gene ontology (GO) [30] analysis is a common method for large-scale functional enrichment studies, including biological process (BP), cell component (CC), and molecular function (MF) terms, and the Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway analysis were carried out to predict the possible courses and pathways. The KEGG [31] is a widely used database storing information on genomes, biological pathways, diseases, and drugs. The R package ClusterProfiler [32] was used to perform GO and KEGG annotation analysis of DMRDEGs. The item screening criteria were p-value < 0.05 and FDR value (qvalue) < 0.25. The padj correction method was Benjamini-Hochberg.

Breast Cancer subtype construction

Consensus Clustering [33] is an algorithm based on resampling to identify each member and its subgroup number and verify the rationality of Clustering. Consensus Clustering is a series of iterations on subsamples of the data set, which provides an indicator of Cluster stability and parameter decision by using subsampling to induce sampling variability. The consensus Clustering method of R package ConsensusClusterPlus [34] was used to identify different disease subtypes of breast cancer. In this process, the maximum number of Clusters is set to 10, and 80% of the total samples are drawn in 100 replicates, with the parameters ClusterAlg = "pam" and distance = "minkowski" set. Then, the expression differences of DMRDEGs in different disease subtypes were verified by group comparison plot.

Gene set enrichment analysis (GSEA)

GSEA [35] is used to evaluate the distribution trend of the genes of a predefined gene set in the gene table ranked by their correlation with the phenotype, so as to judge their contribution to the phenotype. The logFC value of each gene was obtained by differential analysis between different disease subtypes of breast cancer samples in TCGA-BRCA. The R package ClusterProfiler was used to perform GSEA for all genes in TCGA-BRCA based on logFC values. The parameters used in GSEA were as follows: The seed was 2023, the number of computations was 500, and the minimum number of genes contained in each gene set was 1 and the maximum number of genes contained in each gene set was 10. Through MSigDB Database [22] in the access to c2. All. V2023.1. Hs. Symbols. The enrichment of GMT GSEA, The screening criteria of GSEA were padj < 0.05 and FDR value (qvalue) < 0.25.

Gene set variation analysis (GSVA)

GSVA [36] is a non-parametric unsupervised analysis method. It is mainly used to evaluate the enrichment results of gene set of microarray nuclear transcriptome by converting the expression matrix of gene set between different samples into the expression matrix of gene set between samples. Thus, we can evaluate whether different pathways are enriched in different samples. Through MSigDB Database access h.a ll. V2023.2. Hs. Symbols. The GMT gene set and using R package GSVA to TCGA—all the genes in the BRCA GSVA, The functional enrichment differences in different breast cancer disease subtypes were calculated. The screening criterion of GSVA was p-value < 0.05.

Immune infiltration analysis to compare the effects of immunotherapy

CIBERSORT [37] is based on the principle of linear support vector regression to deconvolute the transcriptome expression matrix, so as to estimate the composition and abundance of immune cells in mixed cells. The CIBERSORT algorithm, combined with LM22 feature gene matrix, was used to filter out the data with immune cell enrichment score greater than zero, and the specific results of immune cell infiltration matrix were finally obtained. The proportion of infiltration abundance of immune cells among different disease subtypes in TCGA-BRCA was displayed by stacking bar chart, and then the group comparison diagram between different disease subtypes was drawn to show the expression difference of LM22 immune cells in TCGA-BRCA. Combined with the TCGA-BRCA expression matrix, the correlation between immune cells and the correlation between immune cells and DMRDEGs were calculated, and displayed by the correlation heat map and the correlation scatter plot.

Immunotherapy and IPS analysis

Tumor microenvironment cells and the degree of infiltration of immune cells and stromal cells in the tumor have a significant impact on prognosis. To better understand the impact of genes involved in immunity and stromal cells on prognosis, the R package ESTIMATE [38] uses the unique nature of the transcriptional profile of cancer samples to infer the content of tumor cells as well as different infiltrating normal cells, and the expression profile data of TCGA-BRCA is used to evaluate the immune activity of tumors through the R package ESTIMATE package. Quantitative analysis of immune activity (immune infiltration level) in tumor samples based on gene expression profiling to obtain an immune score for each tumor sample. The differences of immune infiltration characteristics between patients with different disease subtypes were compared.

Immunogenicity refers to the ability of an antigen or its epitopes to act on antigen recognition receptors of T cells and B cells to induce humoral and/or cell-mediated immune responses. This principle of immunity can be called Immunogen. Immunogenicity is determined by several genes, including those associated with effector cells, MHC molecules, immunomodulatory factors, and immunosuppressive cells. The Cancer Immunome Atlas (TCIA) database [39] (https//tCIa.at /home) provides immunophenoscores (IPS) for 20 cancers, It can be a good predictor of CTLA-4 and PD-1 reactivity. The IPS data of TCGA-BRCA were downloaded from TCIA database, and the R package ggplot2 was used to draw group comparison maps based on different disease subtypes of TCGA-BRCA to analyze the differences in IPS.

TMB analysis and TIDE analysis

To obtain the TIDE immunoscore results of TCGA-BRCA, the analysis was performed based on the expression matrix of TCGA-BRCA by the TIDE (Tumor Immune Dysfunction and Exclusion) algorithm [40, 41]. The TIDE immunoscore predicted the potential tumor treatment response of differential genes related to DMRDEGs, and tested the effect of the expression level of each gene in the tumor on patient survival. The differences of TIDE immunoscore in different disease subtypes were calculated according to the analysis results of TIDE immunoscore.

Subsequently, Tumor Mutation Burden (TMB) data were downloaded from the cBioPortal database [42,43,44]. The TIDE Immune score was calculated through the TIDE website. Finally, TMB and TIDE scores of samples with different disease subtypes were compared between groups.

Drug sensitivity analysis

Alterations in the cancer genome strongly influence the clinical response to therapy and are in many cases effective biomarkers of response to drug therapy. The Genomics of Drug Sensitivity in Cancer (GDSC) database [45] (www.cancerRxgene.org) is the largest public resource for information on molecular markers of drug sensitivity and drug response in cancer cells. The cancer Drug Sensitivity Genomics database can be used to find cancer drug response data and genomic sensitive markers. The pRRophetic algorithm [46] was used to predict the sensitivity of patients with different disease subtypes to common anticancer drugs or small molecule compounds by calculating the IC50 value based on the TCGA-BRCA, and the results were displayed by grouping comparison maps.

Protein interaction network

Protein–Protein Interaction Network (PPI Network) is composed of proteins interacting with each other. The STRING database [47] is a database for searching the interactions between known and predicted proteins. Through DMRDEGs, the protein interaction network of DMRDEGs was constructed with the minimum correlation coefficient greater than 0.150 as the standard. Cytoscape software [48] was used to visualize the network model. In addition, the GeneMANIA database [49] (http://genemania.org/) uses a very large set of functional association data to find other genes related to a set of input genes. Association data include protein and genetic interactions, pathways, co-expression, co-localization, and protein domain similarity. GeneMANIA database was used to analyze the protein–protein interaction network of DMRDEGs.

LASSO regression analysis and screening and validation of key genes

Glmnet package [50] was used to identify differentially expressed genes related to DMRDEGs in TCGA-BRCA with random seed number 2023, Ten-fold cross-validation method was used to perform Least absolute shrinkage and selection operator (LASSO) regression to further screen genes. LASSO regression is often used to construct prognostic models. On the basis of linear regression, the penalty term (lambda × absolute value of slope) is added to reduce the overfitting of the model and improve the generalization ability of the model. The results of LASSO regression analysis were visualized by prognostic risk model map and variable trajectory map, and the DMRDEGs contained in them were the key genes. LASSO Risk Score was calculated based on the risk coefficient of LASSO regression analysis, and the risk score was calculated using the following formula:

To validate the LASSO regression model, breast cancer samples were divided into High Risk group and Low Risk group according to the median Risk Score of the model. Group comparison maps and risk factor maps were drawn based on the expression levels of key genes. Time-dependent Receiver Operating Characteristic(ROC) Curve [51] is a coordinate schema-based analysis tool that can be used to select the best model, discard the second-best model, or set the best threshold in the same model. the R package survival ROC was used to generate time-dependent ROC curves based on key genes and Risk scores and calculate the Area Under the Curve (AUC) value to predict 1-year breast cancer. Survival effects at 3 and 5 years.

Survival analysis

Prognostic Kaplan–Meier (KM) curve analysis [52] is a method to analyze and infer the survival time of patients based on data, and to study the relationship between survival time and outcome and many influencing factors and their degree, also known as survival analysis or survival analysis. It is proposed by Kaplan and Meier, so it is called the Kaplan–Meier method, usually referred to as the KM method. The KM method estimates the survival curve by calculating the probability that a patient who has survived for a certain period will survive for the next period (i.e., the survival probability), and then multiplying the survival probabilities one by one, which is the survival rate of the corresponding period. The KM curve was drawn to verify the influence of key genes and Risk Score on survival.

Human protein atlas database and IHC validation

An examination of AIFM1's protein expression in breast cancer and adjacent normal mammary tissue samples was conducted utilizing the Human Protein Atlas online repository (https://www.proteinatlas.org/).

Prognostic clinical correlation analysis

In order to analyze the correlation between the Risk Score and the prognosis of patients, the relationship between the Risk Score and the prognostic pathological characteristics of clinical cases was evaluated. The effect of clinical stage, T stage, M stage, ER and PR immunohistochemical characteristics on the Risk Score was analyzed, and the differences were compared between groups. Then, univariate Cox regression analysis was performed on the key genes and clinicopathological characteristics, and the factors with p-value < 0.05 were selected and included in the multivariate Cox regression analysis to construct the multivariate Cox regression model. Based on the results of multivariate Cox regression analysis, a Nomogram was established to predict the 1 -, 3 -, and 5-year survival of breast cancer patients. A nomogram is a graph that uses a Cluster of disjoint line segments to represent the functional relationship between multiple independent variables in a rectangular coordinate system. Based on multivariate regression analysis, a certain scale is set to represent the situation of each variable in the multivariate regression model, and finally the total score is calculated to predict the probability of the occurrence of an event. Finally, the calibration curve was used to evaluate the accuracy and resolution of the nomogram. The Calibration curve is used to evaluate the prediction effect of the model on the actual outcome by plotting the fitting of the actual probability and the predicted probability of the model under different conditions. It is mainly used for the fitting analysis of the model established by Cox regression method and the actual situation. The horizontal axis of the calibration curve is the survival probability predicted by the model, and the vertical axis is the survival probability shown by the actual data. The lines and points of different colors represent the situation predicted by the model at different time points. The closer the lines of different colors are to the line of the gray ideal situation, the better the prediction effect at this time point. The R package rms was used to construct the nomogram and calibration curve. Decision curve analysis (DCA) [53] is a simple method to evaluate clinical prediction models, diagnostic tests, and molecular markers. The R package ggDCA was used to draw DCA diagram to evaluate the predictive effect of nomogram model on the survival outcome of breast cancer patients at 1, 3, and 5 years. The results can be determined based on the observation that the line of the model can be stable higher than the x value range of All positive line and All negative line. The larger the range of x value, the better the performance of the model.

Cell culture and transfection

The HCC1806 and MDA-MB-231 human breast cancer cell lines were acquired from the American Type Culture Collection (ATCC, Manassas, USA). The cells were propagated in RPMI-1640 (Sigma–Aldrich, Missouri, USA) and DMem (Corning, NY, USA), both enriched with 10% fetal bovine serum (Excell Bio, Suzhou, China) and 1% antibiotic–antimycotic solution, maintained at 37 °C with 5% CO2 in a humidified incubator.

The AIFM1 shRNA were purchased from Tsingke Biotech. (Beijing, China). The sequence of the human AIFM1-specific shRNA was 5’-GCCAAACTATTCAACATTCAT-3’ and 5’-CCTGGAAATAGACTCAGATTT-3’, and 5’-AGATTTCACGGGAAGTCAAAT-3’. The cancer cells, at 30% confluency, were seeded in 6-well plates containing serum-free medium. AIFM1 shRNA or the negative control(NC) shRNA from Tsingke Biotech (Beijing, China) was transfected into these cells using polybrene (Beyotime, Beijing, China). Forty-eight hours post-transfection, the cells were subjected to puromycin selection to establish stable cell lines with the transfection. The efficacy of the knockdown was confirmed through western blot analysis.

Western blotting

Whole cell lysates were prepared with RIPA lysis buffer (NCM Biotech, Suzhou, China) with a protease (NCM Biotech, Suzhou, China) and a phosphatase inhibitor cocktail(NCM Biotech, Suzhou, China). The Pierce BCA portein assay(Thermo Scientific, Massachusetts, USA) was used to measure protein concentration in the supernatant after removing the precipitate. Electrophoresis was used to transfer the protein samples onto PVDF membranes(Millipore, Massachusetts, USA) after boiling the protein. Following the blocking procedure with skim milk, the membranes were incubated overnight at 4℃ with primary antibodies. Subsequently, the membranes were washed thrice with TBST and then incubated with secondary antibodies for 1 h at room temperature. The antibodies used in this study were as follows: anti-AIFM1 (Proteintech, 17984–1-AP), anti-β-Actin (Proteintech, 60009–1-Ig), HRP-conjugated Affinipure Goat Anti-Rabbit IgG (Proteintech, SA00001-2) and HRP-conjugated Affinipure Goat Anti-Mouse IgG(H + L) (Proteintech, SA00001-1).

Colony formation assay

In order to evaluate the cancer cells’ capabilities for colony formation, a plate cloning assay was carried out. The cohort of 600 transfected cells were uniformly distributed across a 6-well plate and subsequently cultured over a span of 14 days with intermittent renewal of the medium. The culture medium was discarded, fixed with 4% paraformaldehype (Biosharp, Beijing, China) for 30 min, stained with crystal violet for 20 min, gently rinsed with running water, air-dried, and photographed for subsequent counting. Utilize the ImageJ software (http://rsb.info.nih.gov/ij/download. html) to ascertain the count of colony formations.

Migration and invasion assays

The scratch assay evaluates the migratory potential of cancer cells. A 6-well plate was seeded with the transfected cells. Upon attaining 80% cell density, a pipette tip was positioned perpendicular to the base of the plate and employed to create a scratch in the cell monolayer. Results were documented with photographs at 0 h, 24 h, and 48 h. Gap distance was quantified using ImageJ software.

The transwell assay evaluates the invasive potential of cancer cells. After pre-coating the Transwell chamber with Matrigel(Corning, NY, USA), 400 transfected cells were resuspended in fresh basal medium and subsequently added to the upper chamber. Full medium was introduced into the lower chamber. After 48 h, cells from the upper chamber were meticulously removed. The cells that remained were fixed with 4% paraformaldehyde, stained with crystal violet, and then imaged under a microscope. The number of invasive cells was quantified using ImageJ software.

Measurement of oxygen consumption

Cancer cells were seeded into each well of a 96-well black-bottom plate(3*10⁴ cells in 100 μl of culture medium per well) and incubated overnight. Oxygen Consumption Rate(OCR) was measured using the OCR Plate Assay Kit (Dojindo, Kumamoto, Japan). The fluorescent signals were measured at 10 min intervals. The calculation of OCR was derived from an analysis of the kinetic profiles obtained from measurements.

Cell viability assay

Cancer cells were seeded into each well of a 96-well plate (2000 cells in 100 μl of culture medium per well) and incubated overnight. After discarding the supernatant, 100 μl of culture medium were added to the control wells, and 100 μl of the culture medium containing various concentrations of Lapatinib(Targetmol,Massachusetts,USA) were introduced into the experimental wells, followed by incubation for 48 h. For the cell growth curve assay, cell Counting Kit-8 (Glpbio,Montclair,USA) was added to the wells for 1 h, then cell numbers were measured at 450 nm.

Statistical analysis

All data processing and analysis in this article were conducted utilizing R software version 4.2.1, or GraphPad Prism 9. Continuous variables are presented as mean ± SD. The Wilcoxon Rank Sum Test was used for comparison between the two groups. If not specified, the correlation coefficient between different molecules was calculated by Spearman correlation analysis, and a p-value of less than 0.05 was used as the criterion for significant difference.

Results

Differentially expressed genes related to drug resistance and mitochondrial energy metabolism

In this study, a total of 1118 BRCA samples from TCGA and 477 BRCA samples from the GEO database (GSE42568, GSE86374, GSE10886) were analyzed. The samples were divided into Tumor group and Normal group. To evaluate differential gene expression between the Tumor and Normal groups, we utilized the R package DESeq2 for differential analysis. Within the dataset, our analysis uncovered 3492 DEGs, among which 2125 were upregulated and 1367 downregulated, in adherence to the criteria of an absolute log fold change (|logFC|) greater than 0.5 and an adjusted p-value (padj) less than 0.05. Subsequently, volcano maps visually represented the differential expression analysis outcomes. To enhance clarity and facilitate reference, the key genes identified during the subsequent screening process were annotated on the volcano plot (Fig. 1A). DEGs were intersected with DRGs and MRGs to obtain the DMRDESGs. 15 genes were identified, and a Venn diagram was generated (Fig. 1B). These 15 genes are: ATP7B, SIRT6, FUS, UCP2, AIFM1, PFKL, IL1B, PTEN, IRS1, PKD2, PTGS2, ALDH1A3, FOXO1, PFKFB3, PPARG. A heatmap compared DMRDEGs expression among different groups within the TCGA-BRCA dataset and the Combined dataset (Fig. 1C-D). A correlation heatmap depivted the relationships among the DMRDEGs within the TCGA-BRCA dataset (Fig. 1E). Finally, the location of 15 DMRDEGs on the human chromosome was analyzed, and a chromosome localization map was drawn (Fig. 1F). The chromosome mapping showed that more DMRDEGs were located in chromosomes 2, 10 and 13, respectively: IL1B, IRS1 in chromosome 2, PFKFB3, PTEN in chromosome 10, FOXO1 and ATP7B in chromosome 13.

Fig. 1

Differential Gene Expression Analysis. A Volcano plot of differentially expressed genes analysis between Tumor group and Normal group in TCGA-BRCA. B Venn diagram of intersection of DEGs with DRGs and MRGs in TCGA-BRCA. C-D Expression heatmap of DMRDEGs in TCGA-BRCA and Combined Dataset. E Correlation heat map of DMRDEGs in TCGA dataset. F Chromosomal mapping of DMRDEGs. *** indicates p-value < 0.001, ** indicates p-value < 0.01, * indicates p-value < 0.05

Mutation analysis and copy number variation analysis

To investigate the somatic mutation profiles of 15 DMRDEGs (ATP7B, SIRT6, FUS, UCP2, AIFM1, PFKL, IL1B, PTEN, IRS1, PKD2, PTGS2, ALDH1A3, FOXO1, PFKFB3, PPARG) within the TCGA-BRCA dataset, mutation analysis was conducted on breast cancer patient samples, and the results were compiled and analyzed. These findings reveal that the primary somatic mutations in TCGA-BRCA are Missense Mutation, Nonsense Mutation, and Frame Shift Del, as well as splice site mutations, Frame Shift Ins mutations, In Frame Del mutations, and others, with missense mutations being the most frequent. In addition, DMRDEGs in BRCA patients mainly feature single nucleotide polymorphisms (SNPs), along with a number of deletions (DELs) and insertions (INS). C > T emerges as the most common single nucleotide variant (SNV) in breast cancer patients, succeeded by C > G and C > A (Fig. 2A). All 15 DMRDEGs in TCGA-BRCA present somatic mutations, with 96 out of 991 somatic mutation samples implicated, accounting for 9.69% of the total. Notably, the PTEN gene exhibits the highest mutation rate, constituting 5% of all mutation samples (Fig. 2B). The mutation status within TCGA-BRCA is depicted using stacked bar charts (Fig. 2C). An association is identified among the 15 DMRDEGs between FOXO1 and PPARG, PFKFB3 and PPARG, and also between PFKFB3 and FUS(p-value < 0.05) (Fig. 2D). From the KM curves of 15 DMRDEGs in TCGA-BRCA, it can be deduced that samples harboring mutations are linked to a worse DSS prognosis (Fig. 2E).

Fig. 2

Mutation Analysis of DMRDEGs in TCGA-BRCA. A Mutation presentation in TCGA-BRCA. B Mutation waterfall plot in DMRDEGs in TCGA-BRCA. C Boxplot and stacked bar plot of TCGA-BRCA mutation types. D Mutation correlation heat map of DMRDEGs in TCGA-BRCA. E KM curve of mutation of DMRDEGs and DSS survival in TCGA-BRCA patients. F CNV of DMRDEGs in TCGA-BRCA

GISTIC2.0 analysis was employed to determine the CNV of TCGA-BRCA, revealing that all 15 DMRDEGs exhibit CNVs. The values span from greatest GAIN to greatest LOSS in the following order: PTGS2, FUS, PFKFB3, PFKL, UCP2, PPARG, ALDH1A3, AIFM1, SIRT6, PKD2, IL1B, FOXO1, ATP7B, PTEN, IRS1 (Fig. 2F).

GO and KEGG analysis

Through enrichment analysis of GO and KEGG, we delved deeper into the connection between BP, CC, MF, and KEGG pathway of 15 DMRDEGs and breast cancer. Then, a combined logFC GO and KEGG enrichment analysis was conducted, which was based on GO and KEGG enrichment analysis. By providing the logFC values of DMRDEGs in differential analysis results, the z-score corresponding to each GO and KEGG entry was calculated. The specific results are shown in Table S3. The results showed that these 15 DMRDEGs were mainly enriched in BP such as negative regulation of transport,regulation of hormone secretion, negative regulation of protein transport,regulation of protein secretion, and negative regulation of establishment of protein localization in breast cancer; in CC, such as caveola, plasma membrane raft, Schmidt-Lanterman incisure,compact myelin, and cytoplasmic side of membrane; in MF, such as NAD + binding, carbohydrate kinase activity, alpha-actinin binding, actinin binding, and NAD binding. Additionally, they were enriched in pathways such as AMPK signaling pathway, Central carbon metabolism in cancer, Longevity regulating pathway,Insulin resistance, and FoxO signaling pathway. The results of the GO and KEGG enrichment analysis were visualized through bar charts and bubble charts (Fig. 3A-B), and the GO and KEGG enrichment analysis results of the joint logFC were displayed through chord and circle plots (Fig. 3C-D).

Fig. 3

GO and KEGG Enrichment Analysis for DMRDEGs. A Bar graph of GO and KEGG enrichment analysis results of joint logFC of DMRDEGs: BP, CC, MF and Pathway. B Bubble plot of GO and KEGG enrichment analysis results of combined logFC of DMRDEGs. C Chordplot of GO and KEGG enrichment analysis results of joint logFC for DMRDEGs. D Circle plot of GO and KEGG enrichment analysis results of combined logFC of DMRDEGs. The screening criteria for GO and KEGG enrichment analysis were p-value < 0.05 and FDR value (q-value) < 0.25

Consensus cluster analysis of DMRDEGs

To investigate the association between breast cancer subtypes and DMRDEGs, we conducted a consensus cluster analysis using 15 genes. As a result, two distinct BRCA subtypes were identified: Subtype 1 (Cluster 1) and Subtype 2 (Cluster 2) (Fig. 4A-C). Among them, subtype 1 (Cluster1) contained 865 samples and subtype 2 (Cluster2) contained 253 samples. To further verify the expression differences of DMRDEGs in BRCA disease subtypes, Boxplots (Fig. 4D) display the expression levels of DMRDEGs and the variations between the two BRCA subtypes. Group comparison plots reveal significant disparities in the expression levels of SIRT6, FUS, PTEN, and PFKFB3 (p-value < 0.001), as well as PKD2 (p-value < 0.01) across the disease subtypes. ATP7B, PTEN, and PFKFB3 exhibited significant differences between the disease subtypes (p-value < 0.001). The differences in the expression of IL1B, IRS1, and FOXO1 among the disease subtypes were statistically significant (p-value < 0.05).

Fig. 4

Consensus Clustering Analysis for DMRDEGs. A Plot of consensus Clustering results for BRCA. B-C. Cumulative distribution function (CDF) plot (B) and Delta plot (C) of concordance Clustering analysis. D Boxplot of differential genes related to DMRDEGs between the two BRCA subtypes. *** denotes p-value < 0.001, ** denotes p-value < 0.01, and * denotes p-value < 0.05

GSEA enrichment analysis

To assess the impact of gene expression levels on breast cancer subtypes in TCGA-BRCA, we initially compared differences among disease subtypes in breast cancer samples and acquired logFC values for each gene. Using GSEA, the expression of all genes in TCGA-BRCA was examined, including BP, CC, and MF. The detailed findings are shown in Table S4. The results showed that all genes in TCGA-BRCA were significantly enriched in NIKOLSKY BREAST CANCER 7Q21 Q22 AMPLICON (Fig. 5A), DACOSTA UV RESPONSE VIA ERCC3 DN (Fig. 5B), WP GPCRS CLASS A RHODOPSINLIKE (Fig. 5C), NIKOLSKY BREAST CANCER 8Q23 Q24 AMPLICON (Fig. 5D) and other biologically related functions and signaling pathways. Then, the above pathways were visualized by drawing a mountain map (Fig. 5E).The findings indicated that enrichment was associated with immune-related functions and pathways.

Fig. 5

GSEA for TCGA-BRCA Dataset. A-D GSEA showed that BRCA subtypes significantly affected NIKOLSKY BREAST CANCER 7Q21 Q22 AMPLICON (A), DACOSTA UV RESPONSE VIA ERCC3 DN (B), WP GPCRS CLASS A RHODOPSINLIKE (C), NIKOLSKY BREAST CANCER 8Q23 Q24 AMPLICON (D). E. 4 biological function mountain map of GSEA for different disease subtypes in TCGA-BRCA. The screening criteria of GSEA were padj < 0.05 and FDR value (q-value) < 0.25

GSVA enrichment analysis

To explore differences in the h.all.v2023.2. Hs.symbols.gmt gene set among breast cancer subtypes, GSVA was conducted on all genes in the TCGA-BRCA dataset. Refer to Table S5 for detailed information. The top 10 logFC positive pathways with padj < 0.05 and the top 10 logFC negative pathways with padj < 0.05 were screened respectively, and the enrichment scores of each sample (Fig. 6A) were shown. Boxplots of the Combined Datasets for group comparison (Fig. 6B) were also displayed. GSVA results showed the top 10 pathways with positive logFC values were: HALLMARK E2F TARGETS, HALLMARK G2M CHECKPOINT, HALLMARK IL6 JAK STAT3 SIGNALING, HALLMARK MYC TARGETS V2, HALLMARK SPERMATOGENESIS, HALLMARK ALLOGRAFT REJECTION, HALLMARK MYC TARGETS V1, HALLMARK TNFA SIGNALING VIA NFKB, HALLMARK MTORC1 SIGNALING, HALLMARK NOTCH SIGNALING; The top 10 pathways with negative logFC values are: HALLMARK TGF BETA SIGNALING, HALLMARK FATTY ACID METABOLISM, HALLMARK HEME METABOLISM, HALLMARK ESTROGEN RESPONSE LATE, HALLMARK PANCREAS BETA CELLS, HALLMARK PEROXISOME, HALLMARK PROTEIN SECRETION, HALLMARK UV RESPONSE DN, HALLMARK BILE ACID METABOLISM, HALLMARK ESTROGEN RESPONSE EARLY.

Fig. 6

GSVA for TCGA-BRCA. A Heat map of enrichment scores for TCGA-BRCA GSVA. B Boxplot of group comparison of GSVA results among BRCA subtypes. *** indicates p-value < 0.001

CIBERSORT immune infiltration analysis

The CIBERSORT algorithm was used to compute the correlation between 22 immune cells and breast cancer subtypes 1 (Cluster1) and 2 (Cluster2). According to the results of immune infiltration analysis, the superimposed bar graph of immune cells in TCGA-BRCA was drawn (Fig. 7A). Then, the expression difference of immune cell infiltration abundance in subtype 1 (Cluster1) and subtype 2 (Cluster2) of TCGA-BRCA was shown by grouping comparison boxplot (Fig. 7B). The results indicated that B cells naive, Plasma cells, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, The abundance of Neutrophils in subtype 1 (Cluster1) and subtype 2 (Cluster2) in TCGA-BRCA showed statistically significant differences (p-value < 0.05). Then, the correlation results of immune cell infiltration abundance were presented by a correlation heat map (Fig. 7C). The results showed the greatest positive correlation between T cells CD8 and T cells regulatory (Tregs)(cor = 0.385, p-value < 0.001) and the greatest negative correlation between NK cells resting and NK cells activated (cor = −0.730, p-value < 0.001). The correlation between DMRDEGs and the abundance of immune cell infiltration in TCGA-BRCA was shown by a correlation dot plot(Fig. 7D). The results showed the greatest positive correlation between IL1B and Mast cells activated(cor = 0.502, p-value < 0.001), and the greatest negative correlation between ATP7B and Macrophages M1 (cor = −0.318, p-value < 0.001), and the relationship between genes and the abundance of immune cell infiltration was shown by a correlation scatter plot (Fig. 7E-F).

Fig. 7

Immune Infiltration Analysis by CIBERSORT Algorithm. A Stacked bar chart of immune cell infiltration abundance in BRCA subtype 1 (Cluster1) and subtype 2 (Cluster2). B Comparison of immune cell infiltration abundance between BRCA subtype 1 (Cluster1) and subtype 2 (Cluster2). C Heat map of infiltration abundance correlation between immune cells. D Dot plot of correlation between differential genes related to DMRDEGs and abundance of immune cell infiltration. E–F Scatter plot of gene-immune cell infiltration abundance correlations with the strongest negative (E, ATP7B and Macrophages M1) and the strongest positive (F, IL1B and Mast cells activated) correlations. *** indicates p-value < 0.001, ** indicates p-value < 0.01, * indicates p-value < 0.05. The absolute value of the correlation coefficient (cor) below 0.3 was weak or no correlation, 0.3–0.5 was weak correlation, 0.5–0.8 was moderate correlation, and above 0.8 was strong correlation

Analysis of immunotherapy responses

The ESTIMATE package leverages the distinct characteristics of cancer sample transcriptome to estimate the composition of tumor cells and various infiltrating normal cells. It primarily computes the immune and stromal scores based on RNA seq data of the sample, and subsequently assesses the tumor purity. The expression profile data of different disease subtypes in TCGA-BRCA were analyzed with the R package ESTIMATE, leading to varied immune and matrix scores, Stromal Score, Immune Score, ESTIMATE Score, and Tumor Purity were obtained for the different disease subtypes of TCGA-BRCA samples. Display the scoring results of Stromal Score, Immune Score, ESTIMATE Score, and Tumor Purity through a group comparison chart (Fig. 8A-D). The results indicated that Stromal Score and Immune Score show significant statistical differences among various disease subtypes(p-value < 0.05). However, There is no significant statistical difference between ESTIMATE Score and Tumor Purity in different disease subtypes (p-value > 0.05).

Fig. 8

ESTIMATE Analysis Between Different BRCA Clusters. A-D Violin plot for comparison of different disease subtypes by Stromal Score (A), Immune Score (B), ESTIMATE Score (C), and Tumor Purity (D) score. *** indicates p-value < 0.001

IPS, TMB, and TIDE analyze

In order to evaluate the predictive power of breast cancer subtypes for immunotherapy efficacy, the IPS related to TCGA-BRCA were retrieved from the TCIA database. The R package ggplot2 was utilized to construct the group comparison plots for various IPS across the breast cancer subtypes (Fig. 9A-D). The analysis disclosed that within the IPS classification, CTLA4(−)PD1( +), CTLA4( +)PD1(−), and CTLA4( +)PD1( +) subtypes showed remarkably significant differences (p-value < 0.001).

Fig. 9

IPS, TMB, TIDE Analysis. D Group comparison of IPS between groups of different BRCA disease subtypes: CTLA4(-)PD1(-) (A), CTLA4(-)PD1( +) (B), CTLA4( +)PD1(-) (C), CTLA4( +)PD1( +) (D). E Group comparison of TMB scores for different disease subtypes of BRCA. F Group comparison of TIDE scores by BRCA disease subtype. *** indicates p-value < 0.001

Furthermore, an examination was conducted into the variance of TMB scores across different breast cancer subtypes, with the results visualized in a group comparison chart (Fig. 9E). This investigation revealed that TMB scores significantly diverged across the subtypes(p-value < 0.001). Subsequently, the TIDE algorithm was applied to assess how sensitive breast cancer patients are to immunotherapy; these specific analytical outcomes were also depicted in a group comparison plot (Fig. 9F). It was found that the TIDE scores, pertaining to immunotherapy potency, did not present notable differences among the breast cancer subtypes (p-value > 0.05).

Drug sensitivity analysis

In the quest to identify suitable mRNA vaccine therapy approaches for breast cancer patients, sensitivity data to prevalent anticancer drugs was harnessed from the GDSC database as a predictive training set for patient responsiveness. A subsequent evaluation spotlighted the disparate drug sensitivities between two breast cancer subtypes. Enumerated here are the top 20 pharmaceutical agents that manifested significant differentiation between the subtypes: MK.2206, Lapatinib, AZD8055, WO2009093972, GDC0941, Temsirolimus, EHT.1864, GW.441756, CCT007093, FH535, PF.4708671, PD.0332991, Elesclomol, AKT.Inhibitor.VIII, Pazopanib, IPA.3, Axitinib, Metformin, NVP.BEZ235, and AMG.706 (Fig. 10A-T). These 20 drugs show statistically significant differences between two subtype groups (p-value < 0.001).

Fig. 10

Drug Sensitivity Analysis. A-T. Comparison of drugs MK.2206 (A), Lapatinib (B), AZD8055 (C), WO2009093972 (D), GDC0941 (E), Temsirolimus (F), EHT.1864 (G) among different disease subtypes based on GDSC database. GW.441756 (H), CCT007093 (I), FH535 (J), PF.4708671 (K), PD.0332991 (L), Elesclomol (M), AKT.inhibitor.VIII (N), Pazopanib (O), The group comparison figure of sensitivity analysis results of IPA.3 (P), Axitinib (Q), Metformin (R), NVP.BEZ235 (S), and AMG.706 (T) is shown. *** indicates p-value < 0.001

Protein interaction network

A PPI network comprising 15 DMRDEGs was constructed utilizing the STRING database (Fig. 11A). The PPI network analysis divulged that there were interactions among 14 out of the 15 DMRDEGs, specifically: ATP7B, SIRT6, FUS, UCP2, AIFM1, PFKL, IL1B, PTEN, IRS1, PTGS2, ALDH1A3, FOXO1, PFKFB3, and PPARG. Additionally, the GeneMANIA database revealed associated genes for these DMRDEGs, which were identified as: ELF5, HTR1B, KLF15, KMT5A, BAK1, GSTA2, PTGS1, SLC25A1, RELA, ATOX1, ATP7A, HAX1, TNNI3, BC1L11, FABP4, NR2F2, NOXA1, TRIM63, AKT1, and PFKFB2 (Fig. 11B).

Fig. 11

PPI Network. A PPI Network of DMRDEGs calculated from STRING database. B PPI interaction network of DMRDEGs in GeneMANIA database

LASSO regression analysis and screening and validation of key genes

Employing clinical data from TCGA-BRCA combined with DMRDEGs, we engaged in a LASSO logistic regression analysis, endeavoring to establish a regression model and to contrive a prognostic risk assessment tool. We present both the LASSO coefficient progression (Fig. 12A) and the constructed model (Fig. 12B) visually. Analysis of the LASSO regression model pinpoints four pivotal genes: ATP7B, FUS, AIFM1, and PPARG, heralded as crucial determinants. The risk score calculation formula for each sample is as follows:

$${\text{Risk Score}}\, = \,\left( { - 0.0{4787}} \right) \, *{\text{ ATP7B}}\, + \,\left( { - 0.{1344}} \right) \, *{\text{ FUS}}\, + \,0.{2892 }*{\text{ AIFM1}}\, + \,\left( { - 0.0{29}0{3}} \right) \, *{\text{ PPARG}}$$

Fig. 12

LASSO Analysis for DMRDEGs. A-B Plots of prognostic risk models (A) and variable trajectories (B) from the LASSO regression model. C-D Comparison of expression groups of key genes in TCGA-BRCA and Combined Dataset. E Risk factor plot of LASSO regression model. F 1 -, 3 -, and 5-year time-dependent ROC curves of Risk scores. The abscissa is the false positive rate (FPR) and the ordinate is the true positive rate (TPR). LASSO, Least absolute shrinkage and selection operator; *** indicates p-value < 0.001

Samples were segregated into high-risk and low-risk categories predicated upon the median value of the Risk Scores. Subsequently, we compared the expression levels of key genes between subjects in the TCGA-BRCA and a Combined Dataset, stratified by risk categorization (Fig. 12C-D). The results show that nearly all key genes display significant expression differences in TCGA-BRCA and Combined Dataset (p-value < 0.001). Among them, ATP7B, FUS, and PPARG show higher expression in low-risk samples, whereas AIFM1 exhibits higher expression in high-risk samples. This suggests that AIFM1 could be a prognostic risk factor in breast cancer, while ATP7B, FUS, and PPARG may serve as protective factors (Fig. 12E). Finally, draw a risk factor map for the LASSO regression results of TCGA-BRCA (Fig. 12F). The results indicate that the LASSO model for breast cancer exhibits high accuracy at one year(AUC > 0.7) but lower accuracy at three and five years(AUC > 0.5).

Survival analysis

Grouping TCGA-BRCA according to the expression levels and Risk Score of four key genes, and plotting the survival prognosis KM curve (Fig. 13A-E). The results showed a significant difference in survival prognosis between the high and low-risk score groups (p-value < 0.001), with the low-risk score group demonstrating a better prognosis than the high-risk score group. There was a highly significant difference in survival prognosis between the high and low expression groups of AIFM1 (p-value < 0.01), with the low expression group having a better prognosis compared to the high expression group; There was a significant contrast in survival prognosis between the high and low expression groups of ATP7B (p-value < 0.05), with the high expression group showing better prognosis than the low expression group.

Fig. 13

Survival Analysis for Key Genes and Risk Score. A-D KM curve of survival prognosis in TCGA-BRCA for different expression levels of key genes. E KM curves of survival prognosis for high and low Risk Score groups in TCGA-BRCA

Clinical prognosis and its corresponding factors

To assess the correlation between risk score and patient prognosis, we evaluated the relationship between risk score and clinical pathological features predictive of outcome. The prognostic impact of higher versus lower risk scores was appraised across several domains: clinical staging, T staging, M staging, ER/PR immunohistochemistry, and additional tumoral pathological characteristics. As we are well aware, there is a significant difference in the Risk Score under these clinical pathological feature groups (p-value < 0.05) (Fig. S1). Visualize the relationships among various clinical pathological characteristics through Sankey diagrams (Fig. 14A). Subsequently, univariate Cox regression analysis was conducted between these clinical pathological characteristics and four key genes, with the results plotted in a forest plot (Fig. 14B). Factors with p < 0.05 were then included in a multivariate Cox regression analysis. Due to the presence of zero and infinite values in the multi-factor hazard ratios for the M stage, this clinical pathological characteristic was also excluded. Therefore, the factors included in the multivariate Cox regression include: Pathlogic stage, ER immunohistochemistry, PR immunohistochemistry, FUS gene, and AIFM1 gene. Draw a column chart to show the contribution of these factors to the prognostic model based on the results of multiple factor regression (Fig. 14C). Subsequently, prognostic calibration curves were plotted to assess the accuracy and precision of the Cox regression model (Fig. 14D). Additionally, decision curve analysis (DCA) was conducted to determine the model's prognostic utility (Fig. 14E). The results showed that the 1-year, 3-year, and 5-year prognostic calibration curves were close to the diagonal of the ideal model, indicating a very high accuracy of the model. In DCA, a model demonstrates superior performance when its line consistently surpasses both the all-positive and all-negative thresholds within a given range. The larger this range, the greater the net benefit achieved by the model. The results show that the line of the 1-year, 3-year, and 5-year models is consistently higher than that of All positive and All negative models within a certain range, and the model has a higher net profit, indicating better performance. The model reveals that the AIFM1 gene is the most significant contributor, indicating poor prognosis for samples with high AIFM1 expression (HR > 1). Based on the TCGA-BRCA database, these four genes manifest differential expression in tumor patients (Fig. 1A). To demonstrate the increased expression of AIFM1 at the tissue protein level in cancer tissues compared to normal breast tissues, we referred to immunohistochemical results from the HPA database and obtained clinical samples for immunohistochemical experiments. The results showed a significantly higher protein expression level in breast cancer tissues compared to normal mammary tissues(Fig. S2). Finally, a baseline data table was drawn to examine the relationship between the high and low expression of the AIFM1 gene and clinical variables, as shown in Table S6.

Fig. 14

Prognostic Model for Cox Regression. A Sankey plot of clinicopathological features. B Clinical characteristics, single factor Cox regression analysis forest plot of four key genes. C-E Prognostic nomogram of multivariate Cox regression model (C), prognostic calibration curve (D), DCA plot (E)

Knocking down AIFM1 mitigates the malignant phenotype of tumor cells and reduces drug resistance

Western blot analysis confirmed the successful knockdown of AIFM1 (Fig. 15A). The silencing of AIFM1 significantly inhibited the proliferation of MDA-MB-231 and HCC1806 cells, as shown by the colony formation assays (Fig. 15B-C). Additionally, scratch assays and Transwell experiments indicated that AIFM1 knockdown led to a considerable decrease in the migratory and invasive capacities of tumor cells (Fig. 15D-G). OCR measurements elucidated an inversely proportional relationship between AIFM1 expression levels and cellular respiration (Fig. 15H). The previous analysis (Fig. 10) identified Lapatinib as one of the top 20 compounds with significant variances and as a potent therapeutic agent for advanced breast cancer. Furthermore, results from the CCK-8 assays have demonstrated that downregulation of AIFM1 enhances the anti-tumor efficacy of Lapatinib (Fig. 15I). In summary, our research demonstrates that knocking down AIFM1 can reduce the virulence and drug resistance of tumor cells.

Fig. 15

Experimental validation of the role of AIFM1 in breast cancer cell lines. A Knockdown efficiency of AIFM1 in MDA-MB-231 and HCC1806 cells. B-C The proliferation ability of MDA-MB-231 and HCC1806 cells was measured by colony formation assay after transfecting AIFM1 shRNAs. D-E Scratch experiment for measuring the migration ability of MDA-MB-231 and HCC1806 cells. F-G. MDA-MB-231 and HCC1806 cell transwell invasion assay were performed in control and sh AIFM1 groups. H The OCR of the MDA-MB-231 and HCC1806 cell lines was assessed following induction by the control and sh AIFM1 group. I Treatment of cancer cells with escalating concentrations of Lapatinib (0.3, 0.6, 0.9, 1.2, 1.5, and 1.8 µM) over a 48-h period resulted in the determination of cell viability, as quantified by CCK-8 assays. *p-value < 0.05; **p-value < 0.01; ****p-value < 0.0001

Discussion

Breast cancer, a heterogeneous disease characterized by widely varying prognoses and therapeutic responses, can be identified through gene or biomarker expression analyses, which are predictive of prognosis. In clinical practice, drug resistance stands as a formidable obstacle, with the examination of tumor resilience through the prism of cellular metabolism eliciting significant interest in the field of cancer research. Previous studies have highlighted that MEM reprogramming stands out as a unique hallmark of cancer. This reprogramming promotes glucose uptake and allows cancer cells to prefer glycolysis as their primary energy source even in normal oxygen conditions. The lactate produced from glycolysis sustains the acidic microenvironment of cancer cells, promoting invasion, metastasis and drug resistance, which correlates with poor prognosis [54,55,56]. Based on this, mitochondrial mutations may correlate with tumor drug resistance and could potentially serve as early breast cancer diagnostic biomarkers. We have revealed that DMRDEGs carry prognostic value, enhancing clinical parameters for predicting prognosis.

Initially, by harnessing public databases and scrutinizing extant literature, we identified a cohort of 15 DMRDEGs. Subsequent functional analyses divulged that these genes are preponderantly engaged in NAD + binding and carbohydrate kinase activity. Moreover, they are overrepresented in pathways including the AMPK signaling pathway and central carbon metabolism in cancer. Within malignant milieus, NAD + conducts a spectrum of crucial functions, encompassing the sophisticated modulation of immune evasion strategies and the exertion of a significant impact upon the tumor microenvironment [57]. Variations in carbohydrate kinase activity are commonly associated with the neoplasm's capabilities concerning growth, survival, proliferation, and migratory potential [58]. The AMPK pathway potentially curbs tumor growth by constraining energy supply while concurrently aiding tumor cell adaptation to metabolic stress, thus bolstering survival [59]. Moreover, central carbon metabolism is frequently reprogrammed to meet the exigent demands of tumor cell growth and division [60]. In summary, these molecular functions and signaling pathways are integral to cancer pathogenesis and progression. Here, the GO and KEGG analyses indicated that 15 DMRDEGs play significant roles in maintaining cellular homeostasis, genomic stability, cell growth and death, and immune response.

Based on 15 DMRDEGs, BRCA has been stratified into two distinct Clusters: Cluster 1 and Cluster 2. GESA and GSVA results show that differentially expressed genes between Cluster 1 and Cluster 2 are distributed among a variety of signaling pathways, collectively encompassing several important carcinogenic trajectories. The emergence of immune checkpoint inhibitors and adoptive cellular immunotherapy has improved the prognosis of breast cancer patients [61, 62]. Determining patient cohorts to immunotherapy is crucial. CIBERSORT results make it evident that, with a reduced count of M0 macrophages, a larger number of macrophages in Cluster 1 are induced to acquire the M2 phenotype. TAMs are correlated with unfavorable prognoses, facilitating treatment resistance, including resistance to immunotherapy [63]. TAMs represent the largest contingent of immune cells within tumors, predominantly comprising two polarized states: the oncogenic M2 macrophages and the antitumoral M1 macrophages [64]. Consequently, Cluster 2 exhibits enhanced responsiveness to immunotherapy, in alignment with subsequent analyses. In contrast, Cluster 1 is characterized by a proliferation of immunosuppressive regulatory T cells (Tregs) relative to Cluster 2. Tregs are instrumental in facilitating tumor immune evasion. By targeting Tregs, especially those residing within the tumor microenvironment, one might amplify the effectiveness of cancer immunotherapies [65]. Moreover, Cluster 1 not only demonstrates a reduced response to immunotherapy but also exhibits resistance to several targeted drugs, thus amplifying the treatment intricacies. Collectively, these findings hold promise for providing more tailored therapeutic strategies for patients with diverse profiles, thereby optimizing personalized treatment protocols and diminishing the overall patient burden.

Following this, we constructed an influential predictive risk model anchored in four hub genes (including ATP7B, FUS, AIFM1, and PPARG), which has emerged as a good prognostic instrument. The risk model had the excellent ability to predict the survival of breast cancer, and possesses the potential to pave the path towards tailored immunotherapies for breast cancer patients susceptible to drug resistance. Survival analysis reveals that amongst the four genes, a disparity in the expression levels of AIFM1 corresponds with the most pronounced prognostic difference in survival outcomes between the groups. Moreover, among the four hub genes, AIFM1 is identified as a poor prognostic factor. The extant literature indicates that AIFM1 is principally implicated in orchestrating the processes of Parthanatos and Oxeiptosis—distinct modalities of programmed cell death. Parthanatos is predicated on the hyperactivation of poly-ADP-ribose polymerase-1 (PARP-1), while oxeiptosis is primarily induced by oxidative stress. however, the expression level of AIFM1 and its role in tumors may vary depending on the type of cancer, the tumor microenvironment, the stage of cancer progression, and other biological factors [66,67,68,69]. Our study has revealed an overexpression of AIFM1 in breast cancer cells. Functional in vitro assays indicate that AIFM1 inhibition decreases the proliferative, migratory, and invasive capabilities of breast cancer cells. Moreover, compared to AIFM1-silenced neoplastic cells, the normal tumor cells exhibit a diminished oxygen consumption rate and an enhanced resistance to Lapatinib. We hypothesize that this observation may be due to metabolic reprogramming, which promotes cellular survival under stress by reducing mitochondrial respiration and increasing glycolysis. Consequently, AIFM1 presents considerable promise as a biomarker for the surveillance of breast cancer recurrence, the assessment of treatment efficacy, and the prediction of patient survival outcomes.

This study has some limitations. The transcriptomic data utilized are exclusively derived from public databases. It is essential to conduct transcriptomic sequencing on breast cancer patients in real-world settings for further analysis and validation. Furthermore, the putative role of AIFM1 in the progression of breast cancer has been provisionally substantiated by a scant number of in vitro experiments, thus its physiological effects in vivo remain obscure and beckon further investigative efforts. In particular, supplementary animal studies designed to investigate the functional role of AIFM1 in cancer would provide robust insights to guide clinical applications. Current literature intimates that AIFM1 is implicated as a tumor suppressor across a range of malignancies, attributable chiefly to its capacity to precipitate apoptosis within cancerous cells [70,71,72]. However, in addition to inducing apoptosis, AIFM1 is pivotal in sustaining the functionality of mitochondrial oxidative phosphorylation (OXPHOS), a metabolic process that contributes to the maintenance of stem-like qualities in tumor cells, thereby endowing them with oncogenic potency and resistance to therapeutic agents [73, 74]. In our study, we found that the upregulation of AIFM1 expression is associated with an enhanced resistance to the therapeutic agent lapatinib, correlating with the concomitant reduction in OCR observed. Considering the integral role of AIFM1 in maintaining functional mitochondrial OXPHOS this could elucidate the oncogenic role of AIFM1 in breast cancer. Taking into account the convoluted nature of human physiology, the prospect of surmounting drug resistance in breast cancer by targeting AIFM1 is contingent upon further scrutiny into the cardinal functions of AIFM1 and the nuances of the tumor microenvironment specific to breast malignancies.

Conclusion

Our team has developed a prognostic model intricately associated with DMRDEGs that possesses the capability to predict with high precision the varied responses of breast cancer patients to multiple therapeutic strategies. Additionally, knockdown experiments of the AIFM1 gene imply that AIFM1 holds potential as an innovative target for the treatment of breast cancer.