Computing the Rearrangement Distance of Natural Genomes

Abstract

The computation of genomic distances has been a very active field of computational comparative genomics over the past 25 years. Substantial results include the polynomial-time computability of the inversion distance by Hannenhalli and Pevzner in 1995 and the introduction of the double cut and join distance by Yancopoulos et al. in 2005. Both results, however, rely on the assumption that the genomes under comparison contain the same set of unique markers (syntenic genomic regions, sometimes also referred to as genes). In 2015, Shao et al. relax this condition by allowing for duplicate markers in the analysis. This generalized version of the genomic distance problem is NP-hard, and they give an integer linear programming (ILP) solution that is efficient enough to be applied to real-world datasets. A restriction of their approach is that it can be applied only to balanced genomes that have equal numbers of duplicates of any marker. Therefore, it still needs a delicate preprocessing of the input data in which excessive copies of unbalanced markers have to be removed. In this article, we present an algorithm solving the genomic distance problem for natural genomes, in which any marker may occur an arbitrary number of times. Our method is based on a new graph data structure, the multi-relational diagram, that allows an elegant extension of the ILP by Shao et al. to count runs of markers that are under- or over-represented in one genome with respect to the other and need to be inserted or deleted, respectively. With this extension, previous restrictions on the genome configurations are lifted, for the first time enabling an uncompromising rearrangement analysis. Any marker sequence can directly be used for the distance calculation. The evaluation of our approach shows that it can be used to analyze genomes with up to a few 10,000 markers, which we demonstrate on simulated and real data.

Keywords: DCJ-indel distance; comparative genomics; genome rearrangements.

FIG. 1.

For genomes A={1 ¯6 5 3, 4 2} and B={1 7 2 3 4 5, 7 ¯8}, the relational diagram contains one cycle, two AB-paths (represented in blue), one AA-path, and one BB-path (both represented in red). Short dotted horizontal edges are adjacency edges, long horizontal edges are indel edges, and top–down edges are extremity edges.

FIG. 2.

(i) A BB-path with four runs. (ii) After an optimal DCJ that creates a new cycle, one $A$-run is accumulated (between edges e₄ and e₃ there is only an adjacency edge) and two $ℬ$-runs are merged (e₂ is in the same run with e₅ and e₆). Indeed, the indel-potential of the original BB-path is three.

FIG. 3.

An optimal recombination with $\Delta d=\Delta \lambda =-2$.

FIG. 4.

Chained recombinations transforming four paths ($2\times A{A}_{Aℬ}+B{B}_A+B{B}_{ℬ}$) into four other paths ($3\times A{B}_{\epsilon }+A{B}_{ℬ}$) with overall $\Delta d=-3$.

FIG. 5.

Natural genomes A=1 3 2 ¯5 ¯4 3 5 4 and $B=1623173413$ can give rise to many distinct pairs of matched singular genomes. The relational diagrams of two of these pairs are represented here, in the top and center. In the bottom, we show the multi-relational diagram $MR\left(A,B\right)$. The decomposition that gives the diagram in the top is represented in red/orange. Similarly, the decomposition that gives the diagram in the center is represented in blue/cyan. Edges that are in both decompositions have two colors.

FIG. 6.

Optimal capping of canonical genomes $A=\left\{21,43\right\}$ and $B=\left\{12,34\right\}$ into $A_{\circ }=\left\{\left(215\right),\left(436\right)\right\}$, and $B_{\circ }=\left\{\left(125\right),\left(346\right)\right\}$. Each pair of AA-+BB-path is linked into a separate AB-cycle.

FIG. 7.

Optimal capping of singular genomes $A=\left\{521,5453\right\}$ and $B=\left\{6162,364\right\}$ into and $B_{\circ }=\left\{\left(761628364\right)\right\}$. This capping shows how to optimally link the four sources of the chained recombinations of Figure 4 into a single AB-cycle.

FIG. 8.

Natural genomes A=1 3 2 ¯5 ¯4 3 5 4 and $B=1623173413$ and their capped multi-relational diagram $M{R}_{\circ }\left(A,B\right)$.

FIG. 9.

Solving times for: (i) genomes with a varying number of duplicate occurrences, totaling 20,000 marker occurrences per genome; (ii) genome pairs with a varying number of linear chromosomes with 20,000 marker occurrences per genome; (iii) varying number of DCJs and indels applied by the simulation to genomes of 35,000 marker occurrences; and (iv) genome pairs with a varying total number of marker occurrences from both genomes.

FIG. 10.

The gene-based distances in Table 9 are used as input to reconstruct Drosophila-Phylogeny: (i) with Neighbor Joining; and (ii) as a Splits diagram.

FIG. 11.

The segmentation-based distances in Table 10 are used as input to reconstruct Drosophila-Phylogeny: (i) with Neighbor Joining; and (ii) as a Splits diagram.