1.1 Protein recognition domains (PRDs) and specificity determining residues (SDRs)

Peptide recognition domains are key elements on proteins with many roles in central signaling pathways [1]. PRDs are involved in many diseases and are a focus of drug discovery efforts [2]. Some viruses co-opt host PRDs via mimicry, emphasizing their relevance to human diseases [3–5]. Representative PRDs include PDZ, SH3 and kinase domains. These domains recognize short (usually around seven to ten amino acids long) motifs in their target proteins. Such motifs often lie in unstructured regions [6] and hence resemble peptides in their conformational flexibility. In particular, PDZ domains recognize hydrophobic C-terminal tails, SH3 domains recognize proline-rich motifs, and kinase domains bind several different classes of motifs; for instance, basophilic kinases bind arginine-rich motifs. As many PRDs are well behaved in solution, they are ideal model systems for studying protein-protein interactions and ligand binding specificity. New approaches, namely phage display [7], protein or peptide arrays [8–10] and oriented peptide libraries [11] have greatly facilitated the measurement of specificities and produced an hitherto unimaginable wealth of binding data. Despite these advances, the residues within the domains that determine their binding specificity have only partly been elucidated. Such specificity determining residues (SDRs) [12, 13] can be identified using dedicated experiments involving site-directed mutagenesis and subsequent measurement of specificity. For example, recent studies have directly tested effects of point mutations in kinases [14] and PDZ domains [7]. However, because of the size of the sequence space to be covered, exhaustive experimental search is infeasible. While co-crystal structures of PRDs with the bound ligand are often used to prioritize residues, identification of SDRs remains a complex problem: on the one hand, close proximity to the ligand does not necessarily implicate a residue in specificity determination. On the other hand, a residue that is far away from the ligand can also affect specificity due to secondary effects on the binding site [13, 15]. Hence, computational methods are needed to select and prioritize the positions to be tested through mutagenesis. Meanwhile, the aforementioned wealth of specificity data offers ample resources to be computationally analyzed for information about SDRs.

Currently, there are two major classes of computational approaches for SDR identification: first, there are structure driven approaches, making use of physical properties from protein structures, such as the hydration site free energies displaced by ligand atoms upon binding [16]. The second class of approaches adopts a statistical approach, and identifies SDRs by looking for sites that are conserved across or within functional groups [17], or more conserved in orthologs than paralogs [12]. Our method also takes on a statistical approach, the systematic nature of which enables the potential identification of non-local interactions between residues that are significant for peptide binding. Due to the lack of large-scale binding data previously, most current statistical methods attempt to detect SDRs based only on conservation signals from multiple sequence alignments of the PRDs. These predictions are noisy as the conservation of a residue could be due to roles other than binding specificity determination. Hence, these approaches do not make optimal use of currently available data. Conversely, utilizing the information about the actual bound peptide profiles from the recent large-scale binding datasets has the potential to boost the power of new predictions. In this paper we demonstrate how these new datasets can be incorporated in the search for SDRs. Previously we explored this idea using a proof-of-principle correlation-based method as part of a larger experimental study [14]. Here, we develop a new approach that is based on a novel formalism of mutual information.

1.2 Traditional and generalized covariation methods

Our approach is based on the notion that SDRs would covary with the binding specificity. The concept of covariation has first been used to study internal residue contacts in proteins [18, 19]. It has subsequently been used in the identification of energetically coupled residues [20], coevolution of proteins with their interaction partners [21–25], protein fold [26], functionally or structurally coupled residues within a protein [27, 28] and protein fusion [29]. The basic idea is to look for two sites in a multiple sequence alignment (MSA), or a pair of MSAs when studying the interaction between two proteins, that exhibit coordinated changes. Such covariation could suggest functional or structural dependency between the two sites, as there exists evolutionary pressure against their independent evolution. In the context of SDR identification, the covariation between a residue in the PRD and a residue in the bound peptide could imply its role in mediating the binding.

Many computational methods have been proposed for identifying covarying sites, including correlation scores [18, 30], Statistical Coupling Analysis (SCA) [20, 31], likelihood methods [32, 33], information theoretic methods [27, 34], independence tests [35, 36] and Bayesian approaches [37]. However, these methods cannot be applied directly to the identification of SDRs, because in a typical protein-peptide binding dataset [13, 38], each PRD does not bind to only one single peptide sequence. Instead, a PRD is associated with a binding profile, often represented in the form of a position weight matrix (PWM). In other words, instead of having two MSAs as in a typical covariation study, SDR identification adopts a more general setting with one MSA of PRDs on the one hand, and a list of respective PWMs on the other hand (Figure 1). Furthermore, a multitude of techniques have been developed in recent studies to handle various issues in the application of the covariation methods [34, 39–42], such as uneven representation of sequences in the MSA and overly diverse sites. These techniques also need to be extended in order to be applied in the current more generalized setting. In this paper we present a novel method to determine SDRs using an information theoretic approach. For each site from the MSA of the PRDs and each site from the aligned peptide PWMs, the method produces a numeric covariation score, which can be used as an indicator of the likelihood that the two sites are involved in binding specificity. Aggregating the scores could also suggest which peptides a PRD is likely to bind. Our approach uses only information from sequences, and thus can be applied to PRDs without known crystal structures.

2.1 Entropy and mutual information

Our method for identifying SDRs is based on an information theoretic measure called mutual information, which has been used in detecting covarying residues in the traditional setting [27, 34]. We first give some general definitions here, and then discuss how they are used in the identification of SDRs in the next subsection. Given a discrete random variable X with a distribution p(x) over the domain of X, , the entropy of X is defined as:

(1)

where log 0 is defined as 0, the asymptotic value of log(x) towards 0. In the general definition of entropy, the base of the logarithm can be set to any value. In analyzing protein sequences, where is the set of 20 amino acid residues, it is convenient to use 20 as the base so that the entropy value is always between zero and one [43].

Similarly, if we have two random variables X and Y with a joint distribution p(x, y), the joint entropy of them is defined as:

(2)

Based on the concepts of entropy and joint entropy, the mutual information between two random variables X and Y is defined as follows:

(3)

2.2 Adapting mutual information to the identification of SDRs

Suppose we are given an MSA A with n rows and m columns. Each row corresponds to a PRD sequence and each column corresponds to a site of the alignment. We use A _ij to denote the residue at site j of sequence i. The residues at a site can be viewed as a sample drawn from a distribution of residues specific to the site. Let A _j be a random variable for the residues at site j. Using Formula 1, we can calculate the entropy of A _j by replacing p(x) with the sample distribution, p _j (x), defined by the frequency of each residue at the site:

(4)

where is the indicator function with and .

The entropy can be interpreted as the uncertainty of which residue we would encounter at the site if we randomly draw a sequence from the MSA, where uncertainty here is mathematically quantified by the number of bits needed to encode the information on average. A completely conserved site has an entropy of zero, and indeed there is no uncertainty as the residue being drawn must always be the same. A site with equal probability of all 20 residues has the maximum possible entropy of one. In general, a more diverse site has a larger entropy.

Similarly, suppose we are given a set of n aligned PWMs W each with w sites. The ith PWM represents the peptides bound by the ith PRD in the MSA (Figure 1). Denote W _ik (y) as the probability of residue y at the kth site of the ith PWM. Let W _k be a random variable for the residues at site k. Again, we can calculate the entropy of W _k using Formula 1 by replacing p(x) with the expected probability of y in the different PWMs, p _k (y), assuming a uniform distribution of the observed sequences:

(5)

Now, we can compute the joint entropy of site j in the MSA and site k in the PWMs using Formula 2, based on the sample distribution of A _j and the probabilities W _ik (y):

(6)

The joint entropy measures the uncertainty of which two residues we would encounter at MSA site j and PWM site k if we randomly draw a sequence from the MSA and get its corresponding PWM.

Finally, we can compute the mutual information between MSA site j and PWM site k using Formula 3. Since H (X, Y) is the uncertainty of the two sites that persists even when we consider them together, subtracting it from H (X) + H (Y) gives the uncertainty that is eliminated by considering the two sites together. In other words, mutual information measures the information shared by the two sites. A larger mutual information indicates a stronger dependency between them. This could indicate a functional or structural relationship between the two sites. For example, it could suggest that the two sites coevolve in the sense that when the residue at the MSA site is changed, binding strength is restored by having a corresponding change at the PWM site.

2.3 Handling uneven sequence representation

In many cases, the input MSA for studying residue covariation has an uneven representation of sequences from different clades. For example, it is common to have more sequences from model organisms or species that are better studied, and fewer sequences from other species. As a result, the MSA could contain many sequences that are highly similar, and few that are significantly different. Each of the highly similar sequences contributes little additional information, but still has an equal amount of influence on the calculation of mutual information under the assumption of a uniform distribution of the observed sequences. Consequently, they could mask the novel information from low abundance sequences. To counteract this effect, we associate weights with the sequences so that the more unique ones receive higher weights. Statistically, it is equivalent to placing a prior distribution to the observed sequences if the weights are normalized to take values between zero and one and have a sum of one. Here we assume there is an external procedure for determining the weights. For instance, one possible way is to construct a phylogenetic tree from the sequences in the MSA, and then recursively distribute the total weight to different branches of the tree, so that each sequence in the crowded branches will receive a smaller share of weights [44]. Suppose the procedure assigns α _i as the weight of sequence i, we can redefine p _j (x), p _k (y) and p_j,k(x, y) as follows:

(7)

(8)

(9)

Entropy, joint-entropy and mutual information can then be calculated using these new definitions of probabilities.

2.4 Handling uneven site conservation

To deal with this problem, various kinds of normalization have been proposed to penalize overly diverse sites [34]. We will focus on the uncertainty coefficient [45], which was found to be one of the best normalized mutual information scores in our preliminary study. For an MSA site A _j and a PWM site W _k , the uncertainty coefficient is defined as follows:

(10)

We have also tried handling the problem using a statistical test. Specifically, we used mutual information as the test statistic to calculate how unlikely it is to get a mutual information at least as large as the observed one under the null hypothesis that the two sites are independent. The distribution of mutual information can be obtained by permuting the residues of a site, or by using a chi-square approximation. It was shown that when n is large, (2 ln 2n)MI(X, Y) tends to have a chi-square distribution with degrees of freedom [46, 47]. It turns out that the results based on this statistical test were not better than using the simple normalization approach. We thus focus on the use of the uncertainty coefficient measure in handling uneven site conservation in the remaining of this paper.

Table 1 also demonstrates the tradeoff between the mutual information between two sites and their individual conservation, both of which are indicators of their functional importance. One may try to derive a measure that takes both into account, similar to what the Sequence Harmony method handles both the conservation and similarity of two groups of residues simultaneously, for the problem of identifying important residues that determine the functional differences of protein subfamilies [48]. We leave the derivation of a new covariation measure to a future study.

2.5 Predicting the PWMs of bound peptides

One important use of the covariation scores is to contribute towards predicting the PWMs of bound peptides. This can be done by aggregating the detected covariation signals. Suppose we are given a new PRD sequence (the (n + 1)th sequence) of the MSA M without the corresponding PWM of its bound peptides. We would like to predict the PWM based on the n + 1 sequences in the MSA and the n known PWMs. We investigate the use of the covariation scores in this problem by comparing a prediction method that considers site covariation with two methods that do not.

A simple prediction method that does not take site covariation into account is to perform a weighted averaged of the known PWMs, where the weights are based on the similarity of the new PRD sequence and the original ones. Specifically, the probability of finding residue y at site k of the bound peptides of the new PRD is predicted by the following formula:

(11)

where s(i, i') is a similarity between sequences i and i' in the MSA, such as their sequence identity:

(12)

Using the covariation scores, we propose an alternative way to define the similarity function. Each MSA site receives a different weight in the calculation, where the weight depends on the covariation score between the site and the target PWM site k. In other words, the similarity score s(i, i') is replaced with a new score s _k (i, i') that is specific to k:

(13)

(14)

where the uncertainty coefficient UC(A _j , W _k ) is computed based on the n original sequences.

We also investigated if prediction accuracy can be improved by using a more complex model. Specifically, we treated each MSA site as a categorical attribute and trained a regression tree model for each probability value in the PWM. The models were then applied to the new sequence to predict the PWM of its bound peptides. We implemented this method using REPTree of the Weka package [49].

To evaluate the effectiveness of the different methods, we performed left-out validation as follows. Each time, we drew a random sample of PRDs to form the testing set. Each sequence in the testing set took turn to take the role of the (n + 1)th sequence. The sequences not included in the sample formed the training set. These sequences were used to compute covariation scores and train prediction models. The procedure was repeated 1,000 times for PDZ and SH3 and 50 times for kinase (due to the long running time) with different random training-testing splits, and the average performance of the trained models on the testing sets was recorded. To eliminate near-identical sequences in the training and testing sets, for sequences with 90% or more identity, we kept only one of them in the dataset before making the training-testing splits. As most of the synthetic PDZ sequences (described below) are highly similar, we excluded this dataset from this part of study.

Each predicted PWM was compared to the actual PWM, and a prediction error was computed as the root-mean-square difference between their distributions per site:

(15)

where and are the predicted and actual PWMs for the bound peptides of the testing sequence, respectively, and the inner summation is taken over all 20 amino acid residues.

3.1 Application of the method to PDZ, SH3 and kinase domains

3.2 Covariation score correlates with physical proximity and reconfirms previous findings

The covariation score between two sites is found to correlate negatively with the physical distance between them, regardless of the exact definition of the distance measure [see Additional file 1 Figure S1, Additional file 2 Figure S2 and Additional file 3 Figure S3 for the results] when the distance between residue centers, alpha carbon atoms and closest atoms minus their van der Waals' radii were used, respectively. All P-values were computed using Fisher transformation [60]. For PDZ, we have also compared the correlations based on several other PDZ structures, and observed similar patterns [Additional file 4 Figure S4].

Since low-scoring pairs are more subject to noise, here for each PRD site, we focus on the peptide site that gives the highest uncertainty coefficient with it (Figure 2). The site pairs with the highest covariation scores are listed in Table 2.

For natural PDZ domains, the highest-scoring pair (circled) is between Leu60 (β3:α1-1, structural nomenclature from [53]) of the PDZ domain in 2H2C and position -1 of the binding motif on the bound peptide. These residues are in physical contact in the dimer structure (Figure 3, all visualization produced using VMD [61]). Interestingly, it has been reported that the side chain at β3:α1-1 can contribute to the recognition of the -1 position of the motif [53], and in the crystal structure from Shank1, a salt bridge is observed between Asp(β3:α1-1) of the PDZ domain and Arg(-1) of the ligand [62]. Our covariation analysis has thus identified these verified SDRs of the PDZ domains in silico.

For SH3 domains, the highest-scoring pair is between Asn71 of the SH3 domain and Leu7 (P+2 residue) of the ligand in 1N5Z. These residues are in close physical proximity in the crystal structure (Figure 4). Interestingly, we found that the MSA residues in the first (Asn71), third (Ile73) and fourth (Tyr72) top-scoring pairs are consecutive in the protein sequence, and two of them are close to Leu7 of the ligand. In a previous study [63], the corresponding residue of Tyr72 on P53BP2, which is an unusual leucine, was hypothesized to cause the protein to bind a peptide very different from its usual ligands. However, mutating the leucine to tyrosine did not affect recognition specificity. As the corresponding residue of Ile73 on P53BP2 is also mutated from the class consensus, and the covariation scores for Asn71 and Ile73 are also high, the three residues may have some combined effects in affecting recognition specificity.

It is also known that the residue corresponding to Glu31 in the RT loop of SH3 domains is a major determinant of the identity of the P-3 residue of the ligand [64, 65]. Since the P-3 residue does not exist in the 1N5Z structure, it is not included in the correlation plots. However, when we checked the covariation scores, indeed the ligand residue having the highest score with Glu31 is the P-3 residue. This observation illustrates the potential of our method to identify SDRs when structural information is not available.

For protein kinases, the top-scoring pair between Tyr229 of the kinase domain in 1ATP and position +1 of the binding motif in the bound peptide is not physically close. However, the next two pairs (between Leu205 and position +1, and between Leu198 and position +1) are both close in proximity (Figure 5). Both Leu198 and Leu205 were previously found to have hydrophobic interactions with position +1 of the peptide [66, 67]. They are two of the three residues (the third being Pro202, which has the highest covariation score with position +1 as compared to other peptide positions) that form the binding pocket for position +1 of the peptide, and contributes to the positioning of a significant portion (-3 to +1 positions) of the bound peptide [67].

We found that the pair between Tyr229 and position +1, besides being the top-scoring pair when uncertainty coefficient was used as the covariation score, was also consistently one of the highest-scoring pairs when covariation was measured by other normalized forms of mutual information. These consistent results made us hypothesize that Tyr229 is also involved in determining binding specificity of the kinase domain, and has long-range coupling with position +1 of the bound peptide. In a previous study, this residue was predicted to be involved in linking nucleotide binding and peptide binding in protein kinases [68]. Interestingly, in the study Tyr229 and Leu205 are predicted to belong to the same special network (termed the θ-shaped network) of related residues. It thus might be the case that Tyr229 acts through the residues in the network to affect the recognition of the +1 position of the peptide.

3.3 PDZ predictions are consistent with mutagenesis results

We further validated our predictions with the natural PDZ domains by using a mutagenesis dataset from a domain specificity study [13] as described in the Materials and Methods section.

In our covariation calculation, among the ten mutated sites, four were filtered as they were too conserved (Ile36 [β2-1], Ile38 [β2-3], His88 [α2-1] and Val92 [α2-5]). Interestingly, in the mutagenesis study, indeed no significant changes to the peptide PWM could be observed for Ile36 and Ile38. The high conservation of them is thus probably caused by a structural or functional role that is independent of peptide binding.

The predicted pair between Phe34 and position -3 has a distance of 15.0 Å in the PDB structure 2H2C. If the SDRs of different class I PDZ domains are similar, this predicted pair is another example that suggests our covariation method could potentially identify physically distant SDR pairs.

3.4 Covariation scores improve prediction of bound peptide profiles

As described in the Materials and Methods section, we compared three different methods for predicting the profiles of bound peptides. The prediction results are shown in Figure 6.

In general, the prediction error is lower with a smaller fraction of PRDs left out for testing. This is expected, as having more PRDs in the training set allows for more accurate computation of covariation scores and the construction of more informative prediction models.

In all cases, the weighted average method using covariation scores outperformed regression tree and the weighted average method using uniform scoring. This suggests that the covariation score provided a meaningful way to weight useful features (that is, PRD sites) for predicting residue distributions of the bound peptides. Interestingly, while the regression tree method also performed feature weighting, in general it performed worse than both weighted average methods. The low performance could be due to over-fitting, as the regression trees are rich in expressive power while the numbers of PRDs in the datasets are small.