MUSCLE: a multiple sequence alignment method with reduced time and space complexity

Abstract

Background: In a previous paper, we introduced MUSCLE, a new program for creating multiple alignments of protein sequences, giving a brief summary of the algorithm and showing MUSCLE to achieve the highest scores reported to date on four alignment accuracy benchmarks. Here we present a more complete discussion of the algorithm, describing several previously unpublished techniques that improve biological accuracy and / or computational complexity. We introduce a new option, MUSCLE-fast, designed for high-throughput applications. We also describe a new protocol for evaluating objective functions that align two profiles.

Results: We compare the speed and accuracy of MUSCLE with CLUSTALW, Progressive POA and the MAFFT script FFTNS1, the fastest previously published program known to the author. Accuracy is measured using four benchmarks: BAliBASE, PREFAB, SABmark and SMART. We test three variants that offer highest accuracy (MUSCLE with default settings), highest speed (MUSCLE-fast), and a carefully chosen compromise between the two (MUSCLE-prog). We find MUSCLE-fast to be the fastest algorithm on all test sets, achieving average alignment accuracy similar to CLUSTALW in times that are typically two to three orders of magnitude less. MUSCLE-fast is able to align 1,000 sequences of average length 282 in 21 seconds on a current desktop computer.

Conclusions: MUSCLE offers a range of options that provide improved speed and / or alignment accuracy compared with currently available programs. MUSCLE is freely available at http://www.drive5.com/muscle.

Figure 1

Progressive alignment. Sequences are assigned to the leaves of a binary tree. At each internal (i.e., non-leaf) node, the two child profiles are aligned using profile-profile alignment (see Figure 2). Indels introduced at each node are indicated by shaded background.

Figure 2

Profile-profile alignment. Two profiles (multiple sequence alignments) X and Y are aligned to each other such that columns from X and Y are preserved in the result. Columns of indels (gray background) are inserted as needed in order to align the columns to each other. The score for aligning a pair of columns is determined by the profile function, which should assign a high score to pairs of columns containing similar amino acids.

Figure 3

Gap penalties in the SP score This figure shows a multiple alignment of three sequences s, t and u. The SP score is the sum over all pairs of sequences of their pairwise alignment score. The contribution to the SP score from the pair (s, t) is computed by discarding columns in which both sequences have indels (arrows). Such indels are said to be external with respect to the pair. Gaps in the remaining columns (gray background) are assessed affine penalties g + λe where g is the per-gap penalty, λ is the gap length, and e is the gap extension penalty.

Figure 4

Position-specific gap penalties. An alignment of two profiles X and Y. Gaps in sequences t and u are embedded in X. Y contains a single sequence w. The gap in w (gray background) is inserted to align the profiles and is not part of Y. Consider the SP score for this alignment. We need not consider pairs of sequences in X as their scores are unchanged under all possible alignments of X to Y, leaving the inter-profile pairs (s, w), (t, w), (u, w) and (v, w). Note that there is no gap penalty for the pairs (u, w) and (v, w) as these pairs do not have gaps relative to each other. The remaining pairs (t, w) and (u, w) are assessed a penalty g + 3e for the gap in Y. The total over all pairs of open or close penalties due to a gap in Y is thus reduced in proportion to the fraction of sequences in X having a gap with the same open or close position. We incorporate this into the PSP score by using position-specific gap penalties b(x) and t(x). For example, b(x) in column 4 of X is half the default value because half of the sequences in X open a gap in that column. Note that there is no open penalty at the N-terminal and no close penalty at the C-terminal. This causes terminal gaps to receive half the penalty of internal gaps.

Figure 5

Tree comparison. Two trees are compared in order to identify those nodes that have the same branching orders within subtree rotation (white). If a progressive alignment has been created using to the old tree, then alignments at these nodes can be retained as the same result would be produced at those nodes by the new tree. New alignments are needed at the changed (black) nodes only.

Figure 6

Neighbor-joining and UPGMA trees for progressive alignment. Here we show the same set of four sequences and the order in which they will be aligned according to a neighbor-joining tree (above) and a UPGMA tree (below). Notice that t and u are the most closely related pair, but (s, t) and (u, v) are evolutionary neighbors. With neighbor joining, t and u are not aligned to each other until the root, in contrast to UPGMA, which aligns s and t as the first pair.

Figure 7

Additive profiles. The profile functions in MUSCLE require amino acid frequencies for each column. Here we show the alignment of two profiles X and Y, giving a new profile Z. Note that the count n^Z_i for amino acid i in a given column of Z is the sum of the counts in the child profiles, i.e. n^Z_i = n^X_i + n^Y_i. In terms of frequencies, this becomes f ^Z_i = N^Xf ^X_i /N^Z + N^Yf ^Y_i/N^Z, where N^X, N^Y, N^Z are the number of sequences in X, Y and Z respectively. Therefore, given a suitable sequence weighting scheme, it is possible to compute frequencies in Z from the frequencies in X and Y. This avoids the step of building an explicit multiple alignment for Z in order to compute frequencies, as done in CLUSTALW and MAFFT.

Figure 8

Occupancy frequencies in additive profiles. Here we show an alignment of profiles X and Y giving Z. A column C of indels (shaded background) has been inserted at position x in order to align X to Y. To compute the number of gap-extensions in column x of Z, three cases must be considered: (1) a gap-extension in the corresponding column of Y, (2) a gap-open in the preceding column of X, and (3) a gap-extension in the preceding column of X. By enumerating all such cases, it is straightforward to compute the occupancy frequencies in Z from the occupancy frequencies of X and Y, plus the alignment path.

Figure 9

E-strings. (1) The effect of the e-string operator <3,-1,2> on the sequence MQTIF. A positive number n skips n letters, a negative number -n insert n indels. (2) The effect of applying two successive e-strings. In the last line, the result is expressed as a new e-string applied to the original string. (3) We define multiplication on two e-strings as yielding the e-string that is equivalent to applying the two e-strings in order. (4) An alignment path is conventionally represented as a vector of edge types (M, D and I). In this example, MDMIMM, shown above a pairwise alignment, is the path that generates that alignment. The alignment can also be generated by a pair of e-strings (shown to the right). An alignment path is therefore equivalent to a pair of e-strings.

Figure 10

Root alignment construction. Here we show the same progressive alignment as Figure 1. Each edge in the tree is labeled with the e-string for its side of the alignment at the parent node. The e-string needed to insert indels into a sequence in the root alignment can be determined by multiplying e-strings along the path to the root. For example, for sequence LSF, the root e-string is <3,-1,1>*<1,-1,2> = <1,-1,1,-1,1>.

Figure 11

Dimers in the {X,-} alphabet. Gap penalties for the sequence pair (s, u) can be computed be considering all aligned pairs of dimers in the alphabet {X,-}, where X is any amino acid and - is the usual indel symbol. Four cases are highlighted. Note that an aligned pair of identical dimers never contribute a gap penalty as any indels in the dimers are necessarily external, as in the left-most example.

Figure 12

Problem dimer pair. The aligned dimer pair -X ↔ -- causes a problem because its gap penalty contribution cannot be computed without additional information. Note that the first column of indels is external; after this column is discarded, different penalties may be needed, as these two examples show.

Figure 13

Dimer substitution matrix. This matrix specifies the contribution to the total gap penalty for a pair of sequences for each possible pair of aligned dimers. Here, g is the per-gap penalty, e is the gap-extension penalty. The problem case -X ↔ -- is approximated as tg, where t is a tunable parameter.

Figure 14

Discrimination plot for PP2. The x axis is the number of true column pairs with scores ≤ S for some value S, as a fraction of the total number of true pairs; the y axis is the number of false column pairs with scores ≤ S, as a fraction of the total number of false pairs. The databases were constructed from the PP2 test set. Shown are discrimination plots for the log-expectation (LE), log-average (LA), Yona-Levitt (YL), LAMA, and profile sum of pairs (PSP) functions. The LE function shows higher discrimination over the entire range of scores than any other function we tested (complete results not shown). The poor performance of the "standard" PSP function is striking. PSP displays negative discrimination over some of its range where it falls below the diagonal (dashed line).

Figure 15

Discrimination plot for PP. This is similar to Figure 13, except that the database was generated from the PP test set. Here we see an ambiguous result as the discrimination plots for LE and PSP intersect.

Figure 16

Execution time as a function of N. This plot shows the execution time as a function of N (number of sequences) for the tested alignment methods. Input data is from 200 to 1,000 sequences in increments of 200. Average sequence length is 282, maximum length 454.