A Humanized Yeast Phenomic Model of Deoxycytidine Kinase to Predict Genetic Buffering of Nucleoside Analog Cytotoxicity

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.


Introduction
Genomics has enabled targeted therapy aimed at cancer driver genes and oncogenic addiction [1], yet traditional cytotoxic chemotherapeutic agents remain among the most widely used and efficacious anticancer therapies [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 [3,4]. A potential advantage of so-called synthetic lethality-based treatment strategies is that they could have efficacy against passenger gene mutation or compensatory 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.
The synthetic genetic array (SGA) method, a way to introduce an allele of interest into the YKO/KD library and recover haploid double mutants [28,29], was used to derive a haploid YKO/KD collection with doxycycline-inducible dCK expression.

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 [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 [9,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) [7]. The CPPs, primarily K and L, were used as phenotypes to measure drug-gene interaction.

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 [6,9,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 [5,6,9,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 (Figure 2C-F; Additional File 2, Figure S2A-D). Thus, shift was calculated in the presence of 5 µg/mL doxycycline without nucleoside analog.
Interaction scores were calculated as previously described [9,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 ): 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 INTmean(KD INT ))/SD(KD INT ).

Recursive Expectation-Maximization Clustering (REMc) and Heatmap Generation
REMc is a probability-based clustering method and was performed as previously described [40]. Clusters obtained by Weka 3.5, an EM-optimized Gaussian mixture-clustering module, were subjected to hierarchical clustering in R (http://www.r-project.org/) to further aid visualization with heatmaps. REMc was performed using L and K interaction z-scores ( 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. 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.

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 http://search.cpan.org/dist/GO-TermFinder/ [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. 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 [42].

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 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.
The GTA analysis is contained in Additional File 6 as tables and interactive plots created using the R plotly package https://CRAN.R-project.org/package=plotly. GTA results were analyzed using both the L and K interaction scores and are included in Additional File 6 (Files A-C).

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 [32]. The PharmacoGx R/Bioconductor package [43] was used to analyze the GDSC1000 (https://pharmacodb.pmgenomics.ca/datasets/5) and gCSI (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 [44,45] was used with the Ensembl database [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 [47].

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 [7,9,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 [29,48] ( Figure 1B). The dependence of gemcitabine and cytarabine toxicity on dCK expression was demonstrated for the reference strain (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 ( Figure 2C-F). (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 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 [6][7][8][9]33,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 [39] (Figure 2).
Growth inhibition was greater for gemcitabine than for cytarabine (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 (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 (Figure 2A-H). Low correlation between the gene-drug interaction profiles suggested differential buffering of these two drugs, consistent with their distinct antitumor efficacies ( Figure 2I,J).

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 [40] (see Figure 3A-F; Table 1; Additional File 5). As described previously, REMc results were assessed with GO Term Finder for gene ontology functional enrichment [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 [40,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 [39], were used as a complement to REMc/GTF for identifying functions that buffer or promote drug effects (Table 2, Figure 4, and 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 ( Figure 3G-H).  . 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 Additional Files 5 and 6 for all REMc/GTF and GTA results, and Additional File 7 for GO term-specific heatmaps.   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.  Elongator holoenzyme complex Gem Table 1 for data descriptions. 'NA' indicates terms identified by GTA only (i.e., not identified by REMc/GTF).
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 (Additional File 7).
Cancer pharmacogenomics data in PharmacoDB were mined using PharmacoGx [43] and biomaRt [44,45] with the Genomics of Drug Sensitivity in Cancer (GDSC) [30,49] or Genentech Cell Line Screening Initiative (gCSI) [31,50] datasets to match yeast drug-gene interaction by homology to differential gene expression in gemcitabine or cytarabine sensitive cancer cell lines ( Figure 3G-H; 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 ( 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 (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 ( Figure 3G,H, Tables 3-5, and Additional File 8, File A).
In summary REMc, GTF, and GTA revealed functional genetic modules that alternatively buffer (deletion enhancing) or promote (deletion suppressing) drug cytotoxicity [5,40,51] and illustrated whether the effects were shared or differential between gemcitabine and cytarabine ( Figure 4). Yeast phenomic information was integrated with pharmacogenomics data results according to yeast-human gene homology to identify correlated differential gene expression associated with drug sensitivity in cancer cell lines ( Figures 5-7). This approach serves to generate hypotheses regarding whether differential expression of a particular gene is causal for increased drug sensitivity [52] and ultimately whether yeast phenomic models can improve the predictive value of cancer pharmacogenomics data in the context of precision oncology [53][54][55][56][57][58].

Genetic Modules that Buffer Cytotoxicity of Both Gemcitabine and Cytarabine
To characterize gemcitabine and cytarabine, which have similar molecular structures and mechanisms of action, yet different spectra of antitumor efficacy, we first surveyed for buffering genes shared in common. Examples of genes with deletion-enhancing interactions for both drugs are displayed in clusters 2-0.2-0, 1-0-14, 2-0.2-2, and 2-0.16-1 ( Figure 3B). GO enrichment was observed in these clusters for the DNA integrity checkpoint; positive regulation of DNA replication; and the Mre11, RecQ helicase-Topo III, CORVET, HOPS, GET, and Ubp3-Bre5 deubiquitination complexes ( Figure 4, Table 1). GTA identified many of the same functions and additionally revealed the terms vesicle fusion with vacuole and checkpoint clamp complex ( Table 2). We mapped yeast gene-drug interactions to respective human homologs in PharmacoDB to find evidence for evolutionary conservation of gene-drug interaction ( Figure 5C,D, Additional File 8, Files B-D) and buffering mechanisms.

DNA Integrity Checkpoint and Repair-Related Complexes
As gemcitabine and cytarabine triphosphate analogs are incorporated into DNA, we anticipated shared interactions with genes functioning in DNA metabolism and repair. Overlap was observed, however there were differential effects between genes assigned to the same gene ontology terms, such that GO TermFinder enrichment in REMc clusters was less than might have been expected. For example, deletion-enhancing gene-drug interaction for the GO term, DNA integrity checkpoint, was enriched in cluster 1-0-14, which displayed deletion enhancement for both gemcitabine and cytarabine (Table 1,  Figures 3 and 5A). However, its child term, intra-S DNA damage checkpoint, was not GO-enriched because of differential clustering among drug-gene interactions associated with the term (Additional File 5, File C). Similarly, intra-S DNA damage checkpoint was not identified by GTA due to variation in interaction between genes assigned to the term, highlighting the utility in displaying the phenomic data for each GO term for manual review ( Figure 5A).

Figure 5.
Drug-gene interaction common to gemcitabine and cytarabine. Genes that similarly influence the cytotoxicity of both gemcitabine and cytarabine suggest common pathways that buffer or promote toxicity, as illustrated by: (A) GO term-specific heatmaps for DNA integrity checkpoint and its child term intra-S DNA damage checkpoint, which buffer gemcitabine and cytarabine, along with (B) genes comprising other DNA checkpoint/repair-related GO terms, such as positive regulation of DNA-dependent DNA replication initiation, and the Mre11, checkpoint clamp and RecQ helicase-Topo III complexes; (C,D) REMc clusters filtered for PharmacoDB results for yeast-human homologs that exhibited (C) deletion enhancement and UES or (D) deletion suppression and OES; and (E) deletion enhancing endosomal transport-related GO terms, including vesicle fusion with vacuole, and the CORVET/HOPS, ESCRT, GET, and Ubp3-Bre5 deubiquitination complexes. Gene labels are color-coordinated with legends in panels B and E, and as described in Figure 3H for panels C and D. Genes in the YKD collection are underlined.
Enriched complexes functionally related to the DNA integrity checkpoint function included the RecQ helicase-Topo III, the checkpoint clamp, and the Mre11 complexes ( Figure 5B). Rmi1, Top3, and Sgs1 form the RecQ helicase Topo III complex, which is involved in Rad53 checkpoint activation and maintenance of genome integrity [59], and together with replication protein A function in DNA decatenation and disentangling of chromosomes [60]. RMI1 and SGS1 deletion-enhancement clustered together in 1-0-14, while TOP3 had a similar, but slightly weaker interaction pattern in cluster 1-0-16 (Additional File 5, File B). The human homolog of SGS1, RECQL5, was UES for cytarabine in lung cancer cells ( Figure 5C; see 1-0-14 in 5C, all cluster heatmaps available in Additional File 5, File B). RECQL5 preserves genome stability during transcription elongation, and deletion of RECQL5 increases cancer susceptibility [61,62]. Human TOP3A was also UES for cytarabine in lung tissue ( Figure 5C; 2-0.16-1). TOP3A is underexpressed in ovarian cancer, and mutations in TOP3A are associated with increased risk for acute myeloid leukemia and myelodysplastic syndromes, suggesting potential cancer vulnerabilities if somatic, but they can also occur in the germline, which would lead to enhanced host toxicity [63][64][65]. . Gemcitabine-specific gene interaction. (A-C) Cellular processes that buffer gemcitabine to a greater extent than cytarabine included: (A) autophagy-related processes; (B) histone modification and chromatin remodeling (particularly for K interaction); and (C) peptidyl-tyrosine dephosphorylation, representing the genes OCA (1-6) (OCA5 was manually added to the panel (see text); OCA3/SIW14 are aliases). (D,E) When comparing gene-drug interactions of homologs across cancer pharmacogenomic and yeast phenomic experiments, human genes are predicted to (D) buffer gemcitabine toxicity if they are UES and deletion enhancing, or to (E) promote gemcitabine toxicity if they are OES and deletion suppressing. (F) Apoptosis-related genes and complexes were observed to promote toxicity of gemcitabine more than toxicity of cytarabine. Gene labels are color-coordinated with legends in panels A, B, and F, and as described in Figure 3H for panels D and E. 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 (  The checkpoint clamp in yeast is comprised of Rad17/hRad1, Ddc1, and Mec3, which function downstream of Rad24/hRad17 in the DNA damage checkpoint pathway [75][76][77] to recruit yDpb11/hTopB1 to stalled replication forks and activate the yMec1/hATR protein kinase activity, initiating the DNA damage response [78]. The human homolog of yeast RAD24, RAD17, was UES for gemcitabine in both lung and hematopoietic and lymphoid tissue ( Figure 5C; 1-0-14), representing a synthetic lethal relationship of potential therapeutic relevance. Consistent with this finding in yeast, depletion of hRAD17 can sensitize pancreatic cancer cells to gemcitabine [79].
Mre11, Xrs2, and Rad50 constitute the Mre11 complex, which participates in the formation and processing of double-strand DNA breaks involved in recombination and repair [126], and clustered together in 1-0-14 ( Figure 5B,C). Deficiency in the Mre1 complex is known to sensitize human cells to nucleoside analog toxicity [127], as also seen in cancer cell lines deficient for other checkpoint-signaling genes, such as Rad9, Chk1, or ATR, [128]. Single nucleotide polymorphisms in DNA damage response (ATM and CHEK1) have been associated with overall survival in pancreatic cancer patients treated with gemcitabine and radiation therapy [129]. Taken together, the results highlight evolutionarily conserved genes that function in DNA replication and recombination-based repair and are required to buffer the cytotoxic effects of both cytarabine and gemcitabine.

Positive Regulation of DNA-Dependent DNA Replication Initiation
The term, positive regulation of DNA-dependent DNA replication initiation, was identified by REMc/GTF and GTA for buffering interactions with both drugs, though stronger for gemcitabine (Tables 1 and 2). Genes representing this term were FKH2, RFM1, and SUM1 ( Figure 5B). The origin binding protein, Sum1, is required for efficient replication initiation [130] and forms a complex with Rfm1 and the histone deacetylase, Hst1, which is recruited to replication origins to deacetylate H4K5 for initiation [131]. HST1 was also a strong deletion enhancer but was observed only for the L parameter and clustered in 2-0.2-2. The forkhead box proteins, Fkh1 and Fkh2, contribute to proper replication origin timing and long range clustering of origins in G1 phase [132] and appear to buffer the cytotoxicity of gemcitabine more so than cytarabine, with FKH2 deletion showing a stronger effect than its paralog ( Figure 5B). Multiple human forkhead box protein homologs (yFKH2/hFOXJ1/FOXG1/FOXJ3/FOXH1) ( Figure 6D) were observed as UES in PharmacoDB, of which FOXJ1 underexpression is a marker of poor prognosis in gastric cancer [133], reduced expression of FOXG1 is correlated with worse prognosis in breast cancer [134], FOXJ3 is inhibited by miR-517a and associated with lung and colorectal cancer cell proliferation and invasion [135,136], and FOXH1 is overexpressed in breast cancer and FOXH1 inhibition reduces proliferation in breast cancer cell lines [137]. Although not UES in PharmacoDB, inhibition of the HST1 homolog, SIRT1, by Tenovin-6 inhibits the growth of acute lymphoblastic leukemia cells and enhances cytarabine cytotoxicity [138], enhances gemcitabine efficacy in pancreatic cancer cell lines, and improves survival in a pancreatic cancer mouse model [139]. Thus, loss of this gene module that positively regulates DNA replication initiation appears to be robustly involved in oncogenesis and is also synthetic lethal with gemcitabine and cytarabine.

'Non-GO-enriched' Homolog Pairs with Corresponding UES and Deletion Enhancement
We next explored yeast-human homologs exhibiting yeast deletion enhancement and underexpression sensitivity in cancer, systematically and regardless of whether their functions were enriched within gene ontology (Table 3, Figure 5C). 'Non-enriched' interaction can be explained by a small total number of genes performing the function, only select genes annotated to a term impacting the phenotype, by genes contributing to a function without yet being annotated to it, by novel functions, and other possibilities.
With regard to the above, human homologs of the yeast type 2C protein phosphatase, PTC1, included PPM1L and PPM1E ( Figure 5C; 1-0-14). PPM1L has reduced expression in familial adenomatous polyposis [74], while PPM1E upregulation has been associated with cell proliferation in gastric cancer [73]. Such differential interactions of paralogs could result from tissue specific expression and functional differentiation of regulatory proteins. Previously, we reported ptc1-∆0 to buffer transcriptional repression of RNR1 [33], which is upregulated as part of the DNA damage response to increase dNTP pools [153].
The microtubule binding proteins, yBIM1/hMAPRE2/hMAPRE3, were deletion enhancing in yeast and UES in cancer for gemcitabine ( Figure 5C; 1-0-14), of which frameshift mutations were reported in MAPRE3 for gastric and colorectal cancers [84], however, MAPRE2 is upregulated in invasive pancreatic cancer cells [83], demonstrating that the yeast phenomic model could help distinguish causal influence in cases of paralogous gene expression having what appear to be opposing effects on phenotypic response of cancer cells to cytotoxic chemotherapy.
NAM7 is a yeast RNA helicase that was deletion-enhancing for both drugs, though slightly stronger for cytarabine, while its human homologs HELZ, HELZ2, and UPF1, were UES only with cytarabine ( Figure 5C; 1-0-14). HELZ has differential influence in cancer, acting as a tumor suppressor or oncogene [66][67][68][69]. UPF1 downregulation is associated with poor prognosis in gastric cancer and hepatocellular carcinoma, and mutations often occur in pancreatic adenosquamous carcinoma [70][71][72]. Thus, it is possible cytarabine could have efficacy for patients with mutational loss of function in members of this helicase family.
ASF1/ASF1B ( Figure 5C; 2-0.16-1) functions in nucleosome assembly as an anti-silencing factor and is one of the most overexpressed histone chaperones in cancer [85]. The yeast phenomic data suggest that anticancer approaches that target ASF1 as a driver [154] could be augmented by combination with gemcitabine or cytarabine.
AVL9/AVL9 ( Figure 5C; 2-0.16-1) functions in exocytic transport from the Golgi [86]. AVL9 knockdown resulted in abnormal mitoses associated with defective protein trafficking and increased cell migration with development of cysts [87], but also reduced cell proliferation and migration in other studies [88]. Regardless, the yeast phenomic model together with pharmacogenomics data would predict that functional loss of AVL9 renders cells vulnerable to cytarabine.
PMR1 is a P-type ATPase that transports Mn++ and Ca++ into the Golgi. Several of its human homologs, ATP1A1, ATP1A2, ATP1A3, ATP1A4, ATP2C1, were UES, either for gemcitabine or cytarabine, in the PharmacoDB analysis ( Figure 5C; 2-0.16-1). Reduced expression of ATP1A1 can promote development of renal cell carcinoma [89], reduced expression of ATP1A2 is associated with breast cancer [90], and mutations in ATP2C1 impair the DNA damage response and increase the incidence of squamous cell tumors in mice [91,92]. Like with FKH2 (described above), PMR1 deletion-enhancement points to multiple human homologs that are both implicated in the cancer literature to promote cancer when underexpressed, yet are also UES in the pharmacogenomics data, suggesting a potentially clinically useful synthetic lethal vulnerability.
KTI11/DPH3 ( Figure 5C; 1-0-14), is a multi-functional protein involved in the biosynthesis of dipthamide and tRNA modifications important for regulation of translation, development and stress response [80,81], and has promoter mutations associated with skin cancer [82]. It was observed to be UES only for cytarabine and in hematopoietic and lymphoid cancer (the context cytarabine is used clinically).
ACB1 binds acyl-CoA esters and transports them to acyl-CoA-consuming processes, which is upregulated in response to DNA replication stress [155]. Human homologs of ACB1 exhibiting UES ( Figure 5C; 2-0.16-1) included: (1) DBI, which is upregulated in hepatocellular carcinoma and lung cancer, and its expression is negatively associated with multidrug resistance in breast cancer [101][102][103]; (2) ACBD4, which promotes ER-peroxisome associations [98] and is upregulated by a histone deacetylase inhibitor, valproic acid, in a panel of cancer cell lines [99]; and (3) ACBD5, which also promotes ER-peroxisome associations, but its link to cancer is unclear [100]. Thus, it appears this gene family may influence epigenetic processes that buffer the cytotoxic effects of gemcitabine and cytarabine.

Deletion Suppression of Toxicity for Both Nucleosides
As opposed to deletion-enhancing interactions, which represent functions that buffer the cytotoxic effects of the drugs, deletion suppression identifies genes that promote toxicity, thus predicting overexpression sensitivity (OES) in pharmacogenomics data that represent causal tumor vulnerabilities. REMc/GTF identified as deletion-suppressing the GO terms glutaminyl-tRNA(Gln) biosynthesis (1-0-3), the nucleoplasmic THO (2-0.8-1), RNA cap binding (1-0-3), and the NuA3b histone acetyltransferase complexes (1-0-7) (Additional File 5, File C), while GTA identified mitochondrial translational elongation and the nuclear cap binding complex (Additional File 6, File A). However, the respective term-specific heatmaps revealed weak effects and high shift for many of the genes (Additional File 2, Figure S3), highlighting the utility of this phenomic visualization tool for prioritizing findings and leading us to shift our focus to individual yeast-human homologs identified in gene deletion suppressing clusters that were OES in the pharmacogenomics analysis, as detailed below.
Of the autophagy-related complexes, Npr2 and Npr3 form an evolutionarily conserved heterodimer involved in mediating induction of autophagy by inhibition of TORC1 signaling in response to amino acid starvation [216], also promoting non-nitrogen starvation-induced autophagy [217] ( Figure 6A). The RAVE complex (RAV1/2) promotes assembly of the vacuolar ATPase [218,219], which is required for vacuolar acidification and efficient autophagy [220]. Gene deletion strains in the term negative regulation of macroautophagy (PHO85, PHO80, and SIC1) [221], which seemed from the automated assessment to suggest an opposing effect, were less compelling following detailed visualization of the data, due to the associated high shift and cytarabine deletion-enhancing interaction ( Figure 6A).
Regarding the term intralumenal vesicle formation, Vps24 and Did4 are components of the ESCRT-III complex (see Figures 5E and 6A), which functions at endosomes, and the ATPase Vps4 is required for disassembly of the complex [222]. Doa4 interacts with Vps20 of ESCRT-III to promote intralumenal vesicle formation, which also requires BRO1 [223]. Pharmacogenomics correlation revealed UES in cancer cell lines for DID4/CHMP2A/CHMP2B ( Figure 6D; 1-0-10). During autophagy, CHMP2A translocates to the phagophore to regulate separation of the inner and outer autophagosomal membranes to form double-membrane autophagosomes [171]. CHMP2B is a member of the ESCRT-III complex required for efficient autophagy and has reduced expression in melanoma [172,173], raising the hypothesis that gemcitabine could have efficacy in that context.
Other genes involved in autophagy-related processes that had human homologs UES in cancer cell lines included: (1) PEP12/STX12 ( Figure 6D; 1-0-10), a t-SNARE required for mitophagy [180], for which underexpression is associated with risk of recurrence [181]; and (2) VPS30/BECN1, knockdown of which enhances gemcitabine cytotoxicity in pancreatic cancer stem cells [170]. Furthermore, gemcitabine treatment has been found to upregulate autophagy in pancreatic or breast cancer, which buffers drug cytotoxicity as inferred by the combination of gemcitabine with autophagy inhibitors' increased killing of cancer cells [224][225][226]. Thus, autophagy-related findings from the yeast model appear consistent with, and to build upon, previous cancer cell models.

Histone Modification and Chromatin Remodeling
GTF/REMc identified the Hda1 and Set1C/COMPASS (1-0-0) complexes as gemcitabine-specific deletion enhancing, which was confirmed by term-specific heatmaps ( Figure 6B). The Set1C complex has been characterized to have a role in cell cycle coordination [227], which may be reflected by greater deletion enhancing interaction for the K than for the L CPP. The Set1C/COMPASS complex catalyzes mono-, di-, and tri-methylation of histone H3K4, which can differentially influence gene transcription depending on the number of methyl groups added [228][229][230][231], and was implicated by BRE2, SWD1, SWD3, SDC1, SPP1, and SHG1 ( Figure 6B). The SWD1 ortholog, RBBP5, which was UES with gemcitabine in lung tissue (Additional File 8, File C; 1-0-4), is upregulated in self-renewing cancer stem cells in glioblastoma and necessary for their self-renewal, is involved in the epithelial-mesenchymal transition in prostate cancer cells via its role in H3K4 trimethylation, and is upregulated in hepatocellular carcinoma [232][233][234][235]. Furthermore, gemcitabine sensitivity of pancreatic cancer cell lines was enhanced by H3K4me3 inhibition with verticillin A [233].

Peptidyl-Tyrosine Dephosphorylation
REMc/GTF identified peptidyl-tyrosine dephosphorylation (1-0-0), for which the term-specific heatmap (Additional File 7) revealed six genes previously characterized for their requirement in oxidant-induced cell cycle arrest and RNA virus replication [239,240], OCA1-6. Two additional tyrosine phosphatases, YMR1 and PTP1, had similar interaction profiles ( Figure 6C). OCA1-3 deletions enhance growth defects associated with reactive oxygen species or caffeine treatment [239,240], and OCA1-4 and OCA6 are deletion suppressors of the cdc13-1 mutation [241]. Although it does not have a tyrosine phosphatase motif, Oca5 deletion also displayed gemcitabine-specific enhancement, consistent with the other genes annotated to this module ( Figure 6C). However, due to the regulatory nature and limited evolutionary conservation of tyrosine phosphorylation, it is not obvious how to predict functionally homologous genetic modules in cancer cells.

Gemcitabine-Buffering by Non-GO-Enriched Yeast-Human Homologs
Homologs with correlated gemcitabine-specific yeast gene deletion enhancement and cancer cell UES (clusters 2-0.2-1, 1-0-10, and 1-0-0) included the family of nucleoside diphosphate kinases (NDKs) ( Figure 6D; Table 4). A single member of the NDK family, YNK1, exists in yeast, while the human genome encodes several paralogs (NME genes) (Additional File 8, File A). The NDKs transfer the γ phosphate of ATP to nucleoside diphosphate as the final step of purine and pyrimidine nucleoside and deoxynucleoside triphosphate biosynthesis and salvage [248,249]. Thus, NDK appears to modulate gemcitabine toxicity by differential activity for endogenous substrates vs. nucleoside analog drugs. In yeast, deletion enhancement by YNK1 was selective for gemcitabine, however the effects in cancer cells are potentially more complex due to multiple NDK genes. In PharmacoDB, NME3 and 5 were UES for gemcitabine, while NME4, 6, and 7 were OES for cytarabine, implicating differential specificity of NME genes for natural and/or medicinal nucleosides as well as possible influences of other kinases, which have, for example, been shown to act on gemcitabine diphosphate [250]. NME5 overexpression was previously associated with gemcitabine-resistant cancer, and its knockdown can increase gemcitabine efficacy [182]. Thus, the anticancer efficacy of gemcitabine could be influenced by differential expression and activity of NDK isoforms across tissues [251], such that NME gene expression could be predictive of response to nucleoside analogs, or perhaps targeted for synergistic antitumor activity.
KEX2 is the yeast member of the calcium-dependent proprotein convertase subtilisin/kexin type serine proteases, which functions in the secretory pathway. Four of the seven human homologs of KEX2 were UES in the pharmacogenomics analysis ( Figure 6D; 1-0-10), including: (1) PCSK1, which can be downregulated by pancreatic cancer derived exosomes [176], (2) PCSK2, which has reduced expression in lung cancer [177], (3) PCSK5, which is also reduced in lung cancer and, furthermore, when reduced in triple negative breast cancer, leads to loss of the Gdf11 tumor suppressor [177,178], and (4) PCSK7, which has been reported both to have reduced expression in lung cancer and increased expression in gemcitabine-resistant cells [177,179]. Thus, loss of this gene family may create cancer-specific vulnerabilities to gemcitabine cytotoxicity.
NMA1 and its human homologs NMNAT1, NMNAT2, and NMNAT3 are nicotinic acid mononucleotide adenylyltransferases involved in NAD biosynthesis and homeostasis, which were found to be UES for both gemcitabine and cytarabine ( Figure 6D, 1-0-0). Loss of function mutations and underexpression of NMNAT1 are associated with increased rRNA expression and sensitivity to DNA damage in lung cancer cell lines [166], consistent with the hypothesis that they could have deletion-enhancing therapeutic benefit in cancers treated with gemcitabine or cytarabine.
RAD54 is a DNA-dependent ATPase that stimulates strand exchange in recombinational DNA repair, which is a known vulnerability of cancer [252]. The human homolog of RAD54, ATRX, was UES by PharmacoDB analysis (Figure 6D, 1-0-0), and loss of ATRX has been associated with improved response to gemcitabine plus radiation therapy in glioma patients with IDH1 mutations [168].
CLB5, a B-type cyclin, is involved in initiation of DNA replication and G1-S progression, for which promoter hypermethylation of the human homolog, CCNA1, is associated with multiple cancers [157], and which was found to be UES with gemcitabine ( Figure 6D, 1-0-0).

Gemcitabine-Specific Gene Deletion Suppression
Representing this class of gene interaction, pharmacogenomics integration is highlighted for clusters 2-0.8-0 and 1-0-7 ( Figure 6E). Although there was limited gene ontology enrichment, the term phosphatidylserine biosynthetic process (UME6 and CHO1) and the GARP (VPS51-54) and Lem3p-Dnf1p complexes were identified ( Figure 6F, Table 2). Ume6 is involved in both positive and negative regulation of the phosphatidylserine synthase, Cho1 [253,254]. Phosphatidylserine exposure to the plasma membrane is a marker of yeast and mammalian apoptosis [255], the latter of which is induced by gemcitabine [256]. In pancreatic cancer cells, addition of the sphingolipid, sphingomyelin, enhances gemcitabine cytotoxicity through increased apoptosis [256,257]. Moreover, GARP complex deficiency leads to reduction of sphingomyelin [258] and accumulation of sphingolipid intermediates, consistent with the hypothesis that reduced sphingolipid metabolism alleviates gemcitabine-mediated apoptosis. Lem3 complexes with Dnf1 or Dnf2 to form phospholipid flippases at the plasma and early endosome/trans-Golgi network membranes and regulate phosphatidylethanolamine and phosphatidylserine membrane content [259,260], potentially further influencing the apoptotic response. The Lem3-Dnf1 and Lem3-Dnf2 flippases are regulated by the serine/threonine kinase Fpk1 [261], which is also a gemcitabine-specific deletion suppressor ( Figure 6F).
Overexpression of ALDH genes is observed in colorectal and pancreatic cancer [183,184] and is a prognostic marker of cancer stem cells [185].

Cytarabine-Specific Gene Deletion Enhancement
Cytarabine-specific deletion enhancement suggests functions that buffer cytotoxic effects of cytarabine to a greater extent than gemcitabine, potentially informing on differential activities of the drugs. There was no notable GO enrichment by REMc/GTF, but four functions of potential relevance were revealed by GTA ( Figure 7A, Table 2). Two of them, the HIR complex (HIR1-3, HPC2) and sphinganine kinase activity (LCB4, LCB5) were relatively weak, being deletion-enhancing only for the L CPP ( Figure 7A). LCB4/5 homologs that were UES in PharmacoDB included: (1) CERKL (Additional File 8, Files B-C; 1-0-6), a ceramide kinase-like gene that regulates autophagy by stabilizing SIRT1 [264], a gene mentioned above for its inhibition being synergistic with cytarabine against acute lymphoblastic leukemia cells [138], and (2) AGK, which is overexpressed in hepatocellular carcinoma, glioma, breast, and cervical squamous cell cancers [265][266][267][268]. Two stronger interaction modules, evidenced by deletion enhancement for both the K and L CPPs, were protein localization to septin ring (HSL1 and ELM1) and the Sec61 translocon complex (SBH1, SBH2, and SEC61) ( Figure 7A, Table 2). In yeast, Hsl1 and Elm1 are annotated as "bud sensors" to recruit Hsl7 to the septin ring at the bud site to degrade the mitotic inhibitor, Swe1 [269]. The HSL1 homologs, BRSK1 and BRSK2, were UES in the cancer data. BRSK1 is mutated in gastric and colorectal carcinoma [270] and its decreased expression is associated with breast cancer [271], but BRSK2 is overexpressed in pancreatic cancer, where it is AKT-activating [272]. PharmacoDB also identified the SEC61 homolog, SEC61A1, which is upregulated in colon adenocarcinoma tissue [273]. The data descriptions are the same as for Table 3.

Human Genes that have Deletion Enhancing Yeast Homologs and Confer Cytarabine UES
We identified human genes that were UES to cytarabine and homologous to yeast genes in REMc clusters (1-0-9 and 2-0.13-0) displaying a pattern of cytarabine-specific deletion enhancement ( Figure 7B; Table 5). Cancer-relevant examples include: (1) Ptm1, which is a protein of unknown function that copurifies with late Golgi vesicles containing the v-SNARE, Tlg2p, but interestingly, its human homologs, TMEM87A and TMEM87B, were UES for cytarabine and identified in a study focused on cytarabine efficacy in acute myelogenous leukemia [281]. (2) NAP1/NAP1L3/NAP1L4, which is a nucleosome assembly protein involved in nuclear transport and exchange of histones H2A and H2B and also interacts with Clb2, is phosphorylated by CK2, and has protein abundance that increases in response to DNA replication stress [155]. NAP1L3 is overexpressed in breast cancer [280]. (3) CCH1, which is a voltage-gated high-affinity calcium channel with several homologs that were UES, including: CACNA1A, underexpressed in breast, colorectal, esophageal, gastric, and brain cancers; CACNA1B, underexpressed in breast and brain cancers; CACNA1C, underexpressed in brain, bladder, lung, lymphoma, prostate, and renal cancers; CACNA1E, underexpressed in breast, brain, gastric, leukemia, lung, and prostate cancers; and CACNA1F, underexpressed in lymphoma [274]. (4) IZH1, a yeast membrane protein involved in zinc ion homeostasis, having a human homolog, PAQR1/ADIPOR1 that encodes the adiponectin receptor protein 1, which is differentially regulated in breast cancers [278,279]. (5) FAT1, a yeast fatty acid transporter and very long-chain fatty acyl-CoA synthetase that corresponds to SLC27A2 (very long-chain acyl Co-A synthetase), which is underexpressed in lung cancer [275], and SLC27A3 (long-chain fatty acid transport), which is hypermethylated in melanoma [276]. (6) FOL2/GCH1, a GTP-cyclohydrolase that catalyzes the first step in folic acid biosynthesis.

Discussion
Informative phenomic models have been developed for multiple human diseases, including cystic fibrosis, neurodegenerative disorders, and cancer [9,[282][283][284]. Molecular models include mutations in conserved residues of yeast homologs of a disease gene and introduction of human alleles into yeast. Complementation of gene functions by human homologs, and vice versa, has demonstrated evolutionary conservation of gene functions [285][286][287]. Like their basic functions, gene interactions are conserved [288,289] and yeast is unique in its capability to address complex genetic interactions experimentally [290]. Here, we model how yeast phenomic assessment of gene-drug interaction could be employed as part of a precision oncology paradigm to predict efficacy of cytotoxic chemotherapy based on the unique cancer genetic profiles of individual patients.
To model the networks that buffer deoxyribonucleoside analogs, we humanized yeast by introducing deoxycytidine kinase into the YKO/KD strain collection, as yeast do not encode dCK in their genomes, and thus cannot activate the unphosphorylated drugs. We hypothesized that gemcitabine and cytarabine would have different buffering profiles, despite their similar mechanisms of action, due to their distinct anticancer efficacies. Results of the unbiased yeast phenomic experiments confirmed this expectation, revealing distinct, though partially overlapping, gene interaction networks. Differential interaction predominated despite the similarity of the molecules, illustrating that distinct mechanisms for buffering anticancer cytotoxic drug responses can be inferred from yeast phenomics and thus applied to predict how an individual's cancer genome could influence responses to treatment [3,5].
Deletion enhancement of both gemcitabine and cytarabine suggested processes that function to buffer nucleoside analog cytotoxicity in common ( Figure 5), in contrast to buffering mechanisms that acted differentially in response to the drugs. Functionally enriched processes that buffered both drugs to a similar extent included the intra-S DNA damage checkpoint; positive regulation of DNA-dependent DNA replication initiation; vesicle fusion with vacuole; and the Mre11, checkpoint clamp, RecQ helicase-Topo III, CORVET, HOPS, ESCRT, GET, Ubp3-Bre5 deubiquitination complexes.
Among the drug-specific deletion enhancing interactions, autophagy, histone modification, chromatin remodeling, and peptidyl-tyrosine dephosphorylation buffered gemcitabine more so than cytarabine ( Figure 6). There were only a few cytarabine-specific deletion enhancing GO-enriched terms, but there were many individual genes with human homologs having cancer relevance that buffered cytarabine relatively specifically (Figure 7). On the other hand, genes that preferentially promote cytotoxicity were observed primarily for gemcitabine, and enriched functions were related to apoptosis, including phosphatidylserine biosynthesis, and the GARP and Lem2/3 complexes ( Figure 6).
The model we constructed incorporates the powerful pharmacogenomics datasets and analysis tools from PharmacoDB, mining them by integration of yeast phenomic drug-gene interaction experiments. We integrated yeast phenomic and PharmacoDB data to identify, across the respective datasets, correlations between deletion enhancement and underexpression sensitivity or deletion suppression and overexpression sensitivity. Deletion enhancement indicates genes that are biomarkers and synergistic targets to augment drug efficacy and expand the therapeutic window, whereas deletion suppression identifies genes that promote drug cytotoxicity, and thus confer sensitivity when hyper-functional and resistance when deficient. A particularly attractive class of drug-gene interaction is overexpression sensitivity involving driver genes, however, anticancer efficacy could also be conferred by lethal drug-gene interactions involving passenger genes, tumor suppressor genes, or components of genetic buffering networks that become altered due to genomic instability ( Figure 1A). The cancer literature revealed many deletion enhancing/UES and deletion suppressing/OES genes to have roles in cancer, suggesting that integration of yeast phenomic models and pharmacogenomics data could have clinical utility for choosing cytotoxic treatments based on gene expression profiles of individual cancers. While predictions sometimes involved GO-enriched processes, often the genes were identified individually. Assessment of conserved buffering genes (number of yeast deletion enhancers with human homologs exhibiting underexpression sensitivity) was estimated to be around 50% (Additional File 8, File A-see 'Conservation_buffering_homologs').
We focused discussion of conserved homologs with cancer relevance based on the integration of yeast and human data, the extent of homology, the annotated functions, and the existing cancer literature. We note that although genes specifically buffering cytarabine were less well annotated and less conserved overall, they were equally relevant in the model, with about half of all buffering gene interactions in yeast evidenced in the PharmacoDB functional genomic data (Additional File 8, File A.-'Conservation_buffering_homologs' worksheet). Thus, the supplemental materials serve as a resource for future analyses of genes that remain to be annotated (Additional File 8, Files A-D).
As shown by Birrell et al., differential gene expression is a poor predictor of which genes are required for response to a drug [52]. Thus, yeast phenomic models (i.e., Q-HTCP of the YKO/KD library) may help clarify the milieu of potentially causal associations between differential gene expression and drug sensitivity observed in cancer cells by individually testing in yeast the influences on cell proliferation of evolutionarily related protein products. As far as we know, this work represents the first application of this fundamental observation from yeast to systems level experimental data from human cells. Literature-based validations of the yeast phenomic model of nucleoside analogs in human cancer cell lines and other cancer models are exemplified in Table 6. These examples illustrate that integrative, systems level drug-gene interaction modeling employing the experimental power of S. cerevisiae phenomics could be applicable to cancer genomic profiling for systems level, precision oncology. Depletion of BECN1 sensitizes pancreatic cancer stem cells to gemcitabine LCB4/5/CERKL sphinganine kinase activity ceramide kinase like [138,264] CERKL stabilizes SIRT1, SIRT1 chemical inhibition sensitizes acute myeloid leukemia cells to cytarabine Deletion-enhancing/UES drug-gene interactions are highlighted; most exemplify loss of buffering functions that lead to increased drug sensitivity; however, there is one instance (KEX2/PCSK7) of overexpression of the buffering gene that increases drug resistance.
In summary, the yeast phenomic model of nucleoside analog toxicity appears to serve as a valuable resource for interpreting cancer pharmacogenomics data regarding gene-drug interaction that could be predictive of patient-specific chemotherapeutic efficacy. Since it's not possible to collect comparable phenomic information from human populations or cancerous tissue alone [5], systems-level yeast phenomic models can help expand and integrate relevant (i.e., evolutionarily conserved) aspects of the extensive cancer literature with regard to cancer-specific vulnerabilities to cytotoxic therapies. A deeper understanding of how genomic instability influences the genetic network that buffers chemotherapeutic agents like nucleoside analogs could guide future research to personalize anticancer therapies based on cancer genomic profiles unique to individual patients. Thus, a future direction for this work should include development of algorithms that prospectively predict chemotherapy response in individual patient cancer cells, which could be tested as part of a prognostic evaluation. In this initial study, we focused on expression data from PharmacoDB, rather than mutation data, because differential gene expression in cancer cells is more analogous to the quantitative changes in gene expression resulting from gene KO/KD than qualitative or indeterminate effects resulting from point mutations, for example. However, an interesting future direction could be to analyze mutation data in conjunction with gene expression and yeast phenomic data, for example to identify eQTLs associated with UES or OES genes [291]. With regard to the number of interactions observed in the yeast model, we note this depends on the inclusiveness of homology. S. cerevisiae, as a single cell organism and evolutionarily distant relative of humans, is informative about gene interactions across many different human cell types. Thus, yeast genetic interactions are merely hypothesis-generating and require appropriate future testing of human homologs and in human cell types. While yeast phenomic models of human disease are powerful engines of discovery, findings can only be truly prioritized and focused at this time in models more directly relevant to the human conditions of interest.

Conclusions
A humanized yeast phenomic model of deoxycytidine kinase was developed to map drug-gene interactions modulating antiproliferative effects of nucleoside analogs in a eukaryotic cell and to investigate the relevance of the resulting networks for precision oncology by integration with cancer pharmacogenomics-derived associations between gene expression and cancer cell line drug sensitivity. The yeast phenomic model revealed gene-drug interaction for the two deoxycytidine analogs, gemcitabine and cytarabine, to be largely different, consistent with the distinct types of cancer for which they are used clinically. The model overall suggested evolutionary conservation of drug-gene interaction that could be used as a resource to predict anticancer therapeutic efficacy based on genetic information specific to individual patients' tumors. Yeast phenomics affords a scalable, high-resolution approach to model, at a systems level, the genetic requirements for sensitivity and resistance to cytotoxic agents and, thus, the potential to resolve complex influences of genetic variation on drug response more accurately. Global and quantitative models of the distinct genetic buffering networks required to maintain cellular homeostasis after exposure to chemotherapeutic agents could aid precision oncology paradigms aimed at identifying composite genomic derangements that create enhanced cancer cell-specific vulnerabilities to particular anticancer drugs. Further in this regard, cytotoxic chemotherapeutic agents are used in combination, so another direction for yeast phenomic analysis of anticancer agents would be to characterize clinically relevant drug combinations.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/10/770/s1, Additional File 1: Supplemental tables: Tables S1-S6. Table S1: Primers used in construction of the synthetic genetic array (SGA) 'query' strain (MIY16) for doxycycline-inducible deoxycytidine kinase expression (see Figure S1), Table S2: Drug-gene interaction data from genome-wide experiments for gemcitabine YKO/KD strains, Table S3: Drug-gene interaction data from genome-wide experiments for gemcitabine reference strains, Table S4: Drug-gene interaction data from genome-wide experiments for cytarabine YKO/KD strains. Table S5: Drug-gene interaction data from genome-wide experiments for cytarabine reference strains, Table S6: Summary statistics of interaction z-scores for the YKO/YKD and reference cultures from the phenomic analyses of gemcitabine and cytarabine drug-gene interaction, including the ranges, averages, standard deviations, and representative interaction z-scores in the tails of the distributions. Additional File 2: Supplemental figures. Figure S1: Construction of tet-inducible dCK allele, Figure S2: Reference strain distributions for rate (r) and area under curve (AUC) with gemcitabine or cytarabine treatment, Figure S3: Gene deletion suppression modules exhibiting high shift or variable interaction, Figure S4: Elongator holoenzyme complex, protein urmylation, and tRNA wobble position uridine thiolation gene modules exhibit variable shift and deletion enhancing interaction profiles. Additional File 3: Gene-gemcitabine interaction plots for (A) YKO, (B) KD, and (C) reference strains. Additional File 4: Gene-cytarabine interaction plots for (A) YKO, (B) KD, and (C) reference strains. The first two pages in Additional Files 3 and 4 display the reference distributions used to calculate gene-drug interactions. Summary data comparing them are in Additional File 1, Table S6. Additional File 5: REMc results: File A is a table providing REMc results and associated gene interaction and shift data. File B is the heatmap representation of each REMc cluster after incorporating shift values and hierarchical clustering. Gene label colors for YKO and YKD are black and red, respectively. File C contains tables of GTF results from the process, function, and component ontologies for each REMc cluster. Additional File 6. Gene ontology term averaging (GTA) results and interactive plots. File A contains all GTA values, cross-referenced with REMc-enriched terms. Files B and C display GTA values associated with above-threshold GTA scores plotted for gemcitabine vs. cytarabine for L and K, respectively. They should be opened in an Internet web browser so that embedded information from Additional File 6A can be viewed by scrolling over points on the graphs. Subsets in each of the plots can be toggled off and on by clicking on the respective legend label. In the embedded information, X1 represents gemcitabine and X2 represents cytarabine information. Additional File 7: GO term-specific heatmaps for REMc/GTF-enriched modules, generated as described in methods and Figure 3. Child terms are presented in subsequent pages of the parent file name. GO terms having more than 100 children, with two or fewer genes annotated to the term, or a file size over 400 KB are not represented. All heatmaps are generated with the same layout (see Figure 3). Additional File 8: Integration of yeast phenomic drug-gene interaction data with cancer cell line pharmacogenomic data (gene expression and drug sensitivity correlations) to predict differential gene expression in cancer cells that is causally associated with enhanced gemcitabine or cytarabine cytotoxicity. File A contains three tables of all UES and OES human genes from the GDSC (gemcitabine and cytarabine) and gCSI (gemcitabine only) databases, also indicating whether deletion-enhancing or deletion-suppressing gene-drug interactions from phenomic analysis of yeast homologs predicted causal associations. Homology relationships and other associated data are coassembled. Files B-D consist of heatmaps and corresponding tables of yeast drug-gene interaction sets that predict causality for UES or OES human genes identified (File B) across all tissue, (File C) in lung, or (File D) in hematopoietic and lymphoid tissue. See Figure 3 for explanation of colors for gene labels and homology type. Note: the teal color, which represents cytarabine-specific UES/OES in the heatmaps in the main manuscript figures corresponds to darker blue in the supplemental heatmaps, while gold, representing gemcitabine-specific UES/OES in the main manuscript corresponds to bright yellow. Black corresponds to UES or OES for both drugs).

List of Abbreviations and Glossary of Terms
AraC-cytarabine; cytosine arabinoside. CPPs-cell proliferation parameters: parameters of the logistic growth equation used to fit cell proliferation data obtained by Q-HTCP. The CPPs used to assess gene interaction in this study were K (carrying capacity) and L (time required to reach half of carrying capacity) [7][8][9]37]. DAmP-decreased abundance of mRNA. Production: refers to a method of making YKD alleles, where the 3 UTR of essential genes is disrupted, reducing mRNA stability and gene dosage [292]. dCK-deoxycytidine kinase. dCMP-deoxycytidine monophosphate. DE-deletion enhancer: gene loss of function (knockout or knockdown) that results in enhancement/increase of drug sensitivity [9]. dFdC-2 ,2 -difluoro 2 -deoxycytidine, gemcitabine. dNTP-deoxyribonucleotide triphosphate. DS-deletion suppressor: gene loss of function (knockout or knockdown) that results in suppression / reduction of drug sensitivity [9]. ESCRT-endosomal sorting complex required for transport. GARP complex-Golgi-associated retrograde protein complex. gCSI-the Genentech cell line screening initiative: one of two pharmacogenomics datasets used in this study (https: //pharmacodb.pmgenomics.ca/datasets/4). GDSC1000-Genomics of Drug Sensitivity in Cancer: One of two pharmacogenomics datasets used in this study (https://pharmacodb.pmgenomics.ca/datasets/5). GO-gene ontology. GTF-gene ontology term finder: an algorithm to assess GO term enrichment amongst a list of genes; applied to REMc (clustering) results [41]. GTA-gene ontology term averaging: an assessment of GO term function obtained by averaging the gene interaction values for all genes of a GO term. GTA value-gene ontology term average value. gtaSD-standard deviation of GTA value. GTA score-(GTA value -gtaSD). HaL-hematopoietic and lymphoid tissue. HDAC-histone deacetylase complex. HLD-human-like media with dextrose [8]: the yeast media used in this study. INT-interaction score. NDK-nucleoside diphosphate kinase. OES-overexpression sensitivity: refers to association of increased gene expression with drug sensitivity in pharmacogenomics data [32]. PharmacoDB-the resource used for cancer pharmacogenomics analysis [32]. PPOD-Princeton protein orthology database. Q-HTCP-quantitative high-throughput cell array phenotyping: a method of imaging, image analysis, and growth curve fitting to obtain cell proliferation parameters [7,37]. Ref-reference: the genetic background from which the YKO/KD library was derived. REMc-recursive expectation maximization clustering: a probabilistic clustering algorithm that determines a discrete number of clusters from a data matrix [40]. RNR-ribonucleotide reductase. SD-standard deviation. SGA-synthetic genetic array. SGD-saccharomyces genome database. UES-underexpression sensitivity: refers to association of reduced gene expression with drug sensitivity in pharmacogenomics data [32]. YKO-yeast knockout. YKD-yeast knockdown: DAmP alleles. YKO/KD-yeast knockout or knockdown.