Overview of NIBR and HIPLAB screens
Because all our comparisons are based on the ability to compare both datasets, we describe each dataset in detail. The data processing strategies of the raw data were fundamentally different between the two research sites (Table 1). In the HIPLAB dataset, the raw data was normalized separately for the strain-specific uptags and downtags, independently for the heterozygous and homozygous strains, creating 4 sets of results: uptag/het, uptag/hom, downtag/het, downtag/hom. For each set, logged raw average intensities were normalized across all arrays using a variation of median polish that incorporates batch effect correction [6]. Because the performance of the two tags in each strain can vary significantly, a ‘best tag’ was identified for each strain, defined as the tag with the lowest robust coefficient of variation across all of the control microarrays. For each array, tags were removed if they did not pass the computed compound and control background thresholds, calculated from the median + 5MADs of the raw signal from the unnormalized intensity values of the used (corresponding to strain tags) and unused (control) features on the array across all arrays. In contrast, in the NIBR dataset, arrays were normalized by “study id”, (a set of ~40 compounds) but were not corrected for batch effects. Rather, tags that performed poorly, based on their correlation values of uptags and downtags across different intensity ranges in the control arrays, were removed and the remaining tags were averaged to obtain strain intensity values.
In the HIPLAB dataset, relative strain abundance was quantified for each strain as the log2 of the median signal in the control condition divided by the signal from the compound treatment. The final fitness defect (FD) score is expressed as a robust z-score where the median of the log2 ratios for all strains in a given screen is subtracted from the log2 ratio of a specific strain and divided by the MAD of all log2 ratios for all strains in that screen. In the NIBR dataset, the inverse log2ratioHIPLAB was used with three differences: 1) average intensities of controls were used (instead of median signals) and 2) because NIBR used replicates for each compound, the average of signals of the compound samples were used instead of a single value (Table 1) and 3) the final gene-wise z-score normalizes for median and standard deviation of each strain across all experiments using quantile estimates (see Methods).
Both laboratories constructed pools of heterozygous and homozygous strains in a similar manner and collected samples robotically for both the HIP and HOP assays as previously described [8]. For NIBR experiments, samples were collected at fixed time points (which served as a proxy for the number of cell doublings), whereas in the HIPLAB experiments cells were collected based on actual doubling time. Notably, in the NIBR pools, ~300 strains fewer homozygous deletion strains were detectable compared to the HIPLAB pools. These strains correlate with strains known to be slow- growers in the absence of drug [9] and their absence is likely due to the fact that the pool was allowed to grow overnight (~16hrs) in the NIBR assays, during which slow-growing strains drop out before the start of the experiment.
Another difference between protocols was that NIBR screened all heterozygous strains, deleted for both essential and nonessential genes, while the HIPLAB screened only the essential heterozygotes. We decided against screening non-essential heterozygotes based on the following logic: because the concept of the HIP assay relies on a fitness defect resulting from gene dosage being decreased from two copies to one in a heterozygous diploid deletion strain, it follows that such fitness defects should not be observed if that gene is not required for growth, as is the case for nonessential genes [10]. Indeed, we find in practice that the HIP profiles of the nonessential heterozygotes do not correlate with HOP profiles for the same drug, nor are these nonessential heterozygote profiles biologically informative. This is illustrated by the profiles for DNA damaging agents. For example, In HOP screens of nonessential deletion strains, RAD genes have high FD scores in the presence of a DNA damaging agent (mechlorethamine), but none of these strains were sensitive as heterozygotes (Figure S1). The exception to this is the small number of nonessential heterozygous strains that exhibit severe fitness defects as homozygotes. As these strains exhibit ‘nearly essential’ phenotypes as homozygotes, they would be expected to exhibit drug-induced haploinsufficiency as heterozygotes and therefore should be included in the HIP assay.
We next compared the depth and breadth of each screening dataset. The NIBR screening library included 1641 propriety compounds and 135 reference compounds with known mechanisms of action. In total there were 2956 HIP and 2923 HOP experiments, for 1776 discrete chemical structures.
However, because NIBR HIP screens included heterozygous strains deleted for both essential and nonessential genes (as mentioned above) when we combined the HIP and HOP NIBR datasets for shared compounds 2725 full HIPHOP screens spanning 1771 distinct compounds remained. ~56% of the NIBR screening library, however (representing 596 compounds) could practically be considered replicate screens because they exhibit correlations on par with true replicates, even though they were screened at different concentrations. For example, we observe such “practical replicates” when a particular compound is screened at a different concentration, yet the level of inhibition is comparable. Supporting this observation, the majority of such “replicates” clustered together (~65%; those with more than one replicate are included if at least one pair is clustered together). Given these experimental caveats, the informative datapoints were reduced from ~30 million to ~15 million due to the nonessential heterozygotes, and to ~9 million unique datapoints if the replicates were excluded.
The HIPLAB screening library comprised 3356 screens and 3250 unique compounds selected from a set of > 50,000 maximally diverse small molecules (~20 million data points) with unknown mechanisms and ~characterized drugs or chemical probes.
The structural diversity of the screening libraries reflects the scale of a large screening effort. While NIBR did not provide the compound structures of their libraries, they reported that 50% of the pairwise comparisons between compounds had Tanimoto coefficients less than 0.1 [7]. In comparison, the HIPLAB compounds were of lower diversity; ~43% of the pairwise comparisons had Tanimoto coefficients less than 0.1. Because the NIBR structures were not provided, however (with the exception of 135 reference compounds and 15 novel inhibitors) this claim is not verifiable [7].
Coinhibition between chemogenomic profiles
To compare the HIPLAB and NIBR screens, we first compared ~150 chemogenomic profiles representing ~50 reference compounds with known MoA that were screened by both NIBR and HIPLAB (Table 2). For many of these compounds, the drug target is well-established in yeast. Chemogenomic profiles were compared individually using ‘coinihibition’ values, where coinhibition is defined as the degree of similarity between two chemogenomic profiles, i.e., the FD scores across all genes in each screen, using Pearson correlation as a metric. The HOP profiles for the mechlorethamine, a DNA damaging agent, identified a similar set of DNA repair genes including RAD1, RAD2, RAD4, RAD5, RAD10, RAD14, RAD18, REV7, REV3, SRS2 and PSO2 (Figure 1A). Likewise, we did a pairwise comparison of four nocodazole chemogenomic profiles exhibiting between-drug correlations of 0.48 and greater across the entire set of deletion strains (Figure 1B). These correlation values increased when comparing only those individual genes exhibiting significant FD scores in the NIBR and HIPLAB in the HIP nocodazole profiles (Figure 1C). In this case, the HIP genes identified are enriched for genes required for tubulin folding (CCT genes). Finally, based on a correlation value of > 0.5 with the nocodazole profiles, we highlight the HIP profiles of two novel compounds, NIBR 2667 and HIPLAB 5790901, both identifying a nearly identical set of genes (Figure 1D). It should be noted that because the screens were performed at different concentrations, a linear correlation of one is not expected.
In addition to measuring the correlation between chemical profiles, a valuable metric is the correlation of gene fitness scores across compounds and between datasets. For this comparison we employed ‘cofitness’; the degree of similarity between fitness profiles in which the FD scores between two genes are measured across all compounds, using Pearson correlation as a metric. Genes that exhibit a high degree of correlation or cofitness between fitness profiles across compounds are often functionally related. In this case, where we have two independent datasets, we expect the same gene to be cofit across the 50 compounds that were shared between the two datasets. Overall, we observed an overall correlation of ~0.15 between the same gene. Because most genes are not perturbed in any given experiment, we expect that those with highly variable scores (therefore more likely to be more biologically informative) to exhibit greater correlation. This is indeed the observation, when only significantly sensitive strains are considered (genes with standard deviations in the top 5%) the correlation between genes increases to ~0.5. As the correlation increases, the mechanistic similarity between drugs that significantly perturb a given deletion strain also increases. For example, while the IDP1 gene profiles exhibit a similar pattern of perturbation (R-value ~0.4) (Figure 2A), RAD5 and HMG1 exhibit higher correlations (R-value ~0.7, ~0.9, respectively) and significant perturbations are seen in mechanistically related compounds such as 1) the DNA damaging agents’ hydroxyurea, mechlorethamine and methyl methanesulfonate (MMS) and 2) the sterol pathway inhibitors fluconazole and fluvastatin (Figure 2B, 2C). Similarly, in the case of TOR1 (R-value 0.85), outlier fitness deviations all arise from the same compound (rapamycin) (Figure 2D). When we examined the genes exhibiting cofitness within the shared 50 compounds, we observe enrichment for pairs in both sets that reflect the mechanistic enrichment of the compounds as a whole. For example, 6 of the 50 compounds were DNA damaging agents, and as a result, several of the top cofit genes were pairs where both genes were involved in DNA damage.
To examine the agreement between the NIBR and HIPLAB datasets at a more comprehensive level, we combined, and then hierarchically clustered the two HIPHOP datasets together for a subset of mechanistically related compounds. In the first case ~100 compounds representing 19 distinct mechanistic classes including TOR signaling, microtubule poisons, FAS1 inhibitors, cell wall inhibitors, statins, ionophores, ion channel blockers, azoles and morpholine antifungals, (Figure 3A), and in the second case ~40 DNA damaging agents representing eight mechanistic classes including the drug and tool compounds doxorubicin, camptothecin, hydroxyurea, mechlorethamine and MMS (Figure3B). In the resulting heatmaps, in both cases, the two identical dendrograms reveal that the screens cluster primarily by the mechanism of drug action and not by the research institute. Screens from NIBR and HIPLAB were interspersed, and all replicates and compounds with the same mechanism clustered together. In the DNA damaging clustergram, one notable exception was observed for two aclarubicin profiles (one from each research site) that did not cluster with the other anthracycline compounds including doxorubicin, daunorubicin and epirubicin. These differences between specific anthracyclines likely reflect true mechanistic differences between these closely related compounds [11, 12]. For example, the individual aclarubicin HIPHOP profiles implicate RPO31 (encoding an RNA polymerase III subunit) as a potential target. In select cases, compounds with similar mechanisms (i.e., part of the same pathway) also clustered together, including the morpholine antifungals, e.g. fenpropimorph and amorolfine, both targeting ERG2, the azoles, e.g. fluconazole and clotrimazole, targeting ERG11, and the statins, e.g. atorvastatin and fluvastatin, targeting HMG1, with all three targets in the sterol biosynthesis pathway). Other examples include clustering of ion channel blockers next to ionophores, amiodarone and nigericin, respectively, and clustering of rapamycin next to caffeine, known to target the TOR pathway in Saccharomyces cerevisiae (Figure 3A).
Common response signatures
Our previous global analysis of the HIPLAB dataset [6] revealed that, despite the complexities of pharmacological inhibition, the cellular response to small molecules is limited and can be described by a network of 45 major response signatures. These responses comprise chemogenomic profiles with; 1) a characteristic gene signature, 2) distinct GO enrichments and 3) enriched chemical sub-structures. In this 2014 study we found that, by subsampling, the majority of these signatures (~80%) could be identified after screening less than 30% of the compounds, suggesting that the cellular response to small molecules is limited. To test if these response signatures are also present in the NIBR dataset, we used the same methodology as in Lee et al., (2014)[6] to hierarchically cluster the NIBR screens using coinhibition as a distance metric and a dynamic branch cutting method [13] to generate discrete clusters. 96 robust clusters were initially identified covering ~41% of the profiles. Compared to the 45 major responses in the HIPLAB dataset covering ~36% of the profiles, the number of NIBR response signatures was two-fold greater. However, many of these signatures were redundant with respect to their GO enrichments and associated gene signatures. Gene signatures were also longer, as would be expected when clusters are small and when compounds within a cluster are replicates. While it is not entirely clear, one explanation for this observation is that the NIBR screening library contains a large number of replicates (56% of all screens) which would produce partially redundant clusters. To identify discrete clusters with minimal redundancy, dynamic branch cutting parameters were modified to be less sensitive to smaller clusters (see Methods), which resulted in a final set of 42 robust NIBR clusters, comparable to the 45 HIPLAB response signatures. The median number of genes in the response signatures was similar between the final NIBR signatures (7 genes) and HIPLAB signatures (8 genes). Using the overlap between gene signatures to measure similarity between response types, we found that ~66.7% of the 45 major HIPLAB response types were detectable in the NIBR clusters. These common signatures include; iron & copper homeostasis, cell wall signaling, mitochondrial stress, and perturbation of the plasma membrane. More specific responses, often including drugs of known mechanism, included the responses: unfolded protein, anthracycline transcription coupled DNA repair, azoles and statins, ERAD & cell cycle, heme biosynthesis & mitochondrial translocase, NEO1-PIK1, tubulin folding & SWR complex, superoxide and DNA damage.
The majority of these conserved chemogenomic response signatures are enriched for biological processes, details of which can be visualized at the accompanying website (Comparative chemogenomics). Taken together, these results provide further support for the concept that the cellular response to small molecules is limited and that it can be defined by chemogenomic signatures.
Because the chemogenomic signature comparison may be impacted by biases in screening library composition, we asked which of the final 42 responses were unique to the NIBR dataset. Response signatures that were not detectable in the HIPLAB responses included those comprising the three TOR signaling clusters, the GPCR inhibitor response as well as the eukaryotic translation initiation factor (eIF) complex inhibitor signature. Other responses were also gene/target-specific and included: inhibitors of VRG4, encoding a Golgi GDP-mannose transporter, RPL15A & SPP41, encoding a ribosomal gene and a regulator of spliceosome components, respectively, and FAS1, encoding fatty acid synthase. The finding that these ‘missing’ responses likely reflect NIBR screening a small number of target- focused compound sets, combined with our initial finding of many, small and highly redundant clusters suggests that the NIBR libraries are enriched for sub-libraries of mechanistically and/or structurally related molecules.
We used the same approach to compare the signatures of the combined dataset to the HIPLAB responses. In this case, of the resulting 47 chemogenomic signatures, ~84% of the original 45 HIPLAB signatures were detected. Of these 38 overlapping signatures, common signatures included all the DNA damage responses, as well as the azole & statin, superoxide, tubulin folding & SWR complex, unfolded protein and mitochondrial-specific stress responses. Interestingly, by combining the two datasets, some of the NIBR signatures that had not previously matched a HIPLAB response were merged into one of the 38 overlapping responses. Only three signatures were comprised solely of HIPLAB profiles: the NEO1, ubiquinone biosynthesis & proteosome, and the RSC complex & mRNA processing signatures. Conversely, the signatures driven by NIBR profiles largely overlapped the target-specific responses unique to the NIBR dataset including: TIM54, RPL15A & SPP41, VRG4, eIF, and GPCR inhibitors as well as the major TOR signaling response.
Target frequency comparison
Compared to the HIPLAB dataset, which focused on screening diverse compounds with unknown mechanisms, NIBR clusters were highly enriched for screens identifying genes as potential drug targets. The most frequently identified targets that dominated specific clusters in the NIBR dataset include, 1) ERG11 (of the sterol biosynthesis pathway) and KOG1, AVO1, and TOR2, encoding subunits of the Targets of Rapamycin (TOR1 and TOR2) complexes, 2) FAS1, encoding fatty-acid desaturase, and 3) the mitochondrial transport gene TIM54. The coherence of these signatures suggests that the contributing compounds represent structural analogs. In the azole & statin and ERG11-GCN responses, ERG11 is identified as the target in 31% and 67% of the screens, respectively. In the three rapamycin clusters, the TOR1 and TOR2 subunits are identified as targets in over half (26) of the 51 screens. In 2012, NIBR published a study of novel Erg11 inhibitors, suggesting these published inhibitors may be present in the NIBR screening library [14]. Similarly, the high frequency of targeting mTOR (mammalian target of rapamycin) complexes (as evidenced by the three responses associated with TOR signaling) suggests an enrichment of rapamycin analogs in the NIBR compound library. This is consistent with the fact that Rapamycin and aging are active areas of inquiry at NIBR [15–17].
HIPHOP profiles presented in studies previously published by NIBR researchers allowed us, in select cases, to infer the structure of blinded screens. For example, a Nature Chemical Biology study published by the group demonstrated TIM23-dependent mitochondrial import as the target of the natural product stendomycin [18]. A HIPHOP profile in the NIBR dataset was nearly identical to the published version, particularly after accounting for differences in concentration (Figure S2). Similarly, the NIBR group also published a novel geranylgeranyltransferase inhibitor (uncovering sensitivities of strains encoding subunits of the CDC43/RAM2 heterodimer) that was highly correlated to the HIPHOP profile for NIBR compound 5692 in the NIBR dataset [19] (Figure S3).
Compounds and mechanism of action inferred by clustering with reference compounds
One of the NIBR clusters revealed the mechanism of NIBR compounds by virtue of its correlation to reference compounds. The HIPLAB amphotericin B HIP screen was highly correlated with the NIBR 4247 and 1020 HIP screens (> 0.7, p-value < 1e-16) (Figure 4A). In another example, the NIBR compounds 1208, 1209, 1210 and 1211 exhibited correlations of > 0.8 (p-value < 1e-16) with the hydroxyurea screens, a correlation value on par with that observed between replicates, suggesting these compounds are most likely structural analogs of hydroxyurea or closely related derivatives (Figure 4B).