DBAC: A simple prediction method for protein binding hot spots based on burial levels and deeply buried atomic contacts
© Li et al; licensee BioMed Central Ltd. 2011
Published: 20 June 2011
Skip to main content
© Li et al; licensee BioMed Central Ltd. 2011
Published: 20 June 2011
A protein binding hot spot is a cluster of residues in the interface that are energetically important for the binding of the protein with its interaction partner. Identifying protein binding hot spots can give useful information to protein engineering and drug design, and can also deepen our understanding of protein-protein interaction. These residues are usually buried inside the interface with very low solvent accessible surface area (SASA). Thus SASA is widely used as an outstanding feature in hot spot prediction by many computational methods. However, SASA is not capable of distinguishing slightly buried residues, of which most are non hot spots, and deeply buried ones that are usually inside a hot spot.
We propose a new descriptor called “burial level” for characterizing residues, atoms and atomic contacts. Specifically, burial level captures the depth the residues are buried. We identify different kinds of deeply buried atomic contacts (DBAC) at different burial levels that are directly broken in alanine substitution. We use their numbers as input for SVM to classify between hot spot or non hot spot residues. We achieve F measure of 0.6237 under the leave-one-out cross-validation on a data set containing 258 mutations. This performance is better than other computational methods.
Our results show that hot spot residues tend to be deeply buried in the interface, not just having a low SASA value. This indicates that a high burial level is not only a necessary but also a more sufficient condition than a low SASA for a residue to be a hot spot residue. We find that those deeply buried atoms become increasingly more important when their burial levels rise up. This work also confirms the contribution of deeply buried interfacial atomic contacts to the energy of protein binding hot spot.
Protein-protein interactions are dominated by hydrogen bonds, salt bridges and hydrophobic contacts across the interface [1, 2]. These local interactions have to be desolvated, densely packed, and hence deeply buried to make contribution to the binding free energy [3–6]. This is why the energetically important hot spot residues in the interface tend to cluster into local regions with low solvent accessible surface area (SASA) values [7, 8].
Identifying these energetically important residues, which can offer useful information to protein engineering and better understanding of protein-protein interaction , is usually done by site-directed alanine mutagenesis. This experimental method mutates the target residue into alanine which only has a C β heavy atom on its side-chain [10, 11]. A residue whose mutation results in a large binding free energy change (≥2.0 kcal/mol, for example) is defined as a hot spot residue .
Many feature-based [13–17] energy-based [18–23] and even feature-and-energy-combined [24, 25] computational approaches have been proposed to address the hot spot prediction problem. Almost all of these feature-based methods use SASA information of the residue as a critical feature in the prediction. A low SASA is necessary for a residue to be a hot spot residue; however, it is not sufficient, as a large number of non hot spot residues also have low SASA values. Therefore, SASA is not effective for distinguishing between slightly buried residues—a large part of which are non hot spot residues—and deeply buried residues that are very likely to be hot spot residues.
In this work, we introduce a new descriptor for protein atoms and residues. It is named “burial level”. In the definition of burial level, the buried immobilized water molecules are treated as an integral part of the protein complex. We show that our definition of residue burial level is nicely correlated to ∆∆G. A high burial level is not only in general necessary for hot spot residues but also more sufficient for them in comparison to SASA. In other words, most hot spot residues tend to have high burial level while most non hot spot residues are exposed or just slightly buried. We also define the burial level of atomic contacts and we determine the number of three types of buried interfacial atomic contacts at different burial level that are directly broken when the residue is substituted by alanine. The number of those deeply buried atomic contacts together with the burial level of the residue itself are further fed into SVM as features to classify interfacial residues into hot spot residues or non hot spot residues. We name this SVM-based model DBAC since the features used are mainly based on the Deeply Buried Atomic Contacts. By applying our method to a data set of 258 mutations, we achieve an F measure of 0.6237 under the hot spot definition of ∆∆G ≥ 2.0 kcal/mol, which is better than other computational methods. We also conduct a detailed analysis of the features used in this work; and we find that hot spot residues tend to have significantly more deeply buried atomic contacts than non hot spot residues.
Our data set is collected by retrieving the experimental alanine mutagenesis data from the alanine scanning energetics database (ASEdb)  and other previously published works [27–31]. We require that: the 3D structure of the wild-type protein complex is solved by X-ray crystallography and is reported in PDB , and the associated solvent information is also included in the PDB file. We do not consider protein-ligand interaction or protein-peptide interaction in this work; thus those interactions without an extended interface are excluded. The reason is that the interfaces of protein-ligand interactions are small and most interfacial residues are exposed in the solvent to a certain degree; thus the burial levels of the atoms are mostly very low and imply little information. The structural similarity of the complexes are tested by the CE algorithm . If the two chains of the two complexes have a significant similarity, their binding interfaces are further examined to ensure that there is no redundancy in the data set. Furthermore, only mutations in the interface are considered. We used another version of this data set in our previous work , where the requirement that the mutations have to be in the interface was not applied.
Our data set in this work consists of 258 mutations distributed in 13 protein complexes. Hot spot residues are usually defined by setting ∆∆G ≥ 1.0 kcal/mol or ∆∆G ≥ 2.0 kcal/mol. We prefer the second choice, as only a higher ∆∆G threshold can reflect the direct influence of the mutation. That is, the interfacial atomic contacts that are directly broken by the mutation are taken into consideration with more weights. Under the ∆∆G ≥ 2.0 threshold, there are 50 hot spot residues and 208 non hot spot residues in our data set. Some researchers even suggested that a residue should have a ∆∆G higher than 4.0 kcal/mol so as to have a strong impact on the binding of the two proteins . In practice, a lower value has to be considered in order to get enough data for statistical analysis .
The data set is available at http://22.214.171.124:8080/DBACdata .
Our definition of burial level is based on atomic contact graph. The atomic contact graph of a protein complex is an undirected graph with heavy atoms as nodes and atomic contacts as edges. The atoms in this graph are labeled as exposed or buried according to its SASA. If the SASA of an atom is not less than 10.0Å 2, it is exposed, otherwise it is buried. SASA is calculated by the NACCESS software based on the Lee-Richards algorithm . All the exposed water molecules, which we consider as a part of the bulk solvent, are removed, while the buried water molecules are kept as a part of the complex. Thus the oxygen atoms of the buried water molecules are a part of the atomic contact graph.
The atomic contact is defined by a distance threshold and the Voronoi diagram. Voronoi diagram decomposes the 3D space into cells by assigning every point in the space to its nearest neighboring input site. Here in this work, the input sites are the positions of the atoms in the complex structure. If two atoms’ Voronoi cells are adjacent to each other, they are somehow “sheltering” each other. We define the atomic contact by adding an extra distance requirement to Voronoi diagram. Two atoms are considered to be in contact if they have a distance less than their Van der Waals radius plus the diameter of a water molecule (2.75 Å) and they share a Voronoi facet. We actually used the Delaunay diagram that is dual to Voronoi diagram. Two facet-sharing Voronoi cells correspond to two connected points in the Delaunay diagram. The Delaunay diagram is calculated by using the qdelaunay program in qhull . This distance threshold, 2.75 Å, is based on a water-free idea and it has been used in .
In an atomic contact graph, the burial level of an atom is defined as the length of the shortest path from this atom to its nearest exposed atom. For example, the burial level of exposed atoms is 0 and the burial level of their immediate buried neighbors is 1. We calculate the burial levels by adding a pseudo node, which represents bulk solvent, to the atomic contact graph. This node is connected to all of the exposed nodes directly. Then the burial level of any atom equals to the length of the shortest path from this atom to the pseudo node minus 1. This is exactly the single-source-shortest-path problem and it can be solved using Dijkstra’s algorithm .
The reason for using a Voronoi-diagram-combined definition of atomic contact is as follows. If only distance information is used, there will be many false atomic contacts in the atomic contact graph whose two atoms cannot contact with each other at all (as they do not share Voronoi facet thus they are spatially blocked by other atoms), and the atomic contact graph will be a trivial discretization of the Euclidean distance between atoms, and the atom burial level will only depend on the distance of the atom to the surface of the complex, especially when a large distance threshold (2.75 Å) is used. Adding Voronoi diagram to the definition makes the burial level depend also on the organization of the atoms inside. Intuitively, the burial level of atoms in a protein complex depends on the size of the protein complex. In general, the larger the protein complex is, the more deeply buried atoms there are. Burial level also depends on the shape of the interacting proteins. For example, globular proteins and protein complexes generally have more deeply buried atomic contacts than those with other shapes.
Note that the calculation of burial level requires information on buried water molecules. In our previous work , we have systematically analyzed the contribution of water molecules to the calculation of burial level as well as to protein binding hot spots.
The burial level of a residue is the average value of the burial levels of all atoms in the residue. For an atomic contact, if the burial levels of the two atoms are the same, the burial level of the atomic contact is taken as the burial level of the two atoms, otherwise it is defined as the smaller one of the two burial levels. The difference of the burial levels of two contacting atoms is at most one.
There are some existing concepts that are related to burial level. In [39–43], the authors defined their concept of depth of an atom as the Euclidean distance to the closest exposed atom or to the closest surface water molecule. There are even some sequence-based methods [44–46] that are capable of predicting its value. This definition is only based on the Euclidean distance and hence it cannot capture the contacts between atoms or the organization of the atoms. And the calculation of the shortest distance contains an exhaustive search among all exposed atoms or surface water molecules for every buried atom, while our graph-based concept of burial level can be easily calculated by transferring the calculation into a single-source-shortest-path problem.
When a residue is mutated into alanine, some interfacial atomic contact are directly broken because of the removal of certain atoms; Some other interfacial atomic contacts may also be broken or distorted due to the conformational change in the local region . Given a residue, we consider its directly broken interfacial atomic contacts, namely those contacts formed by an atom other than C, N, O, C α or C (these atoms are the remainder after a residue is mutated into alanine) and any contact partner from the other chain, even its backbone atoms. For example, for a serine in a complex, the atomic contact between its Oγ and a C α of any residue from the interacting partner is a directly broken atomic contact of this serine, while the contact between its C β and whatever atom from the interacting partner is not. We classify the atomic contacts into three types. If a contact is between a positively charged atom and a negatively charged atom, which usually corresponds to a salt bridge, it is called a Type-I contact. If a contact is between a hydrogen bond donor and a hydrogen bond acceptor, which usually is a hydrogen bond, it is classified as Type-II contact. Contacts that are neither Type-I nor Type-II are classified as Type-III. Here, the definitions for positively charged atoms, negatively charged atoms, hydrogen bond donors and hydrogen bond acceptors are as given in . We do not further divide the Type-III contacts into subtypes such as other polar contact, hydrophobic contact and so on because they are all not as specific as Type-I and -II contacts. Note that the definitions for Type-I and Type-II contacts are not exactly the same as salt bridges and hydrogen bonds in terms of geometrical requirements, yet they can be still very important .
In this work, we use deeply buried atomic contacts whose burial level is not less than 2. We refer to atomic contacts at burial level 0 as exposed atomic contacts and those at burial level 1 as slightly buried atomic contacts. Let C(i, j) denote the number of Type-i directly-broken interfacial atomic contacts at burial level j of a residue. Then our model uses 6 features to describe a residue: C(I, ≥ 2), C(II, 2), C(II, ≥ 3), C(III, 2), C(III, ≥ 3), plus the burial level of the residue. An SVM model based on this feature set is named DBAC (Deeply Buried Atomic Contacts). For comparison, we have also built another model named AC (Atomic Contacts) based on another feature set comprising C(I, 0), C(I, 1), C(I, ≥ 2), C(II, 0), C(II, 1), C(II, 2), C(II, ≥ 3), C(III, 0), C(III, 1), C(III, 2), C(III, ≥ 3), and the burial level of the residue. The maximum value of burial level depends on the size of the protein complexes, the size of the interfaces as well as the shape of the complex and the organization of the atoms. In general, very few contacts have burial level larger than 3, so we do not distinguish further burial levels larger than 3. For Type-I contact, there are very few cases that have burial level larger than 2, thus we do not use C(I, > 3) as a feature but merge it with C(I, 2) into C(I, ≥ 2).
Support Vector Machines (SVMs) are widely used in many classification and regression problems. They have also been adopted in hot spots prediction problems [15, 17, 24] with various feature sets and training-testing protocols. In this work, we use the LIBSVM software , which is a tool for SVM model training and testing available at http://www.csie.ntu.edu.tw/~cjlin/libsvm. We use the radial basis function (RBF) as the kernel. We do not conduct feature selection because our method is straightforward, and the number of features is not large. However, we evaluate performance on two different feature sets: the deeply buried atomic contacts only (by DBAC) and all the atomic contacts (by AC). The latter feature set is evaluated just for comparison.
The performance is evaluated under leave-one-out cross-validation. To avoid over fitting, we have strictly followed a nested-loop cross-validation procedure. There are 258 mutations in our data set, each time one mutation is taken as the test data and the remaining 257 mutations are used to train the model. The two parameters, namely cost and gamma, are optimized on the training data by a grid search. The grid search evaluates the performance, F measure, of SVMs with different parameter values on the training data using 5-fold cross-validation, and the parameter values with the best performance are chosen to build a training model on the training data. This training model is then applied to the test data, that is, the mutation held out in advance. This procedure is repeated 258 times until every mutation in the data set is tested.
Performance is measured by sensitivity, precision, specificity, accuracy and F measure (F1). These measures are defined asfollows: , , , , and , where TP, FP, TN and FN are the number of true positives, false positives, true negatives and false negatives, respectively. A better classifier should predict hot spot residues with less false positives and less false negatives; thus the F measure, which combines sensitivity and specificity, is used to indicate the overall performance.
We also test the significance of the difference in ∆∆G values of predicted hot spot and non hot spot residues. A classifier divides the mutations in the data set into two groups: computational hot spot residues and computational non hot spot residues. The significance of the ∆∆G value difference in these two groups are tested by Mann-Whitney test . The result of a classifier with higher F1 value can be less significant when its false positives have very low ∆∆G values (near 0 kcal/mol or even negative) and the false negatives have high ∆∆G values.
We also examine the value distribution of individual features in hot spot and non hot spot residues. The significance of the difference in the two classes is also tested by Mann-Whitney test.
Performance of our method (DBAC) in comparison with using all atomic contact (AC) and Robetta
The reason we use leave-one-out cross-validation is that we have a small data set and, moreover, there are only a small number of positive samples (hot spot residues). To test the robustness of our method, we evaluate performance using leave-n-out cross-validation under the same training-testing procedure. We find that when n is not large (< 7), the performance change is not significant, sometimes better (F1=0.6304, n=5) and sometimes worse (F1=0.5934. n=3) than that by using leave-one-out. Anyway, as shown later, no matter how large n is, our performance is always better than Robetta and FoldX. For example, the performance of our method by 5-fold (leave-51-out) cross-validation is 0.5870 in F measure. This indicates that our method is robust.
We compare our method with three energy-based methods, Robetta [19, 21], FoldX [20, 50] and EGAD . Robetta is an online service. It can be used to predict the ∆∆G value of interfacial residues by computational alanine scanning based on an energetic function. It can thus be applied to hot spot prediction. Actually, it is a widely recognized gold standard for benchmark comparison in the field. Its performance on our data set is shown in the forth row of Table 1. Our performance is remarkably better than that of Robetta in terms of both F1 and p-value.
FoldX is also available online. It is able to predict the change in both protein stability and affinity. Its energy function contains different sources of contributions such as van der Waals interactions, hydrogen bonds and even water-bridged interactions. We calculated the ∆∆G of mutations in our data set by using FoldX version 3.0 beta 4 with default parameters. From the fifth row of Table 1, although FoldX shows a better performance than Robetta, probably due to the fact that it is being updated, our method still achieves a better performance than that by FoldX.
Comparison of our method with EGAD
These energy-based method are complicated and time-consuming. The energy functions usually contain many terms that represent different kinds of energies. Both binding and folding of proteins can affect the binding free energy between two proteins. But binding and folding are very complicated processes whose details are difficult to capture. When a residue is mutated into alanine, the new structure of the mutated protein and the mutated protein complex must be predicted to get the values of all energy terms of the mutated structure, which is also very difficult. Thus the ∆∆G are hard to be accurately estimated even by these complicated energy functions. From Tables 1 and 2, the performance of these energy-based methods are not very good yet.
We also compare our method with another machine learning method, MINERVA , which uses SVM as well and is based on a larger feature set containing various aspects of information of target residue such as weighted atomic packing density, relative surface area burial, weighted hydrophobicity and so on.
Comparison of our method with MINERVA
Statistical analysis on the features.
5.0897 × 10 -10
C(I, ≥ 2)
1.2419 × 10 -9
C(II, ≥ 3)
3.5031 × 10 -6
1.9061 × 10 -16
C(III, ≥ 3)
1.1621 × 10 -9
We show in Table 5 the performance of our method after one feature is removed under the same training-testing protocol described in Methods. It seems that the removal of residue burial level, C(I, ≥ 2), C(II, 2) or C(II, ≥ 3) has little impact on performance. A reason is that the features we are using are somehow correlated with each other. In general, when an interfacial residue is deeply buried with high residue burial level, it has several deeply buried atomic contacts with the other side. On the other hand, the removal of C(III, 2) or C(III, ≥ 3) reduces performance a lot. The reason is that there are often many Type-III atomic contacts in hot spot residues, as Type-III atomic contacts is not as specific as Types I or II. This suggests that we can further divide Type-III atomic contacts into subtypes. We have actually tried dividing Type-III contacts into other polar-polar contacts and hydrophobic contacts; and it turns out that the performance change is not much and hydrophobic contacts become the new dominant one. However, this does not mean residue burial level, C(I, ≥ 2), C(II, 2) and C(II, ≥ 3) do not contribute to the performance because all the six features have significantly different values in hot spot and non hot spot residues. In fact, we can achieve an F measure of 0.5 with only the residue burial level.
Performances of our method after one feature is removed.
C(I, ≥ 2)
C(II, ≥ 3)
C(III, ≥ 3)
In contrast, as shown in Figure 2(b), hot spot residues tend to have a high burial level, while non hot spot residues do not. More than 60% of hot spot residues have a burial level no less than 2.0, whereas less than 20% of non hot spot residues have such burial levels. Thus, we conjecture that a high burial level is not only necessary but also more sufficient than a low SASA value for a hot spot residue.
Type-I atomic contacts roughly correspond to salt bridges. Some researchers believe that buried salt bridges provide neutral or even negative contribution to protein stability [52, 53] because the desolvation of charged groups requires more energy than the interaction energy of the formation of the salt bridge . But in protein-protein interaction, it is found that interfacial salt bridges are more buried than intra-chain salt bridges, and the salt bridges are found favorable across the interface . Perhaps the two proteins are folded independently with more charged residues exposed and their conformation change during complex formation is very restricted, thus the two proteins prefer to interact in an electrostatic complementary manner.
Hydrogen bonds play a key role in protein-protein interaction . Most interfacial hydrogen bonds are extremely buried; and the more buried a hydrogen bond donor/acceptor is, the more likely it is to form a hydrogen bond . Thus, being buried is favorable for interfacial hydrogen bonds. Figure 3(b) shows the percentages of hot spot residues and non hot spot residues whose C(II, x) are larger than 0. It can be seen that nearly 30% of non hot spot residues have exposed Type-II atomic contacts, but very few of them have deeply buried hydrogen bonds. The case is totally different in hot spot residues. There are more hot spot residues that have deeply buried Type-II atomic contact while a few of them have exposed ones. The number of residues that have extremely buried (burial level≥ 3) atomic contacts is limited by the size of the protein complexes.
Type-III contacts contain all other kinds of contacts that are neither salt bridges nor hydrogen bonds, including hydrophobic contacts and other polar contacts. Actually hydrophobic contacts are not specific contact between atoms but are the packing of groups of hydrophobic side chains. The contribution of hydrophobic contacts to bonding free energy is correlated with the buried surface area . Thus energetically important hydrophobic contacts are those buried ones. Generally, protein-protein interfaces are dominated by salt bridges, hydrogen bonds and hydrophobic contacts; but sometimes other contacts also make contribution to the binding . A hot spot is usually a densely packed region in the interface, thus the number of buried contacts of a hot spot residue tends to be large, which can be reflected by deeply buried Type-III contacts. As shown in Figure 3(c), more than 80% of hot spot residues have Type-III contact at burial level 2 and only about 20% non hot spot residues have Type-III contact at this burial level.
Another example as shown in Figure 4(b) is TYR-29 of barstar. This residue is exposed with an SASA of 64.12 Å2; however, it is a hot spot residue with a ∆∆G of 3.4 kcal/mol. There are only two hot spot residues that have SASA larger than 60Å2 in our data set. We can still successfully identify it as a hot spot residue by using its deeply buried atomic contacts. The side-chain of tyrosine which contains an aromatic ring and a hydroxyl group is capable of forming aromatic π-interactions and hydrogen bonds . As can be seen from the figure, although TYR-29 of barstar is partially exposed, its side-chain stretches into the complex and forms many deeply buried atomic contacts. For the 4 Type-II interfacial atomic contacts shown in the figure, 3 are deeply buried and 1 is slightly buried. There are another 8 deeply buried Type-III contacts, 7 of which are made by the aromatic ring and 5 are atomic contacts with HIS-102, an active site residue of barnase .
We have proposed a feature-based method to predict protein-binding hot spots by using deeply buried interfacial atomic contacts that are directly broken during alanine substitution. The method is based on a graph theoretical definition of burial level of residues, atoms and atomic contacts. We achieved an F measure of 0.6237 when ∆∆G ≥ 2.0 is used as the threshold to define hot spot residues. The burial level of a residue is more intuitive than the concept of SASA. It is nicely correlated with the ∆∆G of a residue. We have shown that a high residue burial level is in general necessary for a residue to be a hot spot residue. Furthermore, it is more sufficient than SASA, a frequently used feature in existing hot spot prediction methods. Our results also reveal that hot spot residues tend to have deeply buried atomic contacts while non hot spots tend to have exposed and slightly buried ones. This is consistent with previous studies that emphasize the energetic contribution of buried salt bridges, hydrogen bonds and hydrophobic contacts.
This research work was funded by two Singapore MOE Tier-2 grants (T208B2203 and MOE2009-T2-2-004).
This article has been published as part of BMC Systems Biology Volume 5 Supplement 1, 2011: Selected articles from the 4th International Conference on Computational Systems Biology (ISB 2010). The full contents of the supplement are available online at http://www.biomedcentral.com/1752-0509/5?issue=S1.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.