Methodology of aiQSAR: a group-specific approach to QSAR modelling

Background Several QSAR methodology developments have shown promise in recent years. These include the consensus approach to generate the final prediction of a model, utilizing new, advanced machine learning algorithms and streamlining, standardization and automation of various QSAR steps. One approach that seems under-explored is at-the-runtime generation of local models specific to individual compounds. This approach was quite likely limited by the computational requirements, but with current increases in processing power and the widespread availability of cluster-computing infrastructure, this limitation is no longer that severe. Results We propose a new QSAR methodology: aiQSAR, whose aim is to generate endpoint predictions directly from the input dataset by building an array of local models generated at-the-runtime and specific for each compound in the dataset. The local group of each compound is selected on the basis of fingerprint similarities and the final prediction is calculated by integrating the results of a number of autonomous mathematical models. The method is applicable to regression, binary classification and multi-class classification and was tested on one dataset for each endpoint type: bioconcentration factor (BCF) for regression, Ames test for binary classification and Environmental Protection Agency (EPA) acute rat oral toxicity ranking for multi-class classification. As part of this method, the applicability domain of each prediction is assessed through the applicability domain measure, calculated on the basis of the fingerprint similarities in each local group of compounds. Conclusions We outline the methodology for a new QSAR-based predictive tool whose advantages are automation, group-specific approach to modelling and simplicity of execution. Our aim now will be to develop this method into a stand-alone software tool. We hope that eventual adoption of our tool would make QSAR modelling more accessible and transparent. Our methodology could be used as an initial modelling step, to predict new compounds by simply loading the training dataset as an input. Predictions could then be further evaluated and refined either by other tools or through optimization of aiQSAR parameters. Electronic supplementary material The online version of this article (10.1186/s13321-019-0350-y) contains supplementary material, which is available to authorized users.


Introduction
The use of quantitative structure-activity relationship (QSAR) methods has expanded significantly in recent decades [1,2]. Consequently-and due to the success of early QSAR models-there is growing interest in these methods from the regulatory perspective, which in turn influences further QSAR development [3,4]. Some of the driving factors are certainly the very high costs-financially, but also timewise [5]-required for standard in vivo and in vitro regulatory tests, as well as ethical concerns related to the use of animals for in vivo tests [6].
Recent developments of QSAR methods show several characteristics. The first is the move towards an integrated approach in modelling [7]. This means that the final output of a model is a consensus of several predictions, each obtained by a distinct QSAR approach [8][9][10]. This is further evidenced by the structure of community-wide modelling efforts, such as the ongoing initiative organized by the acute toxicity workgroup (ATWG) of the Interagency Coordinating Committee on the Validation of Alternative Methods (ICCVAM) [11], which shifted the paradigm from the commonly used approach of selecting a best method, to an integrated approach, combining predictions (with different weights) of all models developed within the initiative.
The second characteristic is streamlining and automation of the workflow. The same shift can be observed in various fields that are emerging from big-data science, such as bioinformatics [12]. Briefly, automation and standardization tools are increasing in number, aiming to streamline some of the most frequently used steps in QSAR development. These range from general tools for model development [13,14], to procedures that automatically manipulate datasets [15] and consequently, tools to simplify data conversions and interactions between various platforms used in different steps [16]. A very good outcome of these streamlining efforts is more transparency in QSAR model development, especially in terms of training set preparation [17] and methods for performance validation [18].
The last characteristic is the shift from simple, linear QSAR methods towards more sophisticated machinelearning algorithms, which were shown to perform significantly better [19,20]. However, this comes at a cost. Linear regression provided a very straightforward way of interpreting the model since it was immediately obvious which descriptors contributed to prediction the most. For the majority of machine learning algorithms, getting the mechanistic interpretation is very hard-if not impossible. A similar issue is seen in many computer science fields, where the analogous "black box" behavior problem of these algorithms makes any backtracking of the mathematical procedure (and therefore the interpretation of results) unfeasible [21].
These are all established though under-explored methodologies [22]. Some examples are the Food and Drug Administration (FDA) method of the T.E.S.T. software [23] and the lazar framework [24,25]. Our aim is to incorporate these developments into a new QSAR-based tool: aiQSAR. The tool should be automated so that it takes as an input the dataset of training compounds, regardless of the endpoint, and predicts values for new, unknown compounds. As a mathematical procedure, we employ an array of machine learning algorithms and integrate their predictions into the final output. The focus of the tool is on building distinct, local-group based models, specific to each compound in the dataset and generated at-the-runtime. Therefore, our models are built "from scratch" for each compound in the dataset, used for a single prediction, and the procedure is repeated from the beginning for the next compound.

Methods
In this section we use the term "training set" (TS) to indicate the dataset of compounds used for model derivation whose endpoint values are known and the term "evaluation set" (ES) to indicate the list of compounds of the same endpoint whose values are unknown. Two distinct modes in which aiQSAR method can be employed are cross-validating (re-testing) the training set and predicting the endpoint values for compounds in the evaluation set. The workflow (Fig. 1, Additional file 1) can be divided into 8 steps: Step 1-Calculating descriptors 3839 descriptors from Dragon 7 software [26] are calculated for all TS and ES compounds. All 1D and 2D available descriptors are selected. An alternative source of descriptors can be used and involves importing a numerical matrix of descriptors (calculated for all TS and ES compounds).
Step 2-Calculating fingerprints Fingerprints are calculated for all TS and ES compounds by the R package "rcdk" [27]. Two types of fingerprints, "pubchem" and "extended", are obtained by the "get.fingerprint" function [28]. Pubchem fingerprints are structural keys defined by the Pubchem database [29]. Extended fingerprints are obtained through a published fingerprinting procedure [30] in which paths up to six chemical bonds of depth are considered and the resulting sequence is folded to the length of 1024 bits.
The following steps are carried out for each compound (referred to as "the target compound") individually, for the training and the evaluation set: Step 3-Calculating fingerprint similarities Tanimoto distance [31] is used as a measure of fingerprint similarities. The distance between the target compound and all TS compounds is calculated by the R package "fingerprint" [32].
Step 4-Selecting the local group An algorithm ( Fig. 2) is used to select between 20 and 50 TS compounds that are most similar to the target compound. This similarity is based on Tanimoto distances between fingerprints of the target compound and all TS compounds (Step 3): First, all compounds that meet initial criteria are selected. Then, if the number of compounds is not met, a further refinement is done based on the average ranks of two similarity measures.
Step 4a-Training set For a target compound in the TS, the local group of similar compounds is selected from the remainder of the TS after the target compound is removed.
Step 4b-Evaluation set For a target compound in the ES, the local group of similar compounds is selected from the whole TS.
Step 5-Filtering descriptors All descriptors with a missing value for any compound in the local group or for the target compound are discarded. Then zero-variance and near-zero-variance descriptors are removed using the "nearZeroVar" function from the R package "caret" [33]. It should be emphasized that this procedure is carried out independently for each target compound, so different descriptors are likely to be selected each time.
Step 6-Predicting the target compound Several mathematical models are built from the local group and each model is used to predict the endpoint value of the target compound. In this step, the R package "caret" is used. Different methods are used depending on the endpoint type. Specific methods were selected based on their performance in an independent evaluation of a large collection of datasets (data not shown). The names given here correspond to the respective methods in the "caret" training function [34].
Step 7-Calculating the consensus value Predictions from all methods are combined into a single output value. Different procedures are used depending on the endpoint type.
Step 7a-Regression First, we discarded predictions from any of the five individual models that are more than 10% higher than the overall highest TS value, or more than 10% lower than the lowest TS value. The final prediction is calculated as a mean of the remaining values predicted by individual models. Fig. 2 Selecting the local group of compounds. T (pubchem) and T (extended) are Tanimoto distances based on "pubchem" and "extended" fingerprints, respectively (see Step 2) Step 7b-Binary classification The final prediction is the majority vote out of the seven individual models.
Step 7c-Multi-class classification The final prediction is the most frequently predicted class out of the five individual models. In case of a tie, the class that is least represented in the TS (out of those that are tied) is selected.
Step 8-Assessing the applicability domain Applicability domain measure (ADM) for the target compound is computed based on similarities between the target compound and compounds in its local group (Fig. 3). The target compound is assigned a rank between 1 and 5. The higher average similarity with compounds from the local group corresponds to the higher ADM rank. This measure can therefore be used to assess the coverage of chemical space relevant to the class of the target compound and, on the basis of the desired ADM cutoff, either to accept the final prediction as reliable or to reject it.

Results and discussion
The workflow was tested on a number of datasets. Three examples will be discussed, one for each endpoint type (regression, binary classification and multi-class classification).

Regression
The bioconcentration factor (BCF) was the endpoint tested for regression. The dataset was obtained from the VEGA software [44] (Additional file 2). BCF is the ratio of the concentration of a chemical in an organism to the concentration of that chemical in the surrounding environment in steady-state conditions [45]. Information on the dataset and performance statistics are given below (Table 1). In general, a higher ADM corresponds to greater average similarity between the target compound and compounds in its local group. Therefore, we would expect the method to perform better on compounds with a higher ADM. If we plot the R 2 statistic against the ADM class (Fig. 4a), we get the expected positive correlation.

Binary classification
Ames mutagenicity was the endpoint tested for binary classification. The dataset was obtained from the VEGA software [44] (Additional file 3). Ames test is a biological assay to assess the mutagenic potential of chemical compounds [46]. Information on the dataset and performance statistics are given below (Table 2). Again, there is a positive correlation between performance measures (Accuracy, Cohen's Kappa -κ -) and ADM (Fig. 4b).
Since the modelling procedure for binary classifications consists of building up seven different models, the default requirement for consensus is four predictions of the same class (Fig. 1). However, we can assess the performance of the method in relation to this requirement by varying the number of predictions required for a given class (in this case-Positives; Fig. 5). In this example our method performs best with the default cutoff of 4. More broadly, this way we can artificially raise or lower the likelihood of certain statistical errors (e.g. false positives, false negatives). We might want to do this if our model has some specific application. For example, if it is used in a regulatory context it would be reasonable to reduce the default requirement on positives in order to minimize the number of false negatives. We also noticed that for some highly unbalanced datasets reducing this requirement in favor of the minority class actually improved the overall performance.

Multi-class classification
Environmental Protection Agency (EPA) ranking for the acute oral toxicity in rats was the endpoint tested for multi-class classification. This dataset was analyzed as part of the acute oral systemic toxicity modelling initiative launched by the NTP Interagency Center for the Evaluation of Alternative Toxicological Methods (NICEATM) [47]. The objective of this initiative was predictive modelling of five endpoints related to acute toxicity and measured as rat oral LD 50 [48]: very toxic (binary classification), nontoxic (binary classification), LD 50 point estimates (regression), EPA ranking (multi-class classification) and GHS ranking (multiclass classification). EPA's hazard classification [49] splits chemicals into four hazard categories: category I (LD 50 ≤ 50 mg/kg) is the highest toxicity category; category II (moderately toxic) includes chemicals with 50 < LD 50 ≤ 500 mg/kg; category III (slightly toxic) includes chemicals with 500 < LD 50 ≤ 5000 mg/kg. Safe chemicals (LD 50 > 5000 mg/kg) are included in Category IV.
To start with, the training set of 8890 compounds (Additional file 4) was released to the community for model development related to the EPA toxicity ranking endpoint. Then the test set of 48,138 compounds was released and predictions for these compounds were requested for each of the five endpoints. Finally, the evaluation set of 2896 compounds (Additional file 5), which were an undisclosed subset of the test set and whose values for all endpoints were known, was used to assess the submitted predictions and rate the corresponding models. All models submitted were rated and ranked by the NICEATM panel on the basis of several performance statistics [47].
The performance statistics for the TS and the ES are given below (Tables 3 and 4 respectively). For both datasets there was a positive correlation between  Table 3 EPA acute oral toxicity ranking TS-statistics, confusion matrix and performance measures Average balanced accuracy was calculated as an average of balanced accuracies for each class performance measures (Accuracy, κ ) and ADM (Fig. 4c,  d). There was also very good concordance between the results for the TS and the ES. This is likely a consequence of the robust framework based on local models, chosen independently for each compound, which can mitigate some faults or errors of the dataset: such flaws will only affect a limited number of compounds rather than the whole ES.

Conclusions
We have proposed a new methodology for predictive QSAR modelling based on local group selection during model development and at-the-runtime execution. The main feature of aiQSAR is compound-specific building of predictive models: for each target compound, only a local group of TS compounds that are structurally similar to the target is considered. We quantified this similarity in terms of the ADM, which is based on Tanimoto distances of molecular fingerprints. Three application examples are discussed: BCF for regression, Ames test for binary classification and EPA acute oral toxicity ranking for multi-class classification. The ADM positively correlated with performance metrics in all of them.
In general, we believe that this QSAR approach based on local group selection and modelling automation can raise the overall quality of in silico predictions. This could be especially important for certain endpoints with under-represented classes, where localization in a natural way, with no sampling required in balancing procedures, might improve the chemical space for QSAR predictions. We plan to fully develop this methodology into a stand-alone software tool that could easily be run on any dataset, without prior knowledge of the intricate details involved in QSAR modelling. This should help open up the field to a wider range of users.

Table 4 EPA acute oral toxicity ranking ES-statistics, confusion matrix and performance measures
Average balanced accuracy was calculated as an average of balanced accuracies for each class