Skip to main content
  • Research Article
  • Open access
  • Published:

COUSCOus: improved protein contact prediction using an empirical Bayes covariance estimator



The post-genomic era with its wealth of sequences gave rise to a broad range of protein residue-residue contact detecting methods. Although various coevolution methods such as PSICOV, DCA and plmDCA provide correct contact predictions, they do not completely overlap. Hence, new approaches and improvements of existing methods are needed to motivate further development and progress in the field. We present a new contact detecting method, COUSCOus, by combining the best shrinkage approach, the empirical Bayes covariance estimator and GLasso.


Using the original PSICOV benchmark dataset, COUSCOus achieves mean accuracies of 0.74, 0.62 and 0.55 for the top L/10 predicted long, medium and short range contacts, respectively. In addition, COUSCOus attains mean areas under the precision-recall curves of 0.25, 0.29 and 0.30 for long, medium and short contacts and outperforms PSICOV. We also observed that COUSCOus outperforms PSICOV w.r.t. Matthew’s correlation coefficient criterion on full list of residue contacts. Furthermore, COUSCOus achieves on average 10% more gain in prediction accuracy compared to PSICOV on an independent test set composed of CASP11 protein targets. Finally, we showed that when using a simple random forest meta-classifier, by combining contact detecting techniques and sequence derived features, PSICOV predictions should be replaced by the more accurate COUSCOus predictions.


We conclude that the consideration of superior covariance shrinkage approaches will boost several research fields that apply the GLasso procedure, amongst the presented one of residue-residue contact prediction as well as fields such as gene network reconstruction.


A multiple sequence alignment (MSA) of orthologous protein sequences not only carries evolutionary sequence information, but also information about functional and structural constraints imposed on the three-dimensional (3D) structure of a protein. Conserved or slightly mutated columns indicate important protein positions for protein stability and function. Additionally, non-conserved positions may also play key roles in maintaining the functionality when accompanied by compensatory mutations at other positions [1, 2]. It is of high interest to develop methods predicting coevolution patterns from MSAs, because coevolving positions mainly involve protein positions proximal in 3D structure [3] and they serve as a valuable source of distance constraints in protein structure [47] as well as in protein complex interface predictions [8, 9].

Due to the substantial increase in sequence data in the post-genomic era, a broad range of methods have been introduced for detecting residue-residue contacts from MSAs in the past decades. Mutual information (MI) was one of the first metrics to be applied for contact prediction from MSAs [10, 11]. An improved MI version that corrects for background noise and phylogenetic effects (MIp) has been introduced by Dunn et al. [12]. Recent methodological improvements are able to distinguish between direct and indirect couplings and have demonstrated enormous accuracy in predicting real couplings and coevolution. Such methods include a Bayesian network approach (BvN) [13], Direct Coupling Analysis (DCA) [14, 15], Protein Sparse Inverse COVariance (PSICOV) [16], and pseudolikelihood approaches implemented in plmDCA [17] and GREMLIN [18].

Most recently, new hybrid methods have been developed, amidst many others such as DNCON [19], PConsC [20], CoinDCA [21] or MetaPSICOV [22], where contact detecting methods are combined along with protein physiochemical features to provide more accurate contact predictions.

In the present study, we developed Contact predictiOn Using Shrinked COvariance (COUSCOus), a residue-residue contact detecting method approaching contact inference in a similar manner as PSICOV, by applying the sparse inverse covariance estimation technique introduced by Meinshausen and Bühlmann [23]. Here, we used a different covariance matrix shrinkage approach, the empirical Bayes covariance estimator, which has been shown by Haff to be the best estimator in a Bayesian framework [24], especially dominating estimators of the form aS, such as the smoothed covariance estimator applied in PSICOV. By analysing the original PSICOV benchmark test set [16] and proteins from the Critical Assessment of techniques for protein Structure Prediction 11 (CASP11) experiments, we show that COUSCOus significantly outperforms PSICOV. Furthermore, we designed a simple random forest (RF) meta-classifier that includes contact detecting techniques and sequence-derived physiochemical features and showed that replacing PSICOV with COUSCOus enhances the prediction outcome.



The benchmark dataset used in this study is the original PSICOV test set introduced by Jones et al. [16]. We used the same alignments without modification as have been made available to ensure comparability. However, for validation we selected 146 out of 150 single domain monomeric proteins because the latest PSICOV version V2.1b3 was unable to provide contact predictions on the remaining four test cases, when used with default parameters. This is due to the insufficient number of effective sequences within the four alignments.

The second test set consists of 37 proteins of the CASP11 experiment (see Additional file 1). We selected only those proteins where the latest version of PSICOV successfully provided predictions, to make a fair comparison. The training set introduced by Jones et al. [22] was used to build the RF meta-classifier.

Coevolution analysis methods

The residue-residue contact prediction metrics applied in this study are MI [10, 11], MIp [12], OMES [25], BvN [13], DCA [15] and PSICOV [16]. The resulting coevolution between pairs of amino-acids using MI, MIp and OMES were calculated using Evol module of Prody [26]. BvN results were generated using Perl scripts and C++ source code kindly provided by the authors [13]. PSICOV results were calculated using the code available online [16]. DCA results were obtained using the fast and free software version FreeContact introduced by Kaján et al. [27]. Methodological details for the different methods may be found in the original studies.



In our approach, we generate a sample covariance matrix S from the input MSA. The MSAs are composed of n orthologous protein sequences where each sequence represents a row. Each protein sequence is made of m amino acids as a result of which we have L columns per alignment row. The size of the covariance matrix S is 21L×21L. This is because we compute the marginal single site frequencies f(A i ) and f(B j ) of observed amino acid types (20 natural occurring amino acids and a gap) in columns i and j and their corresponding pair site frequencies f(A i B j ):

$$ S_{ij}^{ab} = f\left(A_{i}B_{j}\right)-f(A_{i})f(B_{j}) $$

Interestingly the precision matrix Θ which is the inverse of the covariance matrix S will contain the partial correlations of all pairs of variables taking into consideration the effects of all other variables. Hence, the non-zero entries of Θ will provide the extent of direct coupling between any two pairs of amino acids at sites i and j.

Yet, due to the fact that we are generating a covariance matrix S out of MSAs representing homologous protein sequences where not all amino acids are present at each site of the MSA, it is certain that S is singular and not directly invertible. Several approaches have been proposed to approximate the precision matrix in such cases. The most powerful and widely used technique is the sparse inverse covariance estimation using the graphical lasso (GLasso) [28].


We briefly summarise the basic motivation and algorithm. Consider matrix X=[X 1,…,X p ] where X i is a random vector of length n with covariance matrix Σ and precision matrix Θ={θ ij }1≤i,jp . Further, let S denote the empirical covariance matrix obtained from the data. The estimation of the precision matrix Θ is challenging when it is sparse. Interestingly, this task is closely related to selection of graphical models.

Let G=(V,E) be a graph representing conditional independence relations between components of X. G is composed of a set of vertices V with p components {X 1,…,X p } and an edge set E of ordered pairs (i,j), with (i,j)E, if an edge between X i and X j exists. The edge between X i and X j is excluded from the edge set E if and only if X i and X j are independent given all other components {X k ,ki,j}. Assuming that the raw data X is multivariate gaussian (XN(μ,Σ)), the conditional independence between X i and X j given all other components is equivalent to zero in the precision matrix (θ ij =0) as shown in [29]. Hence, for gaussian distributions recovering the structure of graph G is equivalent to the estimation of the support of the precision matrix.

The precision matrix Θ can then be estimated using a L 1 penalised log-likelihood approach. The GLasso algorithm, introduced by Friedman et al. [28], efficiently computes the solution by:

$$ \hat{\Theta}_{\text{GLasso}} := \arg\min_{\Theta \succ 0}\left\{\text{tr}(S\Theta) - \log\det(\Theta) + \lambda\|\Theta\|_{1}\right\}, $$

with tr as trace, Θ1 as the sum of the absolute values of the elements in Θ and λ as a positive tuning parameter to control the sparsity.

Empirical Bayes covariance estimator

The most natural estimator of the covariance Σ is the sample covariance matrix S. The estimator is optimal in the classical settings with large number of samples and fixed low dimensions (n>p). However, it performs poorly in high-dimensional settings (n<<p), see Johnstone [30]. The GLasso approach operates very well in this context, but the computational time required to reach convergence can be large in some cases such as for protein families with low number of sequences. As an alternative to the natural estimator S, several shrinkage estimators have been proposed in the literature [31, 32]. They take a weighted average of the sample covariance matrix S, with a suitable chosen target diagonal matrix. Jones et al. applied a smoothed covariance estimator that shrinks the matrix towards the shrinkage target \(F=diag(\bar {S},\bar {S},\ldots,\bar {S})\) [16]. In this work, we applied the empirical Bayes estimator proposed by Haff [24]:

$$ \hat{\Sigma}=S+\frac{p-1}{n\ tr(S)}I_{p}, $$

where I p represents the identity matrix of order p. In a Bayesian framework, it has been proven by Haff that this estimator is the best estimator of the form a(S+u t(u)C), with 0<a<1/(n−1) and u=1/tr(S −1 C). Here t(·) is non-increasing and C an arbitrary positive definite matrix. It dominates estimators of the form aS by a substantial amount. More precisely, it has been proven that under the L 2 loss, the uniform reduction in the risk function is at least \(100\frac {p+1}{n+p}\)%. In this study, we performed the shrinkage until the adjusted covariance matrix \(\hat {\Sigma }\) is no longer singular, i.e. is a positive-definite matrix. The adjusted covariance matrix \(\hat {\Sigma }\) was finally applied in the GLasso algorithm to obtain the sparse precision matrix, that contains the degree of direct coupling between any pair of amino acids.

APC correction

The coevolution pair list is generated identically to the PSICOV final processing. For each MSA column pair i and j we compute the L 1-norm out of the corresponding 20×20 submatrix in Θ (only contributions of the 20 amino acid types are considered):

$$ C_{ij} = \sum_{ab}\mid \Theta_{ij}^{ab}\mid $$

Furthermore, we adjust the coupling score by the average product correction (APC), previously applied for MI by Dunn et al. [12] to reduce entropic and phylogenetic bias:

$$ {COUSCOus}_{ij} = C_{ij} - \frac{\bar{C}_{(i-)}\bar{C}_{(-j)}}{\bar{C}} $$

with \(\bar {C}_{(i-)}\) as the mean precision norm of column i and all other columns, \(\bar {C}_{(-j)}\) as the corresponding for column j and \(\bar {C}\) as the mean precision norm of all coupling scores.

Random forest meta-classifier

As previously indicated, new hybrid methods that combine coevolution detecting tools with other sources of information such as protein physiochemical features outperforms single methods like PSICOV. However, we are convinced that improvements of single residue-residue contact detecting methods can boost new emerging hybrid techniques. We designed a RF meta-classifier that includes several contact prediction methodologies along with a small number of sequence-derived physiochemical features.

In particular, we built a RF classifier using the training set alignments from the MetaPSICOV study [22]. In total, we used 336 protein alignments where PSICOV was able to successfully provide contact predictions. The RF was trained using the following features:

  • Contact detecting methodologies MI, MIp, BvN, PSICOV or COUSCOus, FreeContact and CCMpred

  • Secondary structure and solvent exposure probabilities derived from PSIPRED [33]

  • Shannon entropy using R [34] package bio3d [35]

  • Hydrophobicity using R package Interpol [36]

  • Amino acid physiochemical properties

The RF meta-classifier was trained using 500 trees with ten features and a max-depth of eight. We performed five-fold cross-validation while training the classifier and optimised the area under the curve (AUC) metric for performance.

Evaluation metrics and distance descriptions

The problem of predicting protein residue-residue contacts is well-known to be an extremely difficult one as on average only 3% of all possible residue pairs in known protein structures are identified to be real contacts. In the latest CASP11 challenge [37], this problem was tackled by dividing the contact prediction task into two categories. First, evaluation of predicted contacts using quality metrics like Accuracy and X d [38] on reduced lists (RL). Second, evaluation of predicted contacts using quality metrics like Matthew’s correlation coefficient (MCC) [39] and area under the precision-recall (A U C pr ) for full lists (FL). The RL are usually defined by considering the top L/n predicted contacts where L is the length of the evaluation target or protein sequence and n is a small integer (e.g. 1,5 or 10). RL metric accuracy is calculated as \(\frac {TP}{TP+FP}\) where TP defines a correctly predicted contact and FP an incorrectly predicted contact. The second RL metric X d represents the difference between the distance distributions of the predicted contacts and all pairs distance distributions in the 3D target structure. It is defined as \(X_{d} = \sum \limits _{i=1}^{15} \frac {Pip-Pia}{di*15}\), with Pia and Pip are the percentages of pairs included in the i th bin for the whole target and predicted contacts respectively. Additional details can be found in [38]. The FL metrics used in this study are A U C pr as it is a robust metric for unbalanced classes and the Matthew’s correlation coefficient [39] to evaluate all residue pairs for contact prediction.

Results and discussion

We first illustrate as an example in Fig. 1 the spatial proximity of the predicted contacts obtained by PSICOV (a and b) and COUSCOus (c and d) for the Immunoglobulin V-set domain (Protein family database [40] (PFAM) ID: PF07686). The upper triangles of the presented contact maps (Fig. 1 a and c) display the native contacts. A residue-residue pair is hereby considered to be in contact if the two amino acids are proximal in the 3D structure, in particular if their C β - C β (C α in the case of glycine) distance is less than 8 Å ngström (Å). The lower triangles show the L/2 contact predictions (in this case 62) obtained by using either PSICOV or COUSCOus. Correctly predicted contacts are coloured in green and incorrect ones in red. Further, we mapped the top L/2 predictions to the structure of the myelin oligodendrocyte glycoprotein (Protein Data Bank [41] (PDB) ID: 1PKO), solved at 1.45 Å resolution (Fig. 1 b and d). Out of 62 possible contacts, COUSCOus correctly predicted 49 (accuracy: 0.79) compared to 39 (accuracy: 0.63) by PSICOV resulting in higher accuracy of COUSCOus. The figure shows (b) that the incorrect identified pairs are mainly located in loop regions at the top and the bottom. Hence, the pairs may still have distances less than 8Å considering that the unordered regions are not static as illustrated by a crystal structure. In contrast, incorrect predicted pairs from PSICOV (Fig. 1 d), are distributed over the entire protein.

Fig. 1
figure 1

The top L/2 (in this case 62) long, medium and short contact predictions for the Immunoglobulin V-set domain family (PFAM ID: PF07686) obtained using PSICOV (a and b) and COUSCOus (c and d) and mapped to the myelin oligodendrocyte glycoprotein 3D crystal structure (PDB ID: 1PKO) (right panel). Correctly predicted contacts are shown in green and incorrect ones in red. Upper triangles of the contact maps display all the native C β C β contacts (left panel). The lower triangles show contacts predicted by PSICOV (a and b) and COUSCOus (c and d)

The performance of COUSCOus was evaluated on the original PSICOV benchmark test set using standard assessment metrics applied in CASP: accuracy, X d ,M C C and A U C pr (see Methods). We distinguished between three types of contacts: short (6 ≤ residue separation < 12), medium (12 ≤ residue separation < 24) and long range contacts (residue separation ≥ 24).

We list in Table 1 the mean accuracies of COUSCOus and PSICOV for the top- L/10, top- L/5 and top-L for long, medium and short range predicted contacts on the original PSICOV benchmark set. For the top- L/10 predictions COUSCOus is 10, 8 and 13% more accurate than PSICOV for long, medium and short range contacts. Similarly, for the top- L/5 predictions we observed the similar gains in accuracy when using COUSCOus. For the top-L predictions we observed different accuracies for the three types of contacts. The gain in accuracy of 11% when using COUSCOus was similar for the long range contacts but the gain dropped to 6 and 5% for medium and short range contacts.

Table 1 Mean accuracies of COUSCOus vs. PSICOV on PSICOV benchmark dataset

The second evaluation metric applied in this study, the X d score, estimates the deviation of the distribution of distance in the RL sets (L/10,L/5 or L) of contacts from the distribution of the distances in all residue pairs within the protein (see Methods). Table 2 summarises the average X d scores for COUSCOus and PSICOV. For the top-L predictions COUSCOus is more accurate than PSICOV on long, medium and short range contacts (16, 10 and 18%). For the top- L/10 and top- L/5 predictions we observed smaller improvements in X d score ranging from 1 to 11%. Moreover, we compared the performance of COUSCOus and PSICOV on FL; considering all possible residue pairs. In Table 3 we summarise the mean A U C pr values of the precision recall (PR) curves for COUSCOus and PSICOV. COUSCOus outperforms PSICOV with gains in accuracy of 14, 11 and 11% for long, medium and short range contacts, respectively.

Table 2 Mean X d values of COUSCOus and PSICOV on PSICOV benchmark dataset
Table 3 Mean AUC pr values

We also performed an exhaustive analysis of COUSCOus and PSICOV w.r.t. MCC for long, medium and short range contact predictions. In Figure 2 we illustrate via box plot the distribution of the MCC values for the two prediction methods. It is apparent from the box plots that COUSCOus is superior in predicting real residue contacts. To further test for significance we performed a t-test, after successfully testing for normal distribution and variance homogeneity, on the MCC distributions for different contact ranges. COUSCOus outperforms PSICOV on all types of contacts significantly, with P-values of 6×10−7,1.2×10−2 and 6×10−4 for long, medium and short range contacts, respectively.

Fig. 2
figure 2

MCC distributions for PSICOV benchmark proteins in case of long, medium and short range contacts predicted by PSICOV and COUSCOus. The stars represent statistical significance where is used to represent P-value < 0.05 and is used to represent P-values < 0.001

Next, we analysed the dependence of the performance of PSICOV and COUSCOus with regard to the size of the protein family. In this case, we used the number of effective sequences in a MSA as comparison metric to account for the fact that highly similar homologous do not provide any additional contact information than a single one. Similar to Ma et al. [21] we grouped the test set members into five categories by l n(N eff ): [4,5),[5,6),[6,7),[7,8),[8,10), and calculated the averaged L/10 accuracies for each group. Figure 3 shows clearly that COUSCOus outperforms PSICOV regardless of the l n(N eff ) on long, medium and short range contacts. In addition, we tested the performance of COUSCOus and PSICOV on an independent test set from the latest CASP11 experiment. In Table 4 we show the mean accuracies of COUSCOus and PSICOV for the top- L/10, top- L/5 and top-L for long, medium and short range predicted contacts. For the top- L/10 predictions COUSCOus is 10, 9 and 13% more accurate than PSICOV for long, medium and short range contacts, respectively. Similarly, for the top- L/5 predictions we observed gains in accuracy of 13, 10 and 9% when using COUSCOus. For the top-L predictions we observed different accuracies for the three types of contacts. A gain in accuracy of 13% is observable when using COUSCOus for the long range contacts but the gain dropped to 4 and 8% for medium and short range contacts.

Fig. 3
figure 3

Dependence of the performance of PSICOV and COUSCOus on the effective number of sequences (N eff ) in the MSAs. The performance is evaluated using accuracies for the top L/10long, medium and short contacts. The solid line represents the averaged accuracies of the test set binned into five different categories of Neff (ln(Neff): [4, 5), [5, 6), [6, 7), [7, 8), [8, 10)). COUSCOus outperforms PSICOV independent of the l n(N eff ) in the test set

Table 4 Mean accuracies of COUSCOus vs. PSICOV on the CASP11 benchmark dataset

Next, we designed two experiments using a RF meta-classifier where we combine contact detecting tools with protein sequence derived features. In the first case we used predictions of PSICOV as a feature in the meta-classifier and in the second case we replaced PSICOV predictions with COUSCOus predictions. The classifier including COUSCOus results as a feature outperforms the classifier including PSICOV results for the top- L/10 and top- L/5 predictions on all types of contacts except for the top-L predictions for long range contacts (see Table 5).

Table 5 Mean accuracies of RF meta-classifier including COUSCOus or PSICOV as a feature on the PSICOV benchmark dataset

In Additional file 2 we illustrate the mean accuracies of different contact detecting techniques. Our newly developed technique COUSCOus (green upper triangle) outperforms PSICOV (black points) and is equally well as FreeContact (red lower triangle) on all contact types. The best single contact detecting tool is CCMpred (magenta rectangle). Our simple RF meta-classifier that combines 6 single residue-residue contact detecting tools along with 4 sequence-derived features outperforms all single methods. However, MetaPSICOV, a multi-stage neural network hybrid method that combines five coevolution techniques along with a broad range of sequence-derived features is still the best performing method.


In this work, we assessed the performance of our newly developed method COUSCOus in predicting residue-residue contacts from MSAs. In particular, the performance was tested in comparison to PSICOV on the original PSICOV benchmark test set as well as on CASP11 targets. On the RL sizes COUSCOus outperformed PSICOV substantially, with on average 10% gain in prediction accuracy for all types of contacts. Moreover, COUSCOus proved to be superior over PSICOV on FL sizes with average A U C pr gains of 12%. With regard to MCC scores COUSCOus is even significantly outperforming PSICOV, illustrated with the help of box plots and hypothesis tests (see also Additional file 3). Further, we reported that COUSCOus’s gain in accuracy is independent of the number of effective sequences in a given MSA.

The main motivation of this work was to highlight that improvements of single residue-residue contact detecting tools, in this case PSICOV, might lead to improvements of new hybrid methods that combine contact detecting techniques with physiochemical and other sequence derived protein features. As proof of concept, we showed with the help of a simple RF meta-classifier that PSICOV should be replaced in hybrid classifiers by COUSCOus.


Jones et al. [16] demonstrated in their initial work that GLasso in principle performs excellently in identifying directly coupled columns within a MSA. In the present study, we highlighted that the application of a different shrinkage approach than the one used in PSICOV, the empirical Bayes covariance estimator, in combination with GLasso substantially increased the contact precision. The theoretically shown superiority of the empirical Bayes covariance estimator over simpler smoothed covariance estimators of the form aS is also valid within this application of contact detection from MSAs.

Furthermore, it is worth mentioning that other research fields that apply the GLasso procedure, such as gene network reconstruction, may also benefit from applying other shrinkage techniques. We are keen to investigate the effect of shrinkage in other graphical inference problems in future work.

Another important application that we are keen to investigate in future is the de novo structure prediction of proteins or protein complexes using COUSCOus or a hybrid classifier, including COUSCOus contact predictions as distance constraints, similarly to what have been applied in EVFold [4] and EVComplex [8], or PconsFold [42].





Å ngström


area under the curve


Bayesian network approach


Contact predictiOn Using Shrinked COvariance


Critical Assessment of techniques for protein Structure Prediction 11


Direct Coupling Analysis


full lists


graphical lasso


Matthew?s correlation coefficient


mutual information


corrected mutual information


multiple sequence alignment


Protein Data Bank


Protein family database


precision recall


Protein Sparse Inverse COVariance


Random forest


reduced lists


  1. Yanofsky C, Horn V, Thorpe D. Protein structure relationships revealed by mutual analysis. Science (New YorkNY). 1964; 146:1593–4.

    Article  CAS  Google Scholar 

  2. Fitch WM, Markowitz E. An improved method for determining codon variability in a gene and its application to the rate of fixation of mutations in evolution. Biochem Genet. 1970; 4(5):579–93.

    Article  CAS  PubMed  Google Scholar 

  3. de Juan D, Pazos F, Valencia A. Emerging methods in protein co-evolution. Nat Rev Genet. 2013; 14(4):249–61.

    Article  CAS  PubMed  Google Scholar 

  4. Marks DS, Colwell LJ, Sheridan R, Hopf TA, Pagnani A, Zecchina R, et al.Protein 3D Structure Computed from Evolutionary Sequence Variation. PLoS ONE. 2011; 6(12):e28766.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Marks DS, Hopf TA, Sander C. Protein structure prediction from sequence variation. Nat Biotechnol. 2012; 30(11):1072–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Hopf TA, Colwell LJ, Sheridan R, Rost B, Sander C, Marks DS. Three-dimensional structures of membrane proteins from genomic sequencing. Cell. 2012; 149(7):1607–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Kosciolek T, Jones DT. De Novo Structure Prediction of Globular Proteins Aided by Sequence Variation-Derived Contacts. PLoS ONE. 2014; 9(3):e92197.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Hopf TA, Schärfe CPI, Rodrigues JPGLM, Green AG, Kohlbacher O, Sander C, et al.Sequence co-evolution gives 3D contacts and structures of protein complexes. eLife. 2014; 3:e03430.

    Article  PubMed Central  Google Scholar 

  9. Ovchinnikov S, Kamisetty H, Baker D. Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information. eLife. 2014; 3:e02030.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Gloor GB, Martin LC, Wahl LM, Dunn SD. Mutual information in protein multiple sequence alignments reveals two classes of coevolving positions. Biochemistry. 2005; 44(19):7156–65.

    Article  CAS  PubMed  Google Scholar 

  11. Martin LC, Gloor GB, Dunn SD, Wahl LM. Using information theory to search for co-evolving residues in proteins. Bioinformatics (Oxford England). 2005; 21(22):4116–24.

    Article  CAS  Google Scholar 

  12. Dunn SD, Wahl LM, Gloor GB. Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction. Bioinformatics (Oxford England). 2008; 24(3):333–40.

    Article  CAS  Google Scholar 

  13. Burger L, van Nimwegen E. Disentangling direct from indirect co-evolution of residues in protein alignments. PLoS Comput Biol. 2010; 6(1):e1000633.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Weigt M, White RA, Szurmant H, Hoch JA, Hwa T. Identification of direct residue contacts in protein-protein interaction by message passing. Proc Natl Acad Sci USA. 2009; 106(1):67–72.

    Article  CAS  PubMed  Google Scholar 

  15. Morcos F, Pagnani A, Lunt B, Bertolino A, Marks DS, Sander C, et al.Direct-coupling analysis of residue coevolution captures native contacts across many protein families. Proc Natl Acad Sci USA. 2011; 108(49):E1293—301.

    Article  PubMed  Google Scholar 

  16. Jones DT, Buchan DWA, Cozzetto D, Pontil M. PSICOV: precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments. Bioinformatics (Oxford England). 2012; 28(2):184–90.

    Article  CAS  Google Scholar 

  17. Ekeberg M, Lȯvkvist C, Lan Y, Weigt M, Aurell E. Improved contact prediction in proteins: using pseudolikelihoods to infer Potts models. Phys Rev E Stat Nonlinear Soft Matter Phys. 2013; 87(1):012707.

    Article  Google Scholar 

  18. Kamisetty H, Ovchinnikov S, Baker D. Assessing the utility of coevolution-based residue-residue contact predictions in a sequence- and structure-rich era. Proc Natl Acad Sci USA. 2013; 110(39):15674–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Eickholt J, Cheng J. Predicting protein residue-residue contacts using deep networks and boosting. Bioinformatics. 2012; 28(23):3066–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Skwark MJ, Abdel-Rehim A, Elofsson A. PconsC: combination of direct information methods and alignments improves contact prediction. Bioinformatics. 2013; 29(14):1815–6.

    Article  CAS  PubMed  Google Scholar 

  21. Ma J, Wang S, Wang Z, Xu J. Protein contact prediction by integrating joint evolutionary coupling analysis and supervised learning. Bioinformatics. 2015; 31(21):3506–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Jones DT, Singh T, Kosciolek T, Tetchner S. MetaPSICOV: combining coevolution methods for accurate prediction of contacts and long range hydrogen bonding in proteins. Bioinformatics. 2015; 31(7):999–1006.

    Article  CAS  PubMed  Google Scholar 

  23. Meinshausen N, Bühlmann P. High-dimensional graphs and variable selection with the Lasso. Ann Stat. 2006; 34(3):1436–62.

    Article  Google Scholar 

  24. Haff LR. Empirical Bayes Estimation of the Multivariate Normal Covariance Matrix. Ann Stat. 1980; 8(3):586–97.

    Article  Google Scholar 

  25. Kass I, Horovitz A. Mapping pathways of allosteric communication in GroEL by analysis of correlated mutations. Proteins Struct Funct Genet. 2002; 48(4):611–7.

    Article  CAS  PubMed  Google Scholar 

  26. Bakan A, Dutta A, Mao W, Liu Y, Chennubhotla C, Lezon TR, et al.Evol and ProDy for bridging protein sequence evolution and structural dynamics. Bioinformatics. 2014; 30(18):2681–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kaján L, Hopf TA, Kalas̆ M, Marks DS, Rost B. FreeContact: fast and free software for protein contact prediction from residue co-evolution. BMC Bioinforma. 2014; 15(1):85.

    Article  Google Scholar 

  28. Friedman J, Hastie T, Tibshirani R. Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 2008; 9(3):432–41.

    Article  PubMed  Google Scholar 

  29. Lauritzen SL. Graphical Models, 1st ed. Oxford: Oxford University Press; 1996.

    Google Scholar 

  30. Johnstone IM. On the Distribution of the Largest Eigenvalue in Principal Components Analysis. Ann Stat. 2001; 29(2):295–327.

    Article  Google Scholar 

  31. James W, Stein C. Estimation with quadratic loss. In: Proc. Fourth Berkeley Symp. Math. Statist. Prob. Berkeley: University of California Press: 1961. p. 361–379.

    Google Scholar 

  32. Ledoit O, Wolf M. Improved estimation of the covariance matrix of stock returns with an application to portfolio selection. J Empir Financ. 2003; 10(5):603–21.

    Article  Google Scholar 

  33. Jones DT. Protein secondary structure prediction based on position-specific matrices. J Mol Biol. 1999; 292:195–202.

    Article  CAS  PubMed  Google Scholar 

  34. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna; 2014.

  35. Grant BJ, Rodrigues APC, ElSawy KM, McCammon JA, Caves LSD. Bio3d: an R package for the comparative analysis of protein structures. Bioinformatics. 2006; 22(21):2695–6.

    Article  CAS  PubMed  Google Scholar 

  36. Heider D, Hoffmann D. Interpol: An R package for preprocessing of protein sequences. BioData Min. 2011; 4(1):16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Park H, DiMaio F, Baker D. CASP11 refinement experiments with ROSETTA. Proteins. 2016; 84:1097–0134.

    Article  Google Scholar 

  38. Izarzugaza JMG, Graṅa O, Tress ML, Valencia A, Clarke ND. Assessment of intramolecular contact predictions for CASP7. Proteins Struct Funct Bioinforma. 2007; 69(S8):152–8.

    Article  CAS  Google Scholar 

  39. Matthews BW. Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochim Biophys Acta Protein Struct. 1975; 405(2):442–51.

    Article  CAS  Google Scholar 

  40. Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, et al.Pfam: the protein families database. Nucleic Acids Res. 2014; 42(D1):D222—D230.

    Google Scholar 

  41. Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al.The Protein Data Bank. Nucleic Acids Res. 2000; 28(1):235–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Michel M, Hayat S, Skwark MJ, Sander C, Marks DS, Elofsson A. PconsFold: improved contact predictions improve protein models. Bioinformatics. 2014; 30(17):i482—i488.

    Article  PubMed Central  Google Scholar 

Download references


We thank Erik van Nimwegen and Lukas Burger for providing code to implement their Bayesian network approach.


Not applicable.

Availability of data and materials

Implementation of COUSCOus is available as R package ( Raw data (alignments) and other applied software can be found as referred in the Methods section.

Authors’ contributions

RR contributed to the concept, programmed COUSCOus, coordinated the project and wrote the manuscript; RM programmed the RF and assisted in the analysis and writing; KK, MA, and EU assisted in the analysis, reviewed and edited the manuscript; ME and HB contributed to the concept, reviewed and edited the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Reda Rawi.

Additional files

Additional file 1

A complete list of all CASP11 protein codes applied in this study. (PDF 72.9 kb)

Additional file 2

Mean accuracies of several contact detecting methods for the top- L/10,L/5 and L predictions for all contact types applying the PSICOV benchmark dataset. (PDF 85.6 kb)

Additional file 3

Scatterplot comparing the accuracies of the top L contacts of PSICOV to COUSCOus, using sequence separation ≥6. (PDF 78 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rawi, R., Mall, R., Kunji, K. et al. COUSCOus: improved protein contact prediction using an empirical Bayes covariance estimator. BMC Bioinformatics 17, 533 (2016).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: