SCUMBLE: a method for systematic and accurate detection of codon usage bias by maximum likelihood estimation

Abstract

The genetic code is degenerate--most amino acids can be encoded by from two to as many as six different codons. The synonymous codons are not used with equal frequency: not only are some codons favored over others, but also their usage can vary significantly from species to species and between different genes in the same organism. Known causes of codon bias include differences in mutation rates as well as selection pressure related to the expression level of a gene, but the standard analysis methods can account for only a fraction of the observed codon usage variation. We here introduce an explicit model of codon usage bias, inspired by statistical physics. Combining this model with a maximum likelihood approach, we are able to clearly identify different sources of bias in various genomes. We have applied the algorithm to Saccharomyces cerevisiae as well as 325 prokaryote genomes, and in most cases our model explains essentially all observed variance.

Figure 1.

(a) Cumulative histogram of the normalized variance for named genes in S. cerevisiae for models with various numbers of trends; actual genome (solid lines) compared to randomized genome (dotted lines). Models with 0 or 1 trend explain the data poorly, as the curve for the real genome is very different from that of a randomized genome, and there are many genes with very high normalized variance. (b) Average (black) and median (red) normalized variance for models with up to 10 trends.

Figure 2.

Experimental values for cellular mRNA/protein levels plotted against the first offset/CAI value of each gene for S. cerevisiae. Several groups of highly expressed genes are plotted in different colors.

Figure 3.

Median normalized variance for 325 prokaryote genomes, using models with 0–10 trends. The different genomes are slightly offset along the abscissa, in alphabetical order. The dotted brown line shows approximate median normalized variance for randomized genomes generated from the models (Supplementary Figure S7). Results for the average normalized variance are very similar, except that in rare but not exceptional cases, individual genes dominate the average due to extremely low estimated probabilities of using a specific codon which is, in fact, used.

Figure 4.

A four-trend model of Helicobacter pylori. (a)–(c) GC3 or GT3 plotted against the first three offsets. Genes for ribosomal proteins are circled in red. The cumulative distributionsof the offsets are shown above each graph, for all genes (black) and for ribosomal genes (red). (d) β₂ plotted against the number of the gene along the genome, with genes on different strands in different colors. The green and blue lines are 50-point running averages for strand 1 and 2, respectively.

Figure 5.

Scatter plot of the first two axes from the four-trend model found by SCUMBLE (a), WCA (b) and CA/RSCU (c) for the genes of Anaeromyxobacter dehalogenans. Genes for ribosomal proteins are circled in red. In (b) and (c), most genes are clustered near the origin; only a small fraction of the genes have significantly negative abscissae.

Figure 6.

Solid lines: number of prokaryote genomes (out of 325) for which the total fraction of the GC (a), GT (b), CT (c) or random (d) preference signal captured by the first n trends exceeds the abscissa, where n is given by the color. Total shaded area of each color is proportional to the average fraction of signal captured by the corresponding trend.

Figure 7.

Scatter plot of the first two offsets for the four-trend model of B. subtilis, with the genes' colors given by their cluster identity given in ref. (17).