Abstract Background Drug resistance is the major cause of the high death rates in Breast cancer (BC), which continues to be the most frequent disease among women. Chemoresistance is significantly mediated by drug efflux transporters, including ATP-binding cassette transporter ABCC1, and glycolytic enzymes, especially lactate dehydrogenase A (LDHA). Improving treatment results requires an understanding of the expression patterns, genetic changes, and prognostic importance of ABCC1 and LDHA. Objective This study aims to elucidate the role of LDHA and ABCC1 in BC prognosis, tumor progression, and treatment resistance using integrated bioinformatics and in-silico approaches. Methodology To investigate the expression and correlation between LDHA and ABCC1, in-silico analysis was done using a range of bioinformatics tools, such as UALCAN, TIMER 2.0, Bc GeneExminer, DISCO, and others, were used. GeneMANIA and STRING databases were used to explore gene–gene and protein–protein interaction networks, while KM Plotter evaluated survival correlations. Functional enrichment and pathway analyses were conducted using Enrichr for Gene Ontology (GO) and KEGG pathways. For therapeutic targeting, structure-based molecular docking was performed using AutoDock Vina, where selected anticancer compounds were docked against LDHA and ABCC1 to identify potential inhibitors. Results Our research indicates that the expression levels of LDHA and ABCC1 are elevated in several malignancies including Breast Cancer. The elevated expression levels of LDHA and ABCC1 significantly correlate with worse overall survival. Expression analysis of ENO1 (Enolase 1) and ESR1 (Estrogen Receptor 1) genes in relation to LDHA and ABCC1 mutations in Breast Cancer samples revealed that higher ENO1 expression is observed in mutant samples, while ESR1 expression is significantly reduced, suggesting an association with altered metabolic and hormonal pathways. Furthermore, NHI-2 and Sulfinpyrazone were found as the potential chemical that targets LDHA and ABCC1 through in-silico studies. Conclusion This study highlights the oncogenic significance of LDHA and ABCC1 in Breast Cancer progression and therapy resistance. The in-silico identification of NHI-2 and sulfinpyrazone as potential inhibitors supports a novel dual-targeting strategy to simultaneously disrupt metabolic and drug efflux pathways, which could enhance therapeutic efficacy and overcome resistance in Breast Cancer. Further experimental validation is warranted to confirm these findings and facilitate clinical translation. Keywords: Breast cancer, Glycolysis, LDHA, ABCC1, Drug resistance, Metabolism Introduction Cancer remains a major global health challenge, accounting for over 19.3 million new cases and approximately 10 million deaths worldwide in 202 [[24]1]. The disease’s increasing global prevalence has increased demand for effective treatments [[25]2, [26]3]. Global Demographic Characteristics predict an increase in cancer incidence in the next decades, with over twenty million new cases expected by 2025 [[27]4]. BC (BC) is the most prevalent form of malignancy among females, accounting for almost 685,000 deaths in 2020 [[28]5–[29]11]. As a highly heterogeneous malignancy, understanding its common molecular mechanisms could aid in identifying targeted therapies for its various subtypes [[30]12, [31]13]. Drug resistance is the main cause of this increased mortality and a significant challenge in the treatment of BC. Drug resistance accounts for up to 90% of cancer-related fatalities, leading to the lack of effectiveness of pharmacological treatment [[32]14–[33]16]. Drug resistance in cancer results from the disease developing a tolerance to medication, with various causes including genetic mutations, epigenetic modifications, conserved but increased drug efflux, and other cellular and molecular pathways contributing to resistance to anticancer medicines [[34]17, [35]18]. Cancer cells exhibit metabolic alterations that sustain rapid ATP production, increased macromolecule biosynthesis, and maintain cellular redox state [[36]19, [37]20]. These alterations are influenced by both intrinsic and extrinsic mechanisms, including oncogenic and tumor suppressor signalling pathways like HIF1, p53, and MYC [[38]21]. The tumor microenvironment, including interaction with surrounding cells and nutrient and oxygen availability variations, also influences metabolism [[39]22, [40]23]. Central carbon metabolism pathways, such as glycolysis and the tricarboxylic acid cycle, are affected. Cancer cells consume glucose and glutamine to meet their metabolic needs [[41]24]. The addiction to specific metabolic pathways has led to the development of novel drugs targeting these metabolic vulnerabilities [[42]25]. This highlights the importance of a rewired metabolism in cancer cells for rapid cell growth and proliferation [[43]26, [44]27]. The question of whether chemoresistance is caused directly by enzymes, metabolites, or transporters, or indirectly by an overall improvement in glycolytic metabolism in cancer cells brought on, by the upregulation of one of these enzymes, is crucial. Both circumstances can be possible: the upregulation of a single enzyme causes an overall increase in glycolytic metabolism, which is the cause of drug resistance [[45]28, [46]29]. The latter occurs because glycolytic enzymes possess nonenzymatic properties which generate drug resistance [[47]30, [48]31]. Various methods have been used to document the direct participation of nonenzymatic processes in drug resistance in tumor cells. A glycolytic enzyme’s post-translational modification was directly responsible for the production of several nonenzymatic activities, including drug resistance [[49]30], and preventing this post-translational alteration inhibited the development of drug resistance [[50]32]. Glycolytic enzyme mutants that have lost their enzymatic function but nonetheless cause chemoresistance can also cause drug resistance [[51]33]. Ultimately, drug resistance brought on by glycolysis can also result from increased efflux of drugs or blockage of drug influx. The activity of transporter proteins in the plasma membrane regulates drug efflux, with the transmembrane proteins that make up the ATP-binding cassette (ABC) transporter family being the most extensively studied. Cellular ATP levels have a significant impact on ABC transporter efflux activity [[52]34], and drug resistance is induced by the ATP produced in tumor cells by glycolytic metabolism. High glycolytic rates cause a distinctive rise in lactate synthesis, which builds up in the extracellular environment and has been linked to immune evasion, invasion, and senescence suppression [[53]35–[54]37]. Lactate levels in the tumor microenvironment can increase by up to forty times and are strongly associated with poor prognosis and distant metastases BC stem cells, which have elevated levels of glycolytic activity and ABCC1 expression, are at risk of metastasis and resistance to treatment. LDHA and ABCC1 play a crucial role in resistance to drugs, metabolic reprogramming, and stress survival [[55]38]. LDHA drives the Warburg effect, creating an acidic environment that promotes invasion and immune evasion. Under hypoxic conditions, both genes increase, allowing tumor cells to adapt to extreme settings. The lactate-rich acidic environment produced by LDHA may lower the effectiveness of medications, increasing ABCC1’s ability to export medications and lessen oxidative stress [[56]39]. Targeting both LDHA and ABCC1 presents a promising strategy to overcome BC aggressiveness and resistance to treatment. In this connection, we selected two drugs namely Methyl 1-hydroxy-6-phenyl-4-(trifluoromethyl)-1 H-indole-2-carboxylate (NHI-2), a strong inhibitor of LDHA, and sulfinpyrazone a sulfoxide and a member of pyrazolidines to target ABCC1. NHI-2 inhibits LDHA, disrupting the glycolytic pathway that cancer cells rely on for energy production, making it a potential candidate for cancer research. Sulfinpyrazone is a platelet inhibitory and uricosuric agent used to inhibit thrombotic and embolic processes and manage chronic gout phases. In the current work, we aim to investigate the correlation between LDHA and ABCC1 expression levels and their impact on BC prognosis, drug resistance, and tumor progression using bioinformatics tools. The expression and prognostic significance of the metabolic enzyme LDHA and the MDR-associated protein ABCC1 in BC were explored using bioinformatics tools. Furthermore, the two genes namely, ENO1 (Enolase 1) and ESR1 (Estrogen Receptor 1) play crucial roles in cancer metabolism and drug resistance. Their overexpression is linked to therapeutic failure, making them important targets for studying the molecular mechanisms underlying chemoresistance in breast cancer. Given their central roles in metabolic reprogramming and chemoresistance, the investigation of these genes, alongside LDHA and ABCC1, is essential for elucidating novel therapeutic strategies targeting the metabolic and efflux pathways critical for overcoming drug resistance in breast cancer.The novelty of this study lies in its integrative computational approach to identify and evaluate the dual targeting of LDHA and ABCC1—two key proteins implicated in metabolic reprogramming and drug resistance in BC. While these proteins have been individually studied in various cancers, our work is among the first to explore their combined prognostic impact and potential co-targeting strategy. Additionally, we repurpose sulfinpyrazone, traditionally used as a uricosuric agent, as a novel inhibitor of ABCC1 in the context of BC, which to our knowledge has not been previously reported. Together with NHI-2, a potent LDHA inhibitor, this dual-targeting approach provides a unique strategy to simultaneously disrupt glycolytic metabolism and drug efflux mechanisms—offering a promising therapeutic avenue for overcoming chemoresistance. Materials and methods GTEX database The GTEx database was utilized to analyze bulk tissue gene expression data for ABCC1 and LDHA [[57]40]. The expression levels were quantified as Transcripts Per Million (TPM) and visualized using violin plots to illustrate distribution across multiple tissue types. Data were normalized to account for technical variability, and statistical analyses were performed to assess significant differences in gene expression. The insights from this analysis help identify tissue-specific patterns that may influence resistance mechanisms and therapeutic responses. LDHA and ABCC1’s expression pattern in PAN cancer An innovative web-based tool called TIMER 2.0 ([58]http://timer.cistrome.org/) was developed specially to examine immune infiltrates and expression patterns across a variety of cancer types. TIMER, which was developed to improve our knowledge of the tumor microenvironment, enables researchers to look at the connections between the makeup of immune cells and cancer traits [[59]41]. This platform is an effective tool for integrative cancer research since it makes use of data from TCGA. TIMER 2.0 sheds light on the expression of LDHA and ABCC1, in a variety of cancer types. Users can investigate the variations in LDHA and ABCC1 expression across tumor types, examine its relationship to immune infiltrates, and evaluate its possible significance in tumor progression and patient outcome by utilizing TCGA data. To compare gene expression between tumor and normal breast tissues, TCGA-BRCA tumor samples via TIMER2.0 portals used data harmonized through the Toil recompute pipeline, which minimizes batch effects arising from different sequencing protocols. UALCAN One useful tool for analysing OMICS data about cancer is the UALCAN online database ([60]https://ualcan.path.uab.edu/). Its main goal is to make publicly available cancer OMICS data easier to use, which will help identify possible biomarkers and allow candidate genes to be validated in silico [[61]42]. The PERL-CGI programming language was used to create the software, which also includes sophisticated JavaScript-generated visuals. Information about patient survival rates and gene expression profiles for families expressing proteins, long non-coding RNAs, and microRNAs is available to users in the form of graphs and statistics (lncRNAs). To compare gene expression between tumor and normal breast tissues, we utilized TCGA-BRCA tumor samples via the UALCAN portal. This platform use data harmonized through the Toil recompute pipeline, which minimizes batch effects arising from different sequencing protocols. Using the TCGA breast invasive cancer dataset, this study explicitly examined the expression of LDHA and ABCC1 in BC and its various subclasses [[62]43]. bc-GenEXMiner The web-based programme bc-GenExMiner v4.5 ([63]http://bcgenex.centregauducheau.fr/BC-GEM) was developed especially for the analysis of transcriptome data related to BC. With the help of this platform, researchers may conduct thorough analyses of molecular characteristics and their clinical consequences thanks to access to a huge library of gene expression profiles obtained from several BC studies [[64]44]. To examine the relationship between LDHA, ABCC1 expression, and other clinical features of BC individuals, including hormone receptor status, nodal status, NPI status, and SBR grade, bc-GenEXMiner v5.0, was used for annotated BC transcriptome data [[65]44]. GeneMANIA To identify novel genes associated with a collection of input genes GeneMANIA ([66]https://genemania.org/) uses a large library of functional association data [[67]45]. Examples of association data include pathways, co-expression, co-localization, and similarities between genetic and protein domains. We used the online tool GeneMANIA to forecast the role and relationship of LDHA and ABCC1 in the gene-gene interaction study. A flexible and user-friendly interface for figuring out the function of gene sets and genes of interest is provided by this free resource. The database was used to anticipate the function of ABCC1 and LDHA as well as their association. PPI investigation A protein-protein interaction (PPI) network of LDHA and ABCC1 with a 0.7 confidence level was built using the online database STRING ([68]http://stringdb.org). To establish and study functional relationships between proteins, the online biological resource STRING was developed [[69]46]. The PPI web was further analysed and visualised using the Cytoscape tool (version) [[70]47]. The MCODE (Molecular Complex Detection) plug-in for the Cytoscape program was used to investigate the important PPI network modules. The Cytohubba plugin option was used to analyse the top 10 HUB nodes on the PPI web [[71]48]. DISCO DISCO represents the Deeply Integrated Human Single-Cell Omics data ([72]https://www.immunesinglecell.org/) which has an annotated atlas containing information about various tissues like breast [[73]49]. Various interactive plots can be generated using DISCO thus presenting an easy-to-use online tool for researchers. The users can further navigate between various parameters available within DISCO. LDHA and ABCC1 analysis in single-cell sequencing database (TISCH) The TISCH was first processed to observe major modifications in LDHA and ABCC1 expression across various datasets in tissues of BC [[74]50]. To determine the connection between tumor development and LDHA and ABCC1 expression patterns, we analyzed the expression of these proteins in primary and metastatic BC datasets. Kaplan-Meier plotter The KM plotter ([75]https://kmplot.com) an online portal, offers information on BC patient survival and gene expression [[76]51]. The mRNA gene chip BC dataset of the KM Plotter contains OS data for 1879 patients. The effect of LDHA and ABCC1 gene expression on the OS of the assigned groups was investigated using KM-survival curves, and the log-rank P-value and hazard ratio ranges were calculated using the KM plotter. Gene ontology (GO) and pathway enrichment of LDHA and ABCC1 One of the most significant bioinformatics portal is Enrichr. This database’s primary goal is to connect the characteristics of genes and gene products across all species. It provides us with a collection of thousands of systematically arranged terminologies for biological activities, cellular components, and molecular processes. It also predicts gene annotations with a P-value of 0.05 and curates them. The GO annotation inquiry was analyzed by Enricher ([77]http://amp.pharm.mssm.edu/Enrichr), a comprehensive online resource that collects biological data for upcoming biological discoveries. Enricher was also used to do further KEGG pathway analysis [[78]52]. TNM plots TNM Plot ([79]https://tnmplot.com/analysis/) is a bioinformatic visualization tool used to examine the expression of genes across normal (N), tumor (T), and metastatic (M) tissue samples, typically from large datasets like TCGA [[80]53]. The expression data for the gene of interest is extracted from these datasets, which include samples from normal tissues, primary tumors, and metastatic lesions. The plots are generated in the form of bar graphs or boxplots, to display differential gene expression across these groups. The significance of expression differences is assessed using statistical tests including the paired t-test, ANOVA, and Wilcoxon rank-sum test. This approach highlights possible indicators and treatment targets by shedding information on how gene expression varies throughout tumor growth and metastasis. muTarget Mutarget is an online database that integrates several kinds of genomic, transcriptomic, and mutational data to make it easier to identify possible therapeutic targets. Mutarget was used in this study to examine mutation data for genes of interest, such as ENO1, ESR1, LDHA, and ABCC1, across several cancer types, with a particular emphasis on BC. In order to get mutation frequencies, gene expression profiles, and related clinical data, the methodology comprised database queries. To determine the importance of the mutations, statistical analysis was done to compare the target gene expression levels in mutant and wild-type samples. In order to identify therapeutic targets and biomarkers for drug repurposing, the data from Mutarget was also utilized to investigate possible correlations between gene mutations and pathways linked to cancer. Human protein atlas An important portal called the Human Protein Atlas ([81]https://www.proteinatlas.org/) aims to map the patterns of protein expression in healthy human cells, tissues, and cancer types. The project intends to offer a thorough resource for comprehending the function of proteins in human biology and illness by examining 11,200 distinct proteins, which account for about 50% of all human protein-coding genes. The study uses a variety of cutting-edge methods, such as mass spectrometry, transcriptomics, and immunohistochemistry, to methodically assess protein quantity and localization [[82]54]. This gives scientists a better understanding of the functional roles of particular proteins by allowing them to see where they are expressed in various tissues, including the liver, brain, and different forms of cancer. The Human Protein Atlas accomplishes several goals by assembling this information into an easily navigable database. Mitelman database The Mitelman Database of Chromosome Aberrations and Gene Fusions in Cancer is a comprehensive resource catalogs chromosomal alterations associated with various neoplastic diseases. It provides valuable insights into genomic rearrangements, gene fusions, and copy number variations (CNVs) linked to cancer progression and drug resistance. Data from the Mitelman Database was analyzed using a Circos plot to visualize chromosomal alterations, focusing on the ABCC1 gene (chr16, 16p13.11). Copy number variations (CNVs) and translocations were identified and correlated with drug resistance mechanisms. Docking of LDHA+NHI-2 and ABCC1+sulfinpyrazone Target selection and its preparation The target proteins LDHA and ABCC1 were obtained in PDB format from the Protein Data Bank using IDs 4OJN and 8VT4. The ‘Prepare protein’ tool in Discovery Studio was utilized to add missing main- and sidechain atoms, standardize the atom’s name, remove disordered alternative conformations, and simplify the structure for ABCC1 and LDHA. Additionally, H2O molecules and heteroatoms were removed. Molecular docking The goal of this study is to provide insight into the binding propensities. of target LDHA having PDB ID (4OJN) with ligand CID_ 51,355,147(NHI-2) and target ABCC1 having PDB ID (8VT4) with ligand CID_5342 (Sulfinpyrazone). Autodock Vina 4.2.6 was used to conduct docking studies on the molecule in question to have a correct docking to forecast the effectiveness of ligand binding. Following the combination of non-polar hydrogens, the receptor and target molecule are both stored in pdbqt format. Ligands are necessary to design grid boxes with particular dimensions and 0.3 Å spacing. The LGA method was used to guide the protein-ligand complex’s docking studies. A population size of 500, a maximum generational number of 27, 2,500,000 evaluations, 50 solutions, and all other parameters were left at their default settings in each of the three replicates of the molecular docking studies. After the ligand and receptor had completed docking, re-clustering with clustering tolerances of 2.0 Å was used to find the optimal cluster with the lowest energy score and the largest populations. The RMSD clustering maps were so developed. Molecular dynamics simulation (MDS) The stability of the complexes produced by 4OJN and CID_51355147 (LDHA and NHI-2) and 8VT4 and CID_5342 (ABCC1 and Sulfinpyrazone) was investigated using MD simulations. The Desmond 2020.1 was used to do the simulations [[83]55]. At a temperature of 27 °C, separate simulations were run for every complex, including the control. The OPLS-2005 force fields, which have been detailed in earlier research, were used to depict the molecular interactions [[84]56–[85]58]. SPC on molecules is employed in an explicit solvent model [[86]59]. A 10 Å x 10 Å x 10 Å periodic boundary salvation box encloses the simulated systems CID_51355147 and CID_5342. The system’s charge was neutralized by adding sodium ions (Na+), and the physiological environment was simulated by adding NaCl solutions. An NVT ensemble simulation CID_51355147 and CID_5342 or 10 ns was used for the equilibration phase to give the protein-ligand complexes time to stabilize. The purpose of this procedure was to guarantee that the system stabilized [[87]60]. After that, a brief period of equilibration and minimization was conducted for 12 ns utilizing an NPT ensemble. The Nose-Hoover chain coupling approach [[88]61] was used to set up the NPT ensemble CID_51355147 and CID_5342, which had a 1.0 ps relaxation time and varied temperature. In every simulation, the pressure was kept constant at 1 bar. The simulations were conducted at a time step of 2 ps. With a relaxation duration of 2 ps, the Martyna-Tuckerman-Klein chain coupling technique [[89]61] was used to regulate the pressure. With a constant Coulomb interaction radius of 9 Å, long-range electrostatic interactions were computed using the particle mesh Edwld methods [[90]62]. To determine the bonded forces for each trajectory, the RESPA integrator used a time step of two. A final 100 ns manufacturing run was then developed for simulation. The RMSD, the Rg, the RMSF, and the number of hydrogen bonds were among the characteristics that were computed to evaluate MD simulations stability. Throughout the simulations, these metrics were used to keep an eye on the system’s stability. Binding free energy (ΔGbind) analysis To determine the ΔGbind of the ligand-protein complexes, the molecular mechanics combined with the generalized Born surface area (MM-GBSA) approach CID_51355147 and CID_5342 is employed. The calculation of MM-GBSA ΔGbind formed using the Python script “thermal mmgbsa.py” on the simulation trajectory’s last 50 frames, with a sampling size of 1 step. The Prime MM-GBSA ΔGbind (in kcal/mol) is estimated by applying the principle of additivity. This involves the summation of individual energy components, including coulombic interactions, covalent bonds, hydrogen bonds, van der Waals forces, self-contacts, lipophilic interactions, and solvation energies of both the protein and ligand. The equation used to calculate ΔGbind is as follows: graphic file with name d33e491.gif The contribution of each CID_51355147 and CID_5342 ic energy component to the total binding free energy (ΔGbind) is represented by a term in this equation. The MM-GBSA method takes into account several energy parameters, enabling a thorough assessment of the binding affinity. Results Tissue-specific expression patterns of ABCC1 and LDHA genes from GTEx database The GTEx database analysis revealed that ABCC1 and LDHA genes exhibit differential expression across various human tissues. ABCC1 shows relatively high expression in the adrenal gland and whole blood, whereas LDHA demonstrates significantly elevated expression in skeletal muscle and heart tissue. These expression patterns indicate tissue-specific regulation and potential functional roles in metabolism and drug resistance, particularly highlighting the metabolic adaptation and efflux mechanisms contributing to multidrug resistance (MDR) in cancer (Fig. [91]1). Fig. 1. [92]Fig. 1 [93]Open in a new tab Violin plots showing tissue-specific expression of ABCC1 (A) and LDHA (B) genes from the GTEx database. Expression levels (TPM) reveal high ABCC1 expression in the adrenal gland and whole blood, and elevated LDHA expression in skeletal muscle and heart tissue LDHA and ABCC1 is highly expressed in many cancers The expression levels of LDHA and ABCC1 in various cancer types were determined using the TIMER 2.0 database. The TIMER 2.0 study shows that both ABCC1 and LDHA are highly expressed in several cancers. Many cancers, such as BC, CESC, COAD, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LUAD, LUSC, PRAD, READ, STAD, and UCEC, have higher levels of LDHA, according to box plots (Fig. [94]2A). Furthermore, many cancers have increased expression of ABCC1, such as BC, CESC, CHOL, COAD, ESCA, HNSC, KICH, KIRC, KIRP, LIHC, LUSC, PRAD, READ, SKCM, STAD, THCA, and UCEC (Fig. [95]2B). Fig. 2. [96]Fig. 2 [97]Open in a new tab A The expression of LDHA and ABCC1 in a variety of malignancies utilizing the TIMER Database’s box plot approach. LDHA. B ABCC1’s expression is significantly increased in the majority of malignancies, including BC LDHA and ABCC1 are highly expressed in BC patients The UALCAN database was used to analyze LDHA and ABCC1 expression in connection to sample type and BC subtypes (Fig. [98]3). The analysis revealed that LDHA and ABCC1 expression levels were greater in primary breast tumors than normal patients (Figs. [99]3A and [100]3C). Compared to TNBC or luminal subtypes, LDHA expression was higher in individuals with HER2-positive BC (Fig. [101]3B). Additionally, ABCC1 expression was augmented in TNBC followed by luminal and HER2-positive BC (Fig. [102]3D). Fig. 3. [103]Fig. 3 [104]Open in a new tab BC patients’ expression profiles for LDHA and ABCC1. The expression was higher in individuals with primary tumors A and in females C. Patients with HER2-positive BC exhibited high expression of LDHA B and ABCC1 expression was elevated in TNBC followed by luminal and HER2-positive BC patients D Overexpression of LDHA and ABCC1 positively correlates with ER negative and PR negative status Using the bc-GenEXMiner portal, we analyzed the relationship between LDHA, ABCC1, and clinicopathological characteristics of BC individuals. LDHA (Fig. [105]4) and ABCC1(Fig. [106]5) were expressed at higher levels in ER-negative and PR-negative BC patients indicating a positive correlation. Additionally, HER2-positive status exhibited higher levels of LDHA and ABCC1 expression than HER2-negative patients. Also, LDHA was significantly expressed in the IDC, SBR3, NPI3, and TNBC grades. In the case of ABCC1, the nodal negative, IDC, SBR3, NPI2, and TNBC grades showed significant expression. Fig. 4. [107]Fig. 4 [108]Open in a new tab The expression of LDHA is related to many clinicopathological characteristics. Elevated expression of LDHA is linked to the following: ER-ve, PR-ve, HER2 + ve, IDC, SBR3, NPI3, and TNBC status Fig. 5. [109]Fig. 5 [110]Open in a new tab The expression of ABCC1 is related to many clinicopathological characteristics. Elevated expression of ABCC1 is linked to the following: Nodal negative, ER-ve, PR-ve, HER2 + ve, IDC, SBR3, NPI2, and TNBC status A strong correlation was seen between LDHA and ABCC1 We looked into the LDHA and ABCC1 gene-gene interaction network using the GeneMANIA portal (Fig. [111]6). As can be shown, LDHA and ABCC1 had robust interactions with LDHB and LMBRD1, followed by TMEM33, ABCCB, LDHC, ABCC3, TCN2, PKM, MDH1, and several other genes (Fig. [112]6A). The relationships within the gene-gene interaction network were based on genetic interaction pathways, co-localization, co-expression, and physical linkages. With an R-value of 0.29, the analysis of the connection between LDHA and ABCC1 using Gepia2 revealed a substantial association between LDHA and ABCC1 in BC patients (Fig. [113]6B). Fig. 6. [114]Fig. 6 [115]Open in a new tab GGI network of LDHA and ABCC1 using GeneMANIA. A LDHA and ABCC1 depicted a high correlation with LDHB, LMBRD1 which was followed by TMEM33, ABCCB, LDHC, ABCC3, TCN2, PKM, MDH1, and many other genes. B With an R-value of 0.29, the Gepia2 portal depicted a substantial relationship between LDHA and ABCC1 in BC PPI of LDHA and ABBC1 Using the STRING portal, PPI was constructed by joining co-expressed 22 genes (nodes) with 58 proteins (edges) for the metabolic enzymes LDHA and ABCC1 protein (Fig. [116]7). The network (protein-protein interaction) had an average local clustering coefficient of 0.754, an average node degree of 5.27, a PPI enrichment p-value of 3.37e-11, and an anticipated number of edges of 21 for the metabolic enzyme LDHA and the ABCC1 protein (Fig. [117]7A). Furthermore, as shown in (Fig. [118]7B), Cytohubba was used to determine the top 10 web-based hub genes based on degree score. KCNJ11, LDHA, CTH, GOT1, LDHC, TST, LDHAL6A, LDHAL6B, MPST, and GOT1L1 were the top 10 genes of the network (Fig. [119]7B). The most important part of the PPI network was extracted by analysis of the Cytoscape program’s MCODE plug-in. Fig. 7. [120]Fig. 7 [121]Open in a new tab PPI of LDHA and ABCC1 using string portal. Using Cytohubba, the top ten hub genes of the network based on degree score were identified. The top 10 genes for LDHA and ABCC1 of the network included KCNJ11, LDHA, CTH, GOT1, LDHC, TST, LDHAL6A, LDHAL6B, MPST, and GOT1L1 Expression of LDHA and ABCC1 in different types of breast cells We explored the DISCO database to find out the expression levels of LDHA and ABCC1 in various breast cells. The results depicted that LDHA and ABCC1 expression is enriched in several different breast tissues/cells (Fig. [122]8). LDHA was highly enriched in KRT6B mammary basal cells, followed by CREB, MT1A, vascular smooth muscle cells, CXCL, pericyte, SAA2 mammary luminal progenitor cells, and many other cells (Fig. [123]8A). The ABCC1 expression was greatly expressed in Treg cells, ILC, CD4T cells, BNC2, ZFPM2, fibroblast, and other cells (Fig. [124]8B). Fig. 8. [125]Fig. 8 [126]Open in a new tab Analysing the expression level of A LDHA and B ABCC1 through the DISCO database Connection of LDHA and ABCC1 with tumor stroma We examined the expression profiles of LDHA and ABCC1 in the TME-focused TISCH database. Expression pattern analysis was performed on several BC datasets that included information on primary and metastatic tumors (Fig. [127]9). According to the data, metastatic BC had much-raised levels of LDHA and ABCC1 than primary BC. Additionally, there were significant differences in the LDHA and ABCC1 expression patterns across the metastatic cancer cell population. B cells, CD4T, and CD8T cells displayed greater levels of LDHA and ABCC1 expression in BC individuals with metastatic cancer (Fig. [128]9B), while epithelial and fibroblast cells did so in BC individuals with primary carcinoma (Fig. [129]9A). Fig. 9. [130]Fig. 9 [131]Open in a new tab LDHA and ABCC1 expression pattern in BC using single-cell RNA-seq database TISCH. A: Primary BC. B: Metastatic BC Overexpression of LDHA and ABCC1 is associated with a worse prognosis in BC individuals To assess the prognostic significance of LDHA and ABCC1 in BC, we utilized the Kaplan–Meier Plotter tool, which correlates gene expression levels with overall survival using publicly available clinical datasets. Figure [132]10A shows that patients with high LDHA expression have significantly poorer overall survival compared to those with low expression (HR = 1.28, 95% CI: 1.06–1.55, p = 0.0094). Similarly, Fig. [133]10B demonstrates that high ABCC1 expression is also associated with reduced survival probability (HR = 1.25, 95% CI: 1.03–1.5, p = 0.023). These findings suggest that overexpression of LDHA and ABCC1 contributes to aggressive tumor behaviour and may serve as negative prognostic biomarkers in BC. The strong correlation between high expression and poor outcome further supports the therapeutic relevance of targeting these genes in future clinical strategies. Fig. 10. [134]Fig. 10 [135]Open in a new tab The OS revealed that increased LDHA (A) and ABCC1 (B) expression correlates with low survival probability Analysis of gene ontology and KEGG The Enrichr database was utilized for GO and KEGG investigation of LDHA and ABCC1 (Fig. [136]11). Among the biological processes, LDHA and ABCC1 were found to be associated with xenobiotic transport across the blood-brain barrier, transepithelial transport, lipid and phospholipid translocation, glycolytic process, carboxylic acid transmembrane transport, carbohydrate catabolic process, pyruvate metabolic process, and carboxylic acid transport (Fig. [137]11A). Among the cellular compartments, LDHA and ABCC1 were more abundant in the basolateral plasma membrane followed by the nucleus and intracellular membrane-bound organelles (Fig. [138]11B). Among molecular functions, LDHA and ABCC1 were involved in ABC-type xenobiotic transporter activity, carboxylic acid transmembrane transporter activity, oxidoreductase activity, activity on the CH-OH group of donors, and acts as NAD or NADP acceptor and cadherin binding (Fig. [139]11C). The KEGG pathway study investigated that LDHA and ABCC1 were related to vitamin digestion and absorption, propanoate metabolism, ABC transporters, pyruvate metabolism, cysteine and methionine metabolism, glycolysis/gluconeogenesis, central carbon metabolism in cancer, glucagon signalling pathway, HIF-1 signalling pathway, and sphingolipid signalling pathway (Fig. [140]11D). Fig. 11. [141]Fig. 11 [142]Open in a new tab GO and KEGG pathway enrichment analysis of LDHA and ABCC1 using the Enrichr portal. A Biological processes, B Cellular components, C Molecular functions, and D KEGG pathways showing key terms associated with LDHA and ABCC1 in BC Elevated expression of LDHA and ABCC1 correlates with BC progression and metastasis LDHA and ABCC1 were significantly upregulated in BC tissues relative to normal tissues, according to the expression analysis, with metastatic samples showing the highest expression. Both genes’ log2-transformed expression levels increased gradually from normal to tumor to metastatic settings, suggesting that they may be involved in the development and spread of tumors. This supports the idea that LDHA and ABCC1 are possible therapeutic targets for BC since higher expression of these proteins is associated with more aggressive disease stages (Fig. [143]12). Fig. 12. [144]Fig. 12 [145]Open in a new tab Expression analysis of LDHA and ABCC1 shows that upward trend in expression levels as the condition progresses from Normal to Tumor to Metastatic, indicating their potential involvement in tumor progression and metastasis Association of LDHA and ABCC1 mutations with ENO1 and ESR1 expression in BC ENO1 and ESR1 are examined about LDHA and ABCC1 because of their roles in drug resistance and cancer metabolism. Similar to ENO1, LDHA enhances glycolysis in cancer cells, which aids in tumor growth and makes it a possible target for treatment. In malignancies that are ESR1-positive, such as BC, ABCC1, a drug transporter, may pump out therapeutic drugs, decreasing their efficacy. This is linked to chemotherapy resistance. ENO1 expression is significantly higher in samples with LDHA and ABCC1 mutations than in wild-type samples (p = 1.63e-03), according to the analysis, although ESR1 expression is significantly lower in mutant samples (p = 1.90e-03). These findings underline their influence on metabolic and hormonal pathways in BC and point to a possible connection between LDHA and ABCC1 mutations (Fig. [146]13). Fig. 13. [147]Fig. 13 [148]Open in a new tab Expression analysis of ENO1 and ESR1 in relation to LDHA and ABCC1 mutations in BC samples. A shows that ENO1 expression is significantly higher in samples with LDHA and ABCC1 mutations compared to wild-type (p = 1.63e-03), suggesting enhanced glycolytic activity in mutant cases. In contrast, B reveals that ESR1 expression is markedly lower in the mutant group (p = 1.90e-03), indicating a potential suppression of estrogen receptor-mediated signaling Immunohistochemical analysis of ABCC1 and LDHA expression in BC Human Protein Atlas results indicate dissimilar expression patterns of immunohistochemical staining between ABCC1 and LDHA in BC tissues. The ABCC1 protein (Fig. [149]14A) shows strong cytoplasmic and membranous staining indicating high expression levels in BC cells. It seems that the protein may play a critical role in drug resistance through the active efflux of chemotherapeutic agents outside the cells. While the LDHA protein (Fig. [150]14B) indicates the weak to moderate cytoplasmic staining, conveying lower expression levels compared to ABCC1. The difference in expression pattern between these two proteins suggests that ABCC1 is probably involved in chemotherapy resistance and that LDHA possibly contributes to metabolic reprogramming in BC cells. All these findings boost the prediction that targeting ABCC1 and LDHA would enhance treatment efficacy in drug-resistant BC. Fig. 14. [151]Fig. 14 [152]Open in a new tab Immunohistochemical staining of BC tissues showing A strong cytoplasmic and membranous expression of ABCC1 and B weak to moderate cytoplasmic expression of LDHA, suggesting ABCC1’s role in chemoresistance and LDHA’s involvement in metabolic reprogramming Genomic amplification of ABCC1 and LDHA and its role in drug resistance The analysis of data from the Mitelman Database revealed a notable amplification of the ABCC1 gene (dark blue bars) on chromosome 16 (16p13.11), indicating an increased copy number. This amplification is commonly observed in drug-resistant cancer samples. It is associated with enhanced drug efflux capacity, as ABCC1 encodes a membrane-bound transporter protein involved in pumping chemotherapeutic agents out of cells. Consequently, this overexpression significantly reduces intracellular drug accumulation, leading to multidrug resistance (MDR). Importantly, no significant homozygous deletions or loss of function mutations were observed at the ABCC1 locus, emphasizing that gene amplification and overexpression are the primary resistance mechanisms rather than gene loss or inactivation. Also, the analysis of data from the Mitelman Database revealed a notable amplification of the LDHA gene (dark blue bars), indicating an increased copy number on its respective chromosome (chromosome 11, 11p15.4). This amplification is commonly associated with metabolic reprogramming in cancer cells, leading to enhanced aerobic glycolysis (Warburg effect). The overexpression of LDHA promotes lactate production and acidification of the tumor microenvironment, which aids in tumor progression and drug resistance. No significant homozygous deletions or loss of function mutations were observed, suggesting that LDHA overexpression primarily contributes to resistance by shifting cellular metabolism towards glycolysis rather than being lost or inactivated (Fig. [153]15). Fig. 15. [154]Fig. 15 [155]Open in a new tab Circos plot showing ABCC1 gene amplification (dark blue bars) on chromosome 16 and LDHA gene amplification on chromosome 12, indicating increased copy numbers associated with multidrug resistance (MDR) and metabolic adaptation, respectively. Inner lines represent interchromosomal translocations Molecular docking showed high binding Molecular docking The molecular docking analysis reveals that the ligand binds to the protein target NHI-2 + LDHA (CID_51355147 + LDHA) and Sulfinpyrazone + ABCC1(CID_5342 + ABCC1) with a binding energy of -7.2 kcal/mol (Fig. [156]16A) and − 6.8 kcal/mol (Fig. [157]16B), indicating a strong and stable interaction. In case of NHI-2 + LDHA, and (CID_5342 + ABCC1) the ligand is positioned within a well-defined binding pocket of the protein, as visualized in the 3D representation. For NHI-2 + LDHA the interaction network is stabilized by multiple interactions including van der Waals forces (with residues such as LEU A:173 and GLN A:233), contributing to general stabilization of the ligand within the pocket. Hydrogen bonding interactions play a critical role, with conventional hydrogen bonds formed with residues like GLU A:240 and SER A:237, and carbon hydrogen bonds with residues such as ARG A:169. Pi-pi stacking interactions (with TYR A:172) and Pi-alkyl interactions (with VAL A:241, ALA A:168, and ALA A:251) are present, highlighting the hydrophobic contributions to the binding. These interactions, particularly the aromatic stacking and alkyl interactions, suggest a strong affinity of the ligand for hydrophobic regions within the binding site. The collective participation of polar residues like GLN A:233 and GLU A:240, alongside hydrophobic residues such as ALA and VAL, stabilizes the ligand in the pocket. Parallel to NHI-2 + LDHA, for Sulfinpyrazone + ABCC1 the Van der Waals forces (with residues such as PHE A:1127, GLN A:1130, ILE A:1119 and LEU A:1098) contribute to overall stability, while hydrogen bonds (with CYS A:1105) play a critical role in specificity. Additionally, Pi-sigma interactions (with LEU A:1247) and Pi-alkyl interactions (with VAL A:1122 and VAL A:1206) further enhance binding through hydrophobic interactions. The key residues involved include hydrophobic amino acids like LEU and VAL, as well as polar residues like TYR A:1126 and GLN A:1130. This docking result highlights the ligand’s favourable fit within the protein’s active site, supported by a combination of hydrophobic, polar, and hydrogen bonding interactions, suggesting its potential efficacy as a modulator or inhibitor of the protein’s activity. Fig. 16. [158]Fig. 16 [159]Open in a new tab A The frequency of populations at the 2.0 tolerance threshold is displayed by the docked posture of LDHA with the ligand NHI-2. The protein’s surface view shows a 2D interaction plot of the ligands’ binding pocket in the corresponding proteins, as well as a deep core of the binding pocket that holds the ligands. B ABCC1 docked posture with ligand Sulfinpyrazone showing population frequency at tolerance level 2.0. The protein’s surface view shows a 2D interaction plot of the ligands’ binding pocket in the corresponding proteins, as well as a deep core of the binding pocket that holds the ligands Molecular dynamics simulation (MDS) Studies using molecular dynamics and simulation, often known as MD, were carried out with the goal of determining the stability and convergence of the LDHA + NHI-2 (4OJN + CID_ 51355147) and ABCC1 + Sulfinpyrazone (8VT4 + CID_5342) system. When comparing the RMSD measurements, the simulation lasting 100 ns showed a stable conformation. While the RMSD of the Cα-backbone of the LDHA showed a deviation of 1.6 Å, the RMSD of the ligand bound to protein showed a deviation of 1.2 Å (Fig. [160]17A). On the other hand, the RMSD Cα-backbone of ABCC1 + Sulfinpyrazone (8VT4 + CID_5342) exhibited more stable 3.0 Å. The RMSD of ligand Sulfinpyrazone also exhibited 3.0 Å (Fig. [161]17B). A stable RMSD plot during simulation signifies good convergence and stable conformations. Thus, it may be inferred that because of the ligand’s greater affinity, proteins and ligands are rather stable in complex. The plot for RMSF of LDHA in NHI-2 bound complex displayed large spikes of fluctuation in protein observed at 130–150 residue positions (Fig. [162]17C) However, there are significant spikes and fluctuations in the ABCC1 protein bound to Sulfinpyrazone at residues 5–45, 165, 250, 410, 498, 598, and 630. This might be because the residues are more flexible. (Fig. [163]17D). Throughout the 100-ns simulation, the majority of the residues showed minimal fluctuation, demonstrating the rigid amino acid conformations. Higher fluctuating residues provide conformational flexibility might be due to loop formation. Since more flexibility enables the protein to become functionally active, fewer fluctuations signify less active protein. The protein’s compactness is gauged by its Rg. In this study, ligand-bound protein Cα-backbone displayed a stable radius of gyration (Rg) of LDHA from 23.55 to 23.50 Å (Fig. [164]17E) while ABCC1 Cα-backbone bound to Sulfinpyrazone displayed a lowering of Rg from 16.9 to 16.54 Å (Fig. [165]17F). A very compact orientation of the protein in a ligand-bound state is indicated by stable gyration (Rg) and a decrease of Rg. The amount of hydrogen bonds that exist between the ligand and protein indicates how well the complex interacts and how stable it is. Between NHI-2 and LDHA protein, hydrogen bonds displayed constant throughout the simulation time and 3 numbers on average (Fig. [166]17G). (Fig. [167]17H) shows that there were a few hydrogen bonds between ABCC1 and Sulfinpyrazone, with an average of up to 100 ns. SASA in both ligand-bound and unbound states showed identical patterns after Rg analysis. SASA exhibited that the plot of the unbound state of ligand to the receptor has a higher surface area (red), and receptor-bound ligand exhibited a lowering of peak (black) (Fig. [168]17I) for LDHA + NHI-2. Furthermore, Fig. [169]17J shows that in every instance, the unbound form of the receptor protein ABCC1 had a large surface area that was solvent accessible. The SASA value lowered as compared to the unbound state while bound with ligands (Fig. [170]17J). A more firmly bound conformation is indicated by a greater lowering of the peak. Rg’s overall analysis shows that ligand binding forces the corresponding proteins to condense. Fig. 17. [171]Fig. 17 [172]Open in a new tab Molecular dynamics simulation analysis of 100 ns trajectories of A Cα backbone RMSD of LDHA. Molecular dynamics (MD) simulation analysis over 100 ns for LDHA–NHI-2 and ABCC1–sulfinpyrazone complexes. (A) RMSD of the Cα backbone for LDHA (black) and its ligand NHI-2 (red), indicating complex stability. B RMSD of the Cα backbone for ABCC1 and sulfinpyrazone, reflecting structural convergence. C, D RMSF plots showing residue-level fluctuations for LDHA (C) and ABCC1 (D) in their respective complexes. E, F Radius of gyration (Rg) profiles for LDHA (E) and ABCC1 (F), assessing compactness and structural integrity. (G–H) Time-dependent hydrogen bond formation between LDHA and NHI-2 (G), and ABCC1 and sulfinpyrazone (H), indicating interaction stability. I, J Solvent-accessible surface area (SASA) comparison between ligand-bound (black) and unbound (red) states for LDHA (I) and ABCC1 (J), reflecting conformational changes upon binding Molecular mechanics generalized born surface area (MM-GBSA) calculations For each protein LDHA + NHI-2 and ABCC1 + Sulfinpyrazone, utilizing the molecular dynamics simulation trajectory, the binding free energy along with other contributing energy in the form of MM-GBSA were determined. According to the results (Table [173]1), ΔGbindCovalent and ΔGbindSolvGB contributed to the instability of the respective complexes, whereas ΔGbindCoulomb, ΔGbindvdW, and ΔGbindLipo were the main contributors to the ΔGbind’s role in the stability of the simulated complexes. The binding free energies of the ABCC1 + Sulfinpyrazone and LDHA + NHI-2 complexes are substantially greater, as seen in Table [174]1. These findings confirmed that NHI-2 and sulfinpyrazone can form stable protein-ligand complexes and have a high affinity for interacting with the proteins LDHA and ABCC1. Table 1. Calculated binding free energy components by MM-GBSA for the LDHA + NHI-2 and ABCC1 + Sulfinpyrazone Energies (kcal/mol) LDHA + NHI-2 ABCC1 + Sulfinpyrazone ΔG[bind] -76.04 ± 2.63 -77.35 ΔG[bind]Lipo -23.96 ± 1.03 -10.75 ΔG[bind]vdW -51.10 ± 2.0 -45.60 ΔG[bind]Coulomb -8.12 ± 1.99 -39.68 ΔG[bind]H[bond] -0.41 ± 0.22 -5.14 ΔG[bind]SolvGB 16.5 ± 1.09 25.56 ΔG[bind]Covalent 1.56 ± 1.2 1.72 [175]Open in a new tab Discussion BC is a common and deadly cancer, and drug resistance is a significant obstacle to effective treatment [[176]3].BC progression and therapy resistance are tightly linked to metabolic reprogramming and drug efflux, two hallmarks of cancer driven by molecular players like LDHA and ABCC1. Our findings revealed a consistent overexpression of these genes across multiple datasets, correlating with poor patient prognosis. Notably, LDHA and ABCC1 expression was especially elevated in triple-negative and HER2-positive subtypes—aggressive forms of BC often associated with resistance to conventional therapies.This is consistent with earlier research showing the roles of ABCC1 and LDHA in chemoresistance and carcinogenesis. A strong correlation was found between LDHA and ABCC1 expression, suggesting a cooperative mechanism driving BC aggressiveness and therapy resistance. Hypoxic conditions, often present in the tumor microenvironment (TME), may simultaneously upregulate both LDHA and ABCC1 via HIF-1α signalling, enhancing the glycolytic flux and drug efflux capacity of cancer cells. This dual upregulation creates an environment conducive to cancer cell survival, proliferation, and metastasis while reducing sensitivity to chemotherapeutic agents [[177]63]. Our integrative analysis revealed that elevated LDHA and ABCC1 expression correlates strongly with poor prognosis and lower overall survival. This prognostic association is based on predictive in silico evidence, and future validation in treatment-annotated, independent clinical cohorts is necessary to confirm these findings. These findings align with previous studies reporting their roles in chemoresistance, metastasis, and disease recurrence. Interestingly, LDHA and ABCC1 expression were not only elevated individually but also showed a positive correlation, suggesting a cooperative function in driving BC progression and resistance. Expression analysis across GTEx, TIMER 2.0, UALCAN, DISCO, and TISCH highlighted their elevation in tumor and metastatic conditions, and their enrichment in immune and stromal cell populations. Notably, while LDHA and ABCC1 are consistently overexpressed at the transcriptomic level across datasets, there is limited evidence of recurrent structural genomic alterations, such as translocations or frequent copy-number gains. This suggests that their dysregulation is likely transcriptional rather than structural, potentially driven by epigenetic changes or signaling from the tumor microenvironment. Gene interaction and enrichment analyses confirmed their roles in metabolic pathways, xenobiotic transport, and immune modulation. Importantly, while these results highlight potential interaction partners and functional clusters, we interpret them as correlative, not causative. However, the observed enrichment of immune-related and redox pathways suggests a novel functional axis through which LDHA and ABCC1 may influence both metabolic adaptation and immune evasion, a direction deserving further experimental exploration. Additionally, ENO1 was upregulated and ESR1 downregulated in LDHA/ABCC1-mutant BC, implicating metabolic-hormonal cross-talk. Immunohistochemistry validated their expression in BC tissues, supporting their roles in chemoresistance and tumor progression. The co-upregulation of LDHA and ABCC1 observed in our study suggests a converging mechanism that reinforces therapy resistance and tumor aggressiveness. LDHA, a glycolytic enzyme, contributes to the acidic tumor microenvironment by converting pyruvate to lactate, thereby facilitating immune evasion and immune checkpoint upregulation. Meanwhile, ABCC1 not only mediates efflux of chemotherapeutic agents, contributing to drug resistance, but also exports sphingosine-1-phosphate (S1P), a bioactive lipid that promotes angiogenesis, metastasis, and tumor survival. Their co-expression in therapy-resistant subtypes, such as TNBC, highlights a functional axis underpinning metabolic adaptation, immune suppression, and resistance to treatment. In light of these findings, combination therapy strategies are gaining traction. Co-targeting LDHA and ABCC1 offers a two-pronged approach: disrupting cancer’s metabolic backbone and weakening its defense against chemotherapeutic agents. Recent research demonstrates that combining LDHA inhibitors with immunotherapy improves anti-tumor responses by enhancing T-cell infiltration and reversing lactate-induced immunosuppression [[178]64]. Our in-silico docking and dynamics simulations support this dual-targeting hypothesis, identifying NHI-2 and sulfinpyrazone as promising candidates for further preclinical validation. Mechanistically, we propose that hypoxic conditions within the tumor microenvironment may upregulate both LDHA and ABCC1 via HIF-1α signaling, enhancing glycolysis and drug efflux. This dual upregulation supports tumor cell survival under metabolic stress and chemotherapy exposure, and further reinforces the need to target both molecules concurrently. To identify potential inhibitors, we employed molecular docking and MD simulation studies, which showed strong binding affinities of NHI-2 (LDHA inhibitor) and sulfinpyrazone (ABCC1 inhibitor). Binding energies, hydrogen bonding profiles, and reduced solvent-accessible surface areas suggest stable protein–ligand interactions and conformational compactness upon binding. These data support the therapeutic potential of these molecules, particularly as part of a dual-targeted therapy approach aimed at disrupting both glycolytic flux and drug efflux in BC. While bulk RNA-seq datasets such as TCGA provide valuable gene-level insights, they often obscure the cellular origin of gene expression. To resolve this, we utilized single-cell RNA-seq datasets (TISCH, DISCO) to determine the cellular localization of LDHA and ABCC1 expression. LDHA was highly expressed in tumor epithelial cells and fibroblasts, aligning with its role in glycolytic reprogramming and stromal interactions that support tumor progression. In contrast, ABCC1 showed a more heterogeneous expression pattern, spanning both tumor cells and immune populations such as CD4 + T cells, CD8 + T cells, and B cells, suggesting a dual role in mediating drug efflux in cancer cells and potentially modulating immune response. These findings contextualize our bulk RNA-seq observations, reinforcing that elevated LDHA and ABCC1 expression reflects the coordinated contribution of multiple cell types within the TME—supporting their relevance as dual therapeutic target. Our study yields relevant information on expression patterns, prognostic significances, and putative targeting of LDHA and ABCC1 in BC through comprehensive bioinformatics and in-silico approaches. It is important to acknowledge some limitations of the current study. Since this is an in-silico study, and although molecular docking and MD simulations predict NHI-2 and sulfinpyrazone to have strong binding affinities with LDHA and ABCC1 respectively, the findings remain predictive. Experimental validation is imperative to ascertain if such interactions hold biological significance and therapeutic potential experimentally. Following such validation, in future, we aim to carry out in vitro studies involving drug treatment in BC cell lines, Western blot, qPCR, and co-immunoprecipitation as checks for expression pattern and protein-protein interaction, and cell viability, apoptosis, and drug resistance as checks for dual targeting. Experimental validation will help translate computationally predicted tools into clinical means. Conclusion Our study displays the importance of LDHA and ABCC1 in BC advancement and drug resistance. The up-regulation of these genes is directly associated with poor patient prognosis and aggressive tumor phenotypes. Molecular docking and simulation studies highlighted the possibility of targeting LDHA via NHI-2 and the ABCC1 via sulfinpyrazone to disrupt cancer metabolism and drug efflux, respectively. These findings describe a promising statement for dual-targeted strategization in overcoming chemoresistance in BC. Study limitation This study is entirely based on bioinformatics and in-silico analyses, which, while powerful for hypothesis generation, have inherent limitations. The predictive nature of molecular docking and simulation does not confirm biological activity or therapeutic efficacy in vivo. Additionally, the reliance on publicly available datasets may introduce variability due to differences in data quality, sample heterogeneity, and batch effects. No experimental validation was performed to confirm the observed gene expression patterns, protein interactions, or drug responses. Future work should include in vitro and in vivo experiments to validate these findings, including functional assays, drug sensitivity tests, and mechanistic studies, to establish clinical relevance and therapeutic potential. Acknowledgements