==== Front Nat Commun Nat Commun Nature Communications 2041-1723 Nature Publishing Group UK London 19777 10.1038/s41467-020-19777-8 Article Scalable multiple whole-genome alignment and locally collinear block construction with SibeliaZ http://orcid.org/0000-0002-0807-347XMinkin Ilia ivminkin@gmail.com 1 http://orcid.org/0000-0003-3143-594XMedvedev Paul 123 1 grid.29857.310000 0001 2097 4281Department of Computer Science and Engineering, The Pennsylvania State University, 506 Wartik Lab University Park, University Park, PA 16802 USA 2 grid.29857.310000 0001 2097 4281Department of Biochemistry and Molecular Biology, The Pennsylvania State University, 506 Wartik Lab University Park, University Park, PA 16802 USA 3 grid.29857.310000 0001 2097 4281Center for Computational Biology and Bioinformatics, The Pennsylvania State University, 506 Wartik Lab University Park, University Park, PA 16802 USA 10 12 2020 10 12 2020 2020 11 632728 9 2019 29 10 2020 © The Author(s) 2020Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.Multiple whole-genome alignment is a challenging problem in bioinformatics. Despite many successes, current methods are not able to keep up with the growing number, length, and complexity of assembled genomes, especially when computational resources are limited. Approaches based on compacted de Bruijn graphs to identify and extend anchors into locally collinear blocks have potential for scalability, but current methods do not scale to mammalian genomes. We present an algorithm, SibeliaZ-LCB, for identifying collinear blocks in closely related genomes based on analysis of the de Bruijn graph. We further incorporate this into a multiple whole-genome alignment pipeline called SibeliaZ. SibeliaZ shows run-time improvements over other methods while maintaining accuracy. On sixteen recently-assembled strains of mice, SibeliaZ runs in under 16 hours on a single machine, while other tools did not run to completion for eight mice within a week. SibeliaZ makes a significant step towards improving scalability of multiple whole-genome alignment and collinear block reconstruction algorithms on a single machine. Multiple whole-genome alignment is a challenging problem in bioinformatics, especially when computational resources are limited. Here the authors present SibeliaZ, an algorithm and software based on analysis of de Bruijn graphs, which provides improved computational efficiency and scalability. Subject terms Comparative genomicsComputational biology and bioinformaticshttps://doi.org/10.13039/100000001National Science Foundation (NSF)CCF-1439057IIS-1453527DBI-1356529IIS-1421908Medvedev Paul https://doi.org/10.13039/100000057U.S. Department of Health & Human Services | NIH | National Institute of General Medical Sciences (NIGMS)R01GM130691Medvedev Paul issue-copyright-statement© The Author(s) 2020 ==== Body Introduction Multiple whole-genome alignments are the problem of identifying all the high-quality multiple local alignments within a collection of assembled genome sequences. It is a fundamental problem in bioinformatics and forms the starting point for most comparative genomics studies, such as rearrangement analysis, phylogeny reconstruction, and the investigation of evolutionary processes. Unfortunately, the presence of high-copy repeats and the sheer size of the input make multiple whole-genome alignment extremely difficult. While current approaches have been successfully applied in many studies, they are not able to keep up with the growing number and size of assembled genomes1. The multiple whole-genome alignment problem is also closely related to the synteny reconstruction problem and to the questions of how to best represent pan-genomes. There are two common strategies to tackle the whole-genome alignment problem2. The first one is based on finding pairwise local alignments3–7 and then extending them into multiple local alignments8–11. While this strategy is known for its high accuracy, a competitive assessment of multiple whole-genome alignment methods1 highlighted several limitations. First, many algorithms either do not handle repeats by design or scale poorly in their presence, since the number of pairwise local alignments grows quadratically as a function of a repeat’s copy number. In addition, many algorithms use a repeat database to mask high-frequency repeats. However, these databases are usually incomplete and even a small amount of unmasked repeats may severely degrade alignment performance. Second, the number of pairwise alignments is quadratic in the number of genomes, and only a few existing approaches could handle more than ten fruit fly genomes1. Therefore, these approaches are ill-suited for large numbers of long and complex genomes, such as mammalian genomes in general and the recently assembled 16 strains of mice12 in particular. Alternatively, anchor-based strategies can be applied to decompose genomes into locally collinear blocks13. These are blocks that are free from nonlinear rearrangements, such as inversions or transpositions. Once such blocks are identified, they can independently be global aligned13–17. The problem of constructing blocks from anchors is known as the chaining problem which had been extensively studied in the past18–20. All of the methods applicable to datasets consisting of multiple genomes are heuristic since the exact algorithms depend exponentially on the number of genomes. Such strategies are generally better at scaling to handle repeats and multiple genomes since they do not rely on the computationally expensive pairwise alignment. A promising strategy to find collinear blocks is based on the compacted de Bruijn graph21–23 widely used in genome assembly. Though these approaches do not work well for divergent genomes, they remain fairly accurate for closely related genomes. For example, Sibelia23 can handle repeats and works for many bacterial genomes; unfortunately, it does not scale to longer genomes. However, the last three years has seen a breakthrough in the efficiency of de Bruijn graph construction algorithms24–27. The latest methods can construct the graph for tens of mammalian genomes in minutes rather than weeks. We therefore believe the de Bruijn graph approach holds the most potential for enabling scalable multiple whole-genome alignments of closely related genomes. In this paper, we describe an algorithm SibeliaZ-LCB for identifying collinear blocks in closely related genomes. SibeliaZ-LCB is suitable for detecting homologous sequences that have an evolutionary distance to the most recent common ancestor (MRCA) of at most 0.09 substitutions per site. SibeliaZ-LCB is based on the analysis of the compacted de Bruijn graph and uses a graph model of collinear blocks similar to the "most frequent paths” introduced by28. This allows it to maintain a simple, static, data structure, which scales easily and allows simple parallelization. Thus, SibeliaZ-LCB overcomes a bottleneck of previous state-of-the-art de Bruijn graph-based approaches17,22, which relied on a dynamic data structure which was expensive to update. Further, we extend SibeliaZ-LCB into a multiple whole-genome aligners called SibeliaZ. SibeliaZ works by first constructing the compacted de Bruijn graph using our previously published TwoPaCo tool27, then finding locally collinear blocks using SibeliaZ-LCB, and finally, running a multiple-sequence aligner spoa29 on each of the foundation blocks. To demonstrate the scalability and accuracy of our method, we compute the multiple whole-genome alignment for a collection of recently assembled strains of mice. We also test how our method works under different conditions, including various levels of divergence between genomes and different parameter settings. Our software is freely available at https://github.com/medvedevgroup/SibeliaZ/. Results Algorithm overview As described in the introduction, the major algorithmic innovation of this paper is the SibeliaZ-LCB algorithm. SibeliaZ-LCB takes as input a de Bruijn graph built on a collection of assembled genomes. An assembled genome is itself a set of contig sequences. SibeliaZ-LCB identifies and outputs all nonoverlapping blocks of homologous subsequences of the input genomes. A block can be composed of two or more sequences from one or more genomes. In this subsection, we will give a high level overview of SibeliaZ-LCB, leaving the more formal and detailed version for the “Methods”. SibeliaZ-LCB relies heavily on the de Bruijn graph of the genomes. In this graph, the vertices correspond to the k-mers (substrings of fixed length k) of the input. A k-mer that appears multiple times in the input is represented using just one node. Then, k-mers that appear consecutively in some input sequence are connected by an edge from the left one to the right one (see Fig. 1a for an example). This way, each genome corresponds to a path in the graph that hops from k-mer to k-mer using the edges.Fig. 1 The de Bruijn graph and an example of a collinear block. a The graph built from strings “GCACGTCC” and “GCACTTCC”, with k = 2. The two strings are reflected by the blue and magenta walks, respectively. This is an example of a collinear block from two walks. There are four bubbles. The bubble formed by vertices “AC” and “TC” describes a substitution within the block, while three other bubbles are formed by parallel edges. The blue and magenta walk form a chain of four consecutive bubbles. b An example of a more complex block, where we have added a third sequence “CACGTTCC” (turquoise) to the input. We can no longer describe the block as a chain of bubbles, as they overlap to form tangled structures. Instead, we consider the path in the graph (dashed black) that shares many vertices with the three collinear walks. This carrying path shares many vertices with the three extant walks, and each walk forms its own chain with it. The task of finding good collinear blocks can then be framed as finding carrying paths that form good chains with the genomic walks. In this graph, two homologous sequences form what is called a chain: an interleaving sequence of parallel edges, which correspond to identical sequences, and “bubbles”, which correspond to small variations like single nucleotide variants or indels. However, the concept of a chain is difficult to extend to more than two homologous sequences because the tangled pattern in the graph is difficult to precisely define (see Fig. 1b for an example). To address this challenge, we introduce the idea that each block has a “carrying path” in the de Bruijn graph that holds the block together. The basic idea is that the homologous sequences forming the block have a lot of shared k-mers and their corresponding genomic paths go through nearly the same vertices. A carrying path is then a path that goes through the most frequently visited vertices, loosely similar to the notion of a consensus sequence from alignment. Each genomic path from the block then forms a chain with this carrying path (see Fig. 1b for an example). We do not know the carrying paths in advance but we can use them as a guiding mechanism to find blocks. We start with an arbitrary edge e in the graph and all other genomic paths that form bubbles with e. We make e the starting point of a carrying path and use it along with the other genomic paths to initiate the collection of sequences making up the block corresponding to this carrying path. To extend the carrying path, we look at the edges extending the genomic paths in the current block and take the most common one. The data structures maintaining the genomic paths in the block and the carrying path are then updated and the extension procedure repeats. Figure 2 shows an example of running this algorithm.Fig. 2 An example of running Algorithm Find-collinear-blocks (Box 1) on the graph from Fig. 1b, starting from edge GC  →  CC as the seed. Each subfigure shows the content of the collinear block P and the carrying path. The collinear walks are solid, the carrying path is dashed, and the rest of the graph is dotted. Subfigure a shows the state of these variables after the initialization; subfigures b–d show the state after the completion of each phase. We continue this process until the scoring function that describes how well a carrying path holds the block together falls below zero. At that point, we consider the possibility that we may have overextended the block and should have instead ended it earlier. To do this, we look at all the intermediate blocks we had created during the extension process and output the one that has the highest score. Once a block is an output, we output all its constituent edges as used so that they are not chosen as part of a future block. In this way, SibeliaZ-LCB finds a single block. Afterwards, we try to find another block by starting from another arbitrary edge. This process continues until all the edges in the graph are either used or had been tried as potential starters for a carrying path. Datasets, tools, and evaluation metrics Evaluation of multiple whole-genome aligners is a challenging problem in its own right and we, therefore, chose to use the practices outlined in the Alignathon1 competition as a starting point. They present several approaches to assess the quality of a multiple whole-genome alignments. Ideally, it is best to compare an alignment against a manually curated gold standard; unfortunately, such a gold standard does not exist. We, therefore, chose to focus our evaluation on real data. We evaluated the ability of SibeliaZ to align real genomes by running it on several datasets consisting of a varying numbers of mice genomes. We retrieved 16 mice genomes available at GenBank30 and labeled as having a “chromosome” level of assembly. They consist of the mouse reference genome and 15 different strains assembled as part of a recent study12 (Supplementary Table 1). The genomes vary in size from 2.6 to 2.8 Gbp and the number of scaffolds (between 2977 and 7154, except for the reference, which has 377). Their GenBank accession numbers are listed in Table 1. We constructed four datasets of increasing size to test the scalability of the pipelines with respect to the number of input genomes. The datasets contain genomes 1–2, 1–4, 1–8, and 1–16 from Supplementary Table 1, with genome 1 being the reference genome.Table 1 Accession numbers of the assembled mice genomes available at GenBank. Strain Accession number C57BL/6J GCA_000001635.8 129S1/SvImJ GCA_001624185.1 A/J GCA_001624215.1 AKR/J GCA_001624295.1 CAST/EiJ GCA_001624445.1 CBA/J GCA_001624475.1 DBA/2J GCA_001624505.1 FVB/NJ GCA_001624535.1 NOD/ShiLtJ GCA_001624675.1 NZO/HiLtJ GCA_001624745.1 PWK/PhJ GCA_001624775.1 WSB/EiJ GCA_001624835.1 BALB/cJ GCA_001632525.1 C57BL/6NJ GCA_001632555.1 C3H/HeJ GCA_001632575.1 LP/J GCA_001632615.1 To measure accuracy, we used several ground-truth alignments (to be described) and employed the metrics of precision and recall used in the Alignathon and implemented by the mafTools package1. For these metrics, alignment is viewed as an equivalence relation. We say that two positions in the input genomes are equivalent if they originate from the same position in the genome of their recent common ancestor. We denote by H the set of all equivalent position pairs, participating in the “true” alignment. Let A denote the relation produced by an alignment algorithm. The accuracy of the alignment is then given by recall(A) = 1 − ∣H⧹A∣/∣H∣ and precision(A) = 1 − ∣A⧹H∣/∣A∣, where ⧹ denotes set difference. To evaluate recall, we compared our results against annotations of protein-coding genes. We retrieved all pairs of homologous protein-coding gene sequences from Ensembl and then computed pairwise global alignments between them using LAGAN31. The alignment contains both orthologous and paralogous genes, though most of the paralogous pairs come from the well-annotated mouse reference genome. We removed any pairs of paralogous genes with overlapping coordinates, as these were likely mis-annotations, as confirmed by Ensembl helpdesk32. We made these filtered alignments as well as the alignments produced by SibeliaZ available for public download from our GitHub repository (see Section “Data availability” for the links). We define the nucleotide identity of an alignment as the number of matched nucleotides divided by the length of an alignment, including gaps. The distribution of nucleotide identities, as well as the coverage of the annotation, is shown in Supplementary Fig. 1. In our analysis, we binned pairs of genes according to their nucleotide identity. Since protein-coding genes only compromise a small portion of the genome, we also computed all-against-all pairwise local alignments between chromosomes 1 of genomes 1–2 and 1–4 using LASTZ6, a reliable local aligner known for its accuracy. We only computed alignments between chromosomes of different genomes, i.e., did not include self-alignments, which excludes duplications such as paralogous genes from the alignment. We used default settings of LASTZ except that we made it output alignments of nucleotide identity at least 90%. We then evaluated the recall and precision of our alignments but restricted our alignments to the sequences of chromosome 1. We then treated the LASTZ alignments as the ground truth. The LASTZ alignments are available for download from our repository’s supplemental data section. Note that because the alignment is represented as a set of positions pairs, it is possible to evaluate a multiple whole-genome alignments using pairwise local alignments. To measure precision, we use the LASTZ alignments on chromosome 1. However, it is computationally prohibitive to compute such alignments with LASTZ for the whole genome. We therefore also use an indirect way to assess precision for the whole genome. For each column in the alignment, we calculate the average number of nucleotide differences33. In an alignment of highly similar genomes that has high precision, we expect these numbers to below (close to 0) for most of the columns in the alignment. Otherwise, it would suggest the presence of unreliable poorly aligned blocks in the alignment. Formally, given a column c of a multiple whole-genome alignments with ci being its ith element, average number of nucleotide differences is given by π(c)=∑1≤i≤∣c∣∑ibor∣q3∣>b). The third case forbids walks (i.e., gives them a score of  −∞) where the hanging ends are too long, and the first case ignores walks (i.e., gives them a score of 0) that weave through pa but are too short. The second case gives a score that is proportional to the length of the part of pa that forms a chain with p. At the same time, it reduces the score if the collinear walks leave hanging ends q1 and q3—the parts of pa not participating in the chain. The penalty induced by these ends is squared to remove spuriously similar sequences from from the collinear block. This form of scoring function was chosen because it performed well during the initial stages of development. We do not penalize for the discrepancy between p and q2 for the sake of simplicity of the scoring function and avoiding extra computation needed to calculate it. Figure 5 shows an example of computing the score.Fig. 5 An example of computing the score of a walk p (solid) relative to a carrying path pa = q1q2q3 (dashed). The path p forms a chain with the subpath q2 of pa, while subpaths q1 and q3 form hanging ends. We count the length of p and subtract lengths of the hanging ends. Thus, the score f(pa, p) = 4 − (1 + 1)2 = 0. The collinearity score of a collinear block is given by 2 f(P)=maxpa ∑p∈Pf(pa,p), where pa can be any path (not necessarily genomic). In other words, we are looking for a path forming longest chains with the collinear walks and thus maximizes the score. The collinear blocks reconstruction problem is to find a set of collinear blocks P such that ∑P∈Pf(P) is maximum and no two walks in P share an edge. Note that the number of collinear blocks is not known in advance. For an example of a complex collinear block in the de Bruijn graph and a carrying path capturing it, refer to Fig. 1b. The collinear blocks reconstruction algorithm Our algorithm’s main pseudocode is shown in Box 1 and its helper function in Box 2. First, we describe the high-level strategy. The main algorithm is greedy and works in the seed-and-extend fashion. It starts with an arbitrary edge in the graph and tries to extend it into a carrying path that induces a collinear block with the highest possible collinearity score Pbest. If the block has a positive score, then Pbest is added to our collection of collinear blocks P. The algorithm then repeats, attempting to build a collinear block from a different edge seed. New collinear blocks cannot use edges belonging to previously discovered collinear blocks. This process continues until all possible edges are considered as seeds. The algorithm is greedy in the sense that once a block is found and added to P, it cannot be later changed to form a more optimal global solution. To extend a seed into a collinear block P, we first initialize the collinear block with a walk for each unused edge parallel to the seed (including the seed) (lines 7 and 8). These parallel edges represent the different occurrences of the seed string in the input and, hence, form the initial collinear block. We then proceed in phases, where each phase is an iteration of the while loop (lines 9–19). During each phase, the carrying path pa is extended using a walk r of length at most b (lines 10–14). Next, we try to extend each of the collinear walks in a way that forms chains with the extended pa (lines 15–19). The extension of a seed into a collinear block is also a greedy process since we only change pa and the walks in P by extending them and never by changing any edges. Finally, we check that the collinearity score for our extended block is still positive—if it is, we iterate to extend it further, otherwise, we abandon our attempts at further extending the block. We then recall the highest-scoring block that was achieved for this seed and save it into our final result P (lines 20–22). To pick the walk r by which to extend pa, we use a greedy heuristic (lines 10–14). First, we pick the vertex t which we want to extend to reach (lines 10–12). We limit our search to those vertices that can be reached by a genomic walk from the end of pa and greedily chose the one that is most often visited by the b-extensions of the collinear walks in P. Intuitively, we hope to maximize the number of collinear walks that will form longer chains with pa after its extension and thereby boost the collinearity score. We then extend pa using the shortest b-extension of the walks in P to reach t. We chose this particular heuristic because it showed superior performance compared to other possible strategies. Once we have selected the genomic walk r by which to extend pa, we must select the extensions to our collinear walks P that will form chains with par. This is done by the function Update-collinear-blocks (Box 2). We extend the walks to match r by considering the vertices of r consecutively, one at a time. To extend to a vertex w, we consider all the different locations of w in the input (each such location is represented by an edge x ending at w). For each location, we check if it can be reached by a b-extension from an existing p ∈ P. If yes, then we extend p, so as to lengthen the chain that it forms with pa. If there are multiple collinear walks that reach w, we take the nearest one. If no, then we start a new collinear walk using just x. Figure 6 shows an example of updating a collinear walk and Fig. 2 shows a full run of the algorithm for a single seed.Fig. 6 Illustration for Algorithm Update-collinear-walks (Box 2). A collinear walk p (solid) requires an update after the carrying path pa is extended with the dashed edge (w0, w). The path pa now ends at the vertex w, which has another incoming edge x. Since x is a part of the b-extension of p (denoted by q), p can be appended with q to form a longer chain and boost the collinearity score. Our description here only considers extending the initial seed to the right, i.e., using out-going edges in the graph. However, we also run the procedure to extend the initial seed to the left, using the incoming edges. The case is symmetric and we, therefore, omit the details. Box 1 Algorithm Find-collinear-blocks Input: strings S, integers k, b and m Output: a set of edge-disjoint subgraphs of G(S, k) representing collinear blocks 1: P←∅ ⊳ Collinear blocks 2: G ← G(S, k)         ⊳ Construct the multigraph 3: for all distinct pairs (u, v) ∈ E(G) do        ⊳ Check possible seeds 4:     Initialize the current-carrying path pa with (u, v) 5:     P←∅ ⊳ Sorted set of collinear walks forming chains with pa 6:     Pbest←∅ ⊳ Highest-scoring collinear block induced by pa 7:     for edges x ∈ E(G) parallel to (u, v) not marked as used do 8:        Add to P a new collinear walk consisting of x 9:     while f(P) ≥ 0 do        ⊳ Extend the carrying path as far as possible 10:        Q ← {q ∣ q is the b − extension of a p ∈ P} 11:        w0 ← last vertex in pa 12:        t ← a vertex, reachable from w0 via a genomic walk, that is visited by the most walks of Q. 13:        Let r ∈ Q be the shortest walk from w0 to t 14:        Denote the vertices of r as w0, w1, …, w∣r∣, with w∣r∣ = t 15:        for i ← 1 to ∣r∣ do 16: Append (wi−1, wi) to the carrying path pa 17: P ← Update-collinear-walks(P, wi) 18: if f(P) > f(Pbest) then 19: Pbest ← P 20:     if f(Pbest) > 0 then 21:        P←P∪{Pbest} 22:        Mark edges visited by walks of Pbest as used 23: return P Box 2 Update-collinear-walks Input: A sorted set of collinear walks P, a vertex w Output: Updated set P 1: for edges x ∈ E(G) ending at w not marked as used do 2:      Let p ∈ P be a walk such that its b-extension q contains x and pos(end(p)) is maximized         ⊳ Find a walk extendable with x 3:      if such p exists then 4:         Truncate q so that end(q) = x 5:         Append p with q         ⊳ Lengthen the chain that p forms with pa 6:     else 7:        Add a new walk consisting of the edge x to P 8: return P Other considerations For simplicity of presentation, we have described the algorithm in terms of the ordinary de Bruijn graph; however, it is crucial for running time and memory usage that the graph is compacted first. Informally, the compacted de Bruijn graph replaces each non-branching path with a single edge. Formally, the vertex set of the compacted graph consists of vertices of the regular de Bruijn graph that have at least two outgoing (or ingoing) edges pointing at (incoming from) different vertices. Such vertices are called junctions. Let ℓ = v1, …, vn be the list of k-mers corresponding to junctions, in the order, they appear in the underlying string s. The edge set of the compacted graph consists of edges {v1 → v2, v2 → v3, …, vn−1 → vn}. We efficiently construct the compacted graph using our previously published algorithm TwoPaCo27. This transformation maintains all the information while greatly reducing the number of edges and vertices in the graph. This makes the data structures smaller and allows the algorithm to fast-forward through non-branching paths, instead of considering each (k + 1)-mer one by one. Our previous description of the algorithm remains valid, except that the data structures operate with vertices and edges from the compacted graph instead of the ordinary one. The only necessary change is that when we look for an edge y parallel to x, we must also check that y and x spell the same sequence. This is always true in an ordinary graph but not necessarily in a compacted graph. An important challenge of mammalian genomes is that they contain high-frequency (k + 1)-mers, which can clog up our data structures. To handle this, we modify the algorithm by skipping over any junctions that correspond to k-mers occurring more than a times; we call a the abundance pruning parameter. Specifically, prior to constructing the edge set of the compacted de Bruijn graph, we remove all high abundance junctions from the vertex set. The edge set is constructed as before, but using this restricted list of junctions as the starting point. This strategy offers a way to handle high-frequency repeats at the expense of limiting our ability to detect homologous blocks that occur more than a times. The organization of our data in memory is instrumental to achieving high performance. To represent the graph, we use a standard adjacency list representation, annotated with position information, and other relevant data. We also maintain a list of the junctions in the paragraph above in the order they appear in the input sequences, thereby supporting next() queries. The walks in the collinear block P are stored as a dynamic sorted set, implemented as a binary search tree. The search key is the genome/position for the end of each walk. This allows performing a binary search in line 2 of Algorithm Update-collinear-walks. Another aspect that we have ignored up until now is that DNA is double-stranded and collinear walks can be reverse-complements of each other. If s is a string, then let s¯ be its reverse complement. We handle double strandedness in a natural way by using the comprehensive de Bruijn graph, which is defined as Gcomp(s,k)=G(s,k)∪G(s¯,k)27. Our algorithm and corresponding data structures can be modified to work with the comprehensive graph with a few minor changes which we omit here. Our implementation is parallelized by exploring multiple seeds simultaneously, i.e., parallelizing the for a loop at line 3 of Algorithm Find-collinear-blocks. This loop is not embarrassingly parallelizable, since two threads can start exploring two seeds belonging to the same carrying path. In such a case, there will be a collision on the data structure used to store used marks. To address this issue, we process the seeds in batches of fixed size. All the seeds within a batch are explored in parallel and the results are saved without modifying the “used” marks. Once the batch is processed, a single arbiter thread checks if there is any overlap in the used marks of the different threads. If there is, it identifies the sources of the conflict and reruns the algorithm at the conflicting seeds serially. Since most seeds do not yield valid carrying paths, such conflicts are rare. Once there is no conflict, the arbiter updates the used main data structures with the results of the batch. This design allows the computation result to be deterministic and independent of the number of threads used. Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article. Supplementary information Supplementary Information Reporting Summary Peer review information Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Supplementary information Supplementary information is available for this paper at 10.1038/s41467-020-19777-8. Acknowledgements We would like to thank Mikhail Kolmogorov for useful suggestions on the empirical evaluation of our algorithm; Robert Harris for his help with running MultiZ; Son Pham for introducing us to the problem; and the Ensembl support team for helping us with retrieving the gene annotations. This work has been supported in part by NSF awards DBI-1356529, CCF-1439057, IIS-1453527, and IIS-1421908 to P.M. The research reported in this publication was supported by the National Institute Of General Medical Sciences of the National Institutes of Health under Award No. R01GM130691. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. Author contributions Conceptualization: I.M., Methodology: I.M. and P.M., Software: I.M., Validation: I.M. and P.M., Writing—original draft: I.M. and P.M., Writing—review and editing: I.M. and P.M., and Funding acquisition: P.M. Data availability Table 1 contains the list of GenBank accession numbers of the mice genomes we used in our experiments (Figs. 3 and 4, Supplementary Figs. 1–3). The nine simulated datasets we generated (Supplementary Figs. 4–6), ground-truth alignments for the mouse data (Fig. 4, Supplementary Figs. 1–3), and alignments produced by SibeliaZ and Progressive Cactus (Fig. 4, Supplementary Figs. 2 and 3) are available for download at https://github.com/medvedevgroup/SibeliaZ/blob/master/DATA.txt. Code availability Our tool is open source and freely available at https://github.com/medvedevgroup/SibeliaZ. Competing interests The authors declare no competing interests. ==== Refs References 1. Earl D Alignathon: a competitive assessment of whole-genome alignment methods Genome Res. 2014 24 2077 2089 10.1101/gr.174920.114 25273068 2. Dewey CN Pachter L Evolution at the nucleotide level: the problem of multiple whole-genome alignment Hum. Mol. Genet. 2006 15 R51 R56 10.1093/hmg/ddl056 16651369 3. Altschul SF Gish W Miller W Myers EW Lipman DJ Basic local alignment search tool J. Mol. Biol. 1990 215 403 410 10.1016/S0022-2836(05)80360-2 2231712 4. Altschul SF Gapped blast and psi-blast: a new generation of protein database search programs Nucleic Acids Res. 1997 25 3389 3402 10.1093/nar/25.17.3389 9254694 5. Schwartz S Human–mouse alignments with blastz Genome Res. 2003 13 103 107 10.1101/gr.809403 12529312 6. Harris, R. S. Improved Pairwise Alignment of Genomic DNA. (The Pennsylvania State University, 2007). 7. Kent WJ Blat—the blast-like alignment tool Genome Res. 2002 12 656 664 10.1101/gr.229202 11932250 8. Blanchette M Aligning multiple genomic sequences with the threaded blockset aligner Genome Res. 2004 14 708 715 10.1101/gr.1933104 15060014 9. Dubchak I Poliakov A Kislyuk A Brudno M Multiple whole-genome alignments without a reference organism Genome Res. 2009 19 682 689 10.1101/gr.081778.108 19176791 10. Angiuoli SV Salzberg SL Mugsy: fast multiple alignment of closely related whole genomes Bioinformatics 2011 27 334 342 10.1093/bioinformatics/btq665 21148543 11. Paten B Cactus: algorithms for genome multiple sequence alignment Genome Res. 2011 21 1512 1528 10.1101/gr.123356.111 21665927 12. Lilue J Sixteen diverse laboratory mouse reference genomes define strain-specific haplotypes and novel functional loci Nat. Genet. 2018 50 1574 10.1038/s41588-018-0223-8 30275530 13. Darling AC Mau B Blattner FR Perna NT Mauve: multiple alignment of conserved genomic sequence with rearrangements Genome Res. 2004 14 1394 1403 10.1101/gr.2289704 15231754 14. Dewey, C. N. Aligning Multiple Whole Genomes with Mercator and MAVID. 221–235 (Humana Press, Totowa, NJ, 2008). 15. Paten B Herrero J Beal K Fitzgerald S Birney E Enredo and pecan: genome-wide mammalian consistency-based multiple alignment with paralogs Genome Res. 2008 18 1814 1828 10.1101/gr.076554.108 18849524 16. Darling AE Mau B Perna NT Progressivemauve: multiple genome alignment with gene gain, loss and rearrangement PloS ONE 2010 5 e11147 10.1371/journal.pone.0011147 20593022 17. Minkin, I., Pham, H., Starostina, E., Vyahhi, N. & Pham, S. C-sibelia: an easy-to-use and highly accurate tool for bacterial genome comparison. F1000Researchhttps://f1000research.com/articles/2-258 (2013). 18. Myers, G. & Miller, W. Chaining multiple-alignment fragments in sub-quadratic time. in Proceedings of the Sixth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’95, 38–47 (Society for Industrial and Applied Mathematics, USA, 1995). 19. Abouelhoda MI Ohlebusch E Chaining algorithms for multiple genome comparison J. Discret. Algorithms 2005 3 321 341 10.1016/j.jda.2004.08.011 20. Ohlebusch, E. & Abouelhoda, M. I. Chaining Algorithms and Applications in Comparative Genomics. (Handbook of Computational Molecular Biology, 2006). 21. Raphael B Zhi D Tang H Pevzner P A novel method for multiple alignment of sequences with repeated and shuffled elements Genome Res. 2004 14 2336 2346 10.1101/gr.2657504 15520295 22. Pham S Pevzner P Drimm-synteny: decomposing genomes into evolutionary conserved segments Bioinformatics 2010 26 2509 2516 10.1093/bioinformatics/btq465 20736338 23. Minkin, I., Patel, A., Kolmogorov, M., Vyahhi, N. & Pham, S. Sibelia: A scalable and comprehensive synteny block generation tool for closely related microbial genomes. in (eds Darling, A. & Stoye, J.) Algorithms in Bioinformatics. 215–229 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2013). 24. Marcus S Lee H Schatz MC Splitmem: a graphical algorithm for pan-genome analysis with suffix skips Bioinformatics 2014 30 3476 3483 10.1093/bioinformatics/btu756 25398610 25. Chikhi R Limasset A Medvedev P Compacting de bruijn graphs from sequencing data quickly and in low memory Bioinformatics 2016 32 i201 i208 10.1093/bioinformatics/btw279 27307618 26. Baier U Beller T Ohlebusch E Graphical pan-genome analysis with compressed suffix trees and the burrows-wheeler transform Bioinformatics 2016 32 497 504 10.1093/bioinformatics/btv603 26504144 27. Minkin I Pham S Medvedev P Twopaco: an efficient algorithm to build the compacted de bruijn graph from many complete genomes Bioinformatics 2017 33 4024 4032 27659452 28. Cleary, A., Kahanda, I., Mumey, B., Mudge, J. & Ramaraj, T. Exploring frequented regions in pan-genomic graphs. in Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics, 89–97 (Association for Computing Machinery, New York, NY, USA, 2017). 29. Vaser R Sović I Nagarajan N Šikić M Fast and accurate de novo genome assembly from long uncorrected reads Genome Res. 2017 27 737 746 10.1101/gr.214270.116 28100585 30. Sayers EW GenBank Nucleic Acids Res. 2019 48 D84 D86 31. Brudno M Lagan and multi-lagan: efficient tools for large-scale multiple alignment of genomic dna Genome Res. 2003 13 721 731 10.1101/gr.926603 12654723 32. Perry, E. Personal communication (2018). 33. Tajima F Statistical method for testing the neutral mutation hypothesis by dna polymorphism Genetics 1989 123 585 595 2513255 34. Armstrong, J. et al. Progressive alignment with cactus: a multiple-genome aligner for the thousand-genome era. Preprint at https://www.biorxiv.org/content/early/2019/10/15/730531 (2019). 35. Paten B Cactus graphs for genome comparisons J. Comput. Biol. 2011 18 469 481 10.1089/cmb.2010.0252 21385048 36. Fiddes IT Comparative annotation toolkit (cat)-simultaneous clade and personal genome annotation Genome Res. 2018 28 1029 1038 10.1101/gr.233460.117 29884752 37. Schwartz AS Pachter L Multiple alignment by sequence annealing Bioinformatics 2007 23 e24 e29 10.1093/bioinformatics/btl311 17237099 38. Sakharkar MK Perumal BS Sakharkar KR Kangueane P An analysis on gene architecture in human and mouse genomes Silico Biol. 2005 5 347 365 39. Pevzner P Tesler G Human and mouse genomic sequences reveal extensive breakpoint reuse in mammalian evolution Proc. Natl Acad. Sci. USA 2003 100 7672 7677 10.1073/pnas.1330369100 12810957 40. Kim J Reconstruction and evolutionary history of eutherian chromosomes Proc. Natl Acad. Sci. USA 2017 114 E5379 E5388 10.1073/pnas.1702012114 28630326 41. Luo H Phylogenetic analysis of genome rearrangements among five mammalian orders Mol. Phylogenet. Evolut. 2012 65 871 882 10.1016/j.ympev.2012.08.008 42. Kolmogorov M Chromosome assembly of large and complex genomes using multiple references Genome Res. 2018 28 1720 1732 10.1101/gr.236273.118 30341161 43. Kim J Reference-assisted chromosome assembly Proc. Natl Acad. Sci. USA 2013 110 1785 1790 10.1073/pnas.1220349110 23307812 44. Kolmogorov M Raney B Paten B Pham S Ragout—a reference-assisted assembly tool for bacterial genomes Bioinformatics 2014 30 i302 i309 10.1093/bioinformatics/btu280 24931998 45. Chen K-T Multi-car: a tool of contig scaffolding using multiple references BMC Bioinform. 2016 17 469 10.1186/s12859-016-1328-7 46. Aganezov, S. & Alekseyev, M. A. Multi-genome scaffold co-assembly based on the analysis of gene orders and genomic repeats. in (eds Bourgeois, A., Skums, P., Wan, X. & Zelikovsky, A.) Bioinformatics Research and Applications. 237–249 (Springer International Publishing, Cham, 2016). 47. Proost S i-adhore 3.0—fast and sensitive detection of genomic homology in extremely large data sets Nucleic Acids Res. 2011 40 e11 e11 10.1093/nar/gkr955 22102584 48. Portwood JL Maizegdb 2018: the maize multi-genome genetics and genomics database Nucleic Acids Res. 2018 47 D1146 D1154 10.1093/nar/gky1046 49. Onodera, T., Sadakane, K. & Shibuya, T. Detecting superbubbles in assembly graphs. in (eds Darling, A. & Stoye, J.) Algorithms in Bioinformatics. 338–348 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2013). 50. Sung W Sadakane K Shibuya T Belorkar A Pyrogova I An \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{O}}(m\mathrm{log}\,m)$$\end{document} O ( m log m ) -time algorithm for detecting superbubbles IEEE/ACM Trans. Comput. Biol. Bioinform. 2015 12 770 777 10.1109/TCBB.2014.2385696 26357315 51. Iliopoulos, C. S., Kundu, R., Mohamed, M. & Vayani, F. Popping superbubbles and discovering clumps: Recent developments in biological sequence analysis. in (eds Kaykobad, M. & Petreschi, R.) WALCOM: Algorithms and Computation. 3–14 (Springer International Publishing, Cham, 2016). 52. Brankovic L Linear-time superbubble identification algorithm for genome assembly Theor. Comput. Sci. 2016 609 374 383 10.1016/j.tcs.2015.10.021 53. Paten, B., Novak, A. M., Garrison, E. & Hickey, G. Superbubbles, ultrabubbles and cacti. in (ed Sahinalp, S. C.) Research in Computational Molecular Biology. 173–189 (Springer International Publishing, Cham, 2017).