Proceedings

Biologists frequently search databases for similar sequences to newly sequenced genes or proteins ([Lesk02], 177). The results, however, are usually a set of close matches whose interrelationships may need further exploration. For example, it is often of interest to trace the evolutionary history of a series of mutation events that gave rise to different versions of a protein. This is the multiple sequence alignment problem ([Edga06]).

Two-Way Global Alignment

A multiple sequence alignment is usually built up from a series of two-way alignments, for which an optimal dynamic programming solution has existed since 1970 ([Need70]). This calculates a ‘global’ alignment, because it aligns every character in the two sequences.

Figure 4: Needleman-Wunsch Algorithm

In Figure 4 a 2-D matrix representing the alignment of the two strings is filled from the bottom right-hand corner to the top left-hand corner, one cell at a time. Each cell is filled with 1 if the letters at those coordinates match, otherwise 0, plus the maximum cell-value of the next row and column down and to the right. Once filled, the matrix can be used to recover a common subsequence with the highest score by following matching pairs from 0,0 down and to the right. There may be more than one such optimal path through the matrix.

Although it is theoretically possible to extend their method to n sequences, the time and space requirements increase exponentially for n > 2. For this reason only the 2-way alignments of all possible pairs of sequences are normally computed, which are then used to construct a close approximation to a global n-way alignment. Even this shortcut turns out to be quite expensive, since the number of possible pairs of n sequences is n(n-1)/2. And the computation of the alignment for each pair of sequences also costs O(MN(M+N)) time, where N and M are the lengths of the two sequences. To give some idea of just how costly that is, the computation of the optimal pairwise alignments of the 36 versions of the Sibylline Gospel ([Schm08]), i.e. 630 pairwise combinations, each around 2,500 characters long, took 6.6 hours on an ordinary desktop computer.

Progressive Alignment

One popular way of aligning multiple sequences is progressive alignment ([Feng87]). In this method the global pairwise alignments of a set of sequences are scored for similarity and the results are used to construct a phylogenetic tree, which records the stemmatic relationships between sequences, to act as a ‘guide’ for the merging process. The multiple-sequence alignment or MSA (Figure 5) is then built up by starting with the most similar pair of sequences, adding gaps where needed (and never removed once inserted), then gradually adding in the other sequences as dictated by the guide-tree. Modern examples of this approach include CLUSTALW ([Thom94]) and T-COFFEE ([Notr00]), although the latter augments the usual progressive technique with iterative refinement.

Figure 5: Multiple Sequence Alignment

Disadvantages of Progressive Alignment

Because this method was originally developed to handle short sequences, it doesn’t take account of transpositions, or repeated segments, which appear to be common features of longer protein sequences ([Pevz04]; [Delc02]; [Gu99]).

Also, it is a greedy algorithm sensitive to the order in which the sequences are added to the MSA. Alignment errors that occur early on get locked in and lead to greater divergences later ([Thom94]).

The initial pairwise alignments are particularly expensive to calculate. This limits the number and length of sequences that can be aligned.

Progressive alignment is based on the concept of edit distance ([Leve66]). This is a count of the number of insertions, deletions and substitutions of individual characters that need to be performed to turn one string into another. This can only ever approximate the real biological processes that gave rise to the differences ([Lass02], 13; [Lesk02], 168).

Partial Order Alignment

A relatively new technique of multiple sequence alignment called POA or Partial Order Alignment has been described by Lee et al. ([Lee02]) and Grasso et al. ([Gras04]). The roots of their method go back to Morgenstern’s Dialign program ([Morg96]). In that method, as in MAFFT ([Kato02]), highly similar regions within two sequences are first identified. They are then resolved into a non-overlapping set of diagonals in n Dimensions, then added to the MSA using a greedy strategy. Unaligned regions are then processed recursively. The difference here is that gaps are no longer explicitly introduced: they are simply what is left over once the aligned portions have been determined.

POA takes this strategy one step further: it replaces the MSA matrix with a directed graph. Instead of pairwise alignment they align each successive sequence – in no particular order – to a partly-built PO-MSA (Partial Order Multiple Sequence Alignment). In their method the initial computation of n(n-1)/2 pairwise alignments is not a preliminary step, but is instead integrated into the alignment process. Because a PO-MSA contains far less data than that contained in all possible pairs of sequences, the computational cost is drastically reduced.

Figure 6 shows their adaptation of the dynamic programming method to the directed graph structure. Instead of aligning two linear sequences within a flat 2-D matrix, in the the PO-MSA method the 2-D matrix becomes a complex planar surface corresponding to the partly merged structure of the directed graph, on the one hand, and the new linear sequence, on the other. At each join between surfaces the program has to decide which surface to align with. Regardless of the order in which the sequences are added to the alignment, each possible pairwise combination of sequences will be considered, since each new sequence is compared to all those already in the graph.

Figure 6: PO-MSA Alignment

The speed of the program is impressive, but the quality of the alignments is somewhat less than could be achieved by traditional progressive alignment techniques ([Lass02]). Grasso et al. ([Gras04]) think the cause is partly that the original algorithm still appears to be sensitive to the order in which sequences are added to the PO-MSA. This is probably because, unlike the progressive method, each new version is aligned only with those versions already in the graph, not with those still to be added.

Disadvantages of POA

In spite of its innovative data model, the POA graph and associated software can’t handle transpositions or repeated segments. The remaining order sensitivity seems to be an inherent deficiency of this approach.

The second area in bioinformatics that provides a close analogy to the current problem is a more radical solution to the multiple sequence alignment problem. It evolved from the related field of contig-assembly. Starting in the 1980s DNA sequencing machines were developed which output fragments of DNA, typically 650 base pairs in length ([Delc02], 2480). These short segments have to be assembled correctly into longer ‘contigs,’ which is an NP-hard problem. Pevzner ([Pevz89]), Idury and Waterman ([Idur95]) and Pevzner ([Pevz01]) developed a heuristic method in which each DNA fragment was broken up into a series of short overlapping substrings of fixed length. For example, if all possible substrings of length 3 are calculated for the string ‘AGTCGTATAC’, then the original string can be reassembled by overlapping these substrings to form the directed graph of Figure 7.

Figure 7: A Sequence Graph

Each node represents two characters shared by its two adjacent arcs, which are labelled with the 3-character substrings. The cycles record the repeats ‘GT’ and ‘TA’. In practice the length of the substrings would be much greater than 3 – around 20 or so ([Pevz01], 9751).

This enabled them to reduce the contig-asssembly problem to one of finding an ‘Eulerian’ path through the graph. An Eulerian path is one that visits each arc exactly once, although Idury and Waterman ([Idur95], 295) interpret this to mean a path that may cross the same arc multiple times as if it were a set of parallel arcs. The problem can then be solved in linear time by Fleury's algorithm ([Aldo00], 202).

The original de Bruijn graph represented all possible strings of length k over an alphabet A and was entirely cyclical ([Debr46]; [Good46]). An example is shown in Figure 8, which shows all the 3-character sequences of the alphabet {0,1}.

Figure 8: A De Bruijn or Good Graph

Raphael et al.'s approach is to add insertions and deletions to the ‘de Bruijn’ graph, calling the result an A-Bruijn graph, and to use it as their data representation ([Raph04]). Although they succeed in representing transpositions (or ‘shuffles’) and repeats, their method only aligns the highly similar regions of protein sequences, while the rest is processed separately using progressive alignment techniques.

The third area of bioinformatics that overlaps with the current problem has developed only recently. As a result of the successful sequencing of entire genomes, biologists became interested in much larger scale alignments than had been possible via traditional multiple sequence alignment. The limitation is imposed by the dynamic programming technique and its essentially O(N²) performance and space requirements per pair of sequences. Also, at the level of millions of base pairs, relocation, reversal and duplication events have become impossible to ignore.

This need has been met by the rapid development of fragment based aligners, using an anchoring strategy, similar to Heckel's ([Heck78]) algorithm. Anchoring is the process of first identifying the most similar parts of two or more sequences. Once those have been identified and merged, the program searches for further alignments in the unmerged portions of the sequences. The early fragment based aligners worked only with two versions, e.g. MUMmer ([Delc02]); MultiPipMaker, based on BLASTZ ([Schw00]); Avid ([Bray03]) and LAGAN ([Brud03a]), and do not handle transpositions, reversals and duplications within the two sequences being compared.

Mauve

Mauve ([Darl04]) handles transpositions, reversals and duplications, and multiple long sequences. Like the other fragment-based alignment programs, Mauve is based on anchors. The alignment process is shown in Figure 9.

Figure 9: Anchor-based Alignment in Mauve

In the first stage the program searches for Maximal Unique Matches across all the versions, the so-called ‘multiMUMs’. The multiMUMs are used to generate a phylogenetic guide tree. Successive multiMUMs that don’t overlap are then merged into local co-linear blocks or LCBs (2). These act as anchors and allow the identification of further multi-MUMs within and between the LCBs. The identification and merging of LCBs then proceeds recursively (3). Finally, progressive alignment is performed on each LCB using the guide tree.

The authors don’t say how fast it is, but it must be slower than Shuffle-LAGAN since it deals with more versions and uses progressive alignment. Its weakest point appears to be that it has no special data representation for the transpositions, and instead uses colour and interconnecting lines to adapt the traditional MSA matrix display.

The practical solutions in this section are not much use for the current problem since they only deal with two versions at a time. Also, since the current problem is NP-complete, it will not be possible to design an algorithm to calculate an optimal solution with respect to edit distance. This is not as bad a result as it seems, since:

A well-designed heuristic algorithm, then, might actually produce better results, as well as take less time.

Froger ([Frog68], 230) describes the execution of what was probably the world's first collation program, running on a Bull GE-55 computer that had only 5K of internal memory. It compared two texts, one word at a time. If there was a mismatch the program gradually increased its lookahead to 2, 5, 25, ... words, comparing all of the words with one another in this window until it found a match, then sliding the window on by one word at a time. This allowed it to recognise insertions, deletions and substitutions over short distances.

This sliding window design, although later much enhanced, appears to have been the basis of all collating programs until recently. Its main disadvantages are that it can’t recognise transpositions, and that it will fail to correctly align inserted or deleted sections larger than its window size. Computational complexity for each pair of versions is O(NW²), where W is the sliding window size, and N is the average length of the two sequences. Typically a modern collation program will compare each version against a single, and possibly idealistic, base text. Examples include COLLATE ([Robi89], 102) and TUSTEP ([Ott00]).

Spencer and Howe ([Spen04]) apply biological progressive alignment to medieval texts of the Canterbury Tales and Parzival. They point out that current automatic collation techniques mimic the original manual process too closely. In particular they object that comparing each version against a base text omits important comparisons between individual witnesses. However, their main objective is to generate a phylogenetic tree or stemma, rather than to produce an apparatus.

JUXTA ([Mcga06]) is an open source program, whose collation algorithm is an adaptation of Heckel's (1978). The adaptation consists in reducing the granularity from a line to a word, and in calling it recursively on the unaligned parts. JUXTA collates plain text with embedded empty XML elements as milestone references, like the COCOA references in COLLATE. JUXTA can also be used to produce an apparatus criticus.

MEDITE ([Bour07]) is a two-way alignment algorithm based loosely on the fragment based aligners Shuffle-LAGAN and Mauve. The advantage of this approach, in which the transpositions and matches are identified first, then substitutions, insertions and deletions, is that a linear performance can be achieved. The algorithm is basically:

Collation algorithms have several potential weaknesses. They are mostly slow, and generally collate plain text, rather than XML. Nowadays transcriptions of documents in linguistics and the humanities are more likely to be in some form of XML. The problem is that if one collates XML, then the generated apparatus will contain unmatched tags. Markup will then have to be ignored, but that means losing valuable information.

The inclusion of COCOA or empty XML references embedded in the text precludes the use of any other kind of markup, such as structural markup. But references are only needed for attaching the apparatus to the main text, and this could be done away with in a purely digital edition.

Comparison needs to be between all versions at once. Currently only Spencer and Howe's ([Spen04]) method achieves this.

Transpositions need to be detected, especially when a work is witnessed by many independent copies. Few of the algorithms do this, and none over n versions.

Taken together, the most advanced algorithm appears to be Bourdaillet's ([Bour07]). However, it can only consider two versions at a time, and so cannot be directly applied to the current problem.

1. Find the longest MUM (Maximal Unique Match) between the new version and the versions of the variant 
graph.
2. Merge the new version and the graph at that point, partitioning them both into two halves.
3. Call the algorithm recursively on the two subgraphs, each spanned by one segment of the new 
version.

When comparing two strings, a Maximal Exact Match, or MEM, is a sequence of consecutive matching characters in both strings that cannot be extended at either end. A MUM, or Maximal Unique Match, is a MEM that occurs only once in the two strings ([Delc02]). In Figure 11 ‘the quick brown fox’ and ‘ brown fox jumps over the ’ are both MUMs, while ‘the ’ is a MEM that also occurs elsewhere.

Figure 11: Maximal Exact Match and Maximal Unique Match

In the current problem MUMs will be sought between the new version being added to the graph and the entire graph, the longest of which is the best MUM. They need to be MUMs and not MEMs because otherwise erroneous alignments would be generated. For example, choosing arbitrary instances of multiple spaces at the start of lines, or repeated sequences of markup tokens as points of equivalence would lead to serious misalignments.

Finding the best MUM is related to the problem of finding the Longest Common Substring (LCS) between two strings. A naive approach to finding the LCS might look like this:

Two problems with this naive approach are immediately clear: firstly it is inefficient, since it will take O(MN) time, where M and N are the lengths of the two strings. Secondly, the number of strings in the table will be enormous, as many as M²N². A much faster and more space-efficient way is to use suffix trees.

A suffix tree stores all the suffixes of a given string. In Figure 12 each path from the root to a leaf spells out one possible suffix of ‘Woolloomooloo’. In the leaves are stored the start position of each suffix in the original string. A node labelled λ indicates that at this point in the tree a suffix occurs with no further text, e.g. for the suffixes ‘oo’ and ‘o’. A suffix tree can be used to find any substring of length M in that string with just M comparisons, however long the string may be.

Figure 12: A Suffix Tree of ‘Woolloomooloo’

A suffix tree can be constructed in linear time using Ukkonen's algorithm ([Ukko95]). The LCS between two strings can also be calculated in linear time by building what is called a ‘generalised suffix tree’ ([Gusf97], 125). A generalised suffix tree represents n strings instead of 1. For two strings this method constructs a suffix tree from one of the strings in linear time, then carries on building the same tree with the second string. The leaf nodes are labelled by the number of the string and the starting index in either or both strings. Any internal node that has only two leaves as children, one from each string, or any leaf shared by the two strings, must be the suffix of a MUM. The MUM itself can be recovered by counting backwards in both strings until a mismatch.

Unfortunately, it is unclear if one can build a suffix tree from a variant graph. How are the suffixes to be identified, and how can one represent starting offsets within arcs that belong to multiple versions? In any case the gains from using a generalised suffix tree appear to be minimal in practice, since its use greatly increases storage requirements.

The method adopted here resembles the streaming technique used in MUMmer ([Delc02]). In this method only one suffix tree is built. The other sequence is then streamed against it to find matches.

Figure 13: Finding Matches

The place of the streamed sequence can be taken by the variant graph, whose characters can be exhaustively listed by breadth-first traversal ([Sedg88], 430f). Each successive character from the graph is traced in a suffix tree built from the text of the new version, as shown in Figure 13. At the end of an arc, for example, at node 3, MATCH calls itself recursively on each diverging arc until a mismatch occurs. At this point the MATCH procedure has detected a raw match, of which there are likely to be very many. Since all that is wanted is the longest MUM, some filtering of matches is clearly required.

Candidate MUMs surviving these tests are passed to the UPDATEMUM procedure. Experimental data suggests that, with a minimum match length of five, only about 15% of matches that make it through the three filters are not MUMs.

The cost of scanning for matches between the graph and the new version is O(M(1+)), where M is the number of characters in the graph, and the average length of a match in the suffix tree. In the case that the graph and the new version have no characters in common, exactly M comparisons will be made. On the other hand, in the worst case that the graph and the new version are identical, the program will require NM/2 comparisons. So the time complexity lies between quadratic and linear, depending on the number and distribution of differences between the graph and the new version.

The candidate MUMs, and their frequencies, are stored in a hashtable. When a new candidate MUM arrives from the MATCH procedure, and is already present in the table, its frequency is incremented, otherwise it is added to the table with frequency 1. After the scan of the variant graph is complete, the longest match in the table whose frequency is still 1 will be the best MUM. The cost of each update to the hashtable is constant, and the cost of the final searching through the entries to find the best MUM is an additional linear cost for the whole operation, which has no effect on the time complexity of the algorithm.

Having found the longest MUM, the next step is to merge the new version with the variant graph at that point. In fact, before searching for the MUM the new version has already been added to the variant graph as a special arc connected to the start and end-nodes, as shown in Figure 14.

Figure 14: Direct Alignment, First Iteration

Definition 1

A special arc represents a new version as yet unassimilated into the variant graph, and is distinguished in some way from other non-special arcs. In all other respects a special arc is still part of the graph, and does not violate its definition because it is part of the unique path of the new version. In Figure 14 the C-arc is a special arc, and is drawn in grey to indicate its status. A special arc is associated with a MUM, which represents its best fit to the variant graph. In Figure 14 the longest MUM, marked in bold, is ‘</l><l>d’un antico acquedotto di sguardi,</l><l>la sua curva sacra e ’ (‘</l><l>of an ancient aqueduct of glances,</l><l>its curve sacred and ’), taken from the B-version. It can be merged with the rest of the variant graph by splitting the special arc and the graph at each end of the MUM, then adding new nodes as required, and finally connecting the new arcs to the right nodes. The result is shown in Figure 15.

Figure 15: Direct Alignment, Second Iteration

Figure 16 shows that there are now two subgraphs on either side of the newly aligned section, and two subarcs left over from the special arc. The terms ‘subarc’ and ‘subgraph’ need to be defined because they will be used in the discussion that follows.

Figure 16: Subgraphs and Subarcs

One Small Problem

In the subgraphs of Figure 16, the right subgraph no longer contains a complete A-version. Since the A and B versions on the right are not parallel, it is impossible to tell which portion of A corresponds to the remaining part of the C-version. This is not ideal and does lead to some misalignments. For example, if the text of version C had ended with ‘magic’ instead of ‘solitari’, then, after merging, ‘magic’ would occur twice in the graph, as shown in Figure 17. Order sensitivity has not been entirely banished, as in the POA algorithm.

Figure 17: A Non-optimal Alignment

One solution is to optimise the graph subsequently. Anomalies like these can be merged away by reconsidering all combinations of arcs between shared endpoints once the first pass of the algorithm has produced a rough version of the graph. In this way the two ‘magic’s could be identified as identical variants of one another and merged. In bioinformatics such iterative refinement techniques have produced high quality alignments, for example TCOFFEE ([Notr00]) DIALIGN ([Morg96]) and MAFFT ([Kato02]). However, this problem does not so often arise in practice, and its solution is left open for future work.

In texts written by the author evidence of actual transpositions could refute any computational heuristic. Figure 18 shows a clear case of transposition. What need could there be to calculate the presence of transpositions when their existence can be precisely verified in this way?

Figure 18: A Clear Case of Transposition

Firstly, even in manuscripts written by the author, the presence of a transposition may not always be immediately obvious. The author may choose to completely rewrite the sentence or to indicate the transposition by some other means, such as insertion and deletion, as in Figure 19.

Figure 19: Deletion and Reinsertion of ‘la sera’

Secondly, in the case of multiple drafts, or multiple witnesses of the one text, the discovery of transpositions between physical versions is extremely difficult using any manual method. Even in the case of original manuscripts with corrections by the author, a good automatic method could produce the correct result most of the time. In that case, the editor would only be obliged to manually intervene in a small number of cases. Exactly how such a manual intervention would work is left for the moment as an open question.

The edit cost of a transposition in biological sequences is considered to be constant, e.g. Lopresti and Tomkins (1997). If this was also true of documentary texts then even single characters that matched at the start and end of a long document would have to be considered as potential transpositions, but intuitively that isn’t likely.

A metric for assessing the validity of transpositions may instead be modelled as a threshold rather than as an edit cost. Whether a match is a transposition or not seems to be a function of the match's size and the distance between its two halves. Over what ranges those quantities of size and distance may vary can be determined by experiment.

A sample of 49 real cases of transposition were manually identified in a wide variety of texts:

The results are shown in Figure 20. ‘Size’ means the length of the transposed text and ‘distance’ is the distance in characters between the insertion and deletion points before the transposition was carried out.

Figure 20: Transpositions: Size vs Distance

It is clear from this graph that real transpositions are mostly short-range affairs. (Not considered here are 'transpositions' produced by the reordering of entire works within a single manuscript, as in [Spen03].) An approximate heuristic is d < |m|^φ, where d is the distance between the two halves of the match, |m| its length, and φ is the golden ratio, approximately 1.618034. This ratio is relevant for two reasons: it has been used in the arts since classical times, particularly in architecture, and also it expresses a ratio of relative length: (a+b)/a = a/b = φ. All 49 examples meet this simple criterion. Using such a threshold yields the double benefit of decreasing the amount of calculation while increasing accuracy.

Before describing the algorithm it is necessary to explain the difference between merging a direct alignment with the variant graph and merging a transposition. Conceptually a direct alignment and a transposition are very similar: both identify aligned parts in the new version and the variant graph, but in the transposition case the two halves of the alignment are not opposite one another.

Figure 21: Merging a Transposition

In Figure 21 ‘parrot’ has been identified as a transposition between versions B and C. Merging has been carried out so that the two instances of ‘parrot’ in versions B and C have been separated from their arcs by the addition of new nodes. This creates two new subarcs ‘the raucous’ and ‘swiftly’ on the left, which are both opposite the same subgraph ‘the dog’. The source instance ‘parrot’ in version B may still participate in further transpositions or even in a direct alignment, for example, if the C-version read ‘a surprised parrot’ instead of ‘a surprised lorikeet’. There is no harm in this, and in fact this is one way to detect repetitions. Apart from the addition of a node or two, merging a transposition does not materially change the structure of a variant graph. All that happens is that one half of the match is marked as the parent (in this case the B-version), and the other as the child (version C).

1.  Set currentArc to a special arc made from the new version
2.  Set currentGraph to the variant graph
3.  Attach currentArc to the start and end of currentGraph
4.  Compute the bestMUM between currentArc and currentGraph
5.  Create a priority queue specials of tuples {Arc,MUM,Graph}
6.  WHILE bestMUM is not NULL DO
7.    Merge bestMUM into currentGraph.
8.    FOR each subarc created by the merge DO
9.      Compute the directMUM, leftTransposeMUM and rightTransposeMUM
10.     Set MUM to the best of the three
11.     IF MUM is not NULL THEN
12.       Insert {subarc,MUM,subgraph} into specials; keep it sorted
13.     END IF
14.   END FOR
15.   POP {currentArc, bestMUM, currentGraph} from specials;
16. END WHILE

Algorithm 2 creates a priority queue of triples A_i,M_i,S_j, where A_i is a special arc, M_i is its best MUM with the graph, and S_j is its opposite subgraph. It has a j subscript, not i because, as explained above, more than one special arc may share the same subgraph. The queue is kept sorted on decreasing length of MUM (line 12). In line 7 the longest MUM is used to create 0-2 new subgraphs and subarcs. Each new subarc, together with its subgraph and its best MUM, are added to the specials queue. When the queue no longer contains a valid MUM the program terminates.

The calculation of the two transpose MUMs uses the same process as the direct MUM, except that the variant graph is searched backwards on the left hand side from the start node of the subgraph, and forwards from its end node, as shown in Figure 22. The best such MUM found is recorded as the left or right transpose MUM for the special arc. In choosing the best of three in line 10, if a transpose MUM and a direct MUM are of the same length, the direct MUM is preferred.

This algorithm is admittedly naive, but it is easy to implement. In reality transpositions are often composed of imperfectly matched sections, that is, literal matches interspersed with small variants, deletions and insertions. Refinement of the algorithm to allow for such fuzzy matching is left for future work. For the moment, this simple atomic approach to transposition produces satisfactory results.

There are three operations that have greater than constant cost:

Thus, summing the direct alignment, transposition and sorting costs, the time complexity of merging one version into a variant graph using Algorithm 2, with block transpositions, is O(MN) in the worst case.

Algorithm 2 has been written up as a test program called MergeTester (http://www.itee.uq.edu.au/~schmidt/downloads.html), and has been run on a wide variety of multi-version example texts. The parameters measured are: author, title, average version size in kilobytes, number of versions, approximate date of manuscripts, total time to merge all versions in seconds and average time to merge one version after the first. Times were computed on a 1.66GHz processor. The test program was single-threaded.