Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models — TEI-tagged view

Dejun JiangpersName, Zhenxing WupersName, Chang-Yu HsiehpersName, Guangyong ChenpersName, Ben LiaopersName, Zhe WangpersName, Chao ShenpersName, Dongsheng CaopersName, Jian WupersName, Tingjun HoupersName (2021date). “Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models.” Journal of Cheminformaticstitle·j.

Object of study — descriptor-based and molecular-graph-based machine-learning models for chemical property prediction and virtual screeningterm

Research question — whether graph neural networksterm·def learn representations that predict molecular properties more accurately and efficiently than strong descriptor-based baselines

Contribution — a multi-dataset, repeated-split comparison of eight algorithms, plus interpretability, computational-cost, data-quality, and HIV virtual-screening analyses

Full text — tagged

Could graph neural networksterm·def learn better molecular representationterm·var for drug discovery? A comparison study of descriptor-based and graph-based modelsterm Abstract Graph neural networksterm·def (GNNterm) has been considered as an attractive modelling method for molecular property prediction, and numerous studies have shown that GNNterm could yield more promising results than traditional descriptor-based methods. In this study, based on 11 public datasets covering various property endpointsterm·var, the predictive capacity and computational efficiencyterm of the prediction models developed by eight machine learningterm (ML) algorithms, including four descriptor-based modelsterm (SVMterm, XGBoostterm, RFterm and DNNterm) and four graph-based modelsterm (GCNterm, GATterm, MPNNterm and Attentive FPterm), were extensively tested and compared. The results demonstrate that on average the descriptor-based modelsterm outperform the graph-based modelsterm in terms of prediction accuracy and computational efficiencyterm. SVMterm generally achieves the best predictions for the regression tasksterm. Both RFterm and XGBoostterm can achieve reliable predictions for the classification tasksterm, and some of the graph-based modelsterm, such as Attentive FPterm and GCNterm, can yield outstanding performance for a fraction of larger or multi-task datasetsterm. In terms of computational cost, XGBoostterm and RFterm are the two most efficient algorithms and only need a few seconds to train a model even for a large dataset. The model interpretations by the SHAP method can effectively explore the established domain knowledge for the descriptor-based modelsterm. Finally, we explored use of these models for virtual screeningterm (VS) towards HIV and demonstrated that different ML algorithms offer diverse VS profilesterm. All in all, we believe that the off-the-shelf descriptor-based modelsterm still can be directly employed to accurately predict various chemical endpoints with excellent computability and interpretability.

Keywords: Graph neural networksterm·def, Extreme gradient boostingterm, Ensemble learning, Deep learningterm, ADME/Tterm prediction Introduction Molecular property modellingterm, which assists in hunting for chemicals with desired pharmacological and ADME/Tterm (absorption, distribution, metabolism, excretion, and toxicity) properties, is one of the most classical cheminformatics tasks [ 1 , 2 ]. A variety of machine learningterm (ML) approaches, such as Naive Bayes (NB) [ 3 – 5 ], k-Nearest Neighbors (k-NN) [ 6 ], logistic regression (LR) [ 7 , 8 ], support vector machine (SVMterm) [ 9 – 13 ], random forestterm·def (RFterm), [ 10 , 14 , 15 ] artificial neural network (ANN) [ 13 ] and more, have been widely employed in property prediction. More recently, the emergence of deep learningterm (DL) methods has revolutionized this traditional cheminformatics task due to their extraordinary capacity to learn intricate relationships between structures and properties [ 16 – 23 ]. The models developed by DL can be roughly classified into two categories: descriptor-based modelsterm and graph-based modelsterm [ 24 ]. As to descriptor-based DL models, molecular descriptorsterm and/or fingerprintsterm commonly used in traditional quantitative structure–activity relationship (QSARterm) models are used as the input, and then a specific DL architecture is employed to train a model [ 25 ]. As to graph-based DL models, the basic chemical information encoded by molecular graphsterm is used as the input, and then a graph-based DL algorithm, such as graph neural networksterm·def (GNNterm), is used to train a model. Similar to the convolutions on the regular data such as images and texts, GNNterm generalizes this operation to the irregular molecular graphterm·def that is a natural representation for chemical structures. More specifically, a graphterm G = ( V , E ) can be defined as the connectivity relations between a set of nodes ( V ) and a set of edges ( E ). Naturally, a moleculeterm can also be considered as a graphterm consisting of a set of atomsterm (nodes) and a set of bondsterm (edges).

Essentially, GNNterm aims to learn the representations of each atomterm by aggregating the information from its neighboring atomsterm encoded by the atom feature vectorterm and the information of the connected bondsterm encoded by the bond feature vectorterm through message passingterm across the molecular graphterm·def recursively (Fig. 1 ), followed by the state updatingterm of the central atomsterm and read-out operationterm. Then, the learned atomterm representations can be used for the prediction of molecular properties through the read-out phase [ 19 , 26 ]. The key feature of GNNterm is its capacity to automatically learn task-specific representationsterm using graph convolutionsterm while does not need traditional hand-crafted descriptorsterm and/or fingerprintsterm. The state-of-the-art accuracy of GNNterm models in property prediction has been well represented [ 17 , 24 , 27 – 32 ]. The representative GNNterm models and their statistical performances on the MoleculeNet benchmark datasets [ 32 ] are summarized in Table 1 . As we can see, their performances on the benchmark datasets vary from one to another, which may be attributed to the discrepancies on the model architectures, evaluation methods, training strategies and so on. Recently, a GNNterm method: Attentive FPterm, has gained increasing attention from the scientific community [ 27 ]. As shown in Table 1 , Attentive FPterm yields the best predictions to 6 out of 11 benchmark datasets, including 2 regression tasksterm (ESOL and FreeSolv) and 4 classification tasksterm (MUVterm, BBBP, ToxCastterm and ClinTox), highlighting its impressive performance in modelling a variety of chemical properties in comparison with several other graph-based methods. A majority of those studies claimed that graph-based modelsterm are typically superior or comparable to traditional descriptor-based modelsterm [ 24 , 30 – 35 ], and only a few studies gave the opposite conclusions [ 36 ]. For example, in 2017, Wu et al. reported MoleculeNet, a large benchmark for molecular machine learningterm, and the evaluation results illustrated that graph-based methods outperformed descriptor-based methods on most datasets [ 32 ]. Similarly, in 2019, Yang et al. introduced a novel GNNterm framework named directed message passingterm neural networks (D-MPNN), and the extensive evaluation on a large dataset collection indicated that D-MPNN consistently matched or outperformed descriptor-based methods on most datasets [ 24 ]. More recently, Korolev et al. reported a universal graphterm convolutional networks (GCNterm) architecture for the predictions of various chemical endpoints [ 33 ], and the application of GCNterm illustrated that its performance was comparable to state-of-the-art ML algorithms such as SVMterm, RFterm, and gradient boostingterm decision treesterm (GBDT).

The general workflow of GNNterm in molecular property prediction The reported GNNterm models in molecular property prediction In most of these reported studies, traditional ML models such as LR, RFterm, SVMterm (especially ‘gold standard’ RFterm) [ 31 , 37 ] were employed to develop the prediction models based on a set of individual fingerprintsterm (especially Extended Connectivity Fingerprintsterm, ECFP) [ 31 – 33 ]. However, it is well known that the performance of descriptor-based modelsterm is highly depending on the descriptors used in training and many previous studies have highlighted that ML models only based on molecular fingerprintsterm are not such well-performing [ 4 , 5 , 38 , 39 ]. In addition, little attention was paid to several newly state-of-the-art ML algorithms, such as XGBoostterm and LightGBM, which have illustrated great potentials for modelling various molecular properties [ 39 – 42 ]. Accordingly, the conclusion that graph-based methods outperform traditional descriptor-based methods still remains controversial.

The present study attempts to give a comprehensive evaluation of descriptor-based and graph-based modelsterm on 11 public datasets with different property endpointsterm·var. Four ML algorithms were used to develop the descriptor-based modelsterm, including SVMterm, extreme gradient boostingterm (XGBoostterm), RFterm and deep neural networks (DNNterm). In order to better represent the chemical and structure features of the moleculesterm for the descriptor-based modelsterm, the combination of one set of molecular descriptorsterm (206 MOEterm 1-D and 2-D descriptors) and two sets of fingerprintsterm (881 PubChem fingerprintsterm and 307 substructure fingerprintsterm) were considered, and such molecular representationsterm·var are also commonly seen and easily accessible. Three typical GNNterm architectures (GCNterm, GATterm and MPNNterm) and a state-of-the-art graph-based modelterm·def (Attentive FPterm) were used as the graph-based modelterm·def baselines, and the informationized molecular graphterm·def using atom-level or bond-level features were taken as the input. Both of the predictability and computability of these models were assessed. The results illustrate that the computational cost of the descriptor-based modelsterm is far less than that of the graph-based modelterm·def baselines, and the descriptor-based modelsterm generally yield more promising predictions than the graph-based methods. More concretely, SVMterm generally performs best on the regression tasksterm. Both RFterm and XGBoostterm are reliable classifiers for the classification tasksterm, but the graph-based modelsterm, such as GCNterm and Attentive FPterm, can also show excellent performance on some tasks. In terms of computational cost, XGBoostterm and RFterm are efficient and they only need a few seconds to train a model even for a large dataset. Moreover, the established descriptor-based modelsterm were interpreted by the Shapley additive explanationsterm (SHAP), and the important descriptors and structural features learned by the prediction models were highlighted. Finally, the developed ML models were used to conduct a virtual screeningterm (VS) study toward human immunodeficiency virus (HIV), and the results indicate that different ML models offer varied performance in identifying potential HIV inhibitors. All in all, we believe that the ready-made and light-weight descriptor-based modelsterm can reach better or comparable accuracy, computability, and interpretability to the highly complicated and specialized graph-based DL models.

Materials and methods Datasets To well compare the performance of descriptor-based and graph-based modelsterm, the dataset collection related to drug discovery used by Attentive FPterm was also adopted in this study [ 27 ]. This dataset collection contains 11 different datasets originally reported in MoleculeNet for a variety of chemical endpoints [ 32 ]. In the study reported by Xiong et al. [ 27 ], the moleculesterm that could not be successfully processed by RDKit [ 43 ] or the Attentive FPterm model were excluded from the original datasets. The details of those datasets are summarized in Table 2 . Here, three datasets were used for the regression tasksterm, including ESOL, FreeSolv, and Lipop, and the remaining eight datasets were used for the classification tasksterm, which can be further divided into the single-task datasetsterm (ESOL, FreeSolv, Lipop, HIV, BACE, and BBBP) and the multi-task datasetsterm (CilnTox, SIDER, Tox21, ToxCastterm, and MUVterm). Notably, we found that, in the ToxCastterm multi-task datasetsterm, some subdatasets are extremely imbalanced (the ratio of two classes is higher than 50) or quite small (the number of compounds is smaller than 500). Apparently, it seems reluctant to include these subdatasets for the development and assessment of ML models because of biased evaluation metric or insufficient training data, especially for traditional ML methods. One of the strengths for graph-based modelsterm is that multi-task learningterm·var can be applied for such highly imbalanced subdatasets and the corresponding statistics may be improved in comparison with traditional ML methods, but the prediction performances for such highly unbalanced subdatasets are not so convinced. Therefore, for the sake of fairness and simplification, such subdatasets were excluded directly, leading to the number of the tasks for ToxCastterm is 182, not the original number of 617. All the assessed ML models were evaluated based on the same remaining 182 ToxCastterm tasks, and we believe that the results can still make sense.

The detailed information of the datasets used in this study Molecular representationterm·var Graph-based methods are capable of learning molecular representationsterm·var by operating the convolutions on the encoded molecular graphsterm directly. In the graphterm representation for a moleculeterm, the connectivity relation between atomsterm is represented by a graphterm G = ( V , E ). Here, the nodes V are represented by the node feature vector X v consisting of a series of atomic featuresterm and the edges E are represented by the edge feature vector E km consisting of a series of bond featuresterm, where the subscript km indicates that atomsterm k and m are bonded. Followed by previous studies [ 27 ], almost all the easily accessible atomterm/bond-level features were exhausted to comprehensively squeeze chemical information into molecular graphterm·def for graph-based modelsterm, where include nine kinds of atomic featuresterm (i.e., atomterm symbol, atomterm degree, formal charge, radical electrons, hybridization, aromaticity, hydrogens, chirality and chirality type) and four kinds of bond featuresterm (i.e., bondterm type, conjugation, ring, and stereo). Most of them were encoded into a molecular graphterm·def in a one-hotterm manner and subsequently the encoded molecular graphterm·def was used as the input. The more details about the molecular representationsterm·var for graph-based modelsterm are available in the publication [ 27 ].

All the moleculesterm were minimized using the MMFF94 force field in MOEterm (Version: 2015.1001) with the default parameters. Then, the expert-crafted descriptors and fingerprintsterm were computed to develop the descriptor-based modelsterm. To comprehensively represent molecular structures, 206 MOEterm 1-D and 2-D descriptors and two sets of fingerprintsterm, including 881 PubChem fingerprintsterm (PubchemFP) and 307 substructure fingerprintsterm (SubFP), were used. The MOEterm descriptors were calculated by MOEterm (Version: 2015.1001), and the two sets of fingerprintsterm were calculated by PaDEL-Descriptor (Version: 2.1). [ 44 ] Prior to the development of the descriptor-based modelsterm, all the molecular features were pretreated as follows: (1) the features with missing values and extremely low variance (< 0.05) were removed; (2) the feature that has a high correlation ( r > 0.95) with another feature was removed; (3) the retained features were normalizedterm to the mean value of 0 and variance of 1.

Machine learningterm algorithms As one of the most classic cheminformatics problems, molecular property prediction has made considerable progress over the last decade due to the application of new ML methods represented by deep learningterm and ensemble learning [ 25 , 40 , 45 , 46 ]. In this study, four representative ML algorithms (i.e., DNNterm, SVMterm, XGBoostterm and RFterm) were used to develop the descriptor-based modelsterm, and four representative graph-based methods (i.e., MPNNterm, GCNterm, GATterm and Attentive FPterm) were employed to develop the graph-based modelsterm.

Deep neural networks (DNNterm) As one of the typical DL algorithms, DNNterm has an input layer, an output layer, and many hidden layers. DNNterm is composed of many individual neuronsterm [ 16 , 25 ]. Each neuron in DNNterm aggregates information from its connected neurons and then the aggregated information is activated by a non-linear activation function. Such manifestations mimick the behavior of biological neural networks. All the operations in DNNterm aim to learn intricate and rapidly-varying non-linear functions and extract a hierarchy of useful features from the input [ 18 ]. In this study, three hidden layers feed-forward neural networks were employed, and the following key hyper-parameters were optimized: L2 regularization (0 to 0.01), dropout rate (0.0 to 0.5) and neurons for each hidden layer (64, 128, 256, 512). The other important hyper-parameters were fixed: ReLU function that has been recommended by many previous studies was used as the activation function [ 25 , 47 ], and the optimizer was set to an adaptive learning rate algorithm: Adadelta [ 48 ].

Support vector machine (SVMterm) SVMterm is one of the most popular ML approaches and it is appropriate for both classification and regression [ 9 , 49 , 50 ]. It is also capable of dealing with both linearly separable and linearly inseparable problems. For linearly inseparable feature space, the kernel trickterm is needed to map the original feature space onto a new higher separable linear space. The basic objective of SVMterm is to find the optimal hyperplaneterm in the feature space that can maximize the distance between the data points and hyperplane, and the discriminant results generated from this optimal hyperplaneterm should be insensitive to small perturbation of training samples. Here, the commonly used radial basis function (RBF) was used as the kernel and the following main hyper-parameters were optimized: C (0.1 to 100) and gamma values (0 to 0.2).

Extreme gradient boostingterm (XGBoostterm) XGBoostterm is one of the most representative ensemble learning ML algorithmsterm under the frame of gradient boostingterm [ 51 ]. Compared with traditional gradient boostingterm, several algorithm optimizations were introduced to XGBoostterm, such as minor improvement in the loss function by penalizing the complexity of the model, introduction of shrinkage and column subsampling for further preventing over-fitting, employment of sparsity-aware split finding technique for efficient training on sparse data, etc. [ 51 ]. XGBoostterm has gained extensive attention in the property prediction due to its significant predictive power and low computational cost [ 42 , 52 , 53 ]. In the training of XGBoostterm, the following hyper-parameters were optimized: learning_rate (0.01 to 0.2), gamma (0 to 0.2), min_child_weight (1 to 6), subsample (0.7 to 1.0), colsample_bytree (0.7 to 1.0), max_depth (3 to 10) and n_estimators (50, 100, 200, 300, 400, 500, 1000).

Random forestterm·def (RFterm) Random forestterm·def is another representative ensemble learning ML algorithmsterm. It constructs a strong classifier or regressor by an ensemble of individual decision treesterm under the frame of bagging and makes predictions by majority vote or averaging of multiple decision treesterm [ 10 , 15 ]. In the implementation of RFterm algorithm, sample perturbation via bootstrap sampling of the training data and feature perturbation via random feature subset selection are introduced to improve the diversity of base learner (decision treesterm), which corrects for the overfitting habit of decision treesterm and subsequently enhances the generalization ability of RFterm. In the training of RFterm, the following hyper-parameters were optimized: n_estimators (10, 50, 100, 200, 300, 400, 500), max_depth (3 to 12), min_samples_leaf (1, 3, 5, 10, 20, 50), min_impurity_decrease (0 to 0.01) and max_features (‘sqrt’, ‘log2’, 0.7, 0.8, 0.9).

Message passingterm neural networks (MPNNterm) MPNNterm is a common framework for GNNterm that was used for chemical prediction in 2017 by Gilmer et al. [ 54 ], and it has shown versatility in many applications such as natural language processing, image segmentation, chemical/molecular graphsterm, and so on. Many recently proposed GNNterm architectures for molecular property prediction can be formulated in this flexible framework [ 24 , 26 , 34 , 37 ]. In theory, MPNNterm operates the convolutions on undirected molecular graphsterm G = ( V , E ) with node features X v and edge features E km . The forward propagation of MPNNterm has two phases: message passing phaseterm and readout phaseterm. The message passing phaseterm transmits information across the molecular graphterm·def to learn a molecular embedding using the message functions M t and node updating functions U t , and the readout phaseterm computes a feature vector for the whole molecular graphterm·def using some readout functions R to model the properties of interest. More mathematical details are available in the study reported by Gilmer et al. [ 54 ] In the training of MPNNterm, the following hyper-parameters were optimized: L2 regularization (0, 10e-8, 10e-6, 10e-4), learning rate (10e-2.5, 10e-3.5, 10e-1.5), dimension of node feature in hidden layers (64, 32, 16), dimension of edge feature in hidden layers (64, 32, 16), and number of set2set layers (2,3,4). The number of message passingterm steps and set2set steps were fixed to 6.

Graphterm convolutional networks (GCNterm) To date, various GCNterm frameworks and variants have been proposed, and the most classical GCNterm model was introduced by Kipf et al. in 2017 [ 55 ]. Mathematically, it follows the propagation rule: H l + 1 = σ D ^ - 1 2 A ^ D ^ - 1 2 H l W l , where H ( l ) and W ( l ) denote the l th neural networks layer and its corresponding learnable parameters, respectively. σ represents a non-linear activation function. Generally, D and A are the degree matrix and adjacency matrix, respectively, A ^ = A + I where I is the identity matrix, and D ^ is the diagonal node degree matrix of A ^ . The design of the D ^ - 1 2 A ^ D ^ - 1 2 term is intended to add a self-connection to each node and keep the scale of the feature vectors. From the message passingterm point of view, it can also be ascribed to the following two step: (1): aggregate neighbors’ information h v to produce an intermediate representation h ^ u ; (2) transform the aggregated representation h ^ u with a linear projection followed by a non-linearity activation: h u = σ W u h ^ u . In this study, the vanilla GCN modelterm proposed by Kipf et al. was used and the following hyper-parameters were optimized: L2 regularization (0, 10e−8, 10e−6, 10e−4), learning rate (10e−2.5, 10e−3.5, 10e−1.5), dimension of FNN classifier (64, 128, 256), and dimension of GCNterm hidden layers ([128, 128], [256, 256], [128, 64], [256, 128]).

Graphterm attention network (GATterm) GATterm is an extension of the vanilla GCN modelterm, and the biggest distinction between vanilla GCNterm and GATterm is the way of neighboring information aggregation. In the vanilla GCN modelterm, the graph convolutionterm operation aggregates the normalizedterm sum of neighboring information. In the GATterm, attention mechanisms by specifying different weights to different nodes are introduced and the corresponding graph convolutionterm operation aggregates the weighed sum of neighboring information in a formulation: H i l + 1 = σ ∑ j ∈ N i α ij l W l H i l , where α ij l is the normalizedterm attention score between node i and node j in the l th graph convolutionterm layer. W , N i and σ are learnable weight matrix, the neighbors of node i , and non-linear activation function respectively. The calculation of the attention score and other details can be reference to the corresponding publication [ 56 ]. The application of attention mechanisms in the graph convolutionterm can force the model to learn the most meaningful parts in neighbors and local environment and it has gained preferable performance in comparison with other usual GCNterm architectures [ 27 , 34 , 56 ]. In the training of GATterm model, the following key hyper-parameters were optimized for each task: L2 regularization (0, 10e−8, 10e−6, 10e−4), learning rate (10e−2.5, 10e−3.5, 10e−1.5), dimension of GATterm hidden layers ([128, 128], [256, 256], [128, 64], [256, 128]), dimension of FNN classifier (64, 128, 256), and the number of attention heads ([2, 2], [3, 3], [4, 4], [3, 4], [2, 3]).

Attentive FPterm Attentive FPterm, proposed by Xiong et al. [ 27 ] is a state-of-the-art GNNterm model for molecular property prediction. In Attentive FPterm, the recursive neural networks (RNN) was employed to agglomerate the structural information encoding in molecular graphterm·def from nearby to distant and update the state of centered atomterm. Moreover, a graph attention mechanismterm was introduced to allow the model to focus on the most relevant parts of the input. The results reported by Xiong et al. illustrated that Attentive FPterm can achieve state-of-the-art predictions to a wide range of molecular properties (Table 1 ) [ 57 ]. The main hyper-parameters for Attentive FPterm include num_layers (the number of attentive layers for atomterm embedding), num_timesteps (the number of attentive layers for moleculeterm embedding), graphterm_feat_size (fingerprintterm dimension), L2 regularization, learning rate, and dropout rate. Here, all those main hyper-parameters were optimized: L2 regularization (0, 10e-8, 10e-6, 10e-4), learning rate (10e-2.5, 10e-3.5, 10e-1.5), num_layers (2, 3, 4, 5, 6), num_timesteps (1, 2, 3, 4, 5), dropout (0, 0.1, 0.3, 0.5), and graphterm_feat_size (50, 100, 200, 300).

For the development of the four descriptor-based modelsterm, the DNNterm algorithm was implemented in the PyTorchterm package (Version: 1.3.1 + cu92) of Python (Version: 3.6.5 × 64), and the XGBoostterm (Version: 0.80), RFterm and SVMterm algorithms were implemented in the scikit-learnterm package (Version: 0.20.1) of Python [ 58 ]. All the four graph-based modelsterm were implemented by the Deep Graph Libraryterm (DGL) package (Version: 0.4.1) using PyTorchterm as the backend of Python [ 59 ].

Model training, optimization and evaluation protocols In the first stage, the same training, validation and test setsterm at a ratio of 8:1:1 used by Attentive FPterm were also used in our study (Additional file 1 generated from the source code provided in the github). For the assessed ML algorithms, the prediction on the validation set was used to guide the optimization of hyper-parameters. The Tree of Parzen Estimatorsterm (TPE) algorithm was used to identify the best hyper-parameters for different ML models in 50 evaluations (Here four graph-based modelsterm on the HIV and MUVterm datasets were in 30 evaluations due to the high computation overhead). The TPE algorithm is an optimization algorithm under the sequential model-based global optimization frame and capable of finding ideal hyper-parameters only through a few objective function evaluations. TPE was implemented by the hyperopt package (Version: 0.2) in Python (Version: 3.6.5 × 64) [ 60 ]. Then, in the second stage, in order to alleviate the effect of the randomness of data splitting, 50 independent runs with different random seedsterm for data splitting (training/validation/test = 8:1:1) were performed to evaluate each ML model in a more reliable way. Similarly, four graph-based modelsterm on the HIV and MUVterm datasets were in a 20 independent runs due to the high computation overhead, and the optimized hyper-parameters determined in the first stage were straightly adopted. For avoiding overfitting and tremendous time consumption, all the neural network (NN)-based model (i.e. DNNterm, GCNterm, GATterm, MPNNterm and Attentive FPterm) were trained in an early stoppingterm way for all tasks if no validation performance improvement was observed in successive 50 epochs, and followed by the previous DNNterm hyper-parameter recommendations [ 25 , 61 ], the maximum epoch was set as an empirical value of 300 for all the task. The additional check of the training logs also proved that this empirical value is enough to learn representative parameters for NN-based models. The training batch for most tasks was set as 128. However, this number was also merely empirical and could change depending on the complexity of model and data volumeterm·var. All the model training and evaluation scripts were available in Additional file 2 .

According to the recommendations of MoleculeNet benchmarks [ 32 ], the classification models were evaluated by the area under the receiver operating characteristic curve (AUC-ROCterm) for the classification tasksterm except the maximum unbiased validation (MUVterm) dataset, which was evaluated by the area under precision-recall curve (AUC-PRCterm) due to its extreme biased data distribution. The regression models were evaluated by root mean square error (RMSEterm). In a more diverse evaluation, we also considered mean absolute error (MAEterm) and R-Squareterm (R2) metrics for regression model. As shown in Table 2 , five datasets contain more than one task. The multi-task learningterm·var was applied in the development of the five NN-based models including DNNterm, GCNterm, GATterm, MPNNterm and Attentive FPterm for each multi-task datasetterm, and the average performance across multiple tasks was reported. However, it is not practical to generalize the multi-task learningterm·var to traditional descriptor-based modelsterm (i.e. SVMterm, XGBoostterm, and RFterm). In this case, each multi-task datasetterm was split into multiple single-task datasetsterm and the individual descriptor-based modelterm·def on each single dataset was trained, and then the average performance was reported in a similar way.

Model interpretation ML algorithms usually function as a “black-box”, and how to interpret these complicated ML models remains a big challenge. Several interpretation methods have been proposed to uncover the “black-box” essence of ML algorithms and they can be roughly classified into two major categories: model-specific and model-agnosticterm strategies. The model-specific strategies are relevant to the specific structure of a model, such as the feature weights for the simplistic linear model and feature importance determined by Gini index for RFterm model. One of the strengths for the model-agnosticterm strategies is that they do not depend on the specific model architecture and can mitigate the necessity to balance model performance and interpretability [ 62 , 63 ]. Some model-agnosticterm strategies such as sensitivity analysis have been applied in model interpretation but it becomes inefficient with the increase of model dimensionality [ 64 , 65 ].

Here, a recently-developed model-agnosticterm interpretation framework called SHapley Additive exPlanationsterm (SHAP) was employed to interpret the ML models due to its both local and global interpretability [ 66 ]. SHAP method was inspired from the game theory and the corresponding SHAP valueterm was employed to quantify the contributions of single players to a collaborative game [ 65 ]. Some published studies have demonstrated that SHAP method has high potential in understanding arbitrary complicated ML models [ 39 , 65 ]. In a more specific way, this method defines an explanation model that belongs to a linear function of binary variables: f x ≈ g z ′ = ∅ 0 + ∑ i = 1 M ∅ i z i ′ , where z ′ ∈ 0 , 1 M denotes the absence (0) or presence (1) of a certain descriptor, and M is the number of molecular descriptorsterm. ∅ i is the so-called SHAP valueterm, and similar to previous descriptions, it measures the impact of the presence or absence of a descriptorterm on the model output, and the sum of all descriptor attributions g z ′ approximates the output f x of the original model. More details about this method can be found in the relevant publications [ 39 , 65 ]. The SHAP method was implemented in the shap package (Version: 0.35.0) of Python software (Version: 3.6.5 x64).

Washing of the benchmark datasets Data qualityterm is one of the fundamental questions in cheminformatics and the incorrect or inappropriate structures contained in datasets would hinder the effort of developing reliable prediction models. Here we found that some salts, inorganics, counterions, solvents, mixtures and even duplicates with inconsistent labels existing in the datasets provided by Xiong et al. [ 27 ], but we do not remove them first for the sake of fairness. The original data was reported by Wu et al. and it is apparently unreasonable to use such datasets for model building. In this regard, we developed a python script based on MOEterm (Version: 2015.1001) and RDKit (Version: 2019.09.1) to automatically eliminate the incorrect or inappropriate structures from the original datasets with the following steps: (1) For the mixtures and compounds containing salts, counterions, and solvents, we used a compromised method of keeping the major component with the largest number of heavy atomsterm and the retained component was neutralized if possible. This step was accomplished by the sdwash module in MOEterm and the compounds that cannot be recognized by MOEterm were eliminated; (2) A moleculeterm was identified as an inorganics if it does not contain any carbon atomterm and then eliminated from the datasets. Similarly, the compounds that cannot be recognized by RDKit were also eliminated in this step; (3) Duplicates were identified by the canonical SMILESterm generated from RDKit. After that, the duplicated records with inconsistent labels were removed.

Results and discussion Performance of descriptor-based and graph-based modelsterm At the outset, the same training, validation, and test sets for the development of the Attentive FPterm models were adopted, and the corresponding statistical results for the six single-task datasetsterm including three regression tasksterm and three classification tasksterm given by the four descriptor-based and four graph-based modelsterm are summarized in Table 3 (regression tasksterm) and Table 4 (classification tasksterm).

The performance comparison (RMSEterm) of the four descriptor-based and four graph-based modelsterm on the three regression datasets (data folds were generated from Attentive FPterm and the top three model were italic for each dataset) The performance comparison (AUC_ROC) of the four descriptor-based and four graph-based modelsterm on the three classification datasets (data folds were generated from Attentive FPterm and the top three model were italic for each dataset) For the regression tasksterm, one of the graph-based modelsterm, Attentive FPterm, achieves the best statistical performance on ESOL with the RMSEterm of 0.471 for the test set, and the performances of SVMterm (RMSEterm = 0.516) and DNNterm (RMSEterm = 0.553) are slightly worse than it. As we can see, the performances of three classical graph-based modelsterm (i.e., GCNterm, GATterm and MPNNterm) and RFterm are obviously unpleasant on this dataset. For FreeSolv, both SVMterm and XGBoostterm offer considerable and comparable performances with RMSEterm = 0.674 and 0.707 for the test set respectively, which are slightly better than that of DNNterm (RMSEterm = 0.724). With regard to Lipop, three methods including one graph-based method (Attentive FPterm) and two descriptor-based methods (SVMterm and XGBoostterm) achieve similar performances with RMSEterm ≈ 0.560 for the test set and this predictive capability is superior to other methods, especially RFterm, GATterm and MPNNterm. Clearly, the RFterm and three graph-based modelsterm (i.e., GCNterm, GATterm and MPNNterm) show disappointing predictive capability to the three regression tasksterm. On average, SVMterm achieves the best predictions on the test sets of the regression tasksterm. XGBoostterm and Attentive FPterm perform similarly but slightly worse than SVMterm.

As for the three classification tasksterm including HIV, BACE and BBBP, it gets puzzled to tell which type of model, i.e. descriptor-based and graph-based, is superior in the light of statistical results only from one random partition. However, it can be observed that three descriptor-based modelsterm (i.e. XGBoostterm, RFterm, and DNNterm) and two graph-based modelsterm (i.e. GCNterm and Attentive FPterm) are more powerful than the other models in general. Concretely, GCNterm and DNNterm offer almost the same predictions to HIV with AUC-ROCterm ≈ 0.857 for the test set, and three models including XGBoostterm, Attentive FPterm and RFterm are slightly worse than them with AUC-ROCterm ≈ 0.847. Besides, XGBoostterm and Attentive FPterm give the same performances on BACE with AUC-ROCterm = 0.889 for the test set, and DNNterm is slightly inferior to them with AUC-ROCterm = 0.883 for the test set. For BBBP, both Attentive FPterm and RFterm offer the same predictive ability for the test set with AUC_ROC = 0.907, and SVMterm gives slightly worse results with AUC_ROC = 0.899 for the test set.

Next, the performances of the descriptor-based and graph-based modelsterm were further compared on the five multi-task datasetsterm including ClinTox, SIDER, Tox21, ToxCastterm, and MUVterm. As shown in Table 5 , it seems also struggling to distinguish which type of model is more promising. Here from the overall level, the models that perform well in the aforementioned three classification tasksterm (i.e. GCNterm, Attentive FPterm, XGBoostterm, RFterm and DNNterm) can still give satisfactory predictions to the five multi-task datasetsterm. More specifically, for ClinTox, two descriptor-based modelsterm (SVMterm and RFterm) and one-graph based model (GATterm) give more promising predictions than the other models. For both SIDER and Tox21, two descriptor-based modelsterm (XGBoostterm and RFterm) and one graph-based modelterm·def (Attentive FPterm) share similar and more powerful predictions on the corresponding test sets. For MUVterm, one descriptor-based modelterm·def (SVMterm) and two graph-based modelsterm (GATterm and Attentive FPterm) offer more promising results on the test set compared with the others.

The performance comparison (AUC_ROC, MUVterm: AUC_PRC) of the four descriptor-based and four graph-based modelsterm on the five multi-task classification datasets (data folds were generated from Attentive FPterm and the top three model were italic for each dataset) To our surprising, five NN-based models including DNNterm, GCNterm, GATterm, MPNNterm and Attentive FPterm yields much better prediction than three descriptor-based modelsterm to the ToxCastterm dataset (average AUC-ROCterm = 0.897 for five NN-based models and 0.760 for three descriptor-based modelsterm). The careful analysis of the Attentive FPterm source code suggests that the unreasonable data splitting for ToxCastterm may attribute to the over-optimistic predictions of five NN-based models where multi-task learningterm·var was applied. More concretely, it is quite possible that all positive samples or negative samples may occur for some columns of the data folds generated from a strongly biased subdataset in ToxCastterm based on the random data splitting and the average AUC-ROCterm cannot be calculated for such data folds accordingly (AUC-ROCterm metric calculation error). In this case, Xiong et al. adopted a compromised splitting strategy where a stratified sampling at a ratio of 8:1:1 was individually applied to each single task of ToxCastterm to generate 182 independent training sets, validation sets and test sets for 182 different tasks [ 27 ]. After that, those independent training/validation/test sets were merged one task by one task in an outer join manner to produce the final training/validation/test set. It is the fact that the aforementioned situation (AUC-ROCterm metric calculation error) was well avoided, but the issue raised by such splitting strategy is the over-estimated statistical results when multi-task learningterm·var was applied because many samples in the final test or validation sets will be included in the final training set. However, such situation (over-estimated statistical results) was well evaded by descriptor-based modelterm·def where each single task was detached to train the model individually and no duplicated samples could occur in the data folds. More details about the data splitting used by Xiong et al. could be found in their webpage [ 67 ]. Besides, it is the same manner for the splitting of the biased MUVterm dataset in Attentive FPterm. To our knowledge, a reasonable way to solve this problem is to change the random seedterm for data splitting if the randomly generated data folds suffer from such situation. Hence, the obvious inferiority of three descriptor-based modelsterm on ToxCastterm compared with five NN-based models may be reasonably explained by the over-optimistic predictions of our NN-based models (what we will discuss later).

Actually, it seems arbitrary to judge which of models is better only based on the statistical results from one-time run because of the randomness in data splitting. To evaluate the ML models in a more reliable way, 50 times independent runs based on different random seedsterm to split data into 50 different folds of training, validation, and test sets at the ratio of 8:1:1 were conducted for each dataset, and the average performance over the 50 folds with the corresponding standard deviation was used to evaluate the ML models. And the splitting strategy for the ToxCastterm and MUVterm datasets was revised. The corresponding statistical results for the 11 studied datasets given by eight assessed models are listed in Table 6 (three regression datasets), Table 7 (three single-task classification datasets) and Table 8 (five multi-task classification datasets). From the Table 8 , it can be observed that the predictions to the randomly split ToxCastterm datasets (Table 5 ) are much worse than those to the data generated by the original splitting strategy used by Attentive FPterm (average AUC_ROC of five NN-based models: 0.897 to 0.770), demonstrating the over-optimistic predictions given by five NN-based models based on the original splitting strategy. Here, it can be found that the average performance of the 50 times independent runs is worse than that of the one-time run for the 11 studied datasets. To our knowledge, many previous studies evaluated the ML models by only averaging the performance from three independent runs and their results may be sensitive to the randomness of data splitting [ 27 , 32 ]. To well illustrate this point, we counted the average performances for the top three runs and the worst three runs among the 50 times independent runs for XGBoostterm (Additional file 3 : Table S1). It can be recognized that the average performances for the top three runs and the worst three runs have big discrepancies for XGBoostterm. Therefore, with the aim of alleviating the randomness of data splitting, it is recommended to conduct sufficient independent runs to evaluate ML models more reliably.

The performance comparison (average RMSEterm) of the 50 times independent runs on the three regression datasets for the eight models. (the top three model were italic for each dataset) The performance comparison (Average AUC_ROC) of the 50 times independent runs on the three classification datasets for the eight models. (the top three model were italic for each dataset) The performance comparison (Average AUC_ROC, MUVterm: Average AUC_PRC) of the 50 times independent runs on the five multi-task classification datasets for the eight models. (the top three model were italic for each dataset) As shown in the Table 6 , it can be recognized that two descriptor-based modelsterm (SVMterm and XGBoostterm) and one graph-based modelterm·def (Attentive FPterm) generally give better performances than the other models, which is consistent, to some extent, with the findings from the previous one random split. Among them, SVMterm gives the best predictions to the ESOL and FreeSolv datasets with average RMSEterm of 0.569 and 0.852 to the test sets, respectively. Attentive FPterm gives the best predictions to the Lipop dataset with average RMSEterm of 0.553 to the test set, and SVMterm and XGBoostterm are slightly worse than Attentive FPterm with RMSEterm ≈ 0.574. Here, XGBoostterm offers satisfactory but slightly worse predictions to all the three regression datasets (average RMSEterm = 0.582, 1.025, and 0.574 to ESOL, Freesolv, and Lipop respectively) compared with SVMterm and Attentive FPterm. In addition, the MAEterm and R2 metrics given by the eight models on the three regression tasksterm were also calculated (Additional file 3 : Tables S2 and S3). As shown in Additional file 3 : Table S2 and S3, similar conclusions could be drawn where SVMterm, XGBoostterm and Attentive FPterm are well-performing regressors and on average SVMterm is the best one. Here what we found from Table 6 is that the descriptor-based modelsterm, especially SVMterm, generally show much better training set performances in comparison with the graph-based modelsterm (especially for two smallest datasets FreeSolv and ESOL). However, some graph-based modelsterm, especially Attentive FPterm, are able to reach comparable prediction results to the descriptor-based modelsterm for the test sets, implying that the descriptor-based modelsterm are more likely to be over-fitted and less generalized compared with the graph-based modelsterm learnt from small and chemically narrow datasets. As for the three single-task classification datasets shown in Table 7 , what we can find is that the four descriptor-based modelsterm are obviously superior to the four graph-based modelsterm on the BBBP dataset, where the average AUC_ROC of the four descriptor-based modelsterm is 0.924 compared with that of 0.891 for the four graph-based modelsterm. Similarly, on average the four descriptor-based modelsterm can give more reliable predictions to the BACE dataset where the average AUC_ROC of the four descriptor-based modelsterm is 0.891 compared with that of 0.875 for the four graph-based modelsterm. However, for the larger HIV, it seems that the graph-based modelsterm are slightly better than the descriptor-based modelsterm, implying that inclusion of more samples may be helpful to train a better graph-based modelterm·def. In some cases, one may need to re-train their ML models with the gradual accumulation of available experimental datasets. Such operations can benefit more to graph-based modelsterm due to their data-hungry essence, but the rapid accumulation of qualitied experimental datasets is not an easy task. On the contrary, regular re-training of ML models by adding a small number of new compounds one time could be some of routine. Generally speaking, the optimization of hyper-parameters is necessary when re-training ML models, especially for NN-based models where their performances are sensitive to the hyper-parameters such as the initial parameters and learning rate. Compared with graph-based modelsterm, descriptor-based modelsterm such as RFterm or SVMterm may be more stable for a long time. With regard to the five multi-task datasetsterm shown in Table 8 , it can be found that the descriptor-based modelsterm, especially XGBoostterm and RFterm, achieve better predictions than the graph-based modelsterm on the ClinTox, SIDER and MUVterm datasets. However, one graph-based modelterm·def, Attentive FPterm, achieves the best predictions to the two relatively large toxicity-relevant datasets including Tox21 and ToxCastterm with average AUC_ROC of 0.852 and 0.794 to the test sets, respectively, which may benefit from the multi-task learningterm·var and larger data volumeterm·var. Numerous studies demonstrated that multi-task models have advantages over single-task models due to their ability to excavate the inconspicuous hidden relations between different subtasks and transparently share the learned features among all the tasks [ 57 , 68 , 69 ]. Nevertheless, the performance of multi-task models is highly related to the favorable correlations of individual tasks but such ready-to-use tasks are not so commonly seen in practical drug discovery campaigns. For the purpose of simplicity, we counted the top three models and the corresponding performances based on the results from 50 times independent runs for each dataset. As can be seen from Table 9 , the descriptor-based modelterm·def achieves the best predictions to six out of 11 datasets including ESOL, FreeSolv, BBBP, ClinTox, SIDER and MUVterm. Moreover, it can be observed that the top three models of all the datasets were mainly occupied by the descriptor-based modelsterm (the ratio is 24/33 = 73%), substantiating the more powerful predictive abilities of the descriptor-based modelsterm compared with the graph-based modelsterm. It is possible that the superiority of the descriptor-based modelsterm for some datasets (ESOL, FreeSolv, and Lipop) may be partially contributed from the descriptors that are highly correlated to the target values (such as the ‘LogS’ descriptor for the ESOL dataset). To systematically check this problem, we removed the top three descriptors that are highly correlated to the target values according to the Pearson’s correlation coefficients (ESOL: ‘logS’, ‘h_logS’, and ‘SlogP’; FreeSolv: ‘vsa_pol’, ‘h_emd’ and ‘a_donacc’; Lipop: ‘SlogP’, ‘h_logD’, and ‘logS’) and then used the remaining descriptors to reconstruct the four descriptor-based modelsterm based on the optimal hyper-parameter configurations determined in the first evaluation stage. The evaluation metrics were also averaged from the 50 times independent runs (Additional file 3 : Table S4). It can be observed that the performance of the models developed based on the remaining descriptors do not show large difference compared with those developed based on the original descriptors. Moreover, we found that the descriptor-based modelsterm without thse high-related descriptors are still superior to the graph-based modelsterm (Additional file 3 : Table S4). Here what we found is that the graph-based modelsterm can outperform the descriptor-based modelsterm on some lager or multi-task datasetsterm such as the HIV, Tox21 and ToxCastterm datasets, which is in well accordance with the previous conclusions where DNNterm excel at larger amounts of data and multi-task learningterm·var [ 68 , 69 ]. However, to build such generalizable and robust deep models requires large-scale high-quality datasets and the datasets in the practical drug discovery campaigns routinely suffer from narrow chemical diversity and insignificant sample sizes [ 70 ]. On the ground, we believe that the descriptor-based modelsterm can be still widely used and give reliable predictions in the drug discovery campaigns.

The top three model and corresponding performances based on the results from 50 times independent runs for each dataset. (the descriptor-based modelsterm were colored as italic and the graph-based modelterm·def were colored as undeline) In conclusion, regardless of the statistical results on the same data folds used by Attentive FPterm or a more reliable 50 times independent runs, what we found is that the traditional descriptor-based modelsterm generally outperform the state-of-the-art graph-based modelsterm. Among them, SVMterm is the best algorithm in modelling regression tasksterm. Both RFterm and XGBoostterm can be well-performing in modelling classification tasksterm, and some graph-based modelsterm, such as Attentive FPterm and GCNterm, can outperform the descriptor-based modelterm·def on some larger or multi-task datasetsterm.

Computational consumption of different ML algorithms It is worthwhile mentioning that an optimal predictive model should have a good balance between prediction accuracy and computational efficiencyterm. As we all know, the run time complexity of SVMterm is quadratic to the number of training data [ 36 ]. As can be seen from Table 10 , it takes a few seconds (average wall-clock timeterm) to fit a model for the tasks whose data size is less than 4000. However, the average wall-clock timeterm is centupled when fitting the largest HIV dataset (data size of 40,748). That is to say, SVMterm is a good choice in dealing with small to medium datasets, but it will be frustrated when dealing with large datasets. To some extent, the same problem exists for the NN-based methods, which highly depend on the acceleration of graphics processing units (GPU) cards. However, XGBoostterm and RFterm provide a parallel tree training with high efficiency, and one of their strengths is the speed [ 40 ].

The mean wall-clock timeterm (seconds) for the six single-task datasetsterm given by the four descriptor-based and four graph-based modelsterm Here, we summarized the training speed of the four descriptor-based and four graph-based modelsterm on the six single-task datasetsterm (Table 10 ), and the training speed was evaluated by the mean wall-clock timeterm (seconds) from five independent runs where each run is to fit one corresponding model using the corresponding optimal hyper-parameters. It is worthwhile that the training speed of ML models can partly depend on the used hyper-parameters, such as the hidden layers of DNNterm, the trees of RFterm model and the graph convolutionterm layers of GNNterm model. In this study, the training speed of all the ML models were evaluated under the corresponding optimal hyper-parameters determined in the first stage of performance comparison. In addition, we shall emphasize that we are not analyzing the time and space complexity of different algorithms theoretically but intend to provide intuitive and touchable elapsed time of different algorithms under the affordable computational resources. All the compared algorithms were implemented by the recognized python packages (i.e., scikit-learnterm, PyTorchterm and PyTorch-based DGL), and more details can be accessed from the footnote of Table 10 . The choice of one-core, multi-cores or GPU largely depends on the inherent nature and common usage scenarios of algorithms, and what we try to present here is more likely a kind of rough users’ experience under the common usage scenarios, not the exactly CPU or GPU-time.

As shown in Table 10 , the training speed for the descriptor-based modelsterm is overwhelmingly faster than that of the graph-based modelsterm. For the three traditional descriptor-based modelsterm, only a few seconds were needed to finish the training of a model to most datasets. Among them, XGBoostterm and RFterm are the two most efficient algorithm and they are also able to manage big data with high proficiency. As expected, SVMterm performs efficiently on the relatively small datasets but its practicability will become much worse for large datasets. The descriptor-based DNNterm models show higher computability than GCNterm, GATterm, MPNNterm and Attentive FPterm, but all the NN-based models are highly dependent on GPU acceleration as mentioned above. Here, the top-performing graph-based algorithm, Attentive FPterm, demonstrates affordable computational efficiencyterm compared with its counterparts. Among the four graph-based modelsterm, the vanilla GCN modelterm is the most efficient algorithm and MPNNterm model is the worst one, which is in line with the common sense where the frameworks of vanilla GCN modelterm are much simpler than that of MPNNterm model. Actually, the total wall-clock timeterm including the hyper-parameter selection for each model was also analyzed but the conclusions are basically similar to the results discussed above (data are not shown).

Briefly, in terms of computational cost, the descriptor-based modelsterm are basically more efficient than the graph-based modelsterm. Among them, XGBoostterm and RFterm give the best computational efficiencyterm and it only needs a few seconds to train a model even for a large dataset. The descriptor-based DNNterm method is the most efficient one in its counterparts including GCNterm, GATterm, MPNNterm and Attentive FPterm, but the training of them largely depends on GPU acceleration.

The interpretation of XGBoostterm Model To check whether the learned knowledge from XGBoostterm is interpretable and reasonable, the SHAP method was used to analyze and interpret the developed models. Here, the XGBoostterm models for a regression dataset (ESOL) and a classification dataset (BBBP) were used as the examples. The top 20 representative molecular descriptorsterm and the corresponding SHAP valuesterm are presented in Fig. 2 .

Importance of the representative molecular descriptorsterm (the top 20) and the corresponding SHAP valuesterm given by XGBoostterm for the a ESOL and b BBBP datasets. One moleculeterm gets one dot on each descriptor’s line and dots stack up to show density ESOL: ESOL is a small regression dataset for aqueous solubility. As can be seen from Fig. 2 a, the most important descriptor given by the XGBoostterm model is h_logS, which represents the logarithm of aqueous solubility (mol/L). The feature value and SHAP valueterm in Fig. 2 a illustrate a clear positive correlation between the values of h_logS and the values of aqueous solubility, that means a higher h_logS will increase the aqueous solubility of a compound and vice versa, which is well in line with the expert knowledge. In Fig. 2 a, h_logD (the octanol/water distribution coefficient at pH = 7), which is related to the hydrophobicity of moleculesterm, is the second most important descriptor, and it presents a clear negative correlation with the value of aqueous solubility. This finding also well accords with the general phenomenon that higher hydrophobicity means lower solubility. In addition, the most significant parameter in the linear regression model for estimating the aqueous solubility of a compound developed by Delaney et al. is also a descriptor highly related to hydrophobicity (logP octanol ) [ 71 ]. Other two significant descriptors, including vsurf_D1 and vsurf_D7 that measure the hydrophobic volume of a moleculeterm, are highly related to hydrophobicity. Similar to h_logD, both of them have negative correlations with aqueous solubility, which is also well explainable where a higher hydrophobic volume will decrease the solubility of moleculesterm.

BBBP BBBP is a classification dataset for the blood–brain barrier (BBB) penetration of compounds. As we can see from Fig. 2 b, a number of the representative descriptors show clearly inverse correlations with BBB permeability, especially the descriptors a_don, SlopP_VSA2, h_ema and PubchemFP659 (2-(methylamino)ethan-1-ol substurcture), implying higher values of such descriptors will block moleculesterm to cross the BBB. Here, compared with the SHAP valueterm distributions of other descriptors, that of opr_leadlike (Oprea’s lead-like test) shows a huge difference due to the clear and successive blue dots on the left part of Fig. 2 b, indicating that opr_leadlike has positive correlations with BBB permeability. That’s to say, compounds with more lead-likeness would be more likely to cross the BBB. Here, most of those descriptors with inverse correlations with BBB permeability are polar-related descriptors, such as a_don (number of hydrogen bondterm donor atomsterm), h_ema (sum of hydrogen bondterm acceptor strengths) and PubchemFP659 (2-(methylamino)ethan-1-ol substurcture). This is consistent with the well-known fact that highly polar compounds have very low BBB permeation.

Virtual screeningterm profile analysis of different ML methods Many efforts have been dedicated to improving the prediction accuracy of different ML algorithms for molecular property prediction. In reality, these models can be served as VS tools to search for potential candidates from large chemical libraries and promote the discovery process. In our opinion, the efforts to improve the predictive accuracy and explore the VS profilesterm of different ML methods have the same priority because different ML models may offer quite different predictions in practical VS campaigns even they have similar predictive accuracy, which may directly determine what kinds of candidates are experimentally tested. To this end, a case study was conducted by identifying potential inhibitors towards HIV replication through the four descriptor-based and four graph-based modelsterm, and the small moleculeterm drugs deposited in DrugBankterm (Version: 5.1.5) were virtually screened by these models. All the explored models were developed based on the training set of the HIV dataset, optimized by the corresponding validation set and validated by the corresponding test set (the data folds were kept the same as those used in the first evaluation stage). The choice of this dataset was considered because of its relatively large data size and a more realistic proportion between inhibitors and noninhibitors. Prior to the screening, the polymers, inorganics, mixtures, salts were removed from the DrugBankterm small moleculeterm drug database. The duplicated compounds between the DrugBankterm database and the training set were also eliminated from the database. Finally, the remaining 1960 small moleculeterm drugs were used for screening. The output probabilityterm given by the optimal model was used as the score to measure the HIV replication inhibition abilityterm (Fig. 3 ). The higher the prediction score is, the greater the likelihood of being a HIV inhibitor is, and vice versa.

The distributions of the prediction scores for the 1960 screened moleculesterm predicted by the four descriptor-based modelsterm including a SVMterm, b XGBoostterm, c RFterm, d DNNterm and the four graph-based modelsterm including e GCNterm, f GATterm, g MPNNterm and h Attentive FPterm It can be observed that the distributions of the prediction scores for the 1960 moleculesterm given by the eight models vary from one to another although these models have similar prediction accuracy (Table 4 ). If an arbitrary threshold of 0.5 was used to classify inhibitors and non-inhibitors, the number of potential inhibitors given by the eight models are 7, 7, 45, 329, 86, 90, 158 and 284 respectively, highlighting the large difference of the predictions among different models. It seems that the conventional descriptor-based modelsterm (SVMterm, XGBoostterm and RFterm) are inclined to give more conservative predictions and the NN-based models are opposite. Among the eight models, SVMterm and XGBoostterm are the most two conservative models where only seven inhibitors were predicted by them. Inversely, the descriptor-based DNNterm model is the most radical one and about 17% compounds (329/1960) in the DrugBankterm database were predicted as inhibitors. Furthermore, the Euclidean distanceterm of the prediction score distributions for the 1960 drugs given by any two models was used to investigate the VS profileterm similarity of model pairs, and the lower this distance is, the more similar between the VS profilesterm of two models is, and vice versa (the minimum and maximum of this distance here are 0 and 44.27, respectively). As shown in Fig. 4 , with the exception of the SVMterm and XGBoostterm model pair, it is apparent that the Euclidean distancesterm of the prediction scores between any two of the eight models are relatively high, demonstrating that different ML models could perform very differently in practical VS campaigns.

The heat map of the Euclidean distancesterm of the prediction scores for different model pairs In order to uncover the structural features of the potential HIV inhibitors predicted by different ML models. The top 20 compounds with the highest scores given by each model were decomposed into different structural fragments and analyzed using Pipeline Pilot 2017. Three types of structural fragments were used, including Murcko Assemblies (contiguous ring systems plus chains that link two or more rings), Ring Assemblies (contiguous ring systems), and Bridge Assemblies (contiguous ring systems that share two or more bondsterm). The generated fragments were counted and the representative fragments whose counts are higher than or equal to four (not consider the common benzene component) for each model are shown in Fig. 5 (descriptor-based modelsterm) and Fig. 6 (graph-based modelsterm). As expected, the structural features of the potential inhibitors given by different models are highly diverse, demonstrating that different ML models are inclined to identify different sets of candidates and their diverse performances may be contributed from the different features used in training and the different principles of the algorithms. In addition, in the top 160 compounds given by the eight ML models (20 compounds for each model), 116 compounds are unique, and only a small fraction of compounds (10) were ranked in the top 20 in any three models, which also supported the aforementioned argument. Among the 10 compounds, it is pleasurable to observe that one compound used to combat HIV/AIDS, zidovudine, was predicted as a promising HIV inhibitor by all the eight models (Fig. 6 e). Here we found that the inhibitors predicted by the eight models share some nitrogen or oxygen heterocyclic components, four models including SVMterm, XGBoostterm, RFterm and GATterm have the tetrahydrofuran component in their predicted inhibitors and two models including GCNterm and MPNNterm have the tetrahydro-2H-pyran component in their predicted inhibitors. The structural features given by the SVMterm and GATterm models are highly overlapping. However, for all the eight ML models, no common structural component was found and the representative structural features given by the Attentive FPterm model show a high diversity. All in all, the structural features of the identified candidates by different ML models are diverse from each other.

The structural features of the potential inhibitors given by the four descriptor-based modelsterm including a SVMterm, b XGBoostterm, c RFterm and d DNNterm The structural features of the potential inhibitors predicted by the four graph-based modelsterm including a GCNterm, b GATterm, c (MPNNterm) and d Attentive FPterm; e the structure of the known HIV inhibitor identified by all the eight models Washing results of the benchmark datasets As described above, three washing steps were developed to automatically eliminate the incorrect or inappropriate structures from the original datasets. The washed datasets containing the original columns coupled with the canonical SMILESterm column were output as the final datasets. All of them are available in Additional file 1 and the detailed information of them are listed in Additional file 3 : Table S5. As shown in Additional file 3 : Table S5, several datasets, including BBBP, ClinTox, SIDER, Tox21, and ToxCastterm, contain relatively large numbers of incorrect or inappropriate structures (the ratio of the number of the removed compounds to its original number is large than 4%). In order to check the effect of the eliminated structures on model performance, two representative algorithms (i.e., XGBoostterm and Attentive FPterm) were used to build the prediction models for the washed datasets of BBBP, Tox21, ToxCastterm, and SIDER. The same hyper-parameters described above were used in model building. Similarly, the models were validated by 50 times independent runs and the statistical results are listed in Additional file 3 : Table S6. It can be observed that the predictions of the models to the washed datasets do not show large difference compared with those to the original datasets. The predictions to the washed datasets of BBBP become slightly better for both models, while those to the washed datasets of ToxCastterm and SIDER become slightly worse for both models. And the predictions to the washed datasets of Tox21 get slightly better for XGBoostterm and slightly worse for Attentive FPterm. However, it should be noted that our purpose is not highlighting the impact of incorrect or inappropriate structures on the predictive accuracy of models but merely points out that the quality of the public datasets should be carefully checked.

Conclusion GNNterm has gained great interest in molecular property prediction due to its ability to learn molecular representationsterm·var automatically. It appears that most studies reported so far have drawn the conclusion that GNNterm is more promising than traditional descriptor-based modelsterm. In this study, we demonstrated that on average the descriptor-based modelsterm outperform the graph-based modelsterm in the predictions of a variety of molecular properties in terms of predictive accuracy and computational efficiencyterm. SVMterm generally gives the best predictions to regression tasksterm. Both XGBoostterm and RFterm can give reliable predictions to classification tasksterm, and graph-based methods, such as GCNterm and Attentive FPterm, can offer outstanding performance on a fraction of larger or multi-task datasetsterm. In terms of computational efficiencyterm, XGBoostterm and RFterm have fast computability and only need a few seconds to train a model even for a large dataset. Moreover, descriptor-based modelterm·def can be well interpreted by the SHAP method. Finally, the ML models were used to conduct a VS study towards HIV, and the results demonstrate that different ML algorithms offer diverse VS profilesterm. In conclusion, our study illustrates that the descriptor-based modelsterm are able to achieve better or comparable predictions to the highly-intricate and specialized graph-based modelsterm in terms of prediction accuracy, computability and interpretability.

Supplementary information Additional file 1. The datasets for the single tasks.

Additional file 2. The python source codes that implement the workflow.

Additional file 3: Table S1. The performance of the top three runs and the worst three runs among the 50 times independent runs given by the XGBoostterm model for the 11 datasets. Table S2 . The performance comparison (MAEterm metric) of the four descriptor-based and four graph-based modelsterm on the three regression datasets. Table S3 . The performance comparison (R2 metric) of the four descriptor-based and four graph-based modelsterm on the three regression datasets. Table S4: The performance comparison (average RMSEterm) of the 50 times independent runs on three regression datasets including ESOL, FreeSolv, and Lipop before/after removing the top three related descriptors given by the four descriptor-based modelsterm (SVMterm, XGBoostterm, RFterm and DNNterm). All the models named with suffix ‘1’ refer to the models developed based on the remaining descriptors. Table S5. The detailed information for the 11 washed datasets. Table S6. The performance comparison of the 50 times independent runs on four datasets including BBBP, Tox21, ToxCastterm, and SIDER before/after washing for the XGBoostterm and Attentive FPterm models.

Acknowledgements Not applicable.

Authors’ contributions DJ and ZW developed the prediction models, analyzed the data, wrote the manuscript; CH, GC, BL, ZW and CS evaluated and interpreted the results, wrote the manuscript; DC, JW and TH conceived and supervised the project, interpreted the results, wrote the manuscript. All authors read and approved the final manuscript.

Funding This work was financially supported by Key R&D Program of Zhejiang Province (2020C03010), National Natural Science Foundation of China (21575128, 81773632), and Natural Science Foundation of Zhejiang Province (LZ19H300001).

Competing interests The authors declare that they have no competing interests.

Footnotes Contributor Information Dongsheng CaopersName, Email: oriental-cds@163.com.

JianpersName Wu, Email: wujian2000@zju.edu.cn.

Tingjun HoupersName, Email: tingjunhou@zju.edu.cn.

Supplementary information Supplementary information accompanies this paper at 10.1186/s13321-020-00479-8.

References Associated Data Supplementary Materials Additional file 1. The datasets for the single tasks.

Additional file 2. The python source codes that implement the workflow.

Additional file 3: Table S1. The performance of the top three runs and the worst three runs among the 50 times independent runs given by the XGBoostterm model for the 11 datasets. Table S2 . The performance comparison (MAEterm metric) of the four descriptor-based and four graph-based modelsterm on the three regression datasets. Table S3 . The performance comparison (R2 metric) of the four descriptor-based and four graph-based modelsterm on the three regression datasets. Table S4: The performance comparison (average RMSEterm) of the 50 times independent runs on three regression datasets including ESOL, FreeSolv, and Lipop before/after removing the top three related descriptors given by the four descriptor-based modelsterm (SVMterm, XGBoostterm, RFterm and DNNterm). All the models named with suffix ‘1’ refer to the models developed based on the remaining descriptors. Table S5. The detailed information for the 11 washed datasets. Table S6. The performance comparison of the 50 times independent runs on four datasets including BBBP, Tox21, ToxCastterm, and SIDER before/after washing for the XGBoostterm and Attentive FPterm models.