Analysis of circular genome rearrangement by fusions, fissions and block-interchanges.

BACKGROUND: Analysis of genomes evolving via block-interchange events leads to a combinatorial problem of sorting by block-interchanges, which has been studied recently to evaluate the evolutionary relationship in distance between two biological species since block-interchange can be considered as a generalization of transposition. However, for genomes consisting of multiple chromosomes, their evolutionary history should also include events of chromosome fusions and fissions, where fusion merges two chromosomes into one and fission splits a chromosome into two.
RESULTS: In this paper, we study the problem of genome rearrangement between two genomes of circular and multiple chromosomes by considering fusion, fission and block-interchange events altogether. By use of permutation groups in algebra, we propose an O (n2) time algorithm to efficiently compute and obtain a minimum series of fusions, fissions and block-interchanges required to transform one circular multi-chromosomal genome into another, where n is the number of genes shared by the two studied genomes. In addition, we have implemented this algorithm as a web server, called FFBI, and have also applied it to analyzing by gene orders the whole genomes of three human Vibrio pathogens, each with multiple and circular chromosomes, to infer their evolutionary relationships. Consequently, our experimental results coincide well with our previous results obtained using the chromosome-by-chromosome comparisons by landmark orders between any two Vibrio chromosomal sequences as well as using the traditional comparative analysis of 16S rRNA sequences. ConclusionFFBI is a useful tool for the bioinformatics analysis of circular and multiple genome rearrangement by fusions, fissions and block-interchanges.

Background

For the past two decades, genome rearrangements have been studied and can be modelled to learn more about the evolution of mitochondrial, chloroplast, viral, bacterial and mammalian genomes [1]. To evaluate the evolutionary distance between two related genomes in gene order, various rearrangement events acting on genes within or among chromosomes have been proposed, such as reversals (also known as inversions) [1-10], transpositions [11,12], block-interchanges [13-15], translocations [16,17], and fusions and fissions [18,19]. Most genome rearrangement studies in computation involve the issue of solving the combinatorial problem to find an optimal series of rearrangements required to transform one genome into another.

Recently, the study on the genome rearrangement by block-interchanges has increasingly drawn great attention, since the block-interchange event is a generalization of transposition and, currently, its computational models measuring the genetic distance are more tractable than those modeled by transposition. Christie [13] first introduced the concept of block-interchange, affecting a chromosome by swapping two non-intersecting blocks containing any number of consecutive genes. Block-interchange can be considered as a generalization of transposition, since any exchanged blocks via transposition must be contiguous in a chromosome, whereas those via block-interchange need not be. As a matter of fact, the occurrence of an exchange of two non-contiguous blocks has been suggested in the previous studies related to the biological processes of bacterial replication [[20], and references therein]. Christie also proposed an (n2) time algorithm, where n is the number of genes, to solve the so-called block-interchange distance problem that is to find a minimum series of block-interchanges for transforming one linear chromosome into another. Later, we [14] designed a simpler algorithm for solving the block-interchange problem on linear or circular chromosomes with time-complexity of (δn), where δ is the the minimum number of block-interchanges required for the transformation and can be calculated in (n) time in advance. We also demonstrated that block-interchange events play a significant role in the genetic evolution of bacterial (Vibrio) species. Very recently, based on this algorithm, we have further implemented a tool, called ROBIN, for analyzing the rearrangements of gene orders via block-interchanges between two linear/circular chromosomal genomes [15]. Not only gene-order data but also sequence data are allowed to be input into the ROBIN system. If the input is the sequence data, ROBIN can automatically search for the common homologous/conserved regions shared by all input sequences.

It should be noted that the above block-interchange studies were dedicated to genomes containing only one chromosome (i.e., uni-chromosomal genomes) for evaluating their evolutionary relationships. However, for biological species with different numbers of chromosomes, the evolutionary history must also consider events of chromosomal fusions and fissions. A fusion occurs when two chromosomes merge into one and a fission takes place when a chromosome splits into two. The reason is that different chromosomes may as well exchange their genetic material with each other and, moreover, this exchange can only be achieved via inter-chromosomal operations such as fusions and fissions, instead of intra-chromosomal operations like block-interchanges. Hence, it is worthwhile to study genome rearrangements considering fusions, fissions and block-interchanges altogether. In this paper, we solve such a genome rearrangement problem by designing an efficient algorithm to compute and obtain a minimum series of all the events involving fusions, fissions and block-interchanges that are required to transform one circular multi-chromosomal genome into another, when both have the same set of genes without repeats. Although most eukaryotic genomes are linear, most prokaryotic (e.g., bacterial) genomes are circular and some of them consist of multiple circular chromosomes and large plasmids1. For example, some important bacterial pathogens like Brucella, Burkholderia, Leptospira and Vibrio species fall into this category. Notably, our approach is based on permutation group in algebra, instead of breakpoint graph, a commonly used approach in the study of genome rearrangement.

Recently, Yancopoulos et al. [21] used breakpoint graph to design an algorithm to solve a genome rearrangement problem in which the considered reversals, translocations (including fusions and fissions) and block-interchanges were given different weights. Unfortunately, their algorithm cannot be applied to solving our problem in which the events we considered are unweighted, because a series of weighted events with minimum weights in total may not be a minimum series of unweighted events, provided the events are given different weights.

Results and discussion

Based on Algorithm Sorting-by-ffbi developed in this study, we have implemented a web server, called FFBI [22], in which biologists or scientists in genomics can conduct comprehensive analyses of circular genome rearrangements by fusions, fissions and block-interchanges for their scientific interests and needs. Furthermore, we used this web server to conduct the rearrangement analyses on the whole genomes of three pathogenic Vibrio species, including V. vulnificus, V. parahaemolyticus and V. cholerae, to infer their evolutionary relationships.

Each of these three Vibrio pathogens consists of two circular chromosomes, and all their genomic sequences have recently been reported in GenBank with protein-coding genes annotated (see Table 1 for their sequence information). As annotated in GenBank (as of April 2006), the genomes of V. vulnificus, V. parahaemolyticus and V. cholerae contain 5098, 4992 and 4008 genes, respectively. From these protein-coding genes, we identified a total of 2393 (one-to-one) orthologous genes that are physically located in different positions on the chromosomes (see the Method section for construction of orthologous genes). Inevitably, there can be a high possibility that some genes with mis-annotated or uncertain protein functions are included in the genome annotation data. We therefore used only those authentic genes whose protein functions are not annotated as hypothetical or putative proteins, or are conserved and not poorly characterized (e.g., not those genes with only general function prediction or unknown function) in the NCBI COGs [23] database of orthologous genes. As a result, there are 1274 authentic orthologous genes in total remained for the further study of genome rearrangement. The relative orders of these orthologous genes along chromosomes, as well as the annotated COGs of their coding proteins, are detailed in the web site of our server.

Table 1

The sequence information of three pathogenic Vibrio species, each with two circular chromosomes.

Accession NO.	Species	Chromosome	Size (Mbps)
[GenBank:NC_005139]	V. vulnificus YJ016	1	3.4
[GenBank:NC_005140]	V. vulnificus YJ016	2	1.9
[GenBank:NC_004603]	V. parahaemolyticus RIMD 2210633	1	3.3
[GenBank:NC_004605]	V. parahaemolyticus RIMD 2210633	2	1.9
[GenBank:NC_002505]	V. cholerae El Tor N16961	1	3.0
[GenBank:NC_002506]	V. cholerae El Tor N16961	2	1.0

For each pair of these pathogenic Vibrio species, the variation in their gene orders has suggested that the genome rearrangement events have occurred and their genomes are closely related in evolution. To evaluate the contribution of fusions, fissions and block-interchanges to these observed rearrangements, we used the server developed in this study to compute the rearrangement distance between the gene orders of any two Vibrio genomes. Consequently, as shown in Table 2, the calculated rearrangement distance between V. vulnificus and V. parahaemolyticus is smaller than that between V. vulnificus and V. cholerae and that between V. parahaemolyticus and V. cholerae, suggesting that V. vulnificus is closer to V. parahaemolyticus than to V. cholerae in evolutionary relationship. Intriguingly, this result of genome-wide experiment well coincides with those we obtained in the previous chromosome-wide experiments [14,15]. Recall that our previous experiments were conducted in a chromosome-by-chromosome style of computing the block-interchange distance between the landmark orders of any two large/small Vibrio chromosomes, where the used landmarks are the maximal unique matches (MUMs) or the locally collinear blocks (LCBs) that are commonly shared by three large/small Vibrio chromosomes.

Table 2

The calculated rearrangement distances among V. vulnificus, V. parahaemolyticus and V. cholerae by fusions, fissions and block-interchanges.

Species Compared	V. vulnificus	V. parahaemolyticus	V. cholerae
V. vulnificus	0	174	364
V. parahaemolyticus	174	0	391
V. cholerae	364	391	0

In fact, the evolutionary relationships of the three pathogenic Vibrio species revealed in our experiment of analyzing their genome rearrangements also confirms that obtained by the biological community on the basis of the traditional comparative analysis of 16S rRNA gene sequences [24-26]. For confirmation, we here repeated this comparative analysis as follows. The 16S rRNA gene sequences of three Vibrios were first aligned using the Clustal W program [27], from which the distance matrix (as shown in Table 3) was then estimated by the algorithm of Kimura's two-parameter model in PHYLIP package [28].

Table 3

The calculated distances among V. vulnificus [GenBank:X76333], V. parahaemolyticus [GenBank:X56580] and V. cholerae [GenBank:X76337] by the traditional comparative analysis of their 16S rRNA gene sequences (accession numbers of 16S rRNAs are given in square brackets).

Species Compared	V. vulnificus	V. parahaemolyticus	V. cholerae
V. vulnificus	0.000000	0.034524	0.050261
V. parahaemolyticus	0.034524	0.000000	0.076739
V. cholerae	0.050261	0.076739	0.000000

Conclusion

In this paper, we studied the genome rearrangement problem between circular genomes with multiple chromosomes by simultaneously considering fusion, fission and block-interchange events. We have shown in the Method section that an optimal series of events required to transform one genome into another can be obtained in a canonical order such that all fusions come before all block-interchanges, which come before all fissions. Based on this property as well as the concept of permutation groups in algebra, we have successfully designed an (n2) time algorithm to obtain the minimum number of fusion, fission and block-interchange events for the transformation and also to generate an optimal scenario of the required rearrangement events. In addition, we have practically implemented this algorithm as a web server and applied it to analyzing by gene orders the whole genomes of three human Vibrio pathogens to infer their evolutionary relationships. As a consequence, our experimental results well coincide with the previous results obtained using the chromosome-by-chromosome comparisons by landmark orders between any two Vibrio chromosomal sequences as well as using the traditional comparative analysis of 16S rRNA sequences. The algorithm, however, should not be applied on linear multi-chromosomal genomes, because as mentioned in the Method section, it is not always possible to have an optimal scenario in a canonical order for linear genomes. Further studies in genome rearrangement can still be pursued to solve the problem for linear multi-chromosomal genomes.

Methods

Permutations versus genome rearrangements

In group theory, a permutation is defined to be a one-to-one mapping from a set E = {1, 2, ..., n} into itself, where n is some positive integer. For example, we may define a permutation α of the set {1, 2, 3, 4, 5, 6, 7} by specifying α(1) = 4, α(2) = 3, α(3) = 1, α(4) = 2, α(5) = 7, α(7) = 6 and α(6) = 5. The above mapping can be expressed using a cycle notation as illustrated in Figure 1 and simply denoted by α = (1, 4, 2, 3) (5, 7, 6). A cycle of length k, say (a1, a2, ..., a), is simply called k-cycle and can be rewritten as (a, a, ..., a, a1, ..., a), where 2 ≤ i 29]. Usually, a cycle of length one in α is not explicitly written and its element, say x, is said to be fixed by α since α(x) = x. Especially, the permutation whose elements are all fixed is called an identity permutation and is denoted by 1 (i.e., 1 = (1) (2) ... (n)).

Figure 1

The illustration of a permutation α = (1, 4, 2, 3) (5, 7, 6) meaning that α(1) = 4, α(2) = 3, α(3) = 1, α(4) = 2, α(5) = 7, α(7) = 6 and α(6) = 5.

Given two permutations α and β of E, the composition (or product) of α and β, denoted by αβ, is defined to be a permutation of E with αβ(x) = α(β(x)) for all x ∈ E. For instance, if we let E = {1, 2, 3, 4, 5, 6}, α = (2, 3) and β = (2, 1, 5, 3, 6, 4), then αβ = (2, 1, 5) (3, 6, 4). If α and β are disjoint cycles, then αβ = βα. The inverse of α is defined to be a permutation, denoted by α-1, such that αα-1 = α-1α = 1. If a permutation is expressed by the product of disjoint cycles, then its inverse can be obtained by just reversing the order of the elements in each cycle. For example, if α = (2, 1, 5) (3, 6, 4), then α-1 = (5, 1, 2) (4, 6, 3). Clearly, α-1 = α if α is a 2-cycle.

Meidanis and Dias [19,30] first noted that each cycle of a permutation may represent a circular chromosome of a genome with each element of the cycle corresponding to a gene, and the order of the cycle corresponding to the gene order of the chromosome. Figure 1, for example, shows a genome with two circular chromosomes, one represented by (1, 4, 2, 3) and the other by (5, 7, 6). Moreover, they observed that global evolutionary events, such as fusions and fissions (respectively, transpositions), correspond to the composition of a 2-cycle (respectively, 3-cycles) and the permutation representing a genome. For instance, let α be any permutation whose cycle decomposition is α1α2 ... α. If ρ = (x, y) is a 2-cycle and x and y are in the different cycles of α, say α= (a1 ≡ x, a2, ..., a) and α= (b1 ≡ y, b2, ..., b) where 1 ≤ p, q ≤ r, then in the composition ρα, αand αare joined into a cycle (a1, a2, ..., a, b1, b2, ..., b), i.e., ρ is a fusion event affecting on α (and also called a join operation of α here). If ρ = (x, y) is a 2-cycle and x and y are in the same cycle of α, say α= (a1 ≡ x, a2, ..., a≡ y, a, ..., a) where 1 ≤ p ≤ r, then in the composition ρα, this cycle αis broken into two disjoint cycles (a1, a2, ..., a) and (a, a, ..., a), i.e., ρ is a fission event affecting on α (and also called a split operation of α here). If ρ = (x, y, z) is a 3-cycle and x, y and z are in the same cycle of α, say α= (a1 ≡ x, a2, ..., a, b1 ≡ y, b2, ..., b, c1 ≡ z, c2, ..., c) where 1 ≤ p ≤ r, then in the composition ρα, the cycle αbecomes (a1, a2, ..., a, c1, c2, ..., c, b1, b2, ..., b), i.e., ρ is a transposition event affecting on α.

Recently, Lin et al. [14] further observed that a block-interchange event affecting on α corresponds to the composition of two 2-cycles, say ρ1 and ρ2, and α under the condition that ρ1 is a split operation of α and ρ2 is a join operation of ρ1α. More clearly, let α= (a1, a2,..., a) be a cycle of α, ρ1 = (a1, a) and ρ2 = (a, a), where 1

As discussed above, any series of fusions, fissions and block-interchanges required to transform one circular multi-chromosomal genome α into another I can be expressed by a product of 2-cycles, say ρρ...ρ1, such that ρρ...ρ1α = I (hence, ρρ...ρ1 = Iα-1). This property implies that Iα-1 contains all information that can be utilized to derive ρ1, ρ2,...,ρfor transforming α into I.

It is well known that every permutation can be written as a product of 2-cycles. For example, (1, 2, 3, 4) = (1, 4) (1, 3) (1, 2). However, there are many ways of expressing a permutation α as a product of 2-cycles [29]. Given a permutation α, let f(α) denote the number of the disjoint cycles in the cycle decomposition of α. Notice that f(α) counts also the non-expressed cycles of length one. For example, if α = (1, 5) (2, 4) is a permutation of E = {1, 2,..., 5}, then f(α) = 3, instead of f(α) = 2, since α = (1, 5) (2, 4) (3). Then the following lemma shows the lower bound of the number of 2-cycles in any product of 2-cycles of a permutation.

Lemma 1 Let α be an arbitrary permutation of E = {1, 2,..., n}. If α can be expressed as a product of m 2-cycles, say α = α1α2...αwith each α being 2-cycle, then m ≥ n - f(α).

Proof. We prove this lemma by induction on m. The lemma is true if m = 0, since α = 1 then, meaning that f(α) = n, and hence m = n - f(α) = 0. Suppose now that the lemma holds for any permutation that can expressed as a product of less than m 2-cycles, where m > 0. Let α' = α2α3, ...,α. Then by the induction hypothesis, we have m - 1 ≥ n - f(α'). Since α = α1α' and α1 is a 2-cycle, α1 operates on α' either as a fusion by joining two cycles of α' into one cycle (i.e., f(α) = f(α') - 1) or as a fission by splitting one cycle of α' into two cycles (i.e., f(α) = f(α') + 1). Whichever α1 operates on α', we have f(α) ≥ f(α') - 1. As a result, m = (m - 1) + 1 ≥ n - f(α') + 1 = n - (f(α') - 1) ≥ n - f(α).     □

Optimal scenario in canonical order

As mentioned previously, each circular multi-chromosomal genome with n genes can be expressed by a permutation of E = {1, 2,..., n}. Given two such genomes G1 and G2 over the same gene set E, the genome rearrangement distance between G1 and G2, denoted by d(G1,G2), is defined to be the minimum number of events needed to transform G1 into G2, where the events allowed to take place are fusions, fissions and block-interchanges. In this section, we shall show that there is an optimal series of events required to transform G1 into G2 such that all fusions come prior to all block-interchanges, which come before all fissions. Here, such an optimal scenario of genome rearrangements is referred as in canonical order.

Lemma 2 d(G1, G2) = d(G2, G1).

Proof. Let Φ = <σ1, σ2,...,σ> be an optimal series of events required to transform G1 into G2. Clearly, Φ' = <σ, σ,...,σ1> is an optimal series of events for transforming G2 into G1 by reversing the role of every event σ, where 1 ≤ i ≤ δ, such that σ is a fission (respectively, fusion) in Φ' if σis a fusion (respectively, fission) in Φ.

Lemma 3 There is an optimal series of events required to transform G1 into G2 such that every fission occurs after every fusion and block-interchange.

Proof. Let Φ = <σ1, σ2,...,σ> be an optimal series of events needed to transform G1 into G2. Of course, if every fission occurs after every fusion and block-interchange in Φ, then the proof is done. Now, suppose that not every fission occurs after every fusion or block-interchange in Φ. Then let i be the largest index in Φ such that σis a fission preceding σthat is either a fusion or a block-interchange. We can then obtain a new optimal series Φ' = <σ1,...,σ, , , σ,...,σ> to transform G1 into G2 such that is a fusion or a block-interchange and is a fission, as discussed below. Suppose that σ splits a chromosome α into α1 and α2. If σis a fusion, then we assume that it joins two chromosomes β1 and β2into β; otherwise, if σis a block-interchange, then assume that it affects β1 such that β1 becomes β through a block-interchange. Clearly, if neither β1 nor β2 is created by σ, then the desired series Φ' is obtained by swapping σ and σin Φ (i.e., = σand = σ). If both β1 and β2 are created by σ, then the net rearrangement of σ (a split operation) followed by σ(a joint operation) either has no effect on α or becomes a block-interchange affecting α. By removing σand σfrom Φ or replacing them with an extra block-interchange, we thus obtain a new optimal series of the events transforming G1 into G2 with strictly less than δ events, a contradiction. Hence, we assume that only one of β1 and β2 is created by σ and without loss of generality, let β1 = α1. Now, we consider the following two cases.

Case 1: σis a fusion. For simplicity of discussion, we let α = (1, 2,..., x - 1, x,..., y - 1), β2 = (y, y + 1,..., z), σ= (1, x) and σ= (1, y), where 1

Case 2: σis a block-interchange. Clearly, the net rearrangement caused by σ and σis to transform α into β and α2, which is equivalent to the rearrangement by first applying σto α and then further splitting it into β and α2 via σ. Then Φ' is obtained by swapping σ and σin Φ.

In other words, we can always obtain Φ' from Φ according to the method described above. Repeating this process on the resulting Φ', we can finally obtain an optimal series of events that are required to transform G1 into G2 such that all fissions come after all fusions and block-interchanges.

Lemma 4 There is an optimal series of events required to transform G1 into G2 in a canonical order such that all fusions come before all block-interchanges, which come before all fissions.

Proof. Let Φ = <σ1, σ2,..., σ> be an optimal series of events required to transform G1 into G2. If there are no fusions or block-interchanges, then the proof is completed. If not, according to Lemma 3, we may assume that all fusions and block-interchanges occur earlier than all fissions. Let i be the index of the last non-fission in Φ and also let G' be the resulting genome after all σ1, σ2,...σhave affected G1. Since Φ is optimal, it is straightforward to see that Φ= <σ1, σ2,...,σ> is an optimal series of fusions and block-interchanges needed to transform G1 into G'. As discussed in the proof of Lemma 2, = <σ, σ,...,σ1> is an optimal series of fissions and block-interchanges for transforming G' into G1. Moreover, by Lemma 3, we can obtain = <, ,...,> from for transforming G' into G1 such that all block- interchanges in occur prior to all fissions. Consequently, <, ,...,> is an optimal series of fusions and block-interchanges needed to transform G1 into G', and all its fusions occur before all its block-interchanges. Therefore, there is an optimal series of events needed to transform G1 into G2 such that all fusions come earlier than all block-interchanges, which come before all fissions.     □

It is worth mentioning that an optimal scenario in a canonical order does not necessarily exist for linear multi-chromosomal genomes. For example, suppose that G1 and G2 are two given linear multi-chromosomal genomes, where G1 = (1, 4, 5) (2, 3) and G2 = (1, 2, 3) (4, 5). Then the optimal scenario between them is a fission, splitting (1, 4, 5) into (1) (4, 5), followed by a fusion, joining (1) and (2, 3) to (1, 2, 3). However, this optimal scenario can not be transformed into another in the canonical order according to the steps as described in Lemmas 3 and 4. Actually, there is no an optimal scenario between such two linear genomes using any two rearrangement events that begin with a fusion.

Algorithm

Let α and I be two given circular multi-chromosomal genomes over the same gene set E = {1, 2,...,n}. Here, we assume that the genes in I are sorted in the order of increasing and consecutive numbers, and that gene i + 1 is on the right side of gene i within the same chromosome, where 1 ≤ i ≤ n - 1. For example, I = (1, 2) (3, 4, 5) (6, 7, 8, 9) if I has three circular chromosomes with two, three and four genes, respectively. In this case, the computation of d(α,I) and its corresponding optimal scenario can be considered as a problem of sorting α using the minimum set of operations, including fusions, fissions and block-interchanges.

Suppose that ρρ... ρ1 is a product of 2-cycles that corresponds to an optimal series Φ of fusions, fissions and block-interchanges for transforming α into I. Then Φ = ρρ...ρ1 = Iα-1 and by Lemma 1, λ ≥ n - f(Iα-1). If Φ' = is another product of λ 2-cycles with Φ' = Iα-1, then the number of 2-cycles in Φ' that function as fusions or fissions must be greater than or equal to that in Φ; otherwise, the total number of fusions, fissions and block-interchanges in Φ' for transforming α into I must be less than that in Φ, a contradiction. The reason is that a fusion or fission requires only one 2-cycle for rearrangement, whereas a block-interchange requires two 2-cycles. In other words, the number of 2-cycles serving as the fusions and fissions is minimum in any optimal series of events.

Based on the above observation as well as Lemma 4, below we design an efficient algorithm for computing d(α,I) and its optimal scenario of rearrangement events in a canonical order. Let χ(α) and χ(I) denote the numbers of chromosomes in α and I, respectively, and let α = α1α2...αand I = I1I2...I. Then, an undirected graph (α,I) = (, , ) is constructed from α and I as follows.

• = {α1, α2,...,α}.

• = {I1, I2,...,I}.

• = {(α, I) | 1 ≤ i ≤ χ (α), 1 ≤ j ≤ χ (I), and αand Ihave at least a common gene}.

For instance, suppose that α = α1α2...α5 = (1, 2, 10) (11, 8, 9, 3, 6) (7, 4, 5,12) (13, 15) (14, 16) and I = I1I2 ... I6 = (1, 2, 3) (4, 5) (6, 7, 8) (9, 10, 11, 12) (13, 14) (15, 16). Then the induced graph (α,I) is shown in Figure 2, which is a bipartite graph since and are independent sets in (α, I) (i.e., no edge between any two vertices in or ). A connected component of (α,I) is defined to be a maximal subgraph of (α,I) such that there exists a path between any pair of vertices in this subgraph. For example, the induced (α,I) as shown in Figure 2 has two connected components. Notice that if in a chromosome Iof I there are two genes that appear in two different chromosomes α and αof α, then (α, I) ∈ and (α, I) ∈ , and hence both α and α belong to the same connected component in (α,I).

Figure 2

The induced bipartite graph (α, I) with two connected components.

Let {1, 2,...,} denote the collection of all connected components in (α,I). For each 1 ≤ i ≤ ω, let βand J denote the chromosomes in α and I, respectively, whose corresponding vertices belong to in (α,I). Let gene(β) and gene(J) be the collections of the genes in all chromosomes of βand J, respectively. Then gene(β) = gene(J) and gene(β) ∩ gene(β) = ∅ for any 1 ≤ j ≠ i ≤ ω. Let n be the number of genes in gene(β). Clearly, n = n1 + n2 + ... + n. In addition, it can be verified that Iα-1 = (J1) (J2) ... (J) and f(Iα-1) = f(J1) + f(J2) + ... + f(J).

According to the properties above, we then find a product Φ of 2-cycles so that Φα = I as follows. We first find a product Φof 2-cycles that corresponds to an optimal series of the rearrangement events required to transform βinto J(i.e., Φβ= J and hence Φ= J) and then let Φ = ΦΦ... Φ1. Clearly, Φ = Iα-1 and hence Φ corresponds to a feasible series of events for transforming α into I. Actually, we shall show later that the number of 2-cycles in each Φis n- f(J), and in Φ the number of 2-cycles functioning as the fusions and fissions is minimum. This causes that the number of 2-cycles in Φ equals to , in which the number of 2-cycles serving as the fusions and fissions is minimum. As a result, Φ is an optimal series of events for transforming α into I. The above description indicates that the original problem can be conquered by independently solving the same problem on the smaller instance whose induced bipartite graph is a connected component of (α, I).

To simplify our discussion, throughout the rest of this section we assume that the induced (α, I) of a given instance α and I has exactly one connected component. Let Φ = <σ1, σ2, ...., σ> be an optimal series of events for transforming α into I in which all fusions precede all block-interchanges that further precede all fissions. Let n, nand ndenote the numbers of fusions, block-interchanges and fissions, respectively, in Φ. Then δ = n+ n+ n. In the following, we shall show that Φ can be expressed by a product of n - f(Iα-1) 2-cycles in which the number of 2-cycles functioning as the fusions and fissions is minimum.

It should be noticed that the chromosomes considered here are disjoint (i.e., without gene duplication). Hence, for any two chromosomes α and αin α with (αI) ∈ and (α, I) ∈ , there must exist a fusion in Φ that joins α and αto one chromosome; otherwise, Ican not be formed from α by a fission later. Since all needed fusions come together in the beginning of Φ, n= χ(α) - 1, which is the lower bound of the number of fusions required in any optimal series of events for transforming α into I. After these nfusions, the resulting α becomes only one chromosome. Since the next nblock-interchanges are intra-chromosomal mutations, we have n= χ(I) - 1. Actually, χ(I) - 1 is the minimum number of the required fissions in any optimal series of events for transforming α into I, since it is the minimum number of the fusions used in the corresponding optimal series of events to transform I into α.

Given any cycle ρ, we use x ∈ ρ to denote that x is a number in ρ. For any two x ∈ ρ and y ∈ ρ, they are said to be adjacent in ρ if ρ(x) = y or ρ(y) = x. Next, we show a way to derive n2-cycles from Iα-1 such that these 2-cycles function as the fusions that join all chromosomes of α to a single one, if α has multiple chromosomes, where n= χ(α) - 1. For simplicity, later in the text we use "cycle in Iα-1" to represent " cycle in the cycle decomposition of Iα-1" in meaning, unless a possible confusion may arise.

Lemma 5 Let αand αbe any two disjoint cycles in α. Then there must exist a cycle in Iα-1 that contains two numbers x and y such that x ∈ αand y ∈ α.

Proof. Since we assume that the induced (α, I) contains exactly and only one connected component, and αand αcontain some numbers u and v, respectively, such that both u and v are in a cycle Iof I. Notice that u ∉ αand v ∉ α. Suppose that there is no cycle in Iα-1 that contains two numbers x and y such that they are in these two different cycles of α, say x ∈ α and y ∈ α. Then all numbers in any cycle of Iα-1 are contained in some cycle of α. Without loss of generality, let I= (u ≡ a1, a2, ..., a≡ v, ..., a) and let p

The following lemma can be easily verified.

Lemma 6 (a1, a2,..., a,..., a) = (a1, a2,..., a) (a,..., a) (a, a), where 1 ≤ i

According to Lemma 5, for any two cycles α and αof α, we can find two numbers x and y in a cycle of Iα-1, say β, such that x ∈ α and y ∈ α. Let β = (a1, a2,..., a), where q ≥ 2. Then we consider the following two cases. Case 1: x and y are adjacent in β. For simplicity, let x = aand y = a. Then by Lemma 6, β = (a1, a2,..., a) (a) (x, y). Case 2: x and y are not adjacent in β. Let x = aand y = a, where 1 ≤ p Lemma 6. In other words, we can derive a 2-cycle (x, y) from β such that it can join α and αto one cycle. After α, and αare joined together via (x, y), the number of the cycles (including 1-cycles) in the resulting Iα-1 increases by one. Repeatedly based on the procedure above, we can derive consecutive n2-cycles from Iα-1, say φ1, φ2,...,, that can join χ (α) cycles in α to a single one, where n= χ(α) - 1. In other words, φ1, φ2,..., function as χ(α) - 1 fusions that transform genome α with χ(α) chromosomes into a genome, denoted by α', with a single chromosome. Clearly, we have α' = ...φ1α, Iα'-1 = Iα-1 φ1φ2 ... , and f(Iα'-1) = f(Iα-1) + n. Hence, we can immediately claim the following.

Claim 1 α' = ...φ1α, Iα'-1 = Iα-1 φ1φ2 ... , and f(Iα'-1) = f(Iα-1) + n, where n= χ(α) - 1.

Without loss of generality, we now suppose that χ(I) > 1. Similarly as the discussion above, we can derive consecutive n2-cycles from α'I-1, say ψ1, ψ2, ..., , such that they serve as the fusions to transform I with χ(I) chromosomes into a genome, denoted by I', with only one chromosome, where n= χ(I) - 1 and α'I-1 is the inverse of Iα'-1 (i.e., α'I-1 = (Iα'-1)-1). Then we have I' = ... ψ1I (hence ψ1ψ2 ... I' = I), α'I'-1 = α'I-1ψ1ψ2 ... and f(α'I'-1) = f(α'I-1) + n. Conversely, we can use , …, ψ1 as fissions to split I' with one chromosome into I with nchromosomes. Since α'I-1 is the inverse of Iα'-1, it can be easily obtained from Iα'-1 by just reversing the order of the numbers in each cycle of Iα'-1 and hence f(α'I-1) = f(Iα'-1), which leads to f(α'I'-1) = f(Iα'-1) + n. As a result, we have f(I'α'-1) = f(α'I-1) = f(Iα'-1) + nsince I'α'-1 = (α'I-1)-1. Therefore, the following claim can be obtained.

Claim 2 ψ1ψ2 ... I' = I, α'I'-1 = α'I-1ψ1ψ2 ... and f(I'α'-1) = f(Iα'-1) + n, where n= χ(I) - 1.

Notice that both α' and I' now are the genomes with only one chromosome. Then based on the algorithm proposed by Lin et al. [14], we can find block-interchanges from I'α'-1 to transform α' into I'. Certainly, these nblock-interchanges can be further expressed by a product of 2n2-cycles, say , such that every two consecutive 2-cycles act as a block-interchange in the process of transforming α' into I', where I'α'-1 = . Hence, we have the following claim immediately.

Claim 3 I' = α'.

Now we let Φ = ψ1ψ2.... Then the result of Φα = I (hence Φ = Iα-1) can be easily verified by Claims 1, 2 and 3 as follows.

In other words, Φ is a product of (n+ (n - f(I'α'-1)) + n) 2-cycles that can transform α into I. More clearly, Φ first uses φ1, φ2,..., (acting as nfusions) to transform α into α', then uses (acting as nblock-interchanges) to transform α' into I', and finally uses , ,..., ψ1 (acting as nfissions) to transform I' into I. By Claims 1 and 2, we can show that (n+ (n - f(I'α'-1)) + n= n - f(Iα-1) as follows.

n+ (n - f(I'α'-1)) + n

= n+ (n - f(Iα'-1) + n) + n    (by Claim 2)

= n+ (n - f(Iα-1) + n- n) + n    (by Claim 1)

= n - f(Iα-1)

As mentioned before, χ(α) - 1 and χ(I) - 1 are the lower bounds of the numbers of fusions and fissions, respectively, required in any optimal series of rearrangement events for transforming α into I. Hence, the number of 2-cycles in Φ that function as the fusions and fissions is minimum. Along with that Φ = Iα-1 can be expressed as a product of n - f(Iα-1) 2-cycles, we thus conclude that Φ is an optimal series of the events that transform α into I with first nfusions, then nblock-interchanges and finally nfissions, where n= χ(α) -1, , and n= χ(I) - 1.

Lemma 7 There is an optimal series of the events needed to transform α into I in a canonical order such that all nn= χ(α) - 1, and n= χ(I) - 1.

Let us take α = (1, 2,10) (11, 8, 9, 3, 6) (7, 4, 5, 12) and I = (1, 2, 3) (4, 5) (6, 7, 8) (9, 10, 11, 12) for an example. It should be straightforward to see that (α,I) is a connected bipartite graph with χ(α) = 3 and χ(I) = 4, and Iα-1 = (1, 11, 7, 9, 6) (3, 10) (4, 8, 12) and hence f(Iα) = 5, since two 1-cycles (i.e., (2) and (5)) are not explicitly shown. First, we are to find two 2-cycles φ1 and φ2 (since n= χ(α) -1 = 2) from Iα-1 to transform genome α with three chromosomes into genome α' with exactly one chromosome. To this purpose, we let φ1 = (3,10) and φ2 = (4,8), since Iα-1 = (1, 11, 7, 9, 6) (4, 12) (4, 8) (3, 10). Then by Claim 1, α' = φ2φ1α = (4, 5, 12, 7, 8, 9, 10, 1, 2, 3, 6, 11) and Iα'-1 = Iα-1φ1φ2 = (1, 11, 7, 9, 6) (4, 12). Next, we need to find three 2-cycles ψ1, ψ2 and ψ3 (since n= χ(I) - 1 = 3) from α'I-1, which is equal to (Iα'-1)-1 = (6, 9, 7, 11, 1) (12, 4) = (1, 7, 11) (1, 9) (1, 6) (12, 4), to transform I into I' with only one chromosome. By letting ψ1 = (12, 4), ψ2 = (1, 6) and ψ3 = (1, 9), we have I' = ψ3ψ2ψ1 = (1, 2, 3, 6, 7, 8, 9, 10, 11, 4, 5, 12) and α'I'-1 = α'I-1ψ1ψ2ψ3 = (1, 7, 11) according to Claim 2. Finally, we will find two 2-cycles and (since n - f(Iα-1) - n- n= 12 - 5 - 2 - 3 = 2) from I'α'-1 that act as a block-interchange to transform α' into I', where I'α'-1 = (α'I'-1)-1 = (11, 7, 1) = (11, 1) (11, 7). By letting = (11, 7) and = (11, 1), we have α' = (11, 1) (11, 7) (4, 5, 12, 7, 8, 9, 10, 1, 2, 3, 6, 11) = (11, 4, 5, 12, 1, 2, 3, 6, 7, 8, 9, 10), which indeed equals I'. Consequently, we find an optimal series of events Φ= ψ1ψ2ψ3φ2φ1 that transform α into I (i.e., Φα = I).

Based on the idea above, we have designed Algorithm Sorting-by-ffbi (meaning sorting by fusions, fissions and block-interchanges) to compute the genome rearrangement distance d(α, I) between two given circular multi-chromosomal genomes α and I, and also to generate an optimal scenario of the required rearrangement events in a canonical order. In Algorithm Sorting-by-ffbi, the purpose of Step 2.3.3 (respectively, Step 2.4.4) is to find two numbers x and y that are both in some cycle of γ = J (respectively, γ = ), but in different cycles in β. By Lemma 5, such x and y exist. In fact, they can be found using the following simple approach. For simplicity, let γ= () be a cycle in γ that contains two numbers x and y such that they are in different cycles in β. Then we only need to check whether and , where 2 ≤ j ≤ l, are in different cycles in βor not. If so, we let x = and y = . The reason is as follows. Suppose that both and for all 2 ≤ j ≤ lare in the same cycle in β. Then all of numbers in γare in the same cycle in β, which contradicts the above assumption that γcontains x and y that are in different cycles in β.

Algorithm sorting-by-ffbi

Input: Two circular multi-chromosomal genomes α and I.

Output: d(α,I) and a minimum series Φ of events required to transform α into I.

1:   Find all connected components , ,...,in graph (α,I);

/* Denote by n the number of genes in . */

2:   for each , 1 ≤ i ≤ ω,do

/* Denote by β(resp. J) the collection of chromosomes in α (resp. I) whose corresponding vertices are in . */

2.1:   Compute J and let γ = J;

2.2:   n= χ(β) - 1, n= χ(J) -1, and δ= n+ n+ n;

2.3:   if χ(β) > 1 then   /* To compute φ1, φ2,..., */

2.3.1:      for each cycle of βdo

Create a set to contain all the numbers in this cycle;

endfor

2.3.2:   /* Let γ1γ2 ... γbe the cycle decomposition of the current γ and

let γ= (), where 1 ≤ q ≤ p and l≥ 2 */

k = 1 and h = 2;

2.3.3:   for j = 1 to ndo

S = find-set ();

while S = find-set()) do

if h

endwhile

x = and y = ;

φ= (x, y) and union(x, y);

endfor

2.3.4:    = ...φ1βand γ = γφ1φ2...;   /* Currently, γ is Ji */

endif

2.4:   if χ(J) > 1 then/* To compute ψ1, ψ2,..., */

2.4.1:      γ = γ-1/* New γ becomes */

2.4.2:   for each cycle of J do

Create a set to contain all the numbers in this cycle;

endfor

2.4.3:   /* Let γ1γ2...γbe the cycle decomposition of the current γ and

let γ= (), where 1 ≤ q ≤ p and l≥ 2 */

k = 1 and h = 2;

2.4.4:   for j = 1 to ndo

S = find-set();

while (S = find-set()) do

if h

endwhile

x = and y = ;

ψ= (x, y) and union(x, y);

endfor

2.4.5:    = ... ψ1Jand γ = γψ1ψ2...;   /* Currently, γ is */

endif

2.5:   /* To compute */

2.5.1:   γ = γ-1;/* New γ becomes */

2.5.2:   ;

2.5.3:   for j = 1 to ndo

Arbitrarily choose two adjacent elements x and y in γ;

/* Let = (a1,a2...,) */

Circularly shift (a1, a2,...,) such that a1 = x and assume y = a;

= (x, y);

for h = 1 to ndo

index (a) = h;

end for

Find two adjacent elements u and v in γ (x, y) such that

index(u) ≤ k - 1 and index(v) ≥ k;

= (u, v);

= and γ = γ;

endfor

3:   Let for each 1 ≤ i ≤ ω;.

4:   Output d(α, I) = and Φ = Φ1Φ2...Φ;

Theorem 1 Given two circular multi-chromosomal genomes α and I over the same gene set E = {1,2,..., n}, the problem of computing the genome rearrangement distance between α and I using fusions, fissions and block-interchanges can be solved and an optimal series of such events in a canonical order can be obtained in (n2) time.

Proof. As discussed above, Algorithm Sorting-by-ffbi transforms α into I using the minimum number of fusions, fissions and block-interchanges. Next, we follow to analyze its time-complexity. Notice that given an undirected graph with p vertices and q edges, all the connected components in this graph can be found in (p + q) time using depth-first search or breadth-first search [31]. As a result, Step 1 can be done in (n2) time for computing the connected components in the induced bipartite graph (α, I), since in the worst case, the number of edges in (α,I) is χ(α) × χ(I) and χ(α) = (n) and χ(I) = (n). In Step 2, there are ω outer iterations, each computing the minimum series of events needed to transform βinto J, where 1 ≤ i ≤ ω. Clearly, Steps 2.1 costs (n) time for computing J and Step 2.2 takes only a constant time. The time cost of Step 2.3 is mostly contributed from that of Step 2.3.3. There are ninner iterations in Step 2.3, each with the purpose of finding two numbers x and y that are both in the same cycle in γ and, however, in the different cycles in β. In the worst case, Step 2.3 needs n find-set operations and nunion operations to finish its overall process. Note that Step 2.3.1 can be implemented by initially creating a set for each number in gene(β) and then performing n- χ(β) union operations to generate χ(β) sets with each corresponding to a cycle in β, where χ(β) = n+ 1. Hence, the total number of union operations is n- 1 in Step 2.3. In fact, these find-set and union operations can be implemented in (n) time using the so-called "static disjoint set union and find" algorithm proposed by Gabow and Tarjan [32]. In other words, Step 2.3 can cost only (n) time. By the same principle, it can be verified that the time cost of Step 2.4 is (n). As for Step 2.5, adopted from our previous work [14], it takes (nn) time, where = n). As a result, the time cost of Step 2 is (ζn), where ζ is the maximum namong all iterations in Step 2 and ζ

Construction of orthologous genes

To analyze the rearrangement of three Vibrio genomes, we identified and constructed a table of orthologous genes that are putatively not involved in horizontal gene transfer (HGT) events by adopting the so-called symmetrical best hits (SymBets for short). In principle, two genes match and give SymBets if they are more similar to each other than they are to any other genes from the compared genomes [33,34]. Detection of such SymBet genes is, yet arguably, the simplest and most suitable method for identification of probable orthologs for closely related genomes [33,34]. Particularly, this prediction of orthologs holds to be applicable even when sequence similarity between the compared proteins is relatively low [33,34]. For the purpose of excluding paralogous genes derived from lineage-specific gene duplications, we here considered only one-to-one orthologous genes, which actually have been demonstrated as a major pattern in prokaryotic genome evolution [34]. Therefore, we used the following steps to identify an HGT-free table of one-to-one orthologous genes from the three complete Vibrio genomes. First, the GenePlot [35] program offered by NCBI was utilized to find and construct a table of SymBet genes between each pair of Vibrio genomes. Next, after removing all one-to-many or many-to-many SymBets, the three resulting tables of SymBet genes were joined to give a new one of one-to-one orthologous genes for all the three Vibrio genomes by using a rule as follows. If genes a (from genome A) and b (from genome B),b and c (from genome C), and c and a are all one-to-one SymBet pairs, then a, b and c are considered as one-to-one orthologous genes for the genomes A,B and C. In other words, the SymBet relationships among a, b and c result in a triangle. Finally, those genes that were involved in the putative HGT events detected and available in the Horizontal Gene Transfer Database [36] were then deleted from the table of one-to-one orthologous genes.

Authors' contributions

CLL and HTC contributed equally to this work. CLL conceived of the study, participated in the design and analysis of algorithm and drafted the manuscript. YLH participated in the algorithm design and software development, and carried out the bioinformatics experiments. TCW participated in the software development. HTC participated in the design and coordination of this study as well as in drafting the manuscript.

Note

1The Institute for Genome Research (TIGR) offers a website of Comprehensive Microbial Resource (CMR) that provides information on all of the publicly available and complete bacterial genomes.