Predicting Anatomical Therapeutic Chemical (ATC) Classification of Drugs by Integrating Chemical-Chemical Interactions and Similarities

1Nowadays, the Anatomical Therapeutic Chemical (ATC) classification system, recommended by the World Health Organization (WHO), is the most widely recognized classification system for drugs. This classification system divides drugs into different groups according to the organ or system on which they act and/or their therapeutic and chemical characteristics. Accordingly, the ATC classification is very helpful for studying utilization of drugs and categorizing them according to different purposes, therapeutic properties, chemical and pharmacological properties (see Report of the WHO Expert Committee, 2005; World Health Organ Tech Rep, Ser:1–119). In the ATC classification system, drugs are classified into 14 main classes (http://www.whocc.no/atc/structure_and_principles/). In order to understand this kind of complicated classification system, some efforts have been made,. In a pioneer study, Gurulingappa et al. proposed a method to study the ATC-classification system by combining the information extraction and machine learning techniques. However, their method can be used to identify the drug compounds only within the class of “Cardiovascular System”, one of the 14 main ATC classes.

2During the past decade, many compound databases, such as KEGG (Kyoto Encyclopedia of Genes and Genomes),, have been established. From these databases many compounds and their properties can be acquired. Such abundant informations provide an opportunity to analyze ATC classification system in greater detail. Encouraged by the successes of using machine learning and data mining methods to investigate complicated problems in a variety of biological areas,,,,, the present study was initiated in an attempt to develop a powerful method by which one can identify query drugs compound among all their 14 posible main classes.

3According to a recent comprehensive review, to establish a really useful statistical predictor for a biological system, we need to consider the following procedures: (i) construct or select a valid benchmark dataset to train and test the predictor; (ii) formulate the samples concerned with an effective mathematical expression that can truly reflect their intrinsic correlation with the target to be predicted; (iii) introduce or develop a powerful algorithm (or engine) to operate the prediction; (iv) properly perform cross-validation tests to objectively evaluate the anticipated accuracy of the predictor. Below, let us describe how to deal with these steps one by one.

4Recently, the information of protein-protein interactions have been used for predicting various attributes of proteins (see, e.g.,,, ), implying that interactive proteins are more likely to share common biological functions than non-interactive ones. Likewise, it is more likely that two interactive drug compounds may have the similar biological function. Actually, it is generally accepted that compounds with similar physicochemical properties often involve in similar biological activities. Accordingly, it is reasonable to assume that the interactive drugs may likely belong to the same ATC-class, and so do those drugs with similar structures. Based on such rational, let us construct the following benchmark to develop a new method for identifying the ATC-classes of drugs.

5The dataset for drugs was obtained from the public available database KEGG, at ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug (June, 2011). There are totally 9,758 drugs. After excluding those without the information of ATC-codes, the remaining are 4,376 drug samples, from which further screening was performed to remove those without the information of both chemical-chemical interactions and chemical-chemical similarities. After the above winnowing procedures, we finally obtained the benchmark dataset containing 3,883 drugs classified into 14 main ATC-classes, as can be formulated bywhere represents the subset for the 1st main ATC class called “Alimentary tract and metabolism”, the subset for the 2nd main ATC class “Blood and blood forming organs”, the subset for the 3rd main ATC class “Cardiovascular system”, and so forth (cf. Table 1 ); while represents the symbol for “union” in the set theory. For convenience, hereafter let us just use C 1, C 2, C 3, …, C 14 as the tags of the 14 classes. A breakdown of the 3,883 drugs into the 14 main ATC-classes is given in Table 1. For the codes of these drugs in each of the 14 classes, see Supporting Information S1. During the course of constructing the benchmark dataset, the information from http://www.genome.jp/kegg-bin/get_htext?br08303.keg was used that collected the drug compounds and their ATC classification information from http://www.whocc.no/atc_ddd_index/ and provided the ATC code for each drug.

6Because some drugs may belong to more than one main ATC-class, like the case in dealing with proteins with multiple location sites,,, it is instructive to introduce the concept of the “virtual drugs” as illustrated as follows. A drug compound belonging to two different ATC-classes will be counted as 2 virtual samples even though they have an identical chemical structure; if belonging to three different classes, 3 virtual samples; and so forth. Accordingly, the total number of the different virtual drug samples is generally greater than that of the total different structural drug samples. Their relationship can be formulated as follows where is the number of total different virtual drug samples in, the number of total different structural drugs, the number of drugs belonging to one ACT-class, the number of drugs belonging to two ATC-classes, and so forth; while is the number of total main ACT-classes (for the current case, (cf. Table 1 ).

7For the current 3,883 drugs in, 3,295 occur in one class, 370 in two classes, 110 in three classes, 37 in four classes, 27 in five classes, 44 in six classes, and none in seven or more classes ( Figure 1 ). Substituting these data into Eq.1, we havewhich is fully consistent with the figures in Table 1 and the data in Supporting Information S1.

8Based on the fact that the interactive compounds often involve in similar biological activities, it is feasible to predict the ATC-class of a query drug using the information of chemical-chemical interactions, as described below.

9STITCH (Search tool for interactions of chemicals) is a large database containing known and predicted interactions between chemicals and between proteins derived from experiments, literature and other databases. We downloaded the information of chemical-chemical interactions from http://stitch.embl.de:8080/download/chemical_chemical.links.v2.0.tsv.gz. Each of these interactions was evaluated by a confidence score, ranging from 1 to 1000, to reflect the likelihood of its occurrence. For any two drugs d 1 and d 2, their interaction confidence score was denoted by. Particularly, if the interaction between d 1 and d 2 does not exist in STITCH, their interaction confidence score was set as zero, i.e.,.

10Suppose that a training dataset consists of n drugs, and that the 14 main ATC-classes are denoted by, where C 1 represents “Alimentary tract and metabolism”, C 2 “Blood and blood forming organs”, and so forth (see Table 1 ). The ATC-classes of any drug di can be formulated aswhereAccording to the chemical-chemical interaction approach, the likelihood for a query drug belonging to Cj, denoted as, can be calculated bywhere means that is an element of the training dataset. According Eq.6, the likelihood that belongs to Cj can be formulated as the maximum of the interaction confidence scores between and those drugs that belong to Cj in the training dataset. Obviously, the larger the score is, the more likely that belongs to. When, it means that the probability for the drug belonging to the class Cj is zero. Given a query drug compound, suppose the outcome derived from Eq.6 iswhich means that the highest probability for the drug belonging to the ATC-class is (“Antineoplastic and immunomodulating agents”), followed by (“Alimentary tract and metabolism”), and so forth (cf. Table 1 ). If there is a tie between two terms in Eq.7, then the probabilities for the drug belonging to the two corresponding classes are the same. But this kind of tie case rarely happened.

11Note that the outcome of Eq.6 might turn out to be trivial, i.e.,indicating that no chemical-chemical interaction exists for the query drug in the training dataset; i.e.,Under such a circumstance, no meaningful result would be obtained by the “interaction-based” method, and we should instead use the “similarity-based method as described in the next section.

12Likewise, based on the fact that the compounds with similar physicochemical properties often have the same biological activities, we can also use the information of chemical-chemical similarities as another feasible avenue to predict the ATC-class for a query drug. To realize this, let us first introduce how to use graphical representation to measure the similarity between two drug compounds.

13Graphical approaches can provide intuitive pictures and useful insights for studying and analyzing complicated biological systems, as demonstrated by many studies on a series of important biological topics (see, e.g.,,,,,,,,,,, ). Here, a special graphic approach was utilized to estimate the similarity of two compounds. Hattori et al. first proposed a means to measure the similarity of two compounds via their graph representations. Since each chemical structure can be easily represented by a 2D (two-dimensional) graph where vertices stand for atoms and edges for bonds between them, the similarity of two compounds can be estimated by the Jaccard coefficient, based on their maximum common subgraph. The similarity scores between compounds by this method can be obtained from the website at http://www.genome.jp/ligand-bin/search_compound. According to the graphical method by Hattori et al., given two drug compounds d 1 and d 2, their similarity score was denoted by. When the similarity score between d 1 and d 2 does not exist in http://www.genome.jp/ligand-bin/search_compound, their similarity was set as zero; i.e.,.

14Thus, the prediction method based on the chemical-chemical similarities can be formulated in a way almost completely parallel to that of the chemical-chemical interactions as done in the preceding section.

15Now, instead of Eq.6, we havewhere the superscript and subscript “s” stands for the 1st letter of “similarity”, implying that the calculation is now based on “chemical-chemical similarity” instead of “chemical-chemical interaction” as done in Eq.6.

16Given a query drug compound, when the integrated method was used to identify its ATC-class, the prediction involved the following two steps.

17The interaction-based method (cf. Eq.6) was first applied to identify its ATC-class.

18If the probabilities thus obtained for the drug belonging to all the 14 ATC-classes were zero as indicated in Eq.8, meaning no meaningful results were obtained at all, then the prediction would continue using the similarity-based method (cf. Eq.10).

19In statistical prediction, the following three cross-validation methods are often used to examine the quality of a predictor: independent dataset test, subsampling (or k-fold crossover) test, and jackknife test. However, of the three test methods, the jackknife test is deemed the least arbitrary that can always yield a unique result for a given benchmark dataset. The reasons are as follows. (i) For the independent dataset test, although all the samples used to test a predictor are outside the training dataset used to train the prediction engine so as to exclude the “memory” effect or bias, the way of how to select the independent samples for testing the predictor could be quite arbitrary unless the number of independent samples is sufficiently large. This kind of arbitrariness might lead to completely different conclusions. For instance, a predictor achieving a higher success rate than the other for a given independent testing dataset might not able to keep so when tested by another independent testing dataset. (ii) For the subsampling (or k-fold crossover) test, the concrete procedure usually used in literatures was the 5-fold, 7-fold or 10-fold cross-validation. The problem with this kind of subsampling test was that the number of possible selections in dividing a benchmark dataset would be an astronomical figure even for a very simple dataset, as elucidated in and demonstrated by Eqs.28–30 in. Therefore, in any practical subsampling cross-validation tests, only an extremely small fraction of the possible selections were taken into account. Since different selections would always yield different results even for a same benchmark dataset and a same predictor, the subsampling test could not avoid the arbitrariness either. A test method unable to generate a unique outcome should not be deemed as a good one. (iii) In the jackknife test, all the samples in the benchmark dataset will be singled out one-by-one and tested by the predictor trained by the remaining samples. During the process of jackknifing, both the training dataset and testing dataset are actually open, and each sample will be in turn moved between the two. The jackknife test can exclude the “memory” effect. Also, the arbitrariness problem as mentioned above for the independent dataset test and subsampling (or k-fold crossover) test can be avoided because the outcome obtained by the jackknife cross-validation is always unique for a given benchmark dataset. Accordingly, the jackknife test has been widely recognized and increasingly adopted by many investigators to examine the quality of various predictors (see, e.g.,,,,,,,,,,,, ). Accordingly, in this study we are to use the jackknife test to examine the prediction quality as well.

20For any given set of query drugs, we can obtain a series of candidate ATC-classes using the aforementioned prediction methods. Ranked by the likelihood according to their descending order, the prediction accuracy can be defined aswhere CPj denotes the number of drugs whose j-th order predicted ATC-class is one of the true ATC-class, and N denotes the total number of query drugs whose ATC-classes are to be identified. According to such a definition, the result of higher ACj with smaller j or lower ACj with larger j indicates that the predicted hits are more concentrated meaning a better prediction. Obviously, the result with high 1st-order prediction accuracy AC 1 always represents a good quality of prediction.

21The average number of ATC-classes for the N query drugs is defined aswhere Ti is the number of ATC-classes for the i-th query drug. Thus, another parameter for measuring the proportion of the true classes successfully identified by the first m-order prediction hits can be calculated as where Pi,m denotes the number of the first m predicted candidate ATC-classes that are the true ATC-classes for the i-th drug in the dataset. Usually, m could take the smallest integer that is equal to or greater than AN; i.e.,where the operator Int means taking the integer part of the quantity right after it. Again, the result of larger Lm with smaller m implies a better prediction with less uncertainty.

22For clarity, the original benchmark dataset of 3,883 drugs (cf. Supporting Information S1) can be separated into two subsets; i.e.,where contains 2,144 drugs that had the chemical-chemical interaction information, while contains drugs that had no chemical-chemical interaction information. Listed in Table 2 are the results obtained by the aforementioned three different prediction methods in identifying the 14 main ATC classes for the drugs investigated. By examining the table, we can observe the following.

23For the 2,144 drugs in we could use Eq.6 to conduct the prediction. The results thus obtained are listed in column 2 of Table 2, from which we can see that the 1st-order prediction by the jackknife test on the 2,114 drugs was 67.72%. The success rates generally followed a descending trend with increasing of the order number, indicating that the predicted ATC-classes were well sorted for each of the samples investigated. The average number of the ATC-classes in was (see Eq.12). Thus, it follows according to Eq.14 that, meaning that the first 2-order predictions should be taken into consideration. Substituting these data into Eq.13, we obtained the overall success rate by the predictions of the first two orders for the 2,144 drugs in was, indicating that the interaction-based method is quite promising in identifying the ATC-classed of drugs. However, this method could only be used to deal with those drugs that had the chemical-chemical interaction information.

24For the remaining 1,739 drugs in the dataset (cf. Eq.15) that did not have the chemical-chemical information, the similarity-based method (cf. Eq.10) was used as a backup, and the results thus obtained are shown in column 3 of Table 2. It can be seen from there that the 1st-order prediction by the jackknife test on the 1,739 drugs was 78.49%. The average number of ATC-classes for the drugs in was (see Eq.12), and hence we have (Eq.14), meaning that the first 2-order predictions should be taken into account. Substituting these data into Eq.13, we obtained the overall success rate by the first two orders predictions for the 1,739 drugs without the chemical-chemical interaction information was 75.31%, indicating that the similarity-based method was quite promising as well.

25At a first glance at Table 2, it looks like that the success rates by the similarity-based method (Eq.10) are higher than those by the interaction-based method (Eq.6). However, since the success rates by the two methods as reported in Table 2 were derived from two different datasets, and (cf. Eq.15) respectively, they might not able to reflect the true superiority between the two methods. To make a comparison between them in a more fair manner, let us construct a new dataset, denoted as. It consists of 2,138 drugs with each containing both chemical-chemical interaction and chemical-chemical similarity informations. The details of such a dataset is given in Supporting Information S2.

26Listed in Table 3 are the results obtained by the methods in identifying the 14 main ATC classes for the 2,138 drugs in the dataset. As we can see from the table, the 1st-order prediction accuracy by the interaction-based method was 67.40%, while that by the similarity-based method was 40.36%.

27The average number of ATC-classes for the drugs in was 1.24 (see Eq.12), and hence we have (Eq.14), meaning that the first 2-order predictions should be taken into account. Substituting these data into Eq.13, we obtained the overall success rate by the 1st two orders predictions for the 2,138 drugs in by the interaction-based method (Eq.6) was 71.26%, while that by the similarity-based method (Eq.10) was only 43.69%, indicating that the interaction-based method is superior to the similarity-based method in identifying the ATC-classes of drugs. That is why in the integrated method the first step was to use the interaction method (Eq.6) to identify the ATC-classes for any query drugs. When, and only when no meaningful result was obtained by the interaction-based method, was the similarity-based method (Eq.10) used as a backup to continue the prediction (see the Section of “Prediction by Integrating the Interaction-Based and Similarity-Based Methods”).

28Shown in the 4th column of Table 2 are the results obtained by the integrated method in identifying the 14 main ATC classes for the 3,883 drugs in the benchmark dataset. As we can see there, the 1st-order prediction accuracy was 72.55%. The average numbers of ATC-classes for the drugs in was (see Eq.12). Thus, it follows according to Eq.14 that, meaning that the first 2-order predictions should be taken into consideration. Substituting these data into Eq.13, we obtained the overall success rate by the first two orders predictions for the drugs in was 73.25%.

29These results indicate that the integrated method performed quite well in identifying drugs among their 14 main ATC-classes, and that more attention should be paid to the results hit by the first two order predictions because they covered more than 70% of the true ATC-classes.

30Finally, it is instructive to point out that although the above demonstrations were given for identifying query drug compounds among their main (or 1st level) classification, the method developed here can be straightforwardly extended to cover the 2nd, 3rd, 4th, 5th or any lower-level classification as long as the corresponding statistically significant datasets for training the predictor are available.