MSARC: Multiple sequence alignment by residue clustering
© Modzelewski and Dojer; licensee BioMed Central Ltd. 2014
Received: 1 December 2013
Accepted: 6 April 2014
Published: 16 April 2014
Progressive methods offer efficient and reasonably good solutions to the multiple sequence alignment problem. However, resulting alignments are biased by guide-trees, especially for relatively distant sequences.
We propose MSARC, a new graph-clustering based algorithm that aligns sequence sets without guide-trees. Experiments on the BAliBASE dataset show that MSARC achieves alignment quality similar to the best progressive methods.
Furthermore, MSARC outperforms them on sequence sets whose evolutionary distances are difficult to represent by a phylogenetic tree. These datasets are most exposed to the guide-tree bias of alignments.
MSARC is available at http://bioputer.mimuw.edu.pl/msarc
Determining the alignment of a group of biological sequences is among the most common problems in computational biology. The dynamic programming method of pairwise sequence alignment can be readily extended to multiple sequences but requires the computation of an n-dimensional matrix to align n sequences. Since the size of such a matrix is exponential with respect to n, the time and space complexity of this method is exponential too.
Progressive alignment offers a substantial complexity reduction at the cost of possible loss of the optimal solution. Within this approach, subset alignments are sequentially pairwise aligned to build the final multiple alignment. The order of pairwise alignments is determined by a guide-tree representing the phylogenetic relationships between sequences.
There are two drawbacks of the progressive alignment approach. First, the accuracy of the guide-tree affects the quality of the final alignment. This problem is particularly important in the field of phylogeny reconstruction, because multiple alignment acts as a preprocessing step in most prominent methods of inferring a phylogenetic tree of sequences. It has been shown that, within this approach, the inferred phylogeny is biased towards the initial guide-tree [2, 3].
Second, only sequences belonging to currently aligned subsets contribute to their pairwise alignment. Even if a guide-tree reflects correct phylogenetic relationships, these alignments may be inconsistent with remaining sequences and the inconsistencies are propagated to further steps. To address this problem, in recent programs [4–8] progressive alignment is usually preceded by consistency transformation (incorporating information from all pairwise alignments into the objective function) and/or followed by iterative refinement of the multiple alignment of all sequences. Moreover, recently several strategies avoiding guide trees altogether were also proposed [9–11].
Experiments on the BAliBASE dataset  show that our approach is competitive with the best progressive methods and significantly outperforms most non-progressive algorithms. Moreover, MSARC is the best aligner for sequence sets with very low levels of conservation. This feature makes MSARC a promising preprocessing tool for phylogeny reconstruction pipelines.
MSARC aligns sequence sets in several steps. In a preprocessing step, following Probalign , stochastic alignments are calculated for all pairs of sequences and consistency transformation is applied to resulting posterior probabilities of residue correspondences. Transformed probabilities, called residue alignment affinities, represent weights of an alignment grapha.
Pairwise stochastic alignment
The concept of stochastic (or probability) alignment was proposed in . Given a pair of sequences, this framework defines statistical weights of their possible alignments. Based on these weights, for each pair of residues from both sequences, the posterior probability of being aligned may be computed.
A consensus of highly weighted suboptimal alignments was shown to contain pairs with significant probabilities that agree with structural alignments despite the optimal alignment deviating significantly. Mückstein et al.  suggest the use of the method as a starting point for improved multiple sequence alignment procedures.
where β corresponds to the inverse of Boltzmann’s constant and should be adjusted to the match/mismatch scoring function s(x,y) (in fact, β simply rescales the scoring function).
Let P(a i ∼b j ) denote the posterior probability that residues a i and b j are aligned.
Here we use the notation for an alignment of the sequence prefixes a1⋯a i and b1⋯b j , and for an alignment of the sequence suffixes a i ⋯a m and b j ⋯b n . Analogously, Zi,j is the partition function over the prefix alignments and is the (reverse) partition function over the suffix alignments.
The reverse partition function can be calculated using the same recursion in reverse, starting from the ends of the aligned sequences.
In this case insertions and deletions must be separated by at least one match/mismatch position. This variant was proposed by Miyazawa  and applied in the Probalign  and MSAProbs  aligners.
Let us regard probabilities P(a i ∼b j ) as a representation of a bipartite graph with weighted edges, i.e. a graph with residues from both sequences as nodes and edges joining each a i with each b j .
Given a set S of k sequences to be aligned, we would like to analogously represent their residue alignment affinity by a k-partite weighted graph. It may be obtained by joining pairwise alignment graphs for all pairs of S-sequences. However, separate computation of edge weights for each pair of sequences does not exploit information included in the remaining alignments. Thus we decided to address this problem with a so called consistency transformation[4, 7], successfully used in progressive methods.
where w x y are weights specifying the relative contribution to the transformation of a sequence pair xy.
where · stands for matrix multiplication.
The idea behind the above formula is that the sum of a row/column of a matrix P a b yields the probability that the corresponding residue is aligned to one in the other sequence (not a gap). If sequences a and b are similar, alignments with fewer gaps are preferred, so (at least for the shorter sequence) most of the sums are close to 1. Consequently, the w a b is close to 1 as well. On the other hand, weights are much closer to 0 for pairs of dissimilar sequences.
Thus w a b measures the similarity of sequences a and b. Therefore sequences c that are similar to a and b contribute to more significantly than others.
The consistency transformation may be iterated any number of times, but excessive iterations blur the structure of residue affinity. Following Probalign  and ProbCons , MSARC performs two iterations by default.
Towards this objective, MSARC applies top-down hierarchical clustering (see Figure 2). Within this approach, the alignment graph is recursively split into two parts until no ambiguous cluster is left. Each partition step results from a single cut through all sequences, so clusterings are conflict-free at each step of the procedure. Consequently, the final clustering represents a proper multiple alignment.
Optimal clustering is expected to maximize residue alignment affinity within clusters and minimize it between them. Therefore, the partition selection in recursive steps of the clustering procedure should minimize the sum of weights of edges cut by the partition. This is in fact the objective of the well-known problem of graph partitioning, i.e. dividing graph nodes into roughly equal parts such that the sum of weights of edges connecting nodes in different parts is minimized.
The Fiduccia-Mattheyses algorithm  is an efficient heuristic for the graph partitioning problem. After selecting an initial, possibly random partition, it calculates for each node the change in cost caused by moving it between parts, called gain. Subsequently, single nodes are greedily moved between partitions based on the maximum gain and gains of remaining nodes are updated. The process is repeated in passes, where each node can be moved only once per pass. The best partition found in a pass is chosen as the initial partition for the next pass. The algorithm terminates when a pass fails to improve the partition. Grouping single moves into passes helps the algorithm to escape local optima, since intermediate partitions in a pass may have negative gains. An additional balance condition is enforced, disallowing movement from a partition that contains less than a minimum desired number of nodes.
Fiduccia-Mattheyses algorithm needs to be modified in order to deal with alignment graphs. Mainly, residues are not moved independently; since the graph topology has to be maintained, moving a residue involves moving all the residues positioned between it and a current cut point on its sequence. This modification implies further changes in the design of data structures for gain processing. Next, the sizes of parts in considered partitions cannot differ by more than the maximum cluster size in a final clustering, i.e., the number of aligned sequences. This choice implies minimal search space containing partitions consistent with all possible multiple alignments. In the initial partition sequences are cut in their midpoints.
Therefore we decided to add a refinement step, following the method used in ProbCons . Sequences are split into two groups and the groups are pairwise re-aligned. Re-alignment is performed using the Needleman-Wunsch algorithm with the score for each pair of positions defined as the sum of posterior probabilities for all non-gap pairs and zero gap-penalty. First each sequence is re-aligned with the remaining sequences, since such division is very efficient in removing superfluous spaces. Next, several randomly selected sequence subsets are re-aligned against the rest.
Figures 6(cd) show the results of refining the alignments from Figures 6(ab). Refinement removed superfluous spaces from the clustering process and optimized the alignment. Note that the final post-refinement alignments turned out to be the same for both Fiduccia-Mattheyses and multilevel method of graph partitioning.
Löytynoja and Goldman argue in  that progressive methods tend to force alignments of non-homologous sequence fragments inserted in corresponding locations of aligned sequences. This tendency leads to systematic errors of the downstream analyses in phylogenetic pipelines, including overestimation of substitution and deletion events. Unfortunately, iterative refinement may be one of possible source of such effects. Therefore the number of iterations in subset re-alignment step in MSARC is adjustable, in particular the whole step may be turned off.
Let n denote a number of sequences to align and let l be their maximum length. Both time and space complexities of stochastic alignment are .
Computations in the other steps use data structures for sparse matrices, so the complexity depends on the number c of non-zero values per row/column. This number depends on the cutoff parameter t c (entries <t c are set to 0), namely c≤1/t c . However, we observe that c tends to be much lower than this bound, e.g. c rarely exceeds 5 for the default t c =0.01.
MSARC implementation of consistency transformation requires time. Space complexity of this and the remaining steps is dominated by sparse matrices and equals .
The time complexity of one pass of the Fiduccia-Mattheyses algorithm on whole sequences is . We observe that the algorithm converges after very few passes, but it is hard to prove a reasonable asymptotic bound. The complexity of the whole clustering is asymptoticly equal to the complexity of the main partition step.
The time complexity of iterative refinement belongs to the class .
Benchmark data and methodology
MSARC was tested against the BAliBASE 3.0 benchmark database . It contains manually refined reference protein alignments based on 3D structural superpositions. Each alignment contains core-regions that correspond to the most reliably alignable sections of the alignment. Alignments are divided into five sets designed to evaluate performance on varying types of problems:
Equidistant sequences with two different levels of conservation very divergent sequences (<20% identity)
medium to divergent sequences (20−40% identity)
Families aligned with a highly divergent “orphan” sequence
Subgroups with <25% residue identity between groups
Sequences with N/C-terminal extensions
BAliBASE 3.0 also provides a program comparing given alignments with a reference one. Alignments are scored according to two metrics. A sum-of-pairs score (SP) showing the ratio of residue pairs that are correctly aligned, and a total column (TC) score showing the ratio of correctly aligned columns. Both scores can be applied to full sequences or just the core-regions.
We decided to present results based on core-region scores only, since the corresponding sections of the reference alignments are most reliable. Moreover, results for full sequence scores are very similar.
Benchmarking MSARC variants
Two steps of MSARC algorithm: stochastic alignment and iterative refinement follow the respective steps in Probalign . Therefore we decided to set a bunch of related parameters to Probalign’s defaults. Namely, MSARC was run with Gonnet 160 similarity matrix , gap penalties of −22, −1 and 0 for gap open, extension and terminal gaps respectively, β=0.2, a cut-off value for posterior probabilities of 0.01 (values smaller than the cutoff are set to 0 and operations designed for sparse matrices are used in order to speed up computations), two iterations of the consistency transformation and 100 iterations of iterative refinement.
On the other hand, we decided to evaluate three parameters that seem to be crucial for steps specific for MSARC approach. First, residue clustering may be performed with basic or multilevel Fiduccia-Mattheyses algorithm. Second, weighted or unweighted consistency transformation may be applied. Third, stochastic pairwise alignment may be based on equations (5)-(8) (i.e. stochastic version of classical Gotoh algorithm) or equations (6) and (7) may be replaced with equations (9) and (10), respectively. The modified formula disallows consecutive insertions and deletions, as is done in Probalign and MSAProbs.
Evaluation of MSARC variants
Comparison to other aligners
MSARC results were compared to CLUSTAL Ω [1, 21] ver. 1.1.0, DIALIGN-T  ver. 0.2.2, DIALIGN-TX  ver. 1.0.2, MAFFT  ver. 6.903, MUSCLE  ver. 3.8.31, MSAProbs  ver. 0.9.7, Probalign  ver. 1.4, ProbCons  ver. 1.12, T-Coffee  ver. 9.02, FSA  ver. 1.15.7 and PicXAA  ver. 1.03. All the programs were executed with their default parameters.
Comparison of multiple sequence alignment methods
1 2:1 5
Significance of differences in BAliBASE 3.0 SP/TC scores
Probalign aligns separately the first half of the sequences (blue and green) and the second half of the sequences (from yellow to red). Next, the prefixes of the second group are aligned with the suffixes of the first group, propagating an error within a yellow sub-alignment.
MSAprobs aligns separately the dark blue, light blue and red sequences. Next the blue sub-alignments are aligned together. The resulting alignment has erroneously inserted gaps near the right ends of dark blue sequences. This error is propagated in the next step, where the suffix of the blue alignment is aligned with the prefix of the red alignment. Finally, the single violet sequence is added to the alignment, splitting it in two.
For both programs, alignment errors introduced in the earlier steps are propagated to the final alignment. On the other hand, the non-progressive strategy used in MSARC yields a reasonable approximation of the reference alignment (see Figure 7(ab)).
The progressive principle has dominated multiple alignment algorithms for nearly 20 years. Throughout this time, many groups have dedicated their effort to refine its accuracy to the current state. Other approaches were omitted due to high computational complexity and/or unsatisfactory quality. However, recently several non-progressive methods were proposed. Two of them, PicXAA and MSARC proved to be competitive with the best progressive approaches. Moreover, both programs outperform progressive methods on sequence sets with evolutionary distances that are difficult to represent by a phylogenetic tree.
Despite the algorithmic novelty, the non-progressive approaches to multiple alignment are interesting preprocessing tools for phylogeny reconstruction pipelines. The objective of these procedures is to infer the structure of a phylogenetic tree from a given sequence set. Multiple alignment is usually the first pipeline step. When alignment is guided by a tree, the reconstructed phylogeny is biased towards this tree. In order to minimize this effect, some phylogenetic pipelines alternately optimize a tree and an alignment [24–26]. The unbiased alignment process of MSARC may simplify this procedure and improve the reconstruction accuracy, especially in the most problematic cases.
MSARC has also the potential for quality improvements. Alternative methods of computing residue alignment affinities could be used to improve the accuracy of both MSARC and Probalign based methods. Other approaches to alignment graph partitioning may also lead to improvements in the accuracy of MSARC, for example a better method of pairing residues for multilevel coarsening than the currently used naive consecutive neighbors merging.
The main disadvantage of MSARC is its computational complexity, especially in the case of the multilevel scheme variant (MSARC is ∼2.5× slower than MSAProbs, the MSARC variant with multilevel scheme is even slower). This is the cost of avoiding the progressive approach.
a Our notion of alignment graph slightly differs from the one of Kececioglu : removing edges between clusters transforms the former into the latter.
This work was supported by the Polish Ministry of Science and Higher Education [N N519 652740].
- Thompson JD, Higgins DG, Gibson TJ:Clustal w: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994, 22 (22): 4673-4680.View ArticlePubMedPubMed CentralGoogle Scholar
- Wong KM, Suchard MA, Huelsenbeck JP:Alignment uncertainty and genomic analysis. Science. 2008, 319 (5862): 473-476. doi:10.1126/science.1151532View ArticlePubMedGoogle Scholar
- Löytynoja A, Goldman N:Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. Science. 2008, 320 (5883): 1632-1635. doi:10.1126/science.1158395.View ArticlePubMedGoogle Scholar
- Notredame C, Higgins DG, Heringa J:T-coffee: A novel method for fast and accurate multiple sequence alignment. J Mol Biol. 2000, 302 (1): 205-217. doi:10.1006/jmbi.2000.4042.View ArticlePubMedGoogle Scholar
- Edgar RC:Muscle: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004, 32 (5): 1792-1797. doi:10.1093/nar/gkh340View ArticlePubMedPubMed CentralGoogle Scholar
- Katoh K, Toh H, Miyata T, :Mafft version 5 improvement in accuracy of multiple sequence alignment. Nucleic Acids Res. 2005, 33 (2): 511-518. doi:10.1093/nar/gki198.View ArticlePubMedPubMed CentralGoogle Scholar
- Do CB, Mahabhashyam MSP, Brudno M, Batzoglou S:Probcons: Probabilistic consistency-based multiple sequence alignment. Genome Res. 2005, 15 (2): 330-340. doi:10.1101/gr.2821705,View ArticlePubMedPubMed CentralGoogle Scholar
- Roshan U, Livesay DR:Probalign: multiple sequence alignment using partition function posterior probabilities. Bioinformatics. 22 (22): 2715-2721. doi:10.1093/bioinformatics/btl472.Google Scholar
- Subramanian AR, Weyer-Menkhoff J, Kaufmann M, Morgenstern B:Dialign-t: an improved algorithm for segment-based multiple sequence alignment. BMC Bioinformatics. 2005, 6: 66-doi:10.1186/1471-2105-6-66.View ArticlePubMedPubMed CentralGoogle Scholar
- Bradley RK, Roberts A, Smoot M, Juvekar S, Do J, Dewey C, Holmes I, Pachter L:Fast statistical alignment. PLoS Comput Biol. 2009, 5 (5): 1000392-10.1371/journal.pcbi.1000392. doi:10.1371/journal.pcbi.1000392.View ArticleGoogle Scholar
- Sahraeian SME, Yoon B-J:Picxaa: greedy probabilistic construction of maximum expected accuracy alignment of multiple sequences. Nucleic Acids Res. 2010, 38 (15): 4917-4928. doi:10.1093/nar/gkq255.View ArticlePubMedPubMed CentralGoogle Scholar
- Thompson JD, Koehl P, Ripp R, Poch O:Balibase 3.0: latest developments of the multiple sequence alignment benchmark. Proteins. 2005, 61 (1): 127-136. doi:10.1002/prot.20527,View ArticlePubMedGoogle Scholar
- Fiduccia CM, Mattheyses RM:A linear-time heuristic for improving network partitions. Proceedings of the 19th Design Automation Conference. DAC ’82. 1982, 175-181.http://dl.acm.org/citation.cfm?id=800263.809204], Piscataway, NJ, USA: IEEE Press,Google Scholar
- Miyazawa S:A reliable sequence alignment method based on probabilities of residue correspondences. Protein Eng. 1995, 8 (10): 999-1009.View ArticlePubMedGoogle Scholar
- Mückstein U, Hofacker IL, Stadler PF:Stochastic pairwise alignments. Bioinformatics. 2002, 18 (Suppl 2): 153-160. 10.1093/bioinformatics/18.suppl_2.S153.View ArticleGoogle Scholar
- Yu YK, Hwa T:Statistical significance of probabilistic sequence alignment and related local hidden markov models. J Comput Biol. 2001, 8 (3): 249-282. doi:10.1089/10665270152530845.View ArticlePubMedGoogle Scholar
- Gotoh O:An improved algorithm for matching biological sequences. J Mol Biol. 1982, 162 (3): 705-708.View ArticlePubMedGoogle Scholar
- Liu Y, Schmidt B, Maskell DL:Msaprobs: multiple sequence alignment based on pair hidden markov models and partition function posterior probabilities. Bioinformatics. 1964, 26 (16): 1958-1964. doi:10.1093/bioinformatics/btq338.View ArticleGoogle Scholar
- Hendrickson B, Leland R:A multilevel algorithm for partitioning graphs. Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM). Supercomputing ’95. 1995, doi:10.1145/224170.224228. [http://doi.acm.org/10.1145/224170.224228], New York, NY, USA: ACM, doi:10.1145/224170.224228.Google Scholar
- Gonnet GH, Cohen MA, Benner SA:Exhaustive matching of the entire protein sequence database. Science. 1992, 256 (5062): 1443-1445.View ArticlePubMedGoogle Scholar
- Sievers F, Wilm A, Dineen D, Gibson TJ, Karplus K, Li W, Lopez R, McWilliam H, Remmert M, Söding J, Thompson JD, Higgins DG:Fast, scalable generation of high-quality protein multiple sequence alignments using clustal omega. Mol Syst Biol. 2011, 7: 539-doi:10.1038/msb.2011.75.View ArticlePubMedPubMed CentralGoogle Scholar
- Subramanian AR, Kaufmann M, Morgenstern B:Dialign-tx: greedy and progressive approaches for segment-based multiple sequence alignment. Alg Mol Biol. 2008, 3: 6-doi:10.1186/1748-7188-3-6.View ArticleGoogle Scholar
- Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O:New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of phyml 3.0. Syst Biol. 2010, 59 (3): 307-321. doi:10.1093/sysbio/syq010.View ArticlePubMedGoogle Scholar
- Redelings BD, Suchard MA:Joint bayesian estimation of alignment and phylogeny. Syst Biol. 2005, 54 (3): 401-418. doi:10.1080/10635150590947041.View ArticlePubMedGoogle Scholar
- Lunter G, Miklós I, Drummond A, Jensen JL, Hein J:Bayesian coestimation of phylogeny and sequence alignment. BMC Bioinformatics. 2005, 6: 83-doi:10.1186/1471-2105-6-83.View ArticlePubMedPubMed CentralGoogle Scholar
- Liu K, Raghavan S, Nelesen S, Linder CR, Warnow T:Rapid and accurate large-scale coestimation of sequence alignments and phylogenetic trees. Science. 2009, 324 (5934): 1561-1564. doi:10.1126/science.1171243.View ArticlePubMedGoogle Scholar
- Kececioglu J:The maximum weight trace problem in multiple sequence alignment. Proceedings of the 4th Symposium on Combinatorial Pattern Matching (CPM). Lecture Notes in Computer Science. Berlin Heidelberg: Springer,1993, 106-119.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.