Estimating true evolutionary distances under rearrangements, duplications, and losses

Background The rapidly increasing availability of whole-genome sequences has enabled the study of whole-genome evolution. Evolutionary mechanisms based on genome rearrangements have attracted much attention and given rise to many models; somewhat independently, the mechanisms of gene duplication and loss have seen much work. However, the two are not independent and thus require a unified treatment, which remains missing to date. Moreover, existing rearrangement models do not fit the dichotomy between most prokaryotic genomes (one circular chromosome) and most eukaryotic genomes (multiple linear chromosomes). Results To handle rearrangements, gene duplications and losses, we propose a new evolutionary model and the corresponding method for estimating true evolutionary distance. Our model, inspired from the DCJ model, is simple and the first to respect the prokaryotic/eukaryotic structural dichotomy. Experimental results on a wide variety of genome structures demonstrate the very high accuracy and robustness of our distance estimator. Conclusion We give the first robust, statistically based, estimate of genomic pairwise distances based on rearrangements, duplications and losses, under a model that respects the structural dichotomy between prokaryotic and eukaryotic genomes. Accurate and robust estimates in true evolutionary distances should translate into much better phylogenetic reconstructions as well as more accurate genomic alignments, while our new model of genome rearrangements provides another refinement in simplicity and verisimilitude.


Introduction
Interest in the evolution of genome structure has been growing steadily in the last 10 years, sustained in part by the ever increasing number of sequenced genomes. In particular much work has been done on rearrangements (see, e.g., [1]), using the convention that each chromosome of the genome is represented by an ordered list of identifiers, each identifier referring to a syntenic block or, more commonly, to a member of a gene family. (In the following, we shall use the word "gene" in a broad sense to denote elements of such orderings and refer to such orderings as "gene orders".) Variations in the placement of homologous genes, as well as variations in gene content and multiplicity, among organisms can then be analyzed. Such data is of great interest to evolutionary biologists, but also to comparative genomicists and to any researcher interested in understanding evolutionary changes in pathogens, crop plants, and, more generally, the biome.
The most fundamental task in the analysis of such data is to estimate the amount of evolutionary change between two genomes-that is, to compute a pairwise evolutionary distance. The true distance, that is, the number of actual evolutionary events (rearrangements, duplications, and losses) that took place during the course of evolution, is what we want to obtain, but is not, of course, something that we can compute. Researchers have thus used a two-stage process, in which a well defined measure is first computed (such as an edit distance, that is, the smallest number of evolutionary events needed to transform one genome into the other), then a statistical model of evolution is used to infer an estimate of the true distance by deriving the effect of a given number of changes in the model on the computed measure and (algebraically or numerically) inverting the derivation to produce a maximum-likelihood estimate of the true distance under the model. This second step is usually called a distance correction and has long been used for sequence (DNA) data (see, e.g., [2]) as well as, more recently, for geneorder data, for which see [3][4][5][6][7].
Evolutionary events that affect the gene order of genomes include various rearrangements, which affect only the order, and gene duplications and losses, which affect both the gene content and, indirectly, the order. (Gene insertion, corresponding to lateral gene transfer or neofunctionalization, can be viewed as a special case of duplication.) Rearrangements themselves include inversion, transposition, block exchange, circularization and linearization, all of which act on a single chromosome, and translocation, fusion, and fission, which act on two chromosomes. All of these operations are subsumed in the double-cut-and-join (DCJ) [8,9], which has formed the basis for much of the algorithmic research on rearrangements over the last few years, including a statistically based method to estimate the true evolutionary distance between two genomes [7]. DCJ makes two cuts, which can be in the same chromosome or in two different chromosomes, producing four cut ends, then rejoins the four cut ends in any of the three possible ways. The DCJ model, however, is unrealistic in two major respects. First, if the two cuts are in the same chromosome, one of the two nontrivial rejoinings causes a fission, creating a new circular chromosome. However, circular chromosomes do not normally arise in organisms with linear chromosomes, and prokaryotic genomes normally consist of a single circular chromosome. Nor can this form of rejoining be forbidden as, without it, DCJ simply reduces to inversion. Secondly, DCJ is a model of rearrangements: it does not take into account evolutionary events that alter the gene content, such as duplications and losses.
Of these two problems, the first has not been seriously addressed: the model we present here is, to the best of our knowledge, the first model that naturally preserves the dichotomy between prokaryotic and eukaryotic genomes. While gene (or segment) duplications and losses have long been studied by geneticists and molecular biologist, their integration with rearrangements in a unified model has seen relatively little work to date. El-Mabrouk [10] gave an exact algorithm to compute edit distances for inversions and losses and also a heuristic to approximate edit distances for inversions, losses, and nonduplicating insertions (all of her results assume that genes cannot be duplicated). More recently, Yancopoulos and Friedberg [11] gave an algorithm to compute edit distances under deletions, insertions, duplications, and DCJ operations, under the constraint that each deletion can only remove a single gene. These and other approaches targeted the edit distance, not the true evolutionary distance. Swenson et al. [12] gave an algorithm to approximate the true evolutionary distance under deletions, insertions, duplications, and inversions for unichromosomal genomes and showed good results under simulations and for small-scale phylogenetic reconstruction. Rearrangements, duplications and losses have also been addressed in the framework of ancestral reconstruction (see, e.g., [13]). All of these approaches have focused on parsimony criteria and have used preassigned weights for the various operations.
In this paper, we propose a new evolutionary model which respects the dichotomy between prokaryotic and eukaryotic genomes and which takes gene duplications and losses into account. Using this new evolutionary model, we develop a statistically based method to estimate the true evolutionary distance in terms of the actual number of rearrangements, gene duplications, and gene losses. Finally, we provide extensive experimental results on a wide variety of genome structures to illustrate the robustness and high accuracy of our estimator.

Genomes as gene-order data
We denote the tail of a gene g by g t and its head by g h . We write +g to indicate an orientation from tail to head (g t g h ), -g otherwise (g h g t ). Two consecutive genes a and b can be connected by one adjacency of one of the following four types: If gene c lies at one end of a linear chromosome, then we also have a singleton set, {c t } or {c h }, called a telomere. A genome can then be represented as a multiset of genes together with a multiset of adjacencies and telomeres. For example, the toy genome in Figure 1, composed of one linear chromosome, (+a, +b, -c, +a, +b, -d, +a), and one circular one, (+e, -f), can be represented by the multiset of genes {a, a, a, b, b, c, d, e, f} and the multiset of adjacencies and telomeres {{a t }, Because of the duplicated genes, there is no one-to-one correspondence between genomes and multisets of genes, adjacencies, and telomeres. For example, the genome composed of one linear chromosome, (+a, +b, -d, +a, +b, -c, +a) and one circular one (+e, -f) would have the same multisets of genes, adjacencies and telomeres as that in Figure 1.
Preliminaries on the evolutionary model We use two parameters: the probability of occurrence of a gene duplication, p d , and the probability of occurrence of a gene loss, p l ; the probability of occurrence of a rearrangement is then just p r = 1 -p d -p l . The next event is chosen from the three categories according to these parameters.
For rearrangements, we will select two adjacencies or telomeres with replacement uniformly from the multiset of all adjacencies and telomeres and then decide which rearrangement event we apply. The four cases are as follows.
Select two different adjacencies, or one adjacency and one telomere, in the same chromosome .. a i-1 a i ... a j a j +1 ... a n ). Reversing all genes between a i and a j yields (a 1 ... a i-1 -a j ... -a i a j+1 ... a n ). from one linear chromosome C = (a 1 ... a i a i+1 ... a n ) and { , from another linear chromosome D = (b 1 ... b j b j+1 ... b n ). Now exchange the two segments between these two chromosomes C and D. There are two possible outcomes, (a 1 ... a i b j+1 ... b n ) and (b 1 ... b j a i+1 ... a n ) or (a 1 ... 1 . This operation causes a translocation (or, if at least one chromosome is circular, a fusion).

Select the same adjacency twice
twice from linear chromosome C = (a 1 ... a i a i+1 ... a n ). Then split C into two new linear chromosomes, (a 1 ... a i ) and (a i+1 ... a n ). For gene duplication, we uniformly select a position to start duplicating a short segment of chromosomal material and place the new copy to a new position within the genome. We set L max as the maximum number of genes in the duplicated segment and assume that the number of genes in that segment is a uniform random number between 1 and L max . For example, select one segment a i+1 ... a i+L to duplicate and insert the copy between one adjacency . Such an operation duplicates L genes and L -1 adjacencies, removes one adjacency, and adds two new adjacencies; thus genes a i+1 ,..., a i+L-1 and a i+L are added to the multiset of genes, the adjacency is removed, and L + 1 new adjacencies, For gene loss, we restrict deletions to genes with at least two copies in the genome and we delete one gene at a time. We uniformly select one gene from the set of all candidate genes and delete it.

Methods
An overview of our technique for estimating the true evolutionary distance The problem of estimating the true evolutionary distance is defined as follows: Input: The original genome G and the final genome F.
Output: An estimate of the actual number of evolutionary events that took place in the evolutionary history to transform G into F.
Based on the multisets of genes and of adjacencies and telomeres of G, for any genome G* of N* genes and l* linear chromosomes, we can build the vector

… …
, where C is the upper bound for the number of copies of one gene, NG i * (i = 1,..., C) is the number of genes with exactly i copies in the genome G*, SA i * (i = 1,..., C) is the number of adjacencies with exactly i copies in G* that also appear in G, DA* is the number of adjacencies in G* that do not appear in G, ST * is the number of telomeres in G* that also appear in G, and DT* is the number of telomeres in G* that do not appear in G. We can write Let G k be the genome obtained from G = G 0 by applying k randomly selected evolutionary operations-under our model, the (i + 1)st evolutionary operation is selected from all possible rearrangements, gene duplications, and gene losses on genome G i according to the parameters p d and p l . We can compute the vector In the section, we show that, given G, we can also produce for the expected vector E(V G (G k )), for any integer k > 0. Our approach for estimating the true evolutionary distance is then to return the integer k that minimizes the 1-norm distance between V G G k ( ) and V G (F).
Estimation of the expected vector after some number of random evolutionary events Given the original genome G, the complete vector for genome G k is defined as , where NG i k is the number of genes with exactly i copies in the genome G k , SA i k (shared adjacencies) is the number of adjacencies with exactly i copies in G k that also appear in G, DA k (distinct adjacencies) is the number of adjacencies in G k that do not appear in G, ST k (shared telomeres) is the number of telomeres in G k that also appear in G, and DT k (distinct telomeres) is the number of telomeres in G k that do not appear in G.
Assume the original genome G has N genes, where each gene has at most C = O(1) copies, and l linear chromosomes, with l = O(1). We thus ignore items NG i k and SA i k for (i >C). The initial vector V G (G 0 ) is then

Rearrangements
We select two adjacencies or telomeres uniformly with replacement, from the multiset of all adjacencies or telomeres.
Theorem 1 Assume all genomes have O(1) linear chromosomes, each gene has at most C = O(1) copies, and Proof In our evolutionary model, each rearrangement operation replaces old adjacencies or telomeres with new ones. Obviously, any rearrangement operation will not change the gene content, so NG i k+1 (i = 1,2,..., C) will be the same. http://www.biomedcentral.com/1471-2105/11/S1/S54 We first ignore the adjacencies or telomeres in the original genome G created after a rearrangement event. Remember two adjacencies or telomeres are selected with replacement uniformly from the multiset of all adjacencies and telomeres, and the number of all adjacencies and telomeres for genome G k is (N k +l k ). Consider the multi-set A i of SA i k adjacencies with exactly i copies in G k that also appear in G.
The probability that exactly one of the two selected , the probability that equivalent adjacencies , and the probability that some adjacency from A i is selected twice is for each item. Consider any adjacency (a, b) in G: we might recover it only if we select two adjacencies or telomeres containing two genes a and b. Since each gene has at most C copies in the genome, there are at most C 2 pairs of adjacencies or telomeres that may lead to recovery of the adjacency (a, b). So, with probability at most

Gene duplication
We select uniformly at random an integer between 1 and L max (the maximum number of genes in the duplicated segment), then select uniformly at random a position in the genome where to start the duplication, then insert the copy at another position selected uniformly at random.
Theorem 2 Assume all genomes have O(1) linear chromosomes, each gene has at most C = O(1) copies, no two duplicate genes or adjacencies are within the segment to be duplicated, and . The probability that the placement of the duplicated segment breaks one adjacency at any specific site is 1/(N k + l k ).
We again first ignore the adjacencies or telomeres in the original genome G created after a duplication event.
Since we assume that no two genes or adjacencies are the same within the duplicated segment, we have

Gene loss
We uniformly select one gene with at least two copies and delete it. Proof In our model of gene loss, one gene with at least two copies is uniformly selected. The number of all possible genes to be deleted is N NG ) > 1 genes with exactly i copies in G k , the probability that one of them is selected and deleted is , the number of genes with exactly i copies decreases by i and the number of genes with exactly (i -1) copies increases by (i -1).
We ignore the adjacencies or telomeres in the original genome G to be created after one gene loss. For SA i k (i > 2) adjacencies with exactly i copies in G k which also appears in G, it is difficult to compute the number f i (del j ) of such adjacencies that each single deletion del j (j = 1,..., N k -NG k 1 ) would affect. But we know that each adjacency with exactly i (i > 2) copies must relate to two genes with more than 2 copies, so we have For SA k 1 adjacencies with exactly 1 copy in G k that also appears in G, it is also difficult to compute the number f 1 (del j ) of such adjacencies that each single deletion ( ) ) = = − ∑ is the count of genes with at least two copies but related to those adjacencies with exactly 1 copy in G k that also appear in G. We consider the effect of rearrangements, gene duplications and losses, and we approximate as follows: Finally, we also approximate the number of adjacencies RSA k+1 that we could thus ignore under rearrangements, gene duplications, and gene losses, and distribute it to the correction of SA i k as follows: Now, given G 0 , we estimate E(V G (G k )) for k > 0 by iterating k times the above formulas (using with p d and p l ); at every step we identify E(V G (G k-1 )) with the actual vector V G (G k-1 ).

Results and discussion
We now present experimental results on the accuracy of our estimation of the expected vector after a given number of random evolutionary events and on the quality of our estimator for the true evolutionary distance (in terms of the actual number of evolutionary events). Our experiments all start with one genome with no duplicated genes and some chosen number of linear and circular chromosomes of various sizes. We first apply some number (usually 10) of duplication events (L max = 10 in all cases) to generate the original genome G with some initial duplicated genes. Then this genome is subjected to a prescribed number k of evolutionary events chosen according to p d and p l to obtain a final genome G k . We vary k from 0 to twice the number of genes. We ran tests on any types of initial genomes designed to resemble actual organismal genomes; we tested different choices of parameters on different genomes; and in each case we generated 10,000 runs to obtain a tight estimate of variance.
We compute the vector representations for all intermediate genomes and then use our method to estimate the evolutionary distance. Due to space limitations, we present results on just three initial genomes: 25,000 genes and 25 linear chromosomes (p d = 0.05, p l = 0.15); 10, 000 genes and 5 linear chromosomes (p d = 0.1, p l = 0.2); and 1, 000 genes and 1 circular chromosome (p d = 0.2, p l = 0.6). The first two examples match large and smaller metazoan genomes, the last matches a small bacterial genome.

Accuracy of the expected vector after k random evolutionary events
We study the behavior of our estimator V G G k ( ) by comparing its prediction to the sample mean for V G (G k ), as computed from our 10,000 trials. In all of our experiments, we find that V G G k ( ) is very close to the sample mean for V G (G k ). Figure 2 shows the values in the vector as a function of the actual number of evolutionary events. SA k 3 and NA k 3 represent the number of adjacencies and genes with at least 3 copies in the original genome G, respectively. The figure shows that our estimation and the sample mean for V G (G k ) are always very close.
Accuracy of the estimation of the actual number of evolutionary events We want to study the accuracy of our estimator for the actual number of evolutionary events; in order to do that, we create simulations with controlled numbers of evolutionary events and set up a threshold for correction in the estimation procedure. Specifically, we vary the actual number of evolutionary events from 0 to twice the number of genes in the original genome and we set 4 times the number of genes as an upper limit on the maximum number of evolutionary events. C is set to 10. Thus our estimated number k is chosen to minimize | V G G k ( ) -V G (F)| 1 , the 1-norm distance between V G G k ( ) and V G (F). Figure 3 shows the mean and standard deviation for the actual number of evolutionary events estimated by our approach. Our approach provides accurate estimates, with very small variance.
We also study the mean absolute difference between the actual number of evolutionary events and our estimator, shown in Figure 4. Table 1 shows that the estimates are quite accurate up to very large numbers of events. Rearrangements, gene duplications, and gene losses fall under the category of "rare genomic events" (in the terminology of [14]), yet our estimator works well even for numbers that would instead indicate common events.

Robustness to unknown model parameters
Up to now we have fixed p d and p l . We now consider the case in which these parameters are unknown-clearly the more common case in practice. We generate 10,000 cases with randomly chosen parameters p d and p l (at 1% resolution, p d < 4p l ) and with actual numbers of evolutionary events varying from 0 to twice the number of genes, setting an upper limit of 4 times the number genes for the maximum number of evolutionary events.
Given the original genome, our estimated vector V G

Conclusion
We propose a new evolutionary model for rearrangements, gene duplications and losses, and a corresponding method for estimating true evolutionary distance. The model is, to our knowledge, the first to preserve the structural dichotomy in genomic organization between most prokaryotes and most eukaryotes, and one of the few to unite rearrangements, duplications, and losses.