|Home | About | Journals | Submit | Contact Us | Français|
Advancements in next-generation sequencing (NGS) technologies are generating a vast amount of data. This exacerbates the current challenge of translating NGS data into actionable clinical interpretations. We have comprehensively combined germline and somatic nonsynonymous single-nucleotide variations (nsSNVs) that affect drug binding sites in order to investigate their prevalence. The integrated data thus generated in conjunction with exome or whole-genome sequencing can be used to identify patients who may not respond to a specific drug because of alterations in drug binding efficacy due to nsSNVs in the target protein's gene. To identify the nsSNVs that may affect drug binding, protein–drug complex structures were retrieved from Protein Data Bank (PDB) followed by identification of amino acids in the protein–drug binding sites using an occluded surface method. Then, the germline and somatic mutations were mapped to these amino acids to identify which of these alter protein–drug binding sites. Using this method we identified 12993 amino acid–drug binding sites across 253 unique proteins bound to 235 unique drugs. The integration of amino acid–drug binding sites data with both germline and somatic nsSNVs data sets revealed 3133 nsSNVs affecting amino acid–drug binding sites. In addition, a comprehensive drug target discovery was conducted based on protein structure similarity and conservation of amino acid–drug binding sites. Using this method, 81 paralogs were identified that could serve as alternative drug targets. In addition, non-human mammalian proteins bound to drugs were used to identify 142 homologs in humans that can potentially bind to drugs. In the current protein–drug pairs that contain somatic mutations within their binding site, we identified 85 proteins with significant differential gene expression changes associated with specific cancer types. Information on protein–drug binding predicted drug target proteins and prevalence of both somatic and germline nsSNVs that disrupt these binding sites can provide valuable knowledge for personalized medicine treatment. A web portal is available where nsSNVs from individual patient can be checked by scanning against DrugVar to determine whether any of the SNVs affect the binding of any drug in the database.
With the development of massively parallel sequencing, also known as next-generation sequencing (NGS), a vast amount of NGS data are being generated with greater throughput and decreased cost compared with its predecessor technology, Sanger sequencing.1, 2, 3, 4 The identification of single-nucleotide variations (SNVs) is one of the most common tasks in NGS data analysis.5 Although most SNVs are found in the intergenic region, many can also be found at more crucial locations such as within protein coding regions. For example, missense SNVs (msSNVs) may play a more direct role in causing or exacerbating disease by changing the protein structure or by other mechanisms.6, 7 Pharmacogenetic and pharmacogenomic studies have shown that SNVs can affect how a patient responds to administered drug treatments.8, 9, 10, 11 The most direct example is where msSNVs exists within the coding region of a gene coding for target protein and these SNVs alter the amino acid of binding site of the drug, resulting in changes in drug binding affinity and consequent therapeutic effect of the drug.12, 13, 14 In addition to msSNVs affecting the binding site, there are several other factors that can cause differences in drug responses including gene expression, drug metabolism and environmental factors.15 In this study we focus on the effects of msSNVs that affect protein–drug binding sites.
The Protein Data Bank (PDB),16 a three-dimensional structure database, contains structure data of proteins complexed with small molecules such as substrates, cofactors, inhibitors and drugs and is widely used in drug discovery research.17, 18, 19 Secondary databases, such as DGIdb, CREDO20 and FireDB,21 use data from PDB and provide value-added information through further analysis. The above databases do not provide comprehensive somatic mutation or polymorphism mapping; neither do they provide protein–drug interaction-centric information. Research has shown that individual's genetic makeup can contribute to differential drug response.8, 22, 23, 24 PharmGKB mines information from this type of research publications,25 and at the time of writing this paper it contained over 5000 variant annotations in more than 900 proteins related to over 600 drugs.
This paper describes identification and integration of amino acid–drug binding sites from PDB and nonsynonymous single-nucleotide variations (nsSNVs) compiled from various sources to create a comprehensive data set called DrugVar. This data set can be used to scan exome or whole-genome sequencing data from patients to see whether a patient has a missense mutation that affects any protein–drug binding site, and also how prevalent this SNV is in tumors or whether they exist as polymorphism in the population.
Amino acid–drug complex structure data was obtained from PDB database.26 The Anatomical Therapeutic Chemical Classification (ATC), a hierarchical representation of drugs, was used for identifying cancer and noncancer-related drugs.26 Data were manually curated in order to separate drugs from other small molecules. PDB amino acid sequence positions were mapped to UniProtKB accessions27 followed by pairwise alignment. For each atom of the amino acid residues in the drug binding pocket, occluded surface (buried) area by the drug in its binding pocket was computed using the program OS28 and a list of drug-interacting residues was generated based on ranking the occlusion percentage. All PDB IDs associated with the same protein and all protein–drug pairs were considered.
A comprehensive non-redundant data set of both germline and somatic nonsynonymous nsSNVs was generated using data from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/), International Cancer Genome Consortium (ICGC),29 IntOGen,30 CSR,31 cancer cell lines from NCI-60 panel,32 dbSNP33 and UniProtKB/Swiss-Prot.34 The integration was performed using methods described earlier.35 The final data set contains 1705286 somatic nsSNVs and 1132832 germline nsSNVs. Drug IDs were retrieved from DrugBank and literature to allow linking to detailed pharmacology and pharmaceutical knowledge and drug target information.36, 37, 38 Figure 1 shows the workflow used to generate the protein–drug interaction data set.
A statistical significance analysis was conducted to investigate the distribution of binding sites between amino acids and drugs. Significance of observed versus expected of amino acids in drug binding sites was calculated using the Binomial statistic described by Mi and Thomas39 and applied in our previous studies.40 Protein Functional Class Enrichment analysis was performed using the PANTHER Classification System.41
For paralog identification, each protein in the data set was used to BLAST42 against the human proteome to identify paralogs43, 44 that might also bind the drug. The binding sites in the paralogs were checked using pairwise alignment to ensure they are conserved. Homology modeling of these paralogs was carried out using MODELLER.45 These homology models of paralogs were energy minimized using the AMBER force field in Molecular Operating Environmen (MOE) (Version 2014.06) package developed by the Chemical Computing Group, Montreal, Canada (www.chemcomp.com). For additional homolog identification, all proteins of non-human mammalian origin–drug complex were retrieved from the PDB database. Similar workflow as mentioned above was applied on these data for identification of human homologs that may serve as drug targets.
Visualization of the amino acids in the binding site for the amino acid–drug complexes was carried out using UCSF CHIMERA.46 To understand the effect of mutations on a drug binding to a protein target, the X-ray structure of carbonic anhydrase 2 (CA2) bound to the drug Lacosamide (PDB ID: 3IEO)47 was chosen. The wild-type amino acid–drug and the mutated amino acid–drug complexes were energy minimized using the AMBER force field in MOE and binding energies computed.
Even if there are no amino acid variations in protein–drug binding sites, the variations existing in a binding site's neighboring region may structurally block the protein–drug binding. Therefore, for each protein–drug binding site, mutations within its upstream and downstream region were identified. The total number of variations within the expanded protein–drug binding region was calculated. In addition, as the gene expression level is another factor that may influence drug efficacy, cancer-related gene expression data were retrieved from BioXpress.48 It analyzes and normalizes raw count data of paired sample (tumor and normal) from TCGA using DEseq R package (http://bioconductor.org/packages/release/bioc/html/DESeq.html), a method that is regarded as one of the most robust RNA-seq normalization methods.49 If a protein–drug pair contained one or more amino acid variations in its binding site and if the variations were somatic mutations associated with a specific cancer type, the corresponding tumor versus healthy tissue differential gene expression was evaluated to rank promising validation targets.
Out of 112722 structures available in the PDB database, 20385 PDB structures are associated with 1712 nonstandard polymeric components. After further filtering of crystallization stabilizers and proteins of non-human origin, a total number of 827 PDB structures were identified with drug bound to them. These PDB structures were mapped to 253 UniProtKB proteins binding to 235 unique drugs with 12993 binding sites (Supplementary Table S1). Out of 253 proteins in the amino acid–drug binding data set, 210 proteins had 1618 drug binding site altering variations (Supplementary Tables S2 and S3). Figure 2 shows the distribution of binding sites and binding sites affected by SNVs across 253 drug target proteins. For each drug target protein, the blue bar indicates the ratio between number of drug binding sites and target protein length, whereas the red bar shows the ratio between number of binding site affecting SNVs and binding sites. In 253 identified drug target proteins, on average 4.8% of the amino acids are identified as drug binding sites and 14.8% of these drug binding sites have variation mapped to them. The number of binding sites for each protein is largely influenced by the number of ‘binding pockets' it contains and the structure of ‘binding pocket'. In all, 1162 binding sites are affected because of germline variations and 456 binding sites are affected by somatic mutations. The number of identified drugs and their target proteins is relatively small compared with the total number of FDA (Food and Drug Administration)-approved drugs (1602) because of the limitations of X-ray crystallography technology. With the advent of large-scale structural genomics projects, we expect to see increased coverage in the coming years. In addition, binding sites based on electrospray ionization mass spectrometry50 and mutation studies51 can be used to model additional sites.
It is interesting to note that 83% of the identified drug target proteins contain mutations on their binding sites. Out of the large amount of identified mutations, 0.12% of germline mutation and 0.035% of somatic mutations are mapped to drug binding sites in our study. Amino acid–drug binding sites and nsSNV sites are available from the DrugVar portal and a snapshot of the web interface is shown in Figure 3. All drugs were classified into 14 groups based on ATC (Table 1).
Out of 12993 drug binding sites, 1408 protein binding sites were found to be associated with 25 antineoplastic drugs (inhibiting or preventing the growth and spread of tumors or malignant cells; ATC Code: L01). Some of them act by disrupting the metabolism of DNA synthesis components. For example, folic acid analogs, such as methotrexate and raltitrexed, inhibit key enzymes that are necessary for synthesis of thymidylate, an essential component of DNA.52, 53 Similarly, purine analogs, such as cladribine, and pyrimidine analogs, such as gemcitabine, inhibit cancer cell growth by disrupting DNA synthesis, whereas nilotinib binds to the catalytic site of BCR-Abl, a protein kinase that only appears in cancer cells of chronic myelogenous leukemia patients, and inhibits the tumor growth.
We further examined 25 unique antineoplastic agent drugs and 35 unique target proteins. Figure 4 shows a Circos plot54 representing the protein–drug binding relationship. Imatinib, a tyrosine kinase inhibitor used in the treatment of multiple cancers has been shown through solved three-dimensional structures to bind to as many as eight proteins (gene symbols: ABL1, ABL2, DDR1, KIT, LCK, MAPK14, NQO2 and SYK).55, 56, 57, 58, 59, 60 As one of the most well-known targeted therapy drugs for cancer treatment, imatinib inhibits the BCR-Abl kinase to cut off the energy source for cancer proliferation, thereby inhibiting the tumor growth.55, 61 It is also known that imatinib can also inhibit mast/stem cell growth factor receptor Kit (KIT)56 to stop tumor growth. However, the effect of imatinib on other proteins in Figure 4 still remains largely unexplored even though structures are available. Sunitinib can bind to six proteins (CDK2, KDR, KIT, PHKG2, ITK and MAPK14). It inhibits cellular signaling by targeting tyrosine kinase including platelet-derived growth factor receptor-β, vascular endothelial growth factor receptor 1, KIT and vascular endothelial growth factor receptor 2.62,63,64,65 However, further studies are needed for the physiological outcome of its binding to rest of the proteins.
For the antineoplastic drug target proteins, 608 out of 3133 variations were found on binding sites between 25 antineoplastic drugs and their 36 target proteins. Out of these, 178 are somatic mutations and 429 are germline variations. Gefitinib is primarily used in treatment of non-small-cell lung cancer.66 In the DrugVar data set, six lung cancer-associated somatic mutations exist within binding sites between epidermal growth factor receptor (EGFR) and gefitinib on amino acid positions 744, 766, 790, 792, 844 and 855. In addition, a total number of seven germline mutations exist within binding sites between EGFR epidermal growth factor receptor and gefitinib on amino acid positions 719, 726, 743, 788, 796, 800 and 854. It is possible that these mutations on EGFR may weaken binding efficacy of gefitinib to its target protein in lung cancer patients.
It is reasonable to assume that the more proteins a drug binds to, the more diverse its effects might be. Imatinib has eight target proteins (ABL1, LCK, KIT, NQO2, ABL2, SYK, DDR1 and MAPK14) in the DrugVar database. Out of them, 7 proteins have a total of 58 amino acid mutations in their binding sites with imatinib. As one of the primary target proteins, KIT is inhibited by imatinib67 and contains 11 possible variations in its imatinib binding sites. These mutations may affect the drug binding affinity and make imatinib less effective under regular dose intake. Tyrosine-protein kinase ABL1 (ABL1) is another imatinib targeted protein and contains 29 variations on binding sites. Out of 29 mutations, 19 are germline mutations and 8 are somatic mutations associated with diffuse large B-cell lymphoma. We report here many such amino acid variations in drug binding sites of proteins. The biological outcomes of these variations have not yet been elucidated. Such examples include epithelial discoidin domain-containing receptor 1 (DDR1) that is suggested to be a pharmacologically available target for cancer treatment,60, 68 tyrosine-protein kinase Lck (LCK), tyrosine-protein kinase SYK (SYK), mitogen-activated protein kinase 14 (MAPK14) and abelson tyrosine-protein kinase 2 (ABL2).
To better understand the binding of the drugs, the binding pocket of one of the drugs was further analyzed. Figure 5a shows the superposition of X-ray crystal structures of imatinib binding to its eight target proteins. The drug binds in the active site cleft of the kinase catalytic domain that is represented by ribbon structures. Figure 5b shows the side chains identified as binding site of amino acids with imatinib. They are marked by a green color in the pocket of target proteins. The structural differences could help to modify imatinib to improve the selectivity for a given protein target. Figure 5c shows the amino acids that are reported as mutations for each target protein in the amino acid–drug binding sites. The drug is located in the active site cleft of the kinases, each of which is shown in different color for clarity.
Our study also emphasizes the importance of proteins that bind to specific drugs while still containing a large number of mutations on their binding sites. EGFR is a well-known protein target for antineoplastic drugs.69 Some drugs, such as lapatinib and erlotinib,70, 71 have been designed specifically as antagonists of EGFR to prevent its activation that is required for tumor proliferation.71, 72 We identified 25 mutations within EGFR binding sites with lapatinib, suggesting that lapatinib might vary in effectiveness based on the genomic profile of the patient.
There are 11124 amino acid–drug binding sites associated with 200 non-antineoplastic drugs. These 200 non-antineoplastic drugs can be further categorized into 13 subclasses. Table 1 shows the number of proteins and drugs associated with each group categorized by ATC Classification System. Out of a total of 3133 amino acid variations, 2704 were found on 200 non-antineoplastic drug target proteins.
There are seven proteins in this data set that have at least seven drugs that can bind to them. The proteins are androgen receptor (AR), deoxycytidine kinase (DCK), aldo-keto reductase family 1 member C3 (AKR1C3), aldo-keto reductase family 1 member C2 (AKR1C2), transthyretin (TTR), CA2 and serum albumin (ALB). Most of them play important roles in drug transportation and in metabolic pathways. For example, ALB, a major plasma protein that functions as drug transporter, has 18 drug–protein complex structures available in PDB and contains 97 amino acid variation sites within the binding sites for these drugs. Similarly, TTR, another major drug transporter, also has 11 drug–protein complex structures available in PDB and contains 39 amino acid variations. As a drug that binds to the plasma protein will remain in circulation and only serve as a reservoir before it is unbound from plasma protein, the binding affinity between the drug and plasma protein directly influences the biological half-life of the drug. In addition, many proteins that are involved in drug metabolic pathways can bind to multiple drugs and several structures are available. Cytochrome P450 3A4 (CYP3A4), the major enzyme that modifies and detoxifies drugs by oxidation reaction in the liver, has 5 drug-bound protein structures in PDB and contains 28 variations within drug binding sites. Figure 6a shows the binding of different drugs (bromocriptine, erythromycin, metyrapone, ritonavir, progesterone) and CYP3A4. The drugs are shown in magenta color in the protein pocket bound to the amino acid residues that are shown in cyan and blue color. The blue colors mean the amino acids are altered because of nsSNV. The heme of CYP3A4 is marked in yellow. Panels on the top of Figure 6a show that bromocriptine and erythromycin chemically bind to heme, whereas panels on the bottom show that metyrapone, ritonavir and progesterone occupy the pocket without directly binding to heme. AKR1C2 and AKR1C3, the major enzymes to detoxify drugs by conjugation reaction, have 10 and 8 target proteins available in PDB and contain 15 and 20 mutations within protein binding sites, respectively. For these key proteins, we believe that the identified mutations on the drug binding sites are one of the determinants for drug binding affinity and, consequently, the drug efficacy and toxicity. As an example, Figure 6b shows the superposition of energy minimized structures for the wild-type (WT) CA2 bound to lacosamide (PDB: 3IEO) and the mutated models (N67K, Q92P and F131L) bound to the same drug. Only the WT (orange color) and mutated residues (cyan color) with the bound drugs (yellow: WT and green color: mutant) are shown. Because of the mutations, the drug has to shift from the WT location to another location in the binding pocket that causes an unfavorable binding energy by 4kcal compared with the WT binding. Using the mutational profile, it is possible to understand the difference between a drug bound to WT and the mutant protein. Therefore, the identification of their mutational profile contributes to their pharmacogenetic and pharmacogenomic characterizations.
In terms of orthlogs of non-antineoplastic drug target proteins, 133 orthologs have been identified with 43 drugs involved.
Within the 36 proteins that bind to 25 antineoplastic drugs, 5 proteins have at least one paralog that have the full binding site conserved. Carbonic anhydrase 5A (CA5A), carbonic anhydrase 5B (CA5B), carbonic anhydrase 7 (CA7) and carbonic anhydrase 13 (CA13) were identified as four paralogs of CA2 that is a target protein of hydroxyurea, an antineoplastic agent. Figure 7 shows a model of CA2 and its paralogs CA13 and CA7 binding to hydroxyurea (magenta color). The three-dimensional structure shows the consistent binding between side chains of CA2, CA13, CA7 and hydroxyurea and suggests the potential binding affinity between CA5B, CA5A and hydroxyurea. In the DrugBank database, the only listed target for hydroxyurea is ribonucleodise-diphosphate reductase large subunit (RRM1).73 In PharmaGKB, the only available target protein of hydroxyurea is mitogen-activated protein kinase kinase kinase 5 (MAP3K5). Therefore, the data we provide here are complimentary to what is available on DrugBank, PharmaGKB and other resources.
In terms of 221 non-antineoplastic drug target proteins, we identified at least one paralog for 42 of them that have the drug binding site conserved. CA2 is also a target for non-antineoplastic drugs. It is part of the CA family and exists along with several of its family members ubiquitously across various tissues.74 It is a well-known target for diuretic and glaucoma drugs, such as ethoxzolamide, acetazolamide and methazolamide.75 However, its existence in vascular tissue also causes hypertension as a side effect by a number of drugs. Our study shows that CA2 has paralogs CA13, CA5A, CA5B and CA7 that have conserved drug binding sites. Search for druggable paralogs may lead to structurally modified alternative drugs that target its more tissue-selective paralogs and show similar pharmacological effect but less side effect. Roflumilast, another example, shows anti-inflammatory and anti-modulatory effects in the pulmonary system by selectively inhibiting phosphodiesterase-4.76 Amino acid–drug binding data retrieved from PDB show the binding between cAMP-specific 3',5'-cyclic phosphodiesterase 4D (PDE4D) and Roflumilast. Our paralogs survey shows that paralogs PDE4A, PDE4B and PDE4C are potential alternative target proteins.
For additional homolog discovery, a total number of 142 human homologs of non-human mammalian proteins were identified as potential drug targets (Supplementary Table S4). Among them, 19 are targeted by antineoplastic drugs and the rest of 123 human homologs are targeted by non-antineoplastic drugs.
To further estimate the SNVs that may structurally affect protein–drug binding, we mapped msSNVs affecting within two and five amino acids upstream and downstream region of all the binding sites for each protein–drug pair. The results of 407 protein–drug pairs are shown in Supplementary Table S5. Generally, the number of binding site affecting SNVs increases proportionally as the binding region extends from two amino acids to five amino acids. Table 2 shows the top 10 protein–drug pairs and their somatic mutations and germline mutations within 2 and 5 amino acids upstream and downstream of binding sites. For germline mutation in Table 2, EGFR and tyrosine-protein kinase ABL1 appear to be two cancer drug target proteins with large number of somatic mutations within the drug binding region. The high variability of these two proteins shows potential risk in therapeutics when they are targeted for cancer treatment.
In addition to SNVs, gene expression changes in specific cancer types can also affect therapeutic effect.15 In the current protein–drug pairs that contain somatic mutations within their binding site, we identified 85 differential gene expression changes associated with specific cancer types (see Supplementary Table S6).48 Table 3 shows differential gene expression profiles for antineoplastic drug targets in specific cancer types. For example, KIT is significantly underexpressed in lung squamous cell carcinoma and colon adenocarcinoma. When sunitinib is used with normal dose in targeted therapy in cancer treatment, the therapeutic effect could be diminished because of the lower expression of target protein and somatic mutations within the binding sites between sunitinib and its target protein.
DrugVar (http://hive.biochemistry.gwu.edu/tools/drugvar) scan allows users to upload a comma/tab-delimited annotation files such as those produced by ANNOVAR77 variation annotation pipeline and see whether any of the variations can alter any of the drug binding sites. The interface also supports both protein-based (UniProtKB accession, PDB ID) and drug-based (DrugBank ID, CAS ID, CID ID) searches.
There are many more drug–protein interactions than what is available in structural databases such as PDB. It is expected that over time, because of structural genomics initiatives, more structural data will be available. We will update DrugVar at least once every year to capture all such new data.
This project is partially funded by NCI/EDRN contract no. 156620 to RM. We want to thank K Smith, S Goldweber and E Thompson for useful comments.
Supplementary Information accompanies the paper on The Pharmacogenomics Journal website (http://www.nature.com/tpj)
The authors declare no conflict of interest.