Integrative Network-Based Transcriptomic Analysis Identifies Niclosamide as a Candidate Repositioned Drug for Breast Cancer.

Breast cancer (BC) is a heterogeneous disease with a wide range of pathogenesis and clinical characteristics. 1 BC is the most common cancer type with an incidence rate of 30%, 2 making it one of the most encountered diseases among women. According to the World Health Organization (WHO), BC is categorized into various histopathologic subgroups, including in situ and invasive breast cancer, fibroepithelial and nipple tumors, mesenchymal and hematolymphoid neoplasms, male breast tumors, and genetic tumor syndromes. 3 This classification is based on cellular morphology, growth characteristics, and architectural patterns. The molecular subtypes of BC are classified by the expression profiles of the estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2). 4 Among these, ER+ breast cancer represents the most prevalent subtype and is characterized by hormone-responsive tumor biology. 5 The MCF-7 cell line is one of the most widely used in vitro models for ER+ BC due to its stable phenotype and well-characterized estrogen-dependent signaling, making it suitable for investigating drug responses in hormone-driven BC. 6 , 7 The regulation of messenger RNAs (mRNA), microRNAs, and long non-coding RNAs are key factors involved in cancer development. 8 , 9 Recent studies investigating the identification of breast cancer biomarkers using targeted transcriptome showed that mRNA levels of HER2, ER, and AR reflect the tumor characteristics. 10 Moreover, miR-182 was identified as an oncogene contributing to breast cancer pathogenesis. 11 miRNA-mRNA interactions in both BC and triple-negative breast cancer (TNBC), highlighting the potential use of these indicators in aggressive breast cancer types by identifying miRNAs like hsa-miR-802, hsa-miR-1258, hsa-miR-548a-3p, and hsa-miR-2053 and mRNAs like MELK, NCAPG, CCNA2, and NUSAP1 that were involved in disease progression and survival outcomes. 12 These findings demonstrated that transcriptomic analysis is a more sensitive and reliable method for identifying BC-specific biomarkers and may provide new therapeutic approaches. Conventional chemotherapy remains a cornerstone in BC treatment; however, the limited selectivity of many antineoplastic agents can lead to off-target effects, potentially impacting healthy tissues and contributing to treatment-related morbidity. Therefore, novel therapeutic strategies are needed. The “drug repositioning” strategy could significantly reduce the time and cost of developing innovative drugs with previously investigated molecules to eliminate risk. 13 One crucial strategy for drug repositioning research is utilizing gene expression data linked to any genetic perturbations. 14 The traditional de novo drug development process is quite expensive, roughly $12 billion, laborious, and time-intensive, lasting approximately 10 to 15 years with high clinical attrition rates. 15 Therefore, drug repositioning has emerged as an alternative strategy for identifying new therapeutic indications for existing compounds while bypassing many early-stage safety evaluations. 16 , 17 Recent advances in breast cancer management have improved patient outcomes through increasingly personalized treatment strategies. 18 Nevertheless, disease heterogeneity, treatment resistance, and recurrence remain major clinical challenges, supporting the continued search for novel biomarkers and therapeutic targets. Several repositioned compounds, including metformin and ritonavir, have previously demonstrated potential anticancer effects in BC. 19 , 20 Niclosamide has been previously reported to exhibit anticancer activity through the inhibition of signaling pathways, such as the Wnt/β-catenin and STAT3 pathways, 21–23 while amitriptyline has been associated with mitochondrial dysfunction and cytotoxic effects in cancer cells. 24–26 However, their role within integrative transcriptomic-driven drug repositioning frameworks remains insufficiently characterized. Given these partially overlapping but mechanistically distinct biological effects, we explored whether these compounds could influence BC cell viability individually and under combinatorial conditions. Despite the growing number of computational drug repositioning studies in BC, 27 , 28 many existing approaches rely predominantly on differential expression-based drug matching and do not sufficiently integrate multilayer biological network prioritization with subsequent experimental validation. 29 , 30 In addition, subtype-specific investigations, particularly focusing on ER+ breast cancer, remain limited, particularly in integrative frameworks combining network-based prioritization with survival-associated biomarker candidates for ER+ BC. This study performed an integrative transcriptome-based microarray analysis to identify candidate repositioned small-molecule compounds associated with biologically prioritized BC signatures. A core regulatory cluster was determined by constructing a three-layered biological network using DEGs, followed by survival-based prioritization (log-rank p < 0.05) and independent expression-pattern assessment. The resulting prognostically relevant hub signatures were subsequently used as the input for L1000CDS 2 -based drug repositioning to identify compounds predicted to reverse the BC-associated transcriptional profile. Among the prioritized candidates, niclosamide emerged as the top-ranked compound and was therefore selected for preliminary in vitro validation. While amitriptyline was also included as a comparative exploratory compound, based on previously reported cytotoxic and mitochondrial effects in BC-related cellular models. 24–26 Rather than assuming a predefined synergistic mechanism, the combined evaluation of these two agents was intended to explore whether distinct but partially convergent biological activities could influence breast cancer cell viability. This study aimed to contribute to the identification of candidate therapeutic strategies for BC based on the integration of transcriptome-level data analysis, multiple-level biological network constructions, and drug repositioning.

The publicly available NCBI Gene Expression Omnibus (NCBI-GEO) database 31 was used to retrieve BC-related microarray datasets. The study design involved the use of separate datasets for (i) the identification of biomarkers associated with BC pathogenesis and (ii) the validation of the consistency of the resulting DEGs and network-derived hub molecules across independent cohorts. The microarray dataset of GSE42568 32 was composed of 121 samples, which were analyzed to elucidate the omics signatures of BC. Among these samples, 104 were taken from patients with primary tumor tissue, and 17 from normal breast tissue. Differential gene expression analysis was performed using the GEO2R online statistical platform 31 based on the submitter-provided processed series matrix expression data available in GSE42568 . According to the original dataset publication, 32 these deposited expression values were generated following GC-RMA normalization, quantile normalization, batch-effect adjustment, and probe-level filtering prior to GEO submission. GEO2R subsequently applied the GEOquery and limma Bioconductor packages to these processed data for differential comparison between user-defined sample groups. No additional raw CEL file normalization or independent probe-level preprocessing was carried out in the present study. Differential expression was performed under the default unpaired linear modeling framework implemented in GEO2R. To identify statistically significant DEGs from the selected dataset, the complete GEO2R output table was first downloaded, and statistical filtering was then applied in the downstream analysis: (i) p-value 1. We considered the cut-off for the upregulated DEGs as log 2 FC ≥ 1 and for the downregulated DEGs as log 2 FC ≤ −1. To further visualize the distribution of gene expression changes, both a volcano plot and an MA (mean–difference) plot were generated. The volcano plot was constructed using log 2 FC and p-values to graphically display significantly differentially expressed genes with substantial fold changes and statistical significance. An MA plot was further obtained from the GEO2R visualization output to display log 2 FC against average log 2 expression values. Moreover, the complete list of identified DEGs, including gene symbols, log 2 FC, p-values, and adjusted p-values, was provided as Supplementary Table 1 . The datasets GSE113865 and GSE22820 were used to assess the consistency of expression patterns of the identified hub molecules across independent cohorts. The detailed information related to the selected datasets is shared in Table 1 . Table 1 The Microarray Datasets for Transcriptomic Analysis of Breast Cancer GEO Reference Series Platform Purpose of Utilization Samples Study Design References GSE42568 Affymetrix Human Genome U133 Plus 2.0 Array Analysis 121 17 control/104 breast cancer samples [ 32 ] GSE113865 Illumina HumanHT-12 V4.0 expression bead chip Independent-validation 6 3 control/3 breast cancer samples No citation GSE22820 Agilent-014850 Whole Human Genome Microarray 4x44K Independent-validation 26 10 control /16 breast cancer samples [ 33–37 ] The Microarray Datasets for Transcriptomic Analysis of Breast Cancer Functional enrichment analysis of DEGs (upregulated and downregulated separately) was performed by the GeneCodis4 web-based bioinformatics tool. 38 The cut-off criterion was determined as adj. p-value < 0.05. For the functional annotation of DEGs, the signaling processes represented by KEGG pathways 39 were used. The Gene Ontology (GO) 40 was used to ensure a comprehensive source of the biological processes that were associated with DEGs. Biological network constructions around all DEGs were performed by collecting the interaction data from protein-protein interactions (PPIs) (BioGrid), 41 miRNA gene interactions (mirTarbase), 42 and transcription factor gene interactions (TRRUST). 43 Networks were visualized by Cytoscape software (v3.9.1), 44 and topological network analysis was performed using degree and betweenness centrality metrics. By using the “cytohubba” package of Cytoscape software, topological local and global metrics such as degree and betweenness centrality were interpreted. 45 PPIs from BioGRID were restricted to Homo sapiens ( H. sapiens ) and included experimentally validated interactions curated in the database without additional filtering by interaction type or confidence score. The miRNA–gene interactions (miRTarBase) and TF-gene interactions (TRRUST) were limited to experimentally validated entries. No tissue-specific filtering was applied; instead, a global human interaction network was constructed, and topological filtering metrics were used to identify the most relevant hub molecules. Degree and betweenness centrality were selected as topological metrics because these capture complementary aspects of network topology, namely local connectivity and global information flow, and are among the most widely applied prioritization metrics in biological network analysis. Betweenness centrality represents the number of times a node is present on the shortest path between other nodes. Degree centrality is simply the number of links held by each node. The node represents the genes of interest, while the edges represent the connections between the nodes. The top 10 hub elements for each interaction network were considered significant by topological metrics, including degree and betweenness centrality. Topological properties of the networks, including node and edge numbers, network density, and clustering coefficient, were calculated using the Network Analyzer tool 46 in Cytoscape. The hub molecules obtained from the TF–gene, miRNA–gene, and PPI interaction networks were subsequently integrated into a combined non-redundant candidate signature pool. This integrated hub set constituted the “core cluster” used for downstream prioritization analyses. To further refine the biological and clinical relevance of these signatures, Kaplan-Meier (KM) survival analysis was subsequently applied, and only molecules showing statistically significant prognostic associations (log-rank p < 0.05) were retained as the final prognostic biomarker set. KM survival analysis was performed using the KM Plotter online tool ( http://kmplot.com ), 47 which integrates gene expression and clinical data from multiple publicly available breast cancer cohorts, including the GEO and TCGA. To evaluate the prognostic significance of a particular gene, patient samples were categorized based on the gene’s median expression (high versus low expression), 48 as determined by the KM Plotter algorithm. All analyses were conducted using the default settings of the KM Plotter breast cancer module. The “all datasets” option was selected, and no additional filtering was applied based on molecular subtype or treatment. The total number of patients included in the analysis was approximately 4929, which is indicated by the platform output. Gene expression values were derived from the user-selected probe set option for each gene. Median expression values used for patient stratification were automatically computed by the KM Plotter algorithm. We performed survival analysis around core cluster elements based on the patient’s clinical data gathered from the TCGA and GEO databases. The KM Plotter 47 was used to conduct survival analysis. The core cluster elements, which have log-rank p-values<0.05, were accepted as prognostic biomarkers for BC. Survival analyses were conducted to evaluate the prognostic relevance of hub genes identified from TF, PPI, and miRNA networks. Overall survival (OS) was used as the primary endpoint. All available patients meeting the default inclusion criteria of the platform were included as a single cohort without further stratification. PCA was performed to understand the clustering ability of the prognostic biomarkers, considering diseased and healthy samples depending on their gene expression values in the GSE42568 , GSE22820 , and GSE113865 datasets. For cross-dataset comparisons, the 11 hubs identified from survival analysis were included. Gene expression values were normalized using z-score transformation (centering and scaling) to ensure comparability across datasets before performing PCA. PCA was carried out using R software (version 4.3.2) with the aid of R Studio (release 2024.12.1) and the utilization of the “factoextra” package. No explicit batch effect correction was applied, and analyses were performed separately for each dataset. PCA was provided to distinguish the samples according to their tissue origin (healthy vs diseased). Besides PCA, factor analysis was also conducted to reduce the variables by extracting all their commonalities into a smaller number of factors. The factor analysis also helped us to reveal the contribution levels of prognostic markers to the discrimination of samples. For the drug repositioning stage, the prognostically significant hub elements obtained after KM survival filtering were retained as the final BC molecular signature. These genes were subsequently categorized into upregulated and downregulated groups according to their expression tendencies in the discovery dataset and submitted to the L1000CDS 2 49 search engine to identify candidate compounds predicted to reverse the BC-associated transcriptional profile ( Supplementary Table 2 ). The L1000CDS 2 platform ranks candidate compounds using a transcriptomic overlap/reversal scoring approach based on cosine similarity metrics. Specifically, the platform evaluates the similarity between the uploaded disease-associated gene signature and perturbation-induced transcriptional profiles derived from the LINCS L1000 database. In reverse mode, compounds are prioritized according to their ability to inversely correlate with the disease signature. The ranking metric is based on the transformation 1−cos(α), where cos(α) represents the cosine similarity between the disease-associated signature vector and the perturbagen-associated expression vector. Higher overlap/reversal scores therefore indicate stronger inverse transcriptomic relationships and a greater predicted capacity to reverse the disease-associated molecular state. The top 50 resultant drugs, which were further evaluated by their association with B and were searched through publicly available databases such as DrugBank, 50 NCATS Inxight Drugs, 51 and ClinicalTrials.gov. 52 The top-scoring 10 drugs were selected for further search on their mechanism of action, approval statuses, and association with BC. After identifying repositioned drug candidates, we performed cell viability assays using niclosamide (N3510, Sigma-Aldrich) and amitriptyline (A8404, Sigma-Aldrich). Breast cancer cells, MCF-7 (ATCC, cat no: HTB-22), were cultured in high-glucose DMEM (GIBCO) with 10% fetal bovine serum (FBS) (Capricorn), 1% penicillin-streptomycin (Sigma), and 1% L-glutamine. MCF-7 cells were maintained at 37°C, in 5% CO 2 . After cells reached 80% confluency, trypsin-EDTA was used to detach cells from the flask and harvested by centrifugation at 4500 rpm for 5 minutes. 5×10 3 cells were seeded in a 96-well plate for the viability assays. MCF-7 cells were treated with niclosamide and amitriptyline at 1, 5, 25, 50, and 100 µM concentrations, and cell viability was determined after 24 h and 72 h using AlamarBlue assay at 560–590 nm wavelength 53 (n=3 replicate). Besides that, to analyze the synergistic effect of the drugs, the concentration of the drugs was determined based on the previous viability results. The drug combination was applied to MCF-7 cells under the same experimental conditions. For the analysis of synergistic effect, Bliss independence model was applied. 54 , 55 According to the Bliss definition formula: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document}$\rm {{Expected}}{\ }{{value}} = {{{M}}_{{A}}} + {\ }{{{M}}_{{B}}}-{\ }\left({{{{M}}_{{A}}}{{x}}{\ }{{{M}}_{{B}}}} \right) = {\ }{{{M}}_{{{AB}}}}$\end{document} , where Niclosamide: M A and Amitriptyline: M B . In order to show cytotoxic effect, the viability percentage on the graph is transformed to mortality percentage. Additionally, the combination effects were analysed using the Chou–Talalay method to quantitatively determine the nature of the interaction between niclosamide and amitriptyline at the corresponding effect level. The combination index (CI) was calculated according to Chou–Talalay equation described as: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document}$$\rm {{CI}} = {{{{{D}}_1}} \over {{{\left({{{{D}}_{{x}}}} \right)}_1}}} + {{{{{D}}_2}} \over {{{\left({{{{D}}_{{x}}}} \right)}_2}}}$$\end{document} D 1 and D 2 represent the doses of each drug used in combination, while (Dₓ) 1 and (Dₓ) 2 correspond to the doses of each drug alone required to produce the same effect level (e.g., 50% inhibition).

The differential gene expression analysis between the test (BC samples) and control groups revealed 4266 DEGs. Their expression distribution was further illustrated by volcano and MA plots ( Figure 1 ). The complete DEG list with corresponding statistical parameters is presented in Supplementary Table 1 . The volcano plot revealed a clear distribution of DEGs ( Figure 1A ), with significantly upregulated and downregulated genes separated based on both statistical significance and fold change thresholds. Consistent with this, the MA plot showed that most genes were centered around a log 2 FC of zero, indicating overall balanced expression between groups ( Figure 1B ). Greater variability was observed among genes with lower average expression levels, whereas genes with higher expression showed more stable fold change patterns. These findings support the reliability of the differential expression analysis and indicate that the observed changes are not driven by systematic bias. Among those, 2171 were upregulated, and 2095 were downregulated genes. Genecodis was used for functional enrichment analysis using GO and KEGG databases by defining a cut-off criterion of the adj. p-value<0.05 to elucidate the functions of the DEGs in the BC case. The functional enrichment analysis was constructed from the 10 most enriched pathways for upregulated and downregulated DEGs. Figure 1 Visualization of differential gene expression analysis. ( A ) Volcano plot showing the distribution of genes based on log 2 FC and −log10(p-value), highlighting significantly upregulated and downregulated genes. ( B ) MA plot illustrating the relationship between log 2 FC and average gene expression levels across all genes. The image A showing GSE42568 . A scatter plot with x-axis label log2(Fold Change) and y-axis label minus log10(p-value). The x-axis shows ticks at negative 5, negative 1, 0, 1, 5. The y-axis shows ticks at 0, 10, 20, 30, 40. A horizontal dashed threshold line is drawn slightly above 0 on the y-axis. Two vertical dashed threshold lines are drawn at x equals negative 1 and x equals 1. Points form a volcano shape: a dense central cluster near x around 0 and y near 0, with left and right wings extending upward. A legend lists Down, Insignificant, Up. The image B showing GSE42568 . A scatter plot with x-axis label log2Exp and y-axis label log2FC. The x-axis shows ticks at 4, 6, 8, 10, 12, 14. The y-axis shows ticks at negative 6, negative 4, negative 2, 0, 2, 4. Points form a wedge centered around y equals 0, with higher spread at mid x values and narrowing toward higher x values. Two scatter plots showing differential gene expression with a volcano plot and an MA plot. Visualization of differential gene expression analysis. ( A ) Volcano plot showing the distribution of genes based on log 2 FC and −log10(p-value), highlighting significantly upregulated and downregulated genes. ( B ) MA plot illustrating the relationship between log 2 FC and average gene expression levels across all genes. GO analysis was performed to elucidate DEGs’ biological processes and molecular functions. The downregulated DEGs were significantly enriched in biological processes such as lipid metabolic processes, brown fat cell differentiation, response to bacteria, and tricarboxylic acid cycle ( Figure 2A ). Whereas the upregulated DEGs were significantly enriched in cell division, migration, and DNA replication ( Figure 2B ). Figure 2 Enrichment analysis of DEGs. ( A ) GO-downregulated DEGs analysis. ( B ) GO-upregulated DEGs analysis. ( C ) KEGG pathway analysis of downregulated DEGs. ( D ) KEGG pathway analysis of upregulated DEGs. The number of genes associated with each pathway are indicated by n. Image A: GO Biological Process Analysis for Downregulated DEGs shows processes like lipid metabolism (91 genes, ~4.9), tricarboxylic acid cycle (18 genes, ~4.8) and response to bacterium (38 genes, ~4.6) with significance measured by -log10(Pval Adj). Image B: Upregulated DEGs include cell division (92 genes, ~6.8), apoptotic process (121 genes, ~3.2) and DNA replication (36 genes, ~3.5). Image C: KEGG Pathway Analysis for Downregulated DEGs highlights metabolic pathways (268 genes, ~9.0), propanoate metabolism (20 genes, ~7.8) and carbon metabolism (39 genes, ~7.0). Image D: Upregulated DEGs involve pathways like polycomb repressive complex (24 genes, ~2.2), human papillomavirus infection (60 genes, ~1.8) and cell cycle (34 genes, ~1.7). Each image uses -log10(Pval Adj) to indicate significance, with varying gene counts across processes. A set of four horizontal bar charts of GO and KEGG enrichment for downregulated and upregulated DEGs. Enrichment analysis of DEGs. ( A ) GO-downregulated DEGs analysis. ( B ) GO-upregulated DEGs analysis. ( C ) KEGG pathway analysis of downregulated DEGs. ( D ) KEGG pathway analysis of upregulated DEGs. The number of genes associated with each pathway are indicated by n. KEGG analysis focuses on high-level biological system functions and provides knowledge about genes, molecular pathways, and their interactions. Down-regulated DEGs were significantly enriched in metabolic and PPAR signaling pathways, propanoate and carbon metabolism, fatty acid degradation, citrate (TCA) cycle, valine, leucine, and isoleucine degradation, and focal adhesion pathways ( Figure 2C ). Upregulated DEGs were enriched for the polycomb repressive complex, tight junction, human papillomavirus infection, Fanconi anemia, estrogen signaling pathways, proteoglycans in cancers, human immunodeficiency virus-1 infection, cell cycle, and pyrimidine metabolism pathways ( Figure 2D ). To construct a TF-gene interaction network, the potential TF-gene interactions related to DEGs were revealed from the TRRUST database. The TF–gene interaction network consisted of 420 nodes and 1312 edges, with a network density of 0.015 and a clustering coefficient of 0.115. Associated interactions were visualized in the Cytoscape software to obtain hub molecules of the TF-gene interaction network based on betweenness centrality and degree metrics. Consequently, 17 hub molecules were determined. The top 10 hubs, based on betweenness centrality and degree metrics, were accepted as significant. The determined significant 17 hub molecules were AR, BCL2, CDH1, CXCL8, E2F1, EGR1, ESR1, FOS, HDAC1, IL6, JUN, NFKB1, PTGS2, RELA, SP1, STAT3 , and TP53 ( Figure 3A ). Figure 3 Three-layered biological network construction. ( A ) TF-gene interactions of DEGs. ( B ) miRNA-gene interactions around DEG. ( C ) Protein-protein interactions of DEGs. The image consists of three sub-images depicting different biological interaction networks. A shows TF-gene interactions with 420 nodes and 1312 edges, highlighting hub molecules like AR, BCL2 and TP53. B illustrates miRNA-gene interactions with 774 nodes and 2171 edges, featuring VAV3 and GATA6 as significant hubs. C presents protein-protein interactions with 9974 nodes and 191422 edges, emphasizing hubs such as CSK and NR3C1. Each sub-image includes a gene interaction spectrum indicating interaction strength, with details on network density and clustering coefficient. Three-layer network: TF-gene, miRNA-gene, protein-protein interactions with node/edge details. Three-layered biological network construction. ( A ) TF-gene interactions of DEGs. ( B ) miRNA-gene interactions around DEG. ( C ) Protein-protein interactions of DEGs. To achieve a miRNA-gene interaction network, upregulated and downregulated DEGs related to potential miRNA gene interactions were investigated through the miRTarBase database, and the H. sapiens species was used as a source. The resulting data were used to construct miRNA-gene interactions in the Cytoscape application. The miRNA-gene interaction network comprised 774 nodes and 2171 edges, showing a network density of 0.003 and a clustering coefficient of 0.000. Using the cytohubba plug-in, the data were evaluated in terms of betweenness centrality and degree centrality. The obtained significant hub molecules were hsa-miR-335-5p, hsa-miR-26b-5p, hsa-miR-124-3p, hsa-miR-16-5p, hsa-miR-92a-3p, NUFIP2, hsa-let-7b-5p, IGF1R, hsa-miR-1-3p, hsa-miR-192-5p, hsa-miR-93-5p, hsa-miR-17-5p, VAV3, hsa-miR-8485, hsa-miR-106b-5p, GATA6 ( Figure 3B ). Since upregulated and downregulated DEGs are assumed to code for the same-named proteins, DEGs were used to construct the PPI network. To detect potential protein-protein interactions within our DEGs, the BioGRID database was utilized, and as a species source, H. sapiens data was selected. Using these interactions, the PPI network was constructed on Cytoscape software, and the cytohubba plug-in was employed to calculate topological metrics such as degree and betweenness centrality. As a result, 14 significant hub molecules were obtained. The top 10 hubs from the betweenness centrality metrics and the top 10 hubs from the degree metrics are selected as protein network signatures for BC. By combining these two groups and removing the replicates, 14 significant hub molecules were revealed, including EGFR, KRAS, ESR1, TRIM25, NR3C1, HDAC1, RECQL4, TRIM28, EGLN3, RHOA, DDX39A, APP, SNCA, and CSK proteins ( Figure 3C ). The PPI network was the most extensive structure, containing 9974 nodes and 19,422 edges, with a network density of 0.000 and a clustering coefficient of 0.224. Due to the high density and complexity of the interaction networks, not all node labels are simultaneously displayed with full readability in the global network visualizations. The primary purpose of these figures is to illustrate the overall topological architecture of the biological interaction networks, while emphasizing the visibility of topologically prioritized hub molecules. KM survival analysis was conducted to understand and analyze the prognostic performance of the hub genes. A survival analysis was carried out among TF, PPI, and miRNA hub genes using the KM Plotter web-based tool. Regarding the log-rank p-value ≤ 0.05, 11 hubs were identified as prognostically significant. Among these, high expressions of ESR1 (HR= 0.64, p = 1×10 −16 ),  FOS (HR = 0.71, p = 1.3×10 −11 ),  BCL2 (HR = 0.73, p = 1.2×10 −9 ),  TRIM25 (HR = 0.73, p = 4.5×10 −5 ),  EGR1 (HR = 0.82, p = 1.1×10 −4 ),  CDH1 (HR = 0.84, p = 7.7×10 −4 ),  KRAS (HR = 0.86, p = 3.5×10 −3 ), PTGS2 (HR = 0.89, p = 2.3×10 −2 ), and IL6 (HR = 0.9, p = 4.6×10 −2 ) associated with improved OS, which may have a protective or tumor-suppressive function ( Figure 4 ). In some molecular circumstances, these biomarkers may have tumor-suppressive functions because they are commonly implicated in hormone signaling (ESR1), cell proliferation (FOS), apoptotic regulation (BCL2), cell adhesion (CDH1), antiviral response (TRIM25), tumor suppressor (EGR1), signal transduction (KRAS), inflammation (PTGS2) and immunological modulation (IL6). Although KRAS is upregulated in breast cancer and is widely recognized as an oncogene, our KM analysis showed that higher expression was associated with improved overall survival (HR = 0.86, log-rank p = 0.0035). This finding contrasts with previous reports linking elevated KRAS expression to poorer prognosis and may reflect differences in cohort composition, subtype distribution, or analytical approaches. These results suggest that the prognostic impact of KRAS expression in breast cancer may be context-dependent rather than uniformly oncogenic. In contrast, worse OS was substantially linked to increased expression of CXCL8 (HR = 1.35, p = 4.9 x 10 −9 ) and RECQL4 (HR = 1.54, p = 1.0 x 10 −16 ), suggesting a role in inflammatory signaling and tumor growth. The prognostic associations observed for certain hub molecules, particularly PTGS2 and IL6, is interpreted cautiously due to modest effect sizes and the absence of multiple-testing correction across survival comparisons. The properties of prognostic biomarkers are shared in Table 2 . Table 2 Prognostic Biomarkers Identified from Three-Layered Network Analysis Gene Name Description Regulation Association with BRCA References BCL2 BCL2 Apoptosis Regulator Down BCL2 family protein interactions regulate apoptosis and contribute to cancer development. [ 56 ] CDH1 Cadherin 1 Up Hypermethylation of CDH1 leads to reduced E-cadherin expression, a common feature in breast cancer progression. [ 57 ] CXCL8 C-X-C Motif Chemokine Ligand 8 Down The CXCL8–CXCR1/2 axis plays a key role in breast cancer initiation and progression. [ 58 ] EGR1 Early Growth Response 1 Down EGR1 expression levels have been associated with prognosis in breast cancer. [ 59 ] ESR1 Estrogen Receptor 1 Up ESR1 is associated with clinical phenotype, disease progression, and therapeutic response in breast cancer. [ 60 ] FOS Fos Proto-Oncogene, AP-1 Transcription Factor Subunit Down FOS is involved in the regulation of proliferation and progression in breast cancer. [ 61 ] IL6 Interleukin 6 Down IL6 promotes breast cancer cell proliferation through activation of the JAK/STAT3 signaling pathway. [ 62 ] KRAS KRAS Proto-Oncogene, GTPase Up Elevated KRAS expression has frequently been associated with adverse clinical outcomes in breast cancer. [ 63 ] PTGS2 Prostaglandin-Endoperoxide Synthase 2 Down Encodes the COX-2 enzyme, a key mediator of inflammation, and is frequently upregulated in aggressive breast cancer subtypes. [ 64–66 ] RECQL4 RecQ Like Helicase 4 Up RECQL4 has been identified as a breast cancer susceptibility gene, with expression levels associated with disease prognosis. [ 67–69 ] TRIM25 Tripartite Motif Containing 25 Up TRIM25 is strongly associated with breast cancer metastasis and plays a role in transcriptional regulation of metastatic processes. [ 70 ] Figure 4 Kaplan-Meier (KM) survival curves of the core cluster expression levels. All elements of the core cluster had log-rank p-values <0.05, highlighting that they might serve as potential prognostic biomarkers for BC. The image displays 11 Kaplan-Meier survival graphs based on expression groups. Each graph shows time (0-250 months) on the x-axis and probability (0.0-1.0) on the y-axis. Key findings include: ESR1 with HR 0.64, logrank P < 1e-16; RECQL4 with HR 1.54, logrank P < 1e-16; FOS with HR 0.71, logrank P 1.3e-11; BCL2 with HR 0.73, logrank P 1.2e-09; CXCL8 with HR 1.35, logrank P 4.9e-09; TRIM25 with HR 0.73, logrank P 4.5e-05; EGR1 with HR 0.82, logrank P 0.00011; CDH1 with HR 0.84, logrank P 0.00077; KRAS with HR 0.86, logrank P 0.0035; PTGS2 with HR 0.89, logrank P 0.023; IL6 with HR 0.9, logrank P 0.046. Each graph compares low and high expression curves. Eleven Kaplan-Meier survival plots comparing low versus high gene expression groups over time. Prognostic Biomarkers Identified from Three-Layered Network Analysis Kaplan-Meier (KM) survival curves of the core cluster expression levels. All elements of the core cluster had log-rank p-values <0.05, highlighting that they might serve as potential prognostic biomarkers for BC. PCA was conducted to determine the ability to discriminate between healthy and BC samples, considering the network biomarkers. The GSE42568 , GSE113865 , and GSE22820 datasets were used in PCA analysis to better understand the gene expression data. The PCA graph of variables for each dataset ( Figure 5 ) demonstrates that the healthy controls and BC samples are separated on the principal components. This indicates that different biological circumstances or cellular states can be observed in the dataset. The principal component 1 (PC1) and PC2 explained 53.4% of the total variance in the GSE42568 dataset, with the highest contribution values for explaining the variance found in the IL6, EGR1 , and FOS genes ( Figure 5A ). Although the GSE113865 dataset contains a limited number of samples, it was included as an independent exploratory cohort to assess the consistency of expression-pattern separation across datasets rather than to provide definitive statistical validation. In the GSE113865 data set, the PCA analysis explained 76.7% of the total variance ( Figure 5C ). This analysis determined the genes with the highest contribution values as KRAS, FOS , and PTGS2 . In the GSE22820 dataset, 55.1% of the total variance was explained, and the genes with the highest contribution values were determined as FOS, EGR1 , and KRAS ( Figure 5B ). The significance of gene contributions for explaining the variance of the datasets is estimated using the cos 2 values. A low cos 2 value means that the variable is not perfectly represented by that component (in our case, hub elements). A high cos 2 value, on the other hand, means a good representation of the variable on that component. The variables presented in red and localized away from the circle’s origin exhibited the higher cos 2 values. In contrast, variables that were presented in blue color and closer to the circle’s origin have lower cos 2 values, which exhibit less significance. Multiple-testing correction was not applied in the present survival analyses because the primary objective was exploratory prioritization of candidate biomarkers for downstream network integration and drug repositioning rather than definitive prognostic modeling. We acknowledge that this approach may increase the risk of false-positive findings; however, applying stringent multiple-testing corrections at this stage could also increase false-negative results and potentially exclude biologically relevant candidates for subsequent analyses. These results indicate that EGR1, FOS, IL6, KRAS , and PTGS2 demonstrated sample-separation capability to discriminate between diseased and healthy samples. It specifically highlights the biological significance of these cancer-related genes and the requirement for them to be assessed as potential diagnostic biomarkers. Figure 5 Principal Component Analysis plots for ( A ) GSE42568 , ( B ) GSE22820 , ( C ) GSE113865 . Variables are colored according to their cos 2 values, where blue indicates low contribution and red indicates high contribution to the principal components. The image contains three sets of graphs for datasets GSE42568 , GSE22820 and GSE113865 . Each set includes a PCA scatter plot and a factor analysis plot. In GSE42568 , the PCA plot shows PC1 (34 percent) and PC2 (19.4 percent), highlighting sample clustering and separation. The factor analysis plot indicates strong contributions from genes like IL6, EGR1 and FOS. In GSE22820 , the PCA plot with PC1 (37 percent) and PC2 (18.8 percent) shows distinct sample separation. The factor analysis plot highlights KRAS and FOS as major contributors. In GSE113865 , the PCA plot with PC1 (55.5 percent) and PC2 (21.2 percent) shows limited sample separation. The factor analysis plot identifies FOS and TRIM25 as key contributors. Color scales represent cos2 values for individuals and contribution levels for variables. Across datasets, GSE113865 captures the most variance in PC1, while GSE42568 shows the most distinct sample clustering. PCA/factor plots for GSE42568 , GSE22820 , GSE113865 show sample separation and gene roles. Principal Component Analysis plots for ( A ) GSE42568 , ( B ) GSE22820 , ( C ) GSE113865 . Variables are colored according to their cos 2 values, where blue indicates low contribution and red indicates high contribution to the principal components. To elucidate whether potential diagnostic biomarkers hold a feature of being treatment targets of the disease, we examined them further by utilizing the L1000CDS 2 search engine. It prioritizes thousands of small-molecule signatures and their pairwise combinations using two techniques to reverse or mimic an input gene expression profile. To reverse the BC disease scenario into a healthy state, the tool uses a reverse mode that oppositely converts input gene expression signatures. The inputs were selected as up- and down-regulated potential diagnostic biomarkers. Based on KM survival filtering, 11 prognostically significant hub genes were retained for the drug repositioning analysis. Among these, CDH1, ESR1, KRAS, RECQL4 , and TRIM25 constituted the upregulated input signature, whereas PTGS2, BCL2, CXCL8, EGR1, FOS , and IL6 formed the downregulated input signature submitted to L1000CDS 2 . The tool gave 50 repositioned drug candidates as a result. The candidate drugs from the query results were eliminated based on their approval status, mode of action, and indications. Among all resultant drugs, 4 of them were approved by the FDA, and the other four drugs were stated as investigational. Table 3 summarizes the detailed properties of selected repositioned drug candidates by giving their mechanisms of action, FDA approval status, and indications. Table 3 Repurposed Drug Candidates Based on the BC-Specific Hub Genes Perturbation Indication MOA Approval Status References Niclosamide Taeniasis due to Taenia solium and Taenia saginata , hymenolepiasis, diphyllobothriasis The regulation of Wnt/β-catenin, mTORC1, STAT3, NF-κB, and Notch signaling pathways, along with oxidative phosphorylation Approved [ 23 , 71–74 ] Emetine Hydrochloride Protozoal infections, acute fulminating amebic dysentery, and amebic hepatitis or abscess Induces vomiting, Inhibits protein synthesis and the synthesis and activities of nucleic acids Approved [ 75 , 76 ] Cycloheximide It is not suitable for human use as a treatment. Protein synthesis inhibitor Investigational [ 77 ] Periplocymarin Cardiac insufficiency and cardiogenic hypotension. Its cardiotonic mechanism involves targeting Na + -K + -ATPase to raise the Ca 2+ concentration in cardiomyocytes Investigational [ 78 ] Narciclasine Melanoma Inhibitor of topoisomerase I Investigational [ 79 , 80 ] Anisomycin An antimicrobial agent Blocks protein synthesis by inhibiting peptidyl transferase activity of the 60S ribosomal subunit Investigational [ 81 ] Penfluridol Schizophrenia T-type Ca 2+ channel blocker. It is thought to be a dopamine receptor blocker Approved [ 82 , 83 ] Ouabain Cardiovascular disease Inhibits the Na + -K + -ATPase membrane pump Investigational [ 84 ] Amitriptyline Major depressive disorder A tricyclic antidepressant by inhibiting norepinephrine and serotonin reuptake Approved [ 85 , 86 ] Repurposed Drug Candidates Based on the BC-Specific Hub Genes MCF-7 cells were treated with varying concentrations (1, 5, 25, 50, and 100 µM) of repositioned drugs, including niclosamide and amitriptyline (n = 3 replicate). Although a broad concentration range of drug was tested to evaluate dose-dependent effects, some of the higher concentrations may be above clinically achievable plasma levels. This concentration range was tested to fully evaluate the cellular response profile and determine the effective dose window of drugs. Both drugs exhibited significant cytotoxicity at both 24 h and 72 h. At 24 h, amitriptyline inhibited cell proliferation with an IC 5 0 value of 43.12 (95% CI: 10.88 to 191.7) ( Supplementary Figure 1 ), corresponding to approximately 50% inhibition at 25 µM ( Figures 6A and B ). Niclosamide exhibited a potent antiproliferative effect, with an IC 5 0 value of 1.910 µM (95% CI: 1.419 to 2.624) at 24 h ( Supplementary Figure 2 ). Notably, treatment with 5 µM niclosamide resulted in a marked reduction in cellular proliferation, consistent with its high cytotoxic potency ( Figure 6C, D and Supplementary Figure 3 ). The results from the 24-h observation were consistent with the findings at 72 h. After 72 h incubation of cells with amitriptyline, the effect on cellular viability remained the same, whereas for niclosamide, cellular proliferation was further decreased. Figure 6 Viability analysis and fluorescence microscopy of MCF-7 cells following treatment with amitriptyline and niclosamide at 24 h and 72 h. ( A ) Cell viability profiles of amitriptyline-treated cells at 24 h and 72 h. ( B ) Representative fluorescence microscopy images of amitriptyline-treated MCF-7 cells at 24 h and 72 h. ( C ) Cell viability profiles of niclosamide-treated cells at 24 h and 72 h. ( D ) Representative fluorescence microscopy images of niclosamide-treated MCF-7 cells at 24 h and 72 h. Error bars represent the standard error of the mean (SEM) from two independent experiments (n = 2), each performed with three technical replicates. Statistical significance was determined using one-way ANOVA (*p < 0.05, **p < 0.01, ***p < 0.001, ****p<0.0001). Scale bar: 50 µM. Image A features two bar charts, Amp.24h and Amp.72h, with x-axis labeled Concentration (micro-M) and categories CTRL, 1, 5, 25, 50, 100. The y-axis shows percent Viability from 0 to 150. Amp.24h: CTRL ~95; 1 ~70; 5 ~60; 25 ~50; 50 ~0; 100 ~0, with significance asterisks. Amp.72h: CTRL ~95; 1 ~85; 5 ~75; 25 ~60; 50 ~0; 100 ~0, with asterisks. Image B displays fluorescence microscopy fields for Amp.24h and Amp.72h across concentrations 100, 50, 25, 5, 1 micro-M and Control. Sparse dots at 100 and 50, increasing at 25 and 5, densest at 1 and Control; 72h row is sparser than 24h at same doses. Image C shows Nic.24h and Nic.72h bar charts with similar axes. Nic.24h: CTRL ~95; 1 ~70; 5 ~25; 25 ~20; 50 ~10; 100 ~10, with asterisks. Nic.72h: CTRL ~95; 1 ~95; 5 ~5; 25 ~5; 50 ~5; 100 ~0, with asterisks. Image D shows microscopy fields for Nic.24h and Nic.72h, similar layout as Image B. Few dots from 100 to 5, more at 1, densest in Control; 72h row is sparse at higher doses, dense in Control. Four-part figure of MCF-7 viability and fluorescence images under amitriptyline and niclosamide treatment doses. Viability analysis and fluorescence microscopy of MCF-7 cells following treatment with amitriptyline and niclosamide at 24 h and 72 h. ( A ) Cell viability profiles of amitriptyline-treated cells at 24 h and 72 h. ( B ) Representative fluorescence microscopy images of amitriptyline-treated MCF-7 cells at 24 h and 72 h. ( C ) Cell viability profiles of niclosamide-treated cells at 24 h and 72 h. ( D ) Representative fluorescence microscopy images of niclosamide-treated MCF-7 cells at 24 h and 72 h. Error bars represent the standard error of the mean (SEM) from two independent experiments (n = 2), each performed with three technical replicates. Statistical significance was determined using one-way ANOVA (*p < 0.05, **p < 0.01, ***p < 0.001, ****p<0.0001). Scale bar: 50 µM. To assess potential synergistic effects between these two drugs, the MCF-7 cells were treated with these drugs in combination at critical concentrations. The combined treatment model was employed to potentiate therapeutic efficacy using niclosamide and amitriptyline at concentrations of 1:25 µM and 5:25 µM, respectively. Under the fluorescence microscope, the morphology of the cells and relative fluorescence units (RFUs) were consistent with the cellular viability results. As a result, the higher treatment potential of the drugs against BC was confirmed. The interaction between niclosamide and amitriptyline was initially evaluated using the Bliss independence model. According to the Bliss definition formula: M A + M B – M A x M B = M AB , where Niclosamide: M A and Amitriptyline: M B , • for the niclosamide on day 1: at 1 µM viability percentage is 72.43%, and M A = 27.57% • for the amitriptyline on day 1: at 25 µM viability percentage is 50.17%, and M B = 49.83% o 0.27 + 0.49 – (0.27 X 0.49) =0.62 (Expected Value) • for the niclosamide on day 1: at 5 µM viability percentage is 23.35%, and M A = 76.65% • for the amitriptyline on day 1, at 25 µM viability percentage is 50.17%, and M B = 49.83% o 0.76 + 0.49 – (0.76 X 0.49) = 0. 88 (Expected Value) for the niclosamide on day 1: at 1 µM viability percentage is 72.43%, and M A = 27.57% for the amitriptyline on day 1: at 25 µM viability percentage is 50.17%, and M B = 49.83% o 0.27 + 0.49 – (0.27 X 0.49) =0.62 (Expected Value) 0.27 + 0.49 – (0.27 X 0.49) =0.62 (Expected Value) for the niclosamide on day 1: at 5 µM viability percentage is 23.35%, and M A = 76.65% for the amitriptyline on day 1, at 25 µM viability percentage is 50.17%, and M B = 49.83% o 0.76 + 0.49 – (0.76 X 0.49) = 0. 88 (Expected Value) 0.76 + 0.49 – (0.76 X 0.49) = 0. 88 (Expected Value) The expected effect value of the two drug concentrations at 1:25 µM is 62%, and the expected effect at 5:25 µM is 60%. The measured cytotoxicity levels were lower than the predicted values for 1:25 and 5:25 µM drug concentrations, indicating that the expected level is not reached. The graphical representation of the synergistic effects of the drugs was shared in Figure 7A and B (Fluorescence intensity of cellular images was indicated in Supplementary Figures 3 and 4 ). The results imply that the combinatorial treatment fails to achieve the anticipated additive cytotoxicity, suggesting a possible antagonistic interaction between niclosamide and amitriptyline at both tested ratios. Figure 7 Bliss Independence analysis of drug combinations in MCF-7 cells. ( A ) Mortality rates of single and combined drug treatments at two different concentrations (5:25 μM and 1:25 μM. ( B ) Representative fluorescence images of MCF-7 cells under combined drug treatment. Scale bar: 50 µm. Two bar charts display percent Effect on the y-axis over 24 hours. The first chart uses x-axis label Concentration (5:25 micromolar) with three bars for A, N and N:A. A dashed Expected Effect line sits near 90 percent. Bar A reads approximately 50 percent, bar N approximately 75 percent and bar N:A approximately 65 percent, all falling below the expected line. The second chart uses x-axis label Concentration (1:25 micromolar) with the same three categories. A dotted Expected Effect line sits near 62 percent. Bar A reads approximately 50 percent, bar N approximately 27 percent and bar N:A approximately 15 percent, all below the expected line. Four circular fluorescence microscopy images of MCF-7 cells are labeled 24h and Nic:Amp. Conditions shown are 5:25 micromolar, 2.5:25 micromolar, 1:25 micromolar and Control. A scale bar of 50 micrometers is visible. The Control image shows visibly more fluorescent cells compared to the treated conditions. Two bar charts and four fluorescence images show MCF-7 cell drug effects at 5:25 and 1:25 micromolar over 24 hours. Bliss Independence analysis of drug combinations in MCF-7 cells. ( A ) Mortality rates of single and combined drug treatments at two different concentrations (5:25 μM and 1:25 μM. ( B ) Representative fluorescence images of MCF-7 cells under combined drug treatment. Scale bar: 50 µm. The interaction profile was quantified further using the Chou–Talalay combination index (CI). The calculated CI value (>1) suggests that the two compounds interact antagonistically at this effect level. The 95% confidence interval ranged from 2.04 to 5.82, indicating the robustness of this finding. Analysis by the Chou-Talalay method. The CI was calculated specifically to determine the nature of the interaction at the corresponding level of effect. 54 , 87 , 88 Using the (Chou-Talalay) Equation: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document}$${\mathrm{CI}} = {{{{\mathrm{D}}_1}} \over {{{\left({{{\mathrm{D}}_{\mathrm{x}}}} \right)}_1}}} + {{{{\mathrm{D}}_2}} \over {{{\left({{{\mathrm{D}}_{\mathrm{x}}}} \right)}_2}}}$$\end{document} D 1 and D 2 represent the doses of each drug used in combination, while (Dₓ) 1 and (Dₓ) 2 correspond to the doses of each drug alone required to produce the same effect level (e.g., 50% inhibition). For the tested condition: D 1 (Amitriptyline) = 25 µM, (Dₓ) 1 = 43.12 µM, D 2 (Niclosamide) = 5 µM, (Dₓ) 2 = 1.91 µM The combination index (CI) was calculated as: CI = (25 / 43.12) + (5 / 1.91) = 3.19 CI = (25 / 43.12) + (5 / 1.91) = 3.19 One possible explanation is that overlapping or competing molecular targets may diminish the efficacy of each compound when applied together. Alternatively, differences in pharmacodynamics, such as variations in drug uptake, metabolism, or intracellular bioavailability, may account for the reduced cytotoxic response. It is also conceivable that compensatory signaling pathways are activated under combined treatment, thereby attenuating the expected cytotoxic impact. These findings highlight the necessity for further mechanistic studies, including pathway-specific analyses and dose response modeling, to clarify the interaction profile and optimize combinatorial regimens.

Integrating omics-level data analysis, network construction, and a drug repositioning framework enabled the identification of BC-associated candidate biomarkers, including CDH1, ESR1, KRAS, RECQL4, TRIM25, PTGS2, BCL2, CXCL8, EGR1, FOS , and IL6 genes, with potential biological and prognostic relevance. These biomarkers were associated with survival outcomes and contributed to expression-based separation between breast cancer and normal samples, suggesting potential biological relevance that warrants further investigation. Niclosamide was prioritized as a repositioned candidate through transcriptomic signature-reversal analysis and demonstrated preliminary antiproliferative activity in MCF-7 cells. Given its reported effects on inflammation, proliferation, and oxidative phosphorylation pathways substantially, it was selected for in vitro evaluation and demonstrated cytotoxic effects on MCF7 cells (24 h, 5 µM). Based on previously reported similarities in mitochondrial and inflammatory pathway modulation, amitriptyline was additionally evaluated as a comparative compound in the same experimental conditions. Amitriptyline also showed dose-dependent cytotoxic activity (24 h, 25 µM), although at a higher IC 50 concentration than niclosamide. Overall, these findings support niclosamide emerged as a highly prioritized candidate and demonstrated preliminary antiproliferative activity in the MCF-7 BC model; however, thorough validation across additional BC subtypes and normal breast epithelial cell models, mechanistic and in-depth in vivo studies are required before translational interpretation.

Over the past decades, the therapeutic landscape for BC has undergone substantial evolution, shifting from localized interventions toward a more integrated, multimodal paradigm aimed at controlling both locoregional disease and distant metastases. More recently, the convergence of high-throughput omics technologies and computational biology has fostered the emergence of drug repositioning as a promising therapeutic avenue. To identify therapeutic targets and transcriptional patterns, we compared 104 BC samples to 17 normal breast biopsies. We detected a total of 4266 DEGs, of which 2171 were upregulated, whereas 2095 of them were downregulated. Around these DEGs, biological network constructions revealed 37 BC-specific network signatures that provide insights into potential disease mechanism(s). Among these 37 network signatures, survival analysis revealed the potential BC-specific prognostic biomarkers, including AGR2, ANLN, AR, BCL2, CALM1, CANX, CDH1, CXCL8, E2F1, EGR1, ERBB2, ESR1, EZH2, FASN, FOS, IGF1R, IL6, KRAS, MIDN, MMP9, MOV10, NFKB1, NPM1, PPARG, PTGS2, RECQL4, RELA, SNCA, SOX4, SQSTM1, STAT3, TP53, TRIM25, TXNIP, UBN2, VAV3 , and XPO1 ( Figure 4 , Supplementary Figure 5 and 6 ). Although several prioritized signatures (ESR1, BCL2, KRAS, PTGS2, and IL6) identified in this study have previously been implicated in BC biology, their recovery through the present multilayer network-based framework supports the biological robustness and internal consistency of the analytical strategy. Importantly, the integration of transcriptomic alterations with multiple regulatory interaction layers may facilitate the prioritization of functionally connected neighboring signatures that might remain undetected under conventional single-threshold differential expression analyses. In addition, the incorporation of survival-based prioritization and drug repositioning analyses extends the utility of the framework beyond mechanistic interpretation toward translational therapeutic candidate identification. Among these network signatures, 11 of them showed a statistically significant impact on the survival probabilities of BC patients. To demonstrate the discrimination ability of these 11 prognostic biomarkers regarding the healthy and disease states of the samples, PCA analysis was carried out and showed separating ability of BC cancer patients and healthy controls above a variance of 50% levels ( Figure 5 ). Independent datasets were primarily used to assess the consistency of the findings rather than to provide full external validation. Differences in platform technologies may also introduce variability that could not be completely accounted for across independent datasets. The 11-gene signature identified in this study should currently be interpreted as an exploratory and biologically prioritized molecular signature rather than a clinically validated diagnostic, prognostic, or treatment-guiding panel. In the present framework, these genes contributed to expression-based sample separation, showed survival-associated patterns, and served as transcriptomic inputs for drug repositioning analysis. After revealing these potential molecular signatures, a drug repositioning analysis was conducted to reverse the scenario of BC pathogenesis by targeting and reversing the gene expression profiles of all diagnostic signatures. Drug repositioning analysis revealed specific small molecules that have the potential to treat BC. In this study, 8 potential drug candidates ( Table 3 ) were identified through drug repositioning based on gene expression profiles associated with breast cancer. Potential drug candidates were identified based on significant overlap values, including niclosamide, emetine hydrochloride, cycloheximide, periplocymarin, narciclasine, anisomycin, penfluridol, and ouabain. Although several identified compounds (cycloheximide, periplocymarin, narciclasine, anisomycin, ouabain) remain investigational or are not currently approved for oncology-related clinical use, previous studies have demonstrated potential anticancer activities for some of these molecules in experimental cancer models. Narciclasine has been reported to induce autophagy-dependent apoptosis and inhibit STAT3-associated signaling in breast cancer-related systems, whereas anisomycin has shown anti-proliferative and pro-apoptotic activity across multiple cancer models. 89 , 90 In addition, several investigational compounds identified in the repositioning analysis have previously demonstrated experimentally supported anticancer activities. Periplocymarin, a cardiac glycoside-derived natural compound, has been reported to inhibit tumor proliferation, induce apoptosis, and modulate glycolysis and mitochondrial oxidative phosphorylation pathways in multiple cancer models through PI3K/AKT and MAPK/ERK signaling regulation. 91 Similarly, ouabain has been shown to exert antiproliferative and pro-apoptotic effects in BC and other malignancies through modulation of Na + /K + -ATPase-associated signaling, ERK1/2 activation, STAT3 suppression, and cell-cycle regulatory pathways. 92 Cycloheximide has also been widely used in experimental oncology studies because of its potent inhibition of protein synthesis, despite its limited clinical applicability due to toxicity concerns. 93 Therefore, these compounds were retained primarily as biologically informative computational outputs rather than immediately translatable therapeutic candidates. Among the identified compounds, niclosamide emerged as a prioritized candidate for further investigation based on its favorable overlap score, known pharmacological profile, and prior FDA approval status. Niclosamide is an anthelmintic drug used to treat parasitic infections, and it was approved by the FDA in 1982. 94 Niclosamide is a member of the weakly acidic lipophilic group which is known as salicylanilides. 95 The salicylanilides disrupt the synthesis of adenosine triphosphate by uncoupling the oxidative phosphorylation in the cell mitochondria; therefore, the motility of parasites and perhaps other functions as well are impaired. 96 The mechanism of action of niclosamide is thought to work against parasites by inhibiting the oxidative phosphorylation of mitochondria and the formation of anaerobic ATP, and glucose uptake. 97 , 98 Niclosamide is the potent mitochondrial uncoupler class 73 with the function of inhibiting various biological processes and signaling pathways, including mTORC1, nuclear factor-κB (NF-κB), Wnt/β-catenin, Notch, and signal transducer/activator of transcription 3 (STAT3). 23 , 71–74 Over the past few years, increasing data suggest that niclosamide is a multi-functional drug, raising the possibility that it could be developed as a new treatment for conditions other than helminthic diseases. 99 Wu et al demonstrated that niclosamide might be repositioned in colorectal cancer treatment. 100 Amitriptyline was included as a literature-supported compound with previously reported cytotoxic and mitochondrial regulatory effects in BC cells to provide an additional comparative perspective in the experimental phase. Amitriptyline, a tricyclic antidepressant, was approved as a drug by the FDA in 1961 under the brand name Elavil ® , 101 and is used for the treatment of major depressive disorder (MDD) in adults. 86 The function of amitriptyline is a reuptake inhibitor of serotonin and norepinephrine, exerting significant effects on the serotonin transporter and mild effects on the norepinephrine transporter. 102 , 103 It has also been used to treat post-COVID headaches. 104 Amitriptyline, a tertiary amine, blocks serotonin and norepinephrine reuptake and has strong binding affinities for muscarinic (M1), histamine (H1), and alpha-adrenergic receptors. 105 Amitriptyline is used to treat depression in most breast cancer patients. 25 The cytotoxic effect of amitriptyline on the viability of MCF7 breast cancer cells was examined in vitro, and as a result, it was determined that amitriptyline had a significant cytotoxic effect on these cells. 106 The computational analyses were performed without subtype-specific stratification and therefore reflect molecular signatures associated with breast cancer at a broader level. MCF-7 cells were selected for experimental validation because they are a well-established ER+ breast cancer model. Consequently, while the identified biomarkers emerged from non-stratified breast cancer datasets, the in vitro validation findings should be interpreted within an ER+ cellular context. MCF-7 cells are characterized by ESR1 expression and hormone-responsive behavior, making it highly suitable for investigating ESR1-associated pathways and therapeutic responses. 107 However, since different molecular subtypes of breast cancer exhibit distinct biological characteristics, the direct generalization of these findings to other subtypes, such as triple-negative or HER2+ breast cancer, may be limited. Therefore, further validation in additional breast cancer models representing different molecular subtypes will be necessary in future studies. Given their reported involvement in mitochondrial regulation and inflammatory signaling, niclosamide and amitriptyline were considered biologically relevant compounds for exploratory combination assessment. According to studies, niclosamide decreases ATP production by interfering with oxidative phosphorylation in the mitochondria. In this connection, amitriptyline is known as non-selective monoamine reuptake inhibitors (NSMRIs), 26 which inhibit ATP production and NADH oxidation, which are traits of uncouplers of oxidative phosphorylation. 108 Amitriptyline significantly inhibits mitochondrial complex I-II-linked respiration activity. 26 Both drugs have been reported to modulate inflammatory responses. Recent studies in experimental animals and human models of acute inflammation have shown that amitriptyline has anti-inflammatory properties in addition to its therapeutic uses as an analgesic and antidepressant. 109 Scheuermann et al indicated that amitriptyline efficiently inhibits inflammation in the mouse sponge model. 110 Niclosamide has been shown to directly inhibit the DNA-binding domain of STAT3, thereby modulating key cellular processes such as proliferation and apoptosis, while also suppressing NFĸB signaling, downregulating Wnt, mTOR, and STAT3 pathway proteins, and consequently attenuating inflammatory responses, reducing macrophage-induced viability and cytokine/chemokine secretion in human endometriotic stromal cells, and markedly inhibiting tumor cell growth in both a mouse model of endometriosis and ovarian cancer models [94–96]. Niclosamide may cause acute myeloid leukemia (AML) blast cells to undergo apoptosis by inhibiting the NFκB pathway and raising the reactive oxygen species (ROS) production. 74  It has been proven as a result of an in vitro study with astroglial cell lines from mice that amitriptyline inhibits NF-κB translocation and decreases IL-1β. 111 Along with the literature, the findings of this study indicate that the combinatorial regimen did not achieve the anticipated additive cytotoxicity, suggesting a potential antagonistic interaction between niclosamide and amitriptyline at both tested ratios. Similar outcomes have been documented in drug combination studies, where non-additive or antagonistic effects arose due to overlapping or competing molecular targets that limited therapeutic efficacy when agents were co-administered. 112 , 113 Both compounds have previously been associated with modulation of mitochondrial oxidative phosphorylation and cellular energy metabolism, raising the possibility that simultaneous mitochondrial stress induction may activate adaptive metabolic compensation mechanisms that partially preserve cellular viability. In addition, mitochondrial dysfunction-induced stress responses may trigger compensatory survival-associated pathways, including PI3K/AKT, MAPK/ERK, autophagy-related signaling, or redox-regulatory processes that attenuate the expected additive cytotoxic effect. Another possibility is that alterations in ATP availability, mitochondrial membrane potential, or ROS dynamics may differentially influence the cellular responses to combined treatment conditions. Another plausible explanation may lie in pharmacodynamic discrepancies, including altered uptake, metabolism, or intracellular bioavailability, all of which are known to critically influence drug–drug interactions in breast cancer and other malignancies. 114 , 115 Furthermore, combined therapies have been reported to activate compensatory signaling cascades such as PI3K/AKT or MAPK pathways, that counteract intended cytotoxic effects and attenuate overall treatment efficacy. 116–118 Collectively, these findings highlight the need for further mechanistic studies such as pathway-specific analyses and dose–response modeling, mitochondrial functional assays, apoptosis profiling, ROS quantification, metabolic flux analyses, and pathway-specific transcriptomic or proteomic studies to clarify the molecular basis of the observed antagonistic interaction and to refine combinatorial strategies for translational application in BC. Also, the selectivity of the combinatorial regimen should be thought about carefully because the best treatment plan should focus on cancer cells and not harm normal cells as much as possible. The observed combination’s lack of selectivity may make it even less useful in the clinic, especially when it comes to antagonistic interactions. The observed IC 50 values indicate that niclosamide exerted a comparatively stronger antiproliferative response than amitriptyline under identical experimental conditions, although these concentration-dependent effects should be interpreted cautiously due to the use of a single MCF-7 cell model. Since MCF-7 cells represent a luminal ER+ subtype, the present in vitro findings remain preliminary and cannot be directly generalized to other molecular forms of breast cancer without further validation. A limitation of this study is the lack of subtype-specific stratification (eg, ER, HER2, and TNBC status, as well as clinical variables), which may mask subtype-dependent molecular differences and influence therapeutic interpretation. The present findings should be interpreted within the context of current BC therapeutic strategies. Niclosamide and amitriptyline are not proposed as replacements for established standard-of-care therapies such as endocrine treatment, CDK4/6 inhibition, or HER2-targeted approaches. Instead, these compounds may represent exploratory repositioning candidates with potential adjunctive or combinatorial relevance, particularly in luminal/ER+ BC models such as MCF-7. Given the known involvement of inflammatory signaling, mitochondrial metabolism, and proliferation-associated pathways in therapeutic resistance, further studies may help determine whether these agents could function as sensitizers or complementary therapeutic components alongside existing treatment regimens. Despite the promising transcriptomic and in vitro findings observed in the present study, important translational limitations should be considered. The clinical implementation of given 11-gene signatures would require a stepwise validation pathway, including ROC/AUC-based diagnostic performance assessment, multivariable prognostic modeling, subtype-specific validation, prospective cohort evaluation, and direct testing of whether the signature predicts therapeutic response to prioritized repositioning candidates. Therefore, the present findings provide a hypothesis-generating basis for future biomarker and therapeutic validation studies rather than an immediately applicable clinical decision tool. In addition, niclosamide is known to exhibit poor aqueous solubility, limited oral bioavailability, and low systemic exposure, 99 all of which substantially restrict its clinical applicability despite its broad-spectrum anticancer activity in preclinical models. Several studies have therefore focused on improving niclosamide pharmacokinetics through nanoformulations, prodrug approaches, and alternative delivery systems designed to enhance systemic absorption and tissue distribution. 119 In parallel, amitriptyline is associated with dose-dependent central nervous system and anticholinergic adverse effects, including sedation, dizziness, dry mouth, constipation, and cognitive impairment, which may limit tolerability at concentrations potentially required for anticancer applications. 105 , 120 Therefore, the present findings should currently be interpreted as exploratory repositioning evidence rather than immediate clinical treatment recommendations. Consequently, additional mechanistic investigations, including pathway-targeted gene/protein expression analyses, metabolic profiling, mitochondrial function assays, and subtype-specific validation experiments will be required to clarify the molecular basis of the observed antiproliferative effects. Future translational studies may require optimized formulation strategies, targeted delivery systems, dose-adjustment approaches, or combination-based regimens to improve therapeutic feasibility and safety profiles.