Genomes and rearrangement operations

We model the evolutionary rearrangement of a genome containing n distinct signed genes through the accumulated operation of number of processes familiar in classical genetics: inversion, reciprocal translocation, transposition, chromosome fusion and fission, operating on linear chromosomes. We will not delve into the details of the operations; formally they can all be subsumed under a single operation called double-cut-and-join (DCJ) which need not be described here. All that is needed for our purposes is a formula due to Yancopoulos et al.[10] that gives d(G₁, G₂), the minimum number of rearrangement operations needed to transform one genome into another in terms of properties of the “breakpoint graph” determined by G₁ and G₂, the initial and final genomes. To calculate D efficiently, we construct and analyze the breakpoint graph as follows.

For each genome, each gene g with a positive polarity is replaced by two vertices representing its two ends, i.e., by a “tail” vertex and a “head” vertex in the order g _t , g _h ; for –g we would put g _h , g _t . Each pair of successive genes in the gene order defines an adjacency, namely the pair of vertices that are adjacent in the vertex order thus induced. For example, if i, j, –k are three neighbouring genes on a chromosome then the unordered pairs {i _h , j _t } and {j _h , k _h } are the two adjacencies they define.

If there are m genes on a chromosome, it has determined 2m vertices by this stage. The first and the last of these vertices are called telomeres. We convert all the telomeres in genome G₁ and G₂ into adjacencies with additional vertices all labelled T₁ or T₂, respectively. The breakpoint graph has a blue edge connecting the vertices in each adjacency in G₁ and a red edge for each adjacency in G₂. We make a cycle of any path ending in two T₁ or two T₂ vertices, connecting them by a red or blue edge, respectively, while for a path ending in a T₁ and a T₂, we collapse them to a single vertex denoted “T”.

Each vertex is now incident to exactly one blue and one red edge. This bicoloured graph decomposes uniquely into κ alternating cycles. If n′ is the number of blue edges, then [10]:

d(G₁, G₂) = n′ – κ. (1)

For a given unrooted binary tree T on N given genomes G₁, G₂, ⋯, G _N (and thus with N – 2 unknown ancestral genomes M₁, M₂, ⋯, M _{N –2} and 2N – 3 branches) as depicted in Fig. 1, the small phylogeny problem is to infer the ancestral genomes so that the total edge length of T, namely ∑ _{XY∈E ( T )}d(X, Y), is minimal.

The computational complexity of the median problem, which is just the small phylogeny problem with N = 1, is known to be NP-hard and hence so is that of the general small phylogeny problem. Our method will be shown to run in linear time, so obviously it is not guaranteed to find an exact solution. One of the goals of this paper is to determine for what range of instances PATHGROUPS leads to accurate solutions, and how the approach may be improved to extend this range.

Paths and fragments

We generalize our definition of a path to be any connected subgraph of a breakpoint graph, namely any connected part of a cycle. Initially, each blue edge in the given genomes is a path.

A fragment is any set of genes connected by red edges in a linear order. The set of fragments represents the current state of the reconstruction procedure. Initially the set of fragments contains all the genes, but no red edges, so each gene is a fragment by itself.

Pathgroups

The objective function for the small phylogeny problem consists of the sum of a number of genomic distances, one distance for every branch in the phylogeny. Each of these distances corresponds to a breakpoint graph. A given genome determines blue edges in one breakpoint graph, while the red edges correspond to the ancestral genome being constructed. For each such ancestor, the red edges are identical in all the breakpoint graphs corresponding to distances to that ancestor. A pathgroup is a set of three paths, all beginning with the same vertex, one path from each partial breakpoint graph currently being constructed. Initially, there is one pathgroup for each non-T vertex. (We do not construct pathgroups for each T vertex separately, to be explained later, though paths ending in T vertices are found in other pathgroups.)

Pathgroups overlap because most paths are in two pathgroups, one associated with its initial vertex and one with its final vertex. With respect to a given path xy, we say the pathgroup determined by vertex x is the partner of the pathgroup determined by y. For the kind of binary (or bifurcating) trees we use, each pathgroup may have up to three distinct partners.

Priorities

A pathgroup may receive no priority, if creating any cycle within the pathgroup necessarily creates a circular chromosome. Note that in adding a red edge xy, this causes not only the disappearance of two partnered pathgroups, but it also changes paths in other pathgroups, which we call secondary pathgroups. Furthermore, each secondary pathgroup may itself have partner pathgroups whose paths, though not affected by the addition of xy, may have changed priorities. We call these tertiary pathgroups.

The pathgroups and the priorities are illustrated in Figure 2.

The makeCycles algorithm

By maintaining a list of pathgroups for each priority level, and a list of fragment endpoint pairs (initial and final), together with appropriate pointers, the algorithm makeCycles requires O(n) running time.

Algorithm makeCycles

input: pathgroups each consisting of three blue one-edge paths

output: ancestral genome

Not all the red fragments output by makeCycles are complete chromosomes of the ancestral genome; they may just be chromosome fragments. We know

Proposition [2]: Adding a red edge xy in a pathgroup creates at most four secondary pathgroups and at most eight tertiary pathgroups.

ensure the O(n) running time of the algorithm. If we had allowed red edges to connect T vertices or allowed pathgroups determined by T vertices, the number of potential secondary and tertiary pathgroups affected by the addition of a red edge would have increased considerably, depending on the number of chromosomes in the genomes, but this would not have affected the O(n) property.

Note that when a red edge is defined, the pathgroup is emptied, either by the creation of cycles, or by the integration of x as a non-endpoint of some path.

An example of the solution to the median problem is given as Additional File 1.

Small phylogeny

To apply PATHGROUPS to the small phylogeny problem, we set up an entire set of pathgroups for every internal (ancestral) node. Initially, in the pathgroups for those ancestral nodes connected to two given genomes, one of the paths will be missing and replaced by a single vertex of the breakpoint graph, as illustrated in Figure 2. Those ancestral nodes connected to only one given genome will have two missing paths in each pathgroup, both replaced by the vertex. Finally those ancestral nodes connected only to other internal nodes will have all paths missing in each pathgroup, all replaced by the vertex.

As the algorithm executes, a pathgroup with the highest priority found among any of the internal nodes is chosen to be processed next. Pathgroups connected to two given nodes will tend to be processed first, building up all three paths and combining the pathgroups one by one. Each time a red edge is added to a path, this becomes a blue edge in the corresponding pathgroup for ancestral genome(s) connected to it.

Eventually even the nodes furthest from any given genomes will accumulate enough edges in their pathgroups so that cycles can be formed and so that fragments of the associated genomes begin their reconstruction.

Refining the priorities

As a first improvement, we refine the priority levels as follows. We retain the number of cycles created as the primary classification, so that for example a two-cycle pathgroup always has higher priority than a one-cycle pathgroup. And we continue to refine this classification by taking account of the best new pathgroup that would be set up by processing the pathgroup under examination. In addition, however, to provide another level of refinement, we check all the pathgroups that would be affected by processing the pathgroup under examination and count the net change, positive or negative, in potential cycles among all these. This leads to 55 priority levels. As is seen in Figure 3, this refinement has a dramatic effect both in decreasing d/τ and in postponing the value of τ/n where d begins to rise much faster than τ.

Two-step look-ahead

As a further refinement, we considered the configuration that would be produced two algorithmic steps beyond the current step. Here, after identifying which one or several of the potential (after one step) pathgroups could produce the largest number of cycles by the addition of a red edge, we check what would happen after processing such a pathgroup, namely of the new (second step) pathgroups created, what is the largest number of cycles (1, 2 or 3) any one of them could produce by adding a red edge. This sub-classifies the 55 priorities three ways, creating a system with 165 priority levels.

As is clear in Figure 3, this additional step improves the accuracy even more than the first refinement did. The effect on reducing d/τ is especially strong, while there is little additional effect in delaying the point at which d begins to rise more quickly than τ.

The increase in run time caused by the two-step look-ahead is substantial, as we will learn from the next experiment. However, given that these results are based on rather highly arranged genomes containing 5000 genes, and a moderately large phylogeny (15-trees), the cost is hardly prohibitive.

Iterative local improvement based on the median

Searching for improvements in ancestral genome reconstruction by iteratively applying the median version at each ancestral node, accepting the changes only if they lower the objective function, is a time-honoured strategy in phylogenetics, including in rearrangement phylogeny [3, 11]. This would seem particularly appropriate in the present context since the median version of pathgroups is rapid and, as we have seen, relatively accurate.

Preliminary trials indicated that this procedure would be susceptible to premature capture by local minima. Experimenting with various regimes of simulated annealing, we settled on simply accepting every median reconstruction that was better than or equal to the existing reconstruction, and stopping after a predetermined number of steps.

Figure 4 depicts the results of this experiment (on 10-trees) up to 50 iterations. It can be seen that there is again a dramatic improvement in accuracy, both in reducing d/τ but also in delaying the point at which d begins to rise more quickly than τ. Simulations up to 100 iterations showed very small increases in accuracy over these results.

A time analysis (Figure 5) show that the iterative procedure exacts a higher cost than the other refinements. Run times for 50 iterations are ten times those shown, but again not enough to impede regular use of this method.

To what extent is the improvement seen here due to the initialization of the tree using PATHGROUPS and to what extent is it due to the iterative use of the median (also making use of PATHGROUPS of course)? After the necessary 50 or 100 iterations, the process no longer “remembers” its initialization, and neither approach seems more susceptible to falling into local optima. What about computing time? Figure 5 also shows what happens when the reconstructions are initialized with random genomes. With less rearranged genomes, there is a distinct time saving with the PATHGROUPS initialization, especially with larger phylogenies, but this disappears with more highly rearranged genomes.