Selection of sequence motifs and generative Hopfield-Potts models for protein families

Statistical models for families of evolutionary related proteins have recently gained interest: In particular, pairwise Potts models as those inferred by the direct-coupling analysis have been able to extract information about the three-dimensional structure of folded proteins and about the effect of amino acid substitutions in proteins. These models are typically requested to reproduce the one- and two-point statistics of the amino acid usage in a protein family, i.e., to capture the so-called residue conservation and covariation statistics of proteins of common evolutionary origin. Pairwise Potts models are the maximum-entropy models achieving this. Although being successful, these models depend on huge numbers of ad hoc introduced parameters, which have to be estimated from finite amounts of data and whose biophysical interpretation remains unclear. Here, we propose an approach to parameter reduction, which is based on selecting collective sequence motifs. It naturally leads to the formulation of statistical sequence models in terms of Hopfield-Potts models. These models can be accurately inferred using a mapping to restricted Boltzmann machines and persistent contrastive divergence. We show that, when applied to protein data, even 20–40 patterns are sufficient to obtain statistically close-to-generative models. The Hopfield patterns form interpretable sequence motifs and may be used to clusterize amino acid sequences into functional subfamilies. However, the distributed collective nature of these motifs intrinsically limits the ability of Hopfield-Potts models in predicting contact maps, showing the necessity of developing models going beyond the Hopfield-Potts models discussed here.

Figure 1

Panel (a) represents the (Hopfield-)Potts model as a statistical model for sequences $\underline{A}\in {\mathcal{A}}^L$, typically characterized by a fully connected coupling matrix $J$ and local fields $h$ (not represented). The model can be transformed into a RBM by introducing Gaussian hidden variables $\underline{x}\in {\mathbb{R}}^p$ with $p$ being the rank of $J$. Note the bipartite graphical structure of RBM, which causes the conditional probabilities $P\left(\underline{A}\right|\underline{x})$ and $P\left(\underline{x}\right|\underline{A})$ to factorize. Panel (b) shows a schematic of persistent contrastive divergence (PCD). Initially, the sample is initialized in the training data (the MSA of natural sequences), and then, $k$ alternating steps of sampling from $P\left(\underline{A}\right|\underline{x})$, respectively, $P\left(\underline{x}\right|\underline{A})$'s are performed. Parameters are updated after these $k$ sampling steps, and sampling is continued using the updated parameters.

Figure 2

Statistics of natural sequences (PF00072, horizontal axes) vs MCMC samples (vertical axes) of Hopfield-Potts models for values of $p\in \{5,10,20,40,80,160\}$ and for a full-rank Potts model inferred using BMDCA. The first column [panels (a.1)–(g.1)] shows the one-point frequencies $f_i\left(a\right)$ for all pairs $(i,a)$ of sites and amino acids; the other two columns show the connected two- and three-point functions $c_{ij}(a,b)$ [panels (a.2)–(g.2)] and $c_{ijk}(a,b,c)$ [panels (a.3)–(g.3)]. Due to the huge number of combinations for the three-point correlations, only the 100 000 largest values (evaluated in the training MSA) are shown. The Pearson correlations and the slope of the best linear fit are inserted in each of the panels.

Figure 3

Panel (a) shows the positive predictive value (PPV) for contact prediction as a function of the number of predictions for various values of the pattern number $p$ and for BMDCA. Panels (b)–(d) show, for $p=20,160$, and BMDCA, the distribution of coupling scores $F_{ij}^{\mathrm{APC}}$. All residue pairs are grouped into contacts (red) and noncontacts (blue). The best contact predictions correspond to the positive tail of the red histogram, which becomes more pronounced when increasing $p$ or even going to BMDCA.

Figure 4

Likelihood contribution of the individual patterns for pattern numbers $p=20,40,80$, and 160.

Figure 5

The five patterns of highest likelihood contribution for $p=20$, the eighth ranking is added since used later in the text. The left panels (a.1)–(f.1) show the patterns in logo representation, the letter size is given by the corresponding element ${\xi }_i\left(a\right)$. The middle panels (a.1)–(f.2) show the distribution of the activities, i.e., the projections of sequences onto the patterns. The blue histogram contains the natural sequences from the training MSA, and the red histogram contains sequences sampled by MCMC from the Hopfield-Potts model. The right-hand side [panels (a.3)–(f.3)] shows the connected two-point correlations of the natural data (horizontal axis) vs data sampled from $P_{-\mu }\left(\underline{A}\right)$, i.e., a Hopfield-Potts model with one pattern removed. Strong deviations from the diagonal are evident.

Figure 6

Patterns with multimodal activity distributions for the set of all MSA sequences can be used to cluster sequences. The rows show combinations of patterns 6–13 [panels (a.1)–(f.1)], 6–14 [panels (a.2)–(f.2)], and 13 to 14 [panels (a.3)–(f.3)]. Each sequence corresponds to a density-colored dot. A strongly clustered structure is clearly visible. When dividing the full MSA into functional subclasses, we can relate clusters to subclasses and, thus, patterns to biological function.

Figure 7

The same as Fig. 2 but for the protein family PF00014.

Figure 8

The same as Fig. 3 but for the protein family PF00014.

Figure 9

The same as Fig. 5 but for the protein family PF00014.

Figure 10

The same as Fig. 2 but for the protein family PF00076.

Figure 11

The same as Fig. 3 but for the protein family PF00076.

Figure 12

The same as Fig. 5 but for the protein family PF00076.

Figure 13

The upper two panels show the statistics (two-point connected correlations) of the training data (PF00014) against an i.i.d. MCMC sample extracted from the inferred model and the CD sample used for inference. The perfect coincidence of the two in the CD case demonstrates that the CD algorithm is converged and contrast is lost despite the fact that the correlations extracted from an i.i.d. sample do not match the empirical ones. To understand this, we have selected those $(i,j,a,b)$'s with substantial deviations, cf. the inset in the first panel and analyzed their location in the protein (first panel) and their amino acid composition (second panel, amino acids in alphabetical order of one-letter code $[-,A,C,...,Y]$), densities are represented via heat map plots. Location at the extremities in the sequences and in gap-gap correlations emerge clearly.

Figure 14

The upper two panels show the statistics (two-point connected correlations) of the training data (PF00014) against an i.i.d. MCMC sample extracted from the inferred model and the PCD sample used for inference. As in the CD case, the MCMC sample shows larger deviations from the empirical observations than the PCD sample, but correlations appear overestimated and contrast in the PCD plot is sufficient to drive further evolution of parameters. To understand these observations, we have selected those $(i,j,a,b)$'s with substantial deviations, cf. the inset in the first panel and analyzed their location in the protein (first panel) and their amino acid composition (second panel, amino acids in alphabetical order of one-letter code $[-,A,C,...,Y]$), densities are represented via heat map plots. Locations at the extremities in the sequences and in gap-gap correlations emerge clearly.

Figure 15

Regularization dependence for CD inference, empirical two-point connected correlations (PF00014) are plotted against those estimated from the model using an i.i.d. MCMC sample. The regularization strength is varied over almost two orders of magnitude with ${\eta }_0={\alpha }_{0}L/qM$ and ${\alpha }_0=0.0525$, going from a zone of overfitting to one of over-regularization. Results are shown for various values of $p$, illustrating a strong $p$ dependence of the optimal coupling strength.

Figure 16

Regularization dependence for PCD inference, empirical two-point connected correlations (PF00014) are plotted against those estimated from the model using an i.i.d. MCMC sample. The regularization strength is varied over almost two orders of magnitude with ${\eta }_0={\alpha }_{0}L/qM$ and ${\alpha }_0=0.0525$; results are shown for various values of $p$. We find a strong robustness of results with respect to regularization.