To identify pathway activities and potential genomic variants associated with TKI-induced cardiotoxicity (TIC), we treated six hiPSC-derived cardiomyocyte cell lines from six healthy individuals (Supplementary Fig. 1) (Supplementary Data 1) with one of 27 TKIs, 4 anthracyclines and 26 other cardiac and non-cardiac-acting drugs for 48 hours, using therapeutically relevant drug concentrations (Supplementary Data 2A, Supplementary Table 116). Untreated control data was generated at the same time and has been previously published as part of characterization of these cell lines16. In this study, we used the untreated control data from each experiment to identify the drug-treatment-induced differentially expressed genes (DEGs). All drugs (except endothelin and TNF-alpha) are FDA-approved. The group of TKIs contained 23 small molecule TKIs and 4 monoclonal antibodies against TKs. They could be separated into 10 cardiotoxic and 17 non-cardiotoxic TKIs based on the results of clinical studies (Supplementary Data 3) and FAERS data analysis (Supplementary Fig. 2). Bulk transcriptomic analysis of control16 and drug-treated cell lines generated 266 lists of DEGs (Supplementary Data 4) (Supplementary Data 5 shows averaged DEGs). All lists of DEGs, each representing the response observed for one sample, i.e., a unique cell line/drug combination (Supplementary Fig. 3A), were subjected to pairwise correlation and hierarchical clustering (Supplementary Fig. 3B).Singular value decomposition to reveal drug-selective expression responsesTo identify drug-selective expression responses we used SVD to search for shared gene expression components in multiple cell lines treated with the same drug (Fig. 1A, Supplementary Fig. 4). Clustering of the unprocessed transcriptomic data only grouped a few samples treated with the same drugs into the same cluster (Fig. 1B, Supplementary Fig. 3B). In contrast, the identity of the cell line, or the amplitude of the drug response, i.e., the number of significant DEGs, mainly determined the clustering results (Supplementary Fig. 5A). To determine if we could identify drug-selective clustering, we calculated one F1 score for each drug that documents how close all cell lines treated with that drug cluster together (Supplementary Fig. 5B). With possible values larger 0 up to 1, the median F1-score of 0.116 documents very low drug-selective clustering efficiencies (Fig. 1C).Fig. 1: Singular value decomposition identifies drug-selective gene expression responses.266 samples, each representing a unique combination of one out of three to six hiPSC-derived cardiomyocyte cell lines treated with one out of 54 drugs, were subjected to bulk RNAseq analysis. 266 lists of differentially expressed genes (DEGs) were calculated using the negative-binomial test implemented in the ‘exactTest’ functionality of the edgeR package. A Our computational pipeline uses singular value decomposition (SVD) to identify drug-selective gene expression responses that are components of the complete responses. Flow chart is used with permission from Mount Sinai Health System, licensed under CC BY. See methods section and Supplementary Fig. 4 for details. B Pairwise correlation analysis followed by hierarchical clustering reveals that most drug responses are dominated by cell-line-selective effects hiding drug-selective effects. Heatmap colors describe drugs used for treatment, as documented in E. See Supplementary Fig. 5A for larger dendrogram. C We used F1 score statistics to document the clustering efficiency, i.e., how close samples treated with the same drug cluster together. Low clustering efficiencies quantitatively describe the finding that only a few complete DEG responses are dominated by drug-selective effects. D Projection of the complete DEG responses into each of the identified 54 drug-selective subspaces greatly increases the clustering efficiencies for all 54 drugs. Numbers of treated cell lines are shown below the bars. Orange: Small molecule kinase inhibitors (KI), red: monoclonal antibodies against KIs, purple: anthracyclines, blue: cardiac-acting drugs, turquoise: non-cardiac-acting drugs. E Pairwise correlation of 266 merged drug-selective responses, followed by hierarchical clustering, documents that SVD allows identification of components induced in all cell lines treated with the same drug. Clusters that contain drugs with similar mechanisms are labeled with gray bars. White insets indicate drugs in those clusters that are not part of the outlined mechanisms. White circles indicate outlier samples that were identified by our pipeline and cluster as outliers in the merged dataset as well. See Supplementary Fig. 13 for a larger dendrogram. #: count of, Ø: without, &: and.SVD (Supplementary Fig. 6A) of all 266 samples identified 266 eigenarrays (Supplementary Fig. 6B) whose linear combination using sample-specific coefficients gives the complete gene expression profiles. The first eigenarray correlates (Supplementary Fig. 6C) with the amplitude of the response (Supplementary Fig. 6D), and its most prominent genes (Supplementary Data 6) enrich for muscle contraction (Supplementary Fig. 6E, Supplementary Data 7). These results indicate that the first eigenarray describes a general response to perturbation that masks drug-selective effects. Removal of the first eigenarray (Supplementary Data 8) disrupts the clustering by the number of significant DEGs (Supplementary Fig. 7). Since the clustering is still dominated by cell-line-selective effects, the removal did not markedly improve the drug-selective clustering efficiencies (Supplementary Fig. 8).We ranked the remaining eigenarrays by their ability to separate samples of each drug or cell line from all other samples (Supplementary Fig. 9A). Our algorithm allowed a clear separation of cell-line and drug-selective effects (Supplementary Fig. 9B). For each drug, we combined unique sets of drug-selective eigenarrays to drug-selective subspaces with the goal of optimizing the drug’s clustering efficiency (Supplementary Fig. 10A), while preserving a maximum of the original information, as quantified by cosine similarity. Changing the relative contribution of F1 scores (Supplementary Fig. 10B) and cosine similarities (Supplementary Fig. 10C), we screened the potential subspaces for outlier responses where one cell line showed a significantly different response to the drug of interest than all other cell lines (Supplementary Figs. 10D, 11, 12). Such outlier responses were linked to our genomic analysis at a later stage. For 24 drugs, we could identify subspaces with outlier responses (Supplementary Fig. 10E), for the other 30 drugs, we selected subspaces based on high relative contribution of the F1 score (95%).Drug-selective gene expression profiles are similar for drugs with similar mechanismsProjection of gene expression profiles into the final drug-selective subspaces generated drug-selective gene expression profiles (Supplementary Data 9) (averaged DEGs in Supplementary Data 10) with high clustering efficiencies (Fig. 1D, Supplementary Fig. 10F). Drug-selective DEGs were combined into a new matrix. Clustering of this matrix identifies clusters that mostly contained samples that were treated with the same drug and documented preservation of outlier characteristics (Fig. 1E, Supplementary Fig. 13). Twelve of the 24 identified outlier samples are clearly separated from the cell lines treated with the same drug (closed circles in Fig. 1E and Supplementary Fig. 13). Another seven outlier samples are grouped together with the samples treated with the same drug in a larger cluster (open circles in Supplementary Fig. 13), but are separated from those, once that cluster is subclustered. Only four identified outliers reside in the same cluster (open rectangles) as the other cell lines treated with the same drug. Comparison of Figs. 1B and 1E in terms of the breadth of the colored bars allows for visualization of this drug-selective clustering.Our algorithm identified drug-selective gene expression profiles for each drug independently. Nevertheless, drugs with similar mechanisms and overlapping targets still cluster together (Fig. 1E, Supplementary Fig. 13), indicating the potential biological validity of our approach. All eight TKIs targeting EGFR signaling, either by inhibiting the receptor (EGFR, ERBB2) or its intracellular ERK signaling cascade (MAP2K1, MAP2K2)17 are part of the same cluster. One of the four non-TKI drugs in this cluster, phenylephrine, and isoprenaline which is part of a slightly bigger cluster, stimulate adrenergic signaling that can cross-activate ERK signaling in the heart18,19, as documented for phenylephrine in the perfused rat heart19. Two additional non-TKI drugs in this cluster, verapamil and amiodarone, influence adrenergic signaling as antagonists20,21,22 and verapamil has been shown to antagonize ERK signaling in rat cardiomyocytes23. Similarly, the JAK inhibitors ruxolitinib and tofacitinib and the proteasome inhibitors bortezomib and carflizomib are part of two independent clusters that contain no other drugs. Two of three antidiabetics, rosiglitazone and saxagliptin are part of a smaller cluster that first merges with a cluster containing decitabine, a drug with hyperglycemia as a major side effect24 and then with a cluster containing cyclosporine and olmesartan. Cyclosporine increases insulin resistance in humans25 and olmesartan −/+ amlodipine ameliorates it in rats26 and patients27, respectively. The anthracyclines epirubicin, daunorubicin and doxorubicin are part of larger cluster, while the anthracycline idarubicin clusters together with their identified outlier samples.In summary, our SVD pipeline uncovers similar expression profiles induced by similar drugs. Since it does not actively search for such similarities, but only aims at maximizing the overlap between DEGs induced by the same drug, this finding serves as internal validation for the extraction of drug-selective components by our computational pipeline.The grouping of drugs with similar mechanisms into higher-level clusters explains eleven non-outlier samples that are separated from the other samples treated with the same drug (closed rhombuses), leaving only four samples that are not part of drug-selective clusters for no obvious biological reason (open rhombuses) (Supplementary Fig. 13).SVD-based resolution of drug-dependent DEGs enables identification of subcellular processes relevant to TKI-induced cardiotoxicityComplete and decomposed DEGs were subjected to pathway enrichment analysis using the Molecular Biology of the Cell Ontology (MBCO)28 (Supplementary Fig. 14A). MBCO subcellular processes (SCPs) are organized in three to four levels where higher-level SCPs (i.e., levels with lower numbers) describe more general, and lower-level SCPs more detailed cell biological processes. Predicted up- or downregulated SCPs of each drug/cell line combination and SCP level were ranked by significance (Supplementary Data 11, 12 and 13). We refer to these ranks as enrichment ranks in the following.Overall comparison of the enrichment results before and after decomposition documents a large increase in the number of overlapping SCPs in different cell lines treated with the same drug (Supplementary Figs. 14B, C, D, E, 15A).In addition, our SVD decomposition allowed identification of SCPs that were masked in the complete dataset. For example, we could document downregulation of SCPs involved in single protein degradation as an effect selective for anthracyclines (Supplementary Fig. 15B), results that we did not obtain using the complete DEGs (Supplementary Fig. 15C). For a discussion of these SCPs, their potential link to supportive treatments and examples that demonstrate agreement with prior knowledge from small-scale experiments see Supplementary Note 2.Subcellular processes and cell types indicative of cardiotoxic response to TKI treatmentTo identify signatures associated with TIC we searched for SCPs that are almost exclusively up- or downregulated at higher enrichment ranks by cardiotoxic TKIs as compared to non-cardiotoxic TKIs, using F1 score and Area under the Curve (AUC) statistics (Fig. 2A, Supplementary Figs. 16 and 17, Supplementary Data 14). SCPs were ranked by decreasing AUC. The top 10, 10, 25 and 10 level-1, -2, -3 and -4 SCPs were grouped by functional similarities (Fig. 2B, Supplementary Fig. 18) and integrated into the MBCO hierarchy (Fig. 2C, Supplementary Fig. 19). In addition, we mapped identified SCPs back to the cardiotoxic TKIs that up- or downregulate them (Supplementary Fig. 20).Fig. 2: Potential subcellular processes indicative of TKI-induced cardiotoxicity.Up- and downregulated genes among the top 600 drug-selective gene expression profiles were subjected to pathway enrichment analysis using MBCO and Fisher’s Exact test. Significantly up- or downregulated SCPs (nominal p-value ≤ 0.05) were ranked separately by significance for each sample and SCP level. A To screen for SCPs that are selectively induced or repressed by cardiotoxic TKIs, we calculated how many cardiotoxic and non-cardiotoxic TKIs upregulate an SCP of interest in any cell line at any rank cutoff from 1 to 30. Definition of cardiotoxic TKIs as true positives allowed calculation of an F1 score (beta = 0.25) at each analyzed enrichment rank and quantification of the Area under the Curve (AUC). Similarly, we calculated F1 scores and an AUC by analyzing TKIs that downregulate the same SCP. To filter for mixed effects, we subtracted half of the other AUC from each AUC. SCPs were ranked by decreasing AUCs. Flow chart is used with permission from Mount Sinai Health System, licensed under CC BY. B Top up- (red) or downregulated (dark blue) 25 level-3 SCPs predicted for the cardiotoxic TKIs were grouped based on the higher-level functions. White numbers indicate AUC ranks. C The same analysis was applied to level-1, -2 and -3 SCPs, except that we only focused on the AUC obtained for enrichment ranks 1 to 20, due to a smaller number of SCPs within these levels. We also repeated the whole analysis, screening for SCPs selectively induced or repressed by non-cardiotoxic TKIs. Identified SCPs up- and downregulated for cardiotoxic (red and dark blue, respectively) and non-cardiotoxic TKIs (light blue and orange, respectively) for all levels were integrated into the MBCO hierarchy. Selected branches are shown. See Supplementary Fig. 19 for all predictions.Besides being elevated or decreased above or below a threshold that determines an SCP’s association with a cardiotoxic response, an SCP’s activity might already be beyond that threshold at baseline level, i.e., before TKI treatment. Candidates for this group could be SCPs that are regulated by non-cardiotoxic TKIs. A non-cardiotoxic TKI might change the baseline SCP activity from sufficient to insufficient for a cardiotoxic response. Consequently, our second focus was the identification of SCPs that non-cardiotoxic TKIs would up- or down-regulate to potentially suppress cardiotoxicity (Supplementary Figs. 18, 19, 21).Single cell RNA sequencing (RNAseq) of four of our six cardiomyocyte cell lines identified at least two major subtypes mapping to adult cardiomyocytes and an epicardium-derived cell type that is similar to fibroblasts16 (Supplementary Fig. 22A–D, Supplementary Data 15, 16A, B). The SCPs we identified by bulk transcriptomics often describe canonical functions of cell types other than cardiomyocytes. We determined how the top SCPs in our data map to our subtypes (Supplementary Fig. 23, Supplementary Data 16C–F), and to the major cell types of the adult human heart analyzed by others using single cell transcriptomics29 (Fig. 3A, Supplementary Fig. 23, Supplementary Data 17, 18A–D).Fig. 3: SCPs can be mapped to cellular subtypes and known cardiomyopathy disease mechanisms.A We subjected marker genes for ventricular and atrial cardiomyocytes (VCM, red fields, and ACM, turquoise fields, respectively), cardiac fibroblasts (CFB, brown fields) and smooth muscle cells (SMC, brown fields) obtained from single nucleus RNAseq of the adult human heart29 to pathway enrichment analysis using MBCO and Fisher’s exact test. Significant SCPs of each cell type (nominal p-value ≤ 0.05) were ranked by significance (numbers in the diagram). Names of SCPs whose higher and lower activities favor a cardiotoxic response are colored red and blue, respectively. B DEGs in heart cells obtained by single cell (SC)13 or nucleus (SN)14 RNAseq from patients with DCM or HCM as well as in hiPSC-derived cardiomyocytes obtained from an infant patient with DCM32 were subjected to pathway enrichment analysis using MBCO and Fisher’s exact test. Significantly up- (red fields) or downregulated (blue fields) SCPs of each cell type (nominal p-value ≤ 0.05) were ranked by significance (numbers in the diagram). Names of SCPs are colored as described in (A).Many of the identified SCPs described well-known biology involved in inherited and acquired cardiomyopathies that develop independently of chemotherapy and map to multiple cell types of the adult human heart. Our results suggest that cardiotoxic TKIs might mimic these mechanisms by early (within 48 hours) transcriptional regulation.Four SCPs involved in muscle contraction and sarcomere renewal were simultaneously identified for the cardiotoxic and non-cardiotoxic drugs (Fig. 2B, Supplementary Fig. 18). Another two SCPs were identified for each group independently. These SCPs that span all four SCP levels (Fig. 2C, Supplementary Fig. 19) are regulated by the two TKI toxicity groups in opposite directions. The difference in directionality suggests insufficient sarcomere renewal as relevant for TIC. All cardiotoxic TKIs, except vandetinib and trastuzumab downregulate (Supplementary Fig. 20), while 13 of 17 non-cardiotoxic TKIs upregulate (Supplementary Fig. 21) at least one of these SCPs, respectively. The SCPs map to cardiomyocyte clusters identified by single cell RNAseq of our cell lines and of the adult human heart29 (Fig. 3A, Supplementary Fig. 23). In agreement with our data, proper heart functioning depends on a continuous turnover of cardiomyocyte sarcomere proteins30. About half of hypertrophic and dilated cardiomyopathies are associated with mutations in sarcomeric proteins31. Further supporting this line of reasoning, we found downregulation of muscle contraction-related SCPs in single cell RNAseq from hiPSC-derived cardiomyocytes from an infant with DCM32 (Supplementary Data 19A–C) and from adult cardiomyocytes13,14 (Supplementary Data 20A, B) (Fig. 3B, Supplementary Fig. 24). These SCPs might serve as a starting point for supportive therapy. Four of the six regulated genes of the SCP ‘Thin myofilament organization’ (Supplementary Data 11C) are involved in blocking the binding of the myosin head to the thin myofilament during muscle contraction. This mechanism is targeted by the new drug mavacamten that was recently approved by the FDA to treat obstructive Hypertrophic Cardiomyopathy (HCM)33.Our data associates the upregulation of SCPs involved in the citric acid cycle and mitochondrial energy generation with a cardiotoxic response (Figs. 2B, C, Supplementary Figs. 18, 19). They are almost exclusively upregulated by pazopanib (Supplementary Fig. 20), a TKI with a high rate of cardiotoxicity (>10%) (Supplementary Data 3). In contrast, many studies document an overall reduction in oxidative phosphorylation during heart failure34 and identified SCPs are downregulated in adult Dilated Cardiomyopathy (DCM)13,14 and HCM13,14 cardiomyocytes as well as in hiPSC-derived DCM cardiomyocytes32, as predicted from single cell and nucleus RNAseq datasets (Fig. 3B, Supplementary Fig. 24, Supplementary Data 21A, B, C). In support of our findings, compensatory upregulation of oxidative phosphorylation was suggested for patients belonging to a large DCM subgroup that is caused by truncating titin variants35,36,37.Increasing evidence links ferroptosis, an iron-dependent accumulation of lipid peroxides that triggers cell death, to heart failure38. Some identified SCPs in our data might converge on ferroptosis as an endpoint. The SCP ‘Desaturation of fatty acids’ was associated with a cardiotoxic response (Fig. 2B, C, Supplementary Figs. 18, 19) due to its upregulation by pazopanib (Supplementary Fig. 20). Upregulated genes mapping to this SCP (Supplementary Data 11C) can generate polyunsaturated fatty acids (PUVAs)39,40 that are precursors of lipid peroxides38. A lower activity of the level-2 SCP ‘Cellular antioxidant systems’ and its level-1 parent ‘Cellular redox homeostasis’, pathways offering protection against ferroptosis38, is associated with a cardiotoxic response (Supplementary Figs. 18 and 19). They were downregulated by the cardiotoxic TKIs vandetinib and bevacizumab (Supplementary Fig. 20). In addition, downregulation of ‘Cellular iron storage’ by the anthracyclines doxorubicin and daunorubicin identified by us (Supplementary Figs. 14D, 15B) and others41 suggests stimulation of ferroptosis by an increase of intracellular free iron41. In agreement, supportive therapy with iron chelators protects against AIC42, though the protective effect might involve mechanisms that are unrelated to iron homeostasis43.Potential clinical relevance of our findings is additionally indicated by identification of SCPs involved in cholesterol metabolism and natriuretic peptide signaling (Fig. 2B, Supplementary Fig. 18). Our SCPs agree with recommended treatment schemes of cardiomyopathy and drug-induced cardiomyopathy involving statins44 and neprilysin45,46, respectively. Prediction of extracellular collagen-crosslinking (Supplementary Fig. 18) can be linked to histopathological observations that correlate with disease progression47,48. Involvement of identified signaling pathways in TIC i.e., PDGF, HIF-1alpha, Oncostatin M and Hippo signaling (Fig. 2B, Supplementary Fig. 18) is supported by their involvement in drug-independent cardiomyopathy. These and other findings regarding the predicted SCPs are discussed in further detail in Supplementary Note 2.Taken together the data from our experiments, when integrated with data from the literature, indicate the TKIs are likely to have their cardiotoxic effects by regulating SCPs in multiple cell types of the heart and mimicking mechanisms involved in drug-independent cardiomyopathy and heart failure.Subcellular processes associated with anthracycline toxicityTo document which pathway activities are specifically associated with AIC in distinction to the cardiotoxicity induced by very cardiotoxic TKIs (frequency >10%), we used the F1 score and AUC algorithm to compare both groups of drugs (Supplementary Fig. 25). Results were integrated into the MBCO hierarchy (Supplementary Fig. 26). As already described, the observed downregulation of ‘Cellular iron storage’ agrees with the central role of ferroptosis in AIC41. Observed upregulation of ‘Nucleotide excision repair’ agrees with the documented involvement of both nucleotide excision repair and homologous recombination in recognition and repair of anthracycline-DNA adducts49. Upregulation of histone genes and other genes involved in chromatin remodeling could be the consequence of anthracycline-induced histone eviction50. As suggested by our data, dysregulation of mitochondrial dynamics is another mechanism involved in AIC41. Similarly, dysregulation of proteasomal degradation and autophagy has been reported previously51, as well as activation of cell death pathways52.TKI transcriptomic signatures are similar in cardiomyoctes cocultured with and without endothelial cellsOur experimental protocol raises the possibility that the signatures we obtain are the consequence of our single cell type culture conditions that are missing other cell types of the human heart. One cell type thought to be relevant for cardiotoxicity is the endothelial cell type of the blood vessels53. To determine if our signatures in cardiomyocytes may be influenced by other cell types, we repeated the stimulation of two selected cell lines with pazopanib and dabrafenib in absence or presence of human coronary arterial endothelial cells (HCAECs) (Supplementary Fig. 27A, Supplementary Data 2B). Cardiomyocytes in the presence of endothelial cells showed increased beating rates. Drug treatment at therapeutic doses did not affect the morphology or viability of either endothelial cells (Supplementary Fig. 27B) or cardiomyocytes. Cardiomyocytes and endothelial cells were cocultured in coculture plates with porous membranes for 48 hours and then treated with drugs for another 48 hours. This set up allows to obtain the cardiomyocytes separately for extraction and bulk transcriptomics under the same conditions used for the original experiments.Pazopanib is a multi-tyrosine kinase inhibitor whose antineoplastic effect is associated with inhibition of angiogenesis and interference with vascular endothelial growth factor (VEGF) and basic fibroblast growth factor (bFGF) signaling in human endothelial cells54. Endothelial dysfunction caused by VEGF blockage is discussed as a potential mechanism involved in pazopanib-induced cardiomyopathy55, suggesting that pazopanib’s high cardiotoxicity (>10%, Supplementary Data 3) could at least partly depend on interactions between endothelial cells and cardiomyocytes. As a second drug, we selected the B-Raf inhibitor dabrafenib, as a representative for TKIs with less evidence for cross-cellular effects56,57.Eight lists of DEGs of the new data (Supplementary Data 22A) were merged with the original 266 lists (Supplementary Data 4). The merged 274 lists of DEGs were projected into the drug-selective subspaces identified from the original 266 lists (Supplementary Data 22B, C), followed by pairwise correlation analysis and hierarchical clustering. Both, pazopanib and dabrafenib (Supplementary Fig. 28A) treated cell lines clustered closely together and there was only a slight change in the overall F1 scores (Supplementary Fig. 28B). A distinct clustering of DEGs generated in presence or absence of the endothelial cells was not observed for either TKI. Clustering of the new combined final drug-selective DEG matrix that now contained 274 instead of 266 lists of DEGs showed only minor rearrangements within a larger cluster of nine drugs that contained pazopanib and dabrafenib (Supplementary Fig. 28C and D). Predicted pathways (Supplementary Fig. 29, Supplementary Data 23–25) and cardiotoxic pathways induced or repressed by pazopanib (Supplementary Fig. 30A) and dabrafenib (Supplementary Fig. 30B) showed only minor differences between the new and old datasets. Overall, this analysis suggests that the effect of endothelial cocultures on TKI-induced signatures in cardiomyoctyes is minimal.Candidates for genomic variants associated with drug-induced cardiotoxicityUsing whole genome sequencing of the six cell lines16, our next analysis step focused on the identification of potential genomic variants that might be associated with a higher or lower risk for drug-induced cardiotoxicity (DIC) or more specifically TIC (Fig. 4A, Supplementary Fig. 31). Cardiotoxic responses to TKI treatment are observed in less than 10% to 20% of the treated patients (Supplementary Data 3). Consequently, we hypothesized that the population-wide frequency of a potential monogenic allele associated with DIC should be in a similar range. Relevance for cardiac tissue was assumed for those variants that map to gene coding regions or are part of cis-expression (e) or -splicing (s) QTLs in the adult human heart58. We hypothesized that potentially cardiotoxic or -protective variants could interfere with a drug’s pharmacokinetics (PK) or -dynamics (PD) or with pathways that are targeted by that drug.Fig. 4: Identification of genomic variants that are potentially associated with a cardiotoxic response.Whole genome sequencing of our six cell lines16 was used to identify alleles in our cell lines at known variant positions. A See text and methods for details of our pipeline for identification of potential genomic variants involved in PK/PD or induced SCP activities for a drug of interest. Flow chart is used with permission from Mount Sinai Health System, licensed under CC BY. B Clustering of DEGs within the daunorubicin-selective subspace reveals an outlier response in cell line MSN09 after daunorubicin treatment. C The identified variant rs2229774 in cell line MSN09 maps to the coding region of the transcription factor RARG regulating the expression of TOP2B and ABCB8, both involved in PK/PD of daunorubicin and doxorubicin that induce an outlier response in this cell line. D Enrichment significance ranks for “WNT-Beta-catenin signaling pathway” obtained by analysis of upregulated genes after daunorubicin or doxorubicin treatment of indicated cell lines. The cell line MSN09 (purple) contains the rs2229774 mutation in the RARG gene. Field colors change from bright to dark yellow with increasing ranks. ‘>’ indicates that the SCP was not predicted or predicted with a rank > 99. E In total, we identified 213 and 201 potential variants associated with TIC or AIC by interference with PK/PD mechanisms, respectively. Variants mapping to multiple gene classes are split equally among them to prevent double counting. Drug names are colored according to their class (orange: Small molecule kinase inhibitors (KI), red: monoclonal antibodies against KIs, purple: anthracyclines). F We compared the overlap of identified SCP genes associated with a cardiotoxic or non-cardiotoxic response to genes associated with inherited DCM or HCM, either within the HuGE Phenopedia database or identified in GWAS. TWAS: transcription-wide-association studies. G Variants that meet our population-wide criteria were mapped to up- and downregulated level -2, -3 and -4 SCPs that we predicted as indicative for TIC. Variants that map to identified SCPs of multiple levels are only counted for the lowest level SCPs (higher level numbers) to prevent double counting. Drug names are colored as described in (E). H Up- and downregulated SCPs associated with a cardiotoxic response were mapped back to the cardiotoxic TKIs that induce them. Numbers in brackets show identified variants for each SCP gene. Blue indicates that the SCP is a level-3 SCP.Variants of the first interference group should induce a different transcriptomic response, allowing us to link cell-line-specific variants to identified 24 outlier responses. Since, based on statistical likelihoods, either none or only one of our six healthy volunteers should suffer from cardiotoxicity induced by a drug of interest, we focused on variants that are overexpressed in an outlier cell line for a drug of interest, if compared to the other five lines. Mechanistic knowledge is added by considering only variants that map to genes with potential involvement in a drug’s PK/PD.Variants of the second group should map to pathways that either contain drug target proteins or are induced or repressed by the drug of interest on the transcriptional level. Using the results of our transcriptomic analysis, we focused on the second subgroup with a particular interest in those SCPs that are associated with a cardiotoxic or non-cardiotoxic response.Identification of genomic variants interfering with a drug’s PK/PDBased on outlier responses to daunorubicin and doxorubicin (Fig. 4B, Supplementary Figs. 11 and 12), our PK/PD algorithm allowed re-identification of the variant rs2229774 (Fig. 4C, Supplementary Data 25) within the coding region of the RARG gene15, one out of three genomic variants with the strongest evidence for AIC59. The suppressive activity of the transcription factor RARG on the expression of the anthracycline target TOP2B is reduced by this variant, leading to increased AIC60. Pathway enrichment analysis suggests missing upregulation of WNT receptor signaling as a potential mechanism of rs2229774-triggered AIC (Fig. 4D, Supplementary Fig. 14D). In agreement, inhibition, or activation of WNT signaling in mice and cell assays enhances or mitigates AIC, respectively61,62,63.In total, our outlier-based screening approach identified 128 and 111 potential variants that map to 100 and 88 genes involved in the PK/PD of three anthracyclines and three cardiotoxic TKIs, respectively (Fig. 4E, Supplementary Fig. 32A, Supplementary Data 26). Identified variant candidates could trigger or protect from AIC or TIC. Three of the variants map to drug target proteins of ponatinib and trametinib (Supplementary Fig. 32B), TKIs with medium (1–10%) and high ( > 10%) rates of TIC, respectively (Supplementary Data 3).Considering all 24 drugs with documented outlier responses for this analysis we identified 464 variants (Supplementary Fig. 32A) mapping to 288 genes (Supplementary Fig. 32C) (Supplementary Data 26).SCP-based identification of genomic variants associated with TKI-induced cardiotoxicityIdentified level-2, -3 and -4 SCPs associated with the response to cardiotoxic and non-cardiotoxic TKIs enrich for genes that are mapped to inherited DCM or HCM, either within the HuGE Phenopedia database64 or obtained from GWAS65,66 (Fig. 4F). We left out level-1 SCPs from this analysis, since they were predicted based on genes mapping to a sub-function of the general cellular functions the level-1 SCPs describe. It needs to be considered that some of these variants could contribute to drug cardiotoxicity by modulating propensity for cardiac insufficiency without primarily effecting cardiomyocytes. Nevertheless, the significant overlap supports our hypothesis that TKIs induce cardiotoxicity by mimicking mechanisms involved in inherited or acquired drug-independent cardiomyopathy, as discussed above and in Supplementary Note 2.To identify a second set of potential variants interfering with TIC we mapped all variants that met our population-wide requirements to identified SCPs. In total, we identified 634, 410 and 80 non-overlapping variants mapping to level-2, -3 and -4 SCPs that are up- or downregulated by cardiotoxic TKIs (Fig. 4G) (see Supplementary Fig. 32D for variants mapping to SCPs regulated by non-cardiotoxic TKIs). Mapping identified SCPs back to the drugs that induce them (Supplementary Fig. 20) enables identification of variants that might interfere with cardiotoxicity of a drug of interest (Supplementary Data 27). For example, we predicted 13 variants mapping to four genes of the SCP ‘Thin myofilament organization’ (Fig. 4H) that might interfere with the cardiac response to sunitinib and ponatinib downregulating this SCP in six (enrichment ranks 2 × 3, 4, 2 × 7) and five (5 × 1) of six treated cell lines, respectively (Supplementary Fig. 20), but also to lapatinib, dabrafenib and bevacizumab. Two of those genes are involved in the same function that is targeted by mavacamten, as described above.Since most variants identified in GWAS of drug-independent cardiomyopathy and mapping to our SCPs occur at high frequency and are not part of cis e/s-QTLs, they do not meet our population-wide requirements. Hence, we do not list them as being involved in drug-induced cardiotoxicity. Some of these variants might not primarily affect cardiomycoytes or other abundant cells in the heart. Additionally, HCM and DCM exist as both monogenic and polygenic diseases. Similarly, single or multiple variants might be related to TIC, requiring development of more sophisticated genomic algorithms. Nevertheless, identified SCPs are a good starting point to reduce the burden of multiple testing hypothesis by restricting the focus on genes with a functional implication in TIC.
