Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Bioinform., 16 September 2025

Sec. Drug Discovery in Bioinformatics

Volume 5 - 2025 | https://doi.org/10.3389/fbinf.2025.1661601

Drug discovery for chemotherapeutic resistance based on pathway-responsive gene sets and its application in breast cancer

Dehua Feng&#x;Dehua Feng1Jingwen Hao&#x;Jingwen Hao2Lingxu Li&#x;Lingxu Li1Jian ChenJian Chen1Xinying LiuXinying Liu1Ruijie ZhangRuijie Zhang1Huirui HanHuirui Han1Tianyi LiTianyi Li1Xuefeng WangXuefeng Wang1Xia LiXia Li1Lei YuLei Yu1Bing Li
Bing Li1*Jin Li
Jin Li1*Limei Wang
Limei Wang1*
  • 1School of Intelligent Medicine and Technology, Kidney Disease Research Institute, Hainan Engineering Research Center for Health Big Data, Hainan Medical University, Haikou, Hainan, China
  • 2Key Laboratory of Brain Science Research and Transformation in Tropical Environment of Hainan Province, Department of Anatomy, School of Basic Medicine, Hainan Medical University, Haikou, Hainan, China

Introduction: Chemotherapy response variability in cancer patients necessitates novel strategies targeting chemoresistant populations. While combinatorial regimens show promise through synergistic pharmacological interactions, traditional pathway enrichment methods relying on static gene sets fail to capture drug-induced dynamic transcriptional perturbations.

Methods: To address this challenge, we developed the Pathway-Responsive Gene Sets (PRGS) framework to systematically identify chemoresistance-associated pathways and guide therapeutic intervention. Comparative evaluation of three computational strategies (GSEA-like method, Hypergeometric test-based method, Bates test-based method) revealed that the GSEA-like methodology exhibited superior performance, enabling precise identification of drug-induced pathway dysregulation.

Results: Key experimental findings demonstrated PRGS’s superiority over conventional Pathway Member Gene Sets (PMGS), exhibiting statistical independence (p < 0.0001) and enhanced detection of chemotherapy-driven pathway dysregulation. Application of PRGS to the GDSC dataset identified 8 resistance-associated pathways. Screening of agents targeting these pathways yielded candidates with predicted anti-resistance activity. An in vitro cellular experiment demonstrated that the bortezomib-bleomycin combination exhibited synergistic cytotoxicity (IDAcomboScore = 0.014) in T47D cells, highlighting the potential of PRGS-guided therapeutic strategies.

Discussion: This study establishes a PRGS-based methodological framework that integrates genomic perturbations with precision oncology, demonstrating its capacity to decode resistance mechanisms and guide therapeutic development through dynamic pathway analysis.

1 Introduction

Chemoresistance remains a critical obstacle in cancer treatment, driven by the multifaceted mechanisms through which cancer cells subvert chemotherapeutic efficacy (Talib et al., 2021). Tumor heterogeneity, particularly intra-tumoral genetic and phenotypic diversity, significantly contributes to therapeutic failure by enabling dynamic adaptation to treatment pressures (Fernandez et al., 2018). To overcome this complexity, combinatorial strategies incorporating pathway inhibitors have emerged as a paradigm shift in overcoming drug resistance. These inhibitors specifically target oncogenic signaling cascades to effectively attenuate chemotherapy resistance. They achieve this by simultaneously blocking adaptive survival pathways that cancer cells activate under treatment stress and disrupting pro-survival signaling networks that confer drug tolerance, thereby enhancing therapeutic efficacy (Neel and Bivona, 2017). For instance, Sui et al. (2010) demonstrated that pharmacological inhibition of the estrogen receptor pathway overcomes resistance to vinca alkaloids in breast cancer models, highlighting the potential for targeted co-therapies. Such synergistic interactions between pathway modulators and conventional chemotherapeutics not only enhance treatment efficacy but also provide mechanistic insights into resistance antagonism, paving the way for precision oncology approaches (Song et al., 2022).

Traditional pathway gene annotations are predominantly curated through subjective methodologies, wherein constituent genes exhibit functionally coordinated interactions to execute specific biological processes or pathways (Chicco and Agapito, 2022). These gene sets are compiled from heterogeneous sources, including curated pathway databases (Kanehisa and Goto, 2000), biomedical literature mining (Venet et al., 2011), and domain expert knowledge systems (Zhang et al., 2011). Classification criteria for pathway membership gene sets (PMGS) vary widely, ranging from genes co-participating in canonical signaling cascades to those sharing homologous functional annotations (Li et al., 2017). Notably, when biological systems are perturbed by pharmacological agents, a subset of genes demonstrating robust transcriptional reprogramming, termed “response signatures”, emerges as critical biomarkers for deciphering context-specific therapeutic mechanisms (Atkinson et al., 2019). This conceptual divergence between PMGS and response signatures lies in their foundational definitions: the former is pre-defined through static pathway frameworks, while the latter arises dynamically from observed transcriptional adaptations to pathological stimuli. Such mechanistic dichotomy underscores the necessity for a dynamically informed model that integrates pathway response signatures into existing pathway frameworks. This approach further leverages treatment-perturbed pharmacodynamic profiles to refine therapeutic target prioritization, effectively decoupling analysis from static gene set definitions.

A substantial body of publicly available perturbation experiments has been systematically curated to identify treatment-responsive genes (Rydenfelt et al., 2020). These studies employ differential expression profiling to compare perturbed and control conditions, enabling the detection of genes exhibiting consistent regulatory responses to specific interventions (Parikh et al., 2010). Functionally, these genes are characterized by their tight association with perturbed biological pathways. This collective set is defined as pathway-responsive gene sets (PRGS).

Emerging evidence underscores the critical association between post-treatment response gene signatures and therapeutic outcomes, including both drug sensitivity and acquired resistance. Functional interrogation of these genes has revealed their dual capacity to delineate pathway perturbations driving resistance while simultaneously identifying actionable therapeutic targets. For instance, Cai et al. (2023) demonstrated in breast cancer models that palbociclib (a CDK4/6 inhibitor) induced significant upregulation of cyclin-dependent kinase regulators (CDK4, Cyclin D1, Cyclin E1), which mechanistically activated the PI3K/AKT/mTOR axis through a feedback loop involving Cyclin D1 overexpression. Notably, co-administration of PI3K/mTOR/AKT inhibitors disrupted this circuit by reducing Cyclin D1 levels, effectively attenuating palbociclib resistance. Similarly, Zhao et al. (2024) elucidated a gefitinib resistance mechanism in non-small cell lung cancer (NSCLC) driven by S6K1 hyperactivation, a process mediated through the ELK1/mTOR/S6K1 signaling axis. Pharmacological inhibition of S6K1 using the selective inhibitor PF-4708671 not only attenuated gefitinib resistance in resistant models, while simultaneously enhancing anti-tumor efficacy through targeted pathway suppression, thereby establishing S6K1 as a master regulator of adaptive resistance.

In SPEED (Parikh, Klinger, Xia, Marto and Bluthgen, 2010) and SPEED2 (Rydenfelt, Klinger, Klunemann and Bluthgen, 2020), numerous experimental perturbation datasets have been utilized. These datasets could be employed to construct the PRGS, which maps treatment-induced gene expression to drug sensitivity/resistance profiles, uncovering the mechanistic relationship between pathway dysregulation and therapeutic outcomes.

The evolution of enrichment analysis methodologies has been driven by the need to address limitations in detecting biologically meaningful gene sets. Over-representation analysis (ORA), a foundational method in enrichment analysis, employs the Hypergeometric test to assess statistical over-representation of differentially expressed genes (DEGs) within predefined pathways (Liu et al., 2022). While computationally efficient, ORA exhibits inherent sensitivity biases, particularly failing to capture genes with subtle expression changes that contribute meaningfully to biological processes (Peng et al., 2023). To overcome these constraints, Gene Set Enrichment Analysis (GSEA) introduced a paradigm shift by evaluating pathway enrichment through Kolmogorov-Smirnov statistics applied to ranked gene lists, prioritizing collective pathway behavior over individual gene metrics (Cai et al., 2022; Hou et al., 2020). Concurrently, the Bates test expanded the analytical framework by modeling gene set mean rank shifts through Bates distribution approximations, enabling robust detection of moderate but biologically relevant pathway dysregulation (Rydenfelt, Klinger, Klunemann and Bluthgen, 2020). This tripartite methodological framework, ORA for hypothesis-driven single-gene prioritization, GSEA for systems-level pathway interrogation, and the Bates test for quantitative ranking analysis, was systematically applied to PRGS, achieving synergistic validation of pathway-dysregulation hypotheses through complementary analytical perspectives.

In this study, we aim to identify potential therapeutic agents to address chemoresistance in cancer treatment. We constructed PRGS by systematically integrating multiple experimental perturbation datasets. Subsequent application of enrichment analysis methods (GSEA-like, Hypergeometric test-based, Bates test-based) to PRGS in GDSC data revealed candidate drugs with the capacity to reduce cancer cell resistance. Finally, experimental validation was conducted to assess the synergistic therapeutic effects of these drug combinations.

2 Materials and methods

The overall process of data preparation, differential expression analysis, dataset construction, enrichment analysis, and cellular experimental validation in this study is shown in the flowchart (Figure 1).

Figure 1
Analysis flowchart and results depicting differentially expressed genes (DEGs) between drug-resistant and drug-sensitive samples. Panels include pathway expression dataset tables, pathway-responsive gene sets, and validation in T47D cell line. Graphs show enrichment framework, pathway connections with drugs and inhibitors, and significantly dysregulated pathways analyzed through different test methods, like GSEA-like and Bates test-based. Sankey diagrams visualize potential drug discovery pathways.

Figure 1. Comprehensive workflow for identifying chemoresistance-associated pathways and therapeutic candidates through PRGS-based enrichment analysis frameworks. (A) DEGs were systematically identified between drug-sensitive and drug-resistant groups using GDSC data. (B) PRGS was constructed by integrating multiple experimental perturbation datasets. (C) A comparative analysis of three PRGS-based enrichment methods (GSEA-like, Bates test-based, and Hypergeometric test-based) was performed to identify pathways associated with chemotherapy resistance. (D) Significant pathways linked to chemoresistance were identified between drug-resistant and drug-sensitive groups. (E) Potential therapeutic candidates capable of mitigating chemoresistance were discovered based on these pathways. (F) The synergistic effect of the bortezomib-bleomycin combination was experimentally validated in T47D cells.

2.1 Data sources

2.1.1 Pathway perturbation-response gene expression datasets

The experimental perturbation-response datasets were obtained from the study previously published by Rydenfelt et al. (2020), containing 15 pathways, 534 experimental datasets, and 29,133 genes. Each pathway contains multiple experimental datasets.

2.1.2 Drug sensitivity data in breast cancer cell line

We curated a panel of 78 chemotherapeutic agents approved by the American Cancer Society (ACS, https://www.cancer.org/) (Narayanan and Honavar, 2017). Chemosensitivity profiles were retrieved from the Genomics of Drug Sensitivity in Cancer (GDSC) database for breast cancer cell lines, with complete datasets available for 19 agents (Yang et al., 2013). To ensure statistical robustness, compounds exhibiting extreme response bias (defined as >90% of cell lines showing uniform resistance or sensitivity) were excluded, eliminating seven agents. This selection yielded 12 drugs with balanced resistance/sensitivity distributions for quantitative analysis: 5-fluorouracil; bleomycin; docetaxel; doxorubicin; epirubicin; etoposide; gemcitabine; paclitaxel; vinblastine; vincristine; vinorelbine; and vorinostat.

The GDSC repository provided chemosensitivity profiles across 51 breast cancer cell lines, with individual drug sensitivity data coverage ranging from 38 to 51 cell lines per agent among the 12 selected agents. Drug sensitivity was quantified using the area under the concentration-response curve (AUC) metric. Drug-sensitive and drug-resistant groups were defined using predefined AUC thresholds: drug-sensitive groups by AUC<0.8 and drug-resistant groups by AUC≥0.8 (Callari et al., 2023). The sensitivity of each chemotherapeutic agent was ranked in descending order based on AUC values. Subsequently, drug-sensitive and drug-resistant groups were identified through threshold-based selection, prioritizing the lowest eight sensitive and highest 8–12 resistant cell lines.

2.1.3 Gene expression in breast cancer cell line

Gene expression profiling data comprising 17,737 genes across 1,018 cancer cell lines were retrieved from the GDSC database. Subsequently, gene expression data were extracted from breast cancer cell lines for drug sensitivity analysis, encompassing 17,419 genes and 51 distinct breast cancer cell lines.

2.2 Differentially expressed genes between drug-sensitive and drug-resistant groups

Breast cancer datasets were categorized into drug-sensitive and drug-resistant groups based on responses to 12 chemotherapeutic agents (5-fluorouracil, paclitaxel, etc.). DEGs analysis was conducted using independent samples t-tests, with statistical significance defined as p-value < 0.01.

2.3 Pathway-responsive gene sets construction

For each experimental perturbation dataset, genes with expression below the median were excluded. When a gene was expressed in multiple datasets, the z-value corresponding to the highest expression was selected; conversely, for genes expressed in a single dataset, their observed z-values were retained. An adjusted z-value reflecting gene expression changes specific to each perturbation type was calculated using Formula 1:

zadjusted=zraw*signtype(1)

where zadjusted represents the adjusted z-value for the gene, zraw represents the raw z-value for gene expression, the sign represents the signum function (returning 1 for input > 0, 0 for input = 0, or −1 for input<0), and type represents the perturbation types (assigning 1 for activation and −1 for inhibition).

For each pathway, z-transformed values were derived by applying Gaussian transformation to all adjusted z-values within each dataset (Hanzelmann et al., 2013) (Supplementary Figure S1).

Subsequently, the mean of these z-transformed values across all experiments was calculated for each gene (hereafter termed z-score), with genes expressed in fewer than two datasets excluded from the analysis.

Pathway-related gene signatures (PRGS) were identified using a stringent statistical threshold |z-score|>2.58, p-value < 0.01).

2.4 PRGS-based enrichment analysis

We conducted systematic pathway enrichment analysis using the PRGS framework to identify dysregulated pathways between chemoresistant and chemosensitive breast cancer subtypes. Pathway significance was determined with p-values (p < 0.01). Three methodologies were utilized in the analysis: GSEA-like, Bates test-based, and Hypergeometric test-based analyses, each implemented with distinct input variable configurations to identify relevant pathways. Within the PRGS framework, variables were categorized into continuous and discrete types based on their characteristics and roles in pathway analysis. Continuous PRGS were defined as genes within pathways that underwent computational processing without statistical thresholding, retaining quantitative measurements of gene expression changes. Discrete PRGS were defined as signature genes meeting stringent statistical criteria and closely associated with pathway regulation. Total genes were derived from the gene expression dataset following independent t-test analyses, retaining their processed were obtained by applying statistical thresholds to these total genes, identifying those significantly associated with drug sensitivity or resistance.

Input variables for the three PRGS-based enrichment analysis methods were as follows: discrete PRGS and total genes were used for the GSEA-like methodology; continuous PRGS and DEGs were employed for the Bates test-based methodology; DEGs and discrete PRGS were used for the Hypergeometric test-based methodology (Figure 2).

Figure 2
A three-panel diagram compares gene expression analysis methods. Panel A shows estrogen-responsive signatures ranked by t-statistic, with results visualized in a running enrichment score plot. Panel B presents similar signatures identified as differentially expressed genes, depicted with a Venn diagram showing overlap with pathway-responsive gene sets. Panel C orders these signatures by adjusted z-score, highlighting their differential expression in various pathways, displayed in a heatmap with upregulated and downregulated distinctions.

Figure 2. Systematic workflow for PRGS-based enrichment analysis methodologies. (A) In the GSEA-like methodology, genes are ranked using t-statistics, and PRGS enrichment is assessed within the generated ranked gene list. (B) The overlap between DEGs and PRGS is evaluated by the Hypergeometric test-based method, and a p-value is calculated to determine the significance of pathway enrichment. (C) Within the Bates test-based method, genes are ranked using adjusted z-scores, and PRGS enrichment is evaluated within the resulting ranked list.

2.4.1 GSEA-like methodology

In GSEA-like method (Figure 2A), genes were ranked based on t-statistics, and the enrichment score (ES) for each pathway was calculated. The distribution of the ES for the jth gene g in the ith position of PRGS was calculated using a modified Kolmogorov-Smirnov statistic. The ES was evaluated by calculating the maximum absolute deviation between Phit and Pmiss, as defined by Formula 24 (Marczyk et al., 2021):

PhitS,i=gjSjirjpNRwhereNR=gjSrjp(2)

Where S represents the PRGS, r is the ranked score, p is a weighting parameter, N is the number of total genes in the ranked list, and NR is calculated by summing the ranks of genes within PRGS.

PmissS,i=gjSji1NNH(3)

Where NH is the number of genes in the PRGS.

ES=PhitPmiss(4)

Where Phit represents the weighted sum of genes within PRGS, Pmiss represents the weighted sum of genes outside PRGS. The ES increases when a gene is part of a gene set and decreases otherwise. This process was repeated for PRGS and total genes in each group.

2.4.2 Hypergeometric test-based methodology

The Hypergeometric test-based method was employed to ascertain the statistical significance of a given number of successes (k) in a sample size (n) observed within a larger total size (N), under the assumption that the total population encompasses a total of (M) occurrences of the specific characteristic being analyzed (Elia Venanzi et al., 2024). In this study, PRGS-based Hypergeometric test-based enrichment analysis was conducted on DEGs to identify significant over-expression within the PRGS (Figure 2B). The probability of observing at least k overlapping genes between DEGs and PRGS was calculated by Formula 5:

PXk=1i=0k-1Mi*NMniNn(5)

Where N represents the total number of genes, n is the number of DEGs, M is the number of overlapping genes between total genes and PRGS, and k is the number of overlapping genes between DEGs and PRGS. A smaller p-value indicates a higher likelihood of DEG enrichment within the PRGS.

2.4.3 Bates test-based methodology

The significance of pathway rankings under the null hypothesis was evaluated through the Bates test-based methodology (Figure 2C), which employs a null distribution modeled as the arithmetic mean of N independent uniform variables spanning the interval [−1, 1]. For each biological pathway, differential expression z-values were computed as the ratio of experimental fold changes (stimulated condition vs. control condition) to LOESS regression-derived standard deviations, effectively normalizing expression variation across experimental datasets. These raw z-values underwent Gaussian transformation to address distributional skewness, yielding z-transformed metrics (z-transformed value) with improved normality characteristics. Subsequent analytical steps included: (i) calculation of pathway mean z-transformed values across all N experiments, (ii) Quantile normalization to mitigate batch effects, and (iii) iterative sorting-rescaling procedures that transformed normalized means into final z-scores bounded within [−1,1]. Gene ranking within each PRGS was ultimately determined by these standardized z-scores, as defined in Formula 6:

rg=zm+1e50;ifzmgUa+Ub2zm1e50;ifzmg<Ua+Ub2(6)

Where zm is the mean z-value of gene, Ua = 1 and Ub = −1 represent the upper and lower bounds of the standardized interval. The second rank was performed using Formula 7, scaling the initial rank to the interval [−1,1].

rg=1rg;ifrg01rg;ifrg>0(7)

The r'(g) values of all genes are sorted in ascending order, and each gene is assigned a unique rank from 1 to n based on its position, resulting in the third rank r''(g). Subsequently, these r''(g) values are scaled to the interval [−1,1], yielding the fourth rank r'''(g). Within the context of a designated set of DEGs, overlapping genes were pinpointed through the intersection of PRGS and DEGs. Subsequently, the corresponding r'''(g) values for these overlapping genes were retrieved. The mean rank R(G) was computed by averaging these r'''(g) values, as delineated in Formula 8:

RG=i=1nrgn(8)

Where n is the number of overlapping genes, r'''(g) is the value of the fourth rank. The Bates p-value was calculated based on the mean rank and the number of overlapping genes using Formula 9:

Px=batesCDFx,n;ifxa+b21batesCDFa+bx2,n;ifx>a+b2(9)

Where x is the R(G), n is the number of overlapping genes, a corresponds to Ua and is assigned a value of 1, while b corresponds to Ub and is assigned a value of −1.

2.5 Discovery of therapeutic inhibitors

Based on the statistical significance determined by PRGS-based analysis, the top 1-2 pathways most strongly associated with chemoresistance were prioritized as therapeutic targets for each chemotherapeutic agent. Targeted inhibitors were identified by screening the Selleck platform (Wojtaszek et al., 2019) (https://www.selleck.cn/index.html) and PubMed, adhering to the following criteria: (i) inhibitors that directly target pathways identified by PRGS; (ii) among pathways with multiple targeted inhibitors, priority should be given to drugs that have been validated through in vitro and in vivo experiments, employed in numerous studies and demonstrate broad applicability; (iii) drugs that modulate specific pathways indirectly were also considered; (iv) relevant inhibitors were sourced from PubMed-indexed preclinical research; (v) extracts and plant-derived compounds were excluded.

2.6 In vitro validation

The predicted synergistic interaction between bortezomib and bleomycin was systematically validated through a two-phase pharmacodynamic evaluation in STR-authenticated T47D breast cancer cells maintained under standardized culture conditions (37 °C, 5% CO2). In the first phase, single-agent dose-response profiling was performed to determine the optimal concentrations for subsequent combination experiments. Bortezomib (0.5–40 nM) and bleomycin (50–400 nM) were tested in eight technical replicates per concentration. Dose-response curves revealed that 200 nM bleomycin corresponded to the optimal IC50 threshold concentration, establishing it as the reference concentration for combination studies. In the second phase, combination experiments were conducted by exposing cells to a fixed concentration of 200 nM bleomycin combined with a gradient of bortezomib (0.5–40 nM). Each treatment arm was performed in eight technical replicates to ensure robustness of the data. After 48 h of co-culture, cell viability was quantified using the Cell Counting Kit-8 (CCK-8) via absorbance measurement at 450 nm. Synergistic efficacy was quantitatively assessed using the IDAcomboScore method as described by Ling and Huang (2020). Full details of all experimental protocols are available in Supplementary Material.

3 Results

3.1 Pathway-responsive gene sets

Heterogeneous data distribution across cross-experimental gene expression datasets poses significant challenges to robust analytical outcomes, particularly when dealing with heavy-tailed distributions that violate normality assumptions. To address this issue, we systematically evaluated three normalization strategies to determine the most suitable preprocessing approach.

The most appropriate method was selected from Gaussian normalization, Logarithmic transformation, and Quantile normalization. Using the GSE13837 dataset as an illustrative example, 70% of the data points were initially distributed within the interval [−20, 20], while the remaining 30% exhibited a heavy-tailed distribution outside this interval (Supplementary Figure S2). After normalization, Gaussian normalization markedly diminished distributional heterogeneity, concentrating 90% of the data points within the narrower interval of [−3, 3], while the ratio of negative to positive values was balanced, changing from 58.9:41.1 to 49.6:50.4 (Supplementary Figure S3).

In contrast, the Logarithmic transformation resulted in data clustering predominantly within the range [−5, 5], while the overall distribution remained skewed. Notably, the ratio of negative to positive values shifted markedly to 25.1:74.9.

Although the distribution resulting from Quantile normalization approximates a normal distribution within the interval [−15, 15], a limited number of extreme values persist beyond these boundaries. Similar patterns were observed in the other two datasets (Supplementary Figures S2 and S3).

Overall, Gaussian normalization emerged as the optimal method. It maps data onto a standard normal distribution, effectively compressing extreme values and improving distributional balance.

Substantial inter-dataset distributional disparities were observed, and we systematically applied Gaussian normalization to all pathway-derived datasets, achieving scale consistency across expression datasets (Hai et al., 2023).

A systematic analysis of these post-normalization datasets spanning multiple pathways identified genes with consistent regulation patterns across experimental conditions, enabling the construction of PRGS (Supplementary Figure S4A). Differences in the number of datasets across pathways were observed, with IL-1 (9 datasets) and Wnt (12 datasets) having relatively small numbers of experimental data, while JAK-STAT (66 datasets) and TGFa (44 datasets) had relatively large numbers. This variation may stem from differences in experimental design, the diversity of available perturbation reagents, or the maturity of research on different pathways. Quantitative comparison revealed that the IL-1 pathway exhibited the highest number of responsive signatures (410), which may be attributed to its role as a key pro-inflammatory factor that triggers extensive transcriptional changes (Supplementary Figure S4B). These changes are associated with numerous biological processes, including immune responses, cell proliferation, and apoptosis, where significant gene expression alterations occur (Xiao et al., 2022). In contrast, the JAK-STAT pathway exhibited the fewest responsive genes, with only 91 identified. This may indicate that fewer genes meet the criteria for consistent response under the experimental conditions used.

Despite variations in the number of responsive genes, they provide a robust foundation for subsequent functional analysis, enabling more effective identification of pathway dysregulation in specific biological samples. For instance, in a study by Gregory et al., their analysis of dysregulated genes in response to imatinib therapy identified the Wnt pathway as a therapeutic target (Gregory et al., 2010).

3.2 Differential expression of drug-sensitive genes across drug groups

DEGs between drug-sensitive and resistant breast cancer cell groups across 12 chemotherapeutic agents were identified through differential expression analysis. To focus on genes with significant directional expression changes, we applied an additional filter (|t-statistic|>2.8) to prioritize biologically impactful candidates. This analysis revealed substantial heterogeneity in t-statistic distributions across agents (Table 1; Supplementary Figure S5A), indicating that chemotherapy resistance-associated pathway dysregulation exhibits distinct molecular signatures under different drug regimens.

Table 1
www.frontiersin.org

Table 1. The number of DEGs of drug-sensitive genes.

In cancer, downregulated tumorous suppressor genes have been linked to chemotherapy resistance (Tallen et al., 2003). Half of the DEGs in the some analyzed group showed a tendency toward downregulation (i.e., a greater number of downregulated genes than upregulated genes). Conversely, the remaining DEGs of sample groups exhibited an opposite expression pattern, suggesting potential differences in their chemoresistance mechanisms. Additionally, significant expression changes of some DEGs were only observed in the paclitaxel- or gemcitabine-treatment group, with minimal overlap between these groups (Supplementary Figure S5B). These findings indicated a potential association between these genes and resistance to specific pharmaceutical agents.

3.3 Evaluation of different enrichment methods

Three enrichment analysis methodologies were systematically applied to the PRGS framework using drug response-associated DEGs. These genes were analyzed to compare drug-sensitive and -resistant subgroups within the GDSC dataset, facilitating the identification of pathways driving chemoresistance. The GSEA-like methodology demonstrated superior performance, identifying eight key pathways (p-value < 0.01) with strong activation of IL-1 and TNFα pathways in multidrug-resistant groups (Figure 3A). Notably, both of which were inflammatory pathways. Functional interrogation further uncovered shared regulatory nodes (e.g., NF-kB activation) between these two pathways, mechanistically linking their coactivation to sustained chemoresistant phenotypes (Feng et al., 2025; Gao et al., 2024).

Figure 3
Circular diagrams (A, B, C) show connections between cytokines, pathways, and drugs like Vorinostat, Vincristine, and Docetaxel. Diagram D is a heatmap displaying p-values for different tests (GSEA-like, Bates, Hypergeometric, PMGS-based GSEA-like) with color gradients representing significance levels in pathways such as PPAR, Notch, and Hypoxia.

Figure 3. Identification of significant pathways by three enrichment analyses. (A–C) Significant pathways were identified using three enrichment analysis methods: GSEA-like, Hypergeometric test-based, and Bates test-based. A chord diagram illustrates the complex relationships between pathways and drug-sensitive/resistant groups, with the width of each band reflecting the strength of these relationships. (D) A comparison of enrichment results across different drug-sensitive/resistant groups is shown. The x-axis represents drug-sensitive/resistant samples, while the y-axis represents 15 pathways. Color intensity indicates the significance of enrichment p-values, with grey indicating pathways not identified by the method. This visualization highlights the performance differences among the three methodologies.

Three chemoresistance-associated pathways were identified through Hypergeometric test-based methodology across two drug-sensitive/resistant subgroup comparisons (Figure 3B). The statistical significance of IL-1 and TNFα pathways was consistently validated by both Hypergeometric test-based and GSEA-like methodologies. Gene-sharing limitations between PRGS and DEGs resulted in undetectable key pathways for 5-fluorouracil, bleomycin, and docetaxel groups (Figure 3D).

Seven chemoresistance-associated pathways were identified through Bates test-based analysis, with Hippo, PPAR, and JAK-STAT pathways showing consistent activation across multiple subgroups (Figure 3C). Functional validation confirmed these pathways as master regulators of chemoresistance acquisition in preclinical models (Nascimento et al., 2017; Zeng and Dong, 2021; Zhang et al., 2024). The statistical comparison revealed that the Bates test-based methodology identified significantly fewer chemoresistance-associated pathways compared to the GSEA-like methodology (Figure 3D).

The comparative analysis demonstrated that pathways identified through the GSEA-like methodology achieved superior statistical significance (p < 0.0001). Significant overlap was observed between pathways identified by the GSEA-like methodology and those detected by alternative approaches, with the GSEA-like methodology identifying the highest number of pathways across all evaluated methods. These findings supported the conclusion that the GSEA-like methodology represents the optimal approach for pathway identification.

3.4 Comparing differences between PRGS and PMGS

Gene sets for ten pathways were retrieved from the KEGG database to evaluate differences between PRGS and PMGS. The Jaccard index and the number of overlaps for each pathway were calculated (Supplementary Figure S6; Supplementary Table S1). The results show extremely low overlap between PRGS and PMGS across all pathways, with an average Jaccard index of only 0.012 (range: 0–0.036).

From a biological perspective, this limited overlap stems from their distinct design objectives. PMGS includes genes that are annotated as components of well-known canonical pathways and signaling cascades (Gerling et al., 2013). These encompass pathway receptors, kinases, transcription factors and co-regulators whose functions remain unaffected by external perturbations (Rashidiani et al., 2024). To illustrate this point, consider the estrogen pathway. In estrogen receptor-positive (ER+) breast cancers that HER2 and EGFR are overexpressed, downstream signaling components can be activated. This activation subsequently enhances the activity of the estrogen receptor and its co-activator AIB1, thereby facilitating the estrogen agonistic activity of tamoxifen in breast cancer (Rahem et al., 2020). This phenomenon results from dysregulated pathway activity due to aberrant gene function within the PMGS. In contrast, PRGS comprises genes whose expression levels are markedly upregulated or downregulated in response to external perturbations affecting the pathway (Wang, Karikomi, et al., 2019). Within the estrogen pathway, 17β-estradiol (E2) acts as a perturbing ligand that binds to estrogen receptor α (ERα) in the cytoplasm. This interaction induces a conformational change in ERα, prompting its translocation to the nucleus, where it modulates gene transcription by binding to estrogen response elements (EREs) located in the promoter regions of target genes (Krolick and Shi, 2022). These target genes, which constitute the PRGS, are not core pathway components but rather represent downstream “response products” elicited following pathway activation.

The GSEA-like methodology demonstrated superior performance when applied to PRGS using GDSC data (p < 0.0001). Subsequent implementation of the GSEA-like approach on PMGS identified four key pathways (Supplementary Figure S7). The Insulin and JAK-STAT pathways have been previously reported to be associated with chemoresistance (Han et al., 2021; Macaulay et al., 2013). The GSEA-like methodology is based on differential expression rankings, which are inherently supported by the continuous gene expression values of PRGS. In contrast, PMGS-based analysis produced statistically significant but less robust enrichment outcomes (p < 0.01), as its pathway annotations lack quantitative expression data required for optimal GSEA-like implementation (Figure 3D). This discrepancy arises because PRGS captures dynamic pathway regulation through expression magnitude, while PMGS only reflects static gene set membership.

3.5 Discovering potential drugs to overcome chemoresistance in breast cancer

Eight key pathways associated with chemoresistance in breast cancer were identified using the PRGS-based GSEA-like approach. These pathways are dysregulated in response to specific chemotherapeutic agents, and can be targeted for therapeutic intervention with pathway-specific inhibitors (Figure 4). In line with the screening criteria, the selected pathway inhibitors were obtained from the Selleck platform and relevant literature indexed in PubMed.

Figure 4
Sankey diagram illustrating relationships among chemotherapy drugs, pathways, and pathway inhibitors. Drugs like 5-Fluorouracil, Vinblastine, and Paclitaxel connect to pathways such as Hippo, Notch, and TNFa, which further connect to inhibitors like CAY10585, Verteporfin, and Tamoxifen. Each pathway and connection is color-coded to indicate the relationships.

Figure 4. Discovery of drug candidates to overcome chemoresistance in cancer treatment. The Sankey diagram illustrates the relationships between chemotherapeutic drugs, dysregulated pathways, and pathway inhibitors. Nodes on the left represent chemotherapeutic agents for cancer treatment, central nodes correspond to dysregulated pathways associated with these drugs, and nodes on the right indicate pathway inhibitors capable of modulating these pathways. Distinct colors and node shapes correspond to various chemotherapy agents, dysregulated pathways, and inhibitors, respectively. The connecting lines illustrate the directional flow between these three elements.

The estrogen pathway was found to be dysregulated in paclitaxel-resistant subgroups. Tamoxifen and fulvestrant, both explicitly classified as “estrogen receptor antagonists” in the Selleck database, function by inhibiting estrogen activity and suppressing the proliferation of breast cancer cells. Previous studies have shown that fulvestrant can make estrogen receptor-negative breast tumors more susceptible to chemotherapy, whereas tamoxifen can mitigate paclitaxel resistance mediated by ERα-positive signaling (Jiang et al., 2014; Nzegwu et al., 2021). Although both clomifene citrate and tamoxifen are estrogen receptor modulators, clomifene citrate has limited applicability and lacks sufficient experimental validation in cancer models. In contrast, tamoxifen is supported by Phase II clinical trial data and exhibits greater practical utility, with multiple in vitro studies confirming its capacity to induce apoptosis in breast cancer cells. Consequently, tamoxifen and fulvestrant are viable inhibitors for the direct targeting of the estrogen pathway.

Bortezomib, while primarily a proteasome inhibitor, also inhibits NF-kB activity. Its combined administration with IL-1 receptor antagonists has been shown to reduce tumor formation and growth in vivo (McLoed et al., 2016). In patients with rheumatoid arthritis, bortezomib suppresses the release of pro-inflammatory cytokines TNFα and IL-1β (Maseda et al., 2014). Therefore, it serves as an indirect inhibitor of the IL-1 pathway.

In addition, drugs that had previously been validated as inhibitors of specific pathways were considered. ST-162 and ST-168 are small-molecule bifunctional inhibitors targeting both the MEK and PI3K pathways, demonstrating dose-dependent suppression of MEK and PI3K signal transduction (Smith et al., 2018). They are being developed as novel anti-tumor agents, but are not currently documented in the Selleck platform. Thus, they represent potential inhibitors for targeting the MAPK+ PI3K pathways.

Promising drug combinations targeting resistant cell lines were identified through this workflow (Supplementary Table S2). A notable example is the HCC1428 paclitaxel-resistant cells (AUC = 0.95), exhibited significant estrogen pathway dysregulation (Yang et al., 2013). Fulvestrant was determined to be an inhibitor of the estrogen pathway (Wakeling et al., 1991). Therefore, they were considered as a potential therapeutic combination. These findings provide diverse combination therapy strategies to overcome chemoresistance in breast cancer.

3.6 Bleomycin-bortezomib synergistic effects validation in T47D cells

Both IL-1 and TNFα pathways were identified as key pathways in both drug-sensitive and drug-resistant subgroups through GSEA-like and Hypergeometric test-based analyses.

Studies have shown that the IL-1 pathway increases drug resistance by activating and reinforcing the release of pro-inflammatory cytokines (Stanam et al., 2016). Moreover, in breast cancer cells exhibiting high sensitivity to IL-1β, stimulation with IL-1β markedly induces the upregulation of BIRC3 expression, consequently conferring resistance to doxorubicin treatment (Rho et al., 2022). In alignment with these findings, BIRC3 is identified as a responsive gene within the IL-1-responsive gene set, thereby strongly suggesting its role as a mediator of IL-1-induced chemoresistance.

Overexpression of TNFα can fragment tumor DNA, thereby contributing to resistance to the cytotoxic effects of chemotherapeutic agents (Adrian et al., 2023). TNFα can also activate the NF-kB signalling pathway and promote the expression of CXCL1/2, which leads to the amplification of the CXCL1/2-S100A8/9 loop and induces chemotherapy resistance (Cruceriu et al., 2020). Consistently, the presence of CXCL1 and CXCL2 genes associated with the TNFα-responsive gene set was observed, suggesting that this mechanism likely plays a significant role in TNFα-mediated drug resistance.

The T47D cell line represents the luminal A subtype of breast cancer, characterized by positivity for estrogen receptor (ER+) and progesterone receptor (PR+) and human epidermal growth factor receptor 2 (HER2-) (Ahmed et al., 2022). This subtype constitutes approximately 70% of all breast cancer cases (Fiorillo et al., 2020). The proliferation of T47D cell was inhibited upon exposure to both native and recombinant human interleukin-1 (IL-1) isoforms, specifically the α and β variants (Gaffney and Tsai, 1986). Importantly, analysis using PRGS has identified the IL-1 pathway as a key therapeutic target in the bleomycin response subgroup. Bortezomib exhibits anti-inflammatory activity and suppresses IL-1 (Mao et al., 2012). These findings suggest the potential for synergistic benefits in combining bleomycin with IL-1 inhibitors for use in patients who do not respond to bleomycin. Consequently, bortezomib was selected for combination treatment with bleomycin in T47D cell.

Under optimal drug concentrations, this combination reduced cellular viability to 20% (Figure 5), with bleomycin and bortezomib demonstrating IC50 values of 200 nmol/L and 40 nmol/L, respectively. An IDAcomboscore of 0.014 (threshold≥0.004 for synergy) confirmed significant interaction between bleomycin and bortezomib.

Figure 5
Box plot comparing cell activity for three groups: Bleomycin, Bortezomib, and their combination. Bleomycin shows lower median activity, Bortezomib has a wider range, and the combination group exhibits intermediate activity.

Figure 5. The effect of bleomycin, bortezomib monotherapy and their combination on T47D cell viability. The x-axis represents the different treatment groups (bleomycin, bortezomib, combination therapy), while the y-axis represents cell viability.

4 Discussion

Despite the increasing application of gene set analysis, systematic investigations into gene set selection criteria remain limited. The impact of distinct gene set definitions on analytical outcomes remains underexplored. Two gene set categories, PRGS and PMGS, were compared in this study. PMGS comprise curated gene sets associated with biological processes, chromosomal locations, disease associations, or pathway membership (Hermida et al., 2013), reflecting their functional annotation purposes. Conversely, PRGS capture pathway-specific expression changes under perturbation conditions, maintaining causal relationships between gene expression and pathway activity. Significant differences in gene composition were observed between PMGS and PRGS, with limited overlap indicating divergent design objectives. While PRGS genes may overlap with PMGS, reciprocal inclusion was not observed. Consequently, PRGS were prioritized to characterize chemotherapy resistance-associated expression patterns.

Three enrichment analysis methodologies based on PRGS were implemented. The GSEA-like approach demonstrated superior performance, identifying a greater number of statistically significant pathways compared to alternative methods. This superiority was attributed to its emphasis on aggregate gene set trends within ranked lists. Notably, predefined gene lists exhibited substantial differential expression levels, contributing to substantial pathway enrichment. The Bates test-based method involves multiple sorting of gene z-values, focusing on the average ranking of genes within a given gene set. Conversely, Hypergeometric test-based analysis was constrained by limited gene overlap between PRGS and DEGs, with overlapping gene quantity prioritized over expression magnitude. This methodological limitation resulted in the lowest number of detected pathways among all evaluated approaches.

Current strategies to address chemotherapeutic resistance include combination therapies, nanomedicine platforms, and gene therapeutic approaches. A self-delivery nanomedicine incorporating α-tocopherol succinate and doxorubicin was developed by Zheng et al. to overcome drug resistance through synergistic chemotherapy (Zheng et al., 2021). Similarly, CRISPR/Cas9-mediated gene editing was employed to disrupt uPAR expression, resulting in enhanced chemosensitivity of HCT 8/T cells (Wang et al., 2019). In this study, pathway inhibitors were combined with chemotherapeutic agents to target dysregulated pathways underlying chemoresistance. The methodology proposed herein exhibits broad applicability beyond chemotherapy resistance, extending to endocrine drug and other drug class-associated resistance mechanisms through combination with pathway inhibitors. As an illustrative example, Dong et al. demonstrated that PI3K/AKT/mTOR pathway inhibition restores estrogen receptor signaling, thereby overcoming resistance to endocrine treatments in preclinical models (Dong et al., 2021). Drug combination of bleomycin and bortezomib was validated in T47D cell line. However, more experiments need to be performed in the following studies.

In summary, this study highlights the importance of PRGS in characterizing chemotherapy resistance-associated expression patterns. By prioritizing PRGS over PMGS, we were able to capture dynamic pathway-specific expression changes and their causal relationships with pathway activity, providing a more accurate representation of biological mechanisms underlying chemoresistance. The GSEA-like approach, with its focus on aggregate gene set trends, outperformed other methodologies in identifying statistically significant pathways, underscoring its utility in pathway analysis. These findings not only advance our understanding of chemoresistance mechanisms but also offer a robust framework for developing targeted therapeutic strategies. The proposed method exhibits broad applicability beyond chemotherapy resistance, extending to other drug class-associated resistance mechanisms through the integration of pathway inhibitors. Thus it provides a foundation for future research aimed at overcoming treatment resistance and improving therapeutic outcomes in cancer patients.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used. Ethical approval was not required for the studies on animals in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

DF: Conceptualization, Formal Analysis, Visualization, Writing – original draft. JH: Investigation, Methodology, Writing – original draft. LL: Resources, Validation, Visualization, Writing – original draft. JC: Investigation, Methodology, Writing – original draft. XL: Investigation, Methodology, Writing – original draft. RZ: Formal Analysis, Visualization, Writing – original draft. HH: Formal Analysis, Visualization, Writing – original draft. TL: Formal Analysis, Visualization, Writing – original draft. XW: Formal Analysis, Visualization, Writing – original draft. XL: Resources, Writing – review and editing. LY: Resources, Writing – review and editing. BL: Resources, Writing – review and editing. JL: Conceptualization, Funding acquisition, Resources, Supervision, Writing – review and editing. LW: Resources, Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work has been supported by the Natural Science Foundation of Hainan Province (Nos. 824RC514, 821QN0894, 621MS041); National Natural Science Foundation of China (No. 32260155); Academic Enhancement Support Program of Hainan Medical University (Nos. XSTS2025048 and XSTX2025030).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbinf.2025.1661601/full#supplementary-material

References

Adrian, K., Ghaib, H., and Ali, I. (2023). Tumour necrosis factor-alpha levels as predictor factor on clinical response of anthracycline-based neoadjuvant chemotherapy in locally advance breast cancer patients: experimental research. Ann. Med. Surg. (Lond). Apr 85, 807–811. doi:10.1097/ms9.0000000000000424

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahmed, E. A., Alkuwayti, M. A., and Ibrahim, H. M. (2022). Atropine is a suppressor of epithelial-mesenchymal transition (EMT) that reduces stemness in drug-resistant breast cancer cells. Int. J. Mol. Sci. 30, 23. doi:10.3390/ijms23179849

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkinson, M. A., Roep, B. O., Posgai, A., Wheeler, D. C. S., and Peakman, M. (2019). The challenge of modulating beta-cell autoimmunity in type 1 diabetes. Lancet Diabetes Endocrinol. 7, 52–64. doi:10.1016/s2213-8587(18)30112-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, Z., Tang, B., Chen, L., and Lei, W. (2022). Mast cell marker gene signature in head and neck squamous cell carcinoma. BMC Cancer 22 (22), 577. doi:10.1186/s12885-022-09673-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, Z., Wang, J., Li, Y., Shi, Q., Jin, L., Li, S., et al. (2023). Overexpressed Cyclin D1 and CDK4 proteins are responsible for the resistance to CDK4/6 inhibitor in breast cancer that can be reversed by PI3K/mTOR inhibitors. Sci. China Life Sci. 66. 94–109. doi:10.1007/s11427-021-2140-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Callari, M., Sola, M., Magrin, C., Rinaldi, A., Bolis, M., Paganetti, P., et al. (2023). Cancer-specific association between Tau (MAPT) and cellular pathways, clinical outcome, and drug response. Sci. Data 10, 637. doi:10.1038/s41597-023-02543-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chicco, D., and Agapito, G. (2022). Nine quick tips for pathway enrichment analysis. PLoS Comput. Biol. 18, e1010348. doi:10.1371/journal.pcbi.1010348

PubMed Abstract | CrossRef Full Text | Google Scholar

Cruceriu, D., Baldasici, O., Balacescu, O., and Berindan-Neagoe, I. (2020). The dual role of tumor necrosis factor-alpha (TNF-alpha) in breast cancer: molecular insights and therapeutic approaches. Cell Oncol. (Dordr) 43, 1–18. doi:10.1007/s13402-019-00489-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, C., Wu, J., Chen, Y., Nie, J., and Chen, C. (2021). Activation of PI3K/AKT/mTOR pathway causes drug resistance in breast cancer. Front. Pharmacol. 12, 628690. doi:10.3389/fphar.2021.628690

PubMed Abstract | CrossRef Full Text | Google Scholar

Elia Venanzi, N. A., Basciu, A., Vargiu, A. V., Kiparissides, A., Dalby, P. A., and Dikicioglu, D. (2024). Machine learning integrating protein structure, sequence, and dynamics to predict the enzyme activity of bovine enterokinase variants. J. Chem. Inf. Model 64, 2681–2694.doi:10.1021/acs.jcim.3c00999

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, B., Zhao, D., Zhang, Z., Jia, R., Schuler, P. J., and Hess, J. (2025). Ligand-receptor interactions combined with histopathology for improved prognostic modeling in HPV-Negative head and neck squamous cell carcinoma. NPJ Precis. Oncol. 9 (9), 57. doi:10.1038/s41698-025-00844-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez, H. R., Gadre, S. M., Tan, M., Graham, G. T., Mosaoa, R., Ongkeko, M. S., et al. (2018). The mitochondrial citrate carrier, SLC25A1, drives stemness and therapy resistance in non-small cell lung cancer. Cell Death Differ. 25, 1239–1258. doi:10.1038/s41418-018-0101-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Fiorillo, M., Toth, F., Brindisi, M., Sotgia, F., and Lisanti, M. P. (2020). Deferiprone (DFP) targets cancer Stem Cell (CSC) propagation by inhibiting mitochondrial metabolism and inducing ROS production. Cells 9, 1529. doi:10.3390/cells9061529

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaffney, E. V., and Tsai, S. C. (1986). Lymphocyte-activating and growth-inhibitory activities for several sources of native and recombinant interleukin 1. Cancer Res. Aug 46, 3834–3837.

PubMed Abstract | Google Scholar

Gao, W., Sun, L., Gai, J., Cao, Y., and Zhang, S. (2024). Exploring the resistance mechanism of triple-negative breast cancer to paclitaxel through the scRNA-seq analysis. PLoS One 19, e0297260. doi:10.1371/journal.pone.0297260

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerling, M., Nousiainen, K., Hautaniemi, S., Kruger, S., Fritzsche, B., Homann, N., et al. (2013). Aneuploidy-associated gene expression signatures characterize malignant transformation in ulcerative colitis. Inflamm. Bowel Dis. 19, 691–703. doi:10.1097/mib.0b013e31827eeaa4

PubMed Abstract | CrossRef Full Text | Google Scholar

Gregory, M. A., Phang, T. L., Neviani, P., Alvarez-Calderon, F., Eide, C. A., O'Hare, T., et al. (2010). Wnt/Ca2+/NFAT signaling maintains survival of Ph+ leukemia cells upon inhibition of Bcr-Abl. Cancer Cell 18, 74–87. doi:10.1016/j.ccr.2010.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Hai, Y., Shi, S., and Qin, G.For The Alzheimer's Disease N (2023). Bayesian and influence function-based empirical likelihoods for inference of sensitivity to the early diseased stage in diagnostic tests. Biom J. 65, e2200021. doi:10.1002/bimj.202200021

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, J., Yun, J., Quan, M., Kang, W., Jung, J. G., Heo, W., et al. (2021). JAK2 regulates paclitaxel resistance in triple negative breast cancers. J. Mol. Med. Berl. 99, 1783–1795. doi:10.1007/s00109-021-02138-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hermida, L., Poussin, C., Stadler, M. B., Gubian, S., Sewer, A., Gaidatzis, D., et al. (2013). Confero: an integrated contrast data and gene set platform for computational analysis and biological interpretation of omics data. BMC Genomics 14, 514. doi:10.1186/1471-2164-14-514

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, L., Jiao, Y., Li, Y., Luo, Z., Zhang, X., Pan, G., et al. (2020). Low EIF2B5 expression predicts poor prognosis in ovarian cancer. Med. Baltim. 99, e18666. doi:10.1097/md.0000000000018666

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, D., Huang, Y., Han, N., Xu, M., Xu, L., Zhou, L., et al. (2014). Fulvestrant, a selective estrogen receptor down-regulator, sensitizes estrogen receptor negative breast tumors to chemotherapy. Cancer Lett. 346 (346), 292–299. doi:10.1016/j.canlet.2014.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi:10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Krolick, K. N., and Shi, H. (2022). Estrogenic action in stress-induced neuroendocrine regulation of energy homeostasis. Cells 11, 879. doi:10.3390/cells11050879

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Wang, X., Xiao, G., and Gazdar, A. (2017). Integrative gene set enrichment analysis utilizing isoform-specific expression. Genet. Epidemiol. 41, 498–510. doi:10.1002/gepi.22052

PubMed Abstract | CrossRef Full Text | Google Scholar

Ling, A., and Huang, R. S. (2020). Computationally predicting clinical drug combination efficacy with cancer cell line screens and independent drug action. Nat. Commun. 11, 5848. doi:10.1038/s41467-020-19563-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Yuan, M., Mitra, R., Zhou, X., Long, M., Lei, W., et al. (2022). CTpathway: a CrossTalk-based pathway enrichment analysis method for cancer research. Genome Med. 14, 118. doi:10.1186/s13073-022-01119-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Macaulay, V. M., Middleton, M. R., Protheroe, A. S., Tolcher, A., Dieras, V., Sessa, C., et al. (2013). Phase I study of humanized monoclonal antibody AVE1642 directed against the type 1 insulin-like growth factor receptor (IGF-1R), administered in combination with anticancer therapies to patients with advanced solid tumors. Ann. Oncol. Mar. 24, 784–791. doi:10.1093/annonc/mds511

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X., Pan, X., Peng, X., Cheng, T., and Zhang, X. (2012). Inhibition of titanium particle-induced inflammation by the proteasome inhibitor bortezomib in murine macrophage-like RAW 264.7 cells. Inflammation 35, 1411–1418. doi:10.1007/s10753-012-9454-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Marczyk, M., Macioszek, A., Tobiasz, J., Polanska, J., and Zyla, J. (2021). Importance of SNP dependency correction and association integration for gene set analysis in genome-wide association studies. Front. Genet. 12, 767358. doi:10.3389/fgene.2021.767358

PubMed Abstract | CrossRef Full Text | Google Scholar

Maseda, D., Bonami, R. H., and Crofford, L. J. (2014). Regulation of B lymphocytes and plasma cells by innate immune mechanisms and stromal cells in rheumatoid arthritis. Expert Rev. Clin. Immunol. 10, 747–762. doi:10.1586/1744666x.2014.907744

PubMed Abstract | CrossRef Full Text | Google Scholar

McLoed, A. G., Sherrill, T. P., Cheng, D. S., Han, W., Saxon, J. A., Gleaves, L. A., et al. (2016). Neutrophil-derived IL-1β impairs the efficacy of NF-κB inhibitors against lung cancer. Cell Rep. 16, 120–132. doi:10.1016/j.celrep.2016.05.085

PubMed Abstract | CrossRef Full Text | Google Scholar

Narayanan, R., and Honavar, S. G. (2017). A tale of two drugs: off and on-label. Indian J. Ophthalmol. 65, 549–550. doi:10.4103/ijo.ijo_586_17

PubMed Abstract | CrossRef Full Text | Google Scholar

Nascimento, A. S., Peres, L. L., Fari, A. V. S., Milani, R., Silva, R. A., da Costa Fernandes, C. J., et al. (2017). Phosphoproteome profiling reveals critical role of JAK-STAT signaling in maintaining chemoresistance in breast cancer. Oncotarget 8 (8), 114756–114768. doi:10.18632/oncotarget.21801

PubMed Abstract | CrossRef Full Text | Google Scholar

Neel, D. S., and Bivona, T. G. (2017). Resistance is futile: overcoming resistance to targeted therapies in lung adenocarcinoma. NPJ Precis. Oncol. 1, 3. doi:10.1038/s41698-017-0007-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Nzegwu, M., Uzoigwe, J., Omotowo, B., Ugochukwu, A., Ozumba, B., Sule, E., et al. (2021). Predictive and prognostic relevance of immunohistochemical testing of estrogen and progesterone receptors in breast cancer in South East Nigeria: a review of 417 cases. Rare Tumors 13, 20363613211006338. doi:10.1177/20363613211006338

PubMed Abstract | CrossRef Full Text | Google Scholar

Parikh, J. R., Klinger, B., Xia, Y., Marto, J. A., and Bluthgen, N. (2010). Discovering causal signaling pathways through gene-expression patterns. Nucleic Acids Res. 38, W109–W117. doi:10.1093/nar/gkq424

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, B., Ling, X., Huang, T., and Wan, J. (2023). HSP70 via HIF-1 alpha SUMOylation inhibits ferroptosis inducing lung cancer recurrence after insufficient radiofrequency ablation. PLoS One 18, e0294263. doi:10.1371/journal.pone.0294263

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahem, S. M., Epsi, N. J., Coffman, F. D., and Mitrofanova, A. (2020). Genome-wide analysis of therapeutic response uncovers molecular pathways governing tamoxifen resistance in ER+ breast cancer. EBioMedicine. 61, 103047. doi:10.1016/j.ebiom.2020.103047

PubMed Abstract | CrossRef Full Text | Google Scholar

Rashidiani, S., Mamo, G., Farkas, B., Szabadi, A., Farkas, B., Uszkai, V., et al. (2024). Integrative epigenetic and molecular analysis reveals a novel promoter for a new isoform of the transcription factor TEAD4. Int. J. Mol. Sci. 13, 25. doi:10.3390/ijms25042223

PubMed Abstract | CrossRef Full Text | Google Scholar

Rho, S. B., Byun, H. J., Kim, B. R., and Lee, C. H. (2022). Snail promotes cancer cell proliferation via its interaction with the BIRC3. Biomol. Ther. Seoul. 30, 380–388. doi:10.4062/biomolther.2022.063

PubMed Abstract | CrossRef Full Text | Google Scholar

Rydenfelt, M., Klinger, B., Klunemann, M., and Bluthgen, N. (2020). SPEED2: inferring upstream pathway activity from differential gene expression. Nucleic Acids Res. 2 (48), W307–W312. doi:10.1093/nar/gkaa236

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., Pawar, M., Van Dort, M. E., Galban, S., Welton, A. R., Thurber, G. M., et al. (2018). Ocular toxicity profile of ST-162 and ST-168 as novel bifunctional MEK/PI3K inhibitors. J. Ocul. Pharmacol. Ther. 34, 477–485. doi:10.1089/jop.2017.0126

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, G., Shang, C., Sun, L., Li, Y., Zhu, Y., Xiu, Z., et al. (2022). Ad-VT enhances the sensitivity of chemotherapy-resistant lung adenocarcinoma cells to gemcitabine and paclitaxel in vitro and in vivo. Invest. New Drugs 40, 274–289. doi:10.1007/s10637-021-01204-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanam, A., Gibson-Corley, K. N., Love-Homan, L., Ihejirika, N., and Simons, A. L. (2016). Interleukin-1 blockade overcomes erlotinib resistance in head and neck squamous cell carcinoma. Oncotarget 7 (7), 76087–76100. doi:10.18632/oncotarget.12590

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, M., Jiang, D., Hinsch, C., and Fan, W. (2010). Fulvestrant (ICI 182,780) sensitizes breast cancer cells expressing estrogen receptor alpha to vinblastine and vinorelbine. Breast Cancer Res. Treat. 121, 335–345. doi:10.1007/s10549-009-0472-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Talib, W. H., Alsayed, A. R., Barakat, M., Abu-Taha, M. I., and Mahmod, A. I. (2021). Targeting drug chemo-resistance in cancer using natural products. Biomedicines 29, 1353. doi:10.3390/biomedicines9101353

PubMed Abstract | CrossRef Full Text | Google Scholar

Tallen, G., Riabowol, K., and Wolff, J. E. (2003). Expression of p33ING1 mRNA and chemosensitivity in brain tumor cells. Anticancer Res. 23, 1631–1635.

PubMed Abstract | Google Scholar

Venet, D., Dumont, J. E., and Detours, V. (2011). Most random gene expression signatures are significantly associated with breast cancer outcome. PLoS Comput. Biol. 7, e1002240. doi:10.1371/journal.pcbi.1002240

PubMed Abstract | CrossRef Full Text | Google Scholar

Wakeling, A. E., Dukes, M., and Bowler, J. (1991). A potent specific pure antiestrogen with clinical potential. Cancer Res. 1 (51), 3867–3873.

PubMed Abstract | Google Scholar

Wang, K., Xing, Z. H., Jiang, Q. W., Yang, Y., Huang, J. R., Yuan, M. L., et al. (2019). Targeting uPAR by CRISPR/Cas9 system attenuates cancer malignancy and multidrug resistance. Front. Oncol. 9, 80. doi:10.3389/fonc.2019.00080

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Karikomi, M., MacLean, A. L., and Nie, Q. (2019). Cell lineage and communication network inference via optimization for single-cell transcriptomics. Nucleic Acids Res. 20 (47), e66. doi:10.1093/nar/gkz204

PubMed Abstract | CrossRef Full Text | Google Scholar

Wojtaszek, J. L., Chatterjee, N., Najeeb, J., Ramos, A., Lee, M., Bian, K., et al. (2019). A small molecule targeting mutagenic translesion synthesis improves chemotherapy. Cell 178, 152–159.e11. doi:10.1016/j.cell.2019.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, K., Liu, P., Yan, P., Liu, Y., Song, L., Liu, Y., et al. (2022). N6-methyladenosine reader YTH N6-methyladenosine RNA binding protein 3 or insulin like growth factor 2 mRNA binding protein 2 knockdown protects human bronchial epithelial cells from hypoxia/reoxygenation injury by inactivating p38 MAPK, AKT, ERK1/2, and NF-κB pathways. Bioengineered 13, 11973–11986. doi:10.1080/21655979.2021.1999550

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, R., and Dong, J. (2021). The Hippo signaling pathway in drug resistance in cancer. Cancers (Basel) 16, 13. doi:10.3390/cancers13020318

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Wang, H., Bathke, A. C., Harrar, S. W., Piepho, H. P., and Deng, Y. (2011). Gene set analysis for longitudinal gene expression data. BMC Bioinforma. 12, 273. doi:10.1186/1471-2105-12-273

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Xiao, B., Liu, Y., Wu, S., Xiang, Q., Xiao, Y., et al. (2024). Roles of PPAR activation in cancer therapeutic resistance: implications for combination therapy and drug development. Eur. J. Pharmacol. 964, 176304. doi:10.1016/j.ejphar.2023.176304

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Wang, Y., Sun, X., Zhang, X., Simone, N., and He, J. (2024). ELK1/MTOR/S6K1 pathway contributes to acquired resistance to gefitinib in non-small cell lung cancer. Int. J. Mol. Sci. 17, 25. doi:10.3390/ijms25042382

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, R. R., Zhao, L. P., Liu, L. S., Deng, F. A., Chen, X. Y., Jiang, X. Y., et al. (2021). Self-delivery nanomedicine to overcome drug resistance for synergistic chemotherapy. Biomater. Sci. 9, 3445–3452. doi:10.1039/d1bm00119a

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: chemotherapeutic resistance, pathway-responsive gene sets, drug discovery, drug combination therapy, precision oncology

Citation: Feng D, Hao J, Li L, Chen J, Liu X, Zhang R, Han H, Li T, Wang X, Li X, Yu L, Li B, Li J and Wang L (2025) Drug discovery for chemotherapeutic resistance based on pathway-responsive gene sets and its application in breast cancer. Front. Bioinform. 5:1661601. doi: 10.3389/fbinf.2025.1661601

Received: 08 July 2025; Accepted: 06 August 2025;
Published: 16 September 2025.

Edited by:

Vikram Dalal, Washington University in St. Louis, United States

Reviewed by:

Sanjaykumar Patel, University of Bridgeport, United States
Kuldeep Giri, Princeton University, United States

Copyright © 2025 Feng, Hao, Li, Chen, Liu, Zhang, Han, Li, Wang, Li, Yu, Li, Li and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Bing Li, YmluZ2xpanBuMjAwM0BhbGl5dW4uY29t; Jin Li, bGlqaW5AbXVobi5lZHUuY24=; Limei Wang, d2FuZ2xtQG11aG4uZWR1LmNu

These authors have contributed equally to this work and share first authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.