Generative power of a protein language model trained on multiple sequence alignments

Abstract

Computational models starting from large ensembles of evolutionarily related protein sequences capture a representation of protein families and learn constraints associated to protein structure and function. They thus open the possibility for generating novel sequences belonging to protein families. Protein language models trained on multiple sequence alignments, such as MSA Transformer, are highly attractive candidates to this end. We propose and test an iterative method that directly employs the masked language modeling objective to generate sequences using MSA Transformer. We demonstrate that the resulting sequences score as well as natural sequences, for homology, coevolution, and structure-based measures. For large protein families, our synthetic sequences have similar or better properties compared to sequences generated by Potts models, including experimentally validated ones. Moreover, for small protein families, our generation method based on MSA Transformer outperforms Potts models. Our method also more accurately reproduces the higher-order statistics and the distribution of sequences in sequence space of natural data than Potts models. MSA Transformer is thus a strong candidate for protein sequence generation and protein design.

Keywords: computational biology; deep learning; none; protein design; protein families; protein language models; protein sequence generation; protein sequences; systems biology.

Figure 1.. Comparison of homology, coevolution, and structure-based scores between natural sequences and sequences generated by MSA Transformer or Boltzmann machine DCA (bmDCA).

For each Pfam family in Appendix 1—table 5, we compare a natural MSA from Pfam and three synthetic MSAs of the same depth. The first synthetic MSA was obtained using MSA Transformer via our iterative masking procedure, and the second and third ones were generated by a Potts model inferred from the natural MSA using bmDCA with two different pairs $\mathrm{(}\lambda \mathrm{,}T\mathrm{)}$ of regularization strength $\lambda$ and sampling temperature $T$. For each of the four scores described in ‘Scoring individual sequences’, we show the distributions of score values among sequences in each MSA as a violin plot. Higher score values are better for all scores except root-mean-squared deviation (RMSD) (bottom panel), where smaller values indicate a closer match to an experimental structure. Top panel: For each Pfam family, HMMER scores are divided by the highest score found in the natural MSA. Note that sequences below HMMER’s default homology detection score (E-value larger than 10), and whose HMMER score is thus 0, are not shown (the median over families of the fraction of such sequences is 2% for bmDCA (10⁻², 1.00)-generated MSAs, while there are no such sequences among the MSA-Transformer–generated ones). Second panel: Statistical energy scores are defined as minus the bmDCA statistical energies. To accommodate the highly family-dependent ranges of these scores, for each Pfam family we show their values after shifting by the mean score in the natural MSA, and normalizing by the standard deviation of natural MSA scores. Third panel: AlphaFold’s predicted local-distance difference test (pLDDT) confidence scores. Bottom panel: RMSD of predicted structures with respect to the experimental structures in Appendix 1—table 5. Structural scores (pLDDT and RMSD) were computed on 200 randomly chosen sequences from each MSA. All kernel-smoothed histograms are normalized such that all violins have the same maximal width. Outliers (less than 1% in all cases) were discarded for legibility.

Figure 1—figure supplement 1.. Multiple sequence alignment (MSA) diversity for each protein family and each generation method.

We show the relative effective depth of each MSA (natural or generated), that is $M_{\mathrm{\text{eff}}}\mathrm{/}M$, where $M_{\mathrm{\text{eff}}}$ is effective MSA depth and $M$ is actual MSA depth, versus the similarity threshold $\mathrm{1}\mathrm{-}\delta$ used to define the effective depth (see Equation 8).

Figure 2.. Homology and coevolution scores versus distance to the natural multiple sequence alignment (MSA), for protein families PF00072 and PF00153.

We show contour plots of the HMMER score and the statistical energy score (defined as minus the DCA statistical energy, shifted by its mean value in the natural MSA) versus the Hamming distance of each sequence to the closest natural sequence (which is not itself, in the case of natural sequences). Results are shown for natural sequences and for sequences generated using MSA Transformer and Boltzmann machine DCA (bmDCA) (the same two $\mathrm{(}\lambda \mathrm{,}T\mathrm{)}$ pairs as in Figure 1 are used for bmDCA). The lightest contours shown include 99% of the cumulative probability mass.

Figure 3.. Application of our sequence generation method based on MSA Transformer to small protein families.

We consider seven small protein families, with natural MSAs that comprise from nine to a few hundreds of sequences, see Appendix 1—table 6. As in Figure 1, for each family, we compare the natural MSA and three synthetic MSAs of the same depth. In all cases, we show violin plots of the same four scores as for large families in Figure 1, as well as of the Hamming distance to the closest natural sequence, which is not itself in the case of natural sequences (‘Distance’). For the three smallest families (left panel; fewer than 40 sequences), we also show the score of each individual sequence as a swarm plot. Note that while we employ the same sampling temperatures $T$ as in Figure 1 for Boltzmann machine DCA (bmDCA), here, we use regularization strength $\lambda ={10}^{-2}$ throughout, due to MSA shallowness (see ‘Sampling sequences from Potts models’).

Figure 4.. Similarity of statistics between synthetic and natural multiple sequence alignments (MSAs).

To compare the statistics of synthetic and natural MSAs at various orders, we compute r20 scores (Haldane et al., 2018; McGee et al., 2021), and plot them versus the number of different MSA columns that are considered (see ‘Analyzing the statistics of MSAs’ for details). All families in Figure 5 are considered. For each of them, the reference MSA comprises either half of the natural MSA (with sequences selected uniformly at random), or 30,000 sequences from it if the natural MSA depth is larger than 60,000. The null model compares the other half of the natural MSA to this reference MSA. It yields an estimate of the expected r20 scores due only to finite-size effects in a model-free, purely data-driven way.

Figure 4—figure supplement 1.. Ability of generated sequences to reproduce one-, two-, and three-body statistics.

In all panels, each marker represents a family in 5. For each statistical or information measure (rows), different scores (columns) comparing generated and natural sequences are the coordinates of these points. All scores are such that being close to 0 is better, and their value for multiple sequence alignments (MSAs) generated by Boltzmann machine DCA (bmDCA) with default parameters is shown versus that for MSAs generated by MSA Transformer. The statistical or information measures considered in each row are defined in ‘Analyzing the statistics of MSAs’ – from the top: one-body frequency, two- and three-body connected correlations, entropy, mutual information, and co-information. For each of them, we consider its values over all MSA columns (or pairs or triplets of columns), and all amino acids if appropriate, for both natural and synthetic MSAs. To obtain the vertical and horizontal coordinates (respectively) of the markers in each panel, we compare these values for each natural MSA with the values from the corresponding synthetic MSAs generated by bmDCA with default parameters or by our method based on MSA Transformer (respectively). We use four different scores for this comparison, and devote each column of the figure to one of these scores – from the left: $|\rho -1|$ where $\rho$ denotes the Pearson correlation; $|\text{Slope}-1|$ where ‘Slope’ means the slope of best linear fit (see Figure 4—figure supplement 2 and Figure 4—figure supplement 3 for illustrations of these first two quantities in the case of two- and three-body connected correlations for families PF00072 and PF00153); the Jensen–Shannon divergence between the distributions of values; the Wasserstein distance between these distributions. For each statistical or information measure (row) and each score (column), and for each family in 5, we have one value of the score comparing the natural and bmDCA-generated MSAs and another one comparing the natural and MSA-Transformer–generated MSAs. We plot the former value versus the latter, yielding one marker per protein family in each plot. Thus, each plot compares the ability of bmDCA and MSA Transformer to reproduce the statistics of the natural data. Blue markers (above the diagonal) mean that the scores for MSA-Transformer–generated MSAs are better, while green markers (below the diagonal) mean the opposite.

Figure 4—figure supplement 2.. Two- and three-body connected correlations estimated from generated multiple sequence alignments (MSAs) versus the natural one, for family PF00072.

Relationships between connected correlations estimated from the MSA generated by MSA Transformer or Boltzmann machine DCA (bmDCA), and those estimated from the natural MSA, are shown as binned scatter plots both for two-body (top row) and three-body (bottom row) statistics. We include a null model (third column) obtained by splitting the natural MSA in half and comparing the statistics of one half with those of the other. It yields an estimate of the expected dispersion in these plots due only to finite-size effects in a model-free, purely data-driven way. Pearson correlation coefficients, $\rho$ and slopes of lines of best fit, are reported in each case. For the comparisons involving MSA Transformer, we used Equations 4 and 5 to estimate correlations. On the other hand, since bmDCA is trained to reproduce frequencies rescaled with phylogenetic weights (w_i in 8), for the comparisons involving bmDCA we rescaled the natural frequencies before using Equations 4 and 5 to estimate correlations.

Figure 4—figure supplement 3.. Two- and three-body connected correlations estimated from generated multiple sequence alignments (MSAs) versus the natural one, for family PF00153.

Same as Figure 4—figure supplement 2, but for family PF00153.

Figure 5.. Distribution of sequences in sequence space, for families PF00072 and PF00153.

We show the distribution of one-hot encoded natural and synthetic sequences projected in the subspace of the first two principal components of the natural multiple sequence alignment (MSA). The same axis limits are used within one family, except for Boltzmann machine DCA (bmDCA) (10⁻³, 0.33) in the case of PF00072. Note that the fraction of the total variance explained by the first two principal components of each MSA is less than 4% for all families and all generation methods.

Figure 5—figure supplement 1.. Distribution of sequences in sequence space for all large protein families in our dataset (part 1).

We show the distribution of one-hot encoded natural and synthetic sequences projected in the subspace of the first two principal components of the natural multiple sequence alignment (MSA), as in 5 for families PF00072 and PF00153, but also including Boltzmann machine DCA (bmDCA) at $\mathrm{(}\lambda \mathrm{,}T\mathrm{)}\mathrm{=}\mathrm{(}{\mathrm{10}}^{\mathrm{-}\mathrm{3}}\mathrm{,}\mathrm{0.66}\mathrm{)}$.

Figure 5—figure supplement 2.. Distribution of sequences in sequence space for all large protein families in our dataset (part 2).

We show the distribution of one-hot encoded natural and synthetic sequences projected in the subspace of the first two principal components of the natural multiple sequence alignment (MSA), as in 5 for families PF00072 and PF00153, but also including Boltzmann machine DCA (bmDCA) at $\mathrm{(}\lambda \mathrm{,}T\mathrm{)}\mathrm{=}\mathrm{(}{\mathrm{10}}^{\mathrm{-}\mathrm{3}}\mathrm{,}\mathrm{0.66}\mathrm{)}$.

Figure 5—figure supplement 3.. Neighbors of natural and synthetic sequences, for families PF00072 and PF00153.

We show the distribution of the number of neighbors of sequences in the natural multiple sequence alignment (MSA), and the distribution of the number of neighbors of the closest natural sequence to each of our generated sequences. Given a sequence in a natural MSA, its number of neighbors is the number of natural sequences that are within a (normalized) Hamming distance $\delta =0.2$ from it. The moving average of the results is shown, using a window representing 5% of the total number of points.

Figure 5—figure supplement 4.. Comparing phylogenies inferred from natural and generated multiple sequence alignments (MSAs) for families PF00072 and PF00153.

We show the averaged spectra of modified graph Laplacian (MGL) matrices computed, for each tree inferred from an MSA, by using the leaves of multiple sub-trees made of 500 randomly sampled sequences from the MSA of interest. Specifically, in each case, we perform an average over 200 different sub-trees of the histograms of counts of the eigenvalues (see ‘Characterizing the distribution of sequences in MSAs’). We compare MSAs generated using either MSA Transformer (left) or Boltzmann machine DCA (bmDCA) (center and right) to the natural MSA.

Figure 6.. Comparison of our generated sequences to those experimentally tested in Russ et al., 2020, for the chorismate mutase family.

Left: The estimated relative enrichment (r.e.) scores of the Boltzmann machine DCA (bmDCA)-generated sequences that are in the top 33% in terms of predicted local-distance difference test (pLDDT) scores are plotted versus their experimentally measured counterparts from Russ et al., 2020. We estimate the expected r.e. of these generated sequences as the r.e. of the closest natural sequence measured in Russ et al., 2020. We observe that high estimated r.e. is associated with high measured r.e., as 71% of sequences with estimated r.e. > 0.4 (green) also have measured r.e. > 0.4. Note that in the top marginals (showing the measured r.e. for bmDCA-generated sequences), the green and yellow histograms are stacked on top of each other. Thus, the stacked histogram shows the distribution of all measured r.e. values for bmDCA-generated sequences that are in the top 33% in terms of pLDDT scores. Top right: Overlaid histograms of estimated r.e. are shown for our MSA-Transformer–generated sequences and for the bmDCA-generated ones from Russ et al., 2020, restricting in both cases to the sequences with top 33% pLDDT scores. Bottom right: Same as top right, but considering all generated sequences.

Figure 6—figure supplement 1.. Homology, coevolution, and structural scores versus distance to the natural multiple sequence alignment (MSA), for the chorismate mutase family.

We assess the performance of our generative method in the case of the chorismate mutase family, by comparing our generated sequences (‘MSA Tr.’) to natural ones, and to the sequences generated in Russ et al., 2020 using Boltzmann machine DCA (bmDCA) at various values of the sampling temperature and of the regularization strength. Our synthetic sequences were generated using the iterative masking procedure based on MSA Transformer, starting from the natural alignment used in Russ et al., 2020. Our sequences score similar to natural sequences and to the bmDCA-generated sequences from Russ et al., 2020, which were experimentally validated there.

Figure 6—figure supplement 2.. Deep mutational scanning (DMS) scores for families PF00595 and PF13354.

In DMS experiments, the fitness effects of all possible one-point mutations from a reference natural sequence are measured. To assign a DMS score to natural and generated sequences, we align each of them to the reference sequence, and sum the experimentally measured fitness effects of the relevant amino acids at each position. Higher values of the DMS score are better, as they mean higher fitness. We show normalized histograms of the DMS scores for sequences in the natural and generated multiple sequence alignments (MSAs) of protein families PF00595 and PF13354, based on the DMS experiments in McLaughlin et al., 2012; Stiffler et al., 2015, respectively. For generated sequences, we restrict to those whose Hamming distance to their closest natural neighbor is larger than $\delta \mathrm{=}\mathrm{0.2}$.

Figure 7.. Iterative masking procedure to generate sequences using MSA Transformer.

Here, the red hashtag (#) stands for a masked amino acid, while blue uppercase letters stand for predicted amino acids at the masked positions.

Figure 7—figure supplement 1.. Evolution of mean scores during the iterative masking procedure, for family PF00153.

Average scores of the generated sequences are reported for different iteration numbers and masking probabilities. The scores employed are: (A) Hamming distances to the closest natural sequence, (B) Mutual information between synthetic and natural columns of the MSAs, (C) HMMER scores – see ‘Scoring individual sequences’, (D) Statistical energy scores (negative DCA statistical energies, shifted by their mean value for natural sequences – see ‘Scoring individual sequences’), (E) Pairwise Hamming distances between the generated sequences, and (F) Mean square deviations (MSD) of the predicted contact maps at each iteration from the initial one (zero iterations). The synthetic MSAs, comprising 5000 sequences, were generated with MSA Transformer starting from natural sequences of family PF00153.

Figure 7—figure supplement 2.. Evolution of mean scores during the iterative masking procedure, for family PF00096.

Same as in Figure 7—figure supplement 1, but for PF00096. This family has the shortest length among those considered here, namely $L\mathrm{=}\mathrm{23}$ (see Appendix 1—table 5).

Figure 7—figure supplement 3.. Evolution of mean scores during the iterative masking procedure, for family PF13354.

Same as in Figure 7—figure supplement 1, but for PF13354. This family has the largest length among those considered here, namely $L\mathrm{=}\mathrm{198}$ (see Appendix 1—table 6).

Figure 7—figure supplement 4.. Evolution of inferred contact maps during the iterative masking procedure, for family PF00153.

Contact maps obtained from MSA Transformer for different iteration numbers and masking probabilities are reported. Probabilities of contacts are computed using the logistic regression on the output attention matrices from Rao et al., 2021a. The input of the model consists of 100 different sequences chosen uniformly at random from the synthetic MSA generated at each iteration of our iterative masking procedure. These contact maps were employed to compute the mean square deviations shown in Figure 7—figure supplement 1F.