SePreSA: a server for the prediction of populations susceptible to serious adverse drug reactions implementing the methodology of a chemical-protein interactome.

Serious adverse drug reactions (SADRs) are caused by unexpected drug-human protein interactions, and some polymorphisms within binding pockets make the population carrying these polymorphisms susceptible to SADR. Predicting which populations are likely to be susceptible to SADR will not only strengthen drug safety, but will also assist enterprises to adjust R&D and marketing strategies. Making such predictions has recently been facilitated by the introduction of a web server named SePreSA. The server has a comprehensive collection of the structural models of nearly all the well known SADR targets. Once a drug molecule is submitted, the scale of its potential interaction with multi-SADR targets is calculated using the DOCK program. The server utilizes a 2-directional Z-transformation scoring algorithm, which computes the relative drug-protein interaction strength based on the docking-score matrix of a chemical-protein interactome, thus achieve greater accuracy in prioritizing SADR targets than simply using dock scoring functions. The server also suggests the binding pattern of the lowest docking score through 3D visualization, by highlighting and visualizing amino acid residues involved in the binding on the customer's browser. Polymorphism information for different populations for each of the interactive residues will be displayed, helping users to deduce the population-specific susceptibility of their drug molecule. The server is freely available at http://SePreSA.Bio-X.cn/.

INTRODUCTION

Drug effect varies among populations. The Japanese population, for instance, exhibit a more rapid response to the lung cancer therapy gefitinib (1), since polymorphisms within the binding pocket of the drug target increase sensitivity to inhibition by the drug (2–4), and these polymorphisms occur more frequently in Asian populations. Serious adverse drug reaction (SADR), an unwanted drug effect, has been an urgent world-wide problem, particularly as tragedies triggered by Vioxx® (5) and Avandia® (6) these years. SADRs, especially type B adverse drug reactions (ADRs) (7), are mainly caused by unexpected interactions of the drug with the SADR targets (8,9), and some polymorphisms within the binding pocket make the population carrying these polymorphisms more susceptible to harmful effects. If drug companies had been able to identify the sensitive population earlier, they should have altered their R&D and marketing strategy beforehand to lower the rate of SADR and to avoid lawsuits. Furthermore, such prediction would suggest candidate polymorphisms for SADR association studies (10), provide hints for interpreting genome-wide association results for SADR and provide primers for functional studies of the SADR mechanism.

Three steps could lead to the prediction of susceptible populations. First, SADR targets which tend to be bound by the compound in question should be prioritized. Second, the conformation of chemical–protein bindings and the interactive residues should be identified. Third, polymorphisms altering drug binding and their minor allele frequency among different populations could then be characterized. Docking programs could be applied to perform the above steps, as the program can not only measure the binding strength of a ligand to a set of SADR targets (11), but can also deliver the binding conformation. Consequently, residues involved in the interaction together with information on their polymorphisms in different populations can be identified. For the first time, these three steps have been integrated into a server named SePreSA. The server computes a relative drug–protein interaction score from a scoring matrix of a chemical–protein interactome (CPI) to prioritize SADR targets which might be affinitive with the user's compound. Considerable specificity and sensitivity in predicting chemical–protein bindings have been achieved by the use of this scoring algorithm. The 3D visualization of the binding pattern is provided, with interactive residues highlighted and the population information presented to the client.

METHODS

Target set for the web server

To make a comprehensive collection of the structural models of the well known SADR targets, we chose all the available structures from PDB which are known to mediate ADRs. The target set includes major phases I and II drug-metabolite enzymes, several types of human MHC I proteins mediating drug hypersensitivity and the pharmacodynamic proteins chosen from DITOP (12) and DART (13) database. All structure models chosen should accord with the following criteria: (i) the species of the proteins is limited to Homo Sapiens; (ii) the protein must contain at least one ligand embedded in it to define the functional site; (iii) no missing residues should be found around this site; and (iv) ligands at the surface of the protein are not acceptable. Now, it contains 91 proteins with 115 ligand-binding site defined. For each ADR target, residues within a 10 Å distance from the ligand were defined as the bioactive pocket of the protein, and balls with a radius ranging from 1.1 to 1.4 Å were generated to fill in the pocket. A grid box was made at 3–5 Å distant from the ‘cloud’ of the balls. Certain water molecules or metal ions play important roles to the protein function, so we used the scientific judgment to decide whether to keep them. The key residue ionization state has been assigned considering the most probable one at the physiological pH, i.e. carboxyl is usually ionized; lysine and arginine residues are protonated; aspartate and glutamate residues are deprotonated; histidine residues are half protonated. We controlled these ionization state using Chimera (14) when preparing the targets. The polymorphism information for each ADR targets is derived from Uniprot (15). We will continue to update ADR targets in SePreSA, and users can subscribe to our updates through RSS feeds.

Dataset for the prediction evaluation

We singled out all co-crystallized ligands embedded in all structures of SePreSA. The antagonists taking up the functional site of the protein were chosen as the probe molecule. Pockets without co-crystallized ligands were excluded, leaving 79 proteins for the construction of the CPI. The molecular probes were submitted through the SePreSA interface to perform an in silico hybridization, generating an interactome of 86 ligands towards 79 protein pockets in the form of a docking-score matrix of 79 × 86 elements. Docking scores ≥ 0 were treated as missing values.

The prediction evaluation

An algorithm named 2-directional Z-transformation (2DIZ) was applied to process the original docking-score matrix. Here, X represents the docking scores of ligand j to protein i. The Z-score was calculated as: where Here, N equaled 86 minus the number of missing values of ligand j. Thus, a Z-score matrix of 79 × 86 elements was generated. The vector for each protein was then normalized to a mean of zero and a standard deviation of one, generating a 79 × 86 Z′-score matrix. These three matrixes allowed us to investigate the distributions of docking scores, Z-scores and Z′-scores on true ligand–protein bindings and the unidentified bindings. Here, we defined the original bindings in PDB structures as the gold standard, and the ability of these three scoring matrixes to predict ligand–protein bindings were presented in ROC curves.

INPUT AND OUTPUT

Users need to upload a drug molecule in mol2 format. A manual at http://sepresa.bio-x.cn/?page=generatemol2file will instruct users in preparing their mol2 files. The server cannot accept molecules in SMILES and mol format, since the ionization state of the drugs have to be specified by users but these two formats cannot include the ionization information of the molecules. The format suitability is checked and its interactome towards all SADR targets was calculated using DOCK (16). The task usually takes up to about 6 h for a molecule, and an email will be sent on completion. The outputs comprise the following three elements.

SADR targets that tend to interact with your molecule will be prioritized.

For each SADR target, binding patterns of the lowest docking score and amino acid residues that interact with your molecule will be highlighted in Jmol (17) applet.

Polymorphism information such as minor allele frequencies among different populations for each of the interactive residues will be displayed.

RESULTS

Prediction of the true and unidentified bindings

We compared the prediction power using different CPI scoring matrixes (Figure 1). The docking-score matrix performed poorest among the three, the area under the curve (AUC) being only 0.62 (Supplementary Table S1), and the lower bound of 95% confidence interval (95% CI) 0.56, which was close to that for a random selection. With the 2DIZ algorithm, however, the AUC reached 0.82 (95%CI: 0.78–0.87). Though performing slightly worse than the Z′-score matrix, the Z-transformation achieved a better performance than simply using the docking-score matrix.

Figure 1.

The ROC curves of different CPI-scoring matrixes in predicting true and unidentified.

The distribution of Z′-scores between true and unidentified bindings are compared in Supplementary Figure S1a. Z′-score for ∼80% of the true bindings, compared with only 30% of the unidentified bindings, were <–0.5 (Supplementary Table S2). Hence, we set a Z′-score threshold of –0.5 to highlight the putative bindings of users’ drugs towards SADR targets, and the sensitivity, specificity and the overall accuracy were 0.80, 0.71 and 0.71, respectively. Docking-score distributions of the true bindings did not seem to be significantly different from those of the unidentified bindings (Supplementary Figure S1b).

The sensitivity of the Z′-score-based prediction was diminished by the high Z′-score of several true bindings. In most circumstances, these false negatives were due to the large size of the probe. Glutathione, a relatively large molecule, could not be docked into 20% of the pockets, resulting in a number of missing values. Consequently, its Z-score vector did not fully reflect its binding profile for all pockets, causing a potential bias in Z′-score of for its co-crystallized enzyme (PDB ID: 11GS). Another reason for the false negative lay in the poor selectivity of the probe. Ethyl dihydrogen phosphate, a small molecule which could possibly crawl into the pocket of all the proteins, did not appear to be selectively affinitive to its co-crystallized protein (1XLV), and thus its Z′-score for 1XLV was not noticeably low.

The specificity of the Z′-score-based prediction might also be much higher than current results, because some of the unidentified bindings, whose Z′-scores were significantly <–0.5, might occur per se, but were regarded as false positives. For example, the Z′-score between catechol O-methyltransferase (3BWY) and S-adenosyl-l-homocysteine (SAH) was –3.1. SAH was structurally similar to S-adenosylmethionine (SAM), which was the original ligand embedded in 3BWY. On the other hand, SAH was originally embedded in nicotinamide N-methyltransferase (2IIP), which belongs to the family of 3BWY. Hence, SAH is very likely to bind to 3BWY i. If this is true, the strong signal indicated by the low Z′-score of the SAH–3BWY interaction could indicate that the 2DIZ scoring algorithm is capable of prioritizing unexpected bindings.

Applying 2DIZ algorithm to the web server

The SePreSA server uses 2DIZ algorithm to prioritize SADR targets of the user's molecule. The prediction mechanism is based on a user-oriented interactome, which calculates the Z′-score of the current molecule from the interactome formed by all molecules submitted by this user, no matter when and where these previous molecules are submitted. Hence, the more molecule a user submits, the more comprehensive CPI profile for each ADR targets he will retrieve.

Case study 1

Serious cutaneous reaction (SCR) triggered by sulfamethoxazole (SMX) might be mediated by the MHC I family members (18). After submitting SMX to SePreSA, we found unexpectedly that the HLA class I histocompatibility antigen, Cw-4 alpha chain (MHC I Cw*4) ranked in the fourth among the total 70 SADR target pockets. By visualizing the binding conformation of the SMX molecule to MHC I Cw*4 in the binding-information page (http://sepresa.bio-x.cn/?reactionid=5069), we found that it tended to ‘root’ at the Y bed of the antigen presentation groove. The identification of HLA-C gene (Cw*4 allele) as the mediator of SCR had been validated in several former studies, from which it was confirmed that the SCR could only be triggered by SMX in presence of MHC I (Cw*4) (19,20).

To our knowledge, no other research has ever disclosed such direct binding model before the results given by SePreSA. Although the identification of this risk allele needs further validation, several ‘wet’ observations support this model. For example, the presentation of SMX parent drug displayed a direct, non-covalent binding fashion to the ‘empty’ presentation groove of MHC I (21). Von Greyerz et al. (19) discovered that most T cell clones exhibited the ‘MHC-allele restricted drug-specific recognition’ stimulated by SMX parent drug. Nassif (20) also uncovered that blister fluid T lymphocytes, which were derived from a patient suffering SMX-induced SCR, could be cytotoxic only when SMX is present in the cells with Cw*4 allele.

Case study 2

The incidence of neuropsychiatric disorder triggered by oseltamivir (Tamiflu®) varies among populations. The Japanese population demonstrates a higher SADR rate than predominantly European-derived populations for this drug (http://www.fda.gov/ohrms/dockets/AC/05/briefing/2005-4180b_06_01_Tamiflu%20AE_reviewed.pdf). The SADR target of oseltamivir (HsNEU2) (8) was also included in our pocket set. We, therefore, submit the active form of the drug, oseltamivir carboxylate, to ascertain whether the HsNEU2 and the susceptible population could be prioritized. The molecule achieved a Z-score of –1.72 with HsNEU2, the second lowest Z-score for 70 pockets, and a Z′-score of –0.92 temporarily, which was much lower than the true binding threshold of –0.5. By then clicking on the ‘Result’ button, the binding conformation and the interactive residues within 6.4 Å of the drug were presented (Figure 2), together with a polymorphism (rs2233385) highlighted among these residues. By clicking on the ‘Show Report’ button, we found that this polymorphism only occurred in Asian population but not in European and African–American populations (Supplementary Table S3). These results suggested that Asian population might be more sensitive to oseltamivir-induced SADR, for which rs2233385 might be responsible.

Figure 2.

Binding conformation of oseltamivir to HsNEU2 and the interactive residues within 6.4 Å of the drug. Among all these residues, the R41Q polymorphism (rs2233385) was highlighted.

Such prediction potential will assist further mechanism studies and genetic association studies on HsNEU2-mediated SADR. If all these facts are approved, the manufacturer of oseltamivir, for example, could improve drug safety through changing the marketing strategy or utilizing pharmacogenomic tests. With the benefit of SePreSA predictions during the early development phase, manufacturers would have the opportunity to redesign or modify the drug in order to weaken such unexpected binding, or they might even give up the Asian market to avoid unexpected lawsuits.

DISCUSSION

Though a lot of downstream events occur when a drug is added to the cell culture, it is undisputed that direct chemical–protein interaction is the primary and the vital factor in drug effects. So, identifying the true bindings of unexpected drug–protein interactions is fundamental in pharmacodynamic research and in the prediction of effects including SADR. Several techniques such as BIACORE® biosensors (22) and drug affinity pull-down (23) can be used to assess such interactions. However, these techniques do not match the dramatic progress achieved by transcriptomics, metabolomics and proteomics. The concept of docking a small molecule into a multi-protein set to prioritize unexpected bindings was first put forward by Chen et al. (11). Several follow-up studies have pursued this logic in prioritizing true targets (24,25), namely that the lower the docking score achieved, the more this binding tends to happen. But this approach does not offer a systematic evaluation of relevant specificity and sensitivity. The docking score might not be sufficient to evaluate the binding strength, e.g. if the docking score of drug A to protein P1 is much lower than A to P2, there is no certainty that P1 is more affinitive to A than P2. However, by considering the mean and the standard deviation of the score vectors upon the two proteins towards multiligands, a more informed judgment can be made.

SePreSA, is the first system to utilize the drug–protein interaction landscape at an interactome level to help users make sound decisions. Although the docking-score matrix of the test CPI now contained 79 × 86 elements, from which the magnitude at either rows or columns did not seem to be very impressive, it already had a total of about 79 × 86 = 6794 ligand–protein pairs to be identified. So, we believe that the classification performances generated at such amount of data can reflect the true performance of Z′-scores to some extent. To our knowledge, this is not only the first, but also the largest evaluation in the target ‘fishing’ methodology using molecular docking in company with clear reported sensitivity, specificity and accuracy data. The experience gathered from using this system also suggests that the use of relative scores from the ‘-omics’ viewpoint can achieve much greater accuracy than simply comparing the docking scores of the two independent interactions. Our algorithm might also inspire the existing virtual screening methodologies. If the interactome profiles of the library molecules towards multiproteins are considered, more accurate results can be achieved.

In this research, our underlying logic was that drug effects would necessarily change when the binding of a drug to its target is altered, and polymorphisms involved in this direct interaction would necessarily change the binding. We have seen that, although drug response is a complex trait (26) mediated by multiple genes, some single polymorphism can also have pronounced effects on drug response. To our knowledge, they all alter the binding conformations of direct drug–protein interactions. Examples include the T790M in the gefitinib binding pocket of EGFR (4); the T164I within the epinephrine binding pocket of β2-adrenergic receptor (27); and the polymorphism within the binding pocket of STI-571 to c-Abl (28). The empirical threshold of 6.4 Å was set to highlight the putative interactive sites according to the distance distribution of drugs to the polymorphism sites that alter drug binding. So SePreSA cannot predict every polymorphism that alters drug binding, but it can predict the interactive residues within the 6.4 Å ‘cloud’, whose polymorphism information are available. Both the premise and the empirical threshold will be more thoroughly evaluated in follow-up research.

CONCLUSION

The core of the SePreSA server is the 2DIZ scoring algorithm. It can accurately predict bindings of a chemical towards multiproteins, and hence could be applied in prioritizing SADR targets.

By using SePreSA, drug enterprises can identify the putative populations that appear sensitive to their drugs, hence early decision during the R&D stage can be made and safety can be promoted in the marketed products.

SADR genetic researchers could find candidate polymorphisms from SePreSA for their SADR association studies. The server could also help to interpret genome-wide association results for SADR and enhance functional studies of the SADR mechanism.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

FUNDING

Major State Basic Research Development Program (2006CB910601) and National Key Technology R&D Program (2006BAI05A05). Funding for open access charge: China Postdoctoral Science Foundation (20070420660) and the Shanghai Postdoctoral Science Foundation (No. 61444).

Conflict of interest statement. None declared.