Structural neighboring property for identifying protein-protein binding sites

Background The protein-protein interaction plays a key role in the control of many biological functions, such as drug design and functional analysis. Determination of binding sites is widely applied in molecular biology research. Therefore, many efficient methods have been developed for identifying binding sites. In this paper, we calculate structural neighboring property through Voronoi diagram. Using 6,438 complexes, we study local biases of structural neighboring property on interface. Results We propose a novel statistical method to extract interacting residues, and interacting patches can be clustered as predicted interface residues. In addition, structural neighboring property can be adopted to construct a new energy function, for evaluating docking solutions. It includes new statistical property as well as existing energy items. Comparing to existing methods, our approach improves overall Fnat value by at least 3%. On Benchmark v4.0, our method has average Irmsd value of 3.31Å and overall Fnat value of 63%, which improves upon Irmsd of 3.89 Å and Fnat of 49% for ZRANK, and Irmsd of 3.99Å and Fnat of 46% for ClusPro. On the CAPRI targets, our method has average Irmsd value of 3.46 Å and overall Fnat value of 45%, which improves upon Irmsd of 4.18 Å and Fnat of 40% for ZRANK, and Irmsd of 5.12 Å and Fnat of 32% for ClusPro. Conclusions Experiments show that our method achieves better results than some state-of-the-art methods for identifying protein-protein binding sites, with the prediction quality improved in terms of CAPRI evaluation criteria.


Introduction
The protein-protein interaction plays a key role in many biological functions, such as drug design and functional analysis. Gaining insights of various binding abilities will deepen our understanding on interaction. Determination of binding sites is widely applied in molecular biology research. Therefore, many efficient methods [1,2] have been developed for identifying binding sites.
Some existing approaches are based on analyzing differences between interface residues and non-interface residues, through machine learning methods or statistical methods. They analyze different features, such as sequence and structural properties or physical attributes. ProMate [3] creates interface or non-interface sphere around each residue. The histograms of many features are statistically obtained from spheres in training proteins. The probability for each sphere of a testing protein can be estimated to be on interface or not. The interface spheres are clustered to identify binding sites. PPI-Pred [4] uses several features to build an SVM model on interface prediction. It generates an interacting patch and a non-interacting patch for each training protein. Seven features are extracted from all interacting and noninteracting patches to predict if a testing patch is an interacting patch. Li et al. [5] divide protein residues into four different classes, which are distinguished by percentage of their neighboring interface residues. The core-SVM model is built over eight features and used to compute whether a residue is a core interface residue. In PINUP [6], an empirical scoring function consists of interface propensity and residue conservation score for predicting binding sites. PINUP takes top scoring patches and ranks residues based on their occurrences in these patches, clustered as predicted interface residues. Burgoyne et al. [7] analyze clefts on surface, that are likely to be binding sites. They can be ranked according to sequence conservation and physical properties. Meta-servers have also been constructed to combine strengths of some existing approaches. The program called meta-PPISP [8] combines three individual servers, namely cons-PPISP, ProMate and PINUP; another program called metaPPI [9] combines five prediction methods, namely PPI-Pred, PINUP, PPISP, Pro-Mate, and Sppider.
In addition, several structural algorithms have also been used to identify binding sites, through analyzing surface structures. SiteEngine [10] recognizes surface regions of a testing protein that are similar to some known binding sites, using geometric hashing triangles. ProBiS [11] predicts interface residues by local surface structure alignment. It compares a testing protein to known binding sites, for detecting structurally similar residues. Ortuso et al. [12] define most relevant interaction areas, based on 3D maps. The GRID program is used to compute on known structural complexes.
Another kind of methods are to examine all possible poses of two protein subunits; that is, how subunits may dock. Docking methods based on fast Fourier transformation (FFT) [13], geometric surface matching [14], as well as intermolecular energy [15] have been proposed. ZRANK [16,17] combines an atom-based potential (IFACE) with five residue-based potentials for ranking docked conformations. It provides fast and accurate rescoring of ZDOCK models [18]. ClusPro [19] develops a fast algorithm for filtering docked conformations with good surface complementarity, and ranks them based on their properties. RosettaDock [20] constructs an energy function using van der Waals energies, orientationdependent hydrogen bonding, implicit Gaussian solvation, side-chain rotamer probabilities and a lowweighted electrostatics energy. HADDOCK [21] makes use of biochemical and biophysical interaction data, such as chemical shift perturbation data resulting from NMR titration experiments. Fernandez-Recio et al. [22] apply docking simulations and analyze interaction energy landscapes to identify interface residues. They use a global docking method based on multi-start energy optimization, and predict low-energy regions as binding sites.
Identifying of protein-protein interface depends on many features, such as sequence, structure, as well as other physicochemical properties. Hydrogen bonds and salt bridges are known to be essential in identifying binding specificity [23]. Most of binding sites are hydrophobic and conserved polar residues at specific locations [24]. Secondary structure composition analysis shows that neither helices nor b-sheets are dominantly populated on interface [25]. Several geometrical features such as weighted atomic packing density, relative surface area burial and weighted hydrophobicity are most effective features for predicting interface residues [26]. Some features only describe properties of current interacting residues, but cannot represent real situation well, thus are insufficient to predict binding sites with high accuracy.
In this paper, we analyze structural neighboring property on protein-protein interface, through Voronoi diagram. Using 6,438 complexes, we study local biases of structural neighboring property on interface. We propose a novel statistical method based on structural neighboring property to extract interacting residues, and interacting patches can be clustered as predicted interface residues. In addition, structural neighboring property can be adopted to limit the search space, for discovering native-like poses. Here, we construct an energy function to evaluate docking solutions, which includes new statistical property as well as existing energy items [27]. Finally, we use trained SVM models to further select best poses for each pair of input proteins.
Experiments show that our method achieves better results than some state-of-theart methods. Here, we use CAPRI evaluation criteria, I rmsd and F nat . Comparing to existing methods for identifying binding sites, our approach improves overall F nat value by at least 3%. On Benchmark v4.0, our method has average I rmsd value of 3.31Å and overall F nat value of 63%, which improves upon I rmsd of 3.89Å and F nat of 49% for ZRANK, and I rmsd of 3.99 Å and F nat of 46% for ClusPro. On CAPRI targets, our method has average I rmsd value of 3.46 Å and overall F nat value of 45%, which improves upon Irmsd of 4.18 Å and F nat of 40% for ZRANK, and I rmsd of Å and F nat of 32% for ClusPro.

Methods
In this paper, we calculate structural neighboring property on protein-protein interface, through Voronoi diagram. We propose a novel statistical method to extract interacting residues, and interacting patches can be clustered as predicted interface residues. In addition, structural neighboring property can be adopted to construct an energy function to evaluate docking solutions, which includes new statistical property as well as existing energy items.

Data set
To obtain statistical property on interface, we adopt a high quality, non-redundant experimental data set. We select 6,438 complexes from Protein Data Bank [28]; each complex consists of two or more subunits. These complexes are determined from X-ray data with resolution less than 2.2Å. Any two complexes share no more than 30% identity.
A complex may contain several subunits and multiple interfaces. Each interface in a complex occurs in a pair of subunits. Two residues between a pair of subunits are called interface residues, if any two atoms, one from each residue, interact. By interact, we mean distance between two heavy atoms is less than 6 Å.

Structural neighboring property
Most of those features only describe current interacting residues, but cannot represent real situation well, thus are insufficient to predict binding sites with high accuracy. Here, we develop a method to calculate structural neighboring property on interface, using Voronoi diagram.

Site features
The physicochemical features are used to characterize potential interacting residues. The most interesting features are described as follows.
• Hydrophobicity: a numerical hydrophobicity of an amino acid [29], • Electrostatic potential: the number of electrostatic charge in an amino acid [30], • Hydrogen bonds: the number of potential hydrogen bonds for all atoms in an amino acid [31].
We use Voronoi diagram to evaluate polygonal face area of each surface residue, and calculate structural neighboring property on these physicochemical features.

Polygonal face area
We use VLDP [32] for geometrically analyzing protein 3D structures, based on Voronoi Tessellation. Voronoi Tessellation is a partition of space into polyhedra, whereas Delaunay diagram builds a graph with vertices at atoms. These graphs define nearest neighbours for each atom of one protein. It calculates Delaunay diagram by using an optimized incremental algorithm. The weights can be interpreted as squared radius. The surface residues appear as a packing of polyhedra, that is necessary to have a reasonable Tessellation throughout entire system. VLDP can be used to evaluate residue contacts and residue volumes, defined as polygonal face area and polyhedral volume for each atom. In particular, contact area of two residues is sum of atomic interface areas on pairs of atoms; surface area of one residue is sum of surface areas exposed to solvent in this residue; total area of one residue is sum of all areas in this residue. We calculate structural neighboring property, based on polygonal face area.

Property function
Given a protein, structural neighboring property of one surface residue x is defined as follow: where p(x) is site feature of each residue, contact(x, y) >0 means that contact area between residues x and y is greater than zero, shows surface area of residue x divided by its total area, as shown in Figure 1.
We use a normal distribution F (x) to estimate probability of structural neighboring property on one side of interface. We also use a bivariate normal distribution F (x1, x2) [33] to estimate probability of structural neighboring property on both sides of interface. Given a pair of proteins, effective free energy of interacting residue pair can be calculated as: where R is a set of all residue pairs on interface.

Extracting interface residues
We propose a statistical method to extract interacting residues, and interacting patches can be clustered as predicted interface residues. The threshold value s th is used to harvest all possible residue pairs between two proteins. The residue pairs with S(x 1 , x 2 ) ≤ s th are called interacting residues.
Considering neighboring residues, we construct a sphere with a radius of 10Å for all interacting residues. The updated statistical property of each interacting residue pair is calculated as follows.
where dis(x, r) is 3D distance of two C a atoms in residues x and r. The residues r i are from protein having residue x 1 , and the residues r j are from protein having residue x 2 .
We rank interacting residues by using updated statistical property. Top interacting residues can be grouped into different regions. All interacting residues served as graph nodes, and each undirected edge is built when two nodes are within distance 10 Å. Strongly connected components are considered as interacting patches. One region, containing a very small number of interacting residues, indicates a weak signal and can be discarded. We cluster interacting patches as predicted interface residues.

Energy function for docking
Given two input proteins, our task is to find the protein-protein interface between them. In first step, we identify docking solutions of two subunits. It performs a large number of rigid transformations to enumerate poses. Top ranking poses are selected through a linear combination of energy items. In second step, we generate possible conformational changes of interface residues from their unbound states to bound states, based on multidimensional scaling method. In third step, we use trained SVM models to further select best poses for input proteins. Structural neighboring property can be effectively applied to identify docking solutions.
Here, we construct a new energy function to evaluate docking solutions, which includes new statistical property as well as existing energy items [27]. The following lists all energy items, and how they are computed: • Structural neighborhood energy is calculated by probability of structural neighboring property on interface.
• Dihedral angle energy is calculated by statistical analysis of dihedral angle frequency and correlation on interface [27].
• Amino acid energy is constructed by probabilities of interface residues.
• Side-chain atoms of interface residues are packed by SCWRL4 [34], and sidechain energy is extracted.
We use a linear combination of these energy items, referred to as initial energy function, to rank poses. The coefficient of each item is optimized by using a linear combination method in [35]. We output top 100 poses with lowest energy values. For conformational changed structures, our method calculates a set of possibly changed conformations of interfaces.
As in [27], we use a training set consisting of 79 complexes from Dockground [36] to produce 79 SVM models, one for each complex, based on these energy items. Finally, we use trained SVM models to further select best 10 poses with lowest energy values for two input proteins.

Assessment of interface prediction
According to CAPRI evaluation criteria [37], three evaluation measures are commonly used in identifying proteinprotein interface. A pair of residues on interface is considered to be in contact if any of their atoms are within 6 Å One is the fraction of native contacts F nat , defined as the number of correct residue-residue contacts in predicted complex divided by the number of contacts in native complex. The other is the fraction of non-native contacts F non−nat , defined as the number of incorrect residuesresidue contacts in predicted complex divided by the total number of contacts in that predicted complex. The third is root-mean-square deviation of interface I rmsd , defined as the rmsd value between all backbone atoms of interfaces in predicted structure and in native complex, after two interfaces are superimposed.
We also calculate P value for binding sites prediction. The calculation of P value should be probability of obtaining not less than n correctly predicted interface residues by randomly picking out N predicted interface residues. The probability that a random method obtains success in one trial is m M , where M is the number of all surface residues, and m is the number of correctly interface residues among them. Therefore, P value for binding sites prediction is given by

Results
In this section, we have done three kinds of experiments. First, we present statistical analysis of structural neighboring property on interface. Then, we compare our method to some existing methods, for identifying binding sites. The results show that our method performs better than other machine learning and statistical approaches. Finally, we examine docking solutions of our method on Benchmark v4.0 and CAPRI targets.
Experiments show that our method outperforms some state-of-the-art methods.

Statistical property
The physicochemical features are used to characterize potential interacting residues. The most interesting features are described by three values, as shown in Table 1.
We present statistical analysis of structural neighboring property on interface. The statistics are carried out on 6,438 complexes. For each feature, we model a bivariate normal distribution. First, we represent an assessment for hydrophobicity on interface, as shown in Figure 2(a). The cluster centered at (1.89, 2.21) can be obtained. The groups of more hydrophobic amino acids often appear on interface. Second, local bias preferences of electrostatic potential on interface are shown in Figure 2b. We observe probability distribution centered at (0.12, −0.09). Many interfaces usually involve a lot of neutral amino acids. Third, we analyze hydrogen bonds on interface, as shown in Figure 2c. The data contains one cluster centered at (8.37, 7.96). The interface residues contain several potential hydrogen bonds, and the number of potential hydrogen bonds on each side of interface must be very similar.
We further investigate whether changing value of s th help to improve interface prediction. Here, we calculate structural neighboring property on 100 interface residues and 100 non-interface residues, randomly extracted from Benchmark v4.0 [38]. The relationship between interface prediction and the value of s th is shown in Table 2. We can observe that when increasing the value of s th , the overall F nat value of prediction is improved; however, the overall F non−nat value of prediction also increases as well.

Binding sites prediction
Some existing methods use machine learning and statistical approaches to predict binding sites. The results show that our method performs better than other existing methods in binding sites prediction.

Comparison to Fernández-Recio's method
In this test, we compare the performance of our method to Fernández-Recio's method. The test data used by this method consists of 43 complexes [22]. The results are reported in Table 3. The overall F nat and F non−nat values for our method are 65% and 32%, respectively. Fernández-Recio method achieves overall F nat and F non−nat values of 62% and 60%, respectively.

Comparison to metaPPI, meta-PPISP and PPI-Pred
In this experiment, we compare our method to metaPPI, meta-PPISP and PPIPred. The test data consists of 41 complexes by metaPPI [9], divided into two categories: enzyme-inhibitor (EI) and others. The overall F nat and F non−nat values for each prediction method are reported in Table 4. The overall F nat values for our method, metaPPI, meta-PPISP and PPI-Pred achieve 62%, 28%, 38% and 38%, respectively. The overall F non−nat values for these four methods achieve 34%, 51%, 54% and 64%, respectively. Our method improves overall F nat value by at least 24%.

Comparison to ProMate and PINUP
Our method is compared to ProMate and PINUP. The test data is originally used by ProMate [3], including 57 unbound proteins and their complexes. The results are reported in Table 5. The overall Fnat values for our method, PINUP and ProMate achieve 61%, 42% and 13%, respectively. The overall F non−nat values for these three methods achieve 45%, 55% and 47%, respectively. Our method improves overall F nat value by at least 19%.

Comparison to core-SVM
We compare our method to core-SVM with 50 dimers [5].
The results are reported in Table 6. The overall F nat values for our method and core-SVM are 63% and 60%, respectively. The overall F non−nat values for these two methods are 36% and 46%, respectively. Our method improves overall F nat value by at least 3%.

Docking result
In this study, we compare our docking solutions with ZRANK [16,17] and external tool, FiberDock [39], specifically designed to handle conformation change after binding. We also compare our docking results with ClusPro [19]. For unbound-unbound docking, experiments show that our method significantly outperforms these existing docking approaches.

Training set
We consider 79 complexes from Dockground [36] as training set. In order to avoid over-fitting, we exclude complexes, which share more than 30 percent identity with cases in testing set. The average I rmsd value is 1.49Å, and the overall F nat and F non−nat values are 85% and 16%.
Evaluation on Benchmark v4.0 On Benchmark v4.0, the average I rmsd values predicted by our method, ZRANK+FiberDock and ClusPro are 3.31Å, 3.89Å and 3.99Å, respectively. The overall F nat values predicted by these three methods are 63%, 49% and 46%, respectively. The results are shown in Table 7.
The complexes are classified into three categories, according to the magnitude of conformational change after binding. In rigid-body group, the average I rmsd values predicted by our method, ZRANK and ClusPro are 2.89Å, 3.31Å and 3.33Å, respectively. The overall F nat values predicted by these three methods are 69%, 56% and 55%, respectively. In medium difficulty group, the average I rmsd values predicted by our method, ZRANK +FiberDock and ClusPro are 3.38Å, 4.46Å and 4.71Å, respectively. The overall F nat values predicted by these three methods are 59%, 39% and 30%, respectively. In difficulty group, the average I rmsd values predicted by our method, ZRANK+FiberDock and ClusPro are 5.41Å,   6.18Å and 6.53Å, respectively. The overall F nat values predicted by these three methods are 36%, 28% and 21%, respectively.

Evaluation on CAPRI
We evaluate docking solutions of our method, ZRANK and ClusPro on CAPRI targets. CAPRI [37] is a community-wide experiment to assess the capacity of docking methods. The average I rmsd values predicted by our method, ZRANK+FiberDock and ClusPro are 3.46Å, 4.18Å and 5.12Å, respectively. The overall Fnat values predicted by these three methods are 45%, 40% and 32%, respectively. The results are shown in Table 8.

Assessment of energy items
To assess effectiveness of energy items, we re-optimize coefficients in each case with only four of five items. We evaluate docking poses of 176 complexes on Benchmark v4.0, by leaving one energy item out. The results are shown in Table 9. The overall F nat value for case without E nb is 58.6%, for case without E pi is 60.5%, for case without E da is 60.2%, for case without E aa is 59.3%, and for case without E sc is 58.1%. The average I rmsd values of five cases are less than that for case with all items. As can be seen, five energy items are all effective, among which structural neighborhood energy item is the most effective one.

Conclusion
In this paper, we calculate structural neighboring property on interface, through Voronoi diagram. We propose a novel statistical method to extract interacting residues, and interacting patches can be clustered as predicted interface residues. Experiments show that our method achieves better results than some state-of-the-art methods. Comparing to existing methods for binding sites    prediction, our approach improves overall F nat value by at least 3%. In addition, structural neighboring property can be adopted to construct an energy function, for evaluating docking solutions. It includes new statistical property as well as existing energy items. On Benchmark v4.0, our method has average I rmsd value of 3.31Å and overall F nat value of 63%. On CAPRI targets, our method has average I rmsd value of 3.46Å and overall F nat value of 45%.