Mutual information and variants for protein domain-domain contact prediction

3D and Reduced Alphabet MI adjustments

To investigate methods that might further enhance the predictive power of MI variants we designed two adjustments. The first adjustment considers triangles of columns rather than pairs, based on the idea that interactions occur in patches[36]. This variant is denoted by the suffix 3D. The second adjustment, suffixed RA, reduces the 20 amino acids to seven categories based on their physical and chemical properties, with the aim of reducing noise.

The idea behind the 3D version is that protein binding involves patches of residues in contact. This idea has been previously used to predict contact residues[12, 36]. Furthermore, the success of MIc lies in its normalising factor, the coevolutionary pattern similarity (CPS) score, which estimates the coevolutionary relationship between the pair of residues currently under consideration and all other residues in the MSA[25]. We thus speculated that adding additional residue information to MI and MIp pair scores may enhance their domain-domain contact predictive capabilities. Hence we created new versions of MI and MIp that consider triangles rather than pairs of columns to identify contacts (MI3D (Equation 13) and MIp3D (Equation 14)). Increasing the dimensionality of MI and MIp in this manner surprisingly worsened performance in both cases; the precision at 20% recall of MI3D and MIp3D are 19.9% and 31.8% respectively as compared to precision of MI and MIp of 24.4% and 42.3% (Figure2 A and B, and Table1). We conjecture that adding an extra dimension to MI and MIp magnifies the noise in the MSA more than it boosts the signal.

Assuming that contact residues mutate in a correlated manner in order to maintain their interaction, it is not evident how much of a change a residue can undergo while still maintaining its contacts. Using a reduced alphabet residue set addresses this point; as it groups residues by their physiochemical properties, under the assumption that residues with the same physiochemical properties will maintain similar interactions. Grouping the 20 amino acids into seven categories only improved the performance of basic MI and MI3D, which rose in precision at 20% recall from 24.4% to 28.4% and 19.9% to 23.5% respectively (Figure2A and Table1). In all other cases the reduced alphabet (RA) appeared to reduce noise as well as signal (Figure2, Table1 and Additional file2: Figure S2).

Case study

The case study Skerker et al.[27] has received a lot of attention for successfully determining inter-protein contact specificity residues with the aid of MI. The authors used original MI (Equation (4)) to determine a subset of contact residues that allow for specific binding of a histidine kinase (HK) with its interacting response regulator (RR). The MSA provided by these authors does not contain the sequence of the structure used in their analysis. Hence we ran MI and MIc on the HK-RR MSA provided by Hamer et al.[12], which does include the sequence of this reference structure.

As Skerker et al. were interested in residue pairs only between the DHp domain (four helix bundle) of the HK and its interacting RR, only these MI and MIc scores were considered when examining performance. In accordance with our evaluation method on the 40 test cases, all buried residues were eliminated, as were residues corresponding to columns that had one or more gaps or an entropy of 0. This leaves us with 46 DHp residues, nine of which are contacts, and 68 RR residues, amongst which 24 are contacts.

As we have no score cut-off for predictions we check the number of correct predictions among the top nine predictions for DHp. If there was no relationship between the MI scores and contact sites, then the number of correct predictions would follow a Binomial distribution with sample size nine and probability of success 9/46. Under this model we would expect 1.76 correct predictions.

For RR there are 24 contact residues. We check the number of correct predictions among the top 24 predictions. If there was no relationship between the MI scores and contact sites, then the number of correct predictions would follow a Binomial distribution with sample size 24 and probability of success 24/68. Under this model we would expect 8.47 correct predictions.

Datasets

For this domain-domain MI investigation we use proteins that have two domains, rather than protein complexes, and treat each domain as a separate protein. In this manner we can be sure that accurate “protein-protein” pairings are used in the MSA. The MSAs are taken from[12], available athttp://www.stats.ox.ac.uk/research/bioinfo/resources, which in turn are based on datasets in[18, 37–39]. The proteins for which each MSA was constructed has a known pdb structure[40] of X-ray resolution 2.5Å or better, and well annotated domain boundaries[12]. This structure is henceforth referred to as the “reference structure” and is used to identify surface, buried, contact and non-contact residue columns within the MSA. The MSA was generated by using the structural protein as a BLAST query[41, 42] against the NCBI-NR database[43]. The homologs identified were made non-redundant at the 90% level using Cd-hit[44]. The final alignment was generated using MUSCLE[45] and MaxAlign[46].

Within the 40 inter-domain MSAs there are non-standard amino acid entries, such as B, Z, X, * and ?. As there is no established method of processing these sequencing uncertainties, we choose to treat them as a gap, while the Brown and Brown pipeline code processes them as additional amino acids[11].

Identifying the surface versus buried residue pairs

For each reference structure protein in the dataset, we calculate the solvent accessibility of the residues using JOY[48]; each domain is treated as a separate entity. In the reference structure, residues that are >7% accessible to a 1.4Å radius water molecule are denoted as “surface” residues[48]. Those that do not meet this criterion are termed “buried.” This information about a residue is then annotated to the entire MSA column to which it belongs.

Employing this criterion on our 40 test cases, along with eliminating residue columns that have an entropy of 0 or contain a gap, leaves us with 5483 surface residues and 2364 buried residues. These numbers decline to 5362 and 2174 respectively when employing the 40 reduced alphabet MSAs, as the reduced alphabet MSAs have a greater number of columns with 0 entropy (Equation (1)). The ratios of surface to buried residues in the reduced and non-reduced alphabet sets are 2.5 and 2.3 respectively. The number of 0 entropy columns in the reduced and non-reduced alphabet set are 668 and 326 respectively, while the number of columns with one or more gaps in both sets are 3479.

Identifying the contact versus non-contact residues

Residues within the binding interface of a pair of interacting domains are labelled as “contact” residues (Figure3). There is not one particular accepted definition for contact residues[26, 49–51]. Here we classified a residue in the representative protein structure as a “contact” residue if:

If not all of the above criteria are met, a residue is denoted as a “non-contact” residue.

Using these criteria, and once again ignoring residue columns having an entropy of 0 or a gap, leaves us with 1342 contact and 4141 non-contact surface residues, over all 40 test cases. These numbers decline to 1306 and 4056 respectively when employing the reduced alphabet on the 40 test cases. The ratio of contact to non-contact residues is 0.32 in both the reduced and non-reduced alphabet sets.

Calculating the Shannon Entropy

We use log _e in our MI calculations so that we may compare our results with other MI investigations[23, 27]. In this equation J is a column in the MSA with probabilities P(J = j) for the discrete set of n amino acids jε{1,…,n}. When P(J = j) = 0 then we set P(J = j)log P(J = j) = 0. The entropy is maximal when all j are equally likely to occur, i.e. P(J = j) = 1/n and$H_{\mathit{\text{unstandardised}}}\left(J\right)=-\sum \frac{1}{n}\text{log}\frac{1}{n}=\text{log}n$[13].

where H _{unstandardised} (J) is the entropy of column J in the MSA, and${\overline{H}}_{\mathit{\text{unstandardised}}}$ and${\sigma }_{H_{\mathit{\text{unstandardised}}}}$ are the average entropy and estimated standard deviation, respectively, over all columns in the MSA combined. Our calculated entropies range from 0.0 to 2.8, while the standardised entropies vary from -4.2 to 3.0.

Calculating MI

where column J has n different residues, and column K has m different residues.

The MI is maximal when residues in columns J and K always covary, i.e. P(J = K) = 1 making the$\mathit{\text{MI}}=-\sum _{j=1}^{n}P(J=j)\text{log}P(J=j)$. The maximum MI that can be observed for protein sequences, which have 20 varying residues, is log 20 ≃ 2.9957[13].

where M I _{unstandardised} (J;K) is the MI of columns J and K in the MSA, and${\overline{\mathit{\text{MI}}}}_{\mathit{\text{unstandardised}}}$ and${\sigma }_{M{I}_{\mathit{\text{unstandardised}}}}$ are the average MI and estimated standard deviation respectively, over all interacting domains’ column pairs in the MSA, excluding pairs with an MI value of 0, involving a 0 entropy residue column or a gapped residue column.

Our calculated M I _{unstandardised} scores vary from 0.0 to 1.6, while the standardised MI scores range from -3.2 to 3.7.

Calculating MIp

where${\overline{\mathit{\text{MI}}}}_{\mathit{\text{unstandardised}}}\left(J\right)$ is the average mutual information for column J,${\overline{\mathit{\text{MI}}}}_{\mathit{\text{unstandardised}}}\left(K\right)$ is the average mutual information for column K, and${\overline{\mathit{\text{MI}}}}_{\mathit{\text{unstandardised}}}$ is the overall average mutual information.

where MI p _{unstandardised} (J;K) is the MIp of columns J and K in the MSA and${\overline{\mathit{\text{MIp}}}}_{\mathit{\text{unstandardised}}}$ and${\sigma }_{\mathit{\text{MI}}{p}_{\mathit{\text{unstandardised}}}}$ are the average MIp and estimated standard deviation respectively, over all calculated column pairs in the protein.

Our MIp scores vary from 0.0 to 0.4, while the standardised MIp scores range from -3.1 to 7.0.

Calculating MIc

where MI c _{unstandardised} (J;K) is the MIc of columns J and K in the MSA and${\overline{\mathit{\text{MIc}}}}_{\mathit{\text{unstandardised}}}$ and${\sigma }_{\mathit{\text{MI}}{c}_{\mathit{\text{unstandardised}}}}$ are the average MIc and estimated standard deviation respectively, over all column pairs being considered in the protein.

The MIc scores calculated on our 40 test cases range from -0.02 to 0.1, while the standardised scores range from -2.4 to 7.7.

3-dimensional (3D) MI and MIp

where MSA column J from domain 1 has n different residues, column K from domain 2 has m different residues, and column L from domain 2 has s different residues. Residues in the representative protein structure, corresponding to columns K and L, should be < 4.5Å from each other in order to be considered as being on the same patch in the domain.

where${\overline{\mathit{\text{MI}}3D}}_{\mathit{\text{unstandardised}}}\left(J\right)$ is the average 3D mutual information for column J,${\overline{\mathit{\text{MI}}3D}}_{\mathit{\text{unstandardised}}}\left(K\right)$ is the average 3D mutual information for column K,${\overline{\mathit{\text{MI}}3D}}_{\mathit{\text{unstandardised}}}\left(L\right)$ is the average 3D mutual information for column L, and${\overline{\mathit{\text{MI}}3D}}_{\mathit{\text{unstandardised}}}$ is the overall average 3D mutual information.

In order to compare the 3D mutual information scores between test cases, MI3D and MIp3D scores were standardised in a manner similar to those described in Equations (5), (8) and (12), respectively. Once again MI3D values of 0 were ignored, as were 0 entropy columns and columns containing one or more gaps.

Reduced Alphabet (RA) MI scores

We grouped the 20 amino acids into the same seven physiochemical categories successfully employed by Hamer et al. in their domain-domain contact predictor, i-Patch[12]. These seven categories include: Small (S,G,A,P), Hydrophobic (V,M,I,L,C), Negatively charged (D,E), Aromatic (F,Y,W), Polar (Q,T,N), Favoured Positively-charged (R,H), and Disfavoured Positively-charged (K). These physiochemical groups are abbreviated to S, H, N, A, P, F and D respectively. Hamer et al. introduced Favoured and Disfavoured categories because Lysine (K) was found to be rare in protein/domain interfaces (propensity 0.66), while Arginine (R) and Histidine (H) were far more common (propensities of 1.05 and 1.11, respectively)[12].

We replaced the amino acid alphabets in each MSA by their corresponding category abbreviation and recalculated MI, MIp, MIc, MI3D and MIp3D as described above. The five new MI variant scores are referred to as MIRA, MIpRA, MIcRA, MI3DRA and MIp3DRA.

We choose to employ this particular set of seven physiochemical categories as it was successfully used by i-Patch[12] in domain-domain contact prediction. We do not expect another grouping to dramatically improve the predictive capabilities of MI and its variants further.

P-ROC curves

For classification, each residue in all 40 test cases was assigned the maximum MI score that its residue column achieved with any other residue column in its MSA. When the average score of each residue was assigned instead, the performance of the MI variants decreased significantly, consequently the maximum score was employed.

TP in these equations denote the number of true positives, FP denotes the number of false positives and FN symbolises the number of false negatives.

Each P-ROC plot contains a flat horizontal line (green) at$\frac{\mathit{\text{TP}}}{\mathit{\text{total}}\mathit{\text{scores}}}$. This line denotes the probability of randomly discriminating positive versus negative cases. For example, in Figure2 the solid green line is at 0.245 because there are 1342 “contact” scores out of 5483 total surface scores in the non-reduced alphabet test set. In this figure the dashed green line is at 0.244 because there are 1306 “contact” scores out of 5362 total surface scores in the reduced alphabet test set. Similarly, in Additional file1: Figure S1 the solid green line is at 0.699 as there are 5483 “surface” scores out of 7847 total scores in the non-reduced alphabet test set, while the dashed green line is at 0.712 for there are 5362 “surface” scores out of 7536 total scores in the non-reduced alphabet test set.

Sub-sampling to test stability of MI scores

To test the stability of the 10 MI variant scores under minor changes in the MSA, for each test case 70% of sequences in the MSA were randomly selected and all 10 MI scores recalculated and 10 respective P-ROC curves were plotted. This sub-sampling and calculation process was repeated 100 times per test case for every MI variant. Then the average and standard error of the precision values for the 100 P-ROC curves were calculated for each MI variant. The precision values at 20% recall for each of the MI variants are listed in Table2.

This sub-alignment creation and MI recalculation process was only carried out on those 24 test cases that had ≥200 sequences to ensure that a minimum of 125 sequences were retained in each sub-alignment, the suggested minimum number of sequences required to reduce the stochastic noise in the MSA[8].

MI scores of 0

Pairing any MSA column with a fully conserved column, i.e. a column with an entropy of 0, results in a joint entropy equivalent to the entropy of the non-fully conserved column and subsequently an MI score of 0 for that pair. Since conserved columns do not give any indication of correlated mutations, MI scores involving these columns are ignored. This is standard procedure; for example,[6]. The relationship between percent of columns in an MSA with an entropy of 0 and percent MI scores of 0 computed can be observed in Figure4. This approximately linear relationship further affirms the direct influence a column with an entropy of 0 has on the MI score of pairs involving that column.