Automated structure–activity relationship mining: connecting chemical structure to biological profiles

Compound sets

We assembled three distinct compound sets with varying structural properties and levels of biological activity annotation. We first selected a collection of 19,637 structurally diverse compounds derived from diversity-oriented synthesis (DOS).²⁰ The library was synthesized by the Center for the Science of Therapeutics at the Broad Institute using a build/couple/pair strategy.²¹ The compounds were built around 23 novel chiral core structures by systematically varying the configuration of core stereocenters and decorating them with various side-chains (see Supplementary Information). Whenever possible, complete sets of stereoisomers were included in our set. No biological activity information was used in the selection of the DOS set. Second, we used PubChem (http://pubchem.ncbi.nlm.nih.gov/) assay data to select 10,162 compounds from the Molecular Libraries Small Molecule Repository (MLSMR) that cover many different chemical structures and biological activities. Third, 2,222 known drugs and probes were added as bioactivity landmarks. Again, these compounds were selected to be structurally diverse and cover a wide range of biological activities. The complete compound set is available as Supplementary Material at http://jbx.sagepub.com/.

Gene-expression profiles (GE)

We followed the protocol of Peck et al. (see Supplementary Methods for a detailed protocol).¹⁰ Briefly, 3,500 U-2 OS cells (ATCC, cat. no. HTB-96) were seeded per well in 384-well plates. After 24h of incubation at 37°C, compounds were added and the cells incubated for another 6h of treatment. Cell lysates were prepared and transferred to oligo-dT capture plates to immobilize mRNAs. cDNA was then generated via reverse transcription and amplified via polymerase chain reaction (PCR). The cDNA was first annealed with a mix of specific probe pairs designed for each of 978 transcripts. All upstream probes contain a universal 20-nucleotide (nt) primer site (complementary to T7 primer), a specific 24-nt barcode sequence, and a transcript-specific sequence of 24nt complementary to the 3′-end of the transcript. Downstream probes consist of a 5′-phosphorylated 20-nt sequence that is complementary to the sequence contiguous to the binding site of the upstream probe on each transcript, as well as a 20-nt universal primer (T3) site. After annealing and removal of free probes, adjacent probes bound to transcript cDNAs were joined by ligation and amplified by PCR using T3 and 5′-biotinylated T7 primers. The resulting amplicons were hybridized with color-coded Luminex microspheres each of which carries capture probes complementary to one of the barcode sequences in the amplicons. After incubation with streptavidin-phycoerythrin, the captured amplicons were quantified by flow cytometry measuring phycoerythrin fluorescence. At the same time, transcript identity was identified by bead color. Compound effects were measured in triplicate treated wells.

Multiplexed cytological morphology profiles (MC)

We followed the protocol of Gustafsdottir et al. (see Supplementary Methods for a detailed protocol).¹¹ Briefly, 1,500–2,000 U-2 OS cells (ATCC, cat. no. HTB-96) were seeded per well in 384-well clear-bottom imaging plates. After 24h of incubation at 37°C, compounds were added and cells incubated for another 48h of treatment. Six fluorescent dyes were added to stain different cell compartments: nucleus (Hoechst 33342), endoplasmic reticulum (concanavalin A/AlexaFluor488 conjugate), nucleoli (SYTO 14 green fluorescent nucleic acid stain), Golgi apparatus and plasma membrane (wheat germ agglutinin/AlexaFluor594 conjugate, WGA), F-actin (phalloidin/AlexaFluor594 conjugate) and mitochondria (MitoTracker Deep Red). WGA and Mitotracker stains were performed on live cells. Afterwards, cells were fixed with 16% paraformaldehyde and the remaining dyes added. Images in 5 fluorescent channels were captured for 9 sites per well at 20× magnification. Images were measured in quadruplicate and CellProfiler software was used to calculate morphological features for each cell.¹¹

Chemical descriptors

Biological descriptors

Frequent pattern mining and association rules

Frequent itemsets were determined using a publicly available implementation (http://mahout.apache.org/) of the Parallel FPgrowth algorithm²³ on a Hadoop platform (version 1.0.4; http://hadoop.apache.org/). Association rules were generated from the list of frequent patterns in R (version using the arules package.²⁴

SAR score

We calculated confidence, purity, and p-values (Fisher’s exact test) for each rule and combined them into a single SAR score. P-values were transformed into negative decadic logarithms. The distributions of all three quality measures were first normalized using a Box-Cox transformation²⁵ to avoid one criterion dominating the others and the resulting values for each quality measure were scaled to a range [0,1]. Purity values were multiplied with a weighting factor of 0.2. The sum of the normalized quality measures yielded the SAR score for each compound. SAR scores were scaled to range [0,1] by division with the theoretical maximum score of 2.2.

Replicate-based hit selection

To determine which compounds led to stable profile changes compared to negative controls (DMSO), we used replicate correlation and connectivity, i.e., the similarity of replicates compared to all other wells on the same set of replicate assay plates. Replicate correlation was calculated as the Pearson correlation coefficient between their GE (or MC) profiles. Replicate connectivity was calculated for each pair of replicates (e.g., R1 and R2) by ranking all wells on the same set of replicate plates by their profile similarity against R1 and identifying the fraction of unrelated wells that rank higher than R2. The reverse calculation (R2 vs. R1) was also performed because connectivity values are not symmetric. To allow individual replicates to fail, we only considered the top 50% of replicate correlation and connectivity values for each compound. We then calculated the two-dimensional distribution of negative-control replicate correlation and connectivity. Compounds were considered hits if they exceeded the 97% confidence interval of this negative-control distribution.

Background

We applied FPM and ARM to small-molecule gene-expression and cell-morphology profiles to derive SAR rules of the form {chemical fragments} → {biological effect pattern}, indicating that compounds containing specific molecular fragments (or chemical substructures) often lead to the same pattern in biological responses (e.g. similar expression changes for a common set of genes). Even though our method can be applied to more abstract chemical attributes (e.g. calculated properties), using fragments makes it easy to interpret the final rules. Rules that connect fragments to biological effects can directly inform synthetic chemistry decisions, especially if these fragments match existing molecular building blocks.

We derived biological attributes from two profiling assays, gene expression (GE) and cell morphology (multiplexed cytological, MC). For GE profiles, the expression levels of 978 protein-coding transcripts were determined by ligation-mediated amplification and Luminex bead-based detection.¹⁰ For MC profiles, changes in 812 cell morphology features were captured by automated microscopy and computational image analysis. We used a “cell-painting” assay with 6 fluorescent dyes in five channels to distinguish cellular compartments and organelles.¹¹

In the following, we follow the notation of Tan, Steinbach, and Kumar¹⁶ to introduce FPM and ARM. FPM and ARM operate on collections of binary attribute vectors (“objects”). Let $I=\left\{i_1,i_2,\dots ,i_d\right\}$ be a set of binary attributes and $T=\left\{t_1,t_2,\dots ,t_N\right\}$ a set of objects where each object is characterized by a subset of all attributes I, i.e., $t_i\subseteq I$. The support count σ of a set of attributes X is defined as the number of objects that have all attributes in X:

Given a minimum support count threshold s, FPM finds all sets of attributes (“itemsets” or “patterns”) that occur at least s times (“frequent patterns”).

For the purpose of SAR rule mining, objects t_i represent compounds tested in a profiling experiment. I then is a set of chemical and biological attributes. Each compound t_i is represented by a specific subset of chemical and biological attributes – the substructures it contains and the biological response parameters it affects (Figure 1). A screening library usually contains only a few compounds with both similar structure and similar biological activity. We therefore used a support count threshold of $s=3$, identifying all patterns that occur at least three times.

From these patterns, ARM derives rules of the form $X\to Y$, indicating that the occurrence of attribute set X (left-hand side, LHS) is associated with the occurrence of attribute set Y (right-hand side, RHS). We limited chemical attributes to the LHS and biological attributes to the RHS to obtain SAR rules of the form {chemical fragments} → {biological pattern}. To reduce redundancy, we only consider rules that consist of closed patterns, i.e., sets of attributes for which no superset exists that has the same support count.

Three measures of rule quality or “interestingness” are then used to filter and prioritize rules. First, the confidence of a rule is defined as the fraction of objects with LHS attributes X for which the rule is true. In other words, the confidence quantifies how likely it is to observe the biological pattern Y when the chemical attributes X are given:

Second, we defined the purity of a rule as the reverse measure, i.e., the fraction of RHS attributes Y for which the rule is true, which gives a measure of how characteristic a biological pattern is for a certain chemotype:

Third, we calculated a p-value for each rule using Fisher’s exact test.²⁶ Low p-values indicate that chemical attributes X and biological pattern Y are unlikely to be independent of each other.

We only kept rules with p-values smaller than 0.01. To rank the remaining rules by their quality, we combined confidence, purity, and p-value into a single SAR score ranging from 0 to 1 (see Methods). A weighting factor can be applied to each individual quality measure to balance their effect on the score. We lowered the influence of the purity with a weighting factor of 0.2 because we considered a high purity less critical for a good rule. This avoids exaggerated down-rating of cases where multiple different structural classes share a biological activity pattern. However, we considered rare and chemotype-specific bioactivities more interesting than common ones and therefore included purity in the score. The set of rules usually contains groups of related rules. We created a “minimal rule set” by removing all rules that consist of attribute subsets or supersets of a rule with a higher SAR score. Unless stated otherwise, all reported rule ranks refer to this minimal set.

We tested our method on GE and MC profiles for two distinct compound libraries. In the following sections, we discuss examples of rules obtained with different chemical structural and biological attributes.

ARM with custom substructure descriptors on planned synthetic libraries leads to intuitive and testable SAR hypotheses

SAR analyses take advantage of structural similarity between compounds. While analog series in advanced projects are tailored to meet this requirement, many screening libraries do not. Early-stage SAR exploration is hence especially useful for planned screening libraries that were designed to include groups of compounds with defined structural relationships.⁴

We therefore included 19,637 compounds in our experiment that were derived from diversity-oriented synthesis (DOS).²⁰ The compounds represent systematic combinations of three basic structural elements: core structures, appendages, and stereochemistry. Compounds were built by modular synthesis, attaching different appendages to 23 chiral core structures and systematically varying the configuration of core stereocenters.²¹ The resulting collection is chemically diverse but, at the same time, contains related structures: groups of compounds that share common cores, stereochemistry, or appendages can be directly compared in SAR analysis. Importantly, this modular design also allows quickly synthesizing structural analogs of identified hits for follow-up studies.

We used as chemical attributes a custom structural fingerprint (“DOS-FP”) that reflects the modular design of DOS compounds to obtain rules that directly relate to compound synthesis. Based on the three DOS diversity elements, DOS-FP specifies for each compound the (a) core structure, (b) appendages attached to the core (as a simple occurrence or as “R-groups” that include information about their attachment position in the core), and (c) configuration of core stereocenters. Stereoisomers of appendages were treated as distinct structures. Note that R-group and stereochemistry features refer to specific positions in the core and can hence only occur in combination with a core feature.

We added a reference set of 2,222 compounds with known bioactivities (BIO) to the DOS library to support the generation of hypotheses about the biological effects of novel DOS compounds. The BIO collection was assembled to contain structurally diverse compounds that cover a wide range of biological activities and often have known direct targets. We measured both GE and MC profiles for most compounds. Exact compound numbers measured with each method differ slightly due to quality-control filters on the experimental data (GE: 17,553 DOS compounds + 1,935 known bioactives; MC: 17,805 DOS compounds + 2,211 known bioactives).

We first aimed to identify SARs on global phenotypic effects. We hierarchically clustered GE and MC profiles to find groups of compounds with similar biological activities. Each node of the resulting dendrogram was considered one cluster. Each compound could therefore be a member of multiple clusters (characterized by superset/subset relationships). Importantly, both DOS and BIO compounds were clustered, but only DOS compounds were subject to ARM for SAR rule mining. While BIO compounds cannot be represented with DOS-FP, this allowed us to evaluate structural rules for DOS compounds in the context of biologically annotated reference molecules.

We found a total of 6,206 rules for GE and 7,098 rules for MC (2,752 and 2,861 for minimal rule sets, respectively). About 70% of the rules for both GE and MC were combinations of core and stereo features and another 17% specified core features alone (Table 1). This distribution is likely to reflect the fact that the DOS library contains many more compounds that match on these two features than on appendages. However, rules that include R-group information are about 5 to 6-fold enriched among the top 10% of rules (Table 1), indicating their importance for describing specific small-molecule effects.

The majority of rules that include R-group information specify analog series, i.e., compounds that vary in one R-group (70% for GE and 80% for MC). Hence, these rules lead to directly testable hypothesis that are easily accessible through DOS. Analogously, 11% (GE) and 14% (MC) of rules with R-group information specify all R-groups and hence apply to series of stereoisomers. The remaining fraction consists of rules that allow two of three R-groups to vary. Overall, chemical feature combinations occur with highly similar frequencies in the rules for both profiling experiments (Table 1). These frequencies did not change appreciably when only profiles with very good replicate agreement were used (see Methods).

We next examined rules for DOS compounds that co-cluster with known bioactive molecules to find rules with interpretable structural and biological parts. The highest-ranking rule identified a cluster that was strongly enriched for compounds carrying an ortho-aminoanilide (OAA) appendage (Figure 2). 50% of all OAA DOS compounds are contained in this cluster (confidence: 0.5) and hence have similar GE profiles. The OAA motif is known as a biasing element for inhibitory activity against class I histone deacetylases (HDACs) due to its ability to interact with zinc and active site resides chains required for hydrolysis of the acetyl moiety on lysine side chains.²⁷ Indeed, the DOS compounds co-cluster with two known HDAC inhibitors that also contain an OAA residue, mocetinostat (MGCD0103) and tacedinaline (CI-994). Furthermore, the only DOS cluster member that does not match the structural rule, BRD4805, has recently been identified as a low micromolar HDAC inhibitor.²⁸

Mocetinostat, tacedinaline, and BRD4805 exhibit selectivity for HDAC1, 2, and 3 over other isoforms.^28,29 Interestingly, less selective inhibitors²⁹ contained in our compound library – trichostatin A (TSA), vorinostat (SAHA), and belinostat (PXD-101) – formed a separate cluster, indicating that GE profiling can distinguish compounds with distinct HDAC isoform selectivity patterns. We therefore hypothesize that the clustered OAA DOS compounds show some degree of specificity for HDACs 1–3. Supporting this hypothesis, stereoisomers of structures 1, 2, and 4 have recently been identified as HDAC1-3 selective inhibitors.³⁰ Members of this structural class shared the highest profile similarity with the known HDAC inhibitors in our experiment (Table S1).

Several of the OAA DOS compounds have further been tested in abiochemical assay measuring inhibition of HDAC2 deacetylase activity (Figure S2). 90% (9 of 10) of the OAA DOS compounds that were in the cluster were hits in this screen. By contrast, OAA DOS compounds that were not part of the profile-based cluster achieved a lower hit rate (57%, 4 of 7). Beyond this enrichment of hits in the cluster, we found no apparent correlation between the profile similarities of DOS compounds to known HDAC inhibitors and their performance in the HDAC2 assay. While HDAC2 serves as a representative class I HDAC, the discrepancies between profile similarity to class I HDAC inhibitors and degree of HDAC2 inhibition may reflect the fact that ortho-aminoanilde-based HDAC inhibitors also target HDAC1 and HDAC3 and/or other in vivo cellular factors that may be influencing the gene-expression signatures.

Taking SAR information into account can increase our confidence in profiling assay results. Only 6 out of the 14 OAA DOS compounds in the cluster had sufficient replicate agreement to be considered ‘hits’ in the GE assay, i.e. compounds that induce a reproducible gene-expression change (see Methods). Relying on replicate agreement as an activity criterion alone would hence disregard the remaining 8 compounds, and thus 55% of the HDAC2 hits in the cluster (5 of 9). This result provides an example of the overall utility of our method for guiding the selection of initial hits for follow-up biological studies and chemical synthesis.

The OAA rule is an example of a single structural motif whose presence is sufficient to elicit a relatively defined bioactivity. Although in principle one could have predicted such compounds to have HDAC inhbitory activity by structure alone, these results help validate the method for prediction of less well-known features. More commonly, combinations of features determine a compound’s biological effect. One such example is a set of related DOS compounds that co-cluster with known microtubule destabilizers (Figure 3). The corresponding rule specifies the core, one of two R-groups, and, importantly, all stereocenters for this compound class, effectively describing a series of analog structures. Stereochemistry is a major determinant of specificity for this compound series. Omitting stereochemistry from the rule causes the confidence to drop from 0.5 (3 rule matches out of 6 structure matches) to 0.065 (3 out of 46). The validity of this rule is further supported by the fact that it occurs at the top of both GE (rank 6) and MC (rank 12) rule lists. Microtubule inhibitors co-cluster with the DOS compounds in both the GE and MC rule (Figure 3). Unlike the OAA case, the known bioactive compounds are not structurally related to the DOS compounds, making this an interesting rule to study for lead-hopping purposes. Using structural features that directly match physical DOS building blocks allows rapidly testing such hypotheses as analogs and stereoisomers can be synthesized quickly (if they do not already exist).

General structure descriptors allow rule generation for arbitrary compound collections

Not all libraries are designed to include compounds with defined structural relationships like the DOS collection. However, arbitrary libraries can be mined by our method by replacing the chemical structural descriptors. We used this approach to analyze MC profile data for a collection of 10,226 compounds from the Molecular Libraries Small Molecule Repository (MLSMR). We used extended connectivity fingerprints (ECFP4²²) as a generic structure descriptor that can be used with any compound collection. While any binary fingerprint can readily be used with our method, we chose ECFPs as a reference for this study because they are widely applied, allow direct mapping of features to structural fragments, and were specifically designed for SAR studies.²² They achieve high structural resolution by encoding a very large number of features and thus present a good test for the ability of ARM to mine large feature spaces.

We used clusters based on MC profiles as biological attributes. Due to the versatility of ECFP4, we were able to include our set of known bioactive compounds in the rule-mining step. We found 14,198 rules (7,328 after redundancy filtering) of which 2,176 contained known bioactive compounds. Unlike the DOS-FP features, ECFP4 substructures can describe arbitrary parts of the molecule through combinations of overlapping substructures (Figure 4). This description allows a more fine-grained and unbiased mapping of relevant substructures. The downside is a reduction of interpretability as rules cannot be directly mapped to existing analogs or straightforward synthesis strategies. In addition, a large number of features are usually needed to achieve high resolution and generality. Even though ARM can handle such large data sets, the number of rules can increase as a result and complicate downstream analyses. Therefore, each compound collection and project requires balancing generality, complexity, and interpretability when choosing a chemical descriptor.

Direct mining of attributes identifies biological signatures for sets of related compounds

Clustering biological profiles captures global similarities and differences between the biological effects of compounds. However, information on individual features (e.g., genes from a GE profile) is lost in this process. Since ARM is designed to operate on high-dimensional data, it can be used to mine biological features without intermediate aggregation. Useful applications include mining compound-target profiles or GE signatures, i.e., sets of genes that change expression levels in response to compound treatment.

We generated a signature of up- and down-regulated genes for each compound by applying a z-score threshold of 2 and −2 to GE profiles, respectively (relative to DMSO-treated controls). The signatures were combined with DOS-FP features and subjected to ARM. The top-scoring rule is again the occurrence of an OAA residue already found with GE clusters as biological attributes (confidence=0.3, purity=1, p=1.68×10⁻²⁶, SAR score=0.98). The associated GE signature consists of only three genes: CASP3:down, IER3:up, TERF2IP:down. This relatively low number of common features leads to important differences in specificity compared to the analogous clustering-based rule. In addition to mocetinostat and tacedinaline, compounds that do not target HDACs are also found to regulate these genes (15D-prostaglandin J2, withaferin A, thioridazine, and aminopurvalanol A). However, more specific rules can be found that achieve only slightly lower SAR scores due to their reduced confidence. For example, the signature CASP3:down, IER3:up, TBPL1:up, TERF2IP:down is only found in benzamide-carrying DOS compounds and mocetinostat and tacedinaline (SAR score 0.95, rank 4 in the full set of rules).

The second-ranked rule is highly specific for a small set of compounds with a signature of 13 up-regulated genes (Figure 5). A functional analysis strongly related the 13-member gene set to sterol and cholesterol synthesis (Table S2). Interestingly, sirolimus and fulvestrant are also found to change the expression of these genes. Irrespective of their distinct mechanisms of action, each has been linked to cholesterol levels and metabolism.^31,32 A corresponding rule found by clustering full-length profiles as described above achieved high ranks for both GE (rank 8) and MC (rank 11). However, neither of the known bioactive compounds co-clusters with the DOS compounds in both assays. This finding supports the view that gene-signatures capture subsets of the global biological effects of perturbagens.³³ We therefore hypothesize that the DOS compounds identified by this rule have related biological effects connected to cholesterol metabolism. This example illustrates how the combination of profiling and automated SAR analysis can aid data mining and quickly lead to testable hypothesis.