
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)02084-4
10.1016/j.isci.2024.110859
110859
Article
Computational pipeline predicting cell death suppressors as targets for cancer therapy
Vinik Yaron 1
Maimon Avi 1
Raj Harsha 1
Dubey Vinay 1
Geist Felix 2
Wienke Dirk 2
Lev Sima sima.lev@weizmann.ac.il
13∗
1 Molecular Cell Biology Department, Weizmann Institute of Science, Rehovot 76100, Israel
2 The Healthcare Business of Merck KGaA, Darmstadt, Germany
∗ Corresponding author sima.lev@weizmann.ac.il
3 Lead contact

30 8 2024
20 9 2024
30 8 2024
27 9 11085915 3 2024
24 6 2024
28 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Summary

Identification of promising targets for cancer therapy is a global effort in precision medicine. Here, we describe a computational pipeline integrating transcriptomic and vulnerability responses to cell-death inducing drugs, to predict cell-death suppressors as candidate targets for cancer therapy. The prediction is based on two modules; the transcriptomic similarity module to identify genes whose targeting results in similar transcriptomic responses of the death-inducing drugs, and the correlation module to identify candidate genes whose expression correlates to the vulnerability of cancer cells to the same death-inducers. The combined predictors of these two modules were integrated into a single metric. As a proof-of-concept, we selected ferroptosis inducers as death-inducing drugs in triple negative breast cancer. The pipeline reliably predicted candidate genes as ferroptosis suppressors, as validated by computational methods and cellular assays. The described pipeline might be used to identify repressors of various cell-death pathways as potential therapeutic targets for different cancer types.

Graphical abstract

Highlights

• TNBC are highly susceptible to ferroptosis, especially the mesenchymal subtype

• A computational pipeline for the detection of targets for cancer therapy is introduced

• The pipeline is based on transcriptomic and vulnerability responses to drugs

• The pipeline can be generalized to other cell-death pathways and to various cancers

Cell biology; Bioinformatics; Cancer; Transcriptomics

Subject areas

Cell biology
Bioinformatics
Cancer
Transcriptomics
Published: August 30, 2024
==== Body
pmcIntroduction

Triple negative breast cancer (TNBC) is an aggressive disease characterized by high metastasis and poor prognosis.1 The disease is defined by the lack of estrogen and progesterone receptors and HER2 amplification, and currently, it has no effective, clinically approved therapeutic targets.2 Although the efficacies of various targeted therapies have been tested in TNBC pre-clinical models and many of them initially eliminated tumor cells through apoptosis, frequently, the tumor cells develop drug resistance and relapse. In recent years, targeting of non-apoptotic cell-death pathways, such as ferroptosis, has been considered as promising strategy for cancer therapy.3 We previously showed that TNBCs are particularly vulnerable to ferroptosis,4 implying that this death module can be exploited for TNBC therapy. However, canonical ferroptosis inducers, such as erastin, which inhibits the System Xc- of cystine/glutamate antiporter, and RSL3, which inhibits glutathione peroxidase 4 (GPX4), the major cellular protector of lipid peroxidation, have little efficacy in vivo and, therefore, limited clinical implications.5 These limitations boosted robust efforts to identify ferroptosis repressors, several have been identified through genetic screens, such as GCH1,6 FSP1,7 DHODH8 and MBOAT1/2.9 Additional genes are continuously discovered and might be potent, targetable suppressors with high anti-cancer efficacy.

Several computational methods have been recently developed to identify therapeutic targets utilizing multiomics profiling and machine learning.10 These methods can predict the outcome of gene targeting,11,12 employing several publicly available transcriptomic resources such as the connectivity map (CMAP), which includes thousands of perturbational datasets. The CMAP can predict the phenotypic outcome of gene targeting based on similarities between the transcriptomic responses induced by numerous pharmacological or genetic perturbations. This resource has been used for gene targets discovery,13 for drug repurposing,14 and for predicting highly potent drug combinations.15 Other approaches integrate gene expression profiles and susceptibility to various drugs to predict drug targets.16 This approach is provided, for example, by the Cancer Dependency Map (DepMap; Broad Institute) portal, which includes large datasets of gene expression, gene essentiality, and drug sensitivity across multiple cancer cell lines.17 Combining the transcriptomic and the vulnerability data of cancer cell lines with the data derived from patients with cancer has been previously used to build prediction models connecting gene expression to drug sensitivity.18 Since both the transcriptomic and vulnerability data rely on large-scale screens, which might yield false positives, their integration may increase the robustness and confidence of the predicted targets.

Here we employ a computation approach, integrating the transcriptomic and the vulnerability responses of TNBC cells to ferroptosis inducers (FINs), to eventually predict ferroptosis repressors as promising candidate genes for cancer therapy. Practically, we measured the vulnerability of multiple TNBC cell lines to FINs and concomitantly profiled the FINs transcriptomic response in the same TNBC cell lines. This experimental data was combined with public resources (CMAP, DepMAP) to assess the potential of a given gene to suppress ferroptosis based on two modules: a transcriptomic similarity module and a correlation module. The output of the pipeline is a ranked list of genes, which can predict ferroptotic cell death in response to their targeting, at least in TNBC. The prediction is based on the proximity of the candidate gene, on a dimension reduction plane, to other key players of ferroptosis, such as GPX4 and GCH1. The power and reliability of the pipeline were validated experimentally by selecting representative genes from the ranked list and demonstrating their ferroptosis suppressing capacity in TNBC. We propose that the integrated drug discovery pipeline can be used as a resource for targeting different cell death pathways in various cancer lineages.

Results

Pipeline design to identify ferroptosis suppressors in triple negative breast cancer

Considering the intrinsic vulnerability of TNBC to ferroptosis,4,19 we constructed a computational pipeline that identified potential suppressors of ferroptosis. We assume that the inhibition of these candidate genes should induce ferroptosis, at least in TNBC. We also assume that some of these targets might be used for cancer therapy, depending on critical considerations (toxicity, druggability, and so forth). The pipeline relies on a set of predictors calculated for each gene (Figure 1). These predictors quantify (i) the similarity between the transcriptomic response induced by targeting the candidate genes to the response induced by the FINs, and (ii) the correlation between the transcription levels of the candidate genes to the vulnerability - of the TNBC cell lines to the FINs. The set of predictors was then integrated into a single score that estimates the potential of each gene targeting to induce ferroptosis. Further computational and experimental validations highlighted the predictive power of the pipeline. The pipeline flow chart is depicted in Figure 1.Figure 1 Pipeline workflow

The pipeline flow chart. The pipeline was generated to predict candidate genes as ferroptosis suppressors using transcriptomic similarity and correlations modules.

Correlations between transcriptomic response and ferroptosis susceptibility

According to the pipeline design, we first measured the correlations between gene expression or gene essentiality to ferroptosis vulnerability, which was quantified by the area under the dose-response curve (AUCs) in TNBC cell lines. High correlations between the level of gene expression and the susceptibility to FINs as measured by the AUCs suggest that a decrease in the expression of these genes might sensitize cells to FINs, thus highlighting possible ferroptosis suppressors. For gene essentiality, we used the gene dependency scores (provided by the Achilles dataset, Broad Institute), which are given on a negative scale, with lower values indicating higher essentiality. Once again, a high correlation of gene dependency score to FINs AUCs might highlight potential repressors.

We previously showed that TNBCs are more susceptible to FINs (erastin and FIN56) compared to cell lines derived from other breast cancer subtypes.4 To expand these findings, we also assessed the susceptibility of TNBC and non-TNBC cell lines to RSL3, thereby analyzing the effects of the three canonical FINs of the major FIN classes (class 1- erastin, class 2-RSL3, and class 3-FIN56)20 using dose-response curves (Figure S1A). We quantified the area under the curves (AUCs), as a measure of ferroptosis susceptibility (Figure 2A), and observed, as expected, that TNBC are more susceptible to ferroptosis (Figure 2B).Figure 2 Assigning the correlation predictors for each gene

(A and B) Area under the curves (AUCs) were quantified from 16 dose-response curves obtained in the indicated breast cancer cell lines in response to erastin, RSL3, or FIN56 treatment. High AUCs represent a low vulnerability to the FINs. T-test (in B) was used to evaluate p-values (∗∗p-value<0.01, ∗∗∗p-value<0.001). The AUCs are averages of 3–4 repeats for each cell line and inducers.

(C–F) Pearson’s correlations were calculated between the expression level of genes (CCLE data, panels C and D) or the essentiality of genes (Achilles data, panels E and F) and the FIN AUCs. The Venn diagrams (C and E) show the number of genes with correlation above 0.65 or below −0.65. The genes with correlation above 0.65 for all 3 FINs are shown in (D) (gene expression vs. AUCs) or (F) (gene dependency vs. AUCs). In both panels, the top plot shows the correlations of the expression (D) or the dependency (F) to the 3 FINs. The bottom plot shows the results of a PubMed search for each gene with 6 ferroptosis related terms. A red box indicates citation(s) of the gene in the related term, while the color hue corresponds to the citation number.

Next, we correlated the AUCs of the 16 BC cell lines to gene expression (CCLE data) and gene dependency scores (Achilles data, Broad Institute). We found 41 genes with a high correlation (r > 0.65) between gene expression levels and the AUCs of the 3 FINs (Figures 2C and 2D). This analysis highlighted DCXR (dicarbonyl and L-xylulose reductase) as the highest correlative gene to erastin AUC, and ASPSCR1 (Tether containing UBX domain for GLUT4) as the highest correlative gene to all three FINs. Neither of them have been reported, thus far, as ferroptosis regulators (Figure 2D). Gene ontology analysis of these 41 genes revealed a significant enrichment of nucleosome-related genes (Figure S1B), which might be relevant for the ferroptosis related effect of HDAC inhibitors.21,22

Similarly, we found 54 genes with a high correlation of gene dependency to the AUCs (Figures 2E and 2F). In total, for each gene we acquired 6 correlation scores (correlation of the gene expression or the gene dependency to the AUC of the three FINs: erastin, RSL3 and FIN56) (Table S1). These scores will be henceforth termed the “correlation predictors”.

Transcriptomic response to ferroptosis inducers in triple negative breast cancer

The assigned 6 predictors described above rely on the correlations to FIN vulnerabilities, suggesting that genes with high predictor values would be associated with increased ferroptosis vulnerability. However, since correlations might not imply causation, we assumed that introducing an additional set of predictors could increase the confidence in the target selection process. We, therefore, introduced the transcriptomic similarity predictors, which reflect the similarity between the transcriptomic response to FINs treatment and the transcriptomic response of gene targeting.

To that end, we profiled the transcriptomic response of FINs, by treating 5 TNBC cell lines of different molecular subtypes, including two basal-like (BL) (MDA-MB-468, HCC70), and three mesenchymal-like (M) (Hs578, SUM-159, BT549), with either erastin or RSL3. We then performed RNAseq to determine the transcriptomic effect of those FINs in the TNBC cell lines. Principal component analysis (PCA) showed that the diversity in the transcriptomic response is mainly associated with differences between the cell lines, rather than the FIN identity (Figure 3A). This was further corroborated by the relatively small number of significant differentially expressed genes (DEGs) between erastin- or RSL3-treated cells vs. DMSO-treated cells (Figure 3B). Among the 10 RNAseq datasets (5 cell lines X 2 FINs), 126 genes were significantly up-regulated in 3 or more datasets (Figures S2A and S2B; Table S2), including 14 genes, which we recently identified as ferroptosis vs. apoptosis biomarkers.23 Analysis of the transcriptomic responses to FINs did not reveal specific expression patterns associated with either the subtype of the TNBC cell lines (BL vs. M) or the identity of the FINs (erastin vs. RSL3) (Figure S2A). However, gene set enrichment analysis (GSEA) showed enrichment of genes related to amino acid deprivation,24 and response to oxidized phospholipids, among other pathways that could be relevant to ferroptosis (Figure S2C). All genes were than ranked by their average fold change in expression in response to FINs versus DMSO. Figures 3C and S2D depict the 75 genes with the highest (Figure 3C) and lowest (Figure S2D) average fold change, representing the most up- and down-regulated genes in response to the applied FINs. Among the 75 up-regulated genes, 20 genes have already been identified as ferroptosis associated genes (FerrDB database,25), including genes with anti-ferroptotic effects, such as PCK226 or SLC1A5,27 among others.Figure 3 Transcriptomic analysis of FINs response in TNBC

(A–C) RNAseq was performed for 5 TNBC cell lines, treated either with erastin (marked with “E”), RSL3 (“R”), or DMSO as control (“C”). Experiments were done in duplicates.

(A) PCA analysis of the RNAseq data.

(B) The number of significant (p-value <0.05) differentially expressed genes was determined for each treatment vs. control (“_E” denotes erastin vs. DMSO, “_R” denotes RSL3 vs. DMSO).

(C) Heatmap showing the 75 genes with the highest average fold change expression between the FINs and DMSO. The points above the heatmap indicate genes that appear in the FerrDB and are known as ferroptosis markers (green), drivers (red), or suppressors (blue). The genes in red are ferroptosis-to-apoptosis biomarkers, as previously described.23

(D) The enrichment of the top 10–75 genes from the heatmap in (C) was assessed in 19 public FIN datasets (Table S3). The enrichment was performed using ssGSEA.

To further characterize this FINs-induced gene set in TNBC and to assess whether it represents a general transcriptomic landscape of ferroptotic response, we collected transcriptomic profiles of 19 publicly available ferroptosis responses to different chemical and genetic perturbations (Table S3) from GEO. Scoring the enrichment of the top-ranking genes (10–75) in the heatmap (Figure 3C) by ssGSEA revealed high enrichment in most of the 19 GEO profiles (Figure 3D), highlighting the importance of these genes in the general transcriptomic response to FINs.

Connectivity map analysis and extraction of transcription similarity scores

The second applied approach for ferroptosis target discovery was associated with transcriptomic similarity. In this approach, we assessed the similarities between the transcriptomic response of FINs obtained in TNBC cell lines (Figure 3C), to that of other chemical or genetic perturbations obtained from the Connectivity Map (CMAP). The CMAP database includes ∼475,000 transcriptomic profiles in 4–9 cancer cell lines in response to ∼19,800 small molecules and ∼7,500 genetic perturbations.28 Usually, the CMAP signatures are compared to a gene set of interest. In our case, we delved into the CMAP signatures database using a gene set composed of the 20 genes with the highest average fold change in their expression in response to the FINs in the 10 datasets as shown in Figure 3C. Indeed, we observed high scores for signatures of canonical FINs such as erastin and sorafenib of shRNA of GPX4. (Figure S2E).

Nevertheless, previous reports highlighted the variable transcriptomic response of a given perturbation across different cell lines.29 Therefore, instead of comparing the combined gene set of the 20 genes derived from all the 10 datasets described in Figure 3C, we used a cell-specific approach and examined the similarities between the CMAP signatures (29,515 signatures = 3932 knockdown genes X 4–9 cell lines) in each individual transcriptomic response of the 10 FINs datasets (Figure 4A).Figure 4 Assigning the transcriptomic similarity scores for each gene

(A) The shRNA consensus transcriptomic profiles were downloaded from the CMAP database. Up- and down-regulated gene sets were extracted, and their enrichment in the 10 RNAseq datasets (5 TNBC cells x 2 FINs) was measured using the CAMERA method.

(B) Plot showing the ranked transcriptomic similarity scores calculated for one representative of the 10 datasets (HCC70 treated with RSL3). In this plot, ∼3900 points are shown, each representing an shRNA from the CMAP database. The x- and y axis are the ranked similarity scores calculated using the upregulated and downregulated gene sets. A potential ferroptosis suppressor is predicted to induce a transcriptomic profile in which the top and the bottom 20 genes are ranked high and low, respectively, in the 10 RNAseq datasets, and therefore should be positioned in the bottom-right corner of the plot (green area).

(C) The final similarity score for each gene in each treatment was calculated as the upregulated rank minus the downregulated rank. The plot shows these scores for all genes in HCC70 treated with RSL3, as a representative.

(D) Predictor selection method. The matrix (left) shows the pairwise Pearson’s correlations between the 10 transcriptomic profiles of FINs in TNBC (shown in Figure 3C). Correlations were calculated using the fold change values of all genes in each profile. Based on these correlations, a metric was evaluated describing how each transcriptomic profile is similar to the consensus of all 10 profiles (right). To calculate this metric, the correlations in each row of the matrix were summed up, excluding the negative correlations and the correlations in the diagonal (which equal 1). These summations were then normalized on a scale of 0–100 and are shown in the plot to the right. The green bars represent the predictors with the highest similarity to the consensus of all profiles, while the gray bars are the predictors with the lowest such similarity, which were filtered out in the next steps.

(E) Validation of predictor selection by text-mining. Sets of genes were delved in PubMed (using the GeneShot tool), to measure the percentage of genes in each set which is cited with the terms shown in the y axis. The pink density plots represent 1023 combinations of predictors, in each the top 50 genes with the highest values of selected predictors were analyzed. In blue, 2000 sets of randomly selected genes were similarly analyzed. The number on the left indicates the mean ± standard deviation of the random sets. The red dot, and the numbers under “selected combination,” represent the actual predictor combination selected for this study (the ∗ indicates p-value <0.05, comparing the selected combination to the random distribution). The black dot represents the predictor combination which includes all the predictors generated.

We used the 29,515 shRNA consensus signatures from the CMAP and extracted from each signature the top 20 and bottom 20 genes, which represent the most upregulated and downregulated genes of each signature. This yielded 59,030 gene sets in total. Next, we used the CAMERA method30 to quantify the enrichment of those gene sets in the 10 transcriptomic datasets of FINs (Figure 4A; Table S4). These scores, henceforth termed the “Transcriptomic Similarity Scores,” represent the similarity between the transcriptomic profiles of ferroptosis induction and the transcriptomic profile of the 3932 knocked-down (KD) genes.

For each KD gene (shRNA), the CMAP database contains 4–9 transcriptomic profiles, one for each cell line in the CMAP project. The 4–9 transcriptomic similarity scores of each shRNA were aggregated into one by selecting the maximum value (Figure S3A). We then ranked the shRNA based on their similarity scores (both for the upregulated and the downregulated gene sets) (Figure S3B). The calculated ranked similarity scores of the 3932 shRNAs to the transcriptomic response of HCC70 cells treated with RSL3 are shown in Figure 4B as an example. The x- and y axis of the plot are the ranked similarity scores calculated using the upregulated (x axis) and downregulated (y axis) gene sets. Promising ferroptosis suppressor candidates should, therefore, be ranked high among the upregulated gene sets similarity scores (Figure S3B, left panel), and low among the downregulated gene sets similarity scores (Figure S3B, right panel). Indeed, GPX4 and GCH1, two canonical ferroptosis repressors, fulfill these criteria, as seen in Figure 4B (both appear at the lower right corner of the plot).

Next, we combine the upregulated and downregulated ranks into a single metric for each shRNA by subtracting the downregulated rank from the upregulated rank. The obtained final scores ranged from ∼4000 (best potential ferroptosis suppressors) to ∼ -4000 (best potential ferroptosis inducers) for each shRNA (Figure 4C). These 10 scores per shRNA (one for each transcriptomic dataset of FINs in TNBC) reflect the similarities between the transcriptomic response of each shRNA from the CMAP to the transcriptomic response of the FINs and were considered to be the 10 transcriptomic similarity predictors. For most predictors, the values are associated with known ferroptosis suppressors. For example, the 40 genes from the FerrDB database of ferroptosis suppressors, which were also analyzed in the CMAP database, were indeed highly enriched compared to other genes (Figure S3C).

To further optimize the pipeline, we introduced a predictor selection step. This step was applied due to the variability in the transcriptomic responses of FINs in TNBC (Figure S2A), which was associated with a high variance between the 10 transcriptomic similarity predictors for each shRNA. To reduce the predictors' variance, we calculated the consensus level between the 10 transcriptomic profiles induced by the FINs in TNBC (Figure 4D, see STAR Methods for details) and then removed the 3 predictors with the lowest similarity to the consensus of all the 10 profiles. This predictor selection step resulted in 7 optimized transcriptomic similarity predictors, which were further combined to the 6 correlation predictors calculated earlier (Figure 2).

To demonstrate the advantage of the predictor selection step, which reduced the number of the transcriptomic similarity predictors from 10 to 7, we used a text-mining analysis. In brief, the 6 correlation predictors were combined with any possible combination of the 10 original transcriptomic similarity predictors (1023 predictor combinations in total), and the top 50 genes for each combination with the highest predictor values were selected. We then calculated the percentage of genes, among these 50 genes, cited with ferroptosis related terms in PubMed (Figure 4E). To evaluate the statistical significance of this analysis, we similarly analyzed randomly selected sets of genes. This analysis indicated that combining the 6 correlation predictors with the 7 selected transcriptomic similarity predictors yields a list of genes with high enrichment for genes cited with those terms, in particular with the term “lipid peroxidation,” compared to most predictor combinations.

Integrating the transcriptomic similarity scores with the correlation scores

Thus far, we acquired 13 predictors for each gene, holding values that measure their potential of being ferroptosis suppressors: 6 correlation scores per gene from the correlation module, and 7 transcriptomic similarity scores from the transcriptomic similarity module. To visualize the 13 predictors, we projected them into a plane using the UMAP (Uniform Manifold Approximation and Projection) method (Figure 5A). Each point in this UMAP represents a single gene (3679 total genes, each with 13 predictors values), with adjacent points sharing similar predictor values. We, therefore, could reduce the 13 predictors into a single metric describing the distance between points in the UMAP projection. In parallel, the overall predictor values for each gene were aggregated (Figure S3D, see STAR Methods), and are marked in the UMAP by a blue hue, with a stronger color representing genes with higher predictor values. Among these, are the canonical ferroptosis regulators GCH1, a major ferroptosis suppressor involved in BH4 synthesis,31 GPX4,32 and CBS, a key regulator of the transsulfuration pathway.33Figure 5 Integrating the transcriptomic similarity and correlation scores to predict potential ferroptosis suppressors

(A) UMAP projection of 13 scores (7 transcriptomic similarity and 6 correlations score) per gene. A selection of known ferroptosis regulators are labeled. The blue scale of the dots (genes) represents an averaging of all the predictor values (see STAR Methods section).

(B) Zoom-in on the GCH1 node neighborhood in the UMAP.

(C) Genes from the GCH1 neighborhood (shown in B), organized by their Euclidean distance from the GCH1 node (top panel). The bottom plot is a PubMed search of each gene with 6 ferroptosis related terms. Red boxes indicate citation(s) with the specified term, while the color hue corresponds to the number of citations. Genes that are not essential pan-cancer, as well as those selected for experimental validation, are marked by green dots.

(D and E) Correlations of gene expression levels (CCLE pan-cancer cell line data) to the dependency score of GPX4 (taken from the Achilles dataset) (in E) or to the AUCs of 4 FINs taken from the CTRP dataset (in D, shown is the average correlation of each gene to the 4 FINs). The red dots indicate the 75 closest neighbors of the GCH1 node in the UMAP shown in Figure (B). Bottom: Gene set enrichment analysis for the top 75 GCH1 closest neighbors, was performed by ranking all the genes based on the correlation shown in the plots above. Normalized enrichment scores (NES) and their p-values are indicated.

(F) Genes selected for experimental validation from the GCH1 neighborhood, showing their Euclidean distance to the GCH1 node, and their values for the correlation (3 out of the 6) and transcriptomic similarity (all the 7) predictors. The boxplots to the right of each row of the table show where these genes (colored dots) are positioned against the distribution of all 3679 genes in the UMAP.

We then focused on GCH1, the 5th gene with the highest average predictor values. We assumed that genes in the GCH1 local neighborhood (Figure 5B; Table S5), would share similar predictor values as GCH1, and therefore could potentially be ferroptosis suppressors. Indeed, many genes cited with ferroptosis related terms in PubMed were found adjacent to the GCH1 node (Figure 5C). Gene ontology enrichment analysis revealed that several genes in the GCH1 neighborhood have oxidoreductase activity (Figure S4A), which can be linked to ferroptosis. To further validate the potential of our prediction, we examined the correlation between all genes in the CCLE dataset to the FIN AUCs (taken from the CTRP drug sensitivity dataset) in pan-cancer (Figure 5D) and to GPX4 dependency score (taken from Achilles dataset, Broad) (Figure 5E). A similar approach was previously used to predict FSP1/AIFM2 as a potential ferroptosis regulator.34 Ranking those correlations, we found a positive, significant enrichment of the 75 genes closest to the GCH1 node in the UMAP (Figures 5D and 5E). These results validate the predictive power of the pipeline outcome. Furthermore, we performed the same analysis for all the 1023 possible combinations of predictors (Figure S4B), and found that most of the combinations gave positive, significant enrichment. These results suggest that the GCH1 cluster could be predicted as a potential ferroptosis suppressors regardless of the predictor selection process. Nevertheless, the selected predictors combinations had one of the highest enrichment scores, suggesting that the selection process indeed optimized the results.

Experimental validation of predicted candidate genes as ferroptosis suppressors

Based on our computation analysis, we predicted that genes in the GCH1 neighborhood would share similar predictor values, and therefore, might suppress ferroptosis. To validate this prediction and highlight the power of the pipeline, we selected several genes from the GCH1 neighborhood for experimental validation. Several considerations were applied for target selection. First and most importantly, we ranked the genes by their Euclidean distance from the GCH1 node in the UMAP (Figure 5C). We filtered out genes which were previously cited in PubMed together with ferroptosis related terms (with few exceptions). To avoid general toxicity, we selected only genes that are not commonly essential in multiple cancer lineages, and therefore their targeting might have fewer side effects. Finally, we prioritized genes which are druggable or ligandable to enable subsequent identification of small molecule inhibitors against those targets.

Among the closest neighbors to the GCH1 node, we proceeded with 7 genes for validation: DECR1, DCXR, NDUFV1, SORD, CANT1, TNFRSF18, and FRAT1. These genes indeed have high values in most of the predictors, especially the correlations of gene expression to AUC (Figure 5F). While the GCH1 neighborhood appeared to be promising, considering the high predictor values of GCH1 (Figure S3D), the GPX4 neighborhood was also examined, given the prominent ferroptotic suppressing activity of GPX4 and its relatively high predictor values (ranked #138 among all genes). Applying similar considerations (such as druggability and toxicity), we selected two additional genes, CTBP1 and PDAP1 from the closest neighborhood of the GPX4 in the UMAP projection (Figures S4C and S4D). As mentioned, most of the 9 selected target genes were not previously associated with ferroptosis (as determined by PubMed search), and are all non-pan essential (Figures S5A and S5B), and druggable targets.

The major consideration in the target selection method is based mainly on the UMAP projection. As the UMAP projection involves a degree of randomness dictated by its initialization seed, we ensured that the genes selection is stable and reproducible regardless the randomness. To this end, we repeated the analysis for 500 UMAP projections, each built with a different random seed (Figure S4E), demonstrating that the selected genes indeed remain in the neighborhood of GCH1 or GPX4. In addition, we further validated that the proximity of the points within the UMAP projection indeed recapitulate proximity in the predictor values (Figure S4F).

For the experimental validation, the 9 selected genes were knocked down by shRNAs in MDA-MB-468 and HCC70 basal-like TNBC cell lines (Figure S5C), and their influence on ferroptotic death was examined by characteristic assays.35 As shown in Figure 6A, knocking down these genes significantly reduced cell viability in the two TNBC cell lines, with the exception of NDUFV1 in MDA-MB-468 cells. In most cases, the inhibitory effect on cell viability was rescued by the ferroptosis inhibitor, ferrostatin (Figure 6B). Next, we assessed the effect of the genes KD on lipid peroxidation using the fluorescent sensor BODIPY-C11. The increased levels of the oxidized form of the sensor upon genes knockdown are demonstrated by representative confocal images (Figure 6C), along with the quantification of the fluorescence intensity (Figure 6D). This increase in lipid peroxidation was accompanied with enhanced 4HNE (4-Hydroxynonenal) levels shown in the Western blots (WB) in Figure 6E. In those experiments, GPX4 targeting by shRNA was used as a positive control. Importantly, one of the selected genes in the GPX4 neighborhood, PDAP1 (PDGFA-associated protein 1), was identified independently in our recent study using a different approach.23 We showed that the depletion of PDAP1 not only induced ferroptosis in basal-like breast cancer cells in vitro but also induced ferroptosis in basal-like breast tumors in a xenograft mouse model and consequently inhibited tumor growth,23 suggesting that the identified ferroptosis suppressors could be used for cancer therapy. Four of the genes (CANT1, DECR1, PDAP1, NDUFV1) had significant prognostic effects, as determined by Kaplan-Meier analysis using relapse-free survival data (Figure 6F). Two of these genes, CANT1 and NDUFV1, were also prognostic using overall survival data as determined by Cox analysis (Figure 6G). The expression level of most of the genes is high in TNBC tumors (Figure 6H), further suggesting that their targeting could be beneficial for TNBC therapy.Figure 6 Experimental validation highlighting the ferroptosis inducing potential of gene targeting

(A and B) Influence of candidate genes on cell viability. The indicated genes were knocked down (KD) in the indicated TNBC cell lines (HCC70, MDA-MB-468). KD efficiency was assessed by qPCR (Figure S5C), while cell viability was evaluated by MTT assay. Cells were seeded in 96 well plates (50% confluency), without (A) or with Ferrostatin (5 μM) (B), and 72 h later, cell viability in the control (pLKO) or genes KD cells was measured. Percent of cell viability (A) and percent rescue with Ferrostatin (B) are shown (mean ± SD of at least two independent experiments). KD of GPX4 was used as a positive control. P-values were measured by one-sample t-test (∗p-value<0.05, ∗∗p-value<0.01, ∗∗∗p-value<0.001, ∗∗∗∗p-value<0.0001).

(C and D) Knock-down of the indicated genes was performed in MDA-MB-468 cells by shRNA. Bodipy-C11 staining was used to visualize lipid peroxidation. Representative confocal images are shown (C) (scale bar, 10μm). Fluorescence measurements are given in (D). pLKO (empty vector) and shGPX4 were used as a negative and positive control, respectively. In (D), means ± SD are shown of two repeats.

(E) Western blot analysis of 4HNE-adduct, for MDA-MB-468 in which the indicated genes were knocked down. Tubulin was used as a housekeeping gene.

(F) Kaplan-Meier plots for genes with significant differences in prognosis.

(G) Hazard ratios (HR) were calculated for the indicated genes in patients with breast cancer (MetaBric data, n = 1904 patients) using the Cox proportional hazards model.

(H) Expression of the indicated genes in patients, taken from TCGA data (n = 114 normal, 879 non-TNBC, and 180 TNBC).

Overall, we validated 9 potential ferroptosis repressors in TNBC, many of which were shown to have an inhibitory effect in BC upon targeting (Table 1). While these genes may enhance ferroptosis through different mechanisms, we found some common gene ontologies, including oxidoreductase activity, which might be associated with many ferroptosis regulators (Figure S5D).Table 1 Key effects of validated genes in breast cancer (BC) or other cancer types

Gene	Effect in cancer	Reference	
CANT1	Prognostic factor, promotes cell proliferation and invasion in lung adenocarcinoma	Yao et al.36	
CTBP1	Elevated in BC, targeting increases sensitivity to chemotherapy	Birts et al.37, Deng et al.38	
DCXR	Promote BC cell proliferation. Targeting inhibits tumor growth in vivo	Jin et al.39	
DECR1	Targeting induces ferroptosis in prostate cancer in vitro and in vivo	Nassar et al.40	
FRAT1	Highly expressed in TNBC	Nam et al.41	
NDUFV1	Increases BC proliferation, however targeting increases the metastatic ability	Kim and Singh42, Santidrian et al.43	
PDAP1	Targeting induces ferroptosis in BC; targeting inhibits tumor growth in vivo	Vinik et al.23	

Discussion

In this study, we established a powerful computational framework applying a unique approach of integrating the correlation and the transcriptomic similarity modules to identify suppressors of cell-death pathways as candidate targets for cancer therapy. Considering the vulnerability of TNBC to ferroptosis and the lack of effective therapeutic targets for patients with TNBC, we focused on ferroptosis as the targeting death pathway.

The integrated pipeline relies on (i) the correlations between gene expression/essentiality and ferroptosis sensitivity, which yielded 6 predictors per gene, and (ii) the similarity between the transcriptomic response of gene targeting and the transcriptomic response to FINs, which yielded 10 predictors (Figure 1). To reduce the variability associated with the diverse transcriptomic response to different FINs, we reduced the number of the transcriptomic similarity predictors, by evaluating the similarities between the 10 profiles and filtering out the transcriptomic profiles most dissimilar to the others. This analysis enabled the reduction of the transcriptomic similarity predictors from 10 to 7, which were eventually combined with the 6 correlation predictors using dimensionality reduction by UMAP. We assumed that the proximity in the UMAP of a gene to canonical ferroptosis repressors (mainly GCH1 and GPX4) is the best metric in the selection process of ferroptosis suppressors. Indeed, experimental validation of 9 genes (7 genes in GCH1 proximity, and 2 from GPX4 neighborhood) demonstrated the predicting power of the pipeline, as their knockdown enhanced lipid peroxidation and reduced TNBC cell viability (Figure 6). Hence, we established a powerful resource for the identification of candidate ferroptosis suppressors.

A major challenge in pipeline optimization was associated with predictor selection, which was a necessary process due to the high variability in the transcriptomic responses to various FINs in different biological models and time points.44 To overcome the high transcriptomic variability, we selected the transcriptomic similarity predictors unbiasedly by assessing their transcriptomic similarity to each other (Figure 4D). We validated the selection process by text-mining analysis of all the possible combinations of predictors (1023 in total) (Figure 4E). This analysis highlights the statistical significance of the predictor selection process and suggests that the text-mining can be used by itself as a method for predictor selection.

Importantly, we also demonstrated that the GCH1 neighborhood (Figure 5B) is statistically significantly enriched with potential targets for ferroptosis (Figures 5D and 5E). A similar analysis for all the possible predictor combinations (Figure S4B) suggests that while the predictor selection step can indeed optimize the results, the pipeline has the potential to predict ferroptosis suppressors regardless of this step. This is an important observation, which highlights the potential of generalizing the pipeline to other death pathways. Furthermore, we showed that the distances between points on the UMAP projection significantly correlate with the distances in the predictor matrix (Figure S4F), confirming the hypothesis that genes with similar predictor values are clustered together. Finally, we made sure that the output of our pipeline is stable for randomness (Figure S4E) considering that the UMAP is a stochastic algorithm. Taken together, we established a statistically significant algorithm that could accurately predict ferroptosis repressors as potential candidates for cancer therapy.

The identified candidate genes may impose their ferroptosis suppression through major known ferroptosis axes,45 including the GSH-GPX4 axis, the FSP1-CoQ10 axis,7,46 the GCH1-BH4 axis,6 the DHODH-CoQH2 axis,8 or by as of yet unknown mechanism. Among the 9 genes that were experimentally validated, DECR1 was previously shown to protect prostate tumor cells from ferroptosis, mainly through the regulation of PUFA oxidation.40 DECR1 was positioned very closely to the GCH1 node, and exhibited high predictor values, similar to GCH1. DCXR was also placed in close proximity to GCH1, had high predictor values, and was found to be the most correlating gene to erastin AUC (Figure 2D). However, its involvement in ferroptosis has not been reported so far. DCXR plays a role in xylitol metabolism,47 which in turn can affect glutathione levels,48 but its role in ferroptosis might not be directly related to the GCH1-BH4 axis and requires further investigation. CTBP1, one of the GPX4's closest neighbors, was also validated to be a good target for ferroptosis induction. This gene, a metabolic sensor of redox status, was recently found to be involved in ferroptosis.49 Another neighbor of GPX4 is PDAP1, which we recently identified as a ferroptosis repressor and demonstrated its potential as a target for TNBC therapy in xenograft TNBC mouse model.23 Importantly, we showed that PDAP1 depletion induced ferroptotic death both in vitro and in vivo, possibly through the AKT-mTOR-SREBP1-SCD1 axis, further demonstrating the power of the described pipeline and the potential of the discovered candidate targets for cancer therapy.

It is important to note that SLC7A11, the actual target of erastin, was ranked very low in this pipeline, suggesting that the output of our approach can produce valid ferroptosis inducing targets, but not all possible targets will be ranked high in the pipeline. In the case of SLC7A11, this can be explained by a negative correlation of its expression to FINs AUCs, which is counter intuitive. While SLC7A11 correlated with FIN AUCs pan-cancer (as shown in Figure 5D), its correlation specifically in breast cancer is poor (also when comparing to CTRP data).

Notably, the CMAP approach was previously used to identify drugs that increased the expression of ferroptosis related genes,50,51 and uncovered a connection between ferroptosis and HDAC inhibition.52 Correlations of gene essentiality to ferroptosis susceptibility identified FSP1 as a major ferroptosis regulator.34 The integration of both methods was proven to be highly beneficial. The integration mode increased the confidence in target selection and reduced the rate of false positives, which are common in many correlation-based drug discovery methods. The described pipeline, combining the correlations to FINs sensitivity together with transcriptomic similarities, followed by the dimension reduction and ranking potential targets by their proximity to major nodes, can potentially be streamlined to produce a computational pipeline revealing new targets to any desired phenotype induced by specific drugs, and not only to ferroptosis. Such a pipeline will require vulnerability data of cancer cells to the drug, and their transcriptomic response, and will produce a ranked list of genes, from which the targets can be selected based on other pharmacological considerations, such as novelty, general toxicity, and druggability. Hence, it can provide a powerful computational tool that can be used to discover targets of many different pathways affecting cancer growth, progression, and patient survival.

Limitations of the study

While the established pipeline has the potential to successfully predict suppressors of different cell-death pathways beyond ferroptosis, we validated its potency only for ferroptosis in the scope of this study. Further target validations of other death modules, in future projects, would greatly increase the applications of this pipeline. For other pathways, we propose to improve the predictor selection process for the pathways exhibiting highly variable predictors. In such pathways, a predictor selection process should be performed based on the predictors’ variance as described in this study and validated with text-mining or with another suitable method. Importantly, for each selected death pathway, the pipeline requires pre-existing information on the pathway in question to assign pathway nodes, such as GCH1 and GPX4 for ferroptosis. Eventually, it would be beneficial to predict death-pathway suppressor candidates based solely on the predictor values, rather than the proximity to major nodes.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Sima Lev (Sima.Lev@weizmann.ac.il).

Material availability

This study did not generate new unique reagents.

Data and code availability

• RNAseq data have been deposited at GEO and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. This article also analyzes existing, publicly available data. These accession numbers for the datasets are listed in the key resources table.

• The code written for this article is available in Github (https://github.com/SimaLevLab/Pipeline-for-discovery-of-ferroptosis-targets), and was also deposited at Zenodo. The code is publicly available as of the date of publication. DOI is listed in the key resources table.

• Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request.

Acknowledgments

Sima Lev is the incumbent of the Joyce and Ben B. Eisenberg Chair of Molecular Biology and Cancer Research. We thank Itay Tirosh for reviewing the article, and Sven Lindermann, Klaus Urbahns and other 10.13039/100004334 Merck members for their valuable input throughout the project.

This study is supported by 10.13039/100009945 Merck KGaA , Darmstadt, Germany.

Author contributions

YV designed and performed the computational analysis and wrote the article. AM designed and performed the experiments in Figure 6. HR and VD prepared the RNA for RNAseq and performed the experiment shown in Figure S1A. FG and DW contributed to the discussion and reviewed the article. SL supervised the study, contributed to the analysis and interpretation of the data, ensured financial support, and wrote the article.

Declaration of interests

Authors declare that they have no competing interests.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
	
anti-4-HNE	Sigma	AB5605; RRID: AB_569332	
Anti-beta-Tubulin	Cell signaling	2146; RRID: AB_2210545	
	
Chemicals, peptides, and recombinant proteins	
	
Erastin	Sigma	E7781	
RSL3	Sigma	SML2234	
FIN56	Sigma	SML1740	
Ferrostatin-1	Sigma	SML0583	
	
Critical commercial assays	
	
C11-BODIPY 581/591 fluorescence reporter	Cayman	27086	
	
Deposited data	
	
RNA-seq for MDA-MD-468 and HCC70	this paper	GEO: GSE235201	
RNAseq for SUM159, HS578 and MDA-MB-231	this paper	GEO: GSE255459	
CMAP signature dataset	Broad Institute	GEO: GSE92742	
CCLE dataset, version 20Q4	DepMap project	https://depmap.org/portal/	
Achilles dataset, version 21Q1	DepMap project	https://depmap.org/portal/	
FerrDB version 1		http://www.zhounan.org/ferrdb/legacy/index.html	
	
Experimental models: Cell lines	
	
MDA-MB-231	ATCC	RRID: CVCL_0062	
MDA-MB-436	ATCC	RRID: CVCL_0623	
SUM159	ATCC	RRID: CVCL_5423	
HS578T	ATCC	RRID: CVCL_0332	
BT549	ATCC	RRID: CVCL_1092	
MDA-MB-468	ATCC	RRID: CVCL_0419	
BT20	ATCC	RRID: CVCL_0178	
HCC1937	ATCC	RRID: CVCL_0290	
HCC70	ATCC	RRID: CVCL_1270	
SKBR3	ATCC	RRID: CVCL_A2GI	
JIMT1	ATCC	RRID: CVCL_2077	
T47D	ATCC	RRID: CVCL_0553	
HCC1143	ATCC	RRID: CVCL_1245	
MDA-MB-453	ATCC	RRID: CVCL_0418	
MCF-7	ATCC	RRID: CVCL_0031	
BT474	ATCC	RRID: CVCL_0179	
HEK293T	ATCC	RRID: CVCL_0063	
	
Oligonucleotides	
	
Primers for qRT-PCR experiment, see Table S6	this paper	Table S6	
	
Recombinant DNA	
	
pLKO-SORD	Sigma	TRCN0000028052	
pLKO-FRAT1	Sigma	TRCN0000062463	
pLKO-NDUFV1	Sigma	TRCN0000025872	
pLKO-CANT1	Sigma	TRCN0000051898	
pLKO-TNFRSF18	Sigma	TRCN0000058373	
pLKO-DECR1	Sigma	TRCN0000046514	
pLKO-DCXR	Sigma	TRCN0000038955	
pLKO-CTBP1	Sigma	TRCN0000285086	
pLKO-PDAP1	Sigma	TRCN0000299988	
pLKO-GPX4	Sigma	TRCN0000046252	
	
Software and algorithms	
	
R version 4.4.0	https://www.r-project.org/		
RStudio version 2024.04.2	Posit Software		
KMplotter	https://kmplot.com/analysis/		
ZEN imaging software, blue edition, version 3.2	Zeiss		
Code generated in this study	This paper	https://doi.org/10.5281/zenodo.13370885	

Experimental model and study participant details

Cell lines

The breast cancer cell lines MDA-MB-231, MDA-MB-436, SUM159, HS578, BT549, MDA-MB-468, BT20, HCC1937, HCC70, SKBR3, JIMT1, T47D, HCC1143, MDA-MB-453, MCF-7, BT474 (all of female origin) and the human embryonic kidney (HEK) 293T cells were originally obtained from the American Type Culture Collection (USA). The cells were grown in RPMI (all breast cancer lines) or DMEM (HEK293 cells) containing 10% fetal bovine serum (Gibco BRL, USA) and penicillin/streptomycin. Cells were cultured at 37°C in a humidified incubator of 5% CO2. Cell lines were routinely (once a month) checked for mycoplasma using a commercially available kit (Biological Industries, Israel).

Method details

Measurement of cell sensitivity to ferroptosis

Breast cancer cell lines were seeded in 96-wells plate (50% confluency), and 24 hr later were treated with erastin, RSL3 or FIN56 using different doses. Cell viability was measured 72 hr later by MTT assay (M2128, SIGMA). Cell viability (taken as % of control untreated (DMSO) cells) was plotted against drug concentration, and dose response curves were fitted to the data by the drm function from the drc package in R. AUCs (area under the curves) were extracted using the computeAUC function from the PharmacoGx package in R.

Generation of the correlation scores

Pearson’s correlations between the AUCs of the 3 FINs (acquired from the dose response curves) and gene expression levels (extracted from the CCLE data, 20Q4 version) and gene dependency scores (extracted from the Achilles dataset, 21Q1 version) was performed using the cor function in R. CCLE and Achilles datasets were downloaded from DepMap. Venn diagrams showing the genes with high and low correlations in all 3 FINs were done using the VennDiagram package.

RNA sequencing

Five TNBC cell lines (MDA-MB-468, HCC70, BT549, SUM159, Hs578) were treated with erastin (9 hours) and RSL3 (4 hours) using the IC50 concentrations which were deduced from the dose response curves created above. DMSO was used as control. Total RNA was extracted using TRI Reagent (Sigma-Aldrich), and its quality was assessed using Agilent 4200 TapeStation System (Agilent Technologies, Santa Clara, CA). RNA-seq libraries were generated by applying a bulk adaptation of the MARS-seq protocol, as previously described.53 Libraries were sequenced by the Illumina Novaseq 6000 using SP mode 100 cycles kit (Illumina). Mapping of sequences to the genome and generation of the count matrix was performed by the UTAP pipeline (Weizmann Institute). Libraries normalization, filtration of low count genes and discovery of differentially expressed genes was performed using the edgeR and Limma packages in R. PCA plot was generated using the factoextra package. UpSet plot was generated using the ComplexUpset package.

Validation of the ferroptosis response signature

19 public datasets, representing the transcriptomic response of ferroptosis induction in different cell lines and models, were downloaded from GEO (Table S3). In each dataset, the log2 fold changes in gene expression between the inducer and its control were calculated using Limma. The enrichment, in those datasets, of genes upregulated by FINs, was measured using the ssGSEA method, implemented by the gsva function from the GSVA package in R. The top 20 ranking genes were uploaded to the connectivity map site (https://clue.io/). The results were analyzed using the campR package in R.

Generation of transcriptomic similarity scores

The CMAP signature dataset was downloaded from GEO (accession number GEO: GSE92742). The level 5 signatures were downloaded, together with all relevant meta-data. The signatures were filtered to include only the consensus shRNA signatures (trt_sh.cgs). From each signature, the top 20 and bottom 20 genes were taken to form the upregulated and downregulated gene sets associated with that signature. The enrichment of all those gene sets in our 10 treatment points was performed by the Camera method from the Limma package in R.30 In CMAP, each shRNA was tested in several cell lines – the scores for each shRNA were aggregated into one score by taking the maximum value of the original scores. The resulting scores were then ranked. For each shRNA, the final transcriptomic similarity score was calculated to be the ranked score of the upregulated set minus the ranked score of the downregulated set. shRNAs with higher score are therefore considered to induce transcriptomic response similar to those of the FINs, i.e. their target genes have a potential to be ferroptosis suppressors.

As a filtration step for the transcriptomic similarity scores, we adapted a method from Smith et al.,54 by which we measured the similarities between the 10 transcriptomic profiles obtained from the RNAseq. In brief, we calculated the Pearson correlation for every pair of transcriptomic profiles, using the log2 fold changes for all genes in each profile. For each profile, the 9 correlation scores to the other 9 profiles were summed up, excluding any negative correlations or the correlation of each profile to itself. These summations were normalized on a scale from 0 to 100. This normalized metric represents the similarity of each transcriptomic profile to their average profile.

PubMed search for ferroptosis publications

For given sets of genes (the top correlating genes to the FIN AUCs, or the GCH1 closest neighbors), we performed a PubMed search with terms related to ferroptosis (“Ferroptosis”, “Gluthathione”, “GPX4”, “Iron”, “Lipid Peroxidation”, “SLC7A11”). This search was automated using the rentrez package in R, and is updated for November 2023. In our method for validation of the predictor selection method, we generated gene sets as follows: (i) Predictor combination sets: we combined the 6 correlation module predictors, together with any combination of the 10 transcriptomic similarity predictors (in total, ∑k=11010!k!(10−k)!=1023 combinations). For each such combination, we took the top 50 genes with the highest predictor values. (ii) Random sets: 2000 sets of 100 randomly selected genes each. Each of those 3023 sets were analyzed for citations with the same terms in PubMed. Due to the large number of sets, for this analysis we used the Geneshot tool.55

Combining the predictor scores

The transcriptomic similarity and correlation predictors were normalized and used to define a single metric that will predict genes with ferroptosis suppressor activity. The UMAP method was used to project the transcriptomic similarity and correlation scores of all genes on a 2D plane. The UMAP was implemented using the UMAP function in R, setting the “min_dist” parameter to 0.05 and keeping the rest of the parameters as default. The prediction metric was determined to be the Euclidean distance in the UMAP projection of each gene from a central node, selected for being a major ferroptosis suppressor (such as GCH1). In another approach, we averaged all the predictor values as follows: We calculated the average of the 3 correlation scores between gene expression levels and the 3 FINs AUCs, as well as the average of the 3 correlation scores between gene dependency scores and the 3 FINs AUCs. Then, the maximum average of the two was taken. This maximum value was ranked, and averaged with the ranking of the average of all transcriptomic similarity scores. This measure was highlighted on the UMAP in blue scale. In addition, we measured the distance of the genes from each other by creating a distances matrix based on the normalized predictor values, and compared those distances to the distances measured on the UMAP, to validate that local neighborhoods in the UMAP indeed correlates with actual similarities in predictor values.

To validate the hypothesis that the genes closest to the GCH1 node in the UMAP have ferroptosis suppression potential, the pan-cancer CCLE data (version 20Q4), Achilles data (version 21Q1) and CTRP drug response data were downloaded from Depmap portal (https://depmap.org/portal/). Gene expression levels were correlated to the GPX4 dependency score (taken from the Achilles dataset) and to FINs AUCs (Erastin, RSL3, ML210, ML162, taken from the CTRP dataset). To verify that the selected set of genes (GCH1 node neighbors) are among the highest correlating genes in both cases, we performed gene set enrichment analysis, using the fgsea function from the fgsea package in R. In the fgsea function, the enrichment score of the selected genes was measured against a ranked list of the genes, ranked by the above correlations of expression to GPX4 essentiality, and expression to FINs AUCs. The same analysis was also done on all 1023 possible combinations of predictors; for each such combination a UMAP was built and the top 75 genes with the closest Euclidean proximity to GCH1 were analyzed in the same way.

Gene knockdown

Lentiviral vectors encoding shRNAs of the validated genes were purchased from Sigma. The catalog number of the shRNAs used are: SORD (TRCN0000028052), FRAT1 (TRCN0000062463), NDUFV1 (TRCN0000025872), CANT1 (TRCN0000051898), TNFRSF18 (TRCN0000058373), DECR1 (TRCN0000046514), DCXR (TRCN0000038955), CTBP1 (TRCN0000285086), PDAP1 (TRCN0000299988) and GPX4 (TRCN0000046252). 3rd generation lenti-viruses were produced in HEK293T cells by co-transfecting the pLKO.1-shRNA constructs with pLP (pVSVG), pLP1 (pMDL) and pLP2 (pREV) expressing the virus envelope protein. Twelve hours later, the medium was changed and viruses were collected 24 and 48 h later. Viral supernatants were filtered through 0.45 μm pore size filters and stored at -80°C. Cells were infected with the viruses in the presence of 8 μg/ml polybrene for 24 h and then subjected to selection with media containing 1 μg/ml puromycin for 72-96 h. Cells infected with the empty virus were used as control. Following knockdown, cell viability assay was performed using MTT as described above, in the presence of the ferroptosis inhibitor Ferrostatin-1 (5μM) or DMSO as control.

Lipid peroxidation measurement

Lipid peroxidation in live cells was detected by the C11-BODIPY 581/591 fluorescence reporter (#27086, CAYMAN). Cells were plated in High-Content Imaging Glass Bottom 96-well Microplates (Cellvis) for ∼24 hr and then treated with the indicated drugs. The lipid peroxidation sensor C11-BODIPY (581/591) (7 μM) was added for 30-40 min together with 1 mM Hoechst 33342 (Sigma-Aldrich) in regular RPMI full media. Cells were gently washed twice with PBS and incubated in live cell imaging solution (Invitrogen, A14291DJ). Fluorescence was measured at 581/590 nm (excitation/emission) for the reduced dye, and at 488/510 nm (excitation/emission) for the oxidized dye, while Hoechst was measured at 350/461 nm (excitation/emission) using the Infinite 200 PRO Tecan microplate reader (Tecan Inc., Switzerland). The green-to-red fluorescence intensity ratio was used to measure lipid peroxidation. The values were normalized to cells number in each well using Hoechst staining values. Confocal microscopy images of live cells stained with C11-BODIPY were acquired using the LSM800 (Zeiss), 40X oil lens and the ZEN Imaging Software (Zeiss).

Lipid peroxidation was also estimated by the level of 4-HNE (4-hydroxynonenal) protein adducts as detected by anti-4-HNE antibody (AB5605, SIGMA) in western blot analysis. For this experiment, cells were lysed in lysis buffer containing 0.5% Triton X-100, 50 mM Hepes pH 7.5, 100 mM NaCl, 1 mM MgCl2, 50 mM NaF, 0.5 mM NaVO3, 20 mM β-glycerophosphate, 1 mM phenylmethylsulfonyl fluoride, 10 μg ml−1 leupeptin, and 10 μg ml−1 aprotinin. Cell lysates were centrifuged at 14,000 rpm for 15 min at 4°C. Protein concentration of the supernatants was measured by Bradford assay (Bio-Rad, Hercules, CA). Equal amounts of total protein (30–50 μg per sample) were analyzed by SDS–PAGE (polyacrylamide gel electrophoresis) and Western Blotting was performed. Anti-Tubulin antibody (Cell Signaling, 2146) was used as loading control.

RNA extraction and real time PCR

Total RNA from cell lines was extracted and purified using TRI Reagent (Sigma-Aldrich). RNA was reverse-transcribed into complementary DNA (cDNA) using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems; Cat. No. 4368814) with random primers according to the manufacturer’s instructions. Real-time PCR analysis was performed in QuantStudio-3 Real-Time PCR system (Applied Biosystems, Thermo Fisher Scientific) using SYBR Green Master Mix reagents (Roche) according to the manufacturer’s guidelines. House-keeping gene β-actin was used for normalization. The relative levels of mRNA were calculated using the ΔΔCT method. Primer sequences are listed in Table S6.

Patient dataset analysis

Expression levels of genes of breast cancer patients was taken either from the METABRIC or TCGA datasets. Univariate Cox regression was done using the coxph() function from the `survival` package, using the overall survival follow-up data included in the MetaBric dataset. Gene ontology analysis was performed by the gost() function in the gProfiler2 package in R. Kaplan-Meier plots were performed using KMplotter,56 using relapse-free survival as the event, and with automatic selection of expression cutoff (in total, 2032 patients examined).

Quantification and statistical analysis

The statistical test used to determine the significance levels for each experiment was done in R and is described in the figure legends. In most cases (unless otherwise mentioned in the figure legends), t-test was used to measure statistical significance. P-values lower than 0.05 were considered as statistically significant.

Supplemental information

Document S1. Figures S1–S5 and Tables S3 and S6

Table S1. Correlation of AUCs to gene expression and essentiality, related to Figure 2

Table S2. Differentially expressed genes based on FIN induced TNBC lines, related to Figure 2

Table S4. Connectivity map enrichment scores, related to Figure 4

Table S5. Gene list - GCH1 nearest neighbors, related to Figure 5

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110859.
==== Refs
References

1 Baranova A. Krasnoselskyi M. Starikov V. Kartashov S. Zhulkevych I. Vlasenko V. Oleshko K. Bilodid O. Sadchikova M. Vinnyk Y. Triple-negative breast cancer: current treatment strategies and factors of negative prognosis J. Med. Life 15 2022 153 161 10.25122/jml-2021-0108 35419095
2 Lev S. Targeted therapy and drug resistance in triple-negative breast cancer: the EGFR axis Biochem. Soc. Trans. 48 2020 657 665 10.1042/BST20191055 32311020
3 Zhang C. Liu X. Jin S. Chen Y. Guo R. Ferroptosis in cancer therapy: a novel approach to reversing drug resistance Mol. Cancer 21 2022 47 10.1186/s12943-022-01530-y 35151318
4 Verma N. Vinik Y. Saroha A. Nair N.U. Ruppin E. Mills G. Karn T. Dubey V. Khera L. Raj H. Synthetic lethal combination targeting BET uncovered intrinsic susceptibility of TNBC to ferroptosis Sci. Adv. 6 2020 eaba8968 10.1126/sciadv.aba8968
5 Hadian K. Stockwell B.R. A roadmap to creating ferroptosis-based medicines Nat. Chem. Biol. 17 2021 1113 1116 10.1038/s41589-021-00853-z 34675413
6 Kraft V.A.N. Bezjian C.T. Pfeiffer S. Ringelstetter L. Müller C. Zandkarimi F. Merl-Pham J. Bao X. Anastasov N. Kössl J. GTP Cyclohydrolase 1/Tetrahydrobiopterin Counteract Ferroptosis through Lipid Remodeling ACS Cent. Sci. 6 2020 41 53 10.1021/acscentsci.9b01063 31989025
7 Doll S. Freitas F.P. Shah R. Aldrovandi M. da Silva M.C. Ingold I. Goya Grocin A. Xavier da Silva T.N. Panzilius E. Scheel C.H. FSP1 is a glutathione-independent ferroptosis suppressor Nature 575 2019 693 698 10.1038/s41586-019-1707-0 31634899
8 Mao C. Liu X. Zhang Y. Lei G. Yan Y. Lee H. Koppula P. Wu S. Zhuang L. Fang B. DHODH-mediated ferroptosis defence is a targetable vulnerability in cancer Nature 593 2021 586 590 10.1038/s41586-021-03539-7 33981038
9 Liang D. Feng Y. Zandkarimi F. Wang H. Zhang Z. Kim J. Cai Y. Gu W. Stockwell B.R. Jiang X. Ferroptosis surveillance independent of GPX4 and differentially regulated by sex hormones Cell 186 2023 2748 2764.e22 10.1016/j.cell.2023.05.003 37267948
10 You Y. Lai X. Pan Y. Zheng H. Vera J. Liu S. Deng S. Zhang L. Artificial intelligence in cancer target identification and drug discovery Signal Transduct. Targeted Ther. 7 2022 156 10.1038/s41392-022-00994-0
11 Roohani Y. Huang K. Leskovec J. Predicting transcriptional outcomes of novel multigene perturbations with GEARS Nat. Biotechnol. 42 2024 927 935 10.1038/s41587-023-01905-6 37592036
12 Lotfollahi M. Klimovskaia Susmelj A. De Donno C. Hetzel L. Ji Y. Ibarra I.L. Srivatsan S.R. Naghipourfar M. Daza R.M. Martin B. Predicting cellular responses to complex perturbations in high-throughput screens Mol. Syst. Biol. 19 2023 e11517 10.15252/msb.202211517
13 Shah I. Bundy J. Chambers B. Everett L.J. Haggard D. Harrill J. Judson R.S. Nyffeler J. Patlewicz G. Navigating Transcriptomic Connectivity Mapping Workflows to Link Chemicals with Bioactivities Chem. Res. Toxicol. 35 2022 1929 1949 10.1021/acs.chemrestox.2c00245 36301716
14 Zhao Y. Chen X. Chen J. Qi X. Decoding Connectivity Map-based drug repurposing for oncotherapy Brief. Bioinform. 24 2023 bbad142 10.1093/bib/bbad142
15 Musa A. Ghoraie L.S. Zhang S.D. Glazko G. Yli-Harja O. Dehmer M. Haibe-Kains B. Emmert-Streib F. A review of connectivity map and computational approaches in pharmacogenomics Brief. Bioinform. 19 2018 506 523 10.1093/bib/bbw112 28069634
16 Qin Y. Conley A.P. Grimm E.A. Roszik J. A tool for discovering drug sensitivity and gene expression associations in cancer cells PLoS One 12 2017 e0176763 10.1371/journal.pone.0176763
17 Tsherniak A. Vazquez F. Montgomery P.G. Weir B.A. Kryukov G. Cowley G.S. Gill S. Harrington W.F. Pantel S. Krill-Burger J.M. Defining a Cancer Dependency Map Cell 170 2017 564 576.e16 10.1016/j.cell.2017.06.010 28753430
18 Li Y. Umbach D.M. Krahn J.M. Shats I. Li X. Li L. Predicting tumor response to drugs based on gene-expression biomarkers of sensitivity learned from cancer cell lines BMC Genom. 22 2021 272 10.1186/s12864-021-07581-7
19 Doll S. Proneth B. Tyurina Y.Y. Panzilius E. Kobayashi S. Ingold I. Irmler M. Beckers J. Aichler M. Walch A. ACSL4 dictates ferroptosis sensitivity by shaping cellular lipid composition Nat. Chem. Biol. 13 2017 91 98 10.1038/nchembio.2239 27842070
20 Nie Q. Hu Y. Yu X. Li X. Fang X. Induction and application of ferroptosis in cancer therapy Cancer Cell Int. 22 2022 12 10.1186/s12935-021-02366-0 34996454
21 Zille M. Kumar A. Kundu N. Bourassa M.W. Wong V.S.C. Willis D. Karuppagounder S.S. Ratan R.R. Ferroptosis in Neurons and Cancer Cells Is Similar But Differentially Regulated by Histone Deacetylase Inhibitors eNeuro 6 2019 1 10.1523/ENEURO.0263-18.2019
22 Oliveira T. Hermann E. Lin D. Chowanadisai W. Hull E. Montgomery M. HDAC inhibition induces EMT and alterations in cellular iron homeostasis to augment ferroptosis sensitivity in SW13 cells Redox Biol. 47 2021 102149 10.1016/j.redox.2021.102149
23 Vinik Y. Maimon A. Dubey V. Raj H. Abramovitch I. Malitsky S. Itkin M. Ma'ayan A. Westermann F. Gottlieb E. Programming a Ferroptosis-to-Apoptosis Transition Landscape Revealed Ferroptosis Biomarkers and Repressors for Cancer Therapy Adv. Sci. 11 2024 e2307263 10.1002/advs.202307263
24 Yang J. Dai X. Xu H. Tang Q. Bi F. Regulation of Ferroptosis by Amino Acid Metabolism in Cancer Int. J. Biol. Sci. 18 2022 1695 1705 10.7150/ijbs.64982 35280684
25 Zhou N. Yuan X. Du Q. Zhang Z. Shi X. Bao J. Ning Y. Peng L. FerrDb V2: update of the manually curated database of ferroptosis regulators and ferroptosis-disease associations Nucleic Acids Res. 51 2023 D571 D582 10.1093/nar/gkac935 36305834
26 Bluemel G. Planque M. Madreiter-Sokolowski C.T. Haitzmann T. Hrzenjak A. Graier W.F. Fendt S.M. Olschewski H. Leithner K. PCK2 opposes mitochondrial respiration and maintains the redox balance in starved lung cancer cells Free Radic. Biol. Med. 176 2021 34 45 10.1016/j.freeradbiomed.2021.09.007 34520823
27 Chen P. Jiang Y. Liang J. Cai J. Zhuo Y. Fan H. Yuan R. Cheng S. Zhang Y. SLC1A5 is a novel biomarker associated with ferroptosis and the tumor microenvironment: a pancancer analysis Aging (Albany NY) 15 2023 7451 7475 10.18632/aging.204911 37566748
28 Lamb J. Crawford E.D. Peck D. Modell J.W. Blat I.C. Wrobel M.J. Lerner J. Brunet J.P. Subramanian A. Ross K.N. The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease Science 313 2006 1929 1935 10.1126/science.1132939 17008526
29 Baillif B. Wichard J. Méndez-Lucio O. Rouquié D. Exploring the Use of Compound-Induced Transcriptomic Data Generated From Cell Lines to Predict Compound Activity Toward Molecular Targets Front. Chem. 8 2020 296 10.3389/fchem.2020.00296 32391323
30 Wu D. Smyth G.K. Camera: a competitive gene set test accounting for inter-gene correlation Nucleic Acids Res. 40 2012 e133 10.1093/nar/gks461 22638577
31 Soula M. Weber R.A. Zilka O. Alwaseem H. La K. Yen F. Molina H. Garcia-Bermudez J. Pratt D.A. Birsoy K. Metabolic determinants of cancer cell sensitivity to canonical ferroptosis inducers Nat. Chem. Biol. 16 2020 1351 1360 10.1038/s41589-020-0613-y 32778843
32 Seibt T.M. Proneth B. Conrad M. Role of GPX4 in ferroptosis and its pharmacological implication Free Radic. Biol. Med. 133 2019 144 152 10.1016/j.freeradbiomed.2018.09.014 30219704
33 Wang L. Cai H. Hu Y. Liu F. Huang S. Zhou Y. Yu J. Xu J. Wu F. A pharmacological probe identifies cystathionine beta-synthase as a new negative regulator for ferroptosis Cell Death Dis. 9 2018 1005 10.1038/s41419-018-1063-2 30258181
34 Bersuker K. Hendricks J.M. Li Z. Magtanong L. Ford B. Tang P.H. Roberts M.A. Tong B. Maimone T.J. Zoncu R. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis Nature 575 2019 688 692 10.1038/s41586-019-1705-2 31634900
35 Lee J.Y. Kim W.K. Bae K.H. Lee S.C. Lee E.W. Lipid Metabolism and Ferroptosis Biology 10 2021 184 10.3390/biology10030184
36 Yao Q. Yu Y. Wang Z. Zhang M. Ma J. Wu Y. Zheng Q. Li J. CANT1 serves as a potential prognostic factor for lung adenocarcinoma and promotes cell proliferation and invasion in vitro BMC Cancer 22 2022 117 10.1186/s12885-022-09175-2 35090419
37 Birts C.N. Harding R. Soosaipillai G. Halder T. Azim-Araghi A. Darley M. Cutress R.I. Bateman A.C. Blaydes J.P. Expression of CtBP family protein isoforms in breast cancer and their role in chemoresistance Biol. Cell 103 2010 1 19 10.1042/BC20100067 20964627
38 Deng Y. Guo W. Xu N. Li F. Li J. CtBP1 transactivates RAD51 and confers cisplatin resistance to breast cancer cells Mol. Carcinog. 59 2020 512 519 10.1002/mc.23175 32124501
39 Jin Y. Zhang M. Tong Y. Qiu L. Ye Y. Zhao B. DCXR promotes cell proliferation by promoting the activity of aerobic glycolysis in breast cancer Mol. Med. Rep. 27 2023 31 10.3892/mmr.2022.12918
40 Nassar Z.D. Mah C.Y. Dehairs J. Burvenich I.J. Irani S. Centenera M.M. Helm M. Shrestha R.K. Moldovan M. Don A.S. Human DECR1 is an androgen-repressed survival factor that regulates PUFA oxidation to protect prostate tumor cells from ferroptosis Elife 9 2020 e54166 10.7554/eLife.54166
41 Nam S.E. Ko Y.S. Park K.S. Jin T. Yoo Y.B. Yang J.H. Kim W.Y. Han H.S. Lim S.D. Lee S.E. Kim W.S. Overexpression of FRAT1 protein is closely related to triple-negative breast cancer Ann. Surg. Treat. Res. 103 2022 63 71 10.4174/astr.2022.103.2.63 36017142
42 Kim S.H. Singh S.V. The FoxQ1 transcription factor is a novel regulator of electron transport chain complex I subunits in human breast cancer cells Mol. Carcinog. 61 2022 372 381 10.1002/mc.23381 34939230
43 Santidrian A.F. Matsuno-Yagi A. Ritland M. Seo B.B. LeBoeuf S.E. Gay L.J. Yagi T. Felding-Habermann B. Mitochondrial complex I activity and NAD+/NADH balance regulate breast cancer progression J. Clin. Invest. 123 2013 1068 1081 10.1172/JCI64264 23426180
44 Vinik Y. Lev S. A snapshot into the transcriptomic landscape of apoptosis and ferroptosis in cancer Cell Death Dis. 15 2024 371 10.1038/s41419-024-06766-8 38811541
45 Wang D. Tang L. Zhang Y. Ge G. Jiang X. Mo Y. Wu P. Deng X. Li L. Zuo S. Regulatory pathways and drugs associated with ferroptosis in tumors Cell Death Dis. 13 2022 544 10.1038/s41419-022-04927-1 35688814
46 Nakamura T. Hipp C. Santos Dias Mourão A. Borggräfe J. Aldrovandi M. Henkelmann B. Wanninger J. Mishima E. Lytton E. Emler D. Phase separation of FSP1 promotes ferroptosis Nature 619 2023 371 377 10.1038/s41586-023-06255-6 37380771
47 Nakagawa J. Ishikura S. Asami J. Isaji T. Usami N. Hara A. Sakurai T. Tsuritani K. Oda K. Takahashi M. Molecular characterization of mammalian dicarbonyl/L-xylulose reductase and its localization in kidney J. Biol. Chem. 277 2002 17883 17891 10.1074/jbc.M110703200 11882650
48 Tomonobu N. Komalasari N.L.G.Y. Sumardika I.W. Jiang F. Chen Y. Yamamoto K.I. Kinoshita R. Murata H. Inoue Y. Sakaguchi M. Xylitol acts as an anticancer monosaccharide to induce selective cancer death via regulation of the glutathione level Chem. Biol. Interact. 324 2020 109085 10.1016/j.cbi.2020.109085
49 Li Y. Hu G. Huang F. Chen M. Chen Y. Xu Y. Tong G. MAT1A Suppression by the CTBP1/HDAC1/HDAC2 Transcriptional Complex Induces Immune Escape and Reduces Ferroptosis in Hepatocellular Carcinoma Lab. Invest. 103 2023 100180 10.1016/j.labinv.2023.100180
50 Chen Z. Wu T. Yan Z. Zhang M. Identification and Validation of an 11-Ferroptosis Related Gene Signature and Its Correlation With Immune Checkpoint Molecules in Glioma Front. Cell Dev. Biol. 9 2021 652599 10.3389/fcell.2021.652599
51 Liu B. Fu X. Du Y. Feng Z. Liu X. Li Z. Yu F. Zhou G. Ba Y. In Silico Analysis of Ferroptosis-Related Genes and Its Implication in Drug Prediction against Fluorosis Int. J. Mol. Sci. 24 2023 4221 10.3390/ijms24044221
52 Zhang T. Sun B. Zhong C. Xu K. Wang Z. Hofman P. Nagano T. Legras A. Breadner D. Ricciuti B. Targeting histone deacetylase enhances the therapeutic effect of Erastin-induced ferroptosis in EGFR-activating mutant lung adenocarcinoma Transl. Lung Cancer Res. 10 2021 1857 1872 10.21037/tlcr-21-303 34012798
53 Jaitin D.A. Kenigsberg E. Keren-Shaul H. Elefant N. Paul F. Zaretsky I. Mildner A. Cohen N. Jung S. Tanay A. Amit I. Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types Science 343 2014 776 779 10.1126/science.1247651 24531970
54 Smith I. Greenside P.G. Natoli T. Lahr D.L. Wadden D. Tirosh I. Narayan R. Root D.E. Golub T.R. Subramanian A. Doench J.G. Evaluation of RNAi and CRISPR technologies by large-scale gene expression profiling in the Connectivity Map PLoS Biol. 15 2017 e2003213 10.1371/journal.pbio.2003213
55 Lachmann A. Schilder B.M. Wojciechowicz M.L. Torre D. Kuleshov M.V. Keenan A.B. Ma'ayan A. Geneshot: search engine for ranking genes from arbitrary text queries Nucleic Acids Res. 47 2019 W571 W577 10.1093/nar/gkz393 31114885
56 Gyõrffy B. Survival analysis across the entire transcriptome identifies biomarkers with the highest prognostic power in breast cancer Comput. Struct. Biotechnol. J. 19 2021 4101 4109 10.1016/j.csbj.2021.07.014 34527184
