Skip to main content

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

Abstract

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, ..., G k on 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 i = 1 k d ( G i , M ) 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 G 1 , G k 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 k d / n 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 S n ± 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, d r and 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

d c j ( Π ) = n + 1 - | c B P ( Π ) | .
(1)

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 t , (i + 1) h ) we say it is positively oriented if during a counter-clockwise traversal we move along it from i t to (i + 1) h . Otherwise it is negatively oriented. We define the orientation function ξ on the edges of BP(Π) to be:

ξ ( e ) = + 1 if  e is positively oriented - 1 if  e is negatively oriented .
(2)

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 dn be a metric on the space of all signed permutations length n. For a set A of these permutations, define

g A d , n : S n ± 0 = { 0 } ,
(3)
g A d , n ( x ) : = y A d n ( x , y ) .
(4)

Then let

m d , n ( A ) : = m i n { g A d , n ( x ) : x S n ± } .
(5)

md,n(A) is called the median value of A under the metric dn. A signed permutation which makes g A d , n minimum is called a median solution of A. Denote by d r and dcj the reversal and DCJ distances on S n ± .

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

X τ κ + n = X τ κ n o ρ ( i , j ) ,
(6)

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

A t ( n ) : = { X t 1 , n , , X t k , n }
(7)

and

ε t d , n : = g A t d , n ( i d ) - m d , n ( A t ) .
(8)

We investigate the time up to which the median value of X t 1 , n ,, X t k , n , namely md,n (A t ), remains a good estimator for the total divergence time, kt, as well as to the total distance of points in A t to id, namely g A t d , n ( i d ) . 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

d r ( π , I ) = n + 1 - c ( π ) + h ( π ) + f ̃ ( π )
(9)

while

d c j ( π , I ) = n + 1 - c ( π ) ,
(10)

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

Although Berestycki and Durrett only proved their theorem for the random transposition r.w. on S n , 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 a n is an arbitrary sequence such that a n → ∞ as n → 0. When it is unambiguous we drop n from A t ( n ) and X t n .

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

d c j ( i d X c n /2 ) = (1 - f ( c )) n + w ( n ),
(11)

where

f ( c ) : = 1 c k = 1 k k - 2 k ! ( c e - c ) k
(12)

and w ( n ) a n n 0in 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

d c j ( i d , X c n / 2 ) - c n 2 = w ( n )
(13)

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 c n , 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 X c n 2 . 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 t be a reversal r.w.starting at id. Define h t := h(X t ) to be the number of hurdles in BP (X t ). Then

h c n / 2 a n n 0 i n p r o b a b i l i t y .
(14)

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 c n 2 . We call such cycles fragmented cycles. The number of fragmented cycles with length more than n is always less than n . But to count all fragmented cycles in X c n 2 with size less than n we need to find an upper bound for the rate of a fragmentation up to time c n 2 . 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 2 n . (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 n n 2 ( n ) n = 2 n . Integrating up to time t, this gives us the expected number of fragmented cycles at time t is 2 n ·t. For t= c n 2 this expectation is c n . Now, dividing by a n n , 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 t be a reversal r.w. on S n ± starting at id and let d r := d r ( n ) denote the reversal distance on S n ± . Then

d r ( i d , X c n / 2 ) = ( 1 - f ( c ) ) n + w ( n )
(15)

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

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

w ( n ) a n n : = w ( n ) + h c n / 2 + f ̃ ( X c n / 2 ) a n n 0
(16)

in probability, by the convergence of w ( n ) a n n and h c n / 2 a n n in Proposition 1 and Theorem 2 and f ~ X c n / 2 1 .

Theorem 3 Let X t 1 , n ,, X t k , n be k independent reversal r.w in S n ± starting at id. Suppose either

a) d := dcj dcj distance

or

b) d:= d r ( n ) reversal distance.

Then for c < 1 4 we have ε c n d , n a n n 0in probability.

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

d r ( n ) ( X t i , n , X t j , n ) d r ( n ) ( x , X t i , n ) + d r ( n ) ( x , X t j , n ) .
(17)

Therefore,

i j d r ( n ) ( X t i , n , X t j , n ) i j ( d r ( n ) ( x , X t i , n ) + d r ( n ) ( x , X t j , n ) ) .
(18)

We conclude

d r n ( X t i , n , X t j , n ) ( k - 1 ) m d , n ( A t ( n ) ) ( k - 1 ) g A t ( n ) ( i d ) .
(19)

Let c 1 4 . Then by Theorem 2 we have for all i, j ij

d r ( n ) ( X c n i , n , X c n j , n ) = 2 c n - w ( n )
(20)

and

d r ( n ) ( i d , X c n i , n ) = c n - w ( n )
(21)

where w ( n ) ( a n n ) 0 in probability. Thus

k 2 ( 2 c n - w ( n ) ) ( k - 1 ) m d , n ( A c n ( n ) ) ( k - 1 ) k ( c n - w ( n ) ) .
(22)

Then

| m d , n ( A c n ( n ) ) - k c n | k w ( n )
(23)

for a constant k'. Also | g A c n ( n ) ( i d ) -kcn|kw ( n ) . Therefore, there exists a constant k* such that

| m d , n ( A c n ( n ) ) - g A c n ( n ) d , n ( i d ) | k * w ( n ) .
(24)

This implies

ε c n a n n = m d , n ( A c n ( n ) ) - g A c n ( n ) d , n ( i d ) a n n 0 i n p r o b a b i l i t y .
(25)

This proves the theorem. ■

Remark 2 The statement of the theorem suggests ignoring the error of order o ( a n n ) for a n → ∞. id remains as the median of leaves of k independent stochastic processes X t 1 , n ,, X t k , n up to time n 4 asymptotically almost surely.

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

k c n - m d , n ( A c n ) a n n 0 i n p r o b a b i l i t y a s n .
(26)

Proof. This follows directly from the fact that

k c n - g A c n d , n ( i d ) a n n 0
(27)

in probability. ■

Now, it is natural to ask whether the statement of Theorem 4 also holds for some time after n 4 . 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 n 4 , but for now can only prove a weaker upper bound for this time.

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

k c n - m d , n ( A c n ) n α c
(28)

where

α c : = k ( 1 - f ( 2 c ) )
(29)

is strictly positive for c> 1 2

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

Proof.

k c n - m d , n ( A c n ) k c n - g A c n d , n ( i d ) = k ( 1 - f ( 2 c ) ) n + w ( n ) ,
(30)

where w ( n ) a n n 0 in probability. Dividing by n, the result follows. ■

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

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

f ( x 2 ) = 1 3 .
(31)

Then for all c > c* there exists β c d such that

ε c n 4 d , n = o ( β c d n ) .
(32)

Proof.

m d , n ( A c n 4 ) d ( X c n 4 1 , n , X c n 4 2 , n ) + d ( X c n 4 1 , n , X c n 4 3 , n ) .
(33)

Computing d ( X c n 4 1 , n , X c n 4 i , n ) for i = 2, 3 is the same as d ( i d , X c n 2 1 , n ) . This is true since the Cayley graph of S n ± w.r.t. reversals is symmetric and regular and so P ( X 0 = i d , X c n 4 = Π ) =P ( X 0 = Π , X c n 4 = i d ) . But therefore by symmetry of the Cayley graph we can just consider d ( i d , X c n 2 1 , n ) . Hence,

m d , n ( A c n 4 ) 2 ( 1 - f ( c ) ) n + 2 w ( n ) .
(34)

Let x > 0 be so that

0 < - 2 ( 1 - f ( x ) ) + 3 ( 1 - f ( x 2 ) ) .
(35)

This means

g A c n 4 d , n ( i d ) > m d , n ( A c n 4 ) .
(36)

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

1 + 2 f ( x ) - 3 f ( x 2 ) > 1 - 3 f ( x 2 )
(37)

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

1 + 2 f ( c ) - 3 f ( c 2 ) > 1 - 3 f ( c * 2 ) = 0 .
(38)

This proves the statement. ■

Now, we would like to measure the volume of that part of the space S n ± 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 n = c n (Π) be the number of cycles in the BP graph of a random Π S n ± . Let d n be a distance (metric) on S n ± . Define

B c n d = B c n d , n : = { Π S n ± , d ( Π , i d ) c n }
(39)

to be the ball of radius cn in S n ± .

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

a ) γ n = | B c n d c j | | S n ± | 0 a s n ,
(40)
b ) γ ' n = | B c n d r | | S n ± | 0 a s n .
(41)

Proof.

  1. a)
    F o r a l l Π B c n d c j , | c B P ( Π ) | ( 1 - c ) n .
    (42)

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

E ( c n i ) ε ( 1 - c ) n i .
(43)

But by Theorem 2.2 in [16], we have

E ( c n i ) n i 0 a s n i .
(44)

That is in contradiction with the above inequality since

ε ( 1 - c ) n i n i ε ( 1 - c ) > 0 .
(45)
  1. b)

    For the second part it suffices to observe that for all Π S n ± we have

    d r ( Π ) d c j ( Π ) .
    (46)

Therefore

B d r ( Π ) B d c j ( Π )
(47)

and the result follows part (a) since

γ n γ n 0 a s n .
(48)

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.

References

  1. Sankoff D, Cedergren RJ, Lapalme G: Frequency of insertion/deletion, transversion and transition in the evolution of 5S ribosomal RNA. Journal of Molecular Evolution. 1976, 7: 133-149. 10.1007/BF01732471.

    Article  CAS  PubMed  Google Scholar 

  2. Blanchette M, Bourque G, Sankoff D: Breakpoint phylogenies. Genome Informatics. Edited by: S. Miyano & T. Takagi. 1997, Tokyo: Universal Academy Press, 25-34.

    Google Scholar 

  3. Sankoff D, Abel Y, Hein J: A Tree - A Window - A Hill; Generalization of nearest-neighbour inter-change in phylogenetic optimisation. Journal of Classification. 1994, 11: 209-232. 10.1007/BF01195680.

    Article  Google Scholar 

  4. Huson D, Nettles S, Warnow T: Disk-covering, a fast-converging method for phylogenetic tree reconstruction. Journal of Computational Biology. 1999, 6: 369-386. 10.1089/106652799318337.

    Article  CAS  PubMed  Google Scholar 

  5. Tang J, Moret B: Scaling up accurate phylogenetic reconstruction from gene-order data. Bioinformatics. 2003, 19: i305-i312. 10.1093/bioinformatics/btg1042.

    Article  PubMed  Google Scholar 

  6. Sankoff D, Blanchette M: The median problem for breakpoints in comparative genomics. Proceedings of Computing and Combinatorics (COCOON). Edited by: T. Jiang and D.T. Lee. 1997, 1276: 251-263. 10.1007/BFb0045092. Lecture Notes in Computer Science

    Chapter  Google Scholar 

  7. Sankoff D, Sundaram G, Kececioglu J: Steiner points in the space of genome rearrangements. International Journal of the Foundations of Computer Science. 1996, 7: 1-9. 10.1142/S0129054196000026.

    Article  Google Scholar 

  8. Bourque G, Pevzner PA: Genome-scale evolution: Reconstructing gene orders in the ancestral species. Genome Research. 2002, 12: 26-36.

    PubMed Central  CAS  PubMed  Google Scholar 

  9. Zhang M, Arndt W, Tang J: An exact solver for the DCJ median problem. Pacific Symposium on Biocomputing. 2009, 138-149.

    Google Scholar 

  10. Tannier E, Zheng C, Sankoff D: Multichromosomal median and halving problems under different genomic distances. BMC Bioinformatics. 2009, 10: 120-10.1186/1471-2105-10-120.

    Article  PubMed Central  PubMed  Google Scholar 

  11. Zheng C, Sankoff D: On the Pathgroups approach to rapid small phylogeny. BMC Bioinformatics. 2011, 12: S4-

    PubMed Central  PubMed  Google Scholar 

  12. Haghighi M, Sankoff D: Medians seek the corners, and other conjectures. BMC Bioinformatics. 2012, 13 (S19): S5-

    PubMed Central  PubMed  Google Scholar 

  13. Berestycki N, Durrett R: A phase transition in the random transposition random walk. Probability Theory and Related Fields. 2006, 136: 203-233. 10.1007/s00440-005-0479-7.

    Article  Google Scholar 

  14. Bollobás B: Random Graphs. 2001, Cambridge University Press, 2

    Book  Google Scholar 

  15. Hannenhalli S, Pevzner PA: Transforming cabbage into turnip: Polynomial algorithm for sorting signed permutations by reversals. Journal of the ACM. 1999, 46: 1-27. 10.1145/300515.300516.

    Article  Google Scholar 

  16. Székely LA, Yang Y: On the expectation and variance of the reversal distance. Acta Univ. Sapientiae, Mathematica. 2009, 1: 5-20.

    Google Scholar 

Download references

Acknowledgements

Research supported in part by grants from the Natural Sciences and Engineering Research Council of Canada. DS holds the Canada Research Chair in Mathematical Genomics. Thanks to Armin Jamshidpey and Leili Rafiee Sevyeri for help in preparation of the manuscript.

Declarations

Publication of this article was supported by the Canada Research Chair in Mathematical Genomics.

This article has been published as part of BMC Bioinformatics Volume 14 Supplement 15, 2013: Proceedings from the Eleventh Annual Research in Computational Molecular Biology (RECOMB) Satellite Workshop on Comparative Genomics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S15.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Arash Jamshidpey.

Additional information

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.

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Jamshidpey, A., Sankoff, D. Phase change for the accuracy of the median value in estimating divergence time. BMC Bioinformatics 14 (Suppl 15), S7 (2013). https://doi.org/10.1186/1471-2105-14-S15-S7

Download citation

  • Published:

  • DOI: https://doi.org/10.1186/1471-2105-14-S15-S7

Keywords