On the family-free DCJ distance and similarity

Structural variation in genomes can be revealed by many (dis)similarity measures. Rearrangement operations, such as the so called double-cut-and-join (DCJ), are large-scale mutations that can create complex changes and produce such variations in genomes. A basic task in comparative genomics is to find the rearrangement distance between two given genomes, i.e., the minimum number of rearragement operations that transform one given genome into another one. In a family-based setting, genes are grouped into gene families and efficient algorithms have already been presented to compute the DCJ distance between two given genomes. In this work we propose the problem of computing the DCJ distance of two given genomes without prior gene family assignment, directly using the pairwise similarities between genes. We prove that this new family-free DCJ distance problem is APX-hard and provide an integer linear program to its solution. We also study a family-free DCJ similarity and prove that its computation is NP-hard.


Background
Genomes are subject to mutations or rearrangements in the course of evolution. Typical large-scale rearrangements change the number of chromosomes and/or the positions and orientations of genes. Examples of such rearrangements are inversions, translocations, fusions and fissions. A classical problem in comparative genomics is to compute the rearrangement distance, that is, the minimum number of rearrangements required to transform a given genome into another given genome [1].
In order to study this problem, one usually adopts a high-level view of genomes, in which only "relevant" fragments of the DNA (e.g., genes) are taken into consideration. Furthermore, a pre-processing of the data is required, so that we can compare the content of the genomes.
One popular method, adopted for more than 20 years, is to group the genes in both genomes into gene families, so that two genes in the same family are said to be equivalent. This setting is said to be family-based. Without gene duplications, that is, with the additional restriction that each family occurs exactly once in each genome, many polynomial models have been proposed to compute the genomic distance [2][3][4][5]. However, when gene *Correspondence: jens.stoye@uni-bielefeld.de 2 Technische Fakultät and CeBiTec, Universität Bielefeld, Universitätsstr. 25, 33615 Bielefeld, Germany Full list of author information is available at the end of the article duplications are allowed, the problem is more intrincate and all approaches proposed so far are NP-hard, see for instance [6][7][8][9][10].
It is not always possible to classify each gene unambiguously into a single gene family. Due to this fact, an alternative to the family-based setting was proposed recently and consists in studying the rearrangement distance without prior family assignment. Instead of families, the pairwise similarity between genes is directly used [11,12]. This approach is said to be family-free. Although the family-free setting seems to be at least as difficult as the family-based setting with duplications, its complexity is still unknown for various distance models.
In this work we are interested in the problem of computing the distance of two given genomes in a family-free setting, using the double cut and join (DCJ) model [5]. The DCJ operation, that consists of cutting a genome in two distinct positions and joining the four resultant open ends in a different way, represents most of large-scale rearrangements that modify genomes. After preliminaries and a formal definition of the family-free DCJ distance, we present a hardness result, before giving a linear programming solution and showing its feasibility for practical problem instances. Finally, we also study the problem of computing the similarity -a counterpart of the distance function -of two given genomes in a family-free setting using the DCJ model and show its NP-hardness. This paper is an extended version of [13], that was presented at the 14th Workshop on Algorithms in Bioinformatics, WABI 2014.

Preliminaries
Each gene g in a genome is an oriented DNA fragment that can be represented by the symbol g itself, if it has direct orientation, or by the symbol −g, if it has reverse orientation. Furthermore, each one of the two extremities of a linear chromosome is called a telomere, represented by the symbol •. Each chromosome in a genome can be represented by a string that can be circular, if the chromosome is circular, or linear and flanked by the symbols • if the chromosome is linear. For the sake of clarity, each chromosome is also flanked by parentheses. As an example, Since a gene g has an orientation, we can distinguish its two ends, also called its extremities, and denote them by g t (tail) and g h (head). An adjacency in a genome is either the extremity of a gene that is adjacent to one of its telomeres, or a pair of consecutive gene extremities in one of its chromosomes. If we consider again the genome A above, the adjacencies in its first chromosome are 3 t , 3 h 1 h , 1 t 4 t , 4 h 2 t and 2 h . Throughout this paper, let A and B be two distinct genomes and let A be the set of genes in genome A and B be the set of genes in genome B.

Adjacency graph and family-based DCJ distance
In the family-based setting the two genomes A and B have the same content, that is, A = B. When there are no duplications, that is, when each family is represented by exactly one gene in each genome, the DCJ distance can be easily computed with the help of the adjacency graph AG(A, B), a bipartite multigraph such that each partition corresponds to the set of adjacencies of one of the two input genomes and an edge connects the same extremities of genes in both genomes. In other words, there is a one-to-one correspondence between the set of edges in AG(A, B) and the set of gene extremities. Vertices have degree one or two and thus an adjacency graph is a collection of paths and cycles. An example of an adjacency graph is given in Figure 1.
The family-based DCJ distance d DCJ between two genomes A and B without duplications can be computed in linear time and is closely related to the number of components in the adjacency graph AG(A, B) [2]: where n = |A| = |B| is the number of genes in both genomes, c is the number of cycles and i is the number of odd paths in AG(A, B). Observe that, in Figure 1, the number of genes is n = 4 and AG(A, B) has one cycle and two odd paths. Consequently the DCJ distance is d DCJ The formula for d DCJ (A, B) can also be derived using the following approach. Given a component C in AG(A, B), let |C| denote the length, or number of edges, of C. From [14,15] we know that each component in AG(A, B) contributes independently to the DCJ distance, depending uniquely on its length. Formally, the contribution d(C) of a component C in the total distance is given by: ifC is an odd path , |C| 2 , ifC is an even path .
The sum of the lengths of all components in the adjacency graph is equal to 2n. Let C, I, and P represent the sets of components in AG(A, B) that are cycles, odd paths and even paths, respectively. Then, the DCJ distance can be calculated as the sum of the contributions of each component:

Gene similarity graph for the family-free model
In the family-free setting, each gene in each genome is represented by a distinct symbol, thus A ∩ B = ∅ and the cardinalities |A| and |B| may be distinct. Let a be a gene in A and b be a gene in B, then their normalized similarity is given by the value σ (a, b) that ranges in the interval [0, 1]. We can represent the similarities between the genes of genome A and the genes of genome B with respect to σ in the so called gene similarity graph [12], denoted by GS σ (A, B). This is a weighted bipartite graph whose partitions A and B are the sets of genes in genomes A and B, respectively. Furthermore, for each pair of genes (a, b), such that a ∈ A and b ∈ B, if σ (a, b) > 0 there is an edge e connecting a and b in GS σ (A, B) whose weight is σ (e) := σ (a, b). An example of a gene similarity graph is given in Figure 2.

Reduced genomes and their weighted adjacency graph
Observe that, for each edge e ∈ M, we have two edges of weight Examples of weighted adjacency graphs are shown in Figure 3.

The family-free DCJ distance
Based on the weighted adjacency graph, in [12] a familyfree DCJ similarity measure has been proposed. We will come back to this measure later in this paper. Before that, to be more consistent with the comparative genomics literature, where distance measures are more common than similarities, here we also propose a family-free DCJ distance. This family-free distance is based on the weighted DCJ distance of reduced genomes. An important design criterion for this definition is that it must be the same as the (unweighted) family-based DCJ distance when all weights are equal to 1.
The first step in our definition is to consider the components of the graph AG σ (A M , B M ) separately, similarly to the approach described previously for the family-based model. Here the contribution of each component C is denoted by d σ (C) and must include not only the length |C| of the component, but also information about the weights of the edges in C. Basically, we need a function f (C) to use instead of |C| in the contribution function d σ (C), such that: (i) when all edges in C have weight 1, f (C) = |C|, that is, the contribution of C is the same as in the family-based version; (ii) when the weights decrease, f should increase, because smaller weights mean less similarity, or increased distance between the genomes.
The simplest linear function f that satisfies both condi- is the sum of the weights of all the edges in C. Then, the weighted contribution d σ (C) of the different types of components is: ifC is an even path .
Let C, I, and P represent the sets of components in AG σ (A M , B M ) that are cycles, odd paths and even paths, respectively. Summing the contributions of all the components, the resulting distance for a certain matching M is computed as follows: since the number of genes in G(M) is equal to the size of M.
Observe that not only the components of the graph, but also the size and the weight of the matching influence the distance above. For example, in Figure 3, matching M 1 gives the weighted adjacency graph with more components, but whose distance d σ (A M 1 , B M 1 ) = 1 + 5 − 2.7 = 3.3 is larger. On the other hand, M 2 gives the weighted adjacency graph with less components, but whose distance d σ (A M 2 , B M 2 ) = 2 + 5 − 3.9 = 3.1 is smaller.
Our goal in the following sections is to study the problem of computing the family-free DCJ distance, i.e., to find a matching in GS σ (A, B) that minimizes d σ . First of all, it is important to observe that the behaviour of this function does not correlate with the size of the matching. Often smaller matchings, that possibly discard gene assignments, lead to smaller distances. Actually, it is easy to see that, for any pair of genomes with any gene similarity graph, a trivial empty matching leads to the minimum distance, equal to zero. Due to this fact we restrict the distance to maximal matchings only. This ensures that no pairs of genes with positive similarity score are simply discarded, even though they might increase the overall distance. Hence we have the following optimization problem: Problem FFDCJ-DISTANCE(A, B): Given genomes A and B and their gene similarities σ , calculate their family-free DCJ distance where M is the set of all maximal matchings in GS σ (A, B).

Complexity of the family-free DCJ distance
In order to assess the complexity of FFDCJ-DISTANCE, we use a restricted version of the family-based exemplar DCJ distance problem [6,8]: Given genomes A and B, where each family occurs at most s times in A and at most t times in B, obtain exemplar genomes A and B by removing all but one copy of each family in each genome, so that the DCJ distance d DCJ (A , B ) is minimized.
We establish the computational complexity of the FFDCJ-DISTANCE problem by means of a polynomial time and approximation preserving (AP-) reduction from the problem (1, 2)-EXDCJ-DISTANCE, which is NP-hard [8].
Note that the authors of [8] only consider unichromosomal genomes, but the reduction can be extended to multichromosomal genomes, since an algorithm that solves the multichromosomal case also solves the unichromosomal case.

Theorem 1. Problem FFDCJ-DISTANCE(A, B) is APXhard, even if the maximum degrees in the two partitions of GS σ (A, B) are respectively one and two.
Before proving the result, we need some definitions and particularly a formal definition of an AP-reduction. These definitions are based on [16].
An optimization problem is defined by three main elements: a set of instances, a set Sol(I) of feasible solutions for each instance I, and a function val that relates a nonnegative rational number val(I, S) to each instance I and solution S in Sol(I). Thus, in a minimization problem, the aim is to find a feasible solution of minimum value. That is, if is an optimization problem with an instance I, then we want to find S ∈ Sol(I) that minimizes val(I, S), called an optimal solution to the optimization problem. For an instance I, the value of an optimal solution is denoted by opt(I).
An AP-reduction from an optimization problem to an optimization problem is a triple (f , g, β), where f and g are algorithms and β is a positive rational number, such that: An AP-reduction from to is frequently denoted by ≤ AP , and we say that is AP-reduced to . An APreduction is a special type of reduction which preserves both the polynomiality property and the approximation factor. Now, we can proceed with the proof of Theorem 1.
(AP1) Algorithm f receives as input a positive rational number δ and an instance (A, B) of (1, 2)-EXDCJ-DISTANCE where A and B are genomes from a set of genes G and each gene in G occurs at most once in A and at most twice in B, and constructs an instance , (A, B)) of FFDCJ-DISTANCE as follows. Let the genes of A be denoted a 1 , a 2 , . . . , a |A| and the genes of B be denoted b 1 , b 2 , . . . , b |B| . Then A F and B F are copies of A and B, respectively, except that symbol a i in A F is relabeled by i, keeping its orientation, and b j in B F is relabeled by j + |A|, also keeping its orientation. Furthermore, the normalized similarity measure σ for genes in A F and B F is defined as σ (i, k) = 1 for i in A F and k in B F , such that a i is in A, b j is in B, a i and b j are in the same gene family, and k = j + |A|. Otherwise, σ (i, k) = 0. Note that the construction is independent of the value of δ. into a solution (A X , B X ) of (1, 2)-EXDCJ-DISTANCE. This is a simple construction: for each edge (i, k) in M, we add symbols a i to A X and b j to B X , where j = k −|A|. The value of δ does not influence the construction. In the example of Figure 4, a matching M = {(1, 7), (2,8), (−3, −10), (4, 6)}, which is a solution to FFDCJ-DISTANCE(A F , B F ), is transformed by g into the genomes A X = { (• a 1 a 2 a 3 a 4   Let A X := A and B X be an exemplar genome of B, such that (A X , B X ) = g(δ, (A, B), M). We want to prove that (A X , B X ) is such that for some fixed positive rational number β. Denote by c AG and i AG the number of cycles and odd paths, respectively, in the adjacency graph AG(A X , B X ), and by c AG σ and i AG σ the number of cycles and odd paths, respectively, in the weighted adjacency graph AG σ (A M F , B M F ). Observe that the way the functions f and g have been defined, we have |A X | = |B X | = |M|, c AG = c AG σ , i AG = i AG σ , and thus Particularly, it is easy to see that we have opt(FFDCJ-DISTANCE(A F , B F )) = opt((1, 2)-EXDCJ-DISTANCE (A, B)) . Therefore, and Equation (3) holds by setting β := 1.

Corollary 2.
There exists no polynomial-time algorithm for FFDCJ-DISTANCE with approximation factor better than 1237/1236, unless P = NP.
Proof. As shown in [8], (1, 2)-EXDCJ-DISTANCE is NPhard to approximate within a factor of 1237/1236 − ε for any ε > 0. Therefore, the result follows immediately from [8] and from the AP-reduction in the proof of Theorem 1.
Since the weight plays an important role in d σ , a matching with maximum weight, that is obviously maximal, could be a candidate for the design of an approximation algorithm for FFDCJ-DISTANCE. However, we can demonstrate that it is not possible to obtain such an approximation, with the following example.
Consider an integer k ≥ 1 and let A = {(• 1 −2 · · · (2k−1) −2k •)} and B = {(• −(2k+1) (2k+ 2) · · · −(2k+2k−1) (2k+2k) •)} be two unichromosomal linear genomes. Observe that A and B have an even number of genes with alternating orientation. While A starts with a gene in direct orientation, B starts with a gene in reverse orientation. Now let σ be the normalized similarity measure between the genes of A and B, defined as follows: This shows that, for any genomes A and B, a matching of maximum weight in GS σ (A, B) can have d σ arbitrarily far from the optimal solution and cannot give an approximation for FFDCJ-DISTANCE(A, B).

ILP to compute the family-free DCJ distance
We propose an integer linear program (ILP) formulation to compute the family-free DCJ distance between two given genomes. This formulation is a slightly different version of the ILP for the maximum cycle decomposition problem given by Shao et al. [10] to compute the DCJ distance between two given genomes with duplicate genes. Besides the cycle decomposition in a graph, as was made in [10], we also have to take into account maximal matchings in the gene similarity graph and their weights.
Let A and B be two genomes with extremity sets X A and X B , respectively, and let G = GS σ (A, B) be their gene similarity graph. The weight w(e) of an edge e in G is also denoted by w e . Let M be a maximal matching in G. For the ILP formulation, a weighted adjacency graph and E(H) has three types of edges: (i) matching edges that connect two extremities in different extremity sets, one in X A and the other in X B , if there exists one edge in M connecting these genes in G; the set of matching edges is denoted by E m ; (ii) adjacency edges that connect two extremities in the same extremity set if they are an adjacency; the set of adjacency edges is denoted by E a ; and (iii) self edges that connect two extremities of the same gene in an extremity set; the set of self edges is denoted by E s . All edges in H are in E m ∪ E a ∪ E s = E(H). Matching edges have weights defined by the normalized similarity σ , all adjacency edges have weight 1, and all self edges have weight 0. Notice that any edge in G corresponds to two matching edges in H. Now we describe the ILP. For each edge e in H, we create the binary variable x e to indicate whether e will be in the final solution. We require first that each adjacency edge be chosen: We require then that, for each vertex in H, exactly one incident edge to it be chosen: Then, we require that the final solution be consistent, meaning that if one extremity of a gene in A is assigned to an extremity of a gene in B, then the other extremities of these two genes have to be assigned as well: We also require that the matching be maximal. This can easily be ensured if we guarantee that at least one of the vertices connected by an edge in the gene similarity graph be chosen, which is equivalent to not allowing both of the corresponding self edges in the weighted adjacency graph be chosen: To count the number of cycles, we use the same strategy as described in [10]. We first give an arbitrary index for We also require that adjacent vertices have the same label, forcing all vertices in the same cycle to have the same label: We create a binary variable z i , for each vertex v i , to verify whether y i is equal to its upper bound i: Since all the y i variables in the same cycle have the same label but a different upper bound, only one of the y i can be equal to its upper bound i. This means that for each cycle there can be only one z i equal to 1, and the sum of all z i variables is the total number of cycles in the adjacency graph.
In fact, it is possible to reduce the number of z i variables. First, notice that each cycle always has vertices from both genomes. That means that if we label all vertices v i starting with vertices of genome A first and then genome B, then the upper bounds for all y i s from genome A are smaller than the upper bounds for the y i s from genome B, and therefore no z i from genome B will ever be 1, since in the same cycle there will be at least one y i from genome A with a smaller upper bound. Then, all z i corresponding to vertices of genome B may be discarded: Finally, we set the objective function as follows: which is exactly the family-free DCJ distance d FFDCJ (A, B) as defined in Equations (1) and (2).

Simulations and experimental results
We performed some initial benchmarking experiments of the proposed ILP formulation. Therefore, we produced datasets using the Artificial Life Simulator (ALF) [17]. Genome sizes varied from 1000 to 3000 genes, where the gene lengths were generated according to a gamma distribution with shape parameter k = 3 and scale parameter θ = 133. A birth-death tree with 10 leaves was generated, with PAM distance of 100 from the root to the deepest leaf. For the amino acid evolution, the WAG substitution model with default parameters was used, with Zipfian indels at a rate of 0.000005. For structural evolution, gene duplications and gene losses were applied with a rate of 0.001 and reversals and translocations with a rate of 0.0025. To test different proportions of rearrangement events, we also simulated datasets where the structural evolution ratios had a 2-and 5-fold increase.
To solve the ILPs, we ran the CPLEX Optimizer on the 45 pairwise comparisons of each simulated dataset. All simulations were run in parallel on a cluster consisting of machines with an Intel(R) Xeon(R) E7540 CPU, with 48 cores and as many as 2 TB of memory, but for each individual CPLEX run only 4 cores and 2 GB of memory were allocated. The results are summarized in Table 1.

The family-free DCJ similarity
For a given matching M in GS σ (A, B), a formula for the similarity s σ of the reduced genomes A M and B M was first proposed in [12] only considering the cycles of AG σ (A M , B M ). Here we extend this formula to consider all components of the weighted adjacency graph. Again, let C, I, and P represent the sets of components in AG σ (A M , B M ) that are cycles, odd paths and even paths, respectively. Furthermore, w(C) = e∈C σ (e) is the sum of the weights of all the edges in a component C. Then the similarity s σ is the normalized total weight of all components: Here our goal is to study the problem of computing the family-free DCJ similarity, i.e., to find a matching in GS σ (A, B) that maximizes s σ . Similarly to the distance, the behaviour of the similarity does not correlate with the size of the matching. In other words, smaller matchings, that possibly discard gene assignments, can lead to higher similarities.
An approach for solving this problem was proposed in [12], following the one in [11] for gene adjacencies. It consists of a parameterized similarity function F α in which the user-controlled parameter α is a real number between 0 and 1: where, as above, w(M) = e∈M w(e) is the sum of the edge weights of the matching M.
Observe that the parameter α can be adjusted in favor of gene similarity when α is closer to 0, or in favor of genome organization similarity, when α is closer to 1. The closer Each dataset has 10 genomes, totalling 45 pairwise comparisons. Maximum running time was set to 60 minutes. For each dataset, the number of runs is shown that found an optimal solution within the allowed time and their average running time in seconds. For the runs that did not finish, the last row shows the relative gap between the upper bound and the current solution. Rate r = 1 means the default rate for ALF evolution, and r = 2 and r = 5 mean 2-fold and 5-fold increase for the gene duplication, gene deletion and rearrangement rates.
the parameter α is to 0, the closer we are to the problem of finding a maximum weighted matching in the gene similarity graph GS σ (A, B). On the other hand, the closer α is to 1, the closer we are to the problem of computing s σ (A M , B M ). A drawback of this model is that the weights of edges actually appear in both terms of the equation. Furthermore, it remains the problem of finding the "best" value for α.
Here, instead of adopting the parameter α, we restrict the similarity to maximal matchings only, ensuring that no pair of genes with positive similarity score is simply discarded, even though it might decrease the overall similarity. We then have the following optimization problem:

Complexity of the family-free DCJ similarity
We have the following result to the family-free DCJ similarity. Proof. We use the Cook reduction, which is a polynomial time transformation, from (1, 2)-EXDCJ-DISTANCE to FFDCJ-SIMILARITY.
Let A and B be any instance of (1, 2)-EXDCJ-DISTANCE and let k be a positive integer, with k ≤ |A|, where |A| is the number of genes of a genome A. We suppose, without loss of generality, that A and B are circular multichromosomal genomes. We must construct a pair of circular genomes A F and B F , a normalized similarity measure σ for genes in A F and B F , and a positive integer k ≤ |A F | such that the family-free DCJ similarity of A F and B F is at least k if and only if the exemplar DCJ distance of genomes A and B is at most k.
The construction of A F , B F , σ , and k is similar to the transformation f in (AP1) of the proof of Theorem 1. Let G be the underlying gene set, such that each gene in G occurs at most once in A and at most twice in B. Let the genes of A be denoted a 1 , a 2 , . . . , a |A| and the genes of B be denoted b 1 , b 2 , . . . , b |B| . Then A F and B F are copies of A and B, respectively, except that symbol a i in A F is relabeled by i, keeping its orientation, and b j in B F is relabeled by j + |A|, also keeping its orientation. The normalized similarity measure σ for genes in A F and B F is defined as σ (i, k) = 1 for i in A F and k in B F , such that a i is in A, b j is in B, a i and b j are in the same gene family, and k = j + |A|. Otherwise, σ (i, k) = 0. It is easy to see that this construction can be accomplished in poynomial time.
Now we must show that the family-free DCJ similarity of A F and B F is at least k if and only if the exemplar DCJ distance of genomes A and B is at most k. Let n = |A|.
Suppose first that M is a matching in the gene similarity graph

Conclusion
In this paper, we have defined a new distance measure for two genomes that is motivated by the double cut and join model, while not relying on gene annotations in form of gene families. In case gene families are known and each family has exactly one member in each of the two genomes, this distance equals the family-based DCJ distance and thus can be computed in linear time. In the general case, however, it is NP-hard and even hard to approximate. Nevertheless, we could give an integer linear program for the exact computation of the distance that is fast enough to be applied to realistic problem instances. Similar theoretical results hold for the family-free DCJ similarity measure, which is NP-hard.
The family-free model has many potentials when gene family assignments are not available or ambiguous, in fact it can even be used to improve family assignments [18]. The work presented in this paper is another step in this direction.