Chapter 5 Epilogue

Multiple Sequence Alignment

7 min read2 code challenges✓ Read

Watch the lectureMultiple Sequence Alignment

Building a three-dimensional Manhattan

Amino acid sequences of proteins performing the same function are likely to be somewhat similar, but these similarities may be elusive in the case of distant species. You now possess an arsenal of algorithms for aligning pairs of sequences, but if sequence similarity is weak, pairwise alignment may not identify biologically related sequences. However, simultaneous comparison of many sequences often allows us to find similarities that pairwise sequence comparison fails to reveal. Bioinformaticians sometimes say that pairwise alignment whispers and multiple alignment shouts.

We are now ready to use pairwise sequence analysis to build up our intuition for comparison of multiple sequences. In our three-way alignment of A-domains from the introduction, we found 19 conserved columns:

1: YAFDLGYTCMFPVLLGGGELHIVQKETYTAPDEIAHYIKEHGITYIKLTPSLFHTIVNTA
2: -AFDVSAGDFARALLTGGQLIVCPNEVKMDPASLYAIIKKYDITIFEATPALVIPLMEYI
3: IAFDASSWEIYAPLLNGGTVVCIDYYTTIDIKALEAVFKQHHIRGAMLPPALLKQCLVSA

1: SFAFDANFESLRLIVLGGEKIIPIDVIAFRKMYGHTE-FINHYGPTEATIGA
2: -YEQKLDISQLQILIVGSDSCSMEDFKTLVSRFGSTIRIVNSYGVTEACIDS
3: ----PTMISSLEILFAAGDRLSSQDAILARRAVGSGV-Y-NAYGPTENTVLS

However, similarities between A-domains are not limited to these 19 columns, as we can find 10 + 9 + 12 = 31 semi-conservative columns, each of which has two matching amino acids:

1: YAFDLGYTCMFPVLLGGGELHIVQKETYTAPDEIAHYIKEHGITYIKLTPSLFHTIVNTA
2: -AFDVSAGDFARALLTGGQLIVCPNEVKMDPASLYAIIKKYDITIFEATPALVIPLMEYI
3: IAFDASSWEIYAPLLNGGTVVCIDYYTTIDIKALEAVFKQHHIRGAMLPPALLKQCLVSA

1: SFAFDANFESLRLIVLGGEKIIPIDVIAFRKMYGHTE-FINHYGPTEATIGA
2: -YEQKLDISQLQILIVGSDSCSMEDFKTLVSRFGSTIRIVNSYGVTEACIDS
3: ----PTMISSLEILFAAGDRLSSQDAILARRAVGSGV-Y-NAYGPTENTVLS

A multiple alignment of strings v1, ..., vt , also called a t-way alignment, is specified by a matrix having t rows, where the i-th row contains the symbols of vi in order, interspersed with space symbols. We also assume that no column in a multiple alignment contains only space symbols. In the 3-way alignment below, we have highlighted the most popular symbol in each column using upper case letters:

Multiple alignment matrix

The multiple alignment matrix is a generalization of the pairwise alignment matrix to more than two sequences. The three arrays shown below this alignment record the respective number of symbols in ATGTTATA, AGCGATCA, and ATCGTCTC encountered up to a given position. Together, these three arrays correspond to a path in a three-dimensional grid:

(0, 0, 0) → (1, 1, 1) → (2, 2, 2) → (2, 3, 3) → (3, 4, 4) → (4, 5, 5) → (5, 6, 5) → (6, 7, 6) → (7, 7, 7) → (8, 8, 8)

As the alignment graph for two sequences is a grid of squares, the alignment graph for three sequences is a grid of cubes. Every node in the 3-way alignment graph has up to seven incoming edges, as shown in the figure below.

Alignment cube

The score of a multiple alignment is defined as the sum of scores of the alignment columns (or, equivalently, weights of edges in the alignment path), with an optimal alignment being one that maximizes this score. In the case of an amino acid alphabet, we can use a very general scoring method that is defined by a t-dimensional matrix containing 21t entries that describes the scores of all possible combinations of t symbols (representing 20 amino acids and the space symbol). Intuitively, we should reward more conserved columns with higher scores. For more details, see "DETOUR: Scoring Multiple Alignments" in the print companion or at Cogniterra.

Multiple Alignment Problem: Find the highest-scoring alignment between multiple strings under a given scoring matrix.

  • Input: A collection of t strings and a t-dimensional matrix Score.
  • Output: A multiple alignment of these strings whose score (as defined by the matrix Score) is maximized among all possible alignments of these strings.

A straightforward dynamic programming algorithm applied to the t-dimensional alignment graph solves the Multiple Alignment Problem for t strings. For three sequences v, w, and u, we define si, j, k as the length of the longest path from the source (0, 0, 0) to node (i, j, k) in the alignment graph. The recurrence for si, j, k in the three-dimensional case is similar to the recurrence for pairwise alignment:

Multiple alignment recurrence

In the case of t sequences of length n, the alignment graph consists of approximately nt nodes, and each node has up to 2t − 1 incoming edges, yielding a total runtime of O(nt · 2t). As t grows, the dynamic programming algorithm becomes impractical. Many heuristics for suboptimal multiple alignments have been proposed to address this runtime bottleneck.

In the Multiple Longest Common Subsequence Problem, the score of a column of the alignment matrix is equal to 1 if all of the column's symbols are identical, and 0 if even one symbol disagrees.

Code Challenge Solve the Multiple Longest Common Subsequence Problem. Solve at Cogniterra or Rosalind.Solve it with the autograderFree to start in the interactive text

  • Input: Three DNA strings of length at most 10.
  • Output: The length of a longest common subsequence of these three strings, followed by a multiple alignment of the three strings corresponding to such an alignment.

A greedy multiple alignment algorithm

Note that the multiple alignment

AT-GTTaTA
AgCGaTC-A
ATCGT-CTc

induces three pairwise alignments:

AT-GTTaTA       AT-GTTaTA       AgCGaTC-A
AgCGaTC-A       ATCGT-CTc       ATCGT-CTc

But can we work in the opposite direction, combining optimal pairwise alignments into a multiple alignment?

STOP and Think

  1. Does an optimal multiple alignment induce optimal pairwise alignments?
  2. Try combining the following pairwise alignments below into a multiple alignment of the strings CCCCTTTT, TTTTGGGG, and GGGGCCCC.
CCCCTTTT----       ----CCCCTTTT       TTTTGGGG----
----TTTTGGGG       GGGGCCCC----       ----GGGGCCCC

Unfortunately, we cannot always combine optimal pairwise alignments into a multiple alignment because some pairwise alignments may be incompatible (as illustrated by the three pairwise alignments from the previous step, which are reproduced below.

CCCCTTTT----       ----CCCCTTTT       TTTTGGGG----
----TTTTGGGG       GGGGCCCC----       ----GGGGCCCC

The first pairwise alignment implies that CCCC occurs before TTTT in the multiple alignment constructed from these three pairwise alignments. The third pairwise alignment implies that TTTT occurs before GGGG in the multiple alignment. But the second pairwise alignment implies that GGGG occurs before CCCC in the multiple alignment. Thus, CCCC must occur before TTTT, which must occur before GGGG, which must occur before CCCC, a contradiction.

To avoid incompatibility, some multiple alignment algorithms attempt to greedily construct a multiple alignment from pairwise alignments that are not necessarily optimal. The greedy heuristic starts by selecting the two strings having the highest scoring pairwise alignment (among all possible pairs of strings) and then uses this pairwise alignment as a building block for iteratively adding one string at a time to the growing multiple alignment. We align the two closest strings at the first step because they often provide the best chance of building a reliable multiple alignment. For the same reason, we then select the string having maximum score against the current alignment at each stage. But what does it mean to align a string against an alignment of other strings?

An alignment of nucleotide sequences with k columns can be represented as a 4 × k profile matrix like the one shown below, which holds the nucleotide frequencies from each column (amino acid alignments are represented by 20 × k profile matrices). The greedy multiple alignment heuristic adds a string to the current alignment by constructing a pairwise alignment between the string and the profile of the current alignment. As a result, the problem of constructing a multiple alignment of t sequences is reduced to constructing t − 1 pairwise alignments.

STOP and Think Design an algorithm for aligning a string against a profile matrix. How would you score the columns in such an alignment?

Alignment profile

Figure: A profile matrix for a multiple alignment of 10 sequences. Each column of the profile matrix sums to 1 minus the frequency of the space symbol. The most popular nucleotides in each column are shown as uppercase colored letters.

Although greedy multiple alignment algorithms work well for similar sequences, their performance deteriorates for dissimilar sequences because greedy methods may be misled by a spurious pairwise alignment. If the first two sequences picked for building a multiple alignment are aligned in a way that is incompatible with the optimal multiple alignment, then the error in this initial pairwise alignment will propagate all the way through to the final multiple alignment.

After learning how to align multiple sequences, you are now ready to solve a challenge problem that is comparable to the task that Marahiel and his collaborators faced in 1999.

Final Challenge In 1999, Marahiel derived the non-ribosomal code for only 14 of the 20 proteinogenic amino acids (Ala, Asn, Asp, Cys, Gln, Glu, Ile, Leu, Phe, Pro, Ser, Thr, Tyr, Val) because he lacked A-domain sequences for the remaining 6 amino acids. With the availability of many more A-domains, you now have the chance to fill in the gaps in Marahiel’s original paper. Construct the multiple alignment of 397 A-domains, reveal the conservative columns in this alignment, and make your best guess about the signatures encoding all 20 proteinogenic amino acids based on Marahiel’s alignment.Solve it with the autograderFree to start in the interactive text

Download A-Domain Data (.csv format)

Bibliography Notes

Edit distance was introduced by Levenshtein, 1966. The local alignment algorithm described in the text was proposed by Smith and Waterman, 1981. When Doolittle et al., 1983 revealed the similarities between oncogenes and PDGF, they did not know about the Smith-Waterman algorithm. The 3-D structures of firefly luciferase and the A-domain from Bacilus brevis were published by Conti, Franks, and Brick, 1996 as well as Conti et al., 1997. The non-ribosomal code was first presented by Stachelhaus, Mootz, and Marahiel, 1999.