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

doi:10.21203/rs.3.rs-4757124/v1

Download PDF

Research Article

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

https://doi.org/10.21203/rs.3.rs-4757124/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background:

Drug resistance constitutes one of the principal causes of poor prognosis in breast cancer patients. Cancer cells can survive independently of the energy provided by mitochondria; however, they are incapable of synthesizing new DNA strands without mitochondrial involvement.This may suggest that mitochondrial energy metabolism could be related to 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 the 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 the immune infiltration and tumor mutation burden of the genes within our model and scrutinized the divergences in the immune microenvironment between high- and low-risk groups. Small hairpin RNA and lentiviral plasmids were designed for the stable transfection of breast cancer cell lines MDA-MB-231 and HCC1806. By conducting clone formation, scratch test and transwell assays, we initiated a preliminary investigation into the mechanistic roles of AIFM1.

Results:

We utilized DMRDEGs to develop a prognostic model that includes four mRNAs for breast cancer, which, by combining various clinical features and critical breast cancer facets, proved to be remarkably effective in forecasting patient outcomes. Additionally, AIFM1 appeared to enhance the proliferation, migration, and invasiveness of the breast cancer cell lines MDA-MB-231 and HCC1806.

Conclusions:

DMRDEGs have the potential to act as diagnostic markers and therapeutic targets for breast cancer. Within the mutated associated genes, ATP7B, FUS, AIFM1, and PPARG could serve as early diagnostic indicators, and notably, AIFM1 may present itself as a promising therapeutic target.

breast cancer

mitochondrial energy metabolism

drug resistance

prognostic model

Breast cancer has emerged as one of the most common malignant tumor and concurrently 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 within cancerous tissues operate in complex ways, influencing both the genesis and the advancement of tumors, making metabolic intervention an important strategy in the management of cancer[5]. Within the sprawling labyrinth of metabolic pathways in tumors lie countless prospective therapeutic targets, each woven into a web of connections that intertwine with pathways involved in immunity, metastasis, drug resistance, among others[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. Due to continuous growth and inadequate vascular perfusion, cancer cells are prone to an increase in ROS levels, which can diffuse through mitochondrial membrane and ultimately inflicting damage upon 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 molecular mechanism of drug resistance in breast cancer has yet to be fully determined. This investigation examined the genes related to mitochondrial energy metabolism in breast cancer that may be associated with drug resistance, as well as the relationship between hub gene and tumor immune invasion and its potential mechanism, in order to screen 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. Cytological experiments confirmed the role of AIFM1 in breast cancer. These findings will provide promising data for breast cancer patient stratification, accurate prognosis evaluation, and personalized treatment.

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:

$$\:Risk\:Score=\underset{i}{?}coefficient\:\left(gene\:i\right)*expression\:\left(gene\:i\right)$$

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.

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 1h 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 minutes, 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 0h, 24h, and 48h. 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 hours, 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.

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.

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 the 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 were constructed to visually represent the outcomes of the differential expression analysis. To enhance clarity and facilitate reference, the key genes identified during the subsequent screening process were annotated on the volcano plot (Fig. 1A). Intersect the DEGs 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 was generated to compare the expression levels of DMRDEGs among different groups within the TCGA-BRCA dataset and the Combined dataset (Fig. 1C-D). A correlation heatmap was then drawn to depict 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 the 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.

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, a mutation analysis was conducted on breast cancer patient samples, and the results were compiled and analyzed. The results indicate that the primary somatic mutations in TCGA-BRCA were: Missense Mutation, Nonsense Mutation, and Frame Shift Del, along with splice site mutations, Frame Shift Ins mutations, In Frame Del mutations, etc. Missense mutations are the most common among them. Moreover, the mutation types of DMRDEGs in BRCA patients primarily include single nucleotide polymorphisms (SNPs), as well as some deletions (DELs) and insertions (INS). C>T is the most common single nucleotide variant (SNV) in breast cancer patients, followed by C>G and C>A(Fig. 2A).All 15 DMRDEGs in TCGA-BRCA exhibited somatic mutations. Out of the 991 somatic mutation samples, 96 samples were implicated in the mutation, representing 9.69% of the total. Among them, the PTEN gene has the highest mutation rate, representing 5% of all mutation samples(Fig. 2B). Visualize the mutation status of TCGA-BRCA using stacked bar charts(Fig. 2C).Among the 15 DMRDEGs, there is a correlation between FOXO1 and PPARG, PFKFB3 and PPARG, as well as PFKFB3 and FUS(p-value < 0.05)(Fig. 2D).Based on the KM curves of 15 DMRDEGs in TCGA-BRCA, it can be inferred that samples with mutations have a worse DSS prognosis(Fig. 2E).

GISTIC2.0 analysis was utilized to examine the CNV in TCGA-BRCA, where 15 DMRDEGs all show CNVs . The values range from greatest GAIN to greatest LOSS: PTGS2, FUS, PFKFB3, PFKL, UCP2, PPARG, ALDH1A3, AIFM1, SIRT6, PKD2, IL1B, FOXO1, ATP7B, PTEN, IRS1, respectively(Fig. 2F).

GO and KEGG analysis

Through enrichment analysis of GO and KEGG, we delve 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 the GO and KEGG enrichment analysis. By providing the logFC values of DMRDEGs in the 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,negative regulation of establishment of protein localization in breast cancer; in CC, such as caveola, plasma membrane raft,Schmidt-Lanterman incisure,compact myelin,cytoplasmic side of membrane; in MF, such as NAD+ binding, carbohydrate kinase activity, alpha-actinin binding, actinin binding, 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).

Consensus Cluster Analysis of DMRDEGs

To ascertain 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, The boxplot (Fig. 4D) was used to display the expression levels of DMRDEGs and the differences between the two BRCA subtypes. Group comparison plots indicated that the expression levels of SIRT6, FUS, PTEN, and PFKFB3 were significantly different between the disease subtypes (p-value < 0.001), and the expression levels of PKD2 were significantly different between the disease subtypes (p-value < 0.01). 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).

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. By 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.

GSVA enrichment analysis

To explore the 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. The boxplots of Combined Datasets for group comparison (Fig. 6B) were also displayed. The results of GSVA showed that 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.

CIBERSORT immune infiltration analysis

The CIBERSORT algorithm was used to calculate 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 showed 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 shown by correlation heat map (Fig. 7C). The results showed that T cells CD8 and T cells regulatory (Tregs) showed the greatest positive correlation (cor = 0.385, p-value < 0.001). There was the greatest negative correlation between NK cells resting and NK cells activated (cor = -0.730, p-value < 0.001). The correlation between differentially expressed genes related to drug resistance and mitochondrial energy metabolism (DMRDEGs) and the abundance of immune cell infiltration in TCGA-BRCA was shown by correlation dot plot (Fig. 7D). The results showed that IL1B showed the greatest positive correlation with Mast cells activated (cor = 0.502, p-value < 0.001), and ATP7B showed the greatest negative correlation with Macrophages M1 (cor = -0.318, p-value < 0.001), and the relationship between genes and the abundance of immune cell infiltration was shown by correlation scatter plot (Fig. 7E-F).

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).

IPS, TMB, and TIDE analyze

To analyze how well breast cancer subtypes predict immunotherapy, TCGA-BRCA-associated IPS were downloaded from the TCIA database, The R package ggplot2 was used to draw the group comparison map of different IPS in breast cancer among different disease subtypes (Fig. 9A-D). The results showed that CTLA4(-)PD1(+), CTLA4(+)PD1(-) of breast cancer subtypes in IPS class, CTLA4(+) pd1 (-), CTLA4(+) pd1 (-), CTLA4(+)PD1(+) showed a highly statistically significant difference (p-value < 0.001).

The difference of tumor mutation burden (TMB) score in different disease subtypes of breast cancer was analyzed and the group comparison diagram was drawn (Fig. 9E). The results showed that the TMB score in breast cancer was highly statistically significant among different disease subtypes (p-value < 0.001). Then, the sensitivity of breast cancer patients to immunotherapy was evaluated by the TIDE algorithm, and the specific analysis results are shown in the group comparison plot (Fig. 9F). The results showed that there were no significant differences in TIDE immunotherapy scores between breast cancer subtypes (p-value > 0.05).

Drug sensitivity analysis

To explore therapeutic strategies suitable for mRNA vaccination in breast cancer patients, drug sensitivity data from the GDSC database were used as a training set to predict the sensitivity of breast cancer patients to common anticancer drugs. Subsequently, the sensitivity difference of two subtypes of breast cancer to different anticancer drugs was evaluated. Here we list the top 20 drugs showing significant variations among different disease 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，AMG.706 (Fig. 10A-T).These 20 drugs show statistically significant differences between two subtype groups(p-value < 0.001).

Protein interaction network

A protein-protein interaction Network (PPI Network) of 15 drug resistance and mitochondrial energy metabolism-related differential genes (DMRDEGs) was constructed using STRING database (Fig. 11A). The results of protein-protein interaction Network (PPI Network) showed that 14 of the 15 drug resistance and mitochondrial energy metabolism-related differential genes (DMRDEGs) had interaction relationships, which were: ATP7B, SIRT6, FUS, UCP2, AIFM1, PFKL, IL1B, PTEN, IRS1, PTGS2, ALDH1A3, FOXO1, PFKFB3, PPARG. Through GeneMANIA database, the genes associated with 15 drug resistance and mitochondrial energy metabolism-related differential genes (DMRDEGs) were obtained, which were: ELF5, HTR1B, KLF15, KMT5A, BAK1, GSTA2, PTGS1, SLC25A1, RELA, ATOX1, ATP7A, HAX1, TNNI3, BCL1L11, FABP4, NR2F2, NOXA1, TRIM63, AKT1, PFKFB2 (Fig. 11B).

LASSO regression analysis and screening and validation of key genes

Using clinical data from TCGA-BRCA along with DMRDEGs for LASSO regression analysis, we will build a LASSO regression model and develop a prognostic risk model. Visualize LASSO variable trajectory maps (Fig. 12A) and LASSO regression model maps (Fig. 12B). The results indicated that the LASSO regression model included four genes: ATP7B, FUS, AIFM1, and PPARG, which are considered key genes. The risk score calculation formula for each sample is as follows:

Risk Score = (-0.04787) * ATP7B + (-0.1344) * FUS + 0.2892 * AIFM1 + (-0.02903) * PPARG

The sample is then divided into high and low-risk groups based on the median Risk Score. Next, comparing the expression levels of key genes between TCGA-BRCA and Combined Dataset by high and low risk grouping(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.

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. Evaluate the influence of high or low Risk Score on clinical staging, T staging, M staging, ER/PR immunohistochemistry, and other tumor pathological features on prognosis. 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 HR 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 decision curve analysis (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). 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.

AIFM1 knockdown inhibited cellular proliferation, migration, and invasion in breast cancer within the HCC1806 and MDA-MB-231 cell lines.

Western blot analysis corroborated the effective silencing of AIFM1(Fig. 15A). The suppression of AIFM1 markedly curtailed the proliferation of MDA-MB-231 and HCC1806 cells, as evidenced by colony formation assay(Fig. 15B-C). Furthermore, scratch assays and Transwell experiments have yielded results indicating that the migratory and invasive abilities of tumor cells exhibited a notable decline subsequent to the knockdown of AIFM1(Fig. 15D-G). Overall, our findings indicated that AIFM1 may promote the proliferation, migration, and invasion of breast cancer.

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. The issue of drug resistance poses a significant challenge in clinical practice, with the exploration of tumor resistance through cell metabolism gaining considerable attention in oncological research. Previous studies have shown 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, thereby facilitating invasion, metastasis and drug resistance, which correlateswith a poor prognosis[54][55][56]. Based on this, the mitochondrial mutations may correlate with tumor drug resistance and could potentially be used as early breast cancer diagnostic biomarkers. we have revealed that DMRDEGs carry prognostic value, enhancing clinical parameters for predicting prognosis.

First of all, through the utilization of public databases and a review of existing literatures, we ascertained a cohort of 15 DMRDEGs. Subsequent functional analysis demystified that molecular functions of these genes are primarily involved in NAD + binding and carbohydrate kinase activity. At the same time, it is also enriched in pathway such as AMPK signaling pathway,Central carbon metabolism in cancer. Within malignant milieus, NAD + orchestrates an array of pivotal functions including the nuanced regulation of immune subterfuge and exerting substantial influence over 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 curbing tumor growth by constraining energy supply while concurrently aiding tumor cell adaptation to metabolic stress, thus bolstering survival[59]. Moreover, the central carbon metabolism is frequently subject to reprogramming to fulfill the exigent demands for tumor cell growth and division[60]. So, in summation, these molecular functions and signaling pathways are intimately linked with cancer pathogenesis and its evolution. 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 the differentially expressed genes between Cluster 1 and Cluster 2 are distributed among a variety of signaling pathways, which collectively encompass several important carcinogenic trajectories. The emergence of immune checkpoint inhibitors and adoptive cellular immunotherapy has improved the prognosis of breast cancer patients[61][62]. It is crucial to determine which patient cohorts will respond effectively to immunotherapy. According to the results of CIBERSORT, it is discernible that based on a lower count of M0 macrophages, a higher quantity of macrophages in Cluster1 are stimulated to differentiate into the M2 phenotype. Tumor-associated macrophages (TAMs) are believed to be associated with the poor prognosis of tumors, attributable to their role in conferring resistance to treatments, including immunotherapy[63]. TAMs constitute the most plentiful immune cell population within neoplasms, consisting chiefly of two polarization states: the tumorigenic M2 macrophages and the tumor-suppressing M1 macrophages[64]. Hence, Cluster 2 demonstrates an enhanced inclination towards immunotherapy, consistent with the findings of subsequent analyses. Cluster 1 exhibits a higher concentration of Immunosuppressive regulatory T cells (Tregs) in comparison to Cluster 2. Tregs provide a main mechanism of tumor immune evasion. Targeting Tregs, particularly those within the tumor microenvironment, may enhance the efficacy of cancer immunotherapy[65]. Apart from exhibiting a poor response to immunotherapy, the drug sensitivity analysis has also revealed that Cluster1 possesses resistance to a number of targeted medications, thus exacerbating the complexities of its treatment. In sum, the aforementioned outcomes can offer more precise therapeutic recommendations for diverse patients, thereby facilitating personalized treatment and reducing the burden of patients.

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. howeve, 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 revealed a significant elevation in the expression of AIFM1 within breast cancer cells. The results of functional experiments in vitro showed AIFM1 knockdown reduces cell proliferation, migration, and invasion in breast cancer.Therefore, AIFM1 holds promise as a potential biomarker for breast cancer recurrence, treatment monitoring, and survival.

This study has some limitations. The transcriptomic data employed in this study is sourced exclusively from public databases. It is essential to conduct transcriptomic sequencing on breast cancer patients in real-world settings for further analysis and validation. Moreover, the role of AIFM1 in the progression of breast cancer has been tentatively corroborated only by a limited number of in vitro experiments; its in vivo effects remain elusive, warranting additional investigation. Further corroboration of drug sensitivity through cellular experimentation, and additional animal studies to investigate the functional role of AIFM1 in breast cancer would provide robust insights for guiding clinical applications.

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.

DMRDEGs, Drug Resistance & Mitochondrial Energy Metabolism-Related

Differentially Expressed Genes;

ATP, Adenosine Triphosphate;

ROS, reactive oxygen species;

MRGs, Mitochondrial Energy Metabolism-Related Genes;

TCGA, The Cancer Genome Atlas;

TPM, Transcripts Per Million;

DRGs, Drug Resistance-Related Genes;

DEGs, Differentially Expressed Genes;

CNV, Copy Number Variations;

GO, Gene Ontology;

KEGG, Kyoto Encyclopedia of Genes and Genomes;

BP, Biological Process;

CC, Cellular Component;

MF, Molecular Function;

CDF, Empirical Cumulative Distribution Function;

GSEA, Gene Set Enrichment Analysis;

GSVA, Gene Set Variation Analysis;

TCIA, The Cancer Immunome Atlas ;

IPS, Immunophenoscore;

TIDE, Tumor Immune Dysfunction and Exclusion;

TMB, Tumor Mutation Burden;GDSC, Genomics of Drug Sensitivity in Cancer;

PPI, Protein-Protein Interaction;

LASSO, Least absolute shrinkage and selection operator;

ROC, Receiver Operating Characteristic;

AUC, Area Under the Curve;

KM, Kaplan-Meier;

DCA, Decision Curve Analysis；

ATCC, American Type Culture Collection;

NC, negative control;

SNPs, single nucleotide polymorphisms;

DELs, deletions;

INS, insertions;

SNV, single nucleotide variant;

FPR, False Positive Rate;

TPR, True Positive Rate;

MEM, mitochondrial energy metabolism;

TAMs, Tumor-associated macrophages;

Tregs, Immunosuppressive regulatory T cells;

PARP-1, poly-ADP-ribose polymerase-1;

Acknowledgements

We extend our profound gratitude to Professor Xia Yunfei for his guidance and assistance in our work.

Author contributions

TX , CC and HL contributed to the study conception and design. TX and SX wrote the manuscript. TX, CC and SX collected and analysed the raw data. TX, WX, TJ,and YW contributed to the cytological experiments. All authors read and approved the final manuscript.

Funding

This study was financially supported by National Natural Science Foundation of China (Grant No.82102838), National Natural Science Foundation of China (Grant No. 82103437), Basic and Applied Basic Research Foundation of Guangzhou (Grant No. 202201010918) and the Natural Science Foundation of Guangdong Province, China (No.2023A1515010685).

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Consent for publication

All authors confirm their consent for publication the manuscript.

Competing interests

The authors declare that there is no conflict of interest.

Author details

¹ Department of Radiation Oncology, State Key Laboratory of Oncology in South China, Guangdong Key Laboratory of Nasopharyngeal Carcinoma Diagnosis and Therapy, Guangdong Provincial Clinical Research Center for Cancer, Sun Yat-sen University Cancer Center, Guangzhou 510060, P. R. China

² Department of Medical Oncology, State Key Laboratory of Oncology in South China, Guangdong Key Laboratory of Nasopharyngeal Carcinoma Diagnosis and Therapy, Guangdong Provincial Clinical Research Center for Cancer, Sun Yat-sen University Cancer Center, Guangzhou 510060, P. R. China

Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, Jemal A: Global cancer statistics 2022: GLOBOCAN estimates of incidence and mor tality worldwide for 36 cancers in 185 countries. CA Cancer J Clin:10.3322/caac.21834.
Loibl S, Poortmans P, Morrow M, Denkert C, Curigliano G: Breast cancer. Lancet, 397(10286):1750-1769.
Harbeck N, Penault-Llorca F, Cortes J, Gnant M, Houssami N, Poortmans P, Ruddy K, Tsang J, Cardoso F: Breast cancer. Nat Rev Dis Primers, 5(1):66.
Anand U, Dey A, Chandel AKS, Sanyal R, Mishra A, Pandey DK, De Falco V, Upadhyay A, Kandimalla R, Chaudhary A et al: Cancer chemotherapy and beyond: Current status, drug candidates, assoc iated risks and progress in targeted therapeutics. Genes Dis, 10(4):1367-1401.
Schmidt DR, Patel R, Kirsch DG, Lewis CA, Vander Heiden MG, Locasale JW: Metabolomics in cancer research and emerging applications in clinical oncology. CA Cancer J Clin, 71(4):333-358.
Li X, Sun T, Jiang C: Intelligent Delivery Systems in Tumor Metabolism Regulation: Exploring the Path Ahead. Adv Mater, 36(11):e2309582.
Nunnari J, Suomalainen A: Mitochondria: in sickness and in health. Cell, 148(6):1145-1159.
Boese AC, Kang S: Mitochondrial metabolism-mediated redox regulation in cancer progressi on. Redox Biol, 42:101870.
Starkov AA: The role of mitochondria in reactive oxygen species metabolism and sig naling. Ann N Y Acad Sci, 1147:37-52.
Rossmann MP, Hoi K, Chan V, Abraham BJ, Yang S, Mullahoo J, Papanastasiou M, Wang Y, Elia I, Perlin JR et al: Cell-specific transcriptional control of mitochondrial metabolism by T IF1γ drives erythropoiesis. Science, 372(6543):716-721.
Vera-Ramirez L, Ramirez-Tortosa M, Perez-Lopez P, Granados-Principal S, Battino M, Quiles JL: Long-term effects of systemic cancer treatment on DNA oxidative damage : the potential for targeted therapies. Cancer Lett, 327(1-2):134-141.
Wang H, Shi W, Zeng D, Huang Q, Xie J, Wen H, Li J, Yu X, Qin L, Zhou Y: pH-activated, mitochondria-targeted, and redox-responsive delivery of paclitaxel nanomicelles to overcome drug resistance and suppress metas tasis in lung cancer. J Nanobiotechnology, 19(1):152.
Egan G, Khan DH, Lee JB, Mirali S, Zhang L, Schimmer AD: Mitochondrial and Metabolic Pathways Regulate Nuclear Gene Expression to Control Differentiation, Stem Cell Function, and Immune Response in Leukemia. Cancer Discov, 11(5):1052-1066.
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I et al: TCGAbiolinks: an R/Bioconductor package for integrative analysis of TC. Nucleic Acids Res, 44(8):e71.
Clarke C, Madden SF, Doolan P, Aherne ST, Joyce H, O'Driscoll L, Gallagher WM, Hennessy BT, Moriarty M, Crown J et al: Correlating transcriptional networks to breast cancer survival: a larg e-scale coexpression analysis. Carcinogenesis, 34(10):2300-2308.
Zhou Q, Liu X, Lv M, Sun E, Lu X, Lu C: Genes That Predict Poor Prognosis in Breast Cancer via Bioinformatical Analysis. Biomed Res Int, 2021:6649660.
Parker JS, Mullins M, Cheang MCU, Leung S, Voduc D, Vickery T, Davies S, Fauron C, He X, Hu Z et al: Supervised Risk Predictor of Breast Cancer Based on Intrinsic Subtypes. J Clin Oncol, 41(26):4192-4199.
Davis S, Meltzer PS: GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioCo nductor. Bioinformatics, 23(14):1846-1847.
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD: The sva package for removing batch effects and other unwanted variatio n in high-throughput experiments. Bioinformatics, 28(6):882-883.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and m icroarray studies. Nucleic Acids Res, 43(7):e47.
Sun X, Zhang Y, Li H, Zhou Y, Shi S, Chen Z, He X, Zhang H, Li F, Yin J et al: DRESIS: the first comprehensive landscape of drug resistance informati on. Nucleic Acids Res, 51(D1):D1263-D1275.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P: The Molecular Signatures Database (MSigDB) hallmark gene set collectio n. Cell Syst, 1(6):417-425.
Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y et al: The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics, 54:1.30.31-31.30.33.
Lin L, Chen L, Xie Z, Chen J, Li L, Lin A: Identification of NAD<sup>+</sup> Metabolism-Derived Gene Signatures i n Ovarian Cancer Prognosis and Immunotherapy. Front Genet, 13:905238.
Ye Z, Zhang H, Kong F, Lan J, Yi S, Jia W, Zheng S, Guo Y, Zhan X: Comprehensive Analysis of Alteration Landscape and Its Clinical Signif icance of Mitochondrial Energy Metabolism Pathway-Related Genes in Lun g Cancers. Oxid Med Cell Longev, 2021:9259297.
Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-seq data wi th DESeq2. Genome Biol, 15(12):550.
Zhang H, Meltzer P, Davis S: RCircos: an R package for Circos 2D track plots. BMC Bioinformatics, 14:244.
Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP: Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res, 28(11):1747-1756.
Li J, Miao B, Wang S, Dong W, Xu H, Si C, Wang W, Duan S, Lou J, Bao Z et al: Hiplot: a comprehensive and easy-to-use web service for boosting publi cation-ready biomedical data visualization. Brief Bioinform, 23(4):bbac261.
Mi H, Muruganujan A, Ebert D, Huang X, Thomas PD: PANTHER version 14: more genomes, a new PANTHER GO-slim and improvemen ts in enrichment analysis tools. Nucleic Acids Res, 47(D1):D419-D426.
Kanehisa M, Goto S: KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res, 28(1):27-30.
Yu G, Wang L-G, Han Y, He Q-Y: clusterProfiler: an R package for comparing biological themes among ge ne clusters. OMICS, 16(5):284-287.
Lock EF, Dunson DB: Bayesian consensus clustering. Bioinformatics, 29(20):2610-2616.
Wilkerson MD, Hayes DN: ConsensusClusterPlus: a class discovery tool with confidence assessmen ts and item tracking. Bioinformatics, 26(12):1572-1573.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES et al: Gene set enrichment analysis: a knowledge-based approach for interpret ing genome-wide expression profiles. Proc Natl Acad Sci U S A, 102(43):15545-15550.
Hänzelmann S, Castelo R, Guinney J: GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics, 14:7.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA: Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol, 1711:243-259.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA et al: Inferring tumour purity and stromal and immune cell admixture from exp ression data. Nat Commun, 4:2612.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z: Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Rela tionships and Predictors of Response to Checkpoint Blockade. Cell Rep, 18(1):248-262.
Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, Liu XS: Large-scale public data reuse to model immunotherapy response and resi stance. Genome Med, 12(1):21.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B et al: Signatures of T cell dysfunction and exclusion predict cancer immunoth erapy response. Nat Med, 24(10):1550-1558.
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E et al: The cBio cancer genomics portal: an open platform for exploring multid imensional cancer genomics data. Cancer Discov, 2(5):401-404.
de Bruijn I, Kundra R, Mastrogiacomo B, Tran TN, Sikina L, Mazor T, Li X, Ochoa A, Zhao G, Lai B et al: Analysis and Visualization of Longitudinal Genomic and Clinical Data f rom the AACR Project GENIE Biopharma Collaborative in cBioPortal. Cancer Res, 83(23):3861-3867.
Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E et al: Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal, 6(269):pl1.
Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR et al: Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeu tic biomarker discovery in cancer cells. Nucleic Acids Res, 41(Database issue):D955-961.
Geeleher P, Cox N, Huang RS: pRRophetic: an R package for prediction of clinical chemotherapeutic r esponse from tumor gene expression levels. PLoS One, 9(9):e107468.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P et al: STRING v11: protein-protein association networks with increased covera ge, supporting functional discovery in genome-wide experimental datase ts. Nucleic Acids Res, 47(D1):D607-D613.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T: Cytoscape: a software environment for integrated models of biomolecula r interaction networks. Genome Res, 13(11):2498-2504.
Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, Franz M, Grouios C, Kazi F, Lopes CT et al: The GeneMANIA prediction server: biological network integration for ge ne prioritization and predicting gene function. Nucleic Acids Res, 38(Web Server issue):W214-220.
Engebretsen S, Bohlin J: Statistical predictions with glmnet. Clin Epigenetics, 11(1):123.
Park SH, Goo JM, Jo C-H: Receiver operating characteristic (ROC) curve: practical review for ra diologists. Korean J Radiol, 5(1):11-18.
Rich JT, Neely JG, Paniello RC, Voelker CCJ, Nussenbaum B, Wang EW: A practical guide to understanding Kaplan-Meier curves. Otolaryngol Head Neck Surg, 143(3):331-336.
Tataranni T, Piccoli C: Dichloroacetate (DCA) and Cancer: An Overview towards Clinical Applica tions. Oxid Med Cell Longev, 2019:8201079.
Hui S, Ghergurovich JM, Morscher RJ, Jang C, Teng X, Lu W, Esparza LA, Reya T, Le Z, Yanxiang Guo J et al: Glucose feeds the TCA cycle via circulating lactate. Nature, 551(7678):115-118.
Sainero-Alcolado L, Liaño-Pons J, Ruiz-Pérez MV, Arsenian-Henriksson M: Targeting mitochondrial metabolism for precision medicine in cancer. Cell Death Differ, 29(7):1304-1317.
Sun S, Li H, Chen J, Qian Q: Lactic Acid: No Longer an Inert and End-Product of Glycolysis. Physiology (Bethesda), 32(6):453-463.
Lv H, Lv G, Chen C, Zong Q, Jiang G, Ye D, Cui X, He Y, Xiang W, Han Q et al: NAD(+) Metabolism Maintains Inducible PD-L1 Expression to Drive Tumor Immune Evasion. Cell Metab 2021, 33(1):110-127.e115.
Guo D, Tong Y, Jiang X, Meng Y, Jiang H, Du L, Wu Q, Li S, Luo S, Li M et al: Aerobic glycolysis promotes tumor immune evasion by hexokinase2-mediated phosphorylation of IκBα. Cell Metab 2022, 34(9):1312-1324.e1316.
Li M, Zhang L, Guan T, Huang L, Zhu Y, Wen Y, Ma X, Yang X, Wan R, Chen J et al: Energy stress-activated AMPK phosphorylates Snail1 and suppresses its stability and oncogenic function. Cancer Lett 2024, 595:216987.
Chae HS, Hong ST: Overview of Cancer Metabolism and Signaling Transduction. Int J Mol Sci 2022, 24(1).
Schmid P, Salgado R, Park YH, Muñoz-Couselo E, Kim SB, Sohn J, Im SA, Foukakis T, Kuemmel S, Dent R et al: Pembrolizumab plus chemotherapy as neoadjuvant treatment of high-risk, early-stage triple-negative breast cancer: results from the phase 1b open-label, multicohort KEYNOTE-173 study. Ann Oncol, 31(5):569-581.
Zhang Y, Schmidt-Wolf IGH: Ten-year update of the international registry on cytokine-induced killer cells in cancer immunotherapy. J Cell Physiol 2020, 235(12):9291-9303.
Xu T, Yu S, Zhang J, Wu S: Dysregulated tumor-associated macrophages in carcinogenesis, progressi on and targeted therapy of gynecological and breast cancers. J Hematol Oncol, 14(1):181.
Xiang X, Wang J, Lu D, Xu X: Targeting tumor-associated macrophages to synergize tumor immunotherap y. Signal Transduct Target Ther, 6(1):75.
Kang JH, Zappasodi R: Modulating Treg stability to improve cancer immunotherapy. Trends Cancer, 9(11):911-927.
Shan P, Yang F, Yu J, Wang L, Qu Y, Qiu H, Zhang H, Zhu S: A novel histone deacetylase inhibitor exerts promising anti-breast can cer activity via triggering AIFM1-dependent programmed necrosis. Cancer Commun (Lond), 42(11):1207-1211.
Sica V, Bravo-San Pedro JM, Izzo V, Pol J, Pierredon S, Enot D, Durand S, Bossut N, Chery A, Souquere S et al: Lethal Poisoning of Cancer Cells by Respiratory Chain Inhibition plus Dimethyl α-Ketoglutarate. Cell Rep, 27(3):820-834.e829.
Zhang Y, Yang Y, Liu R, Meng Y, Tian G, Cao Q: Downregulation of microRNA-425-5p suppresses cervical cancer tumorigen esis by targeting AIFM1. Exp Ther Med, 17(5):4032-4038.
Xiang J, Su R, Wu S, Zhou L: Construction of a prognostic signature for serous ovarian cancer based on lactate metabolism-related genes. Front Oncol, 12:967342.

Download PDF

Editorial decision: Major revision
30 Oct, 2024
Reviewers agreed at journal
09 Aug, 2024
Reviewers invited by journal
07 Aug, 2024
Editor assigned by journal
20 Jul, 2024
First submitted to journal
17 Jul, 2024

You are reading this latest preprint version

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

Status:

Version 1

Abstract

Background:

Methods:

Results:

Conclusions:

Figures

Background

Methods and Materials

Data Download

Differentially expressed genes related to drug resistance and mitochondrial energy metabolism

Mutation analysis and copy number variation analysis

Functional and pathway enrichment analysis

Breast Cancer subtype construction

Gene Set Enrichment Analysis (GSEA)

Gene Set Variation Analysis (GSVA)

Immune infiltration analysis to compare the effects of immunotherapy

Immunotherapy and IPS analysis

TMB analysis and TIDE analysis

Drug sensitivity analysis

Protein interaction network

LASSO regression analysis and screening and validation of key genes

Survival analysis

Prognostic clinical correlation analysis

Cell culture and transfection

Western blotting

Colony Formation Assay

Migration and invasion assays

Statistical analysis

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1