Amino-acid site variability among natural and designed proteins

Data sets

We analyzed two data sets, one of whole yeast proteins and one of protein domains. The yeast-proteins data set was taken from Ramsey et al. (2011) and comprised 38 protein structures homologous to an open reading frame in Saccharomyces cerevisiae. For each of those structures, we had at least 50 homologous natural sequences, also taken from Ramsey et al. (2011). The protein-domain data set was taken from Ollikainen & Kortemme (2013) and comprised 40 protein domains. Only domains with at least one crystal structure in the Protein Database (PDB) and at least 500 sequences in the Pfam Database were selected for this data set. Domains were selected in order to represent several different types of protein folds and domains were also restricted to a length less than or equal to 150 amino acids. For each of these protein domains, we obtained alignments of homologous natural sequences from the Pfam database (Punta et al., 2012), as described (Ollikainen & Kortemme, 2013).

Protein design

For each structure in both data sets, we computationally designed 500 variants each, using multiple design methods. All design methods we used are implemented in the protein-design software Rosetta (Leaver-Fay et al., 2011). First, we used standard fixed-backbone design (Kuhlman et al., 2003). In this method, the protein backbone remains fixed and only amino-acid side chains are allowed to move. Second, we used the flexible-backbone method Backrub (Smith & Kortemme, 2008), which first generates an ensemble of alternative backbones and then designs side chains onto these backbones (Friedland et al., 2009; Smith & Kortemme, 2010). The Backrub method takes as input a temperature parameter that determines the extent of backbone movements that occur during design. Here, we used temperatures spanning from 0.03 to 2.4 corresponding to increasing backbone movements. For the protein-domain data set, we also used one additional design method, called “Soft”. This method keeps the backbone fixed but the energy function used during sequence design dampens the weight of the repulsive Lennard-Jones (LJ) potential term (Ollikainen & Kortemme, 2013). Protein designs for the protein-domain data set have been previously published (Ollikainen & Kortemme, 2013), while the designs for the yeast-proteins data set were newly generated for the present study.

All designs for the yeast-proteins data set were generated with Rosetta Revision 39284. For fixed-backbone design, we used the following command:

./fixbb.linuxgccrelease -database rosetta_database \
    -s input.pdb -resfile ALLAA.res -ex1 -ex2 \
    -extrachi_cutoff 0 -nstruct 1 -linmem_ig 10

Flexible-backbone design was performed by generating a conformational ensemble of 500 structures and then using fixed-backbone design to predict a low energy sequence for each structure in the ensemble. To generate structures for the conformational ensemble, we used the following command:

./backrub.linuxgccrelease -database rosetta_database \
    -s input.pdb -resfile NATAA.res -ex1 -ex2 \
    -extrachi_cutoff 0 -backrub:mc_kt <T> \
    -backrub:ntrials 10000 -nstruct 1 -backrub:initial_pack

where <T> has to be replaced by the desired design temperature.

The design details for the protein-domain data set can be found in Ollikainen & Kortemme (2013).

Data analysis

We quantified the variability of sites in amino-acid alignments using site entropy H_i, defined as $H_i=\sum _j{p}_{ij}ln{p}_{ij}$ . Here, p_ij is frequency of amino acid j in alignment column i, and the sum runs over all amino acids. We compared amino-acid distributions of designed sequences to those of natural sequences using the Kullback-Leibler (KL) divergence. The KL divergence $D_i^{\text{KL}}$ is defined as $D_i^{\text{KL}}=\sum _j{p}_{ij}ln(p_{ij}/q_{ij})$ , where q_ij is the frequency of amino acid j in column i of the reference alignment, and p_ij is the corresponding frequency in the alignment that is being compared to the reference alignment. The sum runs over all amino acids. When calculating frequencies used for the KL divergence we corrected for the presence of frequencies of zero by adding 1/20 to each amino acid count before calculating the frequencies. The KL divergence is inherently an asymmetric distance measure, comparing a probability distribution of interest to a reference distribution. Unless noted otherwise, we always used natural sequence alignments to calculate the reference frequencies q_ij and designed sequence alignments to calculate the frequencies p_ij. Throughout this work, we calculated $D_i^{\text{KL}}$ separately at each site i in a protein, and then averaged the $D_i^{\text{KL}}$ values for all sites in a protein to obtain a mean KL divergence for that protein.

To compare the shapes of amino-acid distributions while disregarding specific amino-acid identities, we performed a second type of KL calculation where we ordered amino-acids by their relative frequencies. Thus, instead of the frequencies p_ij and q_ij we used p_irj and q_isj, where r_j is the rank of the frequency of amino acid j in column i of the alignment being compared to the reference, and s_j is the rank of the frequency of amino acid j in column i of the reference alignment. This way of calculating the KL divergence compares the frequencies of amino acids at equal frequency rank, regardless of which specific amino acids are the most frequent, second-most frequent, and so on in each alignment. As an example, assume that at a given site there are only three different amino acids in the natural alignment, I, L, and V, at frequencies 0.5, 0.35, and 0.15, respectively. At the same site in the designed sequences, there are amino acids A, V, and I, also at frequencies 0.5, 0.35, and 0.15, respectively. In our calculation of KL divergence comparing amino acids at equal frequency rank, we would then compare the frequency of I in the natural alignment with the frequency of A in the designed alignment (the two most frequent amino acids in the two respective alignments) and similarly the frequency of L with the frequency of V and the frequency of V with the frequency of I, respectively. In this example, since the two sets of three frequencies are exactly the same if we disregard amino-acid type, we would obtain a KL divergence of zero.

We calculated Relative Solvent Accessibility (RSA) of residues by first calculating the absolute Solvent Accessibility (ASA) for each residue, using the software DSSP (Kabsch & Sander, 1983). For each protein, we extracted the chain of interest from the PDB structure and ran DSSP only on that chain. We calculated RSA by dividing the ASA value for each residue by the maximum possible ASA value, as given by Tien et al. (2013). Throughout this work, we only calculated RSA on the native PDB structure. We did not perform any RSA calculations on designed structures.

All our data and analysis scripts are available online at: https://github.com/clauswilke/protein_design_and_site_variability.

Overall site variability

We first compared overall amino-acid variability in designed and natural proteins. We assessed amino-acid variability at individual sites by calculating the entropy H_i at each site i in alignments of either designed or natural proteins. We then calculated the mean entropy over all sites in each alignment and used that quantity as a measure of the overall amino-acid variability in the alignment.

We found that protein design using a fixed backbone generally yielded insufficient site variability compared to natural sequences (Fig. 1). This result was magnified in the smaller protein domains. In fact, for the protein domains, the most variable proteins under fixed-backbone design showed only about as much variability as the least variable natural proteins. Overall, there was a significant shift towards higher variability in natural proteins relative to proteins designed with fixed backbone (paired t test, P = 1.4 × 10⁻¹⁰ for the yeast-proteins data set and P < 10⁻¹⁵ for the protein-domain data set). When switching from fixed-backbone design to variable-backbone design, we found that overall site variability increased. Further, site variability increased monotonously with the degree of backbone flexibility allowed during design, as measured by the design temperature (Fig. 1). At the highest temperatures, site variability in designed proteins consistently exceeded that of natural proteins.

Proteins designed at intermediate temperatures had site variabilities that most closely resembled that of natural proteins. For the yeast-proteins data set, the temperature that provided the closest match was T = 0.03, even though the variability of sequences designed at that temperature still exceeded the variability in natural sequences (paired t test, P = 0.0006). For the protein-domains data set, the temperature that provided the closest match was T = 0.9, for which variability was statistically indistinguishable from that found in natural sequence alignments (paired t test, P = 0.353). However, for both data sets, natural sequences generally showed a larger spread in variabilities than did the designed sequences at the closest-matching temperatures (Brown-Forsythe test for equal variances, P = 0.0003 for the yeast-proteins data set at T = 0.03 and P = 7.3 × 10⁻⁶ for the protein-domain data set at T = 0.9).

Amino-acid distributions

We next compared amino-acid distributions between designed and natural sequences. First we looked at overall amino acid frequencies. We found that by-and-large, amino acid frequencies in designed proteins mirrored those in natural proteins (Fig. 2 and Figs. S1–S5). The biggest differences arose in Pro, His, Trp, Phe, and Ala. (We ignore Cys here because Cys is never used in the design algorithm and thus always at frequency 0.) Overall, we observed that hydrophobic residues tended to be under-represented in designed proteins whereas hydrophilic residues tended to be over-represented. This trend was stronger in the protein core than on the surface. We also observed that the longer proteins in the yeast-proteins data set showed larger deviations between designed and natural sequences than the shorter proteins in the protein-domains data set. Finally, when comparing different design methods and design temperatures, we found that differences in amino-acid distributions were relatively minor, see Fig. 2 and Figs. S1–S5.

Even if overall amino-acid distributions are approximately correct, the amino-acid distributions at individual sites can be poorly predicted (Ramsey et al., 2011). Therefore, we next compared, separately at each site, the similarity between amino-acid distributions in natural proteins and those in designed proteins. To carry out this comparison, we employed the Kullback-Leibler (KL) divergence (Wasserman, 2004), which measures how similar one probability distribution is to a reference distribution. A KL divergence of zero implies that the distributions are identical. The higher the KL divergence, the more dissimilar the focal distribution is to the reference distribution. (Note that KL divergence is not symmetric: if we swap the focal and the reference distribution, we will generally obtain a different KL divergence value.) We calculated the KL divergence at each site in each protein, and then averaged over sites within a protein to obtain a mean similarity score for each protein. As a control, we also randomly split the alignment of natural sequences for each protein structure into two halves and calculated the mean KL divergence of natural sequences against themselves.

First, in all comparisons, we found that the KL divergence of designed relative to natural sequences was much bigger than the KL divergence of natural sequences relative to themselves (Fig. 3 and Fig. S6). This finding indicates a substantial discrepancy between designed and natural sequences at individual sites. Second, we found that the mean KL divergence decreased with increasing design temperature (Fig. 3A and Fig. S6A). Thus, according to the KL divergence measure, structures designed with the most flexible backbones had the most similar amino-acid distributions to those found in natural sequences.

However, the result that sequences designed at the highest temperatures are the most similar to natural sequences may be an artifact of the KL divergence measure. As design temperature increases, amino-acid variability increases, and amino-acid distributions become more uniform. A more uniform distribution is generally going to display more overlap with any given distribution than a more localized distribution, if the localized distribution is not correct. Thus, the decrease in KL divergence with increasing temperature may simply reflect the broadening of the distribution, not an actual improvement in reproducing natural amino-acid distributions. To assess whether amino-acid distributions in designed sequences were simply broadening with increasing temperature, or whether they were actually converging on the natural distributions, we carried out a second set of comparisons. We rank-ordered amino acids by frequency at each site in each protein, and then calculated the KL divergence of the rank-ordered distributions.

This comparison considers only the shape of the distribution and does not assess whether the correct amino acids are present at individual sites. This second comparison generally found much lower KL divergence levels, even though still not as low as what was found for the control comparison of natural sequences with themselves (Fig. 3B and Fig. S6B). More importantly, now KL divergence reached a minimum around a temperature of 0.3 (yeast proteins, Fig. 3B) to 1.2 (protein domains, Fig. S6B) and rose again beyond that value. This finding indicates that higher design temperatures do not unequivocally produce more natural amino-acid distributions. Instead, there is an intermediate temperature, approximately coinciding with the temperature at which overall sequence variability matches best, at which amino acid distributions also are most similar.

Site variability and solvent accessibility

The previous analyses demonstrated that while designed proteins overall look similar to natural proteins, there are also important differences. We next wanted to identify whether these differences were present uniformly throughout the structure or could be located to specific structural regions. In our analysis of amino-acid distributions, we had already seen that amino-acid distributions seemed to deviate more at buried sites than at exposed sites (Fig. 2 and Fig. S4).

We first plotted site variability against relative solvent accessibility (RSA, a dimensionless number from 0 to 1 measuring the relative solvent exposure of individual residues) for individual proteins. See Fig. S7 for one example. We generally found that site variability displayed a substantial spread even for sites of very similar RSA. At the same time, there was an overall trend for sites with higher RSA to be more variable than sites with lower RSA. This trend was generally stronger in flexible backbone designs than in fixed backbone designs (Fig. S7). To analyze the relationship between site variability and RSA more systematically, we calculated the correlation between these two quantities for all proteins (Fig. 4 and Fig. S8). On average, natural sequence alignments showed a higher correlation than alignments of designed sequences, regardless of design method.

Intermediate design temperatures showed the highest correlations, but correlations were nevertheless significantly lower in designed proteins than in natural proteins (paired t test, P = 2.96 × 10⁻¹⁰ [T = 0.3, yeast proteins] and P = 1.75 × 10⁻⁵ [T = 0.3, protein domains]). We also investigated whether the designed proteins with the highest correlations corresponded to the natural proteins with the highest correlations, and found this generally to be the case (Fig. 4B and Fig. S8B).

Our finding that correlations between site entropy and RSA are lower in designed proteins than in natural proteins indicates that, in designed proteins, site variability is too uniform across different solvent exposure states. In short, designed proteins are either too variable in the core or too conserved on the surface. To obtain a clearer picture of how exactly designed proteins differed from natural proteins, we once more considered the distributions of mean site entropies, but now calculated separately for buried sites (RSA ≤ 0.05), for partially buried sites (0.05 < RSA ≤ 0.25), and for exposed sites (RSA > 0.25). Figure 5 shows the medians of these distributions. For designed proteins, the mean site variabilities of exposed and of partially buried sites are close in magnitude while the mean site variabilities of buried sites are generally consistently lower. By contrast, in natural sequences exposed sites show much more variability than partially buried sites.

If buried sites are too variable or exposed sites too conserved in designed proteins, we reasoned that hybrid designs, in which buried sites were taken from sequences designed at a lower temperature and exposed sites from sequences designed at a higher temperature, should display correlations more similar to those seen in natural proteins.

According to Fig. 5, for the yeast proteins buried and partially buried sites in designed proteins had site variability most similar to that of natural sequences in proteins designed with a fixed backbone or in proteins with a design temperature of T = 0.03. In the protein-domains data set, that temperature was T = 0.3 to T = 0.6. By contrast, for exposed sites the site variability in designed proteins was most similar to that of natural sequences at a design temperature of T = 0.1 (yeast proteins) and T = 1.2 (protein domains). We thus built our hybrid designs by combining sites from these temperatures. We found that the distribution of the site-entropy–RSA correlations in hybrid designs was comparable to that in natural sequences (Fig. 6). However, predictions for specific proteins lacked accuracy (Fig. S9).