---
abstract: |
    One of the most challenging and long-standing problems in computational biology is the prediction of three-dimensional protein structure from amino acid sequence. A promising approach to infer spatial proximity between residues is the study of evolutionary covariance from multiple sequence alignments, especially in light of recent algorithmic improvements and the fast growing size of sequence databases. In this paper, we present a simple, fast and accurate algorithm for the prediction of residue-residue contacts based on regularized least squares. The basic assumption is that spatially proximal residues in a protein coevolve to maintain the physicochemical complementarity of the amino acids involved in the contact. Our regularized inversion of the sample covariance matrix allows the computation of partial correlations between pairs of residues, thereby removing the effect of spurious transitive correlations. The method also accounts for low number of observations by means of a regularization parameter that depends on the effective number of sequences in the alignment. When tested on a set of protein families from Pfam, we found the RLS algorithm to have performance comparable to state-of-the-art methods for contact prediction, while at the same time being faster and conceptually simpler.

    The source code and data sets are available at <http://cms.dm.uba.ar/Members/slaplagn/software>
author:
- Massimo Andreatta
- Santiago Laplagne
- Shuai Cheng Li
- Stephen Smale
bibliography:
- 'mybib\_plain.bib'
date: 'March 26, 2014'
title: 'Prediction of residue-residue contacts from protein families using similarity kernels and least squares regularization'
---

Introduction {#section:introduction}
============

A major problem in computational biology is the prediction of the 3D structure of a protein from its amino acid sequence. Anfinsen’s dogma suggests that, in principle, the amino acid sequence contains enough information to determine the full three-dimensional structure [@anfinsen1973principles]. However, a few decades on, the mechanisms of protein folding are still not satisfactorily explained [@dill2012protein]. In particular, the space of possible spatial configurations given a certain amino acid 1D sequence is immense (the “Levinthal paradox”), yet an unfolded polypeptide chain is driven to its native 3D structure in a finite time, typically milliseconds to seconds, upon shifting to folding conditions [@rose2006backbone].

Such enormous search space poses important challenges to the development of *ab initio* methods for structure prediction. Therefore, it is essential to exploit different kinds of information that can help reduce the degrees of freedom in the configurational search space. A powerful way of inferring distance constraints is the prediction of residue-residue contacts from multiple sequence alignments (MSA). The underlying assumption is that contacting residues coevolve to maintain the physicochemical complementarity of the amino acids involved in the contact. That is, if a mutation occurs in one of the contacting residues, the other one is also likely to mutate, lest the fold of the protein may be disrupted. Methods based on residue coevolution aim at inferring spatial proximity between residues (contacts) from such signals of correlated mutations (Figure \[fig:contacts\]).

![Illustration of a residue-residue contact. The contact imposes a constraint on the evolution of residues $i$ and $j$. Vice versa, coevolution of $i$ and $j$ can be used to infer their physical proximity.[]{data-label="fig:contacts"}](Contact_diagram-eps-converted-to.pdf){width="45.00000%"}

Thanks to the recent exponential growth in sequence data collected in databases such as Pfam [@PFAM], algorithms for the prediction of contacting residues from MSA have enjoyed increasing attention. Different kinds of approaches have been recently applied for contact prediction, from mutual information (MI) between pairs of positions [@buslje2009correction; @dunn2008mutual; @wang2013predicting], to Bayesian network models [@burger2010disentangling], direct-coupling analysis [@balakrishnan2011learning; @morcos2011direct; @marks2011protein] and sparse inverse covariance matrix estimation [@jones2012psicov]. See also [@marks2012protein] and [@de2013emerging] for recent reviews. In particular, the more sophisticated and successful methods attempt to disentangle direct and indirect correlations, that is the artifactual correlations emerging from transitive effects of covariance analysis [@lapedes1999correlated; @weigt2009identification]. Morcos et al. [@morcos2011direct] and Marks et al. [@marks2011protein] tackle this problem using a maximum-entropy approach, whereas Jones et al. [@jones2012psicov] estimate partial correlations by inverting the covariance matrix. A very recent pseudo-likelihood method based on 21-state Potts models [@ekeberg2013improved] was shown to outperform other approaches for direct-coupling analysis. Kamisetty et al. [@kamisetty2013assessing] systematically analyzed the conditions under which predicted contacts are likely to be useful for structure prediction, and found several hundred families that meet their criteria.

Here, we propose a new approach for computing direct correlations that employs regularized least squares (RLS) regression to invert a sample covariance matrix $S$. We compute the regularized inverse by the formula $$\label{eq:theta}
 \Theta = (S^2 + \eta \operatorname{Id})^{-1} S,$$ with fixed $\eta > 0$. It proves to be a very simple, direct and fast approach, and requires no assumption on probabilities distributions or sparsity in the correlations.

The RLS algorithm described in this paper was applied to three different sets of protein families, and we compared its performance to state-of-the-art methods for contact prediction. The RLS method achieves precision rates superior to PSICOV [@jones2012psicov] and comparable to plmDCA [@ekeberg2013improved] but it is considerably faster than either.

Approach
========

The covariance matrix {#section:covariance}
---------------------

Let ${\mathscr{A}}$ be the set of $20$ amino acids plus the gap symbol $-$ and $\PP = \{p^m = (p_1^m, \dots, p_L^m)\}_{m=1,\dots,M}$ a given Pfam family of $M$ aligned protein sequences, possibly with gaps, where $L$ denotes the length of the protein domains. On this set of proteins, the covariance between any pair of columns $(i,j)$ for the amino acids pair $(a,b)$ is given by $$S^0_{ij}(a,b) = f_{ij}(a,b) - f_i(a) f_j(b)$$ where the corrected frequencies are calculated as $$\label{eq:pseudo}
    f_i (a) = \frac{1}{\lambda + M_{\text{eff}}} \Big( \frac{\lambda}{21} + \sum_p w(p) \delta(a, p_i) \Big)$$ $$f_{ij} (a,b) = \frac{1}{\lambda + M_{\text{eff}}} \Big( \frac{\lambda}{21^2} + \sum_p w(p) \delta(a, p_i) \delta(b, p_j) \Big)$$ The delta kernel takes value $\delta(a,b)=1$ if $a = b$ and $\delta(a,b)=0$ otherwise. $w(p)$ is the weight of protein $p$ and ${M_{\text{eff}}}= \sum_p w(p)$ (see section \[section:measure\] for details on sequence weighting).

The parameter $\lambda$ is the so-called pseudocount, a regularization parameter that accounts for non-observed pairs. We note that the same, or similar, constructions for the corrected amino acid frequencies have been proposed previously by other authors [@ekeberg2013improved; @jones2012psicov; @morcos2011direct].

### Modified covariance matrix

We set $S^0_{ii}(a,b) = 0$ for $a \neq b$, and call $S$ this new matrix. This modification also appears in the code of PSICOV [@jones2012psicov] although it is not stated in their paper. By setting those values to $0$, the resulting matrix contains in general negative eigenvalues (see Figures and ) and hence is not anymore semi-definite positive, but it is still symmetric. We do not fully understand this step, but it is noteworthy that Equation \[eq:theta\] still makes sense for any $\eta > 0$.

In general, working with $S$ instead of $S^0$ gives better results in our experiments. See Table S1 for the effect of this step on predictive performance.

Regularized inverse – the key algorithm
---------------------------------------

As we mentioned in the Introduction, the covariance between our random variables does not distinguish between direct and indirect correlations. To overcome this problem, a technique used by statisticians is to compute the so-called partial correlations, which can be obtained from the inverse of the covariance matrix using its associated correlation matrix.

Since the covariance matrix is usually singular or ill conditioned, regularization techniques must be used to compute a regularized inverse $\Theta$. We achieve this by solving the following optimization problem $$\label{eq:optimization}
 \Theta = \operatorname*{argmin}_{X \in {\mathbb{R}}^{20L \times 20L}} \| SX - \operatorname{Id}\|_2^2 + \eta \|X\|_2^2,$$ where $\| \cdot \|_2$ denotes the Frobenius norm, and $\eta$ is a regularization parameter to be determined. Observe that the first term is minimized by the inverse of $S$ when it exists.

The problem has a unique solution for any $\eta > 0$ as we see in the next proposition.

\[prop:regularization\] For a symmetric matrix $S \in {\mathbb{R}}^{n \times n}$ and a regularization parameter $\eta > 0$, the optimization problem has a unique solution, which is also symmetric and given by equation \[eq:theta\]. When $S$ is semidefinite positive, then the solution also is.

Since the norms involved are coordinate norms, the problem can be decoupled into independent problems for each column of $X$: $$\Theta^{(i)} = \operatorname*{argmin}_{x \in {\mathbb{R}}^{n \times 1}} \| S^tx - e^{(i)} \|_2^2 + \eta \|x\|_2^2,$$ where $\Theta^{(i)}$ is the $i$-th column of $\Theta$ and $e^{(i)}$ is the $i$-th column of the identity matrix.

This is a well studied problem known as regularized least squares (also called Tikhonov regularization or Ridge regression in different areas, see [@tikhonov1943stability] and [@hoerl1962application]). The unique solution is $\Theta^{(i)} = (S^tS + \eta \operatorname{Id})^{-1} S^t e^{(i)}$.

Hence, the solution to our matrix problem is $\Theta = (S^tS + \eta \operatorname{Id})^{-1} S^t$. Since we are assuming $S$ symmetric, we get $$\Theta = (S^2 + \eta \operatorname{Id})^{-1} S.$$

The matrix $S$ is diagonalizable with all of its eigenvalues real. The eigenvalues of $S$ are transformed by the same formula defining $\Theta$. If $\lambda_k$, $1 \le k \le 20L$, are the eigenvalues of $S$ then the eigenvalues of $\Theta$ will be $$\gamma_k = f(\lambda_k) = \frac{\lambda_k}{\lambda_k^2 + \eta}$$ This function is well defined for all $\lambda \in {\mathbb{R}}$ when $\eta$ is positive, which proves that the matrix $S^2 + \eta \operatorname{Id}$ is invertible. The resulting matrix $\Theta$ is symmetric by standard matrix theory. Finally, $f$ preserves the sign of the eigenvalue and hence $\Theta$ will be a semidefinite positive matrix whenever $S$ is.

Note that $\Theta$ can be computed by solving the linear system $(S^2 + \eta \operatorname{Id}) \Theta = S$, which is faster and more accurate than inverting the matrix $S^2 + \eta \operatorname{Id}$.

For a better understanding of our regularization formula, we study the function $f$ in more detail. The derivative of $f$ is $f'(\lambda) = \frac{-\lambda^2 + \eta}{(\lambda^2+\eta)^2}$. Hence $f$ is increasing for $|\lambda| < \sqrt{\eta}$ and decreasing for $|\lambda| > \sqrt{\eta}$, with maximum value at $\lambda = \sqrt{\eta}$ and minimum value at $\lambda = -\sqrt{\eta}$. We show in Figure the plot of this function for $\eta = \eta'/{M_{\text{eff}}}= 1000/3912$ (see Section \[subsec:regularization\] for the choice of $\eta$).

As mentioned in the proof of Proposition \[prop:regularization\], the function is smooth at $0$, so using this regularization formula we deal in a simple way with the conditioning problem of inverting the covariance matrix.

Aggregation
-----------

The matrix $\Theta$ obtained is a $20L \times 20L$ matrix. Its entries are estimates of the partial correlation between pairs of columns $(i,j)$ for [*all*]{} pairs of amino acids $(a,b)$. Since our goal is to detect relations between pairs of columns in the alignment, we compute a coupling score aggregating the values of $\Theta$ using the $l_1$-norm on the $20 \times 20$ sub-matrices, as in [@jones2012psicov]. That is, $$P(i, j) = \sum_{1 \le a, b \le 20} | \Theta_{ij}(a, b) |.$$

The $l_2$-norm for aggregation showed poorer performance than the $l_1$-norm above. Finally, following[@dunn2008mutual] and [@jones2012psicov] we define a corrected score $P_{\text{APC}}(i, j) = P(i, j) - \frac{P(\cdot, j)P(i, \cdot)}{P(\cdot, \cdot)}$, where $\cdot$ stands for the average over all positions.

The prediction of contacts between pairs of residues can now be obtained by ranking the $P_{\text{APC}}(i, j)$, where higher scores identify more likely residue-residue contacts.

Method details
==============

In this section we give more details on the actual implementation of the algorithm described above.

Sequence weighting {#section:measure}
------------------

Families from the Pfam database contain some degree of redundancy. A common strategy to overcome this problem is sequence weighting, which weighs down groups of similar sequences and assigns higher weights to isolated sequences.

We first define a similarity measure between proteins, following [@smale2013introduction]. We start from the BLOSUM90 frequency substitution matrix $B_{90}(a,b)$ defined in [@henikoff1992amino] and call $\hat{B}_{90}(a,b)$ for a pair of amino acids $(a,b)$ the normalized matrix $$\hat{B}_{90}(a,b) = \frac{B_{90}(a,b)}{\sqrt{B_{90}(a,a) B_{90}(b,b) }}$$

We then proceed to construct a similarity kernel between pairs of proteins $$K^3(p, q) = \sum_{k=1}^{10} \left(\sum_{i=1}^{L-k+1} K^2\left((p_i\dots p_{i+k-1}), (q_i\dots q_{i+k-1})\right)\right)$$ where $$K^2\left((p_i\dots p_{i+k-1}), (q_i\dots q_{i+k-1})\right) = \prod_{j=1}^k  \hat{B}_{90}(p_{i+j-1}, q_{i+j-1}),$$ for $p, q \in \PP$, $1 \le k \le L$ and $1 \le i \le L-k+1$;

The normalized version of $K^3$ is obtained using $$\hat K^3(p, q) = \frac{K^3(p, q)}{\sqrt{K^3(p, p)K^3(q,q)}}.$$

Note that, since Pfam families consist of pre-aligned sequences, our $K^3$ kernel definition differs slightly from [@smale2013introduction] as it only compares aligned amino acid $k$-mers. Also, we limit the $k$-mers considered in the construction of $K^3$ to lengths smaller or equal to $10$. This implies a substantial improvement in computation time, with no significant loss in predictive power.

We fix a threshold $\theta$ (in this paper, $\theta = 0.7$) and for any protein $p \in \PP$ we define the equilibrium measure $$\pi(p) = \sum_{\substack{q \in \PP \\ \hat K^3(p, q) > (1-\theta)}} \hat K^3(p, q).$$

The weight of a protein $p$ is then defined as the reciprocal of the equilibrium measure $w(p) = (\pi(p))^{-1}$ , and the effective number of sequences in the alignment is $M_{\text{eff}} = \sum_p w(p)$. The kernel $ \hat K^3(p, q)$ is a measure of similarity between pairs of proteins, and for a given protein $p$ the quantity $\pi(p)$ effectively counts the number of sequences with similarity larger than a threshold $1 -\theta$, thereby weighing down sequences that are over-represented in the data set.

#### Hamming distance weighting

As a term of comparison, we also applied a more traditional sequence weighting scheme based on the hamming distance between pairs of sequences. For each sequence $p$, we count the number of other sequences in the alignment that share more than $\theta$% sequence identity with $p$ $$m^p = | \{ b \in \{ 1, \dots, M \} : \text{\%id}(p,b) > \theta \} |$$ and then assign a weight $w(p) = 1 / m^p$ to sequence $p$.

This approach was used previously by several authors such as [@morcos2011direct], [@jones2012psicov] and [@ekeberg2013improved], but we note that these authors use different values for the theshold $\theta$ (respectively 0.8, 0.62 and 0.9). In this paper we choose $\theta=0.62$. See See Tables S2-S3 for the optimization of the two weighting schemes with respect to the regularization parameter $\eta'$, and Table S4 for a comparison of their performance.

Regularization parameter {#subsec:regularization}
------------------------

The matrix inversion in Equation \[eq:theta\] contains a regularization parameter $\eta$. We observed that families containing few sequences, where the number of sequences $M$ is comparable in size to the number of random variables ($20L$) require a larger regularization parameter compared to bigger families ($M \gg 20L$). We use then a regularization parameter of the form $\eta = \eta' / {M_{\text{eff}}}$, where ${M_{\text{eff}}}$ is the effective number of sequences defined above.

We tried different values of $\eta'$ over the 15 families from [@marks2011protein], and observed that in general roughly the same $\eta'$ appears to be optimal across families with different ${M_{\text{eff}}}$. Thus the normalization $\eta = \eta' / {M_{\text{eff}}}$ appears appropriate.

In Figure we show how the actual eigenvalues of the modified covariance matrix corresponding to PFAM family PF00028 are transformed when computing the regularized inverse.

Pseudocounts {#subsec:pseudocounts}
------------

The pseudocounts parameter $\lambda$ in Equation \[eq:pseudo\] accounts for non-observed pairs of amino acids. Following [@morcos2011direct], we set $\lambda=44$. However, we observed that the performance gain is small compared to setting $\lambda=0$ (see Table S5). In fact, pseudocounts have a similar regularizing effect to the parameter $\eta'$ described in the previous section, and probably for this reason the contribution of $\lambda$ is minimal.

Results and Conclusion
======================

The method and estimation of parameters described above were first applied to the 15 families studied in Marks et al. [@marks2011protein]. Performance was estimated in terms of the fraction of correct predicted contacts among the $L/5$, $L/3$, $L/2$ and $L$ pairs with highest $P_{\text{APC}}$ score, where $L$ is the length of the alignment. We considered as a true contact a pair of amino acids with beta-carbons (C$\beta$) with distance $< 8$ [Å]{} and at least 5 residues apart along the length of the protein. We find that on these 15 families the optimal value for the regularization parameter is around $\eta'=1000$ (see Table S2).

Table \[table:results15\] compares the performance of the RLS algorithm with PSICOV version 1.11 [@jones2012psicov] and the plmDCA method [@ekeberg2013improved]. We observe that on this set our method outperforms both PSICOV and DCA on the majority of families. Additionally, in Table S6 in shown the positive predictive value of the methods with respect to short range ($5 \leq i - j \leq 11$), medium range ($12 \leq i - j \leq 23$) and long range ($>$ 23) interactions for the prediction of the top $L/5$ contacts.

Next, we applied the three methods using the same parameters on an additional set of families from [@ekeberg2013improved]. This set partially overlapped with the families from [@marks2011protein] studied in Table \[table:results15\], and after removing the duplicates we are left with a set of 22 families. Table \[table:results22\] shows the positive predictive value of the RLS algorithm compared to PSICOV and plmDCA. On this set plmDCA obtains the highest average performance on three out of four ranking categories. However, note that plmDCA was optimized on this set of families, so there may be a bias in favor of this method.

Finally, we constructed an independent set of 10 families, selected randomly from Pfam with the only condition of containing at least 1,000 unique sequences. This set had not been used in the optimization of the algorithms, therefore constitutes a fair ground for comparison. The results (see Table \[table:results10\]) show that RLS outperforms the other two methods in the $L/5$ and $L/3$ subsets, whereas DCA obtains highest average performance for the prediction of the top $L/2$ and $L$ contacts.

A comparison of the running times of the two best algorithms (RLS and DCA) shows that RLS is at least one order of magnitude faster than DCA (Table \[table:timings\]). Note that we used the latest fast version of plmDCA [@ekeberg2014fast], termed “asymmetric plmDCA”, which improves considerably on previous pseudolikelihood methods in terms of speed. Our fast regularized inversion of the covariance matrix allows contact prediction on hundreds of amino acids-long domains in a matter of seconds, practically removing the limitations on the length of proteins that can be analyzed. In fact, the slowest step in the predictions is sequence weighting (not accounted in Table \[table:timings\] for either method), in particular the $K^3$-based weighting can be slow for very large families, but we showed that a simpler and faster weighting strategy does not affect too dramatically the performance (Table S4).

In general, we observed that the performance depends on the effective number of sequences ${M_{\text{eff}}}$ in the alignment. For instance, families PF00390 or PF00793 are composed of several thousand sequences, but they contain much redundancy, which brings down ${M_{\text{eff}}}$ to a few hundred units. Roughly, it appears that at least 1000 non-redundant sequences (${M_{\text{eff}}}> 1000$) are necessary to achieve a reasonable precision for contact prediction. This is in agreement with previous estimates [@marks2011protein; @kamisetty2013assessing] which place this number at about $5L$, where $L$ is the length of the alignment.

In conclusion, we demonstrated how our simple regularization scheme for covariance matrix inversion allows the fast and accurate prediction of residue-residue contacts. Currently, a major restriction to this kind of approach is the fairly high number of non-redundant sequences required to infer coevolution from a multiple sequence alignment, limiting the application to a relatively small subset of Pfam. However, as the number of protein sequences deposited in public databases increases, we expect a larger number of protein families to become accessible to our analysis, as well as improved performance on those that are already accessible.

Acknowledgement {#acknowledgement .unnumbered}
===============

This work received funding from City University of Hong Kong grants RGC \#9380050 and \#9041544. Work by S.L. was partially supported by Ministerio de Ciencia, Tecnología e Innovación Productiva, Argentina.
