Connecting Small Molecules with Similar Assay Performance Profiles Leads to New Biological Hypotheses

Data Sets

Computing biological profiles for small-molecule modulators in cell-based phenotypic profiling experiments requires normalization of measurements beyond what is required to choose compounds for follow-up from a single high-throughput screen. Our strategy is to compare each measurement of a small-molecule perturbation with an appropriate negative-control distribution that best reflects biological and technical noise inherent in the assay. We use this distribution to compute a dimensionless score, intuitively similar to a z score, which for each compound treatment relates measured values to the likelihood that a measurement can be explained by noise. We applied a variation of this approach to raw, single-concentration, high-throughput screening data in ChemBank ¹⁹ as well as data in CBIP (Chemical Biology Informatics Platform), an internal HTS database at the Broad Institute. From ChemBank, we extracted a total of 8.36 million normalized assay measurements, including both compound and control well data. We excluded assay development experiments, measurements of compound autofluorescence in the absence of a biological sample, and measurements from which other included measurements (e.g., differences and ratios) were derived. These data represent exposure of 1212 distinct compound stock plates to 1015 distinct biological assay conditions differing in at least one of biological sample, compound exposure time, concentration, or assay readout. These data are very sparse, with an average coverage of the potential stock plate × assay space of 2.16% (26,511 of 1,228,968 possible stock plate × assay combinations). From CBIP, we extracted more than 6.1 million results for 83 observations in 24 high-throughput screens. We combined these results in a custom database (Suppl. Fig. S1) that we query on demand (Suppl. Fig. S2) to generate assay performance profiles.

Normalization of Assay Results

We expressed HTS measurements as a dimensionless score D representing a normalized weighted average of deviations from appropriate negative-control distributions ( Fig. 1 ). Let x₁, . . ., x_m be all measurements for a single compound in one concentration and one assay outcome (we distinguish between assay and assay outcome, as some assays involve multiple measurements). We required that all measurements were obtained in the presence of negative (DMSO-treatment) controls and assumed that these controls behave according to a normal distribution. To accommodate run-to-run (or batch-to-batch) variability, for each measurement x_p, p = 1, . . ., m, we considered a normal distribution with mean µ _l and standard deviation σ _l for negative controls in the run (or batch) corresponding to x_l. We aggregated assay measurements into a weighted average of background-subtracted values,

x ¯ = w 1 ( x 1 − μ 1 ) + … + w m ( x m − μ m ) w 1 + … + w m ,

OPEN IN VIEWER

Figure 1.

Data preprocessing reduces bias in assay performance profiles. Double-sigmoid transformation (blue trace) of compounds scored relative to appropriate DMSO-control distribution (red trace) is performed to suppress activity differences among measurements in the noise and to normalize the contributions of assays with different dynamic ranges.

where the weights w_p, p = 1, . . ., m, are inversely proportional to the variance of the measurement and vary across plates and runs. We scaled this weighted average by an estimate of the uncertainty in the distance between the weighted average and the negative-control distribution to obtain D:

D = x ¯ σ C 2 + σ M 2 ,

where σ C 2 is the variance of the weighted mean x ¯ and σ M 2 is the variance of the background-subtracted measurements from negative-control–treated wells.

We conceptualize assay results as a large (sparse) matrix with rows corresponding to compounds and columns corresponding to assays. To be more precise, since a single assay may provide multiple measurements, columns correspond to individual assay outcomes, and therefore a single assay may be represented by multiple columns. For mining primary high-throughput assay data, we assumed that for a given assay, all compounds were tested at an appropriate single concentration. Let {c₁, c₂, . . ., c_n} be a set of compounds for which we would like to determine biological similarities based on performance. Let A = {A₁, A₂, . . ., A_N} be a set of applicable high-throughput biological assay outcomes. For an assay outcome A_i and a compound c_j, we have activity a_ij = A_i(c_j). Because of the sparse nature of the data matrix, for some assays and compounds, the activities a_ij are undefined. For compounds tested in multiple concentrations, we use the most extreme value of D.

Assay Performance Profile Similarity

We wished to avoid some assays dominating similarity scores due to differences in assay dynamic range. D scores can be arbitrarily large or small, which can have adverse effects on the accuracy of profile correlation. Specifically, assays with larger dynamic ranges will receive greater weight if D score variability is not considered. Moreover, we wished to avoid assigning importance to similarities between compounds because they both fail to score in the same assay. Therefore, we applied a double-sigmoid transformation that, intuitively, will give values near 1 (or −1) to “active” compounds and values near zero to “inactive” compounds ( Fig. 1 ). The transformed values b_ij are given by

b i j = ( a i j α ) K 1 + ( a i j α ) 2 K ,

where the parameter α controls the width of the central flat region and the (odd) integer parameter K controls the slope. In our analysis, we use K = 3 and α = 2.3538, the latter selected in such manner that the central 95th percentile interval of the negative-control distribution is transformed into the b_ij interval (−0.5, 0.5).

To compute assay performance profile similarity for a pair of compounds c_j and c_k, we find all assays A_i such that both a_ij = A_i(c_j) and a_ik = A_i(c_k) are defined, and collect them into a subset A′ _⊆ A of size N′ = N′(c_j, c_k). Let b′ _ij be a double-sigmoid–transformed value for an assay from A′ and let b′ _j be a vector of all such values for compound c_j. Similarly, vector b′ _k consists of transformed values for compound c_k and assays from A′. The level of similarity between c_j and c_k can be expressed by the Pearson correlation coefficient r_jk = cor(b′ _j , b′ _k ). To assess the significance of r_jk, we note that the Fisher transformation ²⁰ of the correlation coefficient F(r) = ln[(1 + r)/(1 – r)] / 2 approximates a normal distribution with standard error (N′ – 3)^−1/2, where N′ is the length of vectors b′ _j and b′ _k – the number of assays in which both compounds c_j and c_k were tested. This relationship allows us to express assay performance profile similarity as a z score, z_jk = F ( r j k ) N ′ − 3 .

The z score formula implies that two compounds must have at least four common assays to get a valid value. However, if one of the compounds has only “inactive” results, that noise may appear correlated to the real signal of the other compound. Since the values b′ _ij are bounded, we can circumvent such a false correlation by requiring that the two compounds are “active” in at least one common assay. Formally, we let H′ = H′(c_j, c_k) be the number of coordinates i in vectors b′ _j and b′ _k such that |b′ _ij | ≥ 0.5 ∧ |b′ _ik | ≥ 0.5, and we require that H′ > 0.

Significance Thresholds

To determine what levels of transformed correlations represent meaningful similarities, we used the familiar chemical similarity principle: chemically similar compounds will often have similar biological behavior. The relationship between assay performance profile similarity and chemical similarity can be expressed as enrichment of chemically similar compounds among biologically similar compounds, and we can determine when such enrichment is maximal. We also explored how N′ and H′ (the number of common assays and the number of assays in which both compounds were called active) affect the enrichment.

Due to the large size of the overall compound collection (more than 400,000 compounds), it is not practical to compute all pairwise chemical structure similarities. We randomly selected more than 3 × 10⁷ random pairs of compounds (satisfying N′ > 3 and H′ > 0) and computed their chemical structure and assay performance profile similarities. For chemical structure similarity, we used a Tanimoto similarity score between nonfolded extended connectivity fingerprints with maximum bond-distance 6 (ECFP6). ²¹ For a given Tanimoto similarity threshold T_C and a variable assay performance profile similarity threshold T_B, we can construct a receiver operating characteristic (ROC) curve by plotting the true-positive rate (TPR) against the false-positive rate (FPR), where TPR = TP/(TP + FN), FPR = FP/(FP + TN); TP = compounds similar by both measures, FN = compounds similar by structure only, FP = compounds similar by performance only, and TN = compounds similar by neither measure. ROC curves are routinely used to evaluate the quality of predictive models. ²² Our situation is slightly different since we do not expect similarity of biological activity to be perfectly explained by chemical structure similarity. Nevertheless, we can compare ROC curves using an area under the curve (AUC) that increases with higher enrichment of chemically similar compounds among biologically similar compounds.

Clearly, the more assays two compounds have in common, the more reliable a prediction of a common mechanism of action will be. Therefore, we explored what effect the values of N′ and H′ have on the ROC curve (Suppl. Table S1). We observe no enrichment for N′ ≤ 11 and H′ ≤ 2 and therefore we require that two compounds have at least 12 assays in common, or are “active” in at least 3 common assays, to compute their assay performance profile similarity.

To determine the optimal threshold T_B for assay performance profile similarity, we compute precision, the proportion of TP among all pairs with transformed correlations above T_B; precision = TP/(TP + FP) as a function of T_B ( Fig. 2 ). We fit a polynomial of degree 3 to the precision values and found its local maximum at T_B = 10.0442. Thus, in practice, we take correlation scores z_jk > 10 to be those most likely to be significant when considering compounds with chemical structure information but no prior biological annotation.

OPEN IN VIEWER

Figure 2.

Chemical similarity enrichment provides guidance to select thresholds for assay performance profiling. We use precision, the ratio of “true” positives to all positives, to measure the enrichment of chemically similar compounds among compounds with assay performance profile similarity. Maximal chemical similarity enrichment was achieved for T_B ~10. A third-degree polynomial y = 0.00001873 x³ − 0.001363 x² + 0.02171 x − 0.03292 has the best fit (R² = 0.9655) to precision values in the desired range.

Bayesian Modeling of Cross-Reactive Compounds

Cross-reactive (“promiscuous”) compounds are frequently reported as hits in assays due to various biological and chemical artifacts rather than genuine biological activity.^23,24 While the “hit ratio” between the number of assays in which a compound was active and the number in which it was tested is a good measure of compound cross-reactivity, it has the drawback of not properly accounting for the number of measurements. Intuitively, 1/4 and 100/400 do not carry the same weight of evidence, even though they have the same hit ratio. We address this shortcoming with a Bayesian approach to cross-reactivity analysis.

Let θ be the probability that a particular compound will be a hit in a new assay, let N be the number of assays in which the compound was tested, and let n be the number of assays in which the compound was already called a hit, so that we have n/N → θ as N → ∞. Let θ₀ be the level at which compounds are qualitatively considered cross-reactive (based on experience we chose θ₀ = 0.25). Let us fix the number, N, of assays in which the compound was tested. For a given θ, we can assume that n has a binomial distribution,

P ( n | θ ) = ( N n ) θ n ( 1 − θ ) N − n .

However, in our case, we know values for n and N and we need to assess how likely it is that θ is larger than θ₀. We can write

P ( θ > θ 0 | n ) = ∫ θ 0 1 P ( θ | n ) d θ

and use Bayes’s formula to express P ( θ | n ) ,

P ( θ | n ) = P ( θ ) P ( n | θ ) P ( n ) .

A natural choice for the prior P ( θ ) (termed a conjugate prior) is the beta distribution ²⁵

B ( x | α , β ) = Γ ( α + β ) Γ ( α ) + Γ ( β ) x α − 1 ( 1 − x ) β − 1 .

It can be shown that combining binomial distribution equation (1) and prior beta distribution equation (4) with Bayes’s formula equation (2) produces a posterior probability P ( θ | n ) that has a beta distribution with parameters n + α and N – n + β.

We turned to historical screening data to determine appropriate values of the parameters α and β for the prior beta distribution. We selected all compounds that were tested in at least 20 assays, determined their hit ratios, and computed the average (0.126) and the standard deviation (0.0115) of the hit ratio distribution. Expressions for the mean and the variance of the beta distribution,

μ = α / ( α + β ) , σ 2 = α β ( α + β ) − 2 ( α + β + 1 ) − 1 ,

allow us to compute values for α and β. The resulting beta distribution accurately approximates the empirical hit ratio distribution (Suppl. Fig. S3). Once the parameters of the prior distribution are determined, we can use equation (2) to compute the probability that θ will be above the threshold θ₀ for a given n and N using the MATLAB function betainc(theta_0, n+alpha, N-n+beta, ‘upper’) ( Table 1 ).

OPEN IN VIEWER

Table 1.

Bayesian Promiscuity (Cross-Reactivity) Probabilities.

N,n	0	1	2	3	4	5	6	7	8	9	10	11	12
1	0.099	0.353
2	0.075	0.289	0.581
3	0.056	0.236	0.508	0.758
4	0.043	0.191	0.440	0.696	0.873
5	0.032	0.154	0.378	0.632	0.828	0.938
6	0.024	0.123	0.322	0.568	0.779	0.910	0.971
7	0.018	0.099	0.272	0.507	0.726	0.878	0.956	0.988
8	0.014	0.078	0.229	0.448	0.671	0.840	0.937	0.980	0.995
9	0.010	0.062	0.191	0.393	0.616	0.798	0.912	0.969	0.991	0.998
10	0.008	0.049	0.159	0.343	0.560	0.752	0.884	0.955	0.985	0.996	0.999
11	0.006	0.039	0.132	0.297	0.506	0.704	0.851	0.937	0.978	0.994	0.998	1.000
12	0.004	0.031	0.108	0.255	0.453	0.654	0.814	0.915	0.968	0.990	0.997	0.999	1.000
13	0.003	0.024	0.089	0.219	0.404	0.604	0.774	0.890	0.955	0.984	0.995	0.999	1.000
14	0.002	0.019	0.073	0.186	0.358	0.554	0.732	0.861	0.938	0.977	0.992	0.998	1.000
15	0.002	0.015	0.059	0.158	0.315	0.505	0.687	0.829	0.919	0.967	0.989	0.997	0.999
16	0.001	0.012	0.048	0.133	0.276	0.457	0.642	0.793	0.897	0.955	0.983	0.995	0.998
17	0.001	0.009	0.039	0.112	0.240	0.412	0.596	0.756	0.871	0.940	0.976	0.992	0.998
18	0.001	0.007	0.032	0.094	0.208	0.369	0.550	0.716	0.842	0.923	0.967	0.988	0.996
19	0.001	0.005	0.025	0.078	0.179	0.329	0.505	0.674	0.810	0.903	0.956	0.983	0.994
20	0.000	0.004	0.020	0.065	0.154	0.291	0.461	0.632	0.776	0.880	0.943	0.976	0.991

Probabilities P(θ > 0.25) for n ≤ 12 (columns) and N ≤ 20 (rows).

We also want to distinguish whether two compounds that connect to each other do so simply because they are each cross-reactive individually. Assume that two compounds with hit ratios θ and ρ were tested in N′ common assays and each was a hit in H′ of them; the overlap enrichment score, E, of the compound pair can be quantified as a ratio of actual and expected hit counts,

E = H ′ θ ρ N ′ .

Connected cross-reactive compounds tend to have lower overlap enrichment scores with values below 3, and connected selective compounds have overlap enrichment scores above 3.

Comparison of Bioactive Compounds

We selected a set of 2222 bioactive compounds from our HTS collection and explored whether historical screening data can be used to provide additional insights into their activities. For each pair of compounds, we computed an assay performance profile similarity score and filtered out pairs that did not achieve the significance threshold. We were left with 5996 connections among 934 compounds and used Cytoscape to visualize them ( Fig. 3 ). We turned to a method using random graphs with given expected degrees ²⁶ to detect compound communities—sets of compounds with relatively dense connections with other nodes in the community and relatively sparse connections to nodes outside the community. We used the algorithm of Farutin et al. ²⁷ that starts from each single connection and grows a community until no additional compound can improve its community score. By its nature, this algorithm produces overlapping communities and also determines community structures within larger communities. Whenever an overlap was sufficiently large, we manually merged two communities to reduce the resulting number of communities to consider. Finally, any unassigned compound that was connected to only one community was assigned to that community.

OPEN IN VIEWER

Figure 3.

Assay performance profiling distinguishes broadly and narrowly bioactive communities of compounds. Network graph showing connections between compounds (nodes), including cross-reactive compounds (red nodes; p ≈ 1) and non–cross-reactive compounds (blue nodes; p ≈ 0). Eight communities of similar assay performance (A–H) are discussed in the text; edges are colored by membership in these communities.

We used Bayesian modeling of cross-reactive compounds to examine the potential mechanisms of action for compounds within highly connected communities. In general, some bioactive compounds can be thought of as “broadly bioactive” ( Fig. 3 , community A, average overlap enrichment score 2.15, and, partially, communities E, F, and G, with average overlap enrichment scores between 2.8 and 3.15), meaning they score in a high fraction of assays and tend to connect to other such compounds. On the other hand, compounds that tend to score in fewer assays can be thought of as “narrowly bioactive” ( Fig. 3 , communities B–D and H, overlap enrichment scores above 4). We observed that some of these latter compounds segregate into communities based on their annotated mechanism of action.

Community A contains many compounds that are toxic at screening concentrations, including anticancer drugs. Community B corresponds to bioactive compounds that do not appear to be toxic and contains, for example, pain relievers, anti-inflammatory agents, and cardiovascular drugs. Evaluation of medium-sized communities revealed remarkable clustering of certain classes of small molecules. For example, community C is composed of antibiotics that work by a variety of mechanisms, including cephalosporins, fluoroquinolones, and tetracyclines. Furthermore, we observed a distinct subcommunity within C, containing mostly sulfonamide antibiotic compounds. Notably, almost all assays with significant hit-rate enrichment among members of community C were microbial assays. Community D consists entirely of steroids, such as hydrocortisone and dexamethasone. Community E contains a set of compounds that are narrowly bioactive but mildly promiscuous, in the sense that they exhibit relatively higher hit ratios, but not at the level that would classify them as promiscuous according to our Bayesian promiscuity model. Many compounds from this community scored highly in a collection of profiling assays using genetically related cell lines. Community F contains several central nervous system (CNS) drugs as well as antihistamines. Community G is enriched for azole antifungal compounds, which we also previously identified as inducers of OXPHOS transcription in muscle cells. ²⁸ We also observed clustering of close chemical analogues, or even independent instances of the same compound, in the smaller clusters. For example, community H has four compounds, all of which contain long hydrocarbon chains: three polyunsaturated fatty acids and farnesylthioacetic acid. These results demonstrate that historical screening data, even when focused on disparate areas of biology, can be used to reassemble known mechanisms of action and suggest new compound connections.

Target Identification

We applied our methodology to BRD7389, a small molecule we reported to induce insulin in murine pancreatic alpha cells. ²⁹ In that study, we compared the assay performance profile for BRD7389 with profiles of 9995 bioactive compounds across 32 assays. We used the list of the top 50 compounds for decision making, since the analysis was performed before we derived a threshold for z_jk. We observed multiple kinase inhibitors among bioactive compounds with high similarity to BRD7389, leading us to hypothesize that BRD7389 also targets protein kinases. We verified this hypothesis through biochemical profiling of a panel of 219 human kinases. In particular, BRD7389 potently inhibited the family of p90^Rsk kinases, and knockdown of these kinases increased insulin expression in the αTC cell line.

Identifying New Compounds

In the previous example, we compared a novel compound with a set of known bioactive compounds. We also applied this method to the reverse situation, in which we started with a few manually selected compounds with desired behavior and looked for novel compounds that exhibit similar assay performance profiles. We used two bioactive compounds: the mitochondrial uncoupler trifluorocarbonylcyanide phenylhydrazone (FCCP) and the natural product piperlongumine, which was recently shown to selectively target cancer cells by elevating reactive oxygen species (ROS) levels. ³⁰ In the case of FCCP, we had previously profiled a collection of bioactive molecules for mitochondrial activity in muscle cells ²⁸ and were interested in identifying additional mitochondrially active compounds. Assay performance profiling of FCCP revealed an expected anticorrelation with the mitochondrial complex I inhibitor rotenone. Accordingly, measurement of the mitochondrial membrane potential in U2-OS osteosarcoma cells using rotenone and the related uncoupler CCCP reflected these profiling results ( Fig. 4a ). Interestingly, piperlongumine was strongly anticorrelated with CCCP, also reflected experimentally. Although the nuclear factor kappa-B (NFκB) inhibitor parthenolide was connected to FCCP, it did not alter the mitochondrial membrane potential (data not shown).

OPEN IN VIEWER

Figure 4.

Assay performance profiling connects new compounds to compounds of interest. (a) Mitochondrial membrane potential change (ΔΨ_m) for compounds with high assay performance profile similarity to trifluorocarbonylcyanide phenylhydrazone (FCCP). U2-OS cells were treated for 90 min with either CCCP (10 µM) or the indicated compounds (20 µM). ΔΨ_m was measured with the fluorescent dye JC-1. (b) Proportion of cells positive for reactive oxygen species (ROS) after treatment with compounds similar to piperlongumine. U2-OS cells were treated for 2 h with the indicated compounds, and ROS levels were measured with the fluorescent dye CMH2-DCFDA (gray box: 1 standard deviation relative to DMSO treatment, NT). Andrographolide (1) and deoxysappanone B (2) were similar in assay performance profile to piperlongumine but had no effect in this assay. (c) Fluorescent micrographs (×4 magnification) of ROS-positive cells quantified in b.

Piperlongumine is an electrophilic compound reported to display selective toxicity to cancer cells relative to normal cells. Piperlongumine causes several phenotypes associated with oxidative stress, including ROS and protein glutathionylation. We sought to identify additional compounds with cellular activities similar to piperlongumine using assay performance profiling. To that end, we compared the assay performance profile of piperlongumine with other compounds in the Broad Institute screening collection to identify a set of similar compounds that was then tested for ROS levels in U2-OS cells. The connected compounds included two electrophilic compounds, parthenolide and BAY 11-7082, that were shown to be strong inducers of ROS levels ( Fig. 4b , c ), consistent with literature reports.^31,32 In this case, rotenone and CCCP, which also had similar profiles to piperlongumine, each increased ROS levels as well.

Screening Hit Prioritization

Assay performance profile comparisons can be used to prioritize HTS hits. The search for compounds chemically similar to HTS hits is a routine step in the analysis of HTS results. Assay performance profile similarity allows us to add an additional dimension by connecting compounds that potentially have the same or similar mode of action but may be different in chemical structure. That is, we anticipate that assay performance profile similarity will reveal connections not available through structure analysis alone. We computed pairwise assay performance profile similarities for a collection of HTS hits from a leukemic stem cell (LSC) HTS assay ³³ and used them as a basis for hierarchical clustering ( Fig. 5a ). We further combined biological and chemical similarities to build a network of interactions ( Fig. 5b ) that revealed one large and four small clusters of related compounds.

OPEN IN VIEWER

Figure 5.

Assay performance profiling aids in prioritization of new high-throughput screening (HTS) hits. (a) Hierarchical clustering of 41 hits from a leukemic stem cell (LSC) screen ³³ and pairwise assay performance profile similarities among the hits. (b) Network of connections among compounds identified as LSC HTS hits (nodes) using chemical similarity (red edges; T_C > 0.25 using ECFP_6 fingerprints), assay performance profile similarity (green edges; z_jk ≥ 10, and red dotted line in dendrogram), or both (black edges). Edge widths are proportional to the number of assays two compounds have in common (range, 14–300 assays).