Genome-scale knockout simulation and clustering analysis of drug-resistant breast cancer cells reveal drug sensitization targets

Abstract

Anticancer chemotherapy is an essential part of cancer treatment, but the emergence of resistance remains a major hurdle. Metabolic reprogramming is a notable phenotype associated with the acquisition of drug resistance. Here, we develop a computational framework that predicts metabolic gene targets capable of reverting the metabolic state of drug-resistant cells to that of drug-sensitive parental cells, thereby sensitizing the resistant cells. The computational framework performs single-gene knockout simulation of genome-scale metabolic models that predicts genome-wide metabolic flux distribution in drug-resistant cells, and clusters the resulting knockout flux data using uniform manifold approximation and projection, followed by k-means clustering. From the clustering analysis, knockout genes that lead to the flux data near that of drug-sensitive cells are considered drug sensitization targets. This computational approach is demonstrated using doxorubicin- and paclitaxel-resistant MCF7 breast cancer cells. Drug sensitization targets are further refined based on proteome and metabolome data, which generate GOT1 for doxorubicin-resistant MCF7, GPI for paclitaxel-resistant MCF7, and SLC1A5 as a common target. These targets are experimentally validated where treating drug-resistant cancer cells with small-molecule inhibitors results in increased sensitivity of drug-resistant cells to doxorubicin or paclitaxel. The applicability of the developed framework is further demonstrated using drug-resistant triple-negative breast cancer cells. Taken together, the computational framework predicts drug sensitization targets in an intuitive and cost-efficient manner and can be applied to overcome drug-resistant cells associated with various cancers and other metabolic diseases.

Keywords: anticancer drug resistance; drug sensitization; genome-scale metabolic model; metabolic reprogramming; single-gene knockout simulation.

Fig. 1.

The overall scheme of predicting drug-sensitizing gene targets for MCF7/DXR. (A) Workflow of reconstructing and simulating cell-specific GEMs to predict gene targets that can sensitize the MCF7/DXR. Proteome data previously obtained for MCF7 and MCF7/DXR (28) were used to build cell-specific GEMs. The MCF7/DXR-specific GEMs were subjected to single-gene knockout simulation. The resulting knockout flux data were subjected to UMAP and k-means clustering. Knockout genes that clustered with that of MCF7 (red in the UMAP plot) were considered drug sensitization targets. (B) PCA (Left) and GO enrichment analysis (Right) of proteome data from MCF7 and MCF7/DXR. Enriched GO terms on biological processes of significant DEPs are shown on the Right. (C) PCA (Left) and KEGG enrichment analysis (Right) of metabolome data from MCF7 and MCF7/DXR. KEGG pathways significantly enriched with differentially abundant metabolites (DAMs) are shown on the Right. (D) Metabolic pathways significantly enriched with different reaction fluxes between MCF7 and MCF7/DXR, identified through flux enrichment analysis using their specific GEMs. The numbers indicate the ratio of reactions with different fluxes to the total number of reactions in each pathway. (E) UMAP plot of the knockout flux data as well as those from the MCF7- and MCF7/DXR-specific GEMs. Knockout genes corresponding to the knockout flux data clustered with the MCF7 flux data were considered potential drug sensitization targets (dotted red ellipse).

Fig. 2.

Selection of drug-sensitizing gene targets for MCF7/DXR. (A) Results of the joint pathway analysis of proteome and metabolome data from MCF7 and MCF7/DXR. Joint pathway analysis identified significantly affected metabolic pathways (in the blue box, surrounded by the dotted red lines; Left) in MCF7/DXR compared to MCF7. The table on the Right lists the pathways included in the blue box, along with their pathway impact (as defined in MetaboAnalyst 6.0) and -log₁₀(FDR adjusted P-value) values. (B) KEGG pathways of high impact scores based on the joint pathway analysis with the number of predicted drug sensitization targets indicated. Among the 92 initially predicted targets, 28 genes appeared to be associated with the presented pathways. (C) Metabolic map showing the selected drug-sensitizing gene targets (in red): IDH1, IDH2, GOT1, and SLC1A5. Red and blue bars across the graphs correspond to the averaged concentration of indicated metabolites in MCF7 and MCF7/DXR, respectively. The error bars indicate SD values. (D) Venn diagram showing the number of potential drug-sensitizing gene targets predicted using MCF7/DXR-specific GEMs generated from proteome and RNA-seq data. (E) Essentiality of the gene targets in the selected pathways (C) according to the 'CRISPRGeneEffect.csv' dataset available in the DepMap Public 24Q4 Files. The numbers indicate cell lines exhibiting a gene effect score of −0.5 or lower among 14 noncancerous and 1,164 cancerous cell lines when each gene was knocked out. Gene targets selected for experimental validation are shown in red.

Fig. 3.

Sensitization of MCF7/DXR to doxorubicin by inhibiting the predicted targets. (A) Dose–response curve of MCF7/DXR treated with doxorubicin and 25 µM of V-9302, an inhibitor of SLC1A5. Pink and gray shades are 95% CI across the doxorubicin concentrations for the mean cell viability from quadruplicate experiments. (B) The effects of treating MCF7/DXR with doxorubicin and V-9302. The synergy score map using the ZIP model for doxorubicin and V-9302 (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The x- and y-axes present the tested concentrations of the two drugs according to the 6×6 matrix configuration. The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (C) Dose–response curve of MCF7/DXR treated with doxorubicin and 100 µM of compound 2c, an inhibitor of GOT1. Pink and gray shades are 95% CI across the doxorubicin concentrations for the mean cell viability from triplicate experiments. (D) The effects of treating MCF7/DXR with doxorubicin and compound 2c. The synergy score map using the ZIP model for doxorubicin and compound 2c (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (E and F) (E) Dose–response curve of MCF7/DXR treated with doxorubicin when IDH1 or IDH2 was suppressed using corresponding siRNAs. Pink and gray shades are 95% CI across the doxorubicin concentrations for the mean cell viability from triplicate experiments. (F) The downregulation efficiency of IDH1 (Left) or IDH2 (Right) using corresponding siRNAs. Values are mean ± s.e.m. from triplicate experiments. Statistical significance was calculated using one-tailed Welch’s t test where ***P < 0.001 and ****P < 0.0001. (G and H) (G) The number of significantly changed metabolites in MCF7 and MCF7/DXR treated with V-9302 in comparison with MCF7/DXR without the treatment. (H) The direction of the concentration changes for commonly affected 21 metabolites from (G).

Fig. 4.

Prediction and selection of drug-sensitizing gene targets for MCF7/PR. (A) PCA (Left) and GO enrichment analysis (Right) of proteome data from MCF7 and MCF7/PR. Enriched GO terms on biological processes of significant DEPs are shown on the Right. (B) PCA (Left) and KEGG enrichment analysis (Right) of metabolome data from MCF7 and MCF7/PR. KEGG pathways significantly enriched with differentially abundant metabolites (DAMs) are shown on the Right. (C) Metabolic pathways significantly enriched with different reaction fluxes between MCF7 and MCF7/DXR, identified through flux enrichment analysis using their specific GEMs. The numbers indicate the ratio of reactions with different fluxes to the total number of reactions in each pathway. (D) UMAP plot of the knockout flux data as well as those from the MCF7-and MCF7/PR-specific GEMs. Knockout genes corresponding to the knockout flux data clustered with the MCF7 flux data were considered potential drug sensitization targets (dotted red ellipse). (E) Results of the joint pathway analysis of proteome and metabolome data from MCF7 and MCF7/PR. The table below lists the pathways included in the blue box (Upper), along with their pathway impact (as defined in MetaboAnalyst 6.0) and -log₁₀(FDR adjusted P-value) values (Lower). (F) KEGG pathways with high pathway impact scores based on the joint pathway analysis with the number of predicted drug sensitization targets indicated. Among the 77 initially predicted targets, 16 genes appeared to be associated with the presented pathways. (G) Metabolic map showing the selected drug-sensitizing gene target, GPI (in red). Red and blue bars across the graphs correspond to the averaged concentration of a metabolite in MCF7 and MCF7/PR, respectively. The error bars indicate the SD values. (H) Venn diagram showing the number of potential drug-sensitizing gene targets predicted using MCF7/PR-specific GEMs generated from proteome and RNA-seq data. (I) Essentiality of the gene targets in the selected pathways (G) according to the 'CRISPRGeneEffect.csv' dataset available in the DepMap Public 24Q4 Files. The numbers indicate cell lines exhibiting a gene effect score of −0.5 or lower among 14 noncancerous and 1,164 cancerous cell lines when each gene was knocked out. Gene targets selected for experimental validation are shown in red.

Fig. 5.

Sensitization of MCF7/PR to paclitaxel by inhibiting predicted drug sensitization targets. (A) Dose–response curve of MCF7/PR treated with paclitaxel and 40 mM of 2-deoxy-D-glucose (2-DG), an inhibitor of GPI. Pink and gray shades are 95% CI across the paclitaxel concentrations for the mean cell viability from triplicate experiments. (B) The effects of treating MCF7/PR with paclitaxel and 2-DG. The synergy score map using the ZIP model for paclitaxel and 2-DG (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The x- and y-axes present the tested concentrations of the two drugs according to the 6×6 matrix configuration. The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (C) Dose–response curve of MCF7/PR treated with paclitaxel and 40 µM of V-9302. Pink and gray shades are 95% CI across the paclitaxel concentrations for the mean cell viability from quadruplicate experiments. (D) The effects of treating MCF7/PR with paclitaxel and V-9302. The synergy score map using the ZIP model for paclitaxel and V-9302 (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The presented values correspond to the mean from quadruplicate experiments across the drug concentrations.

Fig. 6.

Sensitization of drug-resistant MDA-MB-231 cells to doxorubicin and paclitaxel by inhibiting the predicted targets. (A) Dose–response curve of MDA-MB-231/DXR treated with doxorubicin and 20 µM of V-9302. Pink and gray shades are 95% CI across the doxorubicin concentrations for the mean cell viability from quadruplicate experiments. (B) The effects of treating MDA-MB-231/DXR with doxorubicin and V-9302. The synergy score map using the ZIP model for doxorubicin and V-9302 (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The x- and y-axes present the tested concentrations of the two drugs according to the 6×6 matrix configuration. The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (C) Dose–response curve of MDA-MB-231/DXR treated with doxorubicin and 64 µM of compound 2c. Pink and gray shades are 95% CI across the doxorubicin concentrations for the mean cell viability from quadruplicate experiments. (D) The effects of treating MDA-MB-231/DXR with doxorubicin and compound 2c. The synergy score map using the ZIP model for doxorubicin and compound 2c (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (E) Dose–response curve of MDA-MB-231/PR treated with paclitaxel and 75 mM of 2-DG. Pink and gray shades are 95% CI across the paclitaxel concentrations for the mean cell viability from quadruplicate experiments. (F) The effects of treating MDA-MB-231/PR with paclitaxel and 2-DG. The synergy score map using the ZIP model for paclitaxel and 2-DG (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The presented values correspond to the mean from quadruplicate experiments across the drug concentrations. (G) Dose–response curve of MDA-MB-231/PR treated with paclitaxel and 15 µM of V-9302. Pink and gray shades are 95% CI across the paclitaxel concentrations for the mean cell viability from quadruplicate experiments. (H) The effects of treating MDA-MB-231/PR with paclitaxel and V-9302. The synergy score map using the ZIP model for paclitaxel and V-9302 (Left). The dose–response matrix showing the percentage of inhibition presented in the downward direction (Right). The presented values correspond to the mean from quadruplicate experiments across the drug concentrations.