Aligning the unalignable: bacteriophage whole genome alignments

Background In recent years, many studies focused on the description and comparison of large sets of related bacteriophage genomes. Due to the peculiar mosaic structure of these genomes, few informative approaches for comparing whole genomes exist: dot plots diagrams give a mostly qualitative assessment of the similarity/dissimilarity between two or more genomes, and clustering techniques are used to classify genomes. Multiple alignments are conspicuously absent from this scene. Indeed, whole genome aligners interpret lack of similarity between sequences as an indication of rearrangements, insertions, or losses. This behavior makes them ill-prepared to align bacteriophage genomes, where even closely related strains can accomplish the same biological function with highly dissimilar sequences. Results In this paper, we propose a multiple alignment strategy that exploits functional collinearity shared by related strains of bacteriophages, and uses partial orders to capture mosaicism of sets of genomes. As classical alignments do, the computed alignments can be used to predict that genes have the same biological function, even in the absence of detectable similarity. The Alpha aligner implements these ideas in visual interactive displays, and is used to compute several examples of alignments of Staphylococcus aureus and Mycobacterium bacteriophages, involving up to 29 genomes. Using these datasets, we prove that Alpha alignments are at least as good as those computed by standard aligners. Comparison with the progressiveMauve aligner – which implements a partial order strategy, but whose alignments are linearized – shows a greatly improved interactive graphic display, while avoiding misalignments. Conclusions Multiple alignments of whole bacteriophage genomes work, and will become an important conceptual and visual tool in comparative genomics of sets of related strains. A python implementation of Alpha, along with installation instructions for Ubuntu and OSX, is available on bitbucket (https://bitbucket.org/thekswenson/alpha). Electronic supplementary material The online version of this article (doi:10.1186/s12859-015-0869-5) contains supplementary material, which is available to authorized users.


Background
The most abundant, and probably the most diverse, biological entities are bacteriophages, the viruses that infect bacteria. Helped by recent advances in sequencing, comparative studies [1][2][3] use dozens -even hundreds -of genomes from bacteriophages that infect single or related bacteria species. Using dot plots and clustering techniques, these studies produce meaningful clusters that share significant similarity (> 50 %). In order to explore the relations between bacteriophages in the same cluster, it seems natural to turn to whole genome multiple alignments. One of the main features of whole genome aligners is that they take into account genome rearrangement events that scramble the order of large segments of chromosomes (see [4] for a review). All of these approaches are based on the principle that sequence similarity is a good predictor of functional similarity. And it is, most of the time. One notable exception are bacteriophage genomes, in which similar biological function may be encoded by dissimilar sequences -sequences with no detectable similarity, either as nucleotide sequences, or as amino acid sequences -encoding different protein folds, rendering traditional multiple sequence aligners mostly useless [5][6][7].
In this paper, our goal is to construct biologically meaningful multiple alignments of whole bacteriophage genomes from the Siphoviridae, the largest family of tailed bacteriophages [8]. In order to achieve this, we exploit unique structural properties of these genomes. The main one is that even loosely related tailed bacteriophages are often functionally collinear, meaning that different functionalities mostly follow the same order on the genomic sequence, which is, up to a circular permutation: lysogeny, DNA assembly, head morphology and DNA packaging, tail assembly, and lysis [9]. The second one is that genome size is constrained by the fact that it must fit into a capsid whose shape, thus volume, is geometrically determined by a handful of genes. This size constraint implies that segment duplication is a rare event within a genome, and when it occurs, duplicated sequences are short. Last, but not least, bacteriophage genomes are characterized by an "unusually high degree of horizontal genetic exchange in their evolution" [10], resulting in large sequences -up to thousands of base pairs -that are exact or almost exact matches between different strains.
As an example of this last feature, Fig. 1 compares two segments of S. aureus bacteriophages 88 and 92 (see Table 1 for all accession numbers of bacteriophages discussed in the paper). The figure is composed of 5 columns, three of which corresponding to exact matches of length 259, 35 and 49, and the other two columns containing large segments occupying the same position but without any 'detectable' similarity, meaning neither the nucleotide nor the translated sequences produce any significant BLAST hit.
The mosaic patterns exhibited by the comparison of two bacteriophages, as illustrated in Fig. 1, have given rise to the modular theory of phage genome organization [5], which postulates that biological functions are grouped into modules whose order is mostly conserved along the genomic sequence. Each module has variants that perform the same function, possibly encoded by dissimilar sequences. The fact that two variants are aligned in the pairwise comparison may allow the transfer of functional annotation between two collinear phages.
While it is easy to do a pairwise comparison of genomes, upgrading the comparison to multiple genomes is not simple. The standard approaches used by multiple sequence aligners often start by identifying anchors, that are similar segments of significant length shared by all genomes, and then align the sequences between the anchors. Anchors exist in sets of related bacteriophages, but may be very short or very far apart, and the sequences between the anchors may fail to align properly.
The theory behind our framework is based on partial order alignment graphs [4,[11][12][13] which were initially developed for standard multiple sequence alignments. Most applications of these graphs require, in the last phase of the alignment, a linearization of the graph. However, an interesting suggestion appearing in [13] is to skip this last step and work directly with the partial order: this is exactly what is needed for bacteriophage genomes, but apparently the approach has yet to be applied in this context.
Among the multiple sequence aligners, progressive-Mauve [14] is one of the few that recognizes the need to identify 'local' anchors, that are shared by a subset of the target genomes. Unfortunately, Mauve alignments are linearized, blurring the combinatorial properties of the partial order.
Our evolutionary model includes typical mutations characterizing sequence evolution such as substitutions and indels, gains and losses of functions, and recombinations. For sets of genomes that have evolved under these conditions, we will show that it is possible to identify functionally related sequences even when they lack similarity. We also detect large rearrangements events that contradict the functional collinearity hypothesis, such as gene transpositions, and duplications.
In this paper, we report the conception and implementation of Alpha (Alignments of bacteriophage genomes), the first aligner specifically designed for whole bacteriophage genomes. With the help of partial order structures, Alpha captures the unique mosaic structure of bacteriophage genomes, and provides an interactive graphical interface with the generated multiple alignment. We also give a detailed comparison between Alpha and progres-siveMauve [14] alignments.

Partial order alignment graphs and functional collinearity
This section draws heavily on multiple Whole Genome Alignment (WGA) tools and definitions. However, since we have the goal of identifying functional analogs, some classical notions of this field will have a somewhat different meaning: in these cases, we try to be precise and to underline the differences.
A match M between two genome sequences G and H will be denoted M = .v] }, where s and u are the start positions of the match in the respective genomes, and t − s = v − u; it asserts the equalities We next formalize the notion of homologous positions, defining what constitutes a 'column' of a multiple alignment, along the lines of J. Kececioglu's original paper [11].
Let M be a set of matches on genomes G = {G 1 , G 2 , . . . , G n } and let p ij denote position j of genome G i . Two positions p ij and p kl are equivalent if there is a match in M that asserts that through a sequence of such assertions. The equivalence class of p ij is denoted p ij , and we will refer to it as a column. The support of a column sp p ij is the set of genomes that have a position in the column p ij .
Note that there is no requirement that a column span the complete set G, and it usually does not in bacteriophage genomes. The notion of support generalizes the notion of anchor that is used in WGA, which corresponds to columns whose support is equal to G.
The column graph is obtained by linking the different columns according to their order in each genome. Formally, it is the directed (multi)graph whose vertices are the columns p ij , and there are edges from p ij to p kl for all genomes G g that have positions p gs in p ij and p g(s+1) in p kl . An example of the construction is given in Fig. 2, part (B).

Definition 1. Let M be a set of matches between
In most WGA systems, collinearity is meant to represent sequences that 'align well' . Classical collinearity is broken by pairs of dissimilar sequences, and is an indication of gene rearrangements, gains or losses, in most higher organisms. Functional collinearity, on the other hand, excludes rearrangements, but allows dissimilar sequences to be compared. Interestingly, both concepts have the same formal definition, but the use of the column graph is quite different from one application to the other.
A necessary condition for a column graph to be functionally collinear is that each column have at most one The column graph and the alignment graph. a A set of two exact matches on three genomes; the total length of matches is 8. b The matches define columns for which all positions contain the same nucleotide. The column graph is obtained by adding edges corresponding to consecutive positions in the genomes. c Consecutive columns with the same support are merged to form the alignment graph, whose vertices are exact alignments position in a given genome, since the presence of two or more positions would generate a cycle. We have therefore the following definition:

Definition 2. A set M of matches is said to be valid if, for each column p ij , both p ij and sp p ij have the same number of elements.
The most likely sources of invalid sets of matches are duplicated segments, and tandem repeats with a variable number of repeats, which may cause overlapping matches. Most WGA deal with these invalid matches by eliminating matches that are not unique in a genome, and by trimming overlapping matches that result from tandem repeats. Given the simple data structure that we use for the alignment graph -equivalence classes -we eliminate duplicated segments, and trim overlapping matches by post-processing the columns: any column that has more than one position for a genome is split into singletons.
Allowing the length of matches to be very small may cause misalignments, and may also create cycles in the column graph, due to random small transpositions. It is therefore wise to have a minimal length m for matches. In the Results section, we will discuss how the value of m is set, and how it can be changed.
When the column graph has no cycles, we use a condensed version, called the alignment graph, in which consecutive vertices with the same support are merged, as in Fig. 2, part (C).

Construction of the alignment graph
Given a set of genomes G = {G 1 , G 2 , . . . , G n } we use GenomeTools [15], a suffix-array based program, to get a set of exact matches of size at least m. The parameter m is normally chosen as the minimum length such that there is no cycle in the graph. It is set by default at 15, but can be changed to any desired value. Large values of m are sometimes useful in initial explorations of a set of bacteriophages and small values can be used to refine alignments.
The main steps of the algorithm underlying Alpha are described in Algorithm 1. It takes as input a set of matches of minimum length m, constructs the column graph, and the alignment graph when possible. After initializing the graph, it uses the Union-Find algorithm [16] to construct the columns. This step runs in O(α(M )M ), where M is the total length of matches, and α is the inverse of the Ackermann function [17]. The third step visits all classes, determining whether they are valid. For each class, this determination can be made in O(n) steps, where n is the number of genomes, yielding a complexity of O(nN ) for splitting all the nonvalid classes, where N is the total length of genomes. Note that at most N new classes can be created by the splitting process.
In the fourth step, edges between consecutive columns are added to the graph, and it is checked for the presence of cycles that would indicate non-collinearity. This can be done in time proportional to N . Assuming the genomes are collinear, consecutive columns with the same support are then merged to give the alignment graph.
The total complexity is thus O(α(M )M + nN ). We will see in the Results section that the quantities M and nN vary greatly depending on the data.

From the alignment graph to gapless alignments
The vertices of an alignment graph are exact alignments, since all rows are equal. Regions that do not contain exact matches of length at least m remain unaligned. In order to increase the compactness of whole genome alignments, we contract the vertices further with the following requirement: A contracted vertex is obtained by merging into a single vertex all vertices between a pair of contractible vertices U and V, including U and V. Contracted vertices are gapless alignments, since all sequences in them have the same length. Without upper bounds to the length of contracted vertices, there is a small chance that some of these alignments have no biological foundation, since insertions and deletions could conspire in producing sequences with the same length, but low similarity. As we will see in the Results section, bacteriophage genomes are well behaved in this regard. Figure 3 shows an example of a contracted vertex together with its expanded graph. All vertices display the coordinates of aligned segments, except for narrow ones, together with the common length of these segments. In addition, a contracted vertex also displays the percentage of identity, which is the percentage of columns that have the same conserved nucleotide in each genome of its support. Since the contracted version of the alignment graph turned out to be the most useful in practice (see the Results section), we will refer to it as the alignment graph, and use the term 'expanded alignment graph' when the full version is required.

Results and discussion
In this section, we first present examples of multiple alignments and what can be deduced from them. We next assess the validity of the gapless alignments inferred by Alpha by submitting the aligned sequences to three standard aligners. Finally, we compare the results of Alpha and progressiveMauve [14] on whole genome alignments.
Alpha is an interactive tool that allows the manipulation and the visualization of whole genome alignment graphs with hundreds of vertices, involving dozens of species: we can only hope to give a glimpse of the possibilities. Alpha's graph layouts are powered by the open source graph visualization software Graphviz [18].

The alignment graph
The Alpha aligner, in its most basic mode, takes as input a file containing genomes in Fasta format, and produces an alignment graph when the genomes have no large rearranged segments. The vertices of the alignment graph are gapless multiple alignments inferred by Alpha. The vertices contain the positions of the aligned segments in each genome. They also display the length of the alignment and the percentage of columns that have the same nucleotide in each genome -for readability, vertices smaller than 20 bp may be masked. Alignments are connected by colorcoded arrows, one color for each genome. Figure 4 shows an example of a local alignment of four S. aureus bacteriophages, in which two deletions are easily identified. Dotted arrows replace vertices that span less than 20 bp, implying that phage B236 lacks a group of functionalities that spans over a thousand bp, and phage phiETA3, one that spans over 500 bp.
Bacteriophage genomes have been sequenced for decades, and since they can adopt different linear and circular configurations in their complex life cycle, there is no universal consensus as to where the sequence begins. Functional collinearity is biologically defined up to a circular permutation of the sequences, but bacteriophage Fig. 4 A local alignment of four S. aureus bacteriophages. In column 3, three phages are in a 587 bp gapless alignment, with 99 % identical columns. A major deletion in phage B236 spans columns 2, 3, and 4, and the corresponding arrow is dotted to reflect the fact that some basepairs are not shown, 6 bp in this case sequences in databases are linear. If a set of bacteriophage genomes was obtained from different projects, or different laboratories, it may be necessary to synchronize them using a simple procedure that looks for the largest similar sequence shared by all genomes, and sets the beginning there. We call this process normalization, and Alpha checks its input to decide whether normalization is needed. Figure 5 shows a group of S. aureus bacteriophages whose sequences have been normalized: positions of bacteriophage phiNM1 are offset by more than 16 500 bp with respect to the other genomes, while maintaining collinearity with them. In this figure the second and fourth columns identify two pairs of large sequences. These pairs share similar loci within the alignment but they lack detectable similarity, and they are predicted to be variants of the same module. However, the second column splits sequences {A, D} from {B, C}, while the fourth splits sequences {A, B} from {C, D}, illustrating the very peculiar behavior of horizontal transfer in bacteriophages.

The anchor view
Anchors are alignments that span the whole set of n genomes under study. They are maximal exact matches, in the sense of Hohl et al. [19], but they can be as short as 1 base pair. Each anchor is constructed as the intersection of at least n − 1 exact pairwise alignments of length at least m, thus the reliability of anchors increases with both the number of genomes, and the parameter m. The ordered set of anchors forms the backbone of a set of bacteriophage genomes and captures their common core, as defined in Mosaic [20].
The anchor view of Alpha presents a sequence of anchors for the whole genomes, or for selected regions: it is the normal starting point to explore a set of genomes. Figure 6 shows an example of an anchor view for a set of 29 mycobacteriophages computed with m = 175. There are not many anchors for such a large value of m, but they are very well supported because each anchor is defined by at least 29 exact matches whose length is at least 175 bp. Since anchors are articulation points of the alignment graphs -removing them disconnects the graph -each group of sequences spanning from one anchor to the next can be explored separately, as we show in Fig. 7.
Once a pair of anchors is selected in the anchor view, an alignment spanning from one to the other is generated. It is possible to ask the aligner to align the sequences between them, and a new value of m is automatically computed, generally much smaller -or set to any desired value. In Fig. 7, for example, the value of m is 15.

Assessing the validity of alpha gapless alignments
In order to assess the gapless alignments produced by Alpha, we ran it on three datasets. Table 2 Table 1 for details).
Alpha computed 491 gapless alignments of 2 or more sequences on these datasets. All 491 alignments were re-aligned using Clustal Omega [21], T-Coffee [22] and Muscle [23], using the default values for DNA alignments, and we report the number of alignments that include gaps. For each dataset, Table 3 gives the number of vertices in the expanded and contracted alignment graphs, the number of vertices containing at least two sequences, the number of contracted vertices, and the number of gapped alignments obtained by the three aligners.
All four aligners agree on 479 of the 491 = 78 + 154 + 259 gapless alignments proposed by Alpha. The remaining 12 were contested by one or more of the aligners, for a total of 18 gapped alignments. They cover 8 different regions of the genomes since, in four cases, aligners proposed gapped alignments for the same region in both the Myco6 and Myco29 datasets.   . 6 The anchor view. The alignment of 29 Mycobacterium bacteriophage genomes is a challenge to visual approaches. Here we use a large value of m = 175 to display a sequence of anchors, or backbone, in a single screen. The user can then zoom quickly to a specific region, and detect, as Fig. 7 shows, interesting subsets of genomes. Boxes group pairs of anchors bounding a gapless alignment Judging whether nucleotide alignments are in fact gapless is a delicate task. When proteins sequences were available, we used the gapless protein alignments to reject the corresponding gapped nucleotide alignment. This was done by using tblastx with sequences that do not share the same gap patterns -all the gapped alignments showed only two different gap patterns -and by confirming annotations using blastx. Using this method, we could rule out 12 of the 18 alternative gapped alignments, as Table 4 shows, for Regions 1, 2, 3, 4, and 8. All sequences and alignments are available in the Additional file 1.
For Region 5, the three aligners proposed three different gapped alignments, when aligning the 6 genomes of dataset Myco6, but they all switched to gapless alignments when aligning the same region in the 29 genomes of dataset Myco29, which contains Myco6. Using the principle that alignments with more sequences should be more accurate, those three gapped alignments were also ruled out.
Of the remaining three gapped alignments, one was proposed by Muscle, while Alpha, Clustal Omega and T-Coffee proposed gapless alignments (Region 7). Given the overall poor performance of Muscle, we chose the rule of the majority. Finally, for Region 6, there was a tie between the four aligners, and there was no decision.
As the results of Table 4 clearly indicate, Alpha gapless alignments were confirmed in all cases except for the tie in Region 6, where Alpha and Clustal Omega -the two best gapless aligners -predict a gapless alignment against T-Coffee and Muscle, who have the most confirmed misalignments. This is a welcomed and surprising result, since the only parameter of Alpha is m, the minimal length of exact matches: it does not rely on thresholds, and does not maximize any score. Instead, it relies on the transitive properties of the equality relation to provide reliable anchors, and on the unique constraint imposed on bacteriophages who, before traveling, must pack their whole genome in a small suitcase.

Comparison with Mauve alignments
In the preceding section, we showed that Alpha alignments are sound, in the sense that they predict biological meaningful similarities. The next question is to evaluate to what extent Alpha captures all meaningful similarities.
In order to do this, we compared Alpha with progres-siveMauve [14], since it is one of the few aligners that explicitly computes the partial order underlying collinear blocks, storing this information in the .backbone file generated during an alignment. As Fig. 8 shows, Mauve alignments are displayed in a linear way, using colors to show which subsets of segments are in an alignment. With  Fig. 6. Working with a large number of phages allows us to quickly identify subsets of phages that have distinct patterns, such as phages A, C and E (U2, Alvin and Solon) in this example more than a handful of genomes, such as the Myco29 dataset, this type of visualization quickly becomes impractical (see Additional file 2: Figure S1 that shows the Mauve alignment of the region spanned in Fig. 7).
Comparing partial order alignments stemming from two different aligners turned out to be a daunting task: an alignment proposed by one aligner can be broken into several alignments by the other, since alignments may involve  All blocks predicted by both aligners were considered as valid alignments. Furthermore, short blocks (less than 100 bp) that were predicted by only one aligner -mostly by Mauve -were considered as valid: they averaged 24 Fig. 8 Visual comparison. Alpha and Mauve alignments of the same region of the Staph4 dataset. Numbers in red were added to show the correspondence between Alpha alignments and Mauve color blocks. All Alpha alignments are found in the Mauve alignment, but Mauve apparently detects an alignment of the four genomes that is missed by Alpha: the block labeled 'x'. We show in the text that block x is a misalignment bp. Three of the remaining large blocks -two predicted only by Mauve, and one only by Alpha -were validated using protein alignments. Figure 9 shows, in the various green-blue-shaded rectangles, the repartition of valid alignments, measured in total bp length, between common and exclusive alignments for Mauve and Alpha. These results are further partitioned by discriminating whether a block is an alignment, or an extension of an alignment.
Seven alignments predicted only by Mauve, totaling 2803 bp, were suspicious. They correspond to the yellow and red rectangles of Fig. 9. Three of them had non-significant tblastx results, despite the fact that at least one of the two sequences was annotated as coding. The remaining four alignments had a high number of gaps, contradicting the amino acid alignments of the corresponding homolog coding sequences. An example of each of these cases is detailed below, and the evidence for all others is detailed in Additional file 3.
When one excludes Mauve misalignments -that is, considering the green and blue rectangles of Fig. 9 -both aligners generally agree on the topology of the alignment graphs. As could be predicted, due to Mauve's seed and extend strategy, it has wider alignments compared to Alpha, and this accounts for the 778 bp extensions that Alpha could not predict. On the other hand, the 8 bp extensions that Alpha alone predicted, in four different blocks, were in fact verified to be exact matches that Mauve somehow missed. The two aligners also failed to detect a few small alignments, with a little advantage to Alpha: Mauve detected 294 bp that Alpha did not, while Alpha detected 358 bp that Mauve did not.
The high rate of Mauve misalignments on bacteriophage genomes is problematic. More than 20 % of Mauve column alignments on this dataset are either random, or belong to nucleotide alignments whose gaps contradicts the corresponding annotated proteins alignments. It should be noted, however, that the Mauve algorithm uses Muscle to produce alignments between its own anchors, and the way Muscle introduces gaps in nucleotide alignments does not seem to be appropriate for bacteriophage genomes, as we already saw in the preceding section on validating Alpha gapless alignments.
Overall, Alpha's alignment strategy captures the essential features of bacteriophage genomes. Regions that are similar are detected, while more dissimilar regions are not aggressively aligned. Being a conservative aligner, Alpha relies on functional collinearity to predict distant

Conclusion
Due to their peculiar mosaic structure where similar functions do not correspond to similar sequences, bacteriophage genomes are not well-suited to traditional whole genome alignment techniques. On the other hand, they exhibit features that can be leveraged to obtain alignments, most notably functional collinearity, a low duplication rate, and the presence of long shared sequences.
In this paper, we presented a mathematical model based on partial order graphs for performing multiple alignment of bacteriophage whole genomes, along with algorithms to operate on the model. Relying exclusively on the equality relation, the model is almost parameter free, greatly reducing the need to calibrate the aligner, yet delivers biologically meaningful results. The model has been implemented in the form of an interactive aligner that can perform multiple alignments of dozens of genomes and present the result in an attractive format.
We also showed that Alpha, used on bacteriophage genomes, produces biologically meaningful alignments, while avoiding the high rate of misalignments of complex heuristics such as progressiveMauve.
Our model supposes that all genomes under consideration are functionally collinear. This is often the case, but not always. Our program can detect when this condition is not satisfied -indicating the presence of rearrangements -but does not perform the alignment in such a case. A short-term goal is to extend our mathematical model and aligner to deal with rearrangements.
Some bacteriophage genomes present in the online databases are well annotated, other are less so. Another goal is to extend the aligner in order to perform automated transfer of annotations using the generated alignments.
Finally, while we focussed this study on Siphoviridae, we plan to test Alpha on more general viruse families in which horizontal transfer is widespread, and for which the collinearity property may hold.