bSiteFinder, an improved protein-binding sites prediction server based on structural alignment: more accurate and less time-consuming

Motivation Protein-binding sites prediction lays a foundation for functional annotation of protein and structure-based drug design. As the number of available protein structures increases, structural alignment based algorithm becomes the dominant approach for protein-binding sites prediction. However, the present algorithms underutilize the ever increasing numbers of three-dimensional protein–ligand complex structures (bound protein), and it could be improved on the process of alignment, selection of templates and clustering of template. Herein, we built so far the largest database of bound templates with stringent quality control. And on this basis, bSiteFinder as a protein-binding sites prediction server was developed. Results By introducing Homology Indexing, Chain Length Indexing, Stability of Complex and Optimized Multiple-Templates Clustering into our algorithm, the efficiency of our server has been significantly improved. Further, the accuracy was approximately 2–10 % higher than that of other algorithms for the test with either bound dataset or unbound dataset. For 210 bound dataset, bSiteFinder achieved high accuracies up to 94.8 % (MCC 0.95). For another 48 bound/unbound dataset, bSiteFinder achieved high accuracies up to 93.8 % for bound proteins (MCC 0.95) and 85.4 % for unbound proteins (MCC 0.72). Our bSiteFinder server is freely available at http://binfo.shmtu.edu.cn/bsitefinder/, and the source code is provided at the methods page. Conclusion An online bSiteFinder server is freely available at http://binfo.shmtu.edu.cn/bsitefinder/. Our work lays a foundation for functional annotation of protein and structure-based drug design. With ever increasing numbers of three-dimensional protein–ligand complex structures, our server should be more accurate and less time-consuming.Graphical Abstract bSiteFinder (http://binfo.shmtu.edu.cn/bsitefinder/) as a protein-binding sites prediction server was developed based on the largest database of bound templates so far with stringent quality control. By introducing Homology Indexing, Chain Length Indexing, Stability of Complex and Optimized Multiple-Templates Clustering into our algorithm, the efficiency of our server have been significantly improved. What’s more, the accuracy was approximately 2–10 % higher than that of other algorithms for the test with either bound dataset or unbound dataset


Background
Most biological processes involve the interaction of ligands with proteins. Functional characterization of ligand-binding sites of proteins is a key issue in understanding those biological processes [1][2][3][4]. In addition, identifying the location of protein-binding sites is a vital first step in structure-based drug design [5][6][7][8]. However, functional characterization of proteins through experimental method is a labor intensive and time-consuming process. A computational tool to predict the functional binding sites in a protein is therefore of practical importance.
To date, a variety of computational methods have been developed for protein-binding sites prediction, which can be divided into four categories: geometry based methods [9][10][11][12][13][14], energy based methods [15,16], alignment based methods [17][18][19][20] and other miscellaneous methods [21][22][23]. Alignment based methods can be further divided into sequence alignment based and structural alignment based methods. Recently, increasing structural genomics projects have led to the exponential growth of the number of available protein structures. As a consequence, structural alignment based methods exceeded other methods due to its more efficient and more accurate performance.
In 1996, Lichtarge et al. [17] developed the first structural alignment based algorithm for protein-binding sites prediction, entitled evolutionary trace method (ET method). It is based on the extraction of functionally important residues from sequence conservation patterns in homologous proteins, and on their mapping onto the protein surface to generate clusters identifying functional interfaces. In 2007, Brylinski and Skolnick developed a popular structural alignment method called FINDSITE [18]. For a given target sequence, FINDSITE identifies ligand-bound template structures from a set of distantly homologous proteins recognized by the PROS-PECTOR_3 threading approach and superposes them onto the target's structure using the TM-align structural alignment algorithm. Binding pockets are identified by the spatial clustering of the center of mass of templatebound ligands that are subsequently ranked by the number of binding ligands. In 2009, Oh et al. [24] developed LEE, a two-stage template-based ligand binding site prediction method, where templates are used first for protein 3D modeling and then for binding site prediction by structural clustering of ligand-containing templates to the predicted 3D model. Later in 2010, Wass et al. [25] described a new method called 3DligandSite. Structures similar to the query are identified by using MAM-MOTH [26] against a library of protein structures with bound ligands. The structural based alignment of the similar structures and the query superposes ligands onto the query structures. After filtering, the top 25 ligands are retained for analysis and further clustering. In 2012, another comparative approach called COFACTOR was proposed by Zhang group [19]. COFACTOR recognizes functional sites of protein-ligand interactions using lowresolution protein structural models, based on a globalto-local sequence and structural comparison algorithm. The major advantage of COFACTOR over the existing methods is the optimal combination of global and local structural comparisons for identifying protein-binding sites. But, the global comparison can be distracted by structural variations in the regions far away from the binding pockets; meanwhile the local comparison has a high false positive rate since the number of residues involved is too small. Later in 2013, Zhang group published another structural alignment based algorithm, TM-SITE [20]. Different from COFACTOR, TM-SITE compares the structures of a subsequence from the first binding residue to the last binding residue (called SSFL) on the query and template proteins, which solve the problems of global-to-local structural comparison algorithm. These methods provide us valuable choices to predict the binding sites. However, their performance needs to be improved for lack of accuracy or time-efficiency or both since the structural information of protein-ligand complexes (bound protein) are underutilized.
Herein, we built so far the largest database of bound templates with stringent quality control. And on this basis, Stability of Complex as a new criterion and Optimized Multiple-Templates Clustering algorithm are introduced to improve the accuracy. Meanwhile, Homology Indexing and Chain Length Indexing are used to accelerate the efficiency of the structural alignment. Finally, we presented a user friendly protein-binding sites prediction web server (bSiteFinder), at http://binfo. shmtu.edu.cn/bsitefinder/.

Rules of five
The protein data in PDB database are filtered through the rules below: 1. The macromolecule type is protein, no DNA and RNA. 2. Experiment method is set to X-ray. 3. X-ray resolution is between 0 and 3.0. 4. Has free ligands = yes. 5. Sequence length is over 20.

Number of ligand atoms
In the process of building databases, which database a protein finally falls into depends on whether it contains ligands and whether these ligands have enough atoms. For this reason, ligands identification, which is judged by the rules mentioned below, plays a key role. Every HETATM residue is recognized through HET records from the header of PDB files. Notably, some of the residues are modified on normal chains, which are not counted as true ligands because of their present in the MODRES records. Hence, the selected ligands only come from HET records excluding MODRES ones. Water molecule is included in HETATM but not regarded as a ligand. Analyzing the data, we define that a ligand should possess 6 or more atoms as a basic rule to identify a ligand.

Stability of Complex
The binding site check criterion is using as the standard of judging the bound structure's stability. Only if any one of atoms of the ligand has a distance within 4 Å from the geometry center of the calculated binding site, the structure of complex is considered to be stable.

Homology Indexing
Homology Indexing is implemented by using SCOPe, version 2.03 [27]. First, a four-digit classification number is searched based on PDB ID and CHAIN ID of the query chain. After that, all the protein chains with the same classification number are obtained and used to constitute the template database for subsequent structural alignment.

Chain Length Indexing
Only the chains, which have length difference with query chain less than 30 %, are used as candidates for subsequent structural alignment.

Structural alignment
The structural alignment between query and templates in bSiteFinder is implemented by using Combinatorial Extension (CE) algorithm, which is provided by Biojava [28]. Different from traditional dynamic programming algorithm and Monte Carlo algorithm, CE algorithm defines continuous residues in the sequence as aligned fragment pairs (AFPs), which is used in local alignment between query and template. Finally, the optimized alignment results are obtained by expanding or abandoning the local AFPs.

Optimized Multiple-Templates Clustering
After structural alignment, template will be mapped to query. Then, the templates which meet the requirement of Stability of Complex are ranked according to the similarity with query chain, and ligands of the top 20 templates at most will be picked out. After 20 times of structural alignments, all the ligands in templates will be mapped to the query. Further, these ligands are clustered into different clusters. The number of ligand geometric centers, which have a distance less than 3 Å from the certain ligand geometric center, is counted for each ligand. After that, the ligand with the largest number is defined as the center of the Top1 binding site (Fig. 1). Then, this ligand and all the other ligands within 3 Å are removed for searching the centers of the Top2 and Top3 binding site in the same way.

Detection of binding sites
On the condition that protein chains have ligands, we define all residues within the distance of 8 Å from ligands as the components of the binding site. On the condition that binding site is detected by doing structural alignment with templates, all residues within the distance of 10 Å from mapped ligands are defined as the components of the binding site. It should be noted that if the bound proteins' stabilities did not pass the evaluation of Stability of Complex, the bound proteins would be treated as unbound proteins with original ligands removed.

Test and evaluation methods
For comparing with other binding site prediction algorithms, two widespread adopted datasets from LIG-SITEcsc [29] were used for testing our algorithm with the same criteria of evaluating the accuracy of binding site prediction. The first test set contained 210 proteins with ligands (bound dataset). At the suggestion of RCSB, Fig. 1 Workflow of Optimized Multiple-Templates Clustering. Template (b) is mapped to query (a) by structural alignment to form query-template complex (c). Then, the template chain will be removed, and the ligand will be retained (d). After 20 times of structural alignments, the ligands in templates will be mapped to the query (e). The number of ligand geometric centers, which have a distance less than 3 Å from the certain ligand geometric center, is counted for each ligand (f). The ligand with the largest number is defined as the center of the Top1 binding site (g) protein 1B6N was replaced by 1Z1H. The second test set contained 48 proteins with/without ligands (bound/ unbound dataset).
Here, the accuracy and Matthews Correlation Coefficient (MCC) [30] were both used to evaluate our algorithm.

Accuracy
A widely accepted verification method [13] was used. For bound protein, if the protein-ligand's stability has passed the evaluation of Stability of Complex, the accuracy is 100 %. If the protein-ligand's stability did not pass the evaluation of Stability of Complex, the original ligands of bound protein will be removed and in this situation, the bound protein will be regarded as unbound protein and may have a lower accuracy.
For unbound proteins, if the geometric center of a binding site has a distance within 4 Å from any one of the atoms of the predicted ligands, this binding site is regarded as a correctly predicted binding site. Otherwise, this binding site is regarded as an incorrectly predicted binding site.

MCC
Another evaluation index, MCC, was also used to evaluate the accuracy of binding site prediction. For each protein chain, all the residues were divided into four categories: TP: correctly predicted binding site residues; TN: correctly predicted nonbinding site residues; FP: incorrectly predicted as binding site residues; and FN: incorrectly predicted as nonbinding site residues. MCC scores are defined as: For bound proteins that passed the evaluation of Stability of Complex, the MCC is 1. Otherwise, the bound proteins was regarded as unbound proteins and MCC would be lower than 1.
For unbound proteins, the structural alignment between query and template is implemented to map the ligands in bound proteins to the unbound proteins. Then, the mapped pseudo ligands were used to detect the binding site as describe in "Detection of Binding Sites". To evaluate our methods, we divided the residues of query chains into residues of predicted binding site (Res-BS-Pre) and residues of predicted non-binding site (Res-NBS-Pre). At the same time, we also define residues of experimental binding site as Res-BS-Exp and residues of experimental non-binding site as Res-NBS-Exp according to the original ligands of query chains. Therefore,

Create template database
Our algorithm will maximize the information of bound proteins. Herein, we built so far the largest database of bound templates from PDB database with stringent quality control. Figure 2 shows the workflow of creating template database, which include four steps as follow: (1)

Workflow of binding sites detection
When a query protein is submitted by user for binding site prediction, it will be firstly divided into chains. After that, the prediction will be done for each chain. Figure 3 shows the workflow of binding sites detection. Each protein chain will be processed by following steps: 1. Binding sites prediction of high quality bound protein (Part 1) Detection of Binding Sites is employed for binding site detection, when the protein chains meet the requirement of Number of Ligand Atoms and Stability of Complex. Otherwise, enter the following process.

Binding sites prediction of unbound protein with bound templates of same Homology Indexing (Part 2)
If the query chain has a four-digit classification number in SCOPe and has bound template with the same Homology Indexing in template database, the binding site of this query chain will be detected as the following procedure. First, structural alignments between query chain and templates will be done, and the top 20 bound templates which are the most similar to the query will be selected subsequently. The locations of ligands are detected by mapping the ligands in templates to the query, and then the optimization of binding sites was following by using the new developed Optimized Multiple-Templates Clustering method. Finally, Detection of Binding Sites will be employed for binding site detection. Otherwise, enter the following process.  Workflow of binding sites detection. Each protein chain submitted would be processed successively by following steps: 1 Binding sites prediction of high quality bound protein (Part 1), or enter the following process. 2 Binding sites prediction of unbound protein with bound templates of same Homology Indexing (Part 2), or enter the following process. 3 Binding sites prediction of unbound protein with bound templates of Chain Length Indexing (Part 3). Any protein chains submitted into our system could receive the results of binding sites via efficient computation

Binding sites prediction of unbound protein with bound templates of Chain Length Indexing (Part 3)
If the query chain has no satisfactory homologous bound template, the binding site of this query chain will be detected as the following procedure. Chain Length Indexing will be employed to search the bound templates, which have difference with query chain less than 30 % in length, in template database. Then enter the process as the description above (Part 2 of "Workflow of binding sites detection") with top 20 most similar bound templates. Any protein chains submitted into our system could receive the results of binding sites via efficient computation.

Performance of our algorithm and its comparison with others
Two widely adopted datasets including 210 bound and 48 bound/unbound dataset [29] were used for testing our algorithm, and the results are shown in Tables 1 and 2. The accuracy of our algorithm is approximately 2-10 % higher than that of other algorithms for the test with either bound or unbound datasets. In addition, with size of the dataset increased, our algorithm exhibited even more advantage over others regarding accuracy (The accuracy differences between our algorithm and the second highest algorithm in the Top1 increase from 2.4 % with 48 unbound dataset to 11.8 % with 210 unbound dataset).
For bound chain (such as PDB ID: 5p2p, CHAIN ID: A), the binding site is composed of residues within 8 Å from the ligand (Fig. 4a). For unbound chain (such as PDB ID: 3p2p, CHAIN ID: A), unlike bound chain, the binding site is detected with the aid of templates (PDB ID: 1oxr, CHAIN ID: A). First, the ligand in template is mapped to unbound chain. Then the binding site is composed of residues within 10 Å from the ligand (Fig. 4b). See Method part for details.   Table 3 shows the alignment frequency between templates and the query from the 48 unbound dataset after Homology Indexing is used. Without Homology Indexing, 48 unbound dataset should be aligned with each of chains in template database, which means that there are 48 × 101,315 time-consuming structural alignments needed to be done. But, with the Homology Indexing introduced, it can be reduced to 25,127 structural alignments, which only account for only 0.5 % of that without Homology Indexing. It's worth noting that alignment frequencies, in Table 3, reach hundreds or even thousands in practical, which may be due to the uneven distribution of different protein families in template database at present.

Stability of Complex
Examining the bound chain structures in PDB database, it is observed that ligands do not always have a stable binding with protein chains at binding site, such as PDB ID: 2j22, CHAIN ID: A (Fig. 5). For this kind of bound structures, binding sites could not be computed directly based on their ligands. Thus, Stability of Complex is introduced into our algorithm to avoid these situations.
Looking for similar templates by structural alignments is needed for unbound chains which have no ligands to compute the binding site. In the process of structural alignment and ligand mapping successively, ligand in template may not have a stable bind with unbound chain (Fig. 6a, b). Likewise, Stability of Complex is employed here to decide whether ligand from template and unbound chain can form a new stable bound structure.
Similarly, Stability of Complex is introduced to build a template database (see details in Fig. 2), which reduced the number of bound structures from 117,823 to 101,315 with 14 % structures removed. Not only improved the quality of template database, this operation also reduced the number of time-consuming structural alignments.

An Optimized Multiple-Templates Clustering method
Similar to FINDSITE [31], 3DLigandSite [25] and COFACTER [19], the prediction accuracy of our algorithm is improved by Optimized Multiple-Templates Clustering. However, in other works, the cluster number is required in previous algorithms, which actually could not be obtained before computing. In addition, the distances between ligands in each cluster have no reasonable physical meaning. In our algorithm, this deficiency is overcome by defining a new constraint, which restrict that the distances between geometric centers of all the ligands (for one binding site) in the same cluster should be less than a certain threshold (cluster radius). Ligands in multiple templates could be clustered automatically following the constraint with reasonable physical meaning, and there has no need to estimate cluster number before clustering.
Considering the space complexity of bound structure, cluster radius to be used is optimized based on test set. For 48 unbound dataset, threshold is set from 1.0 to 8.0 Å to compute the accuracy of the Top1 and Top3. Table 6 shows the accuracy computed with different cluster radius, and the accuracies of the Top1 range from 72.3 to 85.4 %. It's worth noting that the accuracy of our algorithm with any cluster radius is higher than that of other algorithms (Tables 2, 6).
Result in Table 6 indicates that the Top1 and Top3 have highest prediction accuracies with 48 unbound dataset, when cluster radius is set to 3.0 Å. Thus, 3.0 Å is set as the default parameter by bSiteFinder in Optimized Multiple-Templates Clustering.

Conclusions
bSiteFinder as a protein-binding sites prediction server was developed based on the largest database of bound templates so far with stringent quality control. Each protein chain submitted would be processed by following steps: (1) Binding sites prediction of high quality bound protein; (2) Binding sites prediction of unbound protein with bound templates of same Homology Indexing; (3) Binding sites prediction of unbound protein with bound templates of Chain Length Indexing. Any protein chain submitted could receive the results of binding sites via efficient computation. By introducing Homology Indexing, Chain Length Indexing, Stability of Complex and Optimized Multiple-Templates Clustering into our algorithm, the efficiency of our server have been significantly improved. What's more, the accuracy was approximately 2-10 % higher than that of other algorithms for the test with either bound dataset or unbound dataset. For 210 bound dataset, bSiteFinder achieved high accuracies up to 94.8 % (MCC 0.95). For another 48 bound/ unbound dataset, bSiteFinder achieved high accuracies up to 93.8 % for bound proteins (MCC 0.95) and 85.4 % Fig. 6 a Unbound chain (PDB ID: 1bbs, CHAIN ID: A, blue) and related appropriate template (PDB ID: 1hrn, CHAIN ID: B, yellow). After mapping the ligand (03D, red) in template to unbound chain, a new stable bound structure is formed with the tightly binding between the ligand and unbound chain. The top 20 templates at most ranked according to the similarity would be subsequently clustered. b Unbound chain (PDB ID: 1bbs, CHAIN ID: A, blue) and related appropriate template (PDB ID: 3g6z, CHAIN ID: A, yellow). After mapping the ligand (NAG, red) in template to unbound chain, a new stable bound structure could not be formed. The reason is that there are more residues (see the red circle) in template than unbound chain which have a close connection with the ligand for unbound proteins (MCC 0.72). An online bSiteFinder server is freely available at http://binfo.shmtu.edu.cn/ bsitefinder/, and the source code is provided at the methods page. Our work lays a foundation for functional annotation of protein and structure-based drug design.
With ever increasing numbers of three-dimensional protein-ligand complex structures, our server should be more accurate and less time-consuming.