Phase change for the accuracy of the median value in estimating divergence time.

We prove that for general models of random gene-order evolution of k ≥ 3 genomes, as the number of genes n goes to ∞, the median value approximates k times the divergence time if the number of rearrangements is less than cn/4 for any c <1. For some c* ≥ 1, if the number of rearrangements is greater than c*n/4, this approximation does not hold.

Introduction

The iterative improvement of approximate solutions to the Steiner tree problem by optimizing one internal vertex at a time has a substantial history in the "small phylogeny" problem for parsimony-based phylogenetics, both at the sequence level [1] and the gene order level [2]. It has been generalized to iterative local subtree optimization methods such as "tree-window-hill" [3] and "disc covering" [4,5]. Here we focus on the "median problem" for gene order where we estimate the location of a single point (the median) in a metric space given the location of the three or more points connected to the median by an edge of the tree. Given k ≥ 3 signed gene orders G1, ..., Gon a single chromosome or several chromosomes, and a metric d such as breakpoints [6], inversions [7], inversions and translocations [8], or double-cut-and-join [9], find the gene order M such that is minimized.

Although it plays a central role in gene order phylogeny, the median suffers from several liabilities. One is that it is hard to calculate in most metric spaces. Not only is it NP-hard [10], but exhaustive methods are costly for most instances, namely unless are all relatively similar to each other, which we will refer to generically as the similar genomes condition. Another problem is that heuristics tend to produce inaccurate results unless a suitable similar genomes condition holds [11]. Still another, is the tendency in some metric spaces to degenerate solutions [12] unless the same conditions prevails.

In this paper we add to this litany of difficulties by showing that as k genomes evolve over time, as modeled by any one of several biologically-motivated random walks, there is a phase change after n/4 steps, where n is the number of genes. With u < n/4 steps, the sum of the normalized distances from each of the genomes to the starting point - the ancestor - converges to ku/n in probability, and this is the median value. When u > c*n/4 steps, for a constant c* ≥ 1, the sum of the normalized distances to the median converges in probability to a value less than ku/n, and that the ancestor is no longer the median.

Our proof is inspired by a result of Berestycki and Durrett [13] in showing that the reversal distance between two signed permutations converges in probability to the actual number of steps, after rescaling, if and only if u < n/2. The technique is to construct a graph with genes as vertices and edges added between vertices according to how they are affected by transpositions. Properties of the number of components of random Erdös-Renyi graphs can then be invoked to prove the result.

Definitions

We represent a unichromosomal genome by a signed permutation, where the sign indicates whether the gene is "read" from left to right (tail-to-head) or from right to left (head to tail) on the chromosome. Let be the signed symmetric group of order n, i.e. the space of all signed permutations of length n. A reversal operation applied to a signed permutation reverses the order, and changes the signs, of one or more adjacent terms in the permutation. A DCJ operation, which can apply not only to signed permutations but to more general genomes containing linear and circular chromosomes, cuts the genome in two places and rejoins pairs of the four "loose ends" in one of two possible new ways (one of which may be equivalent to a reversal). We define the reversal and DCJ distances, dand dcj, to be the minimum number of reversal and DCJ operations, respectively, needed to transform one genome to another.

The breakpoint graph BP(Π, Π') of two genomes represented by Π and Π' contains vertices for the head and tail of each gene, black edges edges defined by the adjoining heads or tails of two adjacent genes in the genome Π and grey edges defined by two adjacent genes in the genome Π'. Let id = I, the identity permutation, and BP(Π) = BP(Π, id). It is well-known that

We need to define an orientation for grey (and black) edges of BP(Π). We traverse a cycle c ∈ cBP(Π) in a counter-clockwise manner if we start at the left-most vertex of BP(Π) (in the usual representation), travel along its unique adjacent black edge and end at the same vertex through its unique adjacent grey edge. Then we say a black edge in c is positively oriented if we move along it from left to right in a counter-clockwise traversal. Otherwise we say it is negatively oriented. Similarly, for the grey edge (i, (i + 1)) we say it is positively oriented if during a counter-clockwise traversal we move along it from ito (i + 1). Otherwise it is negatively oriented. We define the orientation function ξ on the edges of BP(Π) to be:

We say the black (grey) edges e, e' are parallel, denoted by e || e' if ξ(e) = ξ(e'). Otherwise we say they are crossing. This is just a reformulation of Hannenhalli and Pevzner's original concept of oriented cycles. An oriented cycle in this definition is a cycle including at least one positively and one negatively oriented black edge. The mechanism by which a reversal affects a genome can easily be seen using the BP graph. Let ρ be a reversal acting on two black edges e, e' in BP(Π). If they are in two different cycles we have a merger of the two to construct a new cycle. But if e, e' are in a same cycle, that cycle either splits, if e ∦ e', or does not split if e || e'.

Limit Behavior of the Median Value

Suppose dbe a metric on the space of all signed permutations length n. For a set A of these permutations, define

Then let

m(A) is called the median value of A under the metric d. A signed permutation which makes minimum is called a median solution of A. Denote by dand dcj the reversal and DCJ distances on .

Let X0 = id, the identity permutation, and let be a stochastic process on , where at random Poisson times τ, with rate 1, we choose two elements of namely i, j and let ρ(i, j) operate on that is

where ρ(i, j) is the reversal acting on i and j. We call a reversal random walk (r.w.) on. Suppose be k independent reversal r.w. all starting at the identity element, id. Define

and

We investigate the time up to which the median value of , namely m(A), remains a good estimator for the total divergence time, kt, as well as to the total distance of points in Ato id, namely . To answer this question we use the fact that the speed of escape of the r.w. up to some particular time, is the same from any point of the space and is close to 1, the maximum value. Berestycki and Durrett studied speed of transposition and reversal random walks with the related edit distances while in the latter they used "approximate reversal distance" instead of reversal itself, ignoring the effect of hurdles and fortresses. This turns out to be the same as DCJ distance on single chromosomes. We have

while

where h(π) and are the number of hurdles and fortresses, respectively.

Although Berestycki and Durrett only proved their theorem for the random transposition r.w. on S, they suggested that same method should carry over to reversal r.w. The following proposition is proved in [13] for approximate reversal distance (i.e., DCJ distance).

In this result and in the ensuing discussion ais an arbitrary sequence such that a→ ∞ as n → 0. When it is unambiguous we drop n from and .

Propostition 1 [Berestycki-Durrett] Let c be fixed and let Xt be a reversal r.w. on starting at id. Then

where

and in probability.

Remark 1 The function 1 - f is linear for c <1, f (c) = 1 - c/2, and sublinear for c >1, 1 - f (c) < c/2 This means that for c ≤ 1

and r.w. travels on an approximate geodesic (or parsimonious path) asymptotically almost surely. f is the function counting the number of tree components of an Erdös-Renyi random graph with n vertices for which the probability of having each edge is , denoted by G(c, n). See Theorem 12 in [14], Chapter V.

We extend the above theorem for the bonafide reversal distance. To do so we need to estimate the number of hurdles of Recall that an oriented cycle in a breakpoint graph is a cycle including an orientation edge, that is a grey edge with two black adjacency edges e, e', where a reversal involving e and e' splits the cycle [15]. As we discussed this is equivalent to saying e ∦ e'. It is not difficult to show

Lemma 1 Let C ∈ cBP(π), then C is oriented if and only if there exists exactly two equivalence classes of black edges, that is there exist at least two black edges with different signs.

Then

Theorem 1 Let c >0 be fixed and let X:= h(X) to be the number of hurdles in BP (X). Then

Proof. Cycles of the BP that have never been involved in a fragmentation event must be oriented, since the two rejoined black edges resulting from an inversion-induced merger of cycles cannot be parallel.

Therefore we need only to count the number of edges that have been involved in a fragmentation event. To do so we apply the method of counting cycles in [13], Theorem 3. Hurdles occur only in those cycles with length more than one that have been involved in a fragmentation up to time We call such cycles fragmented cycles. The number of fragmented cycles with length more than is always less than . But to count all fragmented cycles in with size less than we need to find an upper bound for the rate of a fragmentation up to time . Since a fragmentation occurs when two black edges in one cycle are chosen, to fragment a cycle in BP, for any chosen black edge e we only can pick another black edge e' in the same cycle whose graph distance in the breakpoint graph is less than . (The coefficient 2 arises from the fact that the cycles are alternating in BP.)

Thus the rate of fragmentation at an arbitrary time t is not more than . Integrating up to time t, this gives us the expected number of fragmented cycles at time t is . For this expectation is Now, dividing by , the result follows from Chebyshev's inequality and the fact that hurdles only occurs in fragmented cycles. ■

Theorem 2 let c >0 be fixed and let X

where f is the same function as in the statement of Proposition 1 and w' (n) is a function with in probability.

Proof. Since d(Π) = dcj(Π) + h(Π) + f˜(Π) by the proposition we have d(X) = (1 − f (c))n + w(n) + h+ f˜(X). But

in probability, by the convergence of and in Proposition 1 and Theorem 2 and

Theorem 3 Let be k independent reversal r.w in starting at id. Suppose either

a) d := dcj dcj distance

or

b) reversal distance.

Then for c

Proof. We prove the theorem only for d. The proof of the DCJ case is similar. For all i, j ∈ {1, ..., k} and for a median solution x of

Therefore,

We conclude

Let . Then by Theorem 2 we have for all i, j i ≠ j

and

where in probability. Thus

Then

for a constant k'. Also . Therefore, there exists a constant k* such that

This implies

This proves the theorem. ■

Remark 2 The statement of the theorem suggests ignoring the error of order for a→ ∞. id remains as the median of leaves of k independent stochastic processes up to time asymptotically almost surely.

Theorem 4 Let be fixed. Suppose d is either DCJ or reversal distance. Then by the hypothesis of Theorem 3

Proof. This follows directly from the fact that

in probability. ■

Now, it is natural to ask whether the statement of Theorem 4 also holds for some time after . In other words, is the median value kcn a fair estimator for the total time of divergence? We conjecture not, that the property is lost after time , but for now can only prove a weaker upper bound for this time.

Theorem 5 Let be fixed. Suppose d is either DCJ or reversal distance. Then by the same hypothesis as in Theorem 3

where

is strictly positive for

Remark 3 This theorem shows after time the error is of order n and so the median value is not a good estimate of k times the divergence time.

Proof.

where in probability. Dividing by n, the result follows. ■

In fact, since f (c), c >0 is decreasing and for c <1, , it is easy to see that in the case k = 3, for c >0.75, is of order for some .

Theorem 6 Let k = 3 and d be either dcj or dr. Consider the same hypothesis in Theorem 3. Assume c* be solution of

Then for all c > c* there exists such that

Proof.

Computing for i = 2, 3 is the same as . This is true since the Cayley graph of w.r.t. reversals is symmetric and regular and so . But therefore by symmetry of the Cayley graph we can just consider . Hence,

Let x >0 be so that

This means

So it suffices to prove above inequality for x = c >c∗. Since f (x) >0 for all x >0

in which the right hand side is strictly increasing, Therefore for all c ≥ c*

This proves the statement. ■

Now, we would like to measure the volume of that part of the space for which median does well, compared with the whole space. The ratio of the two converges to 0 as n goes to ∞, showing that the median is only useful in a highly restricted region of the space.. The following theorem is entailed by a theorem in [16]. Let c= c(Π) be the number of cycles in the BP graph of a random . Let dbe a distance (metric) on . Define

to be the ball of radius cn in .

Theorem 7 Let 0 < c <1 be fixed. Then

Proof.

a)

Suppose γdoes not converge to 0. Therefore there exists a subsequence such that for a constant ε >0. This implies

But by Theorem 2.2 in [16], we have

That is in contradiction with the above inequality since

b) For the second part it suffices to observe that for all we have

Therefore

and the result follows part (a) since

■

Conclusion

We have shown that the median value for DCJ and for reversal distance for a reversal r.w.has good limiting properties if the number of steps remains below cn/4, for any c <1, but for some value c > 1, more than this number of steps destroys these limiting properties. The critical value may indeed be c = 1, but for now we can only show that for c > 3 (and c > 2) the median value is no longer a good estimator of the distance between the id and the current position of the r.w. (and k times the divergence time, respectively).

Note that a simulation strategy to estimate c is not available because of the hardness of calculating the median. As n increases even to moderate values all exact methods require prohibitive computing time.

These results imply that the steinerization strategy for the small phylogeny problem may lead to poor estimates of the interior nodes of a phylogeny unless the taxon sampling is sufficient to assure that a "similar genomes condition" holds for every k-tuple of genomes used in the course of of the iterative optimization search. This can be monitored prior to each step in the iterative optimization of the phylogeny through successive application of the median method.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

AJ and DS planned the study, carried out the research and wrote the article.