Computational Prediction of Site of Metabolism for UGT-catalyzed

Dec 26, 2018 - Yingchun Cai , Hongbin Yang , Weihua Li , Guixia Liu , Philip W Lee , and Yun Tang. J. Chem. Inf. Model. , Just Accepted Manuscript...
2 downloads 0 Views 1MB Size
Subscriber access provided by YORK UNIV

Chemical Information

Computational Prediction of Site of Metabolism for UGT-catalyzed Reactions Yingchun Cai, Hongbin Yang, Weihua Li, Guixia Liu, Philip W Lee, and Yun Tang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00851 • Publication Date (Web): 26 Dec 2018 Downloaded from http://pubs.acs.org on December 27, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Computational Prediction of Site of Metabolism for UGT-catalyzed Reactions

Yingchun Cai, Hongbin Yang, Weihua Li, Guixia Liu, Philip W. Lee, Yun Tang*

Shanghai Key Laboratory of New Drug Design, School of Pharmacy, East China University of Science and Technology, Shanghai 200237, China

*Corresponding author, E-mail: [email protected]

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The investigation of metabolically liable sites of xenobiotics mediated by UDP-glucuronosyltransferases (UGT) is important for lead optimization in early drug discovery. However, it is time-consuming and costly to identify potential susceptible sites experimentally. In silico approaches are hence developed to predict site of metabolism (SOM) of UGT-catalyzed substrates. In the present work, four major types of reactions catalyzed by UGT were collected from the book “Handbook of Metabolic Pathways of Xenobiotics” along with their corresponding glucuronidation products. These observed and non-observed SOMs were identified and encoded by atom environment fingerprint. Four resampling methods were performed to treat with data imbalance, while four feature selection methods were utilized to choose appropriate features. Three tree-form machine learning algorithms were conducted to build discriminating models, and optimal models were then obtained for each type of reaction. Results indicated that all chosen best models showed AUC values ranging from 0.713 to 0.869 for two independent external validation sets. Our study could supply valuable information for optimization of pharmacokinetic profile and contribute to metabolism prediction.

ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Introduction Metabolism is one of the major clearance pathways for about 75% of all drugs.1 It produces metabolites with substantially changed physicochemical, pharmacological and toxicological profiles. These properties could affect activation, deactivation, toxification and detoxification, then alter efficiency and safety of drugs. There are two phases of drug metabolism, namely Phase I - modification and Phase II - conjugation. Glucuronidation, catalyzed by UDP-glucuronosyltransferases (UGT), is the main reaction of Phase II. UGT transfers glucuronic acid from UDP-glucuronic acid to diverse substrates, such as alcohols, carboxylic acids, amines, and so on (Table 1). The human genome encodes 21 UGTs, which can be classified into four subfamilies, named UGT1, UGT2, UGT3, and UGT8. These isoforms are widely distributed in liver, intestine, kidney, lung, and skin.2 UGTs mediate the conjugation of UDP-glucuronic acid with lipophilic chemicals that contain a suitable acceptor reaction group (-OH, -COOH, and -NRx), following a serine hydrolase-like catalytic mechanism.3 Actually, it is based on a second-order nucleophilic substitution mechanism.4 The resulted glucuronide is more hydrophilic and easier to eliminate than the original compound. Both endogenous and xenobiotic chemicals can be metabolized by UGTs. Many isoforms of UGTs are capable of conjugating a variety of structurally unrelated substrates and possess a partial overlap in substrate specificity.5 However, different UGTs can conjugate the same substrate at different rates. Substances containing more than one nucleophilic group may be catalyzed by different UGTs.6, 7

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Glucuronidation not only provides an essential detoxification way for xenobiotics, but also result in short duration of therapy and loss of pharmacological activity for certain drugs. For that reason, medicinal chemists need optimize the pharmacokinetic properties of lead compounds to obtain the expected therapeutic efficacy. For example, chemical group substituents8 or bioisosteric modifications9 were reported to overcome the problem of rapid glucuronidation and design more metabolically stable drugs. In contrast, Millan’s group designed inhaled p38 inhibitors prone to be metabolized by UGTs for minimizing systemic exposure to reduce adverse events.10 Therefore, it is vital to understand UGT metabolism in lead optimization, and it is urgent to identify which site to undergo glucuronidation in a molecule. In general, experimental identification of site of metabolism (SOM) in one UGT-catalyzed reaction is cost-intensive and time-consuming. Hence, there is an increasing interest to develop reliable computational methods for prediction and understanding of the UGT-mediated reactions. In an early study, Sorich et al. exhibited classification models for rapid prediction of chemical metabolism by human UGT isoforms based on quantum chemical descriptors.11 Although 84% accuracy in prediction was obtained, they only focused on reaction phenotypes and the models could only discriminate substrates and non-substrates. Afterwards, they first reported in silico models to predict glucuronidation sites by analyzing each functional group metabolized by UGTs.12 They figured out the varied susceptibility of different types of nucleophilic reactions by calculation of frequency distribution in local chemical structures. Ghemtio et al. developed a novel approach to combine SVM classification

ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

and CoMSIA together on UGT1A6 interacting molecules, and obtained satisfactory results.13 However, this method needs a high-quality 3D protein structure and cannot extend to those without 3D structures. SOM-UGT described by Peng et al. could perform prediction on four types of glucuronidation sites individually14 by neglecting less common or atypical substructures. However, SOM-UGT did not combine their four independent models into a single one to evaluate a given compound. SOMP uses a naïve Bayes classifier to rank all sites in one compound by their UGT-mediated likelihood.15 A major insufficiency of Top-N metric used by SOMP is that it did not consider the bias of SOM prediction for some molecules.16 The solo metric IAP (AUC) in SOMP might not thoroughly judge the models. The latest strategy XenoSite17 presents two models to predict UGT metabolism with higher accuracy than SOM-UGT and SOMP, but it requires calculation of quantum chemical-derived descriptors, which is technically complex and difficult to select proper features to avoid over-fitting of models. In the present study, we developed SOM prediction models for four types of UGT-mediated reactions (AlOH, ArOH, COOH, and Nitrogen) using machine learning methods. Each retrieved reaction was described by atom environment fingerprint and labeled with reaction type by comparing its reactant and product structures. Before model building, resampling strategies were utilized to balance the data. Feature selection methods were also employed. Two independent datasets were collected to validate the performance of models. Optimal models were selected on the basis of classification performance of different models.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Materials and Methods Data Collection and Preparation All the data were retrieved from Lee’s book “Handbook of Metabolic Pathways of Xenobiotics” along with their corresponding glucuronidation products in human body.18 Each drug and its product were redrawn and restored in SMILES format. The observed UGT-catalyzed SOMs were identified by analyzing the structures of reactants and products. All the glucuronidated sites in each molecule were categorized into one of the four substructure groups: aliphatic hydroxyl groups (AlOH), aromatic hydroxyl groups (ArOH), carboxyl groups (COOH) and amino groups (Nitrogen). Here, we only considered these four types of reactions as listed in the first four lines of Table 1. Among them, the nitrogen one can be aromatic and aliphatic amino/ heterocyclic nitrogens. Each SOM was regarded as positive one (labeled as ‘1’) if it was experimentally observed, otherwise it was tagged as negative one (labeled as ‘0’). Finally, each of the clean four types of SOM datasets was randomly divided into training set and test set in the ratio of 8 : 2.

Table 1. Examples of UGT-catalyzed reactions. The first four ones were used in this study. Reaction Types

Reaction Group

Description*

ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

O

O N

AlOH

N OH

N

O O

N

O

OH

HO

OH OH

Antipyrine

OH

OH OH

O-Glucuronidation ArOH

N H

N H

OH

O

OH

O

O

HO

OH

OH

OH

Albuterol O

O

O O

OH

COOH

O

OH

O HO

O

OH OH

Aniracetam

O N+

N-Glucuronidation

N

Nitrogen

N

N

O

O-

HO

OH OH

Nicotine N

S-Glucuronidation

OH N

N

O N H

Sulfur

HO O

N

S

N

S

OH

O

Tanaproget

N N

N N

C-Glucuronidation

Carbon

O

OH

O

O

S O

O S HO O HO

O O

O OH

OH

Sulfinpyrazone * Substructure in black circle is experimentally regarded as site of metabolism.

For performance validation of classification models, we collected two more datasets as external validation sets. One was merged by the test sets from Peng et al.14

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and Rudik et al.15, respectively. Another one was assembled from KEGG database,19, 20

which contained amount of available metabolism resources. We termed the two

datasets as Ext1 and Ext2. Atom Environment Fingerprints Atom environment fingerprint was derived from circular fingerprint,21 which was generated by Morgan algorithm on a set of atom invariants. For a given atom in one molecule, its environment feature vectors was described as the occurrence counts of Daylight atomic invariants at different topological distances from questionable atom. The size of atom environment fingerprint varied based on different bond depths, so the radius should be provided when calculating the fingerprint. In this work, we defined an SOM as an atom where glucuronidation occurred in a metabolic reaction. The SOM was described by various size of atom environment fingerprint. All fingerprints were generated using the open-source toolkit RDKit with size = 1024 and radius = 3.22 Feature Selection Feature selection is a necessary step in machine learning workflow. The aim is to discard those redundant, noisy, or irrelevant features to improve prediction accuracy of models.23 Four feature selection methods were used here, namely variance threshold (VT), select percentile of feature (SPF), recursive feature elimination (RFE), and principle component analysis (PCA). Detailed instruments of these methods can be found at URL: http://scikit-learn.org/stable/user_guide.html. VT was first used to remove those whose values were one or zero in more than

ACS Paragon Plus Environment

Page 8 of 29

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

80% of the samples with 0.8 * 0.2 as the variance threshold. SPF was then used to keep those at the highest scores according to user-defined percentile. ANOVA F-value between label array and feature array was chosen as the score function. Here, label arrays are one-dimensional data, while feature arrays are two-dimensional data. RFE recursively used estimator to train on the smaller and smaller sets of features based on the importance of each feature. The least important features were pruned from current set of features by scores. In another way, PCA was used to compress initial sets of features after VT filter. For assessing the impact of parameters of the four methods on classification performance, we defined different percentiles and components for SPF and PCA approaches, respectively. Data Resampling Data resampling aims to balance the numbers between positives and negatives. For that purpose, three under-sampling methods and one weight method were used here. The weight method imposed different weights on two classes of samples, where a larger weight was set for the majority and a smaller one was set for the minority. The three under-sampling methods included random under-sampling (RU), edited nearest neighbor (EN), and similarity-based sampling (SS). Detailed descriptions of RU and EN see: http://imbalanced-learn.org/en/stable/index.html. SS was a method that we developed to reduce the number of the majority based on cosine similarity. A threshold of 0.8 was set to remove samples. All the resampling operations were only performed on training set. Model Building

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Three machine learning methods were used to build the models for SOM prediction, namely decision tree (DT), random forest (RF) and AdaBoost (AB). All the methods were employed in an open-source machine-learning toolkit scikit-learn (version 0.18).24 The whole workflow was presented in Figure 1. Both local and global models were developed. A local model was built for one of the four types of reactions, whereas a global model was constructed for all the four types of reactions. The initial parameters and a grid of parameter values we used for each classifier were listed in Table S1.

ACS Paragon Plus Environment

Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 1. The whole workflow employed in this work.

Decision Tree. DT is a flowchart-like structure where each internal node represents a “test” on an attribute, each branch represents the outcome of the test, and each leaf node represents a label.25 DT uses non-parametric supervised learning algorithm for decision and predicts the class by learning simple decision rules inferred from the data features. DT is simple to understand and interpret. The constructed trees

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

can be visualized for analysis. Random Forest. RF is an ensemble learning algorithm for classification, which constructs a multitude of DTs and outputs the label.26 In RF classification task, a small set of input features was selected at random to build each tree in an ensemble. Because of the randomness, the bias of the forest usually slightly increases. However, the averaging procedure results in the decrease of variance, to compensate for the bias increase. Hence, RF method can produce an overall better classification performance. AdaBoost. AB classifier is an ensemble method and a meta-estimator. It first fits a base classifier on the original dataset, and trains additional copies of the classifier on the same dataset. The weights of incorrectly predicted instances are then adjusted so that subsequent classifiers concentrate more on difficult cases.27 The predictions from all of them are combined through a weighted majority vote to produce the final result. Validation of Model Performance In our SOM classification task, classification accuracy (ACC), sensitivity (SE), specificity (SP), and area under the curve (AUC) were employed to evaluate the performance of a classifier with 10-fold cross validation. The first three metrics were calculated by the following equations, respectively. ACC =

𝑇𝑃 + 𝑇𝑁

(1)

𝑇𝑃 + 𝑇𝑁 + 𝐹𝑃 + 𝐹𝑁 𝑇𝑃

SE = 𝑇𝑃 + 𝐹𝑁 𝑇𝑁

SP = 𝑇𝑁 + 𝐹𝑃

(2) (3)

AUC value was the area under the receiver operating characteristic curve (ROC) that could provide further information about the performance of the models. By

ACS Paragon Plus Environment

Page 12 of 29

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

iteratively setting different thresholds, ROC was plotted to show the separation ability of a binary classifier.

Results In this work, we collected data from unique sources first. Two series of models were then built to identify positive SOMs from negative ones: four local models (each for one reaction type) and one global model. All the models were assessed by two external validation sets. The optimal models were selected for each type of reaction by comparing metrics among these models. Data Set Analysis In total 400 drugs metabolized by UGT were retrieved from Lee’s book. After data processing, each SOM was collected into our dataset for model building. As shown in Table 2, for Lee’s AlOH dataset, the numbers of SOM+ and SOM- were 93 and 82, respectively. However, as to the other two datasets (ArOH and COOH), positives were much more than negatives, whereas the opposite status was found in Nitrogen dataset. For Ext1 dataset, the four types of reactions contain dozens of SOM+ data, but only a few SOM- data except Nitrogen (126 SOM- data). The total numbers of compounds for Lee, Ext1 and Ext2 were 400, 74 and 112, respectively. Analyzing these three datasets, we realized that each dataset encountered imbalance problem and the maximum ratio of positive to negative reached to 20 and minimum was 0.09. This issue indicated that resampling methods should be employed first.

Table 2. Statistics information for three datasets of SOMs

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Reactions AlOH ArOH COOH Nitrogen Total

Lee

Page 14 of 29

Ext1

Ext2

SOM+

SOM-

NO.C*

SOM+

SOM-

NO.C*

SOM+

SOM-

NO.C*

93 132 84 55 364

82 18 23 545 668

133 131 100 271 400

30 30 20 25 105

7 6 1 126 140

29 26 21 60 74

37 39 14 5 95

21 44 1 53 119

42 43 12 32 112

*NO.C indicated the number of substrates changing in reactions.

Performance of Local Models After data resampling and feature selection, baseline models were built by machine learning methods without parameter optimization, to provide a general view on prediction ability of models for the four types of reactions. All the metrics of the baseline models were averaged with standard deviation. Herein, we took AlOH as an example to exhibit these results as shown in Figure S1 for training and test sets, and in Figure S2 for ext1 and ext2. Other baseline models for rest three types of reactions were shown in Figures S3-S8. From Figure S1, we obviously found that AB performed the worst, followed by DT and RF. Its ACC value was lower in common for training set and test set. A similar observation was found in Figure S2 for ext1 and ext2. As to the three resampling methods, EN was the attractive one compared to the others for training set. However, in test set, the performance of EN was a little worse than RU and SS. As to the two external validation sets, resampling methods could perform better than non-sampling ones generally. For two feature selection methods, SPF exhibited better results than PCA on test set as well as ext1 and ext2 (Figure S2) although its performance on training set was a little worse.

ACS Paragon Plus Environment

Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

According to these baseline models, we chose four better ones (each for one reaction) for parameter tuning to get the optimal models. As shown in Table 3, the ACC values of training set ranged from 0.782 to 0.892, and the AUC values were from 0.870 to 0.965, respectively, among which AlOH reaction was the worst. For test set, the ACC and AUC values were just a little lower than the corresponding ones in the training set, but still very good. Detailed information about these optimal models was shown in the Method Combination column of Table 3. Except for AlOH reaction, the other three models were built by DT methods. These results indicated that our models were effective for identification of positive SOMs from negative SOMs.

Table 3. Performance of optimal models for training and test sets and method combinations Training Set

Reactions AlOH ArOH COOH Nitrogen

Method Combination*

Test Set

ACC

SE

SP

AUC

ACC

SE

SP

AUC

0.782 0.892 0.882 0.861

0.760 0.888 0.941 0.913

0.790 0.923 0.647 0.854

0.870 0.951 0.920 0.965

0.771 0.867 0.818 0.817

0.867 0.880 0.875 0.889

0.700 0.800 0.667 0.811

0.820 0.928 0.865 0.932

RF_SS_40_False DT_0.8_35_True DT_0.7_50 DT_EN_10_True

* The method combination a_b_c_d or e_f_g is described as following: a or e is machine learning method: RF, DT or AB; b or f is data resampling method: RU, SS or EN, a float means weights for negative class; c is percentile of unselected features for SPF method; d: whether to use RFE selection algorithm or not; g is component for PCA.

In order to show our models in detail, we exhibited the four feature selection procedures for these models in Figures S11-S14. As shown in Figure S11, VT could effectively filter low variance feature vectors and keep information-full ones to be filtered further by SPF or RFE. As we expected, these two methods (Figures S12 and

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

S13) employed more advanced algorithms to hold on useful features and reduce dimension space. In Figure S12, we could see that there were 77, 34 and 46 features selected for AlOH, ArOH and Nitrogen, respectively, based on the SPF method. RFE was conducted to further filter useless features (Figure S13). Finally, 20 and 25 features were maintained for ArOH and Nitrogen, separately. In Figure S14, we could see that the percentage of variance explained by PCA method increased along with the component for COOH. To further ensure the performance of our optimal models, Ext1 and Ext2 were utilized. ACC and AUC values for these models were shown in Figure 2. ACC values were from 0.622 to 0.762, while AUC values ranged from 0.713 to 0.869. Ext1 consisted of two test sets from Peng’s work14 and Rudik’s work15. Here, we also did a prediction on Peng’s test set alone by our optimal models. The results were shown in Table 4, where we could see that ACC values ranged from 0.636 to 0.798, and SE values ranged from 0.600 to 0.688. As to SP values, 1 for AlOH, ArOH and COOH, while 0.816 for Nitrogen. These results demonstrated that our optimal models owned satisfactory prediction capability.

ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 2. Model assessment results for two external validation sets. ACC and AUC were shown to illustrate the performance of these models.

Table 4. Performance of our local models on Peng’s test set Reactions AlOH ArOH COOH Nitrogen

TP 14 15 6 11

FN 8 7 4 5

TN 2 2 1 84

FP 0 0 0 19

ACC

SE

SP

0.667 0.708 0.636 0.798

0.636 0.682 0.600 0.688

1 1 1 0.816

Some examples were shown in Figure 3 to demonstrate the performance of our local models. For example, fenoldopam28 is a dopamine D1 agonist that is used as an antihypertensive agent. It could lower blood pressure through arteriolar vasodilation. Its elimination is largely done by conjugation, and glucuronidation is one of the major routes of conjugation. By our models, three out of the four potential sites were correctly

identified.

Another

example

is

atorvastatin29:

an

inhibitor

of

3-hydroxymethy-glutary coenzyme A reductase. Glucuronidation of atorvastatin by

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

UGT helps in lactonization reaction and is associated with toxicity. By our models, three experimentally not observed sites and one experimentally observed site were correctly predicted. These results could provide valuable information for structural optimization to improve the pharmacokinetic profile of drugs.

ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 3. Some examples of predicted SOMs within several molecules.

Performance of Global Models Besides the local models, a global model was also constructed to identify SOMs from non-reactive sites. Similarly, many baseline models were first built (Figures S9 and S10), and the best one was selected for parameter optimization. The results were shown in Figure 4, where ACC, SE, SP, and AUC for the training set were all greater than 0.90. As to test set, except SE = 0.708, the others were over 0.80. ACC and AUC were 0.804 and 0.850 for Ext1, 0.649 and 0.682 for Ext2, respectively. As to Peng’s data set, ACC, SP and AUC were higher than 0.85, while SE was lower than 0.7. According to the above analysis, the optimal global model could effectively identify SOMs of UGT-catalyzed substrates, although the validation results were a little low for Ext2 with SP = 0.535 but SE = 0.810.

Figure 4. Performance of the global model. Method combination includes AdaBoost classifier

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

without resampling operation and 255 components for PCA method.

The combination of global model with local models was used to improve the predictive capability of UGT-mediated SOMs. The combination means that both local and global models were used to predict a certain site and marked by a label, separately. Then the two labels were taken together to execute exclusive OR, which led to the final prediction. The results were shown in the bottom side of Table 5. From Tables 4 and 5, we could see that combination of global model with local models could effectively help to improve performance of local models, especially for models of ArOH and COOH. Table 5. Comparison of our combination results with Peng’s results ACC SE SP TP FN TN FP AlOH 17 4 2 1 0.792 0.810 0.667 ArOH 18 5 1 1 0.760 0.783 0.500 Peng’s Results a COOH 10 0 1 0 1.000 1.000 1.000 Nitrogen 12 2 68 6 0.909 0.857 0.919 AlOH 15 6 2 1 0.708 0.714 0.667 ArOH 21 1 2 0 0.958 0.955 1.000 Our b Results COOH 10 0 1 0 1.000 1.000 1.000 Nitrogen 11 5 84 19 0.798 0.688 0.816 a Peng’s test set results calculated from TP, FP, TN, and FN supplied in the paper. b Our corresponding results by combination of local and global models for Peng’s 54 molecules.

Comparison with Others To compare the performance of our models with others, we did a statistic analysis on the 54 molecules from Peng’s test set14. The results were also shown in the left side of Table 5. As we can see in this table, ACC and SE were slightly lower than Peng’s for AlOH and Nitrogen. However, SP was all equal to or higher than

ACS Paragon Plus Environment

Page 20 of 29

Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Peng’s except Nitrogen. Moreover, for COOH, the prediction was extremely satisfactory (ACC = 1, SE = 1, and SP = 1). Except for the above comparison, feature selection was our advantage over Peng’s paper, which could reduce the models’ complexity. Moreover, tree-form machine learning methods were easily interpretable than SVM method. This comparison demonstrated that our models had satisfactory classification capability and could be comparable to other models.

Discussion In this work, SOM prediction models were established for four types of UGT-mediated reactions with three machine learning methods and atom environment fingerprint. The data used in this work were brand new, and provided important recognition patterns for model building. Four data resampling strategies were applied to reduce the gap between two classes for instance balance in order to make the prediction results more scientific. Before model building, feature selection approaches were also performed on calculated site fingerprints to avoid information overlapping of bits, to improve the robustness of models, and to decrease the complexity of models. Then the optimal classification model for each type of reaction was chosen according to four assessment metrics from amounts of trained models. Here, tree-form machine learning methods (DT, RF, and AdaBoost) were chosen according to interpretation. As we all know, DT can be highly interpretable and simple, because it involves a series of splits based on logical decisions, starting from the most important top-level node. DT visually looks like an upside down tree. As to RF and AdaBoost,

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

they are ensemble of DT. Thus these two methods can also be interpretable. A little limitation was that not very large dataset (Table 2) impeded us adequately exploring the capacity of machine learning methods. In machine learning tasks, feature selection is used to retain a set of special representation of original feature vectors, which is extremely important for model building. This operation aims at choosing those ones strongly related to classes with weak correlation between each other. In real world, however, the retrieved data are usually sophisticated, redundant and varied. These unpopular properties of datum lead to difficulties for human to deal with big data. It is necessary for us to find out useful characteristics from original data. However, manual selection of features depends on amounts of human labor and professional knowledge, and is hard to extend. So in this paper, we employed special techniques to learn and extract features, making our classification task speedier and more efficient. As exhibited above, four feature selection methods (VT, SPF, RFE, and PCA) were used to discard useless features and to improve performance of prediction models. Actually, PCA method does not belong to feature selection technology but a dimensionality reduction strategy. It still can be used to condense features and is extendedly utilized in dataset transformation. As we can see in Figure 5, the best models for AlOH, ArOH and Nitrogen were built using SPF method with 40%, 35% and 10% percentile, respectively. RFE method was also applied for ArOH and Nitrogen models, except for AlOH one. As to COOH model, PCA approach generated the optimal model, whose component was set to 50. Meanwhile, best global model

ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(Figure 6) was also obtained by PCA with component of 255. Based on above discussion,

we

concluded

that

feature

selection

well-performance model building.

ACS Paragon Plus Environment

indeed

contributed

to

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Different percentiles (5% to 100% with 5% as step), components with 5 as step, and class_weight for negative SOM (0.1-0.9 with 0.1 as step and ‘balanced’ as default) were employed to train local models. AUC was used to assessment and compare all these models which were generated under 1024bit fingerprint and radius of three.

Figure 6. Different components (for ‘coarse’ with 50 as step and for ‘detail’ with 5 as step), and class_weight for negative SOM (0.1-0.9 with 0.1 as step and ‘balanced’ as default) were employed to train global models. AUC was used to assessment and compare all these models. All the models were generated under 1024 bit fingerprint and radius of three.

As to SOM prediction, we usually encounter such an issue that the number of positive instances is always much lower than that of negative ones. In such situations, conventional machine learning methods treat with these classes of samples equally,

ACS Paragon Plus Environment

Page 24 of 29

Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

which could result in low recognition accuracy. Moreover, the subsequent prediction can be affected by the imbalance issue. In present study, four resampling approaches (RU, EN, SS, and Weight) were utilized to balance binary classes of datasets and to reduce the effect of imbalance issue on accuracy of models. We illustrated the effect of a classifier with different level of class balance. As indicated in Table 3 and Figure 5, SS and EN methods generated the best models for AlOH and Nitrogen, retaining default ‘balanced’ of class_weight level for Weight method. In Figure 5, for Ext1 we could see that most of the nine imbalanced levels (0.1-0.9) produced lower AUC values than that of ‘balanced’ level. But for Ext2 of Nitrogen data, the ‘balanced’ level didn’t give optimal AUC value. Here we need demonstrate that the best model for each type of reaction was selected based on three validation sets (one internal test set and two external validation set), not based on one validation set. For ArOH and COOH in Figure 5, we found that class_weight of 0.8 and 0.7 of Weight method produced the best models. As to global model in Figure 6, the ‘balanced’ level generated the highest AUC value. These results proved that our resampling methods improved performance of discrimination models effectively.

Conclusions In this work, we developed a generalized workflow to predict SOMs of four types of UGT-catalyzed reactions (AlOH, ArOH, COOH, and Nitrogen) with tree-form classification algorithms. New data extracted from Lee’s book were employed to build the models. Atom environment fingerprint was used to generate

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

feature vectors of each SOM. Four resampling methods were employed to reduce the number gap between positives and negatives, while four feature selection methods were utilized to discard redundant features. The optimal model was obtained for each type of reaction in terms of assessment indexes. Compared with other published studies, our workflow and data are novel, and the predictive accuracy is satisfactory. The employment of resampling, feature selection and tree-from algorithms make these models interpretable, low complexity, and scalability. Moreover, combination of local and global models provides more accurate prediction. A big limitation of this study is that the dataset is small. Hence, what we need to do next is to collect more substrates to enlarge our dataset. This study would add values to further direction and make rational decisions for UGT-mediated lead optimization.

Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant 81872800) and the 111 Project (Grant B07023).

Supporting Information Supporting Information Available: Tables S1, initial parameters and a grid of parameter values; Figures S1-S10, baseline models for the four types of reactions; Figures

S11-S14, feature selection score.

ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

References 1.

Kirchmair, J.; Goller, A. H.; Lang, D.; Kunze, J.; Testa, B.; Wilson, I. D.; Glen, R. C.; Schneider,

G., Predicting Drug Metabolism: Experiment and/or Computation? Nat. Rev. Drug Discovery 2015, 14, 387-404. 2.

Rowland, A.; Miners, J. O.; Mackenzie, P. I., The UDP-Glucuronosyltransferases: Their Role in

Drug Metabolism and Detoxification. Int. J. Biochem. Cell Biol. 2013, 45, 1121-1132. 3.

King, C. D.; Rios, G. R.; Green, M. D.; Tephly, T. R., UDP-Glucuronosyltransferases. Curr. Drug

Metab. 2000, 1, 143-161. 4.

Radominska-Pandya, A.; Czernik, P. J.; Little, J. M.; Battaglia, E.; Mackenzie, P. I.,

STRUCTURAL AND FUNCTIONAL STUDIES OF UDP-GLUCURONOSYLTRANSFERASES*. Drug Metab. Rev. 1999, 31, 817-899. 5.

Magdalou, J.; Fournel-Gigleux, S.; Ouzzine, M., Insights on Membrane Topology and

Structure/Function of UDP-Glucuronosyltransferases. Drug Metab. Rev. 2010, 42, 159-166. 6.

Zhurova, E. A.; Zhurov, V. V.; Chopra, D.; Stash, A. I.; Pinkerton, A. A., 17α-Estradiol·1/2 H2O:

Super-Structural Ordering, Electronic Properties, Chemical Bonding, and Biological Activity in Comparison with Other Estrogens. J. Am. Chem. Soc. 2009, 131, 17260-17269. 7.

Itäaho, K.; Mackenzie, P. I.; Ikushiro, S.-i.; Miners, J. O.; Finel, M., The Configuration of the

17-Hydroxy Group Variably Influences the Glucuronidation of β-Estradiol and Epiestradiol by Human UDP-Glucuronosyltransferases. Drug Metab. Dispos. 2008, 36, 2307-2315. 8.

Tanaka, R.; Rubio, A.; Harn, N. K.; Gernert, D.; Grese, T. A.; Eishima, J.; Hara, M.; Yoda, N.;

Ohashi, R.; Kuwabara, T.; Soga, S.; Akinaga, S.; Nara, S.; Kanda, Y., Design and Synthesis of Piperidine Farnesyltransferase Inhibitors with Reduced Glucuronidation Potential. Bioorg. Med. Chem. 2007, 15, 1363-1382. 9. M.;

Kawada, H.; Ebiike, H.; Tsukazaki, M.; Nakamura, M.; Morikami, K.; Yoshinari, K.; Yoshida, Ogawa,

K.;

Shimma,

N.;

Tsukuda,

T.;

Ohwada,

J.,

Lead

Optimization

of

a

Dihydropyrrolopyrimidine Inhibitor Against Phosphoinositide 3-Kinase (PI3K) to Improve the Phenol Glucuronic Acid Conjugation. Bioorg. Med. Chem. Lett. 2013, 23, 673-678. 10. Millan, D. S.; Bunnage, M. E.; Burrows, J. L.; Butcher, K. J.; Dodd, P. G.; Evans, T. J.; Fairman, D. A.; Hughes, S. J.; Kilty, I. C.; Lemaitre, A.; Lewthwaite, R. A.; Mahnke, A.; Mathias, J. P.; Philip, J.; Smith, R. T.; Stefaniak, M. H.; Yeadon, M.; Phillips, C., Design and Synthesis of Inhaled p38 Inhibitors for the Treatment of Chronic Obstructive Pulmonary Disease. J. Med. Chem. 2011, 54, 7797-7814. 11. Sorich, M. J.; McKinnon, R. A.; Miners, J. O.; Winkler, D. A.; Smith, P. A., Rapid Prediction of Chemical Metabolism by Human UDP-glucuronosyltransferase Isoforms Using Quantum Chemical Descriptors Derived with the Electronegativity Equalization Method. J. Med. Chem. 2004, 47, 5311-5317. 12. Sorich, M. J.; McKinnon, R. A.; Miners, J. O.; Smith, P. A., The Importance of Local Chemical Structure for Chemical Metabolism by Human Uridine 5‘-Diphosphate−Glucuronosyltransferase. J. Chem. Inf. Model. 2006, 46, 2692-2697. 13. Ghemtio, L.; Soikkeli, A.; Yliperttula, M.; Hirvonen, J.; Finel, M.; Xhaard, H., SVM Classification and CoMSIA Modeling of UGT1A6 Interacting Molecules. J. Chem. Inf. Model. 2014, 54, 1011-1026.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 29

14. Peng, J.; Lu, J.; Shen, Q.; Zheng, M.; Luo, X.; Zhu, W.; Jiang, H.; Chen, K., In Silico Site of Metabolism Prediction for Human UGT-Catalyzed Reactions. Bioinformatics 2014, 30, 398-405. 15. Rudik, A.; Dmitriev, A.; Lagunin, A.; Filimonov, D.; Poroikov, V., SOMP: Web Server for in Silico Prediction of Sites of Metabolism for Drug-Like Compounds. Bioinformatics 2015, 31, 2046-2048. 16. Kirchmair, J.; Williamson, M. J.; Tyzack, J. D.; Tan, L.; Bond, P. J.; Bender, A.; Glen, R. C., Computational Prediction of Metabolism: Sites, Products, SAR, P450 Enzyme Dynamics, and Mechanisms. J. Chem. Inf. Model. 2012, 52, 617-648. 17. Dang, N. L.; Hughes, T. B.; Krishnamurthy, V.; Swamidass, S. J., A Simple Model Predicts UGT-Mediated Metabolism. Bioinformatics 2016, 32, 3183-3189. 18. Lee, P. W., Handbook of Metabolic Pathways of Xenobiotics (Volumes 1-5). 2014, John Wiley & Sons Ltd.. 19. Kanehisa, M.; Goto, S.; Sato, Y.; Kawashima, M.; Furumichi, M.; Tanabe, M., Data, Information, Knowledge and Principle: Back to Metabolism in KEGG. Nucleic Acids Res. 2014, 42, D199-D205. 20. Kanehisa, M.; Furumichi, M.; Tanabe, M.; Sato, Y.; Morishima, K., KEGG: New Perspectives on Genomes, Pathways, Diseases and Drugs. Nucleic Acids Res. 2017, 45, D353-D361. 21. Rogers, D.; Hahn, M., Extended-Connectivity Fingerprints. J. Chem. Inf. Model. 2010, 50, 742-754. 22. Czarnecki, W. M.; Podlewska, S.; Bojarski, A. J., Robust Optimization of SVM Hyperparameters in the Classification of Bioactive Compounds. J. Cheminf. 2015, 7, 38. 23. Danishuddin; Khan, A. U., Descriptors and Their Selection Methods in QSAR Analysis: Paradigm for Drug Design. Drug Discovery Today 2016, 21, 1291-1302. 24. Pedregosa F, V. G., Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay, E., Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 2011, 12, 2825--2830. 25. Podgorelec, V.; Kokol, P.; Stiglic, B.; Rozman, I., Decision Trees: An Overview and Their Use in Medicine. Journal of Medical Systems 2002, 26, 445-463. 26. Svetnik, V.; Liaw, A.; Tong, C.; Culberson, J. C.; Sheridan, R. P.; Feuston, B. P., Random Forest:  A Classification and Regression Tool for Compound Classification and QSAR Modeling. J. Chem. Inf. Model. 2003, 43, 1947-1958. 27. Zhu, J.; Zou, H.; Rosset, S.; Hastie, T., Multi-class AdaBoost. Statistics and Its Interface 2009, 2, 349-360. 28. Uchaipichat, V.; Winner, L. K.; Mackenzie, P. I.; Elliot, D. J.; Williams, J. A.; Miners, J. O., Quantitative Prediction of In Vivo Inhibitory Interactions Involving Glucuronidated Drugs from In Vitro Data: the Effect of Fluconazole on Zidovudine Glucuronidation. Br. J. Clin. Pharmacol. 2006, 61, 427-439. 29. Riedmaier, S.; Klein, K.; Hofmann, U.; Keskitalo, J. E.; Neuvonen, P. J.; Schwab, M.; Niemi, M.; Zanger,

U.

M.,

UDP-Glucuronosyltransferase

(UGT)

Polymorphisms

Lactonization In Vitro and In Vivo. Clin. Pharmacol. Ther. 2009, 87, 65-73.

ACS Paragon Plus Environment

Affect

Atorvastatin

Page 29 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

For Table of contents Only

ACS Paragon Plus Environment