Bioactive Molecule Prediction Using Extreme Gradient Boosting

Following the explosive growth in chemical and biological data, the shift from traditional methods of drug discovery to computer-aided means has made data mining and machine learning methods integral parts of today’s drug discovery process. In this paper, extreme gradient boosting (Xgboost), which is an ensemble of Classification and Regression Tree (CART) and a variant of the Gradient Boosting Machine, was investigated for the prediction of biological activity based on quantitative description of the compound’s molecular structure. Seven datasets, well known in the literature were used in this paper and experimental results show that Xgboost can outperform machine learning algorithms like Random Forest (RF), Support Vector Machines (LSVM), Radial Basis Function Neural Network (RBFN) and Naïve Bayes (NB) for the prediction of biological activities. In addition to its ability to detect minority activity classes in highly imbalanced datasets, it showed remarkable performance on both high and low diversity datasets.


Introduction
Recent advancement in technology has been crucial to the explosive growth in the amount of chemical and biological data available in the public domain. Hence, data driven drug discovery and development process has attracted increased research interest in the last decade with a view not only to design and analyze but apply effective learning methodologies to the rapidly growing data. By leveraging one of the important principles of chemical/molecular similarity [1], where similar biological activities and properties are expected of structurally similar compounds, approaches to drug design through screening of large chemical databases have increased over the years. Virtual Screening (VS), the use of computational approaches and tools through the search of large databases for target or activity prediction, has notably witnessed a shift in trend from the traditional similarity searching, through reference compounds, to the use of machine learning tools to learn from the massive big data by training and prediction of unknown activity. In particular, the compound classification, in which compound label prediction is based on knowledge acquired from a training set, has gained increased research interest and many machine learning tools have been proposed to exploit the increasing big data in drug discovery. Support Vector Machines (SVM) [2,3], DT [4], Random Forest [5], K Nearest Neighbors (K-NN) [6], Naïve Bayes Classifier [7] and Artificial Neural Networks (ANN) [8] are some of the most popular machine learning methods used for activity prediction in compound classification [9].
Despite records of successful application of these methods in cheminformatics and computer aided drug discovery, each method has its peculiar shortcomings and practical constraints; such as predictive accuracy, robustness to high dimensionality and irrelevant descriptors, model interpretability, and computational efficiency, that hinders its optimal performance. For example, DT is a method that performs fairly well when it comes to most of the afore-stated criteria; however, its low predictive accuracy has inspired methods involving an ensemble of trees to improve this shortcoming. One of such efforts produced Random Forest which has been shown to be a reliable machine learning tool for compound classification as reported in [5]. In the same vein, while bearing in mind the No free Lunch Theorem [10]; that there is no best algorithm for all problems, we present herein the findings on another impressive ensemble of tree method called Extreme Gradient Boosting (Xgboost) for bioactive molecule prediction.
Xgboost is an efficient and scalable variant of the Gradient Boosting Machine (GBM) [11] which has been a winning tool for several Machine learning competitions [12,13] in recent years due to its features such as ease of use, ease of parallelization and impressive predictive accuracy. In addition to the obvious fact that alternative approaches to target prediction gives a wider perspective of the data rather than a single approach [14], we show in this paper that Xgboost not only produces comparable or even better predictive accuracy than the state of art in bioactivity prediction, but possess the intrinsic ability to handle the highly diverse and complex feature space of descriptors, especially in situations where the class distribution is highly imbalanced.

Tree Ensemble
As described by Chen and Guestrin [15], Xgboost is an ensemble of K Classification and Regression Trees (CART) {T 1 (x i , y i ) . . . ..T K (x i , y i )} where x i is the given training set of descriptors associated with a molecule to predict the class label, y i . Given that a CART assigns a real score to each leaves (outcome or target), the prediction scores for individual CART is summed up to get the final score and evaluated through K additive functions, as shown in Equation (1): where f k represents an independent tree structure with leaf scores and F is the space of all CART. The regularized objective to optimize is given by Equation (2): The first term is a differentiable loss function, l, which measures the difference between the predictedŷ and the target y i . The second is a regularization term Ω which penalizes the complexity of the model to avoid over-fitting. It is given by Ω p f q " γT`1 2 λ ř T j´1 w 2 j Where T and w are the number of leaves and the score on each leaf respectively. γ and λ are constants to control the degree of regularization. Apart from the use of regularization, shrinkage and descriptor subsampling are two additional techniques used to prevent overfitting [15].
Training. For a training dataset of molecules with vectors of descriptors and their corresponding class labels or (e.g., active/inactive) or activity of interest, the training procedure in Xgboost is summarized as follows; i For each descriptor, ‚ Sort the numbers ‚ Scan the best splitting point (lowest gain) ii Choose the descriptor with the best splitting point that optimizes the training objective iii Continue splitting (as in (i) and (ii)) until the specified maximum tree depth is reached iv Assign prediction score to the leaves and prune all negative nodes (nodes with negative gains) in a bottom-up order v Repeat the above steps in an additive manner until the specified number of rounds (trees K) is reached.
Since additive training is used, the predictionŷ at step t expressed aŝ And Equation (2) can be written as And more generally by taking the Taylors expansion of the loss function to the second order q are respectively first and second order statistics on the loss function. A simplified objective function without constants at step t is as follows The objective function can be written by expanding the regularization term as where I j " ti|qpx i q " ju is the instance set of leaf j, for a given structure q pxq the optimal leaf weight, wj , and the optimal objective function which measure how good the structure is are given by Equations (8) and (9) respectively where G j " ř iPI j g i G j " ř iPI j g i and H j " ř iPI j h i . Equation (10) is used to score a leaf node during splitting. The first, second and third term of the equation stands for the score on the left, right and the original leaf respectively. Moreover, the final term, γ, is regularization on the additional leaf.

Machine Learning Algorithms
The performance of Xgboost was compared with four machine learning algorithms that have been used in the previous studies for activity prediction (Lavecchia 2015):The Support Vector Machine LibSVM (LSVM) [16], Random Forest (RF) [5], Naïve Bayes (NB) [17], and the Radial Basis Function Network (RBFN) [18] Classifiers.

Datasets
This work was evaluated on seven carefully selected datasets that have been used to validate fingerprint based molecule classification and activity prediction in the past. A description of COX2 cyclooxygenase-2 inhibitors (COX2) (467 samples), benzodiazepine receptor (BZR) (405 samples) and estrogen receptor (ER) (393 samples) datasets [19,20] is shown in Table 1. The compounds are classified as active or inactive, and divided into training (70%) and validation (30%) sets for the purpose of this work. The table shows the mean pairwise Tanimoto similarity that was calculated based on ECFC_4 across all pairs of molecules for both active and inactive molecules. The fourth dataset utilized as a part of this study is Directory of Useful Decoys (DUD), which was presented by [21]. Although recently compiled as a benchmark data, its use in virtual screening can be found in [22,23].The decoys for each target have been chosen to fulfill a number of criteria to make them relevant and as unbiased as possible. Only 12 subsets of the DUD with only 704 active compounds were considered and divided into training (70%) and validation (30%) set in this study as shown in Table 2. The last three datasets (MDDR1-3), selected from the MDL Drug Data Report MDDR [24], have been previously used for LBVS [22,25] and activity prediction [26]. The MDDR data sets contain well defined derivatives and biologically relevant compounds that were converted to Pipeline Pilot's ECFC_4 fingerprints and folded to give 1024 element fingerprints. A detailed description of each dataset showing the training (70%) and validation (30%) sets, activity classes, number of molecules per class, and their average pairwise Tanimoto similarity across all pairs of molecules is given in Tables 3-5. The active molecules for each dataset were used. For instance, the MDDR1 (Table 3) contains a total of 8294 active molecules, which is a mixture of both structurally homogeneous and heterogeneous active molecules (11 classes). The MDDR2 (5083 molecules) and MDDR3 (8568 Molecules) in Tables 4 and 5 respectively, contain 10 homogeneous activity classes and 10 heterogeneous ones respectively [27]. The datasets were divided into training (70%) and validation (30%) sets for the purpose of this experiment. Ten-fold cross-validation was used for the Training set. In this cross-validation, the data set was split into 10 parts; 9 were used for training and the remaining 1 was used for testing. This process is repeated 10 times with a different 10th of the dataset used to test the remaining 9 parts during every run of the 10-fold cross validation. Figure 1 pictorially illustrates the various stages involved in the work under study.

Xgboost and Machine Learning Algorithms Parameters
Identifying the optimal parameters for a classifier can be time consuming and tedious and Xgboost is not an exception. This is even more challenging in Xgboost due to the wide range of tuneable parameters for optimal performance; a few of which, using the R [28] implementation of Xgboost, we have restricted our scope to in this work. Thus, by using brute force, we obtained the best performance for Xgboost when eta, gamma, minimum child weight and maximum depth were 0.2, 0.16, 5 and 16 respectively. Where; eta is the step size shrinkage meant to control the learning rate and over-fitting through scaling each tree contribution, gamma is the minimum loss reduction required to make a split, minimum child weight is the minimum sum of instance weight needed in a child and max depth is the maximum depth of a child. Other tree booster parameters like maximum delta step, subsample, column sample and the number of trees to grow per round are left at their default values of 1 respectively. For LSVM, WEKA workbench offers a way to automate the search for optimal parameters. By using grid search, a peak performance with the radial basis kernel was obtained when gamma and cost were 5.01187233627273 × 10 4 and 20 respectively. RF performed best when the maximum depth of tree was not constrained and the number of iteration set to its default value of 100. The NB classifier achieved best performance when kernel estimator parameter is used instead of normal distribution. For RBFN, we converted numeric attributes to nominal and set the minimum standard deviation to 0.1 to get the best performance.

Evaluation Metrics
The choice of performance evaluation for both model building and validation have been carefully selected from the most commonly used metrics in the literature. The selected evaluation metrics includes the accuracy, area under curve (AUC), sensitivity (SEN), specificity (SPC) and F-measure (F-Sc). The one run definition of AUC (Equation (11)) also known as balanced accuracy which is given by the average of the sum of sensitivity and specificity has been used in this work. AUC = ((SEN + SPC))/2 (11) while sensitivity (SEN) (Equation (12)) and specificity (SPC) (Equation (13)

Xgboost and Machine Learning Algorithms Parameters
Identifying the optimal parameters for a classifier can be time consuming and tedious and Xgboost is not an exception. This is even more challenging in Xgboost due to the wide range of tuneable parameters for optimal performance; a few of which, using the R [28] implementation of Xgboost, we have restricted our scope to in this work. Thus, by using brute force, we obtained the best performance for Xgboost when eta, gamma, minimum child weight and maximum depth were 0.2, 0.16, 5 and 16 respectively. Where; eta is the step size shrinkage meant to control the learning rate and over-fitting through scaling each tree contribution, gamma is the minimum loss reduction required to make a split, minimum child weight is the minimum sum of instance weight needed in a child and max depth is the maximum depth of a child. Other tree booster parameters like maximum delta step, subsample, column sample and the number of trees to grow per round are left at their default values of 1 respectively. For LSVM, WEKA workbench offers a way to automate the search for optimal parameters. By using grid search, a peak performance with the radial basis kernel was obtained when gamma and cost were 5.01187233627273ˆ10 4 and 20 respectively. RF performed best when the maximum depth of tree was not constrained and the number of iteration set to its default value of 100. The NB classifier achieved best performance when kernel estimator parameter is used instead of normal distribution. For RBFN, we converted numeric attributes to nominal and set the minimum standard deviation to 0.1 to get the best performance.

Evaluation Metrics
The choice of performance evaluation for both model building and validation have been carefully selected from the most commonly used metrics in the literature. The selected evaluation metrics includes the accuracy, area under curve (AUC), sensitivity (SEN), specificity (SPC) and F-measure (F-Sc). The one run definition of AUC (Equation (11)) also known as balanced accuracy which is given by the average of the sum of sensitivity and specificity has been used in this work.
AUC " ppSEN`SPCqq{2 (11) while sensitivity (SEN) (Equation (12)) and specificity (SPC) (Equation (13)) show the ability of the model to correctly classify true positive as positive and true negative as negative respectively, AUC simply describes the tradeoff between them.
SEN " tp{ptp`fnq (12) SPC " tnptn`fpq (13) where tp, tn, fp and fn are true positive, true negative, false positive and false negative respectively. In addition to the accuracy (Equation (14)) which is the sum of the correctly classified divided by the total number of classes, F-measure (FSc) (Equation (15)), which is the harmonic mean of precision and recall is included to serve as measure the model's accuracy.
ACC " pptp`tnq{ptp`tn`fn`fpqq F Sc " 2 pprecisionˆrecallq{pprecision`recallq (15) This work aims to introduce Xgboost for activity prediction through its performance on known datasets in drug discovery. To achieve this aim, the performance of Xgboost was compared with four state of the art machine learning algorithms used in drug discovery based on the afore-stated evaluation metrics. The prediction performances of the different machine learning algorithms on the datasets under study are tabulated in Tables 6-12. The best values for each metric is shaded.   The classification performance of the MDDR1-3, DUD, COX2, BZR and ER datasets are reported in Tables 6-12 respectively.
The experimental results on MDDR1-3 Validation datasets (Tables 6-8) shows that Xgboost produced the best accuracy, sensitivity, specificity, AUC and F-Sc across all the activity classes compared to the other machine learning methods (RF, LSVM, RBFN and NB) despite the obvious imbalance distribution of activity classes in the most of the datasets. Hence, the Xgboost method performed well for the high diverse dataset (MDDR3), and these results are particularly interesting since the MDDR3 is made up of heterogeneous activity classes which are more challenging for most machine learning algorithms.
For DUD Validation dataset (Table 9), Xgboost and RF produced the best accuracy (0.9471) compared to the other methods. In addition, Xgboost produced the best specificity across all DUD sub datasets. However, NB obtained the best sensitivity, AUC and F-Sc results.
For COX2, ER and BZR Validation datasets (Tables 10-12), it is shown that Xgboost performed well and produced the best accuracy and AUC for COX2 and ER datasets. In addition, it obtained the best F-Sc results for COX2 dataset compared to the other state-of-art methods.
Visual inspection of the results shows that Xgboost produced the best accuracy for all used datasets (except for BZR dataset which produced the second best accuracy). While the performance of Xgboost on most activity classes in terms of accuracy and AUC remains the best, it still produces the best average performance across all evaluation metrics. In addition, the good performance of Xgboost is not only restricted to homogenous activity classes since it also performed well on the heterogeneous dataset.
Moreover, a quantitative approach using Kendall W test of concordance was used to rank the effectiveness of all used methods as shown in Table 13. This test shows whether a set of raters make comparable judgments on the ranking of a set of objects. Hence, the XGB, RF, LSVM, RBFN and NB methods were used as the raters, and the accuracy measure (using MDDR1-3, DUD, COX2, BZR and ER datasets respectively) were used as the ranked objects. The outputs of this test are the Kendall coefficient (W) and the associated significance level (p value). In this paper, if the value is significant at a cutoff value of 0.01, then it is possible to give an overall ranking for the methods.
The results of the Kendall analysis for the seven datasets are shown in Table 13. The columns show the evaluation measure, the value of the Kendall coefficient (W), the associated significance level (p value), and the ranking of prediction methods. The overall rankings of the four methods show that Xgboost significantly outperforms the other methods using accuracy measure across all datasets.

Conclusions
This paper investigated the performance of Xgboost on bioactivity prediction and found out that Xgboost is indeed a robust predictive algorithm. Experimental results show that Xgboost is not only effective as a predictive model for homogeneous dataset but can replicate such effectiveness on structurally heterogeneous dataset. Experimental results show that Xgboost produces an impressive predictive accuracy, ranging from 94.47% accuracy in the heterogeneous data to 98.49% in the homogeneous one. In addition to the obvious fact that Xgboost has been shown in this work to be a good predictive tool for bioactive molecule, we are hopeful that by this Xgboost would be seen as an invaluable addition to already known computational approaches to target prediction and thus leading to a wider perspective of the data rather than a single approach.