skbio.alignment.multi_align#

skbio.alignment.multi_align(sequences, /, sub_score=(1.0, -1.0), gap_cost=2.0, free_ends=True, guide_tree=None, ids=None, atol=1e-05, keep_tree=False, keep_distmat=False)[source]#

Perform progressive alignment of multiple sequences.

Added in version 0.7.4.

A versatile workflow for multiple sequence alignment (MSA) by integrating the pairwise alignment engine (pair_align), alignment distance calculation, and UPGMA tree construction following the progressive alignment approach.

Parameters:

sequencesiterable of Sequence, str, or sequence of scalar: Sequences to be aligned. Must be non-empty and ungapped. At least two sequences must be provided.
sub_scoretuple of (float, float), SubstitutionMatrix, or str, optional: Score of a substitution. May be two numbers (match, mismatch), a substitution matrix, or its name. See pair_align for details. Default is (1.0, -1.0).
gap_costfloat or tuple of (float, float), optional: Penalty of a gap. May be one (linear) or two numbers (affine). See pair_align for details. Default is 2.0.
free_endsbool, optional: If True (default), gaps at the sequence terminals are free from penalization.
guide_treeTreeNode, optional: A guide tree determining the merging order of sequences. Must be strictly bifurcating. Tip names must match sequence IDs. If not provided, the function will align all sequence pairs, calculate score-based distances, and compute a guide tree using UPGMA. A provided guide tree will skip these procedures.
idsiterable of str, optional: Unique identifiers in input order to match tips in the provided guide tree. tree. If not provided, the function will use sequence metadata 'id' if present in every sequence and unique, or use ['0', '1', ...] if none is present.
atolfloat, optional: Absolute tolerance in comparing scores of alternative alignment paths. See pair_align for details. Default is 1e-5.
keep_treebool, optional: If True, include the guide tree in the returned object. Default is False.
keep_distmatbool, optional: If True, and if the guide tree is not provided, include the constructed distance matrix in the returned object. Default is False.

Returns:

pathAlignPath: One alignment path in original input sequence order, starting at position zero and consuming every sequence in full.
treeTreeNode or None: Constructed or provided guide tree determining the merging order of sequences (if keep_tree is True).
distmatDistanceMatrix or None: Distance matrix constructed based on pairwise alignment scores and used to compute the guide tree (if keep_distmat is True).

Raises:

ValueError: If inputs are empty or gapped, scores or costs are invalid, identifiers do not match, or a guide tree is not binary with exactly the required tips.
TypeError: If free_ends is not Boolean or guide_tree is not a TreeNode.

See also

pair_align
align_score
AlignPath

Notes

This function implements the classic progressive alignment method for multiple sequence alignment, originally introduced in [1], with later improvements described in [2] and [3]. Compared with the historical method, this implementation represents a refined form of progressive alignment that is often described in educational materials. Specifically, the algorithm consists of the following steps, detailed in reverse order:

Procedures

An alignment of all sequences is constructed by iteratively merging sub-alignments containing one or more sequences. This function adopts the profile alignment approach [3], which aligns two sub-alignments (“profiles”) using the same dynamic programming (DP) algorithm for pairwise sequence alignment (pair_align). Refer to the latter’s documentation for settings and considerations.

The substitution score \(S\) for two matching columns between profiles is calculated as the average substitution score \(s\) across all pairs of characters from the two columns:

\[S = \frac{1}{mn}\sum_{x\in A}\sum_{y\in B}s(x,y)\]

where \(x\) and \(y\) are characters in the two columns \(A\) and \(B\), which contain \(m\) and \(n\) rows (sequences), respectively.

Under the “once a gap, always a gap” rule [1], existing gaps within each profile are treated as neutral characters and assigned a substitution score of 0 with any character. Gap penalties are calculated only for gaps introduced during the DP alignment.

The order of merging is determined by a guide tree. The program traverses the tree in postorder and merges the two child sub-alignments at each internal node. If the guide tree is not explicitly supplied, the program computes one using the UPGMA method (see upgma) from a distance matrix containing all pairwise sequence distances. Each distance is calculated by aligning the two sequences \(a\) and \(b\) using DP (see pair_align) and normalizing the alignment score \(S_{a,b}\) as follows [1]:

\[D = -\ln S_{\mathrm{eff}} = -\ln\left(\frac{S_{a,b} - S_{\mathrm{rand}}} {S_{\mathrm{iden}} - S_{\mathrm{rand}}}\right)\]

where \(S_{\mathrm{iden}} = (S_{a,a} + S_{b,b})/2\) is the average score of the two sequences aligned to themselves. \(S_{\mathrm{rand}}\) is the random score, calculated as [2]:

\[S_{\mathrm{rand}} = \frac{1}{L}\sum_{x\in a}\sum_{y\in b}s(x,y)N_a(x)N_b(y) - G\]

where \(L\) is the length of the alignment, \(G\) is the total gap penalty, and \(N\) is the number of occurrences of a character in the corresponding source sequence.

There are two notes regarding this calculation. First, a factor of 100 is omitted from the effective score \(S_{\mathrm{eff}}\) compared with the original work. Second, to guard against edge cases that would produce undefined, infinite, or negative distances (e.g., when aligning two homopolymers), \(S_{\mathrm{eff}}\) is clipped to the range [1e-6, 1] in this implementation.

Solution quality

The sum-of-pairs (SP) score is the optimality criterion for multiple sequence alignment. This metric can be calculated by applying the align_score function to the resulting alignment. It should be noted that progressive alignment is a heuristic algorithm and the resulting alignment is not guaranteed to be optimal.

Computational efficiency

The algorithm is dominated by the all-vs-all pairwise alignment step, which takes O(n² L²) time for n sequences of comparable length L. Peak memory usage is O(n² + L²) for storing the distance matrix and each DP matrix (the algorithm reuses memory for DP matrices). When a guide tree is supplied, time reduces to O(nL² + n² L), and memory to O(L² + nL).

Terminal gap policy

The function defaults to free_ends=True which prevents terminal gaps from being penalized. This setting is broadly applicable to homologous sequences with incomplete coverage, different domain boundaries, or terminal extensions. However, it can favor short overlaps between weakly related sequences. When sequences are expected to span the same homologous region with defined boundaries (e.g., the coding sequence of a gene), setting free_ends=False is often preferable.

References

[1] (1,2,3)

Feng, D. F., & Doolittle, R. F. (1987). Progressive sequence alignment as a prerequisite to correct phylogenetic trees. Journal of Molecular Evolution, 25(4), 351-360.

[2] (1,2)

Feng, D. F., & Doolittle, R. F. (1996). [21] Progressive alignment of amino acid sequences and construction of phylogenetic trees from them. In Methods in enzymology (Vol. 266, pp. 368-382). Academic Press.

[3] (1,2)

Corpet, F. (1988). Multiple sequence alignment with hierarchical clustering. Nucleic Acids Research, 16(22), 10881-10890.

Examples

>>> from skbio.alignment import multi_align

Align three DNA sequences using default parameters and obtain an alignment path.

>>> from skbio.sequence import DNA
>>> seqs = [DNA('CATTAACGT'),
...         DNA('CGTTACGGT'),
...         DNA('AGTTAACGG')]
>>> path = multi_align(seqs).path
>>> path
<AlignPath, sequences: 3, positions: 11, segments: 7>

Print the aligned sequences.

>>> for seq in path.to_aligned(seqs):
...     print(seq)
CA-TTAACGT-
-CGTTA-CGGT
-AGTTAACGG-

Or convert the alignment path and sequences into a TabularMSA object.

>>> from skbio.alignment import TabularMSA
>>> msa = TabularMSA.from_path_seqs(path, seqs)
>>> msa
TabularMSA[DNA]
----------------------
Stats:
    sequence count: 3
    position count: 11
----------------------
CA-TTAACGT-
-CGTTA-CGGT
-AGTTAACGG-

The quality of the alignment can be evaluated using the align_score function, which calculates the sum-of-pairs (SP) score, and has the same default parameter settings as multi_align does.

>>> from skbio.alignment import align_score
>>> align_score((path, seqs))
7.0

Under the hood, this function performs pairwise alignments, calculates a distance matrix, then infers a guide tree which determines the merging order. The tree and distance matrix can be retained for diagnostic and educational purposes.

>>> path, tree, dm = multi_align(seqs, keep_tree=True, keep_distmat=True)
>>> print(tree.ascii_art())
          /-1
---------|
         |          /-0
          \--------|
                    \-2

>>> print(dm)
3x3 distance matrix
IDs:
'0', '1', '2'
Data:
[[ 0.          0.70444674  0.40215932]
 [ 0.70444674  0.          0.40215932]
 [ 0.40215932  0.40215932  0.        ]]

One can supply a custom guide tree to skip the costly automatic pairwise alignment and tree building process. An accurate tree may improve the alignment quality.

>>> from skbio.tree import TreeNode
>>> tree = TreeNode.read(['((1,2),0);'])
>>> path = multi_align(seqs, guide_tree=tree).path
>>> for seq in path.to_aligned(seqs):
...     print(seq)
CATTAACGT-
CGTTA-CGGT
AGTTAACGG-

>>> align_score((path, seqs))
9.0

By default, sequences match taxa (tip names) of the tree by incremental indices ‘0’, ‘1’, ‘2’… Alternatively, explicit sequence IDs can be defined using the 'id' key of sequence metadata or supplied by the ids parameter of this function.

>>> for seq, id_ in zip(seqs, 'abc'):
...     seq.metadata['id'] = id_
>>> tree = TreeNode.read(['((b,c),a);'])
>>> res = multi_align(seqs, guide_tree=tree)

One can customize the alignment parameters, including substitution scores, gap penalties, and terminal gap policy. Refer to pair_align for details of the parameters.

>>> params = dict(sub_score=(2, -3), gap_cost=(2, 5), free_ends=False)
>>> path = multi_align(seqs, **params).path
>>> for seq in path.to_aligned(seqs):
...     print(seq)
CATTAACGT
CGTTACGGT
AGTTAACGG

Supply the same parameters when calculating the alignment score.

>>> align_score((path, seqs), **params)
4.0

The entire process of reading a multi-FASTA file of original sequences, performing multiple sequence alignment, and writing the aligned sequences into a multi-FASTA file is:

>>> from skbio.io import read as sk_read
>>> seqs = list(
...     sk_read('input.fa', format='fasta', constructor=DNA)
... )
>>> path = multi_align(seqs, **params).path
>>> msa = TabularMSA.from_path_seqs(path, seqs)
>>> msa.write('output.fa')