ERpred: a web server for the prediction of subtype-specific estrogen receptor antagonists

Data compilation and curation

Two data sets of inhibitors against ERα and ERβ were compiled from the ChEMBL database, version 25 (Gaulton et al., 2017), which consisted of an initial set of 11,618 compounds for ERα and 7,810 compounds for ERβ. The IC₅₀ was selected as the bioactivity unit for further investigation and after the data curation process, this led to a final data set of 1,593 and 1,281 compounds for ERα and ERβ, respectively. As this study sets out to develop a classification model for both ERα and β we therefore, defined thresholds of <1 and >10 µM (corresponding to pIC₅₀ values of 6 and 5, respectively) for distinguishing actives from inactives, respectively. Moreover, the intermediate biological activity with IC₅₀ values ranging between 1 and 10 µM were not selected for this study. A final set of non-redundant and curated compounds consisting of 1,194 and 997 inhibitors were obtained for ERα and ERβ, respectively.

Molecular descriptors

Fingerprint descriptors for compounds in the data sets were computed using the PaDEL-descriptor software (Yap, 2011). SMILES notation was used for the calculation of molecular descriptors. Structures were pre-processed so as to remove salt and standardize tautomers using the built-in function of the PADEL-descriptor software. In general, molecular descriptors are important for QSAR studies as they characterize molecular properties and chemical structure information in quantitative or qualitative forms. As previously described in Malik et al. (2020) 12 molecular fingerprints belonging to 9 classes consisting of AtomPairs 2D, CDK fingerprinter, CDK extended, CDK graph only, E-state, Klekota–Roth, MACCS, PubChem and Substructure were used for describing the chemical structures. Furthermore, Substructure, Klekota–Roth fingerprint and 2D atom pairs consisted of two versions: (1) 1 or 0 denotes the presence or absence of the descriptor and the (2) representing the frequency value of the descriptor with a count version.

Data filtering

In order to remove inherent complexity and bias that may be introduced to the model building process, constant and near constant variables were removed to select the fingerprint descriptor sets. Particularly, near constants were identified using a standard deviation (SD) threshold of 0.1 whereby variables with SD values less than 0.1 were selected for further analysis.

Data splitting

The Kennard–Stone algorithm (Kennard & Stone, 1969) was applied for splitting the data into an 80/20 split where, the internal set comprised of 80% of the entire data set while the external set consisted of the remaining 20%. The internal set is used as the training set and also subjected to 5-fold CV. The external set is used as the testing set whereby the trained model will be applied to this data set to make predictions and thus determine the model’s robustness. For the internal set of ER alpha (comprising 528 compounds), active compounds accounted for 53.59% while inactive compounds are comprised of 46.41%. Similarly, for the external set (comprising of 131 compounds), active and inactive compounds were split into sets of 53.44% and 46.56%, respectively. For ERβ, the percentage of active and inactive in the internal data set (composed of 572 compounds) was 78.15% and 21.85%, respectively. Similarly, in the external data set (composed of 142 compounds), the percentage for active and inactive was 78.12% and 21.88%, respectively.

Statistical analysis

As described in our previous work (Malik et al., 2020), trends in individual descriptors of active and inactive compounds were determined through 6 common descriptive statistical parameters, encompassing the minimum (Min), first quartile (Q1), median, mean, third quartile (Q3) and maximum (Max) parameters. The outcome was visualized in the form of a box plot using ggplot2, a package of the R program. In addition, the Mann–Whitney U test (also known as the Wilcoxon Rank Sum test) was conducted to determine the statistical significance in terms of the p-value.

Multivariate analysis

Random forest (RF) is used for the building of classification models owing to its robust model performance and interpretability. RF was successfully used in our recent work for modeling the bioactivity of hepatitis C virus inhibitors (Malik et al., 2020). In essence, RF is an ensemble classifier that employs an N number of decision trees (specified by the ntree parameter) to learn the inherent patterns from the input data (Breiman, 2001; Breiman et al., 1984). In this study, a five-fold cross-validation (5-fold CV) procedure was applied for tuning the ntree parameter (100, 1,000, 100) and the mtry parameter (5, 30, 5) via the use of the tuneRF function from the randomForest package (Liaw & Wiener, 2002). In order to provide a better understanding of the biochemical activity of the inhibitors, feature selection was estimated using the built-in importance estimator of the RF model. The mean decrease of the Gini index (MDGI) was utilized to estimate the important descriptors (Weidlich & Filippov, 2016). Descriptors affording the largest value of MDGI represents the most important features as that descriptor contributes most significantly to the model performance.

Model validation

Parameters commonly used for evaluating the model performance of binary classification problems are typically based on true positives (TP), true negatives (TN), false positives (FP) and false negatives (FN). Particularly, the fitness of the model was assessed using various statistical parameters including the overall prediction accuracy (Ac), sensitivity (Sn), specificity (Sp) and Matthews correlation coefficient (MCC) (Song & Tang, 2004).

(1) $$Ac=TP+TN/\left(TP+TN+FP+FN\right)\times 100$$

(2) $$Sn=TP/\left(TP+FN\right)\times 100$$

(3) $$Sp=TN/\left(TN+FP\right)\times 100$$

(4) $$MCC=\left(TP\times TN\right)-\left(FP\times FN\right)/\sqrt{\left(TP+FP\right)\left(TP+FN\right)\left(TN+FP\right)\left(TN+FN\right)}$$where TP, TN, FP and FN represent the instances of true positives, true negatives, false positives and false negatives, respectively.

Applicability domain analysis

The main purpose of the applicability domain (AD) is to estimate the boundaries within which the model can make reliable and accurate predictions for compounds on the basis of similarity with the compounds on which the model was constructed. The compounds that satisfy the scope of the model are within the AD while the rest are outside the AD. In this study, we used the principal component analysis (PCA) bounding box to assess the AD of compounds from the training (internal) and testing (external) sets.

Web server development

The predictive model was exported as an RDS file (i.e., the model.RDS file) and subsequently deployed as a web server. Particularly, the web server is coded in R via the use of the Shiny R package (i.e., a web framework for the R environment). Technically, the web server is comprised of two major components: (1) user interface and (2) server, which are saved as ui.R and server.R, respectively. The ui.R file accepts input values (i.e. the SMILES notation of query compounds) and transfers this information to the server.R file where the SMILES notation is submitted to the PaDEL-Descriptor software. After descriptor calculation, the computed descriptors are used as input to the predictive model (i.e., the model is exported as model.RDS file) which will then classify query compounds as either being active or inactive (i.e., the bioactivity class label). Such predicted class labels are then printed out onto the web server whereby users can also download the predicted results as a CSV file.

Reproducible research

The data and code used in this study are publicly available on GitHub at https://github.com/chaninlab/ERpred/.

Chemical space analysis

Chemical space analysis is employed to explore the characteristic differences between the active and inactive compounds. The general chemical space was first visualized as a function of the molecular weight (MW) vs the Ghose–Crippen–Viswanadhan octanol-water partition coefficient (ALogP). In addition, the active and inactive compounds were further compared using the Lipinski’s rule-of-five (Ro5) descriptors. Briefly, the Ro5 describes the drug likeness of compounds on the basis of their molecular properties namely molecular weight (<500), octanol–water partition coefficient (ALogP < 5), the number of hydrogen bond acceptors (<10) and the number of hydrogen bond donors (<5) (Lipinski et al., 2001). Visualization of the MW chemical space as a function of ALogP is shown in Fig. 3. As can be observed for both ERα and ERβ, most of the compounds are clustered within the MW range of 200–500 Da with an ALogP in the range of 1 and 6. In addition, Fig. 4 shows the distribution of active and inactive compounds according to the Ro5 descriptors. It is observed that both ER subtypes contain compounds following the Ro5 criteria such as, MW of less than 500 Da, ALogP value of less than 5 and nHBDon and nHBAcc values of less than 10. It can also be seen that for ERα, some of the active compounds have an ALogP of greater than 5, but the number is very minimal. Furthermore, the results from statistical analysis displays a significant difference between the active and inactive compounds using the Mann–Whitney U test. Most of the active compounds (422.68 ± 91.52) were larger (i.e., higher MW) than the inactive compounds (350.35 ± 79.82), which was observed from the mean values of box plots. Similarly, the ALogP values of the active compounds (4.36 ± 1.37) were greater than the inactive compounds (3.17 ± 1.53). However, it was observed that both active and inactive compounds had similar nHBDon values while the active compounds had nHBAcc values that were lower than the inactive compounds. On the other hand, for ERβ, the MW between the active (356.94 ± 92.43) and inactive compounds (351.69 ± 94.80) was not statistically significant as determined using the Mann–Whitney U test. Nonetheless, the ALogP was very statistically significant with the active group (3.82 ± 1.6) displaying higher values than the inactive group (2.91 ± 1.5). Similar to the ERα subtype, the nHBDon values of both the active and inactive groups were on par while the nHBAcc for the active compounds was seen to be a lot lower than the inactive compounds.

Moreover, PubChem fingerprints were utilized for estimating the AD of the CSAR model developed herein which were further used as input values for PCA analysis. The resulting PCA scores plot can be seen in Fig. 5. The data set for ERα and ERβ comprised of 1,194 and 997 compounds, respectively were further divided into internal (80%) and external sets (20%) using the Kennard–Stone algorithm. Of note, the training set is composed of the internal set which is utilized to build the model and thus make predictions on the external set. In addition, both the internal and the external sets are also subjected to 5-fold CV. Furthermore, the chemical space distribution as observed from Fig. 5 shows that the external set (i.e., the testing set, represented using blue dots) lies within the boundaries of the internal set (i.e., the training set, represented using red dots). Thus, the AD is well defined for the CSAR model developed herein, as shown through these results.

In order to develop a deeper understanding of the ER chemical space, the active and inactive datasets for both ERα and ERβ were sorted according to their pIC₅₀ values. The top 10 active and bottom 10 inactive compounds were taken from each set and applied to the Scaffold Hunter software for further analysis (Schäfer et al., 2017). Particularly, the major scaffolds identified were then observed in further detail in terms of the number of nodes in each parent scaffold lineages and were also subjected to extensive analysis as will be discussed in paragraphs hereafter. Fig. 6 shows the schematic process of the scaffold analysis.

As mentioned above, further in-depth exploration of ERα and ERβ inhibitors led to the identification of top actives and inactives as shown in Figs. S1 and S2, respectively. The top active compounds of ERα had a bioactivity range of 8.398–8.698 pIC₅₀ while the top active compounds of ERβ had a bioactivity range of 8.045–8.522 pIC₅₀. Likewise, the bioactivity range for the top inactive compounds pertaining to the ERα and ERβ groups were observed to have pIC₅₀ in the range of 3.008–3.252 and 3.000–3.587 pIC₅₀, respectively.

In addition, a rigorous analysis of the chemical space was conducted by investigating the underlying scaffold structures as presented in actives and inactives using the Scaffold Hunter software (Schäfer et al., 2017). Particularly, chemical structure clouds were created for actives and inactives of ERα and ERβ as shown in Fig. 7. Analysis of ERα actives led to the identification of two major scaffolds namely 2,3-dihydro-1,4-benzoxathiine (with a frequency of 74) and diaryltetrahydronaphthalene (with a frequency of 27). On the other hand, ERα inactives consisted of five major scaffolds namely 4-hydroxypyrimidine (with a frequency of 26), 3-phenyl-4-(3H)quinazolinone (having 4-hydroxypyrimidine as the parent with a frequency of 24), chromane (with a frequency of two), pyran-4-one (with a frequency of 8) and leucoline (with a frequency of 17). Similarly, analysis of ERβ actives led to the identification of five top scaffolds whereby one is in common with ERα actives (i.e, lasofoxifene having a frequency of three). The other 4 scaffolds identified included 2-phenyl-1H-inden-1-one (with a frequency of nine), 2-phenylnaphthalene (with a frequency of 59), 2-phenylbenzofuran (with a frequency of 21), and 1,2,9,9a-tetrahydrofluoren-3-one (with a frequency of 25). Furthermore, three major scaffolds were observed for ERβ inactives, which comprises of coumarone, pyrazolopyrimidine combined with a furan and 1,4-dioxa-8-azaspiro[4.5]dec-8-yl(2-furyl)methanone. The importance of these above-mentioned scaffolds will be discussed in the Structural interpretation section below.

QSAR modeling

This study follows the Organisation for Economic Co-operation and Development (OECD) (OECD, 2014) guidelines for the development of robust QSAR models. These guidelines are applied in all our work as previously mentioned (Malik et al., 2020) and comprises of the following main points: (i) the data set has a defined endpoint, (ii) uses an unambiguous learning algorithm, (iii) the applicability domain of the QSAR model is well defined, (iv) appropriate measures of goodness-of-fit, robustness and predictivity and (v) mechanistic interpretation of the QSAR model. Following these aforementioned guidelines to develop interpretable QSAR models, this study makes use of molecular fingerprints that are interpretable which are computed using the PaDEL-Descriptor software. As shown in Table 1, three out of the 12 fingerprints (i.e., PubChem, Substructure and Klekota–Roth) are readily interpretable. In addition, Table 2 provides details on the model performances for all 12 fingerprints.

In this study, we developed a QSAR model based on the random forest algorithm in order to differentiate the active and inactive inhibitors for ERα and ERβ subtypes. Table 2 shows the results from the RF model with 12 different types of fingerprints over an internal validation test, CV and an external validation test. The best averaged values were observed as Ac of 94.65% and 92.25% and MCC of 0.89 and 0.76 for ERα and ERβ, respectively which was achieved for the PubChem fingerprint descriptors as evaluated by CV. Concurrently, Klekota–Roth and Substructure descriptors also performed well harbouring the second and third highest averaged values for Ac and MCC in which models built using Klekota–Roth fingerprints afforded Ac and MCC values of 90.83% and 0.81%, respectively for ERα and Ac and MCC values of 94.36% and 0.82%, respectively for ERβ. Similarly, models built using the Substructure fingerprints afforded Ac and MCC values of 93.89% and 0.87%, respectively for ERα and Ac and MCC values of 94.36% and 0.82%, respectively for ERβ. Although the Ac and MCC values of models built using the PubChem fingerprints for ERβ were not superior to the models built using fingerprints of Klekota–Roth and Substructure, they are quite comparable. In addition, taking into account the Ac and MCC values as well as the overall external and CV, for both ERα and ERβ and the interpretability of the features, we determined that the PubChem fingerprints were the ideal choice for interpretation of the model.

Mechanistic interpretation of feature importance

In order to provide a better understanding of the mechanistic details governing ERα and β subtypes, an analysis of the feature importance on selected informative descriptors was conducted. Owing to the inbuilt ability of feature importance estimation of the RF model and its great prediction performance, this study utilized it for analysis. Generally, two measures are used to rank the important features, namely the mean decrease of the Gini index and the mean decrease of the accuracy. As reported by Calle & Urrea (2011) that the Gini index had more robust results compared to those from the accuracy, we utilized the mean decrease of the Gini index to rank the importance of the PubChem feature descriptors. The 20 top-ranked PubChem descriptors deduced from the Gini index as derived from the RF model can be found in Fig. 8 and the contributions of their substructure towards the overall functioning of compounds as shown in Tables 3 and 4, will be discussed in the following section.

For ERα, the Gini index pertaining to the top 20 ranked features are shown in Fig. 8; and described in Table 3, which consisted of descriptors pertaining to the following classes: 7 aromatic (2 of which contain sulfur), 6 nitrogen containing features (consisting of amine and amide), 2 non-aromatic sulfur containing compound and 5 aliphatic hydrocarbons or atom counts. In addition, the Gini index pertaining to the top 20 ranked features for ERβ are shown in Fig. 8 and described in Table 4, which consisted of descriptors pertaining to the following classes: 10 nitrogen containing features (consisting of amine, amide and atom count), 3 aromatic, 3 alcohol, 4 aliphatic hydrocarbons or atom counts.

Aromatic fingerprints

For ERα, the maximum number of PubChem fingerprints as obtained from the Gini index, with 7 out of 20 top-ranked features (i.e., PubChemFP199, PubChemFP193, PubChemFP714, PubChemFP777, PubChemFP259, PubChem804 and PubChemFP741) was seen to pertain to the aromatic group. On the other hand, for ERβ, 3 out of the 20 top-ranked Gini features belonged to the aromatic group (i.e., PubChemFP777, PubChemFP714 and PubChemFP259). Surprisingly, all of the 3 aromatic Gini features for ERβ were also present for ERα. Therefore, we can infer that these 3 aromatic features are important to the functioning of the compounds. Diving deeper into the substructure description of these PubChem features, it can be seen that the first, second and ninth ranked features, (i.e., PubChemFP199, PubChemFP193 and PubChemFP259) correspond to aromatic rings of size ≥ 3 or 4. This point is in accordance to the fact that the natural agonist of ER (i.e. estradiol or E2) has only 1 aromatic ring while the most common antagonists of ER (i.e. tamoxifen and fulvestrant) have 3 and 1 aromatic rings, respectively (Bafna et al., 2020). In addition, it is the phenyl group in E2 that forms hydrogen bonds with Glu353 and Arg394 in the active site. Similarly, tamoxifen with its triphenylethylene core forms the same hydrogen bonds as E2 (Bafna et al., 2020). Thus, the phenyl moiety of compounds are important for ER inhibition. In addition, the third and seventh ranked features i.e., PubChemFP714 and PubChemFP777, correspond to 4-methylphenol which is an organic compound used as a precursor or intermediate for the manufacture of other chemicals. Furthermore, 4-methylphenol is also vital in the production of butylated hydroxytoluene (BHT) undergoing coupling to give an extensive family of diphenol antioxidants. These antioxidants are valued because they are relatively low in toxicity (Fiege, 2000). Although BHT has also been postulated as an antiviral drug, it has not yet been approved by any drug regulatory agency for use as an antiviral (Pirtle, Sacks & Nachman, 1986; Lanigan & Yamarik, 2002).

Furthermore, the tenth and fifteenth ranked features (i.e. PubChemFP804 and PubChemFP741) correspond to 3-sulfonyl phenol which according to the SMILES from its substructure description, seems to fit as a part of 4,4′-sulfonyldiphenol (Bisphenol S) (National Center for Biotechnology Information, 2020). As described in the report, Bisphenol S is an organic compound that has many functions, one of them being to act as an endocrine disruptor and it could thus modulate hormone receptors such as ERα and ERβ (Rochester & Bolden, 2015). For example, Viñas and Watson (Viñas & Watson, 2013) studied the nongenomic effects of Bisphenol S since it is an analogue of Bisphenol A, a well-known endocrine disruptor that imperfectly mimics the effects of physiologic estrogens via membrane-bound estrogen receptors. The authors concluded that Bisphenol S disrupts E2-induced cell signaling, leading to altered cell proliferation and cell death. Hence, the presence of such features in the Gini index top 20 is valid.

Nitrogen-containing fingerprints

According to the Gini index for ERβ, 10 out of the top 20 features pertain to the nitrogen containing class which includes PubChemFP392, PubChemFP451, PubChemFP439, PubChemFP393, PubChemFP299, PubChemFP450, PubChemFP645, PubChemFP15, PubChemFP375 and PubChemFP646. On the other hand, features containing nitrogen as obtained from the Gini index for ERα constituted of 6 features, namely, PubChemFP345, PubChemFP540, PubChemFP365, PubChemFP453, PubChemFP391 and PubChemFP418. Taken together, features with nitrogen (i.e., amines and amides) constituted the maximum number of features, spanning both ER subtypes. Interestingly, there were no overlapping features found between the two groups for this category. Furthermore, the substitution of the CH group with a N atom in compounds containing aromatic and heteroaromatic ring systems, is a common bioisosteric transformation conducted to mimic the binding of natural ligands while exerting antagonistic effects (Kumar et al., 2011). Of note, six out of the 10 nitrogen containing features for ERβ (Table 4) belonged to the top 10 of Gini index, emphasizing their importance. In addition, most of the aforementioned features pertain to methanamine, N-methylmethanamine, ethylamine, N,N-dimethylmethanamine etc. which are all precursors of many significant chemical compounds such as Tamoxifen, 4-hydroxy tamoxifen, Raloxifene and many of their derived analogues. A review by Sharma, Kumar & Narasimhan (2018) highlights the substructures of ERα antagonists and their analogues which were analyzed in silico using molecular docking. Through this review, the authors emphasized the need for selective estrogen receptor antagonists for the treatment of breast cancer.

Aliphatic hydrocarbons

From the analysis of the Gini index, for ERα and ERβ, respectively, 2 and 4 features out of the top 20 ranked features belonged to this group. PubChemFP697 and PubChemFP696 corresponding to 2-methlyheptane (which is also isomeric to octane) and octane, respectively, were observed to be a common feature for both ER subtypes (Tables 3 and 4). In addition, Shoda et al. (2015) observed that long alkyl side chains of tamoxifen derivatives acted as inducers of ER degradation. Furthermore, the same group designed a series of diphenylalkane derivatives bearing several long alkyl chains on the hydroxyl group and evaluated their biological properties such as ER degradation, binding affinity, transcriptional activity and anti-proliferation activity (Shoda et al., 2015; Misawa et al., 2017). Among all the compounds, one emerged as a novel ERα downregulator with a binding affinity of IC₅₀ = 4.9 nM and ERα antagonistic activity of IC₅₀ = 45 nM. Moreover, upon conducting computational docking analysis of the novel compound, the authors observed the interactions between the hydrogen atom on the amino group of the compound and the carboxylic acid of Glu351 of ERα which leads to the binding of the long alkyl chain to the hydrophobic groove of ERα. Thus, the amino group and the optimal length of the long alkyl chain in the diphenylheptane skeleton are considered important for ERα downregulation (Nanjyo et al., 2019).

Structural interpretation

As previously mentioned, the LBD of ERα and ERβ show a consensus of around 55% and are composed of 12 α helices. Upon interaction with its natural ligand (i.e., E2), both ER subtypes form hydrogen bonds with residues (ERα/β numbering) Glu353/305, His524/475, Arg394/346 and a water bridge for Arg394/346 which connects to the A-ring of E2 (Salentin et al., 2015). In addition, both ERs also form perpendicular pi-stacking with the phenol rings of residue Phe404/356 and E2. Similarly, hydrophobic interactions were also observed with the interactions of both ER subtypes to E2 (consisting of residues M343, Leu346, Leu349, Ala350, Leu387 and Leu391 for ERα and Leu298, Leu301, Ala302, M336, Leu339, Leu343 and Leu476 for ERβ) however, ERβ possesses an additional interaction with residue Leu476 as compared to ERα (Salentin et al., 2015). Fig. 9 shows these aforementioned interactions, although some interactions are not shown for simplicity.

Generally, a good ER antagonist should possess two OH groups that are linked by a lipophilic central scaffold, which places them at a distance of about 11 Å. At least one of these hydroxyls should be a phenol or a phenol-bioisostere (i.e. replacement of the phenol group with another group that can act as a similar hydrogen bond donor) (Sessler et al., 2017). One of these OH groups form a strong hydrogen bond with residues (ERα/β numbering) Glu353/305, Arg394/346 and a water molecule and the other OH group, instead, mimics the OH group of estradiol and forms an additional hydrogen bond with residue His524/475 in the estradiol D-ring pocket. At the other end of the binding cavity, the D-ring makes non-polar contacts with Ile424/376, Gly521/472 and Leu525/476 (Brzozowski et al., 1997). Beyond these key, energetically important interactions, which appear to be comparable between the two ERs, subtype-selectivity needs to arise from the different shape and hydrophobicity of the central scaffold. This has previously been demonstrated by Pike et al. (1999), where the authors study the interactions of ERβ-Genistein (a ERβ-specific partial agonist) complex and observed that the H12 helix does not adopt the distinctive ‘agonist’ conformation but instead, lies in the ‘antagonistic’ position. In addition, it was further observed that two amino acid changes within the binding cavity (i.e., residues Leu384/Met336, and Met421/Ile373) may be responsible for specificity towards ERβ as well as differences in the propensity of the two ERs to form pockets of different sizes. In any case, the selective activation of ERα or ERβ may depend not only on a selective receptor binding affinity, but also on selective activation of each receptor subtype (Leitman et al., 2010).

Several studies have reported a novel class of Selective Estrogen Receptor Modulators (SERMs) based on dihydrobenzoxathiine (i.e., the top scaffold for ERα actives) as shown in Fig. 7 (Kim et al., 2004; Chen et al., 2004; Viglianisi & Menichetti, 2010). A SAR evaluation of several analogs was carried out for identifying the best structural features for SERM both in vitro and in vivo. It appears that the sulfur atom of the dihydrobenzoxathiine scaffold interacts with discriminative residues (Leu384 for ERα and Met336 for ERβ) in the binding pocket of the two receptor isoforms and thus plays a crucial role in maintaining the subtype-selectivity. This particular scaffold was observed in 5 of the top 10 active compounds (i.e., compounds 4, 7, 8, 9 and 10 of Fig. S1) and therefore, its importance in preserving high activity towards ERα is evident. Another study revealed the diastereomerism of dihydrobenzoxathiine through molecular modeling (Zhuang et al., 2011) and discovered that the full antagonistic activity was achieved through hydrogen bonding with Glu353 and His524 of ERα LBD while van der Waals interactions were most predominant in the binding. In addition, lasofoxifene (derived from diaryltetrahydronaphthalene) was seen to be another major scaffold found in both ERα and ERβ active groups. Lasofoxifene selectively binds to both ERs with high affinity that is found to be similar to that of estradiol as well as other reported SERMs (i.e., raloxifene and tamoxifen) (Gennari et al., 2006). Furthermore, the crystal structure of lasofoxifene bound to ERα revealed features of ERα/SERM recognition whereby the C-terminal AF-2 helix is displaced at the LXXLL motif of coactivator proteins (occupying the space normally filled by residue Leu 540) as well as modulating the conformation of helix 11 residues (His 524, Leu 525) (Vajdos et al., 2007).

Additionally, the 2-phenylnaphthalene scaffold was seen to have a high frequency number for ERβ (Fig. S2). Genistein is a well known ERβ selective inhibitor that is derived from the 2-phenylnaphthalene scaffold. Furthermore, genistein is a key compound in the further identification of ERβ specific inhibitors (Mewshaw et al., 2005). The authors introduced substitutions at the appropriate positions of the phenylnaphthalene scaffold allowing a single orientation to predominate, which accounted for higher ERβ selectivity. Furthermore, Wilkening et al. (2006) noted that by creating analogues of tetrahydrofluorenone through substitution, potent ERβ subtype selective ligands were formed. It should be noted that tetrahydrofluorenone accounted for 3 of the top 10 from ERβ actives (i.e., compounds 23, 25 and 26) as shown in Fig. S2. Furthermore, the authors had also reported several analogs possessing ERβ binding affinities that is comparable to that of 17β-estradiol but with greater than 75-fold selectivity over that of ERα. Therefore, this analysis could help guide the process of novel ER inhibitor development that are subtype-selective.

Model deployment as a web server

In order to make the prediction model presented herein a practical tool that can be widely used by the scientific community, we have constructed a web server called the ERpred using the model as described in previous sections. Briefly, the SMILES notation of the chemical compound of interest is used as the input which is fed into the ERpred web server. The server performs descriptor calculation using the PaDEL-Descriptor software and outputs predictions in the form of the bioactivity class label based upon the constructed random forest model. In addition, other web servers established for bioactivity prediction do so for just one class however, the ERpred web server is unique since it can predict the bioactivity for 2 proteins (herein ERα and ERβ), simultaneously. The Shiny package under the R programming environment was used to establish the web interface which has been made freely accessible at http://codes.bio/erpred/. Fig. 10 provides screenshots of the ERpred web server where the web server prior to input submission (i.e., panel A) and after prediction (i.e., panel B) are clearly depicted.

Briefly, a step-by-step guide on using the web server is given below:

• Step 1. The following URL (http://codes.bio/erpred/) should be entered into the web browser.

• Step 2. SMILES notation for common compounds of interest can be obtained from public databases such as ChEMBL, PubChem or ChemSpider whereas custom molecules can also be drawn into ChemDraw or ChemAxon MarvinSketch so as to generate the SMILES notation.

• Step 3. Input molecules as SMILES notation should be entered into the Input box or uploaded as a file containing the SMILES notation by clicking on the “Choose file” button.

• Step 4. Prediction process can be initiated upon clicking the “Submit” button.

• Step 5. Prediction results are automatically displayed in the grey box found below the “Status/Output” heading. Typically, it takes a few seconds for the server to process the task. Users can also download the prediction results as a CSV file by pressing on the “Download CSV button”.