mutLBSgeneDB: mutated ligand binding site gene DataBase.

Mutations at the ligand binding sites (LBSs) can influence protein structure stability, binding affinity with small molecules, and drug resistance in cancer patients. Our recent analysis revealed that ligand binding residues had a significantly higher mutation rate than other parts of the protein. Here, we built mutLBSgeneDB (mutated Ligand Binding Site gene DataBase) available at http://zhaobioinfo.org/mutLBSgeneDB We collected and curated over 2300 genes (mutLBSgenes) having ∼12 000 somatic mutations at ∼10 000 LBSs across 16 cancer types and selected 744 drug targetable genes (targetable_mutLBSgenes) by incorporating kinases, transcription factors, pharmacological genes, and cancer driver genes. We analyzed LBS mutation information, differential gene expression network, drug response correlation with gene expression, and protein stability changes for all mutLBSgenes using integrated genetic, genomic, transcriptomic, proteomic, network and functional information. We calculated and compared the binding affinities of 20 carefully selected genes with their drugs in wild type and mutant forms. mutLBSgeneDB provides a user-friendly web interface for searching and browsing through seven categories of annotations: Gene summary, Mutated information, Protein structure related information, Differential gene expression and gene-gene network, Phenotype information, Pharmacological information, and Conservation information. mutLBSgeneDB provides a useful resource for functional genomics, protein structure, drug and disease research communities.

INTRODUCTION

Molecular recognition plays a fundamental role in all biological processes (1). Mutation-induced conformational change and induced fit with the ligand are the key factors of protein–ligand interactions in cancer cells (2,3). Point mutations at spatially distinct sites lead to conformational changes and exert hinge effects (4). Some point mutations at ligand binding sites may dramatically change the binding affinities of the ligands (5,6). Studies also reported that mutations at ligand binding sites could link to the resistance to small molecule drugs in patient care (7,8). Recently, we also found a significantly higher mutation rate at ligand binding residues than in other parts of the protein sequence across 16 cancer types (9). Therefore, comprehensive annotations of all ligand binding site mutations in pan-cancer will allow for investigators to better understand cancer mechanisms and identify targetable mutations at ligand binding sites.

Many researchers have identified mutation-induced molecular modifications in ligand-protein interactions. For example, mutations in epidermal growth factor receptor (EGFR) in glioblastoma increased ligand binding affinity for EGF (10). A point mutation in neuraminidase 1 gene (NEU proto-oncogene) conferred high ligand binding affinity (6). Moreover, a few studies reported the roles of ligand binding domain mutations. The association between the ligand binding sites and disease related mutations in the type I collagen was observed (11) and the ligand-binding-domain mutations of androgen receptor (AR) gene led to the disruption of interaction between N- and C-terminal domains (12). Recently, several studies showed that ligand binding site mutations could lead to drug resistance. For example, ligand-binding domain mutations in estrogen receptor 1 gene (ESR1) were found in hormone-resistant breast cancer (7). Two major ligand binding site mutations in ESR1 can confer partial resistance to the currently available endocrine treatments (13). Consequently, the cancer and drug research community has recognized the importance of ligand binding site mutations and called for systematic and comprehensive analyses of genes with ligand binding site mutations (14), which are still largely not done yet despite the exponential growth of cancer and other biomedical data recently.

This paper introduces mutLBSgeneDB (mutated Ligand Binding Site gene DataBase), the web interface, and its applications. As the first database encompassing all human ligand binding site mutations with bioinformatics analyses, it provides unique and useful information for functional genomics, protein structure, disease and drug research communities.

DATABASE OVERVIEW

mutLBSgeneDB contains over 2300 genes with ligand binding site mutations that are annotated with seven categories (Figure 1). (i) Gene summary category provides basic gene information with diverse hyperlinks and the literature evidence in ligand binding site mutations for each gene. (ii) Ligand binding site mutation information category presents detailed information of somatic mutations that occur at the ligand binding sites only. The current version of mutLBSgeneDB includes 11 873 non-synonymous mutations at 10 108 ligand binding sites that were extracted from The Cancer Genome Atlas (TCGA) (15) and a semi-manually curated database for biologically relevant ligand-protein interactions (BioLiP) (16). (iii) Protein structure related information category shows relative stability of proteins encoded by all mutLBSgenes and ligand binding affinity changes with their drugs after the occurrence of mutation at the ligand binding site of carefully selected 20 genes. (iv) Differential gene expression and gene-gene network category shows expressional differences between mutated and non-mutated samples based on co-expressed protein interaction network (CePIN). (v) Phenotype information category includes disease related information on the genetic and mutation level. (vi) Pharmacology information category provides heat map using the top 20 correlated drugs between mutLBSgene expressions and 138 anti-cancer drug responses across 790 cell-lines from Cancer Cell Line Encyclopedia (CCLE) (17). mutLBSgeneDB shows the druggable features of mutLBSgenes covering a total of 1324 drugs from DrugBank. (vii) Conservation information category offers conserved sequences for each ligand binding site residue across eight species.

Figure 1.

Overview of mutLBSgeneDB. mutLBSgeneDB is composed of seven annotation categories that cover the genotypic, transcriptomic, proteomic, phenotypic, pharmacological data and network analyses of 2372 mutLBSgenes with manual curation of 744 targetable_mutLBSgenes.

Table 1 summarizes the statistics of 2372 genes having mutations at ligand binding sites (mutLBSgenes) and 744 drug targetable mutLBSgenes (targetable_mutLBSgenes) for each annotation category. mutLBSgeneDB can be used to explore and predict cancerous features and possible drug repurposing. All aforementioned entries and annotation data are available for browsing and searching on the mutLBSgenes web site (http://zhaobioinfo.org/mutLBSgeneDB).

Table 1.

Annotation entry statistics for mutLBSgenes and targetable_mutLBSgenes

Data type	# Entries	# mutLBSgenesa	# targetable_mutLBSgenesb
		Total 2372 (%)	Total 744 (%)
Targetable genes	# genes
Human Kinomec	267	267 (11.3%)	267 (35.9%)
TRANSFACd	216	216 (9.1%)	216 (29.0%)
IUPHARe	579	579 (24.4%)	579 (77.8%)
Cancer driver genesf	179	179 (7.5%)	179 (24.1%)
Ligand binding site	# LBS
BioLiPg	10 108	2372 (100.0%)	744 (100.0%)
Mutation	# nsSNV
TCGAh	11 873	2372 (100.0%)	744 (100.0%)
Expression	# genes
TCGAi	20 502	2372 (100.0%)	744 (100.0%)
Expression with drug treatment	# genes
CCLEj	19 931	2372 (100.0%)	744 (100.0%)
Molecule	# molecules
DrugBankk	8206 drugs	865 (36.5%)	378 (50.8%)
UniProtl	2374 proteins	2372 (100.0%)	743 (99.9%)
BioLiPm	6108 ligands	1780 (75.0%)	572 (76.9%)
Phenotype	# phenotype
DisGeNetn	6761 disease ID	1449 (61.1%)	662 (85.5%)
ClinVaro	107 phenotype ID	80 (3.4%)	49 (6.6%)
Conservation	# LBS
MUSCLE resultsp	27 269	2371 (100.0%)	744 (100.0%)

aNumber of genes having ligand binding site mutations (mutLBSgenes).

bTargetable mutLBSgenes.

cKinases from The Human Kinome database.

dTranscription factors from TRANSFAC database.

eDrug target genes from IUPHAR database.

fSignificantly mutated genes in 18 TCGA cancer types.

gLigand binding sites from Ligand-protein binding database (BioLiP).

hSomatic non-synonymous single nucleotide variants (nsSNVs) from TCGA in 16 cancer types.

iExpression values from TCGA.

jAnti-cancer drug treated cell-line's gene expression data.

kmutLBSgenes related drug IDs from DrugBank database.

lProtein accession data from UniProt database.

mLigand binding to mutLBSgenes.

nGene-level disease annotation from DisGeNet database.

oMutation-level pathogenic information from ClinVar.

pConservation information across 8 species from MUSCLE.

DATA INTEGRATION

mutLBSgenes and targetable_mutLBSgenes

A total of 145 531 ligand–protein binding interactions for 2874 proteins of UniProt (18) were downloaded from BioLiP (data version January 2016) (16). Somatic point mutation data was retrieved from TCGA (March 2016). Mutations that occur on the direct protein–ligand binding site residue or on its immediate two flanking residues at both upstream and downstream sides were considered to be ligand binding site mutations. There were 4660, 4472 and 2741 nsSNVs located at the direct protein–ligand binding site residues, the immediate flanking residues (±1 aa), and the immediate second flanking residues (±2 aa), respectively. After this data processing, 2372 genes with 11 873 non-synonymous mutations at 10 108 ligand binding sites were obtained. Furthermore, 744 drug targetable mutLBSgenes were identified by incorporating kinase genes from the Human Kinome (19), transcription factors from TRANSFAC (20), all drug target genes from the concise guide to pharmacology (IUPHAR, International Union of Basic and Clinical pharmacology) (21), and cancer driver genes from cancer type specific, significantly mutated genes that we collected and curated previously (22). As a result, the targetable_mutLBSgenes are composed of 220 human kinases, 216 human transcription factors, 579 IUPHAR target genes, and 101 cancer-type specific significantly mutated genes (Supplementary Table S1). Ten common genes among the five gene sets were CREBBP, EP300, ESR1, EZH2, FGFR1, HDAC3, PGR, RXRA, SMARCA4 and SMO.

Manual curation of PubMed articles

For the 744 targetable_mutLBSgenes, a literature query of PubMed was performed in June 2016 using the search expression that applied to each mutLBSgene (using BRAF as an example here: ‘((BRAF[Title/Abstract]) AND mutation[Title/Abstract]) AND ligand[Title/Abstract])’. The abstracts of over 1000 articles were manually reviewed. We found literature evidence (138 articles) for 98 genes (∼4.0%) that support the role of these ligand binding site mutations in cancer or drug response. For the 301 genes annotated as kinase or cancer driver genes in mutLBSgenes, we added 3D structure images by searching the Protein Data Bank (PDB) (23). For the most recurrent mutation in each targetable_mutLBSgene, we added related clinical information from genetically informed cancer medicine (My Cancer Genome) (24). Using this curation method, we created a classification system for the genes in the database to show reliability. Class A consists of genes with literature evidence and is part of the targetable_mutLBSgenes. Class B consists of only targetable_mutLBSgenes without additional evidence. The remaining genes belong to Class C.

Expression data preparation

Gene expression data was downloaded from TCGA (January 2015). Normalized gene expression data from RNASeqV2 was extracted using the R package TCGA-Assembler (25). In addition, microarray gene expression data from over 790 cancer cell lines was extracted from CCLE (October 2012).

Co-expressed protein interaction network (CePIN)

The protein interaction network (PIN) reported in our previous study included 113 473 unique protein-protein interactions connecting 13 579 protein-coding genes (26,27). It was used in conjunction with the Pearson Correlation Coefficient (PCC) calculated for each gene-gene pair to build a CePIN. Co-expressed network figures were drawn using the igraph package in R (28). For each gene, the top 20 neighbor genes with the highest PCC values were kept in the network to reflect the genetic signals.

Gene-drug and gene-ligand interaction networks

Drug–target interactions (DTIs) were extracted from DrugBank (29) and the duplicated DTI pairs were excluded. All drugs were grouped using Anatomical Therapeutic Chemical (ATC) classification system codes (30). Two-dimensional chemical structure images of all drugs were generated using the chemical toolbox, OpenBabel (v2.3.1) (31). Ligand–target interactions (LTIs) were extracted from BioLiP and the duplicated LTI pairs were excluded.

Calculating drug binding affinity for top 20 ranked genes and their drugs

We selected 20 genes ranked by the following information: recurrences in the samples, targeted by the drugs of ‘approved and investigational’ status from DrugBank, and number of mutated ligand binding sites. We further selected the most studied drugs (2 or 3) for each of these genes by searching DrugBank and PubMed and downloaded PDB structure file of drugs and proteins. Using these data sets, we searched the drug binding affinities for these 20 genes. We downloaded the crystal structures of genes and three-dimensional structures of drugs from PDB and a free database of commercially available compounds for virtual screening (ZINC) in mol2 format (32). Individual mol2 files were converted into pdbqt files using the python script prepare_ligand4.py available in the Autodock Tools package (33). Using Autodock package, we computed the free energy and studied the docking. Lastly we searched the optimal fit of each drug into targets. The details about this method were described in previous studies (34,35).

Correlation between drug response and gene expression using CCLE data

Drug response data in 714 cell lines on 142 drugs was extracted from Genomics of Drug Sensitivity in Cancer (http://www.cancerrxgene.org/) (36) (October 2012). Pearson Correlation Coefficient (PCC) between drug response and gene expression was calculated for each drug-gene pair.

Conservation information

All sequences used in comparative alignment and specific positions of amino acid were downloaded from the Conserved Domain Database (CDD) of NCBI. Comparison of homologous sequences was obtained by using the multiple sequence alignment tool with high accuracy and high throughput (MUSCLE) (37).

Database architecture

The mutLBSgeneDB system is based on a three-tier architecture: client, server, and database. It includes a user-friendly web interface, Perl's DBI module, and MySQL database. This database was developed on MySQL 3.23 with the MyISAM storage engine.

WEB INTERFACE AND APPLICATIONS

Ligand binding site mutation information category

This category presents detailed information of non-synonymous mutations (i.e. nsSNVs) located at the ligand binding sites (Figure 2A) such as the lollipop-style plot showing the mutations that only occurred at the ligand binding sites, cancer type specific mutLBS table giving the sorted mutation frequency information across cancer types, and clinical information table showing the specific clinical information for the most frequently recurrent mutations. We obtained clinical information for 74 genes among 744 targetable_mutLBSgenes using My Cancer Genome. For example, the most frequently observed nsSNV of v-raf murine sarcoma viral oncogenes homolog B1 (BRAF) is the V600E driver mutation (BRAF), activates the MAPK pathway in 50% of melanoma patients (38). This mutation is located near the ligand binding site (A598) (39). The cancer type specific mutLBS table shows the consistent results with the previous studies that the two most frequently mutated cancer types of BRAF are thyroid carcinoma (THCA) and skin cutaneous melanoma (SKCM) (40–42). Another example of ESR1 shows the possible usage of cancer type specific mutLBS table for user to examine whether and how the mutLBS present in different cancer types (43) (Supplementary Figure S1).

Figure 2.

Ligand binding site mutation information category and protein structure related information category. (A) Ligand binding site mutation information category. The lollipop plot shows non-synonymous single nucleotide variants (nsSNVs) at the ligand binding sites only in amino acid sequence of protein B-Raf. Cancer type specific mutLBS table shows the frequency of each LBS mutation in each cancer type. (B) Protein structure related information category. The dot plot shows the relative protein stability change after occurrence of each mutation. Red dots represent ligand binding site mutations and grey dots for other nsSNVs. This category also shows the predicted protein structure of wild type and mutant protein with their drugs. For the top 20 selected genes, we also calculated the binding affinity.

To provide a weighted gene list, we sorted mutLBSgenes based on the number of mutated ligand binding site. Among the total 2372 mutLBSgenes, 1891 genes and 203 had more than two and ten ligand binding sites of mutation, respectively. Among these, the top 20 genes were ERGR, ABL, TP53, BCHE, CTNNB1, VHL, CSNK2A1, FOLH1, KRAS, THBS2, CD1D, F2, EP300, HGF, RUNX1T1, ABL1, AGO2, XDH, CD1B and CES1 (Supplementary Table S2). Gene set enrichment tests were performed for the 203 genes to infer the active pathways of mutLBSgenes (WebGestalt, adjusted P-value (i.e. q-value) <0.05, hypergeometric test followed by multiple test correction using Benjamini–Hochberg's method) (44). There were 40, 53 and 54 genes that were enriched in ‘negative regulation of cell death’ pathway, ‘response to endogenous stimulus’ pathway and ‘protein phosphorylation’ pathway with q-value 2.80 × 10−15, 3.58 × 10−16 and 6.35 × 10−15, respectively (Supplementary Table S3). From these pathways, we could infer that the top genes with the most ligand binding site mutations were significantly involved in tumorigenesis and phosphorylation.

Protein structure related information category

This category shows protein stability changes after occurrence of mutation at the ligand binding site for all proteins encoded by mutLBSgenes (Figure 2B). One study comparing protein mean square deviation (MSD) between wild type and mutant proteins of B-Raf V600E showed that even a single mutation of the protein could lead to much different molecular characteristics (45). To this end, we calculated the relative stability of protein structure after one ligand binding site mutation using MuPro1.1, a computational tool that predicts the protein stability changes for single-site mutation using support vector machines and neural networks (46). Our annotation results showed that five mutations (G466V, G466E, G466R, F468L and F595S) at three ligand binding sites of B-Raf may cause the change of protein structure toward a more stable form with a positive stability change value (Figure 2B).

Furthermore, to annotate mutation-induced modifications on protein-drug binding, we selected top 20 genes ranked by recurrences, targeted by the drugs of ‘approved and investigational’ status, and number of mutated ligand binding sites. These genes are BRAF, CDK2, CPS1, CYP11B2, CYP2B6, CYP2C19, CYP2C8, CYP3A4, EGFR, ERBB2, FGFR2, IDE, ITK, KEAP1, KIT, MET, RET, SULT1E1, VDR and XDH. Binding affinity (kcal/mol) between wild type and mutant proteins with their respective drugs (Supplementary Table S4) was calculated for each of these genes. For example, mutated protein encoded by BRAF has a lower free energy of binding to Vemurafenib, a FDA-approved BRAF kinase inhibitor in the treatment of melanoma, compared to other drugs such as Dabrafenib and Regorafenib (Figure 2B).

Differential gene expression and gene–gene network category

This category provides differential gene expression and gene-gene network of mutated versus wild type samples to show expressional differences in each cancer type. First, a violin plot shows the differential gene expression between mutated and wild type samples for each mutLBSgene (Supplementary Figure S2A). Second, a co-expressed protein interaction network (CePIN) was created by calculating the Pearson Correlation Coefficient (PCC) using gene expression values from TCGA. Using this method, we found significantly different gene-gene networks of BRAF between mutated and non-mutated samples in colon adenocarcinoma (COAD). Gene set enrichment analysis using CePIN gene elements of mutated samples showed that two co-expressed genes (BRAF and PAK2) enriched in ‘positive kinase activity’ pathway with q-value 0.0116. In the wild type samples, 10 co-expressed genes were enriched in ‘regulation of defense response’ and ‘activation of immune response’ pathways with q-value 9.86 × 10−06 and 2.04 × 10−05, respectively (Supplementary Figure S2A, Supplementary Table S5). This result suggested that the occurrence of point mutation at the ligand binding sites of BRAF in COAD may not activate the immune response; instead, protein kinases were activated for the tumorigenesis.

Phenotype information category

This category includes two phenotype information tables. The first table shows the related disease information for each gene retrieved from a database of gene-disease associations (DisGeNet) (47). As shown in Supplementary Figure 2B, the most studied disease name for BRAF is ‘Melanoma’ and less frequently studied diseases include ‘Colon cancer’, ‘Thyroid cancer’, and ‘Cardiofaciocutaneous syndrome’. This is consistent with the previous finding that BRAF mutations are present in approximately 50% of melanoma, 60% of thyroid, and 10% of colorectal carcinomas and are less prevalent in other tumor types (48). On the other hand, the mutation-level pathogenic information table shows pathogenic mutation information from a public archive of relationships between sequence variation and human phenotype (ClinVar) (49).

Pharmacological information category

This category provides pharmacological information such as the correlation between drug response and gene expression and the network visualization of genes with their interacting small molecules (Supplementary Figure S3A). Each gene expression profile in cell lines treated with anti-cancer drugs summarizes the correlation between drug response and altered gene expression by integrating microarray gene expression data from the CCLE database. For example, the expression of BRAF was positively correlated with the treatment effect of drug NVP-TAE684. From the network and information table for the relating drugs and ligands of each mutLBSgene, user can retrieve more detailed information including drug structures. Overall, mutLBSgeneDB includes 1198 FDA-approved drugs targeting 961 mutLBSgenes (Supplementary Table S6).

Conservation information category

This category presents the homologous protein sequences of its flanking region for each ligand binding site obtained from MUSCLE to indicate if the ligand binding site is conserved among different species (Supplementary Figure S3B). For example, protein B-Raf has three mutated ligand binding sites: A481, A598, and C532, all of which are conserved in Homo sapiens (common name: human), Mus musculus (mouse), Gallus (chicken), Caenorhabditis elegans (roundworm), and Drosophila melanogaster (fly). In comparison, all ligand binding site mutations in human EGFR have shown conservation in mice, but not in chicken.

DISCUSSION AND FUTURE DIRECTION

This study introduces a unique resource, mutLBSgeneDB, for the systematic annotation of genes having ligand binding site mutations. To serve functional genomics, protein structure, and drug research communities and advance precision medicine research, we will continuously update mutLBSgeneDB in the following directions. (i) Update routinely by checking the new data of mutations and ligand binding sites from TCGA and BioLiP. (ii) Collect high-quality drug pharmacological data from high-throughput screening and drug resistance studies. (iii) Continue to collect articles on ligand binding site mutations. (iv) Add more protein-ligand 3D structures highlighting ligand binding site mutations with their drugs. (v) Collect and curate germline mutations at ligand binding sites and make the data interactive to somatic mutations. (vi) Perform additional integrative analysis by using other omics data such as methylation, microRNA, and proteomics data. mutLBSgeneDB will be useful for many investigators in functional genomics, protein structure, and drug and therapeutic research.