Abstract Knowledge about synthetic lethality can be applied to enhance the efficacy of anticancer therapies in individual patients harboring genetic alterations in their cancer that specifically render it vulnerable. We investigated the potential for high-resolution phenomic analysis in yeast to predict such genetic vulnerabilities by systematic, comprehensive, and quantitative assessment of drug–gene interaction for gemcitabine and cytarabine, substrates of deoxycytidine kinase that have similar molecular structures yet distinct antitumor efficacy. Human deoxycytidine kinase (dCK) was conditionally expressed in the Saccharomyces cerevisiae genomic library of knockout and knockdown (YKO/KD) strains, to globally and quantitatively characterize differential drug–gene interaction for gemcitabine and cytarabine. Pathway enrichment analysis revealed that autophagy, histone modification, chromatin remodeling, and apoptosis-related processes influence gemcitabine specifically, while drug–gene interaction specific to cytarabine was less enriched in gene ontology. Processes having influence over both drugs were DNA repair and integrity checkpoints and vesicle transport and fusion. Non-gene ontology (GO)-enriched genes were also informative. Yeast phenomic and cancer cell line pharmacogenomics data were integrated to identify yeast–human homologs with correlated differential gene expression and drug efficacy, thus providing a unique resource to predict whether differential gene expression observed in cancer genetic profiles are causal in tumor-specific responses to cytotoxic agents. Keywords: yeast phenomics, gene–drug interaction, genetic buffering, quantitative high throughput cell array phenotyping (Q-HTCP), cell proliferation parameters (CPPs), gemcitabine, cytarabine, recursive expectation-maximization clustering (REMc), pharmacogenomics 1. Introduction Genomics has enabled targeted therapy aimed at cancer driver genes and oncogenic addiction [[40]1], yet traditional cytotoxic chemotherapeutic agents remain among the most widely used and efficacious anticancer therapies [[41]2]. Changes in the genetic network underlying cancer can produce vulnerabilities to cytotoxic chemotherapy that further influence the therapeutic window and provide additional insight into their mechanisms of action [[42]3,[43]4]. A potential advantage of so-called synthetic lethality-based treatment strategies is that they could have efficacy against passenger gene mutation or compensatory gene expression, while classic targeted therapies are directed primarily at driver genes ([44]Figure 1A). Quantitative high throughput cell array phenotyping of the yeast knockout and knockdown libraries provides a phenomic means for systems level, high-resolution modeling of gene interaction [[45]5,[46]6,[47]7,[48]8,[49]9], which is applied here to predict cancer-relevant drug–gene interaction through integration with cancer pharmacogenomics resources ([50]Figure 1B). Figure 1. [51]Figure 1 [52]Open in a new tab Experimental model of gemcitabine and cytarabine drug–gene interaction networks. (A) The strategy of cytotoxic anticancer drug–gene interaction is illustrated in the context of driver gene-mediated oncogenesis. Driver genes promote cancer and influence the expression of passenger genes (black arrows), which also leads to genomic instability and alterations in the genetic buffering network. The genetic buffering network (blue arrows) maintains cellular homeostasis and is altered in cancer cells by genomic instability, thereby creating the potential for drug–gene interaction that increases the therapeutic window of anticancer agents (red arrows). Drug–gene interaction can either involve driver or passenger genes directly, or the compromised genetic buffering network, which are systematically characterized by the quantitative yeast phenomic model. (B) The synthetic genetic array (SGA) method was used to introduce tet-inducible human deoxycytidine kinase (dCK) expression in the yeast knockout and knockdown (YKO/KD) collection. The phenomic model incorporates treatment of individually grown cultures of the YKO/KD collection, and 768 replicate reference (Ref) strain cultures, with increasing gemcitabine (0, 5, 10, 20, and 30 μg/mL) or cytarabine (0, 10, 25, 50, and 100 μg/mL) in a dextrose (HLD) media, with dCK induced by addition of doxycycline. Drug–gene interaction profiles were analyzed by recursive expectation-maximization clustering (REMc) and gene ontology (GO) term analysis to characterize phenomic modules with respect to drug–gene interaction for gemcitabine or cytarabine, and integrated with pharmacogenomics data to predict evolutionarily conserved drug–gene interactions relevant to precision oncology. (C) Structures and metabolism of deoxycytidine analogs. Nucleoside analogs include a diverse group of compounds with anticancer, antiviral, and immunosuppressive efficacy [[53]10]. The anticancer agents have tissue-specific efficacy ranging from solid tumors to leukemias, yet details about how these agents confer differential activity are unknown [[54]10,[55]11]. Gemcitabine (2′,2′-difluoro 2′-deoxycytidine, dFdC) and cytarabine (Ara-C) are deoxycytidine analogs that undergo the first step of conversion to their active triphosphate forms by deoxycytidine kinase (dCK) ([56]Figure 1C). The nucleoside triphosphate analogs can be incorporated into DNA and inhibit the functions of polymerases and other enzymes involved in DNA metabolism. For example, gemcitabine inhibits ribonucleotide reductase (RNR), which limits the production of deoxyribonucleotides (dNTPs) that are needed for DNA synthesis and repair [[57]11]. Gemcitabine has been used as a single agent in the treatment of some cancers, such as pancreatic, and in combination with platinum-based drugs in non-small-cell lung, breast, and ovarian cancers [[58]12,[59]13,[60]14,[61]15]. Cytarabine, on the other hand, has been an important agent in treatments for acute myeloid leukemia and acute lymphoblastic leukemia [[62]16]. Deoxycytidine kinase (dCK) phosphorylates deoxycytidine to deoxycytidine monophosphate (dCMP), similarly phosphorylating gemcitabine and cytarabine to dFdCMP and AraCMP, respectively. UMP/CMP kinase and the nucleoside diphosphate kinase are subsequently involved in conversion to the triphosphate form ([63]Figure 1C). Reduced expression of dCK or high expression of RNR subunits RRM1 and RRM2 is associated with increased gemcitabine resistance [[64]10,[65]12,[66]17,[67]18,[68]19,[69]20,[70]21]. Genomic analyses have suggested genetic influences on the efficacy of gemcitabine or cytarabine [[71]22,[72]23,[73]24,[74]25,[75]26], which we model here at a systems level by surveying gene–drug interaction to elucidate biology underlying differential anticancer efficacies of the respective drugs, and thereby aid in predicting treatment outcomes based on individual patient cancer genetic profiles. Saccharomyces cerevisiae does not have a dCK homolog and is thus naturally resistant to gemcitabine and cytarabine. To examine the gene–drug interaction networks for gemcitabine and cytarabine in yeast, we introduced human dCK into the yeast knockout and knockdown (YKO/KD) library by the synthetic genetic array (SGA) method [[76]27,[77]28,[78]29] and conducted phenomic analysis on the resulting double mutant library by quantitative high-throughput cell array phenotyping (Q-HTCP) [[79]6,[80]7,[81]8,[82]9] using multiple growth inhibitory concentrations of gemcitabine or cytarabine ([83]Figure 1B). Cell proliferation parameters (CPPs) obtained by Q-HTCP were used to quantify and compare drug–gene interaction for gemcitabine vs. cytarabine. The unbiased results provide a systems-level resource of genetic and biological information about the cytotoxicity of these drugs, incorporating knowledge about genes that either buffer or promote their effects [[84]3,[85]5]. Recent advances in cancer pharmacogenomics have provided gene expression and drug sensitivity data from hundreds of cancer cell lines, establishing associations between gene expression and anticancer efficacy for many compounds, including gemcitabine and cytarabine [[86]30,[87]31,[88]32]. We investigated the potential utility of a yeast phenomic model of chemotherapy sensitivity and resistance for predicting causality in correlations between differential gene expression and drug sensitivity by generating a network-level drug–gene interaction resource. The resource integrates cancer pharmacogenomic and yeast phenomic data, using the results to query the cancer genetics literature in order to obtain systems-level biological insights about how yeast phenomic models help predict cytotoxic chemotherapy efficacy based on unique genetic alterations specific to each individual patient’s cancer ([89]Figure 1A). 2. Materials and Methods 2.1. Strains, Media, and Drugs We obtained the yeast gene knockout strain library (YKO) from Research Genetics (Huntsville, AL, USA) and the knockdown (KD) collection, also referred to as the decreased abundance of mRNA production (DAmP) library, from Open Biosystems (Huntsville, AL, USA). The YKO library is in the genetic background of BY4741 (S288C MATa ura3-∆0 his3-∆1 leu2-∆0 met17-∆0). Additional information and strains can be obtained at [90]https://dharmacon.horizondiscovery.com/cdnas-and-orfs/non-mammalian -cdnas-and-orfs/yeast/#all. Some mutants appear multiple times in the library and they are treated independently in our analysis. HLD is a modified synthetic complete medium [[91]8] and was used with 2% dextrose (HLD) as the carbon source. Doxycycline hydrochloride (BP26535) was obtained from Fisher Scientific. Gemcitabine (Gemzar) was obtained from Eli Lilly and Company (0002-7502-01). Cytarabine was obtained from Bedford Laboratories (55390-131-10). A tet-inducible dCK query allele was constructed in the SGA background in the following way: An integrating plasmid for doxycycline-inducible gene expression was constructed by subcloning 3′UTR and 5′ORF targeting sequences from the LYP1 locus into pJH023 [[92]33], creating pJH023_UO_lyp1, and the reverse VP16 transactivator (Tet-ON), obtained by PCR from pCM176 [[93]34], was fused to the ACT1 promoter by overlap PCR and subcloned into pJH023_UO_lyp1, replacing the VP16 transactivator (Tet-OFF) and creating the “Tet-ON” construct, pML1055 [[94]35]. pML1055 was digested with NOT1 and transformed into strain 15578-1.2b_LYP1 (MATα his3∆1 leu2∆0 ura3∆0 can1∆0::P[GAL1]-T[ADH1]-P[MFA1]-his5^+[sp] hmr∆0::URA3ca), which was derived by backcrossing 15578-1.2b (MATα his3∆1 leu2∆0 ura3∆0 can1∆0::P[GAL1]-T[ADH1]-P[MFA1]-his5^+[sp] lyp1∆0 hmr∆0::URA3ca) to restore the LYP1 locus. The resulting chromosomal integration of pML1055 between the promoter and ORF at the LYP1 locus was selected with nourseothricin, giving rise to DWY1 (MATα his3∆1 leu2∆0 ura3∆0 can1∆0::P[GAL1]-T[ADH1]-P[MFA1]-his5^+[sp] hmr∆0::URA3ca Pact1-revTetR-VP16-natMX-PtetO7-LYP1). Tet-inducible LYP1 in DWY1 was verified phenotypically by doxycycline-dependent SAEC sensitivity [[95]35]. Overlap PCR was performed to fuse deoxycytidine kinase (from a plasmid, gift of Bo Xu and William Parker, Southern Research) and the HPH gene (from pFA6a-HBH-hphMX4) [[96]36], introducing flanking sequences for replacement of the LYP1 ORF (see [97]Additional File 1, Table S1 for primers). The PCR product was transformed into DWY1 ([98]Additional File 2, Figure S1) and transformants selected on hygromycin were confirmed by doxycycline-induced sensitivity to gemcitabine and cytarabine, yielding MIY16 (MATα his3∆1 leu2∆0 ura3∆0 can1∆0::P[GAL1]-T[ADH1]-P[MFA1]-his5^+[sp] hmr∆0::URA3ca lyp1-∆0::Pact1-revTetR-VP16-natMX-PtetO7-DCK). The synthetic genetic array (SGA) method, a way to introduce an allele of interest into the YKO/KD library and recover haploid double mutants [[99]28,[100]29], was used to derive a haploid YKO/KD collection with doxycycline-inducible dCK expression. 2.2. Quantitative High Throughput Cell Array Phenotyping (Q-HTCP) Q-HTCP, an automated method of collecting growth curve phenotypes for the YKO/KD library arrayed onto agar media, was used to obtain phenomic data [[101]37]. A Caliper Sciclone 3000 liquid handling robot was used for cell array printing, integrated with a custom imaging robot (Hartman laboratory) and Cytomat 6001 (Thermo Fisher Scientific, Asheville, NC, USA) incubator. Images of the 384-culture arrays were obtained approximately every 2–3 hours and analyzed as previously described [[102]9,[103]37]. To obtain CPPs, image analysis was performed in Matlab and data were fit to the logistic equation, G(t) = K/(1 + e^−r(t−l)), assuming G(0) < K, where G(t) is the image intensity of a spotted culture vs. time, K is the carrying capacity, r is the maximum specific growth rate, and l is the moment of maximal absolute growth rate, occurring when G(t) = K/2 (the time to reach half of carrying capacity) [[104]7]. The CPPs, primarily K and L, were used as phenotypes to measure drug–gene interaction. 2.3. Quantification of Drug–Gene Interaction Gene interaction was defined as departure of the corresponding YKO/KD strain from its expected phenotypic response to gemcitabine or cytarabine [[105]6,[106]9,[107]38]. The expected phenotype depends, in part, on the difference in cell proliferation phenotypes between the mutant and reference strain without gemcitabine or cytarabine, which we have termed shift. This difference is applied to the entire data series (hence the term, ‘shift’), when assessing differential response in cell proliferation phenotypes between the mutant and reference strain in the presence of escalating concentrations of gemcitabine or cytarabine [[108]5,[109]6,[110]9,[111]33]. The concentrations of gemcitabine or cytarabine (μg/mL) were chosen based on phenotypic responses being functionally discriminating in the parental strain. Gemcitabine, cytarabine, or doxycycline, alone, did not alter cell proliferation ([112]Figure 2C–F; [113]Additional File 2, Figure S2A–D). Thus, shift was calculated in the presence of 5 μg/mL doxycycline without nucleoside analog. Figure 2. [114]Figure 2 [115]Open in a new tab Phenomic analysis of drug–gene interaction for gemcitabine and cytarabine. Average growth curves (from fitting pixel intensity data of 768 replicate cultures to a logistic function) for the reference (RF) strain, treated with the indicated concentrations of (A) gemcitabine or (B) cytarabine. (C–F) Cell proliferation parameter (CPP) distributions from data depicted in panels A and B, also with and without induction of deoxycytidine kinase (0 or 5 μg/mL doxycycline respectively), for (C–D) gemcitabine and (E–F) cytarabine in μg/mL for (C,E) L and (D,F) K. (G,H) Comparison of drug–gene interaction scores calculated for L vs. K for (G) gemcitabine and (H) cytarabine, where score distributions of yeast knockout/knockdown (YKO/KD, black) and non-mutant parental (Ref, red) strain cultures are indicated along with thresholds for deletion enhancement and suppression (dashed lines at +/− 2). (I–J) Differential drug–gene interaction using L (I) or K (J) as the CPP for gemcitabine vs. cytarabine, classified by specificity of gene-drug interaction, where ‘G’, ‘C’, and ‘B’ indicate gemcitabine, cytarabine, or both, respectively. Deletion enhancement or suppression is indicated by ‘_Enh’ or ‘_Sup’. Interaction scores were calculated as previously described [[116]9,[117]39], with slight modifications, as summarized below. All media conditions used for interaction score calculation had 5 μg/mL doxycycline to express dCK. Variables were defined as: * D[i] = concentration (dose) of gemcitabine or cytarabine; * R[i] = observed mean growth parameter for parental reference strain at D[i]; * Y[i] = observed growth parameter for the YKO/KD mutant strain at D[i]; * K[i] = Y[i] − R[i], the difference in growth parameter between the YKO/KD mutant (Y[i]) and reference (R[i]) at D[i]; * K[0] = Y[0] − R[0], the effect of gene KO/KD on the observed phenotype in the absence of gemcitabine or cytarabine—this value is annotated as ‘shift’ and is subtracted from all K[i] to obtain L[i]; * L[i] = K[i] − K[0], the interaction between (specific influence of) the KO/KD mutation on gemcitabine or cytarabine response, at D[i]; For cultures not generating a growth curve, Y[i] = 0 for K and r. Note, however, that for cultures with extremely low growth, the L parameter asymptote is infinity and thus it was assigned Y[i] max, defined as the maximum observed Y[i] among all cultures exhibiting a minimum carrying capacity (K) within 2 standard deviation (SD) of the parental reference strain mean at D[i]. Y[i] max was also assigned to outlier values (i.e., if Y[i] > Y[i] max). Interaction Was Calculated by the Following Steps: * (1) Compute the average value of the 768 reference cultures (R[i]) at each dose (D[i]): * (2) Assign Y[i] max (defined above) if growth curve is observed at D[0], but not at D[i], or if observed Y[i] is greater than Y[i] max. * (3) Calculate K[i] = Y[i] − R[i]. * (4) Calculate L[i] = K[i] − K[0] * (5) Fit data by linear regression (least squares): L[i] = A + B*D[i] * (6) Compute the interaction value ‘INT’ at the max dose: INT = L[i]-max = A + B*D[max] * (7) Calculate the mean and standard deviation of interaction scores for reference strains, mean(REF[INT]) and SD(REF[INT]); mean(REF[INT]) is expected to be approximately zero, with SD(REF[INT]) primarily useful for standardizing against variance ([118]Additional File 1, Tables S2–S5; Additional Files 3–4). * (8) Calculate interaction z-scores: z-score(YKO[INT]) = (YKO[INT] − mean(REF[INT]))/SD(REF[INT]) z-score(YKO[INT]) > 2 for L or < −2 for K are referred to as gene deletion enhancers of gemcitabine or cytarabine cytotoxicity, and conversely, L interaction score < −2 or K interaction scores >2 are considered gene deletion suppressors. Due to the fact that the CPP distributions for KD strains were different from the reference strain, we used the mean and standard deviation from the KD plates only as a conservative measure of variance where z-score(KD[INT]) = (KD[INT] – mean(KD[INT]))/SD(KD[INT]). 2.4. Recursive Expectation-Maximization Clustering (REMc) and Heatmap Generation REMc is a probability-based clustering method and was performed as previously described [[119]40]. Clusters obtained by Weka 3.5, an EM-optimized Gaussian mixture-clustering module, were subjected to hierarchical clustering in R ([120]http://www.r-project.org/) to further aid visualization with heatmaps. REMc was performed using L and K interaction z-scores ([121]Figure 3A). The effect of gene deletion on the CPP (in the absence of drug), termed ‘shift’ (K[0]), was not used for REMc, but was included for visualization in the final hierarchical clustering. [122]Additional File 5, Files A–B contain REMc results in text files with associated data also displayed as heatmaps. In cases where a culture did not grow in the absence of drug, 0.0001 was assigned as the interaction score, so that associated data (‘NA’) could be easily indicated by red coloring in the shift columns of the heatmaps. Figure 3. [123]Figure 3 [124]Open in a new tab Prediction of drug–gene interaction in cancer cells by integration of yeast phenomic and human pharmacogenomic data. Recursive expectation-maximization clustering results were classified visually by their associated gene interaction profiles (see methods). (A) The order of data columns, which is consistent for all heatmaps, is indicated. K-derived interactions are in columns 2 and 4, with L-derived interactions in columns 6 and 8, for gemcitabine and cytarabine, respectively. To the left of each interaction value (indicated by ‘+’), is the corresponding ‘shift’ value (indicated by ‘−‘), referring to the ∆CPP for the respective YKO/KD culture relative to the reference culture average in the absence of gemcitabine or cytarabine (i.e., the effect of the YKO/KD on cell proliferation independent of drug treatment; see methods). (B–F) Within each panel, clusters in the respective categories are displayed, left to right, in descending order, by relative strength of drug–gene interaction effects evident by the heatmaps. (B) Enhancing gene–drug interactions for both drugs. (C) Gemcitabine-specific enhancement. (D) Cytarabine-specific enhancement. (E) Suppressing gene–drug interactions for both drugs. (F) Gemcitabine-specific suppression. (G) The algorithm for integrating yeast phenomic and cancer pharmacogenomics data: For all cell lines from the Genomics of Drug Sensitivity in Cancer (GDSC) database (either lung, hematopoietic and lymphoid, or across all tissues) with increased drug sensitivity, underexpressed (UES) genes were highlighted by yeast homologs that were deletion enhancing, while overexpressed (OES) genes were highlighted by yeast homologs that were deletion suppressing. (H) Yeast–human homologs identified as described in G. The category of homology from BiomaRt is indicated in the left column of each heatmap (see homology color key) shown at right. The gene label color (at far right) indicates whether the human homolog was found in PharmacoDB for both drugs (black), cytarabine (teal), or gemcitabine (gold). [125]Additional Files 5 (File B) and 8 (Files B–D) contain all REMc heatmaps of the types indicated to the left and right, respectively, in panel H. [126]Additional File 8 includes information for all yeast–human homologs from each category suggested by the study to exhibit functionally conserved gene–drug interaction. 2.5. Gene Ontology Term Finder (GTF) A python script was used to format REMc clusters for analysis with the command line version of the GO Term Finder (GTF) tool downloaded from [127]http://search.cpan.org/dist/GO-TermFinder/ [[128]41]. GTF reports on enrichment of gene ontology (GO) terms by comparing the ratio of genes assigned to a term within a cluster to the respective ratio involving all genes tested. [129]Additional File 5, File C contains GTF analysis of all REMc clusters. GO-enriched terms from REMc were investigated with respect to genes representing the term and literature underlying their annotations [[130]42]. 2.6. Gene Ontology Term Averaging (GTA) Analysis In addition to using GTF to survey functional enrichment in REMc clusters, we developed GTA as a complementary workflow, using the GO information on SGD at [131]https://downloads.yeastgenome.org/curation/literature/ to perform the following analysis: 1. Calculate the average and SD for interaction values of all genes in a GO term. 2. Filter results to obtain terms having GTA value greater than 2 or less than −2. 3. Obtain GTA scores defined as |GTA value| - gtaSD; filter for GTA score > 2. The GTA analysis is contained in [132]Additional File 6 as tables and interactive plots created using the R plotly package [133]https://CRAN.R-project.org/package=plotly. GTA results were analyzed using both the L and K interaction scores and are included in [134]Additional File 6 (Files A–C). 2.7. Prediction of Human Homologs that Influence Tumor Response to Gemcitabine or Cytarabine PharmacoDB holds pharmacogenomics data from cancer cell lines, including transcriptomics and drug sensitivity [[135]32]. The PharmacoGx R/Bioconductor package [[136]43] was used to analyze the GDSC1000 ([137]https://pharmacodb.pmgenomics.ca/datasets/5) and gCSI ([138]https://pharmacodb.pmgenomics.ca/datasets/4) datasets, which contained transcriptomic and drug sensitivity results. A p-value < 0.05 was used for differential gene expression and drug sensitivity. For gene expression, the sign of the standardized coefficient denotes increased (+) or decreased (−) expression. The biomaRt R package [[139]44,[140]45] was used with the Ensembl database [[141]46] to match yeast and human homologs from the phenomic and transcriptomic data, classifying yeast–human homology as one to one, one to many, and many to many. The Princeton Protein Orthology Database (PPOD) was also used to manually review homology in further consideration of whether or not to include particular genes in the results sections and whether or not cancer literature searches related to the homologs would be worthwhile [[142]47]. 3. Results 3.1. Quantitative Phenomic Characterization of Differential Gene–Drug Interaction The Q-HTCP workflow incorporates high-throughput kinetic imaging and analysis of proliferating 384-culture cell arrays plated on agar media to obtain CPPs for measuring gene–drug interaction, as previously described [[143]7,[144]9,[145]37]. To apply it for analysis of dCK substrates, a tetracycline-inducible human dCK allele was introduced into the complete YKO/KD library by the synthetic genetic array method [[146]29,[147]48] ([148]Figure 1B). The dependence of gemcitabine and cytarabine toxicity on dCK expression was demonstrated for the reference strain ([149]Figure 2A–F), as the two nucleosides exerted cytotoxicity only if dCK was induced by the addition of doxycycline. Induction of dCK had no effect on proliferation in the absence of gemcitabine or cytarabine ([150]Figure 2C–F). Interaction scores were calculated by departure of the observed CPP for each YKO/KD strain from that expected based on the observed phenotypes for the reference strain treated and untreated with drugs and the YKO/KD strain in the absence of drugs, incorporating multiple drug concentrations, 768 replicate reference strain control cultures, and summarized by linear regression as z-scores [[151]6,[152]7,[153]8,[154]9,[155]33,[156]37]. Gene interaction scores with absolute value greater than two were selected for global analysis and termed deletion enhancers (z-score_L ≥ 2 or z-score_K ≤ −2) or deletion suppressors (z-score_L ≤ −2 or z-score_K ≥ 2) of drug cytotoxicity, revealing functions that buffer or promote drug cytotoxicity, respectively [[157]39] ([158]Figure 2). Growth inhibition was greater for gemcitabine than for cytarabine ([159]Figure 2A–F), however, the phenotypic variance was also less for cytarabine, such that interactions of smaller effect size were detectable and the range of scores was greater ([160]Additional File 1, Table S6). The CPP, ‘L’, (the time at which half carrying capacity is reached), is most sensitive to growth inhibitory perturbation, while ‘K’ (carrying capacity) reports on more extreme growth differences ([161]Figure 2A–H). Low correlation between the gene–drug interaction profiles suggested differential buffering of these two drugs, consistent with their distinct antitumor efficacies ([162]Figure 2I,J). 3.2. Functional Analysis of Gene Interaction Modules Recursive expectation-maximization clustering (REMc) was used to identify modules of genes that shared similar profiles of buffering or promoting nucleoside toxicity of gemcitabine or cytarabine [[163]40] (see [164]Figure 3A–F; [165]Table 1; [166]Additional File 5). As described previously, REMc results were assessed with GO Term Finder for gene ontology functional enrichment [[167]41] and heatmaps generated by first adding data regarding the main effect of the gene knockout or knockdown (i.e., no drug) on cell proliferation, termed ‘shift’ (see methods), followed by hierarchical clustering [[168]40,[169]41]. GO term average (GTA) scores, which are based on the average and standard deviation of drug–gene interaction for all genes of each GO term [[170]39], were used as a complement to REMc/GTF for identifying functions that buffer or promote drug effects ([171]Table 2, [172]Figure 4, and [173]Additional File 6, Files A–C). Yeast–human homologs were judged, regarding causality of differential gene expression associated with sensitivity to gemcitabine or cytarabine, by the correspondence of yeast phenomic and cancer pharmacogenomics results, thus establishing a model resource to test the utility of yeast phenomics to inform cancer genetic profiling for predicting drug-specific, antitumor efficacy ([174]Figure 3G–H). Table 1. GO terms enriched in REMc clusters. GO Term Drug INT O Cluster Genes in Term p-Value Genes Fig. GTA Gem L GTA Cyt L Ubp3-Bre5 deubiquitination complex Both Enh C 2-0.2-0 2/2 2.57 × 10^−5 UBP3:BRE5 [175]Figure 5D 19.8 14.32 positive regulation of DNA-dependent DNA replication initiation Both Enh P 1-0-2 3/4 2.09 × 10^−4 RFM1:FKH2:SUM1 [176]Figure 5B 15.7 4.9 Mre11 complex Both Enh C 2-0.14-1 2/3 5.66 × 10^−4 RAD50:XRS2 [177]Figure 5B 13.7 26.6 HOPS complex Both Enh C 2-0.14-1 2/7 3.94 × 10^−3 PEP3:VPS33 [178]Figure 5D 12.0 4.8 CORVET complex Both Enh C 2-0.14-1 2/7 3.94 × 10^−3 PEP3:VPS33 [179]Figure 5D 10.4 4.3 RecQ helicase-Topo III complex Both Enh C 1-0-14 2/3 3.31 × 10^−3 SGS1:RMI1 [180]Figure 5B 7.5 14.6 GET complex Both Enh C 2-0.14-0 2/3 4.68 × 10^−4 GET1:GET2 [181]Figure 5D 3.3 18.6 DNA integrity checkpoint Both Enh P 1-0-14 4/40 3.85 × 10^−3 DUN1:RAD17:RAD24:SGS1 [182]Figure 5A 4.8 4.8 α-glucoside transmembrane transporter activity Cyt Enh F 2-0.17-3 2/2 5.98 × 10^−3 MAL31:MAL11 Figure 7A −0.7 2.2 intralumenal vesicle formation Gem Enh P 1-0-10 3/7 2.90 × 10^−3 DOA4:VPS24:BRO1 [183]Figure 6A 9.0 1.6 HDA1 complex Gem Enh C 1-0-0 2/3 7.08 × 10^−2 HDA1:HDA3 [184]Figure 6B 4.8 0.3 Swr1 complex Gem Enh C 1-0-11 3/12 3.46 × 10^−2 SWC3:VPS71:SWR1 [185]Figure 6B 2.9 −1.6 peptidyl-tyrosine dephosphorylation Gem Enh P 1-0-0 5/20 2.18 × 10^−3 OCA2:SIW14:OCA1:OCA4:OCA6 [186]Figure 6C 1.5 0.5 Set1C/COMPASS complex Gem Enh C 1-0-0 3/6 5.74 × 10^−3 SDC1:SWD3:BRE2 [187]Figure 6B 1.0 0.6 phospholipid-translocating ATPase activity Gem Sup F 1-0-8 3/7 9.70 × 10^−3 DRS2:LEM3:DNF2 [188]Figure 6D −1.6 −0.9 [189]Open in a new tab For each GO term, the table indicates which drugs interact with it, the interaction type (enhancing or suppressing), the ontology (‘O’) it derives from (cellular process or component, or molecular function), the REMc cluster ID from which the term was most specific, the fraction of the genes in the term that were observed in the cluster, and the p-value for enrichment of the genes. Relevant figures and associated GTA data are also given. Table 2. GO terms identified by gene ontology term averaging (GTA). Term Drug INT_Type Ont Cluster p-Value Genes Fig. Gem GTA_K Gem GTA_L Cyt GTA_K Cyt GTA_L checkpoint clamp complex Both Enh L/K C NA NA RAD17 | MEC3 [190]Figure 5B −7.3 13.8 −23.5 15.4 HOPS complex Both Enh L/K C 2-0.14-1 3.94 × 10^−3 VPS16 | VPS8 | PEP3 | VPS41 | VPS33 | PEP5 [191]Figure 5D −6.3 12.0 −11.4 4.8 Mre11 complex Both Enh L/K C 2-0.14-1 5.66 × 10^−4 MRE11 | RAD50 | XRS2 [192]Figure 5B −8.8 13.7 −39.3 26.6 RecQ helicase-Topo III complex Both Enh L/K C 1-0-14 3.31 × 10^−3 RMI1 | SGS1 | TOP3 [193]Figure 5B −7.7 7.5 −24.7 14.6 Ubp3-Bre5 deubiquitination complex Both Enh L/K C 2-0.2-0 2.57 × 10^−5 UBP3 | BRE5 [194]Figure 5D −9.2 19.8 −16.9 14.3 vesicle fusion with vacuole Both Enh L/K P NA NA VAM3 | VPS33 [195]Figure 5D −7.4 13.3 −11.4 7.1 Sec61 translocon complex Cyt Enh K C NA NA SEC61 | SBH2 [196]Figure 7A −0.4 1.1 −5.1 1.9 HIR complex Cyt Enh L C NA NA HIR1 | HIR2 | HPC2 | HIR3 [197]Figure 7A −1.0 1.0 −0.6 2.5 sphinganine kinase activity Cyt Enh L F NA NA LCB4 | LCB5 [198]Figure 7A −0.1 0.3 −1.2 3.9 protein localization to septin ring Cyt Enh L/K P NA NA ELM1 | HSL1 [199]Figure 7A −1.3 2.5 −17.8 21.9 autophagosome maturation Gem Enh K P NA NA VAM3 | CCZ1 [200]Figure 6A −5.6 7.7 −1.6 2.5 Elongator holoenzyme complex Gem Enh K C NA NA TUP1 | ELP4 | ELP2 | IKI3 | IKI1 | ELP3 | ELP6 [201]Figure S4C −3.6 3.4 −2.6 2.5 ESCRT I complex Gem Enh K C NA NA STP22 | VPS28 | SRN2 | MVB12 [202]Figure 5D −6.9 9.1 −0.8 2.5 negative regulation of macroautophagy Gem Enh K P NA NA PHO85 | PHO80 | KSP1 | PCL5 | SIC1 [203]Figure 6A −5.8 9.4 −4.1 1.8 protein urmylation Gem Enh K P NA NA ELP2 | UBA4 | NCS2 | URM1 | URE2 | ELP6 [204]Figure S4C −3.7 1.5 1.0 1.2 CORVET complex Gem Enh L/K C 2-0.14-1 3.94 × 10^−3 VPS16 | VPS8 | PEP3 | VPS41 | VPS33 | VPS3 | PEP5 [205]Figure 5D −6.6 10.4 −10.4 4.3 ESCRT-0 complex Gem Enh L/K C NA NA VPS27 | HSE1 [206]Figure 5D −5.7 10.4 −3.9 2.6 HDA1 complex Gem Enh L/K C 1-0-0 7.08 × 10^−2 HDA3 | HDA1 | HDA2 [207]Figure 6B −4.8 4.8 −0.6 0.3 GATOR (Iml1) complex Gem Enh L/K C NA NA NPR2 | NPR3 [208]Figure 6A −4.4 6.4 1.0 2.2 intralumenal vesicle formation Gem Enh L/K P 1-0-10 2.90 × 10^−3 VPS20 | VPS24 | BRO1 | DOA4 | VPS4 | SNF7 [209]Figure 6A −5.7 9.0 −1.8 1.6 positive regulation of DNA-dependent DNA replication initiation Gem Enh L/K P 1-0-2 2.09 × 10^−4 SUM1 | FKH2 | RFM1 | FKH1 [210]Figure 5B −8.1 15.7 −2.4 4.9 RAVE complex Gem Enh L/K C NA NA RAV1 | RAV2 [211]Figure 6A −4.2 3.5 0.6 −0.2 GARP complex Gem Sup L C NA NA VPS51 | VPS53 | VPS54 | VPS52 [212]Figure 6D 1.7 −3.4 1.5 −1.0 Lem3p-Dnf1p complex Gem Sup L C NA NA DNF1 | LEM3 [213]Figure 6D 1.6 −3.4 −0.1 0.2 phosphatidylserine biosynthetic process Gem Sup L C NA NA DEP1 | CHO1 | UME6 [214]Figure 6D 2.6 −3.7 0.8 −0.3 [215]Open in a new tab See [216]Table 1 for data descriptions. ‘NA’ indicates terms identified by GTA only (i.e., not identified by REMc/GTF). Figure 4. [217]Figure 4 [218]Open in a new tab GO annotations associated with deletion enhancement or suppression of gemcitabine and/or cytarabine cytotoxicity. Representative GO terms are listed, which were identified by REMc/GTF (orange), GTA (purple), or both methods, for enhancement (above dashed line) or suppression (below dashed line) of gemcitabine (left, red), cytarabine (right, blue), or both media types (black). Term-specific heatmaps were manually reviewed to decide which terms should be included. Distance above or below the horizontal dashed line reflects the average interaction score for genes identified by REMc/GTF or the GTA score (see methods). See [219]Additional Files 5 and 6 for all REMc/GTF and GTA results, and [220]Additional File 7 for GO term-specific heatmaps. Heatmaps were also produced systematically to visualize drug–gene interaction profiles for all genes assigned to GO terms identified by REMc/GTF or GTA; these are referred to as term-specific heatmaps and are grouped by GO term parent–child relationships ([221]Additional File 7). Cancer pharmacogenomics data in PharmacoDB were mined using PharmacoGx [[222]43] and biomaRt [[223]44,[224]45] with the Genomics of Drug Sensitivity in Cancer (GDSC) [[225]30,[226]49] or Genentech Cell Line Screening Initiative (gCSI) [[227]31,[228]50] datasets to match yeast drug–gene interaction by homology to differential gene expression in gemcitabine or cytarabine sensitive cancer cell lines ([229]Figure 3G–H; [230]Additional File 8). Yeast gene deletion enhancers identified human homologs underexpressed in gemcitabine- or cytarabine-sensitive cells, termed UES, while yeast gene deletion suppressors identified human homologs overexpressed in drug-sensitive cells, termed OES ([231]Figure 3G). The analysis was focused on the GDSC database, because it had expression data available for both gemcitabine and cytarabine; however, analysis of the gCSI data was also conducted for gemcitabine ([232]Additional File 8, File A). Differential expression was analyzed: (1) across all tissue types, to consider interactions that might be applicable in novel treatment settings; (2) in hematopoietic and lymphoid tissue (HaL); and (3) in lung tissue, as cytarabine and gemcitabine are used to treat HaL and lung cancers, respectively. Gemcitabine is also used for pancreatic cancer; however, the number of cell lines tested (30) was lower than for lung (156) or HaL (152). Thus, yeast genes that were deletion enhancing or suppressing were catalogued with human homologs that were UES or OES in PharmacoDB ([233]Figure 3G,H, [234]Table 3, [235]Table 4 and [236]Table 5, and [237]Additional File 8, File A). Table 3. Yeast–human homologs predicted to similarly buffer or promote both gemcitabine and cytarabine toxicity. yGene hGene H Drug Cluster Tissue Gem K Cyt K Gem L Cyt L Ref Description (Human) NAM7 HELZ 2 Cyt 1-0-14 L −6.5 −16.7 1.1 13.6 [[238]66,[239]67,[240]68,[241]69] helicase with zinc finger NAM7 HELZ2 2 Cyt 1-0-14 A, H −6.5 −16.7 1.1 13.6 helicase with zinc finger 2 NAM7 UPF1 2 Cyt 1-0-14 L −6.5 −16.7 1.1 13.6 [[242]70,[243]71,[244]72] UPF1, RNA helicase and ATPase PTC1 PPM1E 2 Both 1-0-14 L −8.8 −12.7 7.9 15.7 [[245]73] protein phosphatase, Mg2+/Mn2+ dependent 1E PTC1 PPM1L 2 Both 1-0-14 A, H −8.8 −12.7 7.9 15.7 [[246]74] protein phosphatase, Mg2+/Mn2+ dependent 1L RAD24 RAD17 1 Gem 1-0-14 H, L −7.4 −27.6 14.2 8.3 [[247]75,[248]76,[249]77,[250]78,[251]79] RAD17 checkpoint clamp loader component SGS1 RECQL5 2 Cyt 1-0-14 L −8.4 −33.4 3.4 19.3 [[252]61,[253]62] RecQ like helicase 5 KTI11_2 DPH3 1 Cyt 1-0-14 H −7.7 −10.3 6.5 9.1 [[254]80,[255]81,[256]82] diphthamide biosynthesis 3 BIM1_2 MAPRE2 2 Gem 1-0-14 A −7.7 −15.4 16.0 20.0 [[257]83] microtubule associated protein RP/EB family member 2 BIM1_2 MAPRE2 2 Both 1-0-14 L −7.7 −15.4 16.0 20.0 [[258]83] microtubule associated protein RP/EB family member 2 BIM1_2 MAPRE3 2 Gem 1-0-14 A −7.7 −15.4 16.0 20.0 [[259]84] microtubule associated protein RP/EB family member 3 ASF1 ASF1B 2 Cyt 2-0.16-1 L −6.1 −9.5 4.1 8.3 [[260]85] anti-silencing function 1B histone chaperone AVL9 AVL9 1 Cyt 2-0.16-1 H −4.3 −2.5 0.2 2.9 [[261]86,[262]87,[263]88] AVL9 cell migration associated PMR1 ATP1A1 2 Cyt 2-0.16-1 A, H −3.8 −9.8 3.6 10.1 [[264]89] ATPase Na+/K+ transporting subunit α 1 PMR1 ATP1A2 2 Gem 2-0.16-1 A −3.8 −9.8 3.6 10.1 [[265]90] ATPase Na+/K+ transporting subunit α 2 PMR1 ATP1A3 2 Cyt 2-0.16-1 L −3.8 −9.8 3.6 10.1 ATPase Na+/K+ transporting subunit α 3 PMR1 ATP1A4 2 Gem 2-0.16-1 A −3.8 −9.8 3.6 10.1 ATPase Na+/K+ transporting subunit α 4 PMR1 ATP1A4 2 Both 2-0.16-1 H −3.8 −9.8 3.6 10.1 ATPase Na+/K+ transporting subunit α 4 PMR1 ATP2C1 2 Cyt 2-0.16-1 A −3.8 −9.8 3.6 10.1 [[266]91,[267]92] ATPase secretory pathway Ca2+ transporting 1 PMR1 ATP2C1 2 Both 2-0.16-1 H −3.8 −9.8 3.6 10.1 [[268]91,[269]92] ATPase secretory pathway Ca2+ transporting 1 TOP3 TOP3A 2 Cyt 2-0.16-1 L −5.2 −4.0 3.3 3.4 [[270]63,[271]64,[272]65] DNA topoisomerase III α VPS21 RAB21 3 Cyt 2-0.16-1 A, H −7.2 −4.1 −0.4 2.4 [[273]93,[274]94] RAB21, member RAS oncogene family VPS21 RAB22A 3 Gem 2-0.16-1 A −7.2 −4.1 −0.4 2.4 [[275]95,[276]96,[277]97] RAB22A, member RAS oncogene family ACB1_2 ACBD4 2 Gem 2-0.16-1 H −5.4 −4.8 4.5 0.6 [[278]98,[279]99] acyl-CoA binding domain containing 4 ACB1_2 ACBD5 2 Cyt 2-0.16-1 H −5.4 −4.8 4.5 0.6 [[280]100] acyl-CoA binding domain containing 5 ACB1_2 DBI 2 Cyt 2-0.16-1 A, H −5.4 −4.8 4.5 0.6 [[281]101,[282]102,[283]103] diazepam binding inhibitor, acyl-CoA binding protein CPR3 PPIA 3 Cyt 2-0.8-1 A, H 2.1 1.6 −4.1 −2.8 [[284]104,[285]105,[286]106] peptidylprolyl isomerase A CPR3 RGPD4 3 Gem 2-0.8-1 A 2.1 1.6 −4.1 −2.8 RANBP2-like and GRIP domain containing 4 ELO3 ELOVL1 3 Both 2-0.8-1 L 2.2 1.3 −3.4 −4.0 [[287]107,[288]108] ELOVL fatty acid elongase 1 ELO3 ELOVL2 3 Cyt 2-0.8-1 H 2.2 1.3 −3.4 −4.0 [[289]109] ELOVL fatty acid elongase 2 ELO3 ELOVL4 3 Cyt 2-0.8-1 H 2.2 1.3 −3.4 −4.0 ELOVL fatty acid elongase 4 ELO3 ELOVL5 3 Cyt 2-0.8-1 H 2.2 1.3 −3.4 −4.0 ELOVL fatty acid elongase 5 ELO3 ELOVL6 3 Both 2-0.8-1 A, L 2.2 1.3 −3.4 −4.0 [[290]110,[291]111] ELOVL fatty acid elongase 6 MDL2 ABCB10 3 Gem 2-0.8-1 H 2.5 1.5 −3.0 −3.0 [[292]112] ATP binding cassette subfamily B member 10 MDL2 TAP1 3 Cyt 2-0.8-1 L 2.5 1.5 −3.0 −3.0 transporter 1, ATP binding cassette subfamily B member PIF1 PIF1 2 Gem 2-0.8-1 A 2.2 1.5 −4.5 −3.4 [[293]113] PIF1 5’-to-3’ DNA helicase RPS1B RPS3A 1 Both 2-0.8-1 A 2.3 0.9 −3.9 −2.3 [[294]114,[295]115] ribosomal protein S3A SAC3 MCM3AP 2 Gem 2-0.8-1 H 2.2 1.5 −5.2 −3.8 [[296]116] minichromosome maintenance complex component 3 associated protein SAC3 SAC3D1 2 Cyt 2-0.8-1 H 2.2 1.5 −5.2 −3.8 [[297]117,[298]118] SAC3 domain containing 1 YTA7 ATAD2 2 Both 2-0.8-1 A, H 1.8 1.0 −6.0 −3.6 [[299]119,[300]120,[301]121,[302]122,[303]123,[304]124,[305]125] ATPase family, AAA domain containing 2 YTA7 ATAD2B 2 Both 2-0.8-1 H 1.8 1.0 −6.0 −3.6 ATPase family, AAA domain containing 2B [306]Open in a new tab Genes selected for discussion in the results were included in the table. The homology types (H) are one to one (1), one to many (2), and many to many (3). Drugs (Gem, Cyt, or Both) with which the genes interacted in a UES or OES manner in the GDSC database are indicated. The REMc clusters 1-0-14 and 2-0.16-1 are deletion enhancing and 2-0.8-1 is deletion suppressing (see [307]Figure 5C,D). Tissue types from which genes were UES or OES in the PharmacoDB data are indicated for across all tissue (A), lung (L), and hematopoietic (H). Related references cited (Ref), and gene descriptions are given.