text:
# Genetic Algorithms for the Imitation of Genomic Styles in Protein Backtranslation

## 1 Introduction

The main components of the cell are nucleic acids (DNA and RNA) and proteins. Both are polymers, long words written in alphabets of 4 and 20 letters: 4 nucleotides for DNA and RNA, and 20 amino acids, for proteins. The “fundamental dogma” of molecular biology describes the usual flow of information in the cell, from DNA to mRNA to protein. The first step, transcription, preserves the sequence read from DNA, which is reversed and complemented in the mRNA (in addition, the alphabet is slightly changed). It is straightforward to obtain the DNA from a given mRNA (it is called then complementary DNA, or cDNA); in fact, Nature does it: retrotranscription is performed by viruses and several small “selfish” units of information.

The second step, translation, is more complicated: the mRNA is read, three nucleotides at a time, and an amino acid encoded by them is added to the forming protein, according to the well known genetic code (see Table 1). This nearly universal code associates to each triplet (codon) an amino acid, or the “stop” meaning.

Table 1: The (standard) Genetic Code

| aaa | K | aga | R | caa | Q | cga | R | gaa | E | gga | G | taa | stop | tga | stop |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| aac | N | agc | S | cac | H | cgc | R | gac | D | ggc | G | tac | Y | tgc | C |
| aag | K | agg | R | cag | Q | cgg | R | gag | E | ggg | G | tag | stop | tgg | W |
| aat | N | agt | S | cat | H | cgt | R | gat | D | ggt | G | tat | Y | tgt | C |
| aca | T | ata | I | cca | P | cta | L | gca | A | gta | V | tca | S | tta | L |
| acc | T | atc | I | ccc | P | ctc | L | gcc | A | gtc | V | tcc | S | ttc | F |
| acg | T | atg | M | ccg | P | ctg | L | gcg | A | gtg | V | tcg | S | ttg | L |
| act | T | att | I | cct | P | ctt | L | gct | A | gtt | V | tct | S | ttt | F |

Unlike retrotranscription, the reversal of this second step (called backtranslation) is ambiguous, due to the degeneracy of the genetic code: as can be seen in Table 1, amino acids are encoded by 1, 2, 3, 4 or 6 different codons. Backtranslation does not occur in natural systems111 Though  suggests that it did occur at the origin of life, and even proposes an in vitro device for backtranslation., but is required for several purposes in genomics and biotechnology. The problem is not trivial, since different species have different “genomic styles” that determine which of the many preimages is used to code for a protein. Thus it may happen that we know the DNA for a given protein produced by, for instance, a plant, but we want to synthesize the protein in a bacterium. We will need to backtranslate the protein into the genomic style of this kind of bacteria. In other cases, the protein is known but no DNA is known for it at all; this may happen with artificial proteins, or with proteins from unsequenced organisms. Other applications, like degenerate primers (for “gene fishing”) and sequence analysis, will be discussed in the last section.

The best known statistical feature of coding sequences is the presence of a periodicity of period 3, which is caused by the structure of the genetic code and the asymmetry of the different codon positions. This property is very important for distinguishing coding from non-coding sequences; however, it is not important for backtranslation, since it is shared by all organisms. On the other hand, we know that codon usage (the degree of preference for the different codons inside each synonymous class) does distinguish one species from another; it is the best known feature of the different “genomic styles”.

The common approach to backtranslation relies on the imitation of the codon usage of the target species (the species whose style we want to imitate). This is the solution currently given by all commercial and non-commercial software, like GCG, EMBOSS, VectorNTI, EditSeq, AiO, and the online tools of Molecular Toolkit and Entelechon. The only different approach we know is , where a neural network was trained to perform backtranslation. However, it was done at the single amino acid level, and thus it cannot account for anything but codon usage.

This current solution can be improved; there are more features peculiar to the different coding styles, which are in part or completely independent from codon usage. In the present article, we consider different possible statistics that may be associated to genomic styles, and then we apply a genetic algorithm to perform backtranslation, taking these features into account. Our approach considers DNA only as a symbolic sequence, ignoring chemical properties or biological features. Furthermore, we will not use biological considerations to decide whether or not a statistical property needs to be imitated: we assume that any property distinguishing the style of a species must be considered in backtranslation (after all, in some cases the origin of known features remains obscure). All the statistics we consider were taken from the literature on sequence analysis, where their possible interpretations are discussed.

## 2 Notation, Materials

Let $`A`$={A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y} and $`B`$={$`a`$,$`c`$,$`g`$,$`t`$} be the alphabets for amino acids and nucleotides, respectively, and denote $`B^3=(B^3)^{}`$. Let $`\tau :B^3(A\{stop\})^{}`$ be the translation of a sequence according to the genetic code. In fact, $`\tau `$ may depend on small variations to the code which do occur in some species and organelles; however, here we will assume the code to be universal. Furthermore, we will consider the sequences without the start and stop signals, i.e., cutting the $`atg`$ codon that initiates a protein and the stop codon that marks its end.

We will say that a function (or stochastic procedure) $`\beta :A^{}B^3`$ is a guess iff $`\tau \beta =id_A^{}`$. If $`CA^{}`$, we will denote $`\beta (C)=\{\beta (u):uC\}`$. A particular guess that will be used for comparison purposes is the canonical backtranslation procedure, which backtranslates each amino acid using the empirical frequencies of its codons as probabilities; we will denote it as $`\beta _{sp}^{cu}`$, with the subindex indicating the species whose codon usage table was used.

Given a sequence $`wB^3`$, $`w=w_0,w_1,\mathrm{}`$ and $`i=1,2,3`$, we will talk about the letters in codon position i to refer to $`w_{i1}`$, $`w_{i+2}`$, $`w_{i+5}`$, …. We will denote with $`\pi _{ry}`$, $`\pi _{ws}`$ and $`\pi _{mk}`$ the three most usual projections of $`B`$ into $`\{0,1\}`$, as follows. We will use the same symbols to refer to the extensions of these functions to $`B^3`$ (projecting each letter).

|  | $`a`$ | $`c`$ | $`g`$ | $`t`$ | refers to: |
| --- | --- | --- | --- | --- | --- |
| $`\pi _{ry}`$ | 0 | 1 | 0 | 1 | purine/pyrimidine |
| $`\pi _{ws}`$ | 0 | 1 | 1 | 0 | weak/strong |
| $`\pi _{mk}`$ | 0 | 0 | 1 | 1 | amino/keto |

It is important to notice that many characters in $`\beta (u)`$ are almost or completely determined by $`u`$. Amino acid K, for instance, is coded by $`aaa`$ and $`aag`$; the first and the second position will be $`a`$ in any backtranslation, and the third one will be either $`a`$ or $`g`$ (and will have $`\pi _{ry}=0`$, so that for any $`\beta `$, $`\pi _{ry}(\beta (K))=001`$). The next table shows the number of amino acids for which characters are fixed in the different codon positions for the different binary alphabets. Most of the ambiguity of backtranslation is in the third position.

|  | $`\pi _{ry}`$ | $`\pi _{ws}`$ | $`\pi _{mk}`$ |
| --- | --- | --- | --- |
| Cod. Pos. 1 | 18 | 18 | 18 |
| Cod. Pos. 2 | 19 | 20 | 19 |
| Cod. Pos. 3 | 11 | 2 | 2 |

Materials

We extracted coding sequences from Genbank release 131 (August 2002), belonging to the following species: Methanosarcina acetivorans C2A ($`A1`$), Sulfolobus solfataricus ($`A2`$), Escherichia coli ($`B1`$), Bacillus subtilis ($`B2`$), Streptomyces coelicolor A3(2) ($`B3`$), Mesorhizobium loti ($`B4`$), Nostoc sp. PCC 7120 ($`B5`$), Saccharomyces cerevisiae ($`E1`$), Arabidopsis thaliana ($`E2`$), Drosophila melanogaster ($`E3`$), Caenorhabditis elegans ($`E4`$) and Homo sapiens ($`E5`$). The selection of species was done trying to have abundant sequences and a rather good representation of the tree of life. All coding sequences (“CDS” features in Genbank) were extracted, provided that they were complete, univoque, and longer than 1029 nucleotides. The average length of the sequences varies between 1500 for $`A1`$ and 2456 for $`E3`$. Please notice that introns -intervening sequences- were removed from the sequences; this may affect the coding statistics that depend on relations between distant nucleotides. We will use the abbreviation of a species to refer to the set of its coding sequences, or to the set of the corresponding proteins, depending on the context. Thus, an expression like $`\beta _{B1}^{cu}(E5)`$ denotes a set of backtranslations obtained for all proteins encoded by the coding sequences of $`E5`$, obtained by the standard backtranslation method, considering the codon usage of $`B1`$.

## 3 Coding Statistics

Here we discuss the results of computations performed on our set of species for several features that have been studied in coding sequences, “generally known as coding statistics, since their behavior is statistically distinct on coding and non-coding regions”. Discussions about the most common coding statistics, their relations, and their use for gene finding, can be found in  and . However, we are not interested in the difference between coding and non-coding regions; rather, we want those statistics that contribute to the “genomic style” of a species.

The notion of genomic style has been around since the “genome hypothesis” of Grantham , who first recognized the idiosyncratic nature of codon usage. Later, Karlin used the bias in dinucleotide usage as the “genomic signature” of a species . Forsdyke suggests that the species “broadcast” their genes in different $`g+c`$ frequencies , and that this could be crucial for speciation; in this way, genomic styles could be the first line of an immune system222Indeed,  shows that some viruses may mimic the genomic style of their host, in order to be expressed.. There have been other proposals, usually for phylogenetic purposes. The reasons for the existence of different styles are debatable: for instance, changes in the molecular machinery, tRNA abundance, environmental temperature, different biases in the mutation rates, the requirements of messages other than the protein sequences, etc. The exact causal relations are subject to discussion.

In order to improve the profile of genomic styles, we want to choose those statistics which: (1) have typical and statistically sound values for each species, with small variability, (2) have different values in different species, and (3) do not depend (exclusively) on the amino acids encoded by a sequence (i.e., they do depend on backtranslation). Because of space limitations, we will not give the values of all computations; in the graphics, not all the species will be displayed, if it is not required. Moreover, we will dispense from data in the case of well known facts. All computations and data sets can be found at .

### 3.1 Nucleotide frequencies

The most natural computation is the frequency of the four nucleotides in the sequences, as well as their frequencies in the different codon positions. For each sequence $`wB^3`$, $`w=w_0,\mathrm{},w_{3N1}`$, and each nucleotide $`\alpha `$, we compute

$$\rho _\alpha (w)=\frac{1}{3N}\underset{i=0}{\overset{3N1}{}}\delta _\alpha (w_i),\rho _\alpha ^j(w)=\frac{1}{N}\underset{i=0}{\overset{N1}{}}\delta _\alpha (w_{3i+j1}),j=1,2,3$$

where $`\delta _\alpha (x)`$ is $`1`$ if $`x=\alpha `$ and 0 otherwise. Our computations confirm a number of facts already known in the literature, like “Chargaff’s second law”, which states that $`\rho _a\rho _t`$ and $`\rho _c\rho _g`$ as can be observed in Graphic 1a. Since, in addition, $`_{\alpha B}\rho _\alpha =1`$, Chargaff’s law implies positive correlation between complementary nucleotides ($`a`$ with $`t`$, and $`c`$ with $`g`$) and negative correlation between non-complementary ones. Thus we can reduce the study to a single value; the usual choice is $`\rho _{g+c}=\rho _c+\rho _g`$. It is well known that $`\rho _{g+c}`$ has different values in different species, and that all the genes in a species have similar values; this can be seen in Graphic 1b, with histograms showing the number of sequences of each species in different $`\rho _{g+c}`$ ranges. Some qualifications are due: First, it is also known that eukaryotic genomes are organized in large “islands” called isochores , with different $`\rho _{g+c}`$ values but each of them relatively homogeneous. Moreover, in a set of closely related species $`\rho _{g+c}`$ may depend more on the genes than on the species. However, the general pattern holds, and it is used both for the detection of genes (since genes tend to be $`\rho _{g+c}`$-richer than non-coding regions) and in the detection of horizontally transferred genes (see section 5).

Graphic 1c shows the values of $`\rho _{g+c}^j=\rho _c^j+\rho _g^j`$ for the different species, together with $`\rho _{g+c}`$. We notice the existence of wide variations in the $`\rho _{g+c}`$ composition depending on the codon position. In addition, extreme values of $`\rho _{g+c}`$ are usually supported by extreme values of $`\rho _{g+c}^3`$; this shows that the sequences were adapted to get a certain $`\rho _{g+c}`$ level, and that the third -usually synonymous- codon position was used for this purpose. As can be seen in Table 2, $`\rho _{g+c}^1`$ and $`\rho _{g+c}^2`$ are almost entirely determined by the encoded amino acids.

### 3.2 Codon usage

The frequency of a given codon $`C=c_0,c_1,c_2B^3`$ in a sequence $`w=w_0,\mathrm{},w_{3N1}B^3`$ is defined as $`\frac{1}{N}_{i=0}^{N1}\delta _{c_0}(w_{3i})\delta _{c_1}(w_{3i+1})\delta _{c_2}(w_{3i+2})`$. For each codon $`CB^3`$, we define its synonymous class $`\theta (C)=\{C^{}B^3:\tau (C)=\tau (C^{})\}`$. Then the synonymous codon usage and the relative synonymous codon usage  of $`C`$ are defined as

$$SCU_C=\frac{\rho _C}{{\displaystyle \underset{C^{}\theta (C)}{}}\rho _C^{}},RSCU_C=\frac{|\theta (C)|\rho _C}{{\displaystyle \underset{C^{}\theta (C)}{}}\rho _C^{}}=|\theta (C)|SCU_C$$

As we mentioned above, the codon choice pattern was noted very early to be a signature of the species, and our data confirm this. We will dispense with extensive SCU tables, since they are well known in the literature, and available in public databases. As we said before, the common approach to backtranslation uses SCU as the probability of choosing a certain codon, given the amino acid. RSCU is used for comparisons between codons from different synonymous classes.

### 3.3 Dinucleotides

Most published results on dinucleotide frequencies consider long DNA sequences, including both coding and non-coding regions . Our own computations, in spite of being limited to coding sequences, confirm most of the facts already noted by the different authors. This accounts for the fact that dinucleotide frequencies are not considered as “coding statistics”: their behavior is similar in coding and in non-coding sequences. However, they do exhibit characteristic patterns according to the different species and groups. Karlin  even used them to define the genome signature of a species as the collection $`\{\varrho _{\alpha \beta }^{}\}`$, with $`\alpha `$ and $`\beta `$ ranging over $`B`$. Here $`\varrho _{\alpha \beta }=\rho _{\alpha \beta }/\rho _\alpha \rho _\beta `$ (with $`\rho _{\alpha \beta }`$ being the frequency of the dinucleotide $`\alpha \beta `$) and $`\varrho ^{}`$ is the computation of $`\varrho `$ over the sequence concatenated to its inverse complement (in order to get the information about both DNA strands).

IDH.  There is an interesting set of indices which can be computed from dinucleotide frequencies. The so called index of DNA homogeneity (IDH) was proposed by Miramontes et al  and is defined for a binary sequence as $`d=\frac{\rho _{00}\rho _{11}\rho _{01}\rho _{10}}{\rho _0\rho _1}`$. We define $`d_{ry}(w)=d(\pi _{ry}(w))`$, $`d_{ws}(w)=d(\pi _{ws}(w))`$, and $`d_{mk}(w)=d(\pi _{mk}(w))`$. This index expresses the degree of local homogeneity of the sequence: long stretches of 0 or 1 will cause $`d`$ to be near 1, while strong alternation will push it toward -1. The three indices $`d_{ry}`$, $`d_{ws}`$ and $`d_{mk}`$ are not independent, and since $`\pi _{mk}`$ is the least meaningful of the binary projections, the choice in  was to plot the species in the $`(d_{ry},d_{ws})`$ plane. The corresponding map with our own data is in Graphic 2a. Graphic 2b displays the distribution of the values in the sequences of some species. Both the specificity and the classificatory power of IDH can be clearly noted.

### 3.4 Fourier harmonics and Periodicities

Another common tool for DNA analysis is the discrete Fourier transform. For a binary sequence $`w=w_0,\mathrm{},w_{N1}`$, we define the spectrum and its $`m`$-smoothed version:

$$S_n(w)=\frac{1}{N^2}\left|\underset{k=0}{\overset{N1}{}}w_ke^{\frac{2\pi ink}{N}}\right|^2\text{ , }\stackrel{~}{S}_n^m(w)=\underset{k=nm}{\overset{n+m}{}}\frac{S_k(w)}{2m+1}$$

$`S_n(w)`$ measures the frequency content of ‘frequency’ $`n`$, which corresponds to a period $`\frac{N}{n}`$; the smoothed value helps to remove the dispersion that appears for small data sets.

The main and better known periodicity in DNA sequences is of period 3; it can be explained by the asymmetry in the codon positions , though its presence in tRNA genes suggests some other origin. Another well documented periodicity is of period 10.5 $`\pm `$ 0.5; it has been attributed to requirements from the structure of both DNA and proteins, and the exact contribution of each is unclear. Some periodicities of higher periods have been shown, but they are not statistically significant for the typical lengths of genes.

We divided each sequence in non-overlapping windows of length 256, and used the fast Fourier transform (FFT) algorithm to compute $`\stackrel{~}{S}^5\pi _{ry}`$, $`\stackrel{~}{S}^5\pi _{ws}`$ and $`\stackrel{~}{S}^5\pi _{mk}`$ for all the species. The results were averaged and are shown in Graphics 3a, 3b and 3c for some of the species; only part of the ordinate axis is used, in order to highlight their differences. The two periodicities mentioned before are present: there is a big peak at $`n=85`$ for the three projections in almost all the species (the top of the peaks is outside the graphics); this corresponds to a period of $`\frac{256}{85}3`$. There is also a minor peak around $`n=24`$, present for most species and for most projections, corresponding to the period $`\frac{256}{24}10.5`$; there are some differences between species, a fact that has been observed before and is related to the various origins of this periodicity.

To show the specificity of the spectrum, we chose a set of 20 collections of sequences, each set selected at random to be 1% of $`E5`$. We computed the average of spectra for each set; the results for $`\pi _{ws}`$ are shown in Graphic 4a.

Position dependent spectra. To take into account the asymmetry of the different codon positions, we computed the spectra for the three subsequences $`w_n^{(i)}=w_{3n+i}`$, $`i=0,1,2`$, using windows of length 64 (data not shown). In absence of period 3, the most notorious feature is a peak at $`n=18`$, corresponding to a period $`\frac{64}{18}3.5`$ in the subsequence, and hence 10.5 in the sequences; it is by far stronger for the middle codon position, a fact that hints for dependence on the amino acid sequence.

### 3.5 Autocorrelation functions

Correlation functions  measure the excess or defect of nucleotides at different distances; if $`\rho _{\alpha ,\beta }(d)`$ is the frequency with which we find a ‘$`\beta `$$`d`$ positions after a ‘$`\alpha `$’, then what we compute is $`\rho _{\alpha ,\beta }\rho _\alpha \rho _\beta `$. More precisely, what we compute for a sequence $`w=w_0,\mathrm{},w_{N1}`$ is

$$\mathrm{\Gamma }_{\alpha ,\beta }(d)[w]=\frac{1}{Nd}\underset{i=0}{\overset{Nd1}{}}\delta _\alpha (w_i)\delta _\beta (w_{i+d})\rho _\alpha (w)\rho _\beta (w)$$

We computed $`\mathrm{\Gamma }_{0,0}`$ for $`\pi _{ry}`$, $`\pi _{ws}`$, $`\pi _{mk}`$. The most notorious result of this computation is the strong oscillation due to period 3; this can be removed by considering the smoothed version, $`\stackrel{~}{\mathrm{\Gamma }}_{\alpha \beta }(d)=\frac{1}{3}_{i=d1}^{d+1}\mathrm{\Gamma }_{\alpha \beta }(i)`$; when this was done, the periodicity of period 10.5 could also be seen. To give an idea of the shape of the curves, and to show their specificity, Graphic 4b shows the results for $`\pi _{ws}`$, for $`B1`$, $`E5`$, and the same subsets of $`E`$ used in Graphic 4a. In general, $`\mathrm{\Gamma }`$ behaves very similar to the Fourier transform, in specificity and in the dependencies on alphabet and/or codon position. This is no surprising, since both express the same information (if $`\mathrm{\Gamma }`$ is computed for a circular sequence, then it can be recovered form the spectra, and vice versa, by the Wiener-Khinchin theorem). Position dependent autocorrelation functions were also computed, with no unexpected results.

## 4 Backtranslation strategy

### 4.1 Genomic style beyond codon usage

We will consider all of the coding statistics reviewed in the previous section as features defining the genomic style of a species. It is important to notice that they are not (or not directly) dependent on the codon usage; if this were the case, then genomic style would reduce to RSCU, and the current approach to backtranslation would be already optimal.

It is clear that $`\rho _\alpha `$ and $`\rho _\alpha ^j`$ are recovered by RSCU, if the amino acid composition is kept constant (this is the case in $`\beta _{B1}^{cu}(B1)`$ and $`\beta _{E5}^{cu}(E5)`$); in general, since amino acid composition is rather similar in all the different species (data not shown), we can expect nucleotide frequencies to be conserved.

For dinucleotides, this is not so clear, even if the amino acid frequencies are kept: in spite of recovering the number of dinucleotides starting at the first and second codon positions, RSCU will not recover those starting at the third. This is important, since most of the degeneracy is in this position, and “genomic style” depends strongly on it; moreover, mutation rates tend to be affected by the neighboring nucleotides , in ways that are species-dependent. In particular, Miramontes et al  show that their indices (IDH) are not determined by codon usage, even when the amino acid frequency was conserved. Our data (not shown) confirm it.

As for the Fourier spectra, Guigó  shows that it is rather independent from $`\rho _{g+c}`$. To discard dependence on RSCU, we computed the spectra on $`\beta _{B1}^{cu}(B1)`$, $`\beta _{B1}^{cu}(E5)`$, $`\beta _{E5}^{cu}(E5)`$ and $`\beta _{E5}^{cu}(B1)`$; results for $`\stackrel{~}{S}_n^5\pi _{ws}`$ are displayed in Graphic 5a. We can see that all the sets of guesses lie between the real spectra, with codon usage being a bit more relevant than the amino acid sequences (the species); this was also the case for $`\pi _{ry}`$ and $`\pi _{mk}`$ (data not shown). Although the autocorrelation function contains the same information as the spectrum, the details of each one are the main lines of the other, and thus, each may be considered apart. Graphic 5b displays computations of $`\stackrel{~}{\mathrm{\Gamma }}_{0,0}\pi _{ws}`$ over the same sets; it can be noticed that in this case the species (amino acid sequences) are the major contribution, with only a small effect of RSCU.

### 4.2 Genetic Algorithms for Backtranslation

We want to obtain a backtranslation that imitates the genomic style of a target species as close as possible; thus, we will look for a backtranslation for which the coding statistics listed above are close to those of the target species, i.e., their distance is minimum. We choose, for $`wB^3`$,

$$\begin{array}{c}f_1(w)=|\rho _{g+c}(w)\rho _{g+c}^{}|f_2(w)=_{CB^3}|RSCU_C(w)RSCU_C^{}|\\ f_3(w)=|d_{ry}(w)d_{ry}^{}|+|d_{ws}(w)d_{ws}^{}|+|d_{mk}(w)d_{mk}^{}|\\ f_4(w)=_{k=3}^{125}a_k|\stackrel{~}{S}_k^5(w)\stackrel{~}{S}_k^{5,}|f_5(w)=_{d=2}^{99}b_k|\stackrel{~}{\mathrm{\Gamma }}_k(w)\stackrel{~}{\mathrm{\Gamma }}_k^{}|\end{array}$$

where the values with “\*” are obtained averaging over the known coding sequences of the target species, and $`a_k`$ and $`b_k`$ are weights, incorporated in order to give more importance to some parts of the curves, e.g. to encourage a uniform convergence. The indices in the sums of $`\stackrel{~}{S}`$ and $`\stackrel{~}{\mathrm{\Gamma }}`$ follow our particular choices of window lengths 256 and 30, respectively.

With these definitions, what we want, for a given $`uA^{}`$ and a given target species, is to minimize $`\stackrel{}{f}(w)`$, with $`w\tau ^1(u)`$. There are two main difficulties involved. First, we have a non-convex problem, in a vast search space, with terms depending on several scales of the sequences. Moreover, it is a problem of multiobjective optimization. For these reasons, we propose the use of genetic algorithms (GA), specially suited for problems with these characteristics. Our particular implementation of a genetic algorithm for backtranslation follows here.

* for $`1in`$ initialize $`w^i=\beta ^{cu}(u)`$
* while not stop condition

+ for $`1j5`$, $`\overline{f}_j=\mathrm{max}_if_j(w^i)`$
    + for $`1in`$, $`1j5`$, $`N_j^i={\displaystyle \frac{\overline{f}_jf_j(w^i)}{n\overline{f}_j_kf_j(w^k)}}`$
    + for $`1in`$, $`N^i=_j\lambda _jN_j^i`$
    + Update $`P`$ using $`\{N^i\}`$ \[stoch. univ. sampling\]
    + Apply genetic operators: crossover and mutation

For a given $`uA^{}`$, we iterate on a population of $`n`$ guesses of $`\tau ^1(u)`$, denoted by $`\{w^i\}`$. As seen in the scheme, our initial condition is the usual backtranslation (imitation of RSCU); the GA is iterated then to optimize coding statistics. $`N_j^i`$ are the expected number of copies of a guess in the next generation; ponderating them with $`\{\lambda _j\}`$ we combine the different objective functions, without needing to make their numeric values comparable. The genetic operators used are crossover and mutation, both adapted to maintain the encoded amino acid sequence $`u`$. In addition, the probability of crossover between two guesses $`w^i`$ and $`w^j`$ depends on the Hamming distance between them, making crossing between distant guesses less probable (this is introduced in order to encourage the exploration of a bigger region in search space).

A special feature of this approximation is the use of the candidate solutions (guesses) as their own encodings for the GA. Of course, this is made possible by the sequential and digital nature of genetic sequences, which were the very inspiration of GA and other forms of evolutionary computation. Obvious as it may seem, this is the only application we know about in which genetic algorithms are applied to genetic sequences.

### 4.3 Results of GA application

The genetic algorithm was run several times for randomly selected sequences of $`B1`$ and $`E5`$ (with the other species as target, in each case), in order to find the best values for its parameters (mutation and crossover rates, population size, etc.), for the ponderations, etc.; this was done first for each $`f_i`$, and then for the combined optimization (detailed data can be found at ). Even when a single function was optimized, we computed all the statistics on the resulting guesses, in order to see the effect of each statistics on the rest. Optimization of spectra and autocorrelation functions, for instance, do not have the same effect on the sequence, in spite of working with the same information. Optimization of $`\stackrel{~}{S}`$ causes strong oscillations in $`\stackrel{~}{\mathrm{\Gamma }}`$, whereas optimization of $`\stackrel{~}{\mathrm{\Gamma }}`$ alone tends to cause a flattening of $`\stackrel{~}{S}`$. In general, imitation of $`\stackrel{~}{\mathrm{\Gamma }}`$ is the most difficult, followed by $`\stackrel{~}{S}`$, with $`\rho _{g+c}`$, $`RSCU`$ and specially IDH being the easier. The joint optimization of the $`f_i`$ arrived at values of each $`f_i`$ only slightly worse than those obtained in single function optimization, with the exception of $`f_4`$, which was actually better. Optimization of $`\rho _{g+c}`$ and $`RSCU`$ appeared to be almost unnecessary: when only $`f_3`$, $`f_4`$ and $`f_5`$ were considered (with $`\beta ^{cu}`$ as initial condition), the final $`\rho _{g+c}`$ and $`RSCU`$ were still closer to the target species than the original sequence was to its own. In general, all $`f_i`$ are optimized by the genetic algorithm; it is even possible to make the periodicity of period 10.5 appear in sequences from which it was absent.

To remove the differences due to the amino acid sequences (which can strongly influence any coding statistic in a sample with just a few sequences), we constructed a test set with sequences encoding homologue proteins in $`B1`$ and $`E5`$. To do this, we extracted from the euGenes database  the list of homologies between these species, chose the cases with a higher identity percentage, and cut the segment of each sequence corresponding to the alignment. Thus we obtained a set $`W^B=\{w_1^B,\mathrm{},w_{20}^B\}`$ of sequences from $`B1`$, and another set $`W^E=\{w_1^E,\mathrm{},w_{20}^E\}`$ from $`E5`$, with each pair $`w_i^B,w_i^E`$ encoding very similar amino acid sequences. We performed a canonical backtranslation on $`\tau (W^E)`$, obtaining $`\beta _{B1}^{cu}(W^E)`$; we perform also a backtranslation by means of our genetic algorithm, obtaining what we will call $`\beta _{B1}^{}(W^E)`$. The computation of the diverse coding statistics allows us to see how this procedure gets the backtranslation closer to the average style of $`B1`$; moreover, since we do have $`W^B`$, we can compare with the values of that particular set of $`B1`$. For instance, for IDH, we can compute a distance between two sets of sequences $`S_1`$ and $`S_2`$ as $`d_{idh}(S_1,S_2)=|d_{ry}(S_1)d_{ry}(S_2)|+|d_{ws}(S_1)d_{ws}(S_2)|+|d_{mk}(S_1)d_{mk}(S_2)|`$. We obtain that $`d_{idh}(W^E,W^B)=0.275`$, while $`d_{idh}(\beta _{B1}^{cu}(W^E),W^B)=0.104`$, and $`d_{idh}(\beta _{B1}^{}(W^E),W^B)=0.049`$. Something similar happens with the other statistics. Graphic 6 shows the graphs of $`\stackrel{~}{S}^{11}\pi _{ry}`$ for the different sets; we can see again how $`\beta ^{}`$ builds a preimage for the image of $`W^E`$ (which is a typical $`E5`$ subset) which is far more similar to $`B1`$ and $`W^B`$ than the usual backtranslation procedure, $`\beta ^{cu}`$. For $`\stackrel{~}{\mathrm{\Gamma }}`$ the results are similar, but not so easy to observe in the graphics; instead of that, Table 3 displays the average difference between the curve $`\stackrel{~}{\mathrm{\Gamma }}(W^B)`$, and those for $`W^E`$, $`\beta _{B1}^{cu}(W^E)`$ and $`\beta _{B1}^{}(W^E)`$. Again, $`\beta ^{}`$ improves with respect to $`\beta ^{cu}`$.

Table 3: Average distance of curves $`\stackrel{~}{\mathrm{\Gamma }}`$

| Projection | $`d(W^E,W^B)`$ | $`d(\beta _{B1}^{cu}\left(W^E\right),W^B)`$ | $`d(\beta _{B1}^{}\left(W^E\right),W^B)`$ |
| --- | --- | --- | --- |
| $`\pi _{ws}`$ | 0.0018 | 0.0013 | 0.0008 |
| $`\pi _{ry}`$ | 0.0016 | 0.0019 | 0.0011 |

## 5 Discussion

The purpose of this article is to propose an improvement of the current procedures of protein backtranslation, through the inclusion of coding statistics other than RSCU which contribute to characterize the different genomes; this can be accomplished by the use of genetic algorithms. We first presented several known coding statistics, showing their idiosyncratic nature. Then we proposed a particular implementation of genetic algorithms, for a small set of coding statistics; this is only an example, since other choices of the statistics, or other implementations of evolutionary computation, may give better results. Our implementation, which is available at , does already produce backtranslations which mimic the coding statistics of the target species, in ways that are not automatically reproduced by RSCU imitation.

The definitive test for our approach would be the use of our procedure for the in vitro generation of actual artificial genes: we expect it to have a higher success frequency than the canonical backtranslation. Meanwhile, the in silico experiment consisting in the backtranslation of a human protein into “bacterial” style, and the comparison of the statistics of the resulting gene to those of an homologue bacterial gene (see section 4.3), suggest that our approach is correct. In fact, the “optimized” preimages had more exact matches with the bacterial genes (at the aligned codon positions) than the simple RSCU-based backtranslation; this happened when human proteins were optimized for “bacterial style”, and also when bacterial proteins were translated into “human”. Though small, the systematic increase in exact matches is surprising: we did not expect the imitation of coding statistics to have this effect, since the number of preimages satisfying a given profile is still huge.

This increase in exact matches suggests that the algorithm could be also applied to the problem of “gene fishing” through PCR reactions primed by degenerate primers, or “guessmers”. This is a particular case of backtranslation, limited to short sequences selected for their minimal ambiguity. Thus, coding statistics are hard to evaluate (sequences are short) and hard to optimize (sequences are rigid). In spite of these difficulties, preliminary in silico experiments seem to support this application.

Another field of application for the ideas presented here is the analysis of sequences: discussions on the relations and origins of coding statistics can be illuminated by massive backtranslation of sequences under some criteria, like we did in 4.1 with RSCU to study its relation to spectra and autocorrelation functions. Of special interest are the comparisons between genes suspected, or known, to be related by horizontal transfer. Values of RSCU and/or $`\rho _{g+c}`$ divergent from the style of a genome have been used to detect horizontally transferred genes; the degree of their divergence has been used as a clock to determine when a gene was acquired. Some authors have done this through a “reverse amelioration” which is a kind of backtranslation, and could be enriched by the results and procedures given here.

## 6 Acknowledgments

This research was supported by CONICYT through the FONDAP program in Applied Mathematics, and was started during a visit to the GREG (Group de Recherche et d’Etude sur les Genomes) at the Institut of Mathematics at Luminy, University of Marseille, France. Special thanks to Alejandro Maass for his lasting support.