Differential bicodon usage in lowly and highly abundant proteins

Data sources

In this work, I have used two kinds of data: (i) genome-wide PA across nine model organisms, and (ii) nucleotide sequences associated with the proteins indicated above. The absolute PA data from three prokaryotes (Microcystis aeruginosa, Bacillus subtilis, and E. coli), one unicellular fungus (S. cerevisiae), one plant (Arabidopsis thaliana), two multicellular eukaryotes (Drosophila melanogaster and Caenorhabditis elegans), and two mammals (Mus musculus and Homo sapiens) were downloaded from the PaxDb web site (http://pax-db.org/) on May 2015 (Wang et al., 2012). From these comprehensive data sets, I selected one sample of the most abundant proteins, and another sample of the less abundant ones. When more than one isoform was present in the comprehensive data set, only one was included in the samples. It is important to point out that PA correlates negatively with coding-sequence length in yeast (Coghlan & Wolfe, 2000). As PA distributions are generally biased, i.e., short proteins should be more abundant than larger ones (see Fig. 1), I selected two sets of 500 sequences, in which the sequence length distribution of both sets was similar. The procedure for sampling sequences with similar length distributions, consisted in ordering all the sequences which corresponded to a given organism in a PA increasing order. To select the sequences of the low PA samples, I began from the lowest PA extreme of the sequence list, and I compared $r_l=exp[-{(l-l_o)}^2/2\sigma ]$ (where l is the length of the sequence, l_o and σ are the mean length and standard deviation (SD), respectively, of the target distribution) with a random number uniformly distributed r. If r_l > r, the sequence was added to the set of low PA sequences. Then, I tested the second sequence in a similar manner and so on, until 500 sequences were selected. To select the high PA sequence set, I performed the same procedure, but beginning from the highest PA extreme of the sequence list. The resulting lists of coding sequences corresponding to these unbiased samples are given in Tables S1 and S2. In Fig. 1 (and also in Figs. S1–S8), I have plotted the distributions of the whole PA for all organisms used in this study, and the PA and the sequence length distribution of the selected data sets. The nucleotide coding sequences corresponding to the selected proteins were downloaded from Ensembl web sites (four eukaryote organisms from ftp://ftp.ensembl.org/pub/ and four prokaryote organisms from http://bacteria.ensembl.org), while A. thaliana coding sequences were downloaded from www.arabidopsis.org.

I also used the cumulative ribosome occupancy computed by Gamble et al. (2016) from the yeast ribosome profiling data set obtained by Jan, Williams & Weissman (2014). The cumulative ribosome occupancy is defined as the ratio between the sum of joint counts at positions with the bicodon in the ribosomal P-, A-sites and E-, P-sites and the joint counts sum across all surrounding window positions, which were extracted from Table S5 of Gamble et al. (2016).

Statistical analyses

Bicodon bias was studied in the context of the low and high PA samples. I counted all consecutive pairs of codons on the same reading frame of the coding sequences belonging to a given sample, which allowed us to compute the occurrence of each bicodon ij for all the sequences of each sample. The index i indicates the codon corresponding to P-site, while j indicates the one corresponding to the A-site. The occurrence of the codon pair ij will be denoted by o_ij. I also computed the number of single codons f_i for all the sequences of each sample.

In summary, I analyzed the bias of bicodon usage in the two samples using three complementary measures: (i) the pause propensity score, which is based on the differential bicodon usage in both samples; (ii) the Fisher’s exact test, which establishes whether the bicodon usage bias is significant; and (iii) the residual score proposed in Gutman & Hatfield (1989), which establishes whether the bias in bicodons can be explained by the codon usage bias, or not.

The pause propensity score, denoted by π, is defined as the difference between the relative synonymous bicodon usage computed over the low PA sequences (RSBU^L), and the one computed over the sequence sample associated with high PA (RSBU^H). Mathematically, ${\pi }_{ij}={\text{RSBU}}_{ij}^L-{\text{RSBU}}_{ij}^H=q_{ap}(f_{ij}^L-f_{ij}^H)/N_{ap}.$ (1)

Here, f_ij^X is the frequency of the bicodon ij computed over the sequence sample X, q_ap is the number of bicodons encoding for the same amino acid pair, and N_ap is the frequency of that amino acid pair for both samples. Thus, a large, or small, value of π_ij indicates the preference of bicodon ij for encoding low, or high, PA sequences, respectively. These values were clustered using a hierarchical average linkage over P-site and A-site codons according to patterns of similar bicodon preference.

Further, I use Fisher’s exact test to examine whether the number of occurrences of bicodon o_ij^L, observed in the sample of sequences associated with lowly abundant proteins was significantly different to the number of occurrences observed in the sample of sequences associated with highly abundant proteins o_ij^H. Thus, I constructed a 2 × 2 contingency table for each bicodon, as shown for an illustrative purposes in Fig. S9 for the particular case of bicodon AAGAAG (Agresti, 1992). In order to compute the p-value, I approximated the factorial operator with Stirling’s formula, $n!\approx \sqrt{2\pi n}{(n/e)}^n$ for n ≥ 25. To improve visualization, the colors in the heat maps are related to the quantity −Slog₁₀(p-value) where S takes value +1 or −1, when the bicodon has preference for sequences with low or with high PA, respectively.

Furthermore, I computed the expected number of occurrences of each codon pair, as $e_{ij}=f_i{f}_j{N}_p/N_{\text{tot}}^2$, where N_tot is the total number of codons in the set of sequences and N_p is the number of bicodons. Following (Gutman & Hatfield, 1989), I removed the contribution due to the nonrandomness of amino acid pairs by normalizing the former expected values as: ${\widehat{e}}_{ij}=e_{ij}\times \frac{{\displaystyle {\sum }_{kl}^{*}{O}_{kl}}}{{\displaystyle {\sum }_{kl}^{*}{e}_{kl}}},$ (2)where the * indicates that the sum is only over codon pairs kl encoding the same amino acid pair encoded by the bicodon ij. From the observed and normalized expected bicodon counts recorded in a given sample X, I computed the residual scores χ^Xij for each codon pair as: ${\chi }_{Xij}^2=\frac{{\left(o_{ij}-{\widehat{e}}_{ij}\right)}^2}{{\widehat{e}}_{ij}},$ (3)where X indicates the sequence samples, i.e., X = L for low PA sample, or X = H for high PA sequence sample. These residual scores can be used to assess whether the bias in a given codon pair can be explained, or not, by the bias in codons and amino acids.

In order to statistically assess the residual scores, I performed a random shuffling control. From each sample of sequences, I generated a second random set of sequences by shuffling the order of codons (but preserving the stop codon at the end of sequence). This procedure removed the codon correlation but not the codon usage. Then, I computed the residual score of bicodons associated with this random sample. I repeated the above procedure 200 times and, finally for each bicodon ij I computed the mean value 〈χ_ij²〉_ran, and associated standard deviation SD_ran(χ_ij²). Thus, the residual scores of bicodon ij can be expressed as the number of SDs from the mean, i.e., $({\chi }_{ij}^2-{\langle {\chi }_{ij}^2\rangle }_{\text{ran}})/{\text{SD}}_{\text{ran}}({\chi }_{ij}^2)$. This procedure was performed for the low and high PA samples of sequences independently.

The above analysis was performed for all possible bicodons (61 × 64 = 3904), excluding stop:sense bicodons and the stop:stop bicodons. The values of frequencies, π, p-value, and χ² for all bicodons and the nine organism are listed in Table S3. The scripts for each analysis in the manuscript are available in Code S1.

Bicodon preferences

In the light of recent findings (Gamble et al., 2016), I am interested here in studying bicodons as determinants of translational elongation rate modulation in nine organisms. The working hypothesis was that sequences encoding abundant proteins should be optimized, in the sense of translation efficiency. This optimization could be reflected in the usage frequency of both codons and bicodons which encode sequences associated with high and low PA. Thus, I expected to capture this bias by statistical analysis of these sequence samples. To check this hypothesis, I selected a set of 500 coding sequences associated with highest abundance proteins, and another 500 coding sequences associated with lowest abundance proteins, in nine model organisms from different kingdoms. Before showing the whole analysis across several organisms, I begin with an illustrative example. Figure 2A shows the histogram of the bicodons (red bars) which encode for the amino acid pair VR, and which were obtained from 500 low PA sequences from S. cerevisiae.

One can observe that bicodon usage is not uniform, i.e., it is biased, and this could be due to the known bias observed at the codon level. However, the expected frequency associated with such bicodons (black bars, obtained by the product of each codon frequency) shows that, although some bicodon frequencies can be explained by the bias in codon usage (for example, bicodon GTGAGA), other bicodons have an associated usage frequency that is underrepresented (such as bicodon GTTAGA), or overrepresented (as bicodon GTACGA). This means that two consecutive codons used to encode a given amino acid pair can be correlated. A similar analysis can be performed with sequences associated with high PA, as shown in Fig. 2B, and in all other amino acid pairs. Evidence for nonrandom associations between codon pairs, even once codon bias and bias against specific amino acid pairings were subtracted, was previously reported in E. coli (Gutman & Hatfield, 1989), and in many other genomes (Buchan, Aucott & Stansfield, 2006). However, what is a new remarkable fact in Fig. 2, is the strong difference between the histograms computed for the low PA and high PA samples. Fisher’s exact test allows one to reject, with a high significance level, the null hypothesis that bicodons are equally used in sequences from the low PA and high PA samples. In the particular case of bicodon GTTAGA, the p-value is 1.73 × 10⁻⁹. By comparison of the histograms corresponding to both low PA (red bars) and high PA (orange bars), it can be seen that some bicodons, such as GTTAGA or GTCCGT, are much more frequently used in sequences with high PA than in sequences with low PA, while the frequency usage of bicodon GTAAGA has an inverse relationship (see Fig. S10 for a direct comparison between both samples). More interestingly, the author found that some bicodons, such as GTACGG, GTACGA, or GTGCGA, are poorly or never used to encode VR amino acid pair in high PA sequences. In particular, the last two bicodons have been identified as inhibitors of translation elongation in yeast (Gamble et al., 2016). Figure 2 only illustrates the particular case of VR pair which is not a special case. There are many other bicodons with low frequency usage in high PA, and the 20 bicodons which mediate strong inhibition of translation elongation reported in Gamble et al. (2016) are a subset. Figure 3 depicts a raster plot of relative synonymous bicodon usage computed over the low PA sequences (RSBU^L) vs. relative synonymous bicodon usage computed over the high PA sequences (RSBU^H). The inset illustrates the region where inhibitor bicodons are. All of them, with the exception of CTTCTG, have an associated RSBU^H value which is very small or zero. Besides these known inhibitor bicodons, many other bicodons exhibit usage bias between both samples, and these bicodons are far from diagonal line in the RSBU^L–RSBU^H plane (dashed line in Fig. 3). Bicodons with similar features could be expected in other organisms and, in order to assess the existence of this bicodon bias observed in S. cerevisiae in other organisms, I devised two alternative heat maps: one based on the pause propensity score, π, and another based on the p-value provided by the Fisher’s exact test.

Pause propensity score across nine organisms

Figure 4 depicts a heat map for S. cerevisiae associated with the pause propensity score of 3904 bicodons. This score quantifies bicodon preference for coding low or high PA sequences. The position of each bicodon in the grid map was determined by clustering codons with similar pairing preference on both axes to facilitate the data visualization. Thus, bicodons with a clear preference for sequences associated with high PA (blue cells) have been grouped on the right-bottom corner of the grid, while red cells (high π score) indicate bicodons with preference for coding sequences associated with low PA. The author noticed that there are more bicodons with preference for encoding low PA sequences, than bicodons with preference for encoding high PA sequences, but the latter group has a stronger bias. Inhibitor bicodons identified in Gamble et al. (2016) are indicated by asterisks. All of them are associated with positive π scores, but not the highest π scores, which reflects the fact that these bicodons also have a smaller frequency associated with low PA sequences than other bicodons with preference for low PA. Another interesting feature of the heat map depicted in Fig. 4, is that some bicodons have an opposite preference depending on the order of the codons. For example, this is the case, to mention a few, for bicodons AGA-TGT and AGA-GCA (indicated by black circles), which have a preference for encoding low PA sequences (yellow cell), while the reverse ordered codon pairs TGT-AGA and GCA-AGA, respectively, have a preference for encoding high PA sequences (light blue cell). This fact cannot be explained solely in terms of codon usage bias.

I also elaborated similar heat maps to the one depicted in Fig. 4, for the other eight organisms which are shown in Fig. 5. For the sake of comparison, the position of each cell has been preserved over all organisms. These maps reveal that discrepancy in bicodon usage between low and high PA samples is not the same across the studied organisms, but is a particular feature of each organism, the same as codon bias. However, C. elegans and S. cerevisiae have many bicodons with the same preference, a feature that is also shared to a lesser extent by D. melanogaster and A. thaliana. These similarities and discrepancies are more apparent when clustering bicodons for all nine organisms (see Fig. S11). Regarding the inhibitor bicodons identified in yeast (indicated by an asterisk in Fig. 5), one can see that many of them share a preference for low PA sequences with other organisms, in particular in the cases of C. elegans, A. thaliana, and D. melanogaster. On the other hand, both studied mammals (H. sapiens and Mus musculus), and E. coli, have less apparent preferences. Interestingly, it was in E. coli where bicodon bias was first reported (Gutman & Hatfield, 1989). However, it is clear that bicodon usage bias across the ORFeome does not imply a different usage preference between highly and lowly expressed proteins. The low preference of bicodon usage observed in E. coli can be a consequence of the fact that, in this particular case, there is not a clear distinction between PA distributions for both samples (see Fig. S6).

The above heat maps show that bicodon usage bias exists in lowly and highly abundant proteins, but they do not indicate whether this bias is significant. For this, I applied Fisher’s exact test to compute the p-value against the null hypothesis of equal distribution in both low and high PA samples. Figure 6 depicts the corresponding heat map associated with S. cerevisiae, where the color of each cell in the grid is determined by the quantity −Slog₁₀(p-value), where S takes the values +1 or −1 when the bicodon has preference for sequences with low or with high PA, respectively. This means that red cells indicate bicodons with significant bias towards low PA sequences, while blue cells indicate bicodons with significant bias to high PA sequences. Again, to facilitate data visualization, the position of each bicodon in the map was determined by clustering codons with similar pairing preference on the P-site and A-site axes. Inhibitor bicodons are indicated by asterisks. Some of them (like bicodons CGAGCG, CTCCCG, and CGACGG) are associated with a not significant bias (p-value < 0.01). This is due to the fact that these bicodons have very small frequency in low PA sequences and almost zero in high PA sequences. Something similar occurs for the sense:stop bicodons, which are grouped in three almost white rows at the center of the plot. For the sake of comparison, I also elaborated similar heat maps for other organisms where the positions of each cell has been preserved over all organisms (Fig. S12). A couple of white rows can be observed in all organisms. These cells correspond to bicodons which do not exhibit any preference or they are usually poorly used in both sequence samples and have poor associated statistics.

The ribosome profiling method has been used for obtaining valuable information about protein synthesis, such as the location and strength of translational attenuation. Thus, sites associated with high levels of ribosome occupancy can be linked to translational pauses (Ingolia et al., 2009; Ingolia, Lareau & Weissman, 2011). To determine if ribosome occupancy correlates with the proposed pause propensity score, I used the cumulative ribosome occupancy for each bicodon, which was computed using the yeast ribosome profiling data (Table S5 of Gamble et al., 2016). Figure 7 depicts a scatter plot of the cumulative ribosome occupancy vs. the pause propensity score. Despite the fact that a clear correlation between these descriptors is not apparent in these plots, one can see that in addition to the inhibitor bicodons identified in Gamble et al. (2016), indicated by red circles, there are some other bicodons which have high levels of cumulative ribosome occupancy (> 0.03) and a significant preference for encoding low PA sequences (π > 0.75 and p-value < 0.01). These bicodons could be also candidates for modulating transcript elongation rate in yeast (they have been indicated by ** in Table S3).

Codon bias cannot explain bicodon bias

The heat maps shown above are very useful to see some common features among organisms, but they do not show whether bicodon preference is explained or not by the preference of the codon in the pair for sequences associated with low or high PA. In order to study this, I have computed residual scores for each bicodon over sequences with low PA, χ_L², and over sequences with high PA, χ_H². When the residual score is high, bicodon usage cannot be explained by codon usage in the same sample of sequences, as it has been previously established in Gutman & Hatfield (1989) and Buchan, Aucott & Stansfield (2006). The value of these residuals scores for all codon pairs and organisms are listed in Table S3. In Fig. 8, a raster plot of these residual scores can be seen for all bicodons in B. subtilis, yeast, humans and E. coli (residual plots associated with the other five organisms are displayed in Fig. S13). As there are two sequence samples for each bicodon, it is convenient to take the quantity ${\chi }^2={\chi }_L^2+{\chi }_H^2$ as a whole residual score. For all organisms, I have considered that those bicodons with χ ≥ 3 SD above the mean are bicodons whose preference for encoding low or high PA sequences is not explained by the preference of the individual codons in the pair (at a significance level of 0.003).

Thus, four types of bicodons can be distinguished in Fig. 8: (i) codon pairs that are significantly more used in sequences associated with low PA than in sequences associated with high PA (p-value < 0.01), and whose preferences cannot be explained by the codon usage bias (red dots); (ii) codon pairs which are significantly more used in sequences associated with high PA than in sequences associated with low PA (p-value < 0.01), and whose preferences cannot be explained by codon usage bias (blue dots); (iii) codon pairs with a significantly different usage frequency in low and high PA samples, but whose preferences can be explained by codon usage bias, i.e., χ² < 3 × SD (green dots) and, finally, (iv) codon pairs whose usage frequencies in low and high PA samples are not significantly different, i.e., p-value > 0.01 (black dots).

These plots indicate that while there are many bicodons with evident preference for low and high PA sequences in C. elegans and S. cerevisiae, there are only few in H. sapiens at this significance level. Humans, mouse and, surprisingly, also E. coli, have few bicodons with evident preference. However, even in these organisms, there are numerous bicodons with similar usage frequencies in both samples, but whose associated frequency usage cannot be explained by codon usage bias (black dots out of the gray quadrant). This fact reveals that, even in these organisms, there is a complex correlation between two consecutive codons.

On the other hand, S. cerevisiae has more codon pairs than B. subtilis, with a significant different usage frequency in low and high PA samples. Nevertheless, the bias of many bicodons is explained by codon usage bias (green dots in S. cerevisiae panel).

Regarding the inhibitor bicodons, I found that only five of them have a significantly different usage frequency in low PA sample and only one in high PA sample (see Table S3). I searched for bicodons with the same preferences shared by S. cerevisiae, C. elegans, and D. melanogaster, by selecting those with p-value < 0.01 and χ² ≥ 3 × SD above the mean. Table 1 lists 16 shared bicodons with preference for low PA, while Table 2 lists 40 bicodons with preference for high PA sequences. It is noteworthy that, even within each species, the number of bicodons with preference for encode high PA sequences (blue cells in Fig. 5) are smaller than the number of bicodons with preference for encoding low PA sequences (yellow cells in Fig. 5), and the number of shared bicodons in Table 2 is substantially greater than those listed in Table 1. This fact seems to indicate that bicodons with preference for high PA sequences are more conserved across these organisms than bicodons with preference for low PA sequences.

Pause propensity score and cotranslational protein folding

The above statistical analysis is able to determine which bicodons are associated with lowly or highly abundant proteins. I hypothesized that bicodons associated with lowly abundant proteins could have a key role in programming translation pauses of the ribosomal machinery. As a proof of principle, I have identified one example that illustrates how even small biases as those observed in H. sapiens, could explain documented protein misfolding and a pathological condition in humans. It was previously reported that the single polymorphism (SNP) rs1045642 in the gene MDR1, which encodes the ABCB1 drug-transport protein, alters the structure and function of the protein with the consequent change in its substrate specificity (Kimchi-Sarfaty et al., 2007).

The SNP in exon 26 at position 3435 changes codon ATC to the synonymous ATT, which reduces codon usage from 47 to 35%. It was argued that the presence of a rare codon affects the timing of co-translational folding and insertion of P-glycoprotein into the membrane. Although it is difficult to consider the ATT codon as rare, it is clear that the SNP alters the timing of ribosomal transit. Here, I offer an alternative cause for translational attenuation. In this sense, I have observed that this SNP is also associated with a large change in the bicodon pause propensity score π. Specifically, bicodon ATCGTG has preference for low PA sequence with π = 0.1, while bicodon ATTGTG has preference for high PA sequence π = −0.8. This means a large change in comparison with the change in codon usage. Further, the other synonymous bicodon ATAGTG has a low pause propensity score, π = −0.55, in agreement with Kimchi-Sarfaty et al. (2007) observation that associates this haplotype to a larger decrease in the inhibitory effect. In addition to the SNP mentioned above, there are other synonymous SNPs related to human diseases that could be explained by a large change in the pause propensity. Among these, one can mention SNP rs34533956 in the gene CFHR5, which is associated with age-related macular degeneration (Narendra, Pauer & Hagstrom, 2009). In this case, the mutation changes bicodon GACGTG to GATGTG, and the associated change in pause propensity is from 0.03 to −0.78, while the change in the relative synonymous codon usage is only 13%. Another example corresponds to SNP rs11615 in the gene ERCC1, which was associated with colorectal cancer (Liang et al., 2010), where the pause propensity change is also large, against a small change in the relative synonymous codon usage. These relationships suggest that some pathological synonymous mutations could be understood in terms of the change in timing needed for co-translational folding programmed by bicodons.