Sequence alignment

Concept

Sequence alignment

T-106Home BU-401Threads information · systems

Statement

Detecting homology by dynamic programming.

Why it matters

Every later result in this unit that compares sequences depends, directly or indirectly, on being able to line two sequences up in a way that reflects shared ancestry rather than coincidence. blast-database-search is a heuristic approximation to exactly this problem, built for speed against enormous databases at some cost in guaranteed optimality; phylogenetic-tree-construction requires a multiple-sequence alignment as its starting input before any distance or likelihood calculation can even begin; and hidden-markov-models offer a probabilistic generalisation of the same underlying alignment problem. Sequence alignment, treated exactly and rigorously here via dynamic programming, is therefore the foundational algorithmic tool this entire unit builds on.

The central difficulty alignment must solve is that homologous sequences accumulate not only substitutions but insertions and deletions since diverging from a common ancestor, so a simple position-by-position comparison of two sequences of possibly different length cannot, by itself, detect true homology; an alignment procedure must be able to introduce gaps in a principled way.

Hypotheses

Homologous sequences diverge from a common ancestor by an accumulation of substitutions, insertions, and deletions, so a biologically meaningful alignment must permit gaps, not just character-by-character comparison of equal-length sequences.Without allowing gaps, two truly homologous sequences that differ in length by even a single inserted or deleted residue could never be correctly aligned across the point of that indel, and every position downstream of it would be spuriously misaligned. A fixed scoring scheme (a reward for matches, a penalty for mismatches, a penalty for gaps) can be chosen such that, on average, alignments preserving true evolutionary homology score more highly than alignments that do not.The entire algorithm searches only for the mathematically optimal alignment under whatever scoring scheme is supplied; the biological usefulness of that optimum depends entirely on how well the chosen scores approximate the real, unknown evolutionary substitution and indel process. The optimal alignment under a given, fixed scoring scheme is not guaranteed to be the true evolutionary alignment. Scoring parameters (substitution matrices such as PAM or BLOSUM for proteins, and gap-penalty values) are simplifying approximations, calibrated against known alignments, and different reasonable parameter choices can occasionally yield different "optimal" alignments for the same sequence pair.

Proof

F(i,0) = -i\cdot\text{gap}, \qquad F(0,j) = -j\cdot\text{gap}

Represent the two sequences (lengths \(m\) and \(n\)) as the axes of an \((m{+}1)\times(n{+}1)\) score matrix \(F\); the first row and column are initialised as the cost of aligning a prefix of one sequence entirely against gaps in the other, since there is only one way to do that. A

F(i,j) = \max\big\{\, F(i{-}1,j{-}1) + s(x_i,y_j),\ \ F(i{-}1,j) - \text{gap},\ \ F(i,j{-}1) - \text{gap} \,\big\}

The Needleman–Wunsch recurrence: the best score of aligning prefixes \(x_1..x_i\) and \(y_1..y_j\) is the best of three options — extend by aligning \(x_i\) to \(y_j\) (a match or mismatch, scored \(s(x_i,y_j)\)), or extend by a gap in \(y\) (consuming \(x_i\) against a gap), or extend by a gap in \(x\) (consuming \(y_j\) against a gap). B

\text{The optimal alignment of the full sequences has score } F(m,n), \text{ by the optimal-substructure principle of dynamic programming.}

Because every cell's value is computed as the best extension of already-optimal smaller sub-alignments, no globally optimal alignment can contain a suboptimal sub-alignment of any prefix pair; this guarantees \(F(m,n)\) is the true maximum score over all possible alignments, not merely a good one, without needing to enumerate every possible alignment explicitly. B

\text{Traceback from } F(m,n) \text{ to } F(0,0), \text{ recording which of the three recurrence terms achieved the maximum at each cell, reconstructs an explicit optimal alignment.}

Each step backward corresponds to one column of the final alignment: a diagonal step is a match/mismatch column, a step up or left is a gap column in one sequence; ties at a given cell (two or more terms achieving the same maximum) correspond to multiple, equally optimal alignments. A

\text{Smith\textendash Waterman (local alignment): } F(i,j) = \max\{\,0,\ F(i{-}1,j{-}1)+s(x_i,y_j),\ F(i{-}1,j)-\text{gap},\ F(i,j{-}1)-\text{gap}\,\}\text{, traceback from the matrix maximum.}

Adding the option of resetting to \(0\) prevents the score from going negative, effectively allowing an alignment to "restart" at any point; starting traceback from wherever the overall maximum occurs in the matrix, rather than forcibly from \((m,n)\), finds the best-scoring matching subsequence rather than requiring the two full sequences to align end-to-end. B

Result

F(i,j) = \max\{\,F(i{-}1,j{-}1)+s(x_i,y_j),\ F(i{-}1,j)-\text{gap},\ F(i,j{-}1)-\text{gap}\,\} \ \Rightarrow\ F(m,n) = \text{optimal alignment score}

Reading. Dynamic programming decomposes the alignment problem into optimal smaller alignments of shorter prefixes, building the full solution up cell by cell; the final cell gives the guaranteed-optimal score, and traceback recovers one (or more, if tied) explicit optimal alignment achieving it.

Scope. The direct algorithm requires \(O(mn)\) time and space; it finds the optimal alignment only relative to a chosen scoring scheme (substitution scores and gap penalty), whose biological appropriateness is a separate question from the algorithm's mathematical correctness (Hypotheses).

Corollaries & converses

blast-database-search trades this algorithm's guarantee of finding the true optimal alignment for enormous speed, using a heuristic seed-and-extend strategy to approximate the same scoring idea across databases far too large to align exhaustively against every query.
phylogenetic-tree-construction requires a multiple sequence alignment, an extension of this pairwise procedure to more than two sequences simultaneously, as the necessary input before any distance matrix or likelihood calculation can be computed.
hidden-markov-models generalise this same alignment problem probabilistically, replacing fixed match/mismatch/gap scores with explicit match/insert/delete state probabilities, and replacing a single optimal alignment with a full probability distribution over possible alignments.

Fails without

Drop the optimal-substructure/dynamic-programming property (Hypotheses' scoring-scheme assumption still holding, but computed by brute force): the number of possible ways to align two sequences of length \(m\) and \(n\), once gaps are permitted, grows combinatorially with sequence length, so checking every possible alignment directly becomes computationally infeasible for anything beyond very short sequences — dynamic programming's guarantee of optimality without exhaustive enumeration is what makes exact alignment of realistic sequences tractable at all.
Choose a badly miscalibrated gap penalty: with a gap penalty of zero, the algorithm can insert arbitrarily many gaps to maximise the count of matched positions trivially, producing a biologically meaningless alignment; with an excessively large gap penalty, genuinely homologous sequences containing real insertions or deletions are instead forced into a run of mismatches rather than a gap, since a mismatch becomes cheaper than the (overly punitive) true gap — either extreme produces an alignment that is mathematically optimal under the chosen scores but biologically wrong.

Common errors

Confusing global alignment (Needleman–Wunsch, forcing both full sequences to align end-to-end) with local alignment (Smith–Waterman, finding the best-scoring matching subsequence) and applying the wrong variant to sequences that share only a partial domain of similarity rather than full-length homology.
Assuming the mathematically optimal alignment under a chosen scoring matrix must automatically be the biologically or evolutionarily correct alignment, rather than merely the highest-scoring one under that particular, approximate choice of parameters (Hypotheses, third assumption).
Forgetting that ties are common: multiple distinct alignments can achieve the identical optimal score, so "the" optimal alignment returned by a particular implementation may be only one of several equally valid optima.
Using a single, fixed gap penalty regardless of how evolutionarily distant the two sequences being compared are, rather than adjusting scoring parameters to the expected divergence (a closely related pair and a distantly related pair are not well served by identical scoring assumptions).

Discussion

Saul Needleman and Christian Wunsch published the global dynamic-programming alignment algorithm in 1970, the first systematic method guaranteed to find an optimal alignment between two sequences under a given scoring scheme. Temple Smith and Michael Waterman published the local-alignment variant in 1981, adapting the same recurrence to find the best-scoring matching subsequence rather than requiring the sequences to align end-to-end — the two algorithms remain the standard reference points for, respectively, global and local pairwise sequence alignment.

Protein alignment commonly uses empirically derived substitution matrices (PAM or BLOSUM families) rather than a simple match/mismatch score, reflecting the observation that not all amino acid substitutions are equally likely or equally functionally tolerated during evolution — substituting one hydrophobic residue for a chemically similar one is scored more favourably than substituting it for a charged residue, even though both are "mismatches" in the crude sense.

Common misconception: that finding the mathematically optimal alignment under some scoring scheme settles the biological question of homology outright. It does not: the optimal alignment is only ever optimal relative to the specific, approximate scoring parameters supplied (Hypotheses' third assumption), and different reasonable parameter choices can occasionally favour different alignments for the same sequence pair.

Worked examples

\text{Align } X = \text{"AC"} \text{ and } Y = \text{"AGC"} \text{ with match} = +1,\ \text{mismatch} = -1,\ \text{gap} = -2.

Initialise \(F(0,0..3) = 0,-2,-4,-6\) and \(F(0..2,0) = 0,-2,-4\). Filling the remaining cells via the recurrence (Step 2) gives \(F(1,1){=}1\) (A/A match), \(F(1,2){=}{-1}\), \(F(1,3){=}{-3}\), \(F(2,1){=}{-1}\), \(F(2,2){=}0\), \(F(2,3){=}0\). A

\text{Traceback from } F(2,3)=0: \text{ diagonal (C/C match) from } F(1,2), \text{ then left (gap) from } F(1,1), \text{ then diagonal (A/A match) from } F(0,0).

Reading the traceback path forward gives the alignment A–C against AGC, i.e. \(X\): A – C, \(Y\): A G C; scoring this directly, \((+1) + (-2) + (+1) = 0\), matching \(F(2,3)=0\) exactly and confirming the traceback reconstruction is consistent with the matrix computation. A

X:\ \text{A}\,\text{--}\,\text{C} \qquad Y:\ \text{A}\,\text{G}\,\text{C} \qquad \text{score} = 0

Reading. Even for two very short sequences, dynamic programming correctly identifies that inserting a single gap to preserve both the leading and trailing matches scores better than any alternative alignment without a gap (e.g. aligning A/A and C/G directly would score \((+1)+(-1)+(-2\text{ for the trailing gap}) = -2\), worse than the gapped alignment's score of \(0\)).

Scope. The identical matrix-filling and traceback procedure, at a larger scale, is exactly what real alignment software performs on sequences hundreds or thousands of residues long.

Problems

Using the same scoring scheme as the Worked example (match \(+1\), mismatch \(-1\), gap \(-2\)), compute \(F(1,1)\) for aligning \(X=\text{"GA"}\) against \(Y=\text{"CA"}\), showing all three terms of the recurrence.
Solution
\(F(1,1) = \max\{F(0,0)+s(G,C),\ F(0,1)-2,\ F(1,0)-2\}\). \(s(G,C)=-1\) (mismatch), so the diagonal term is \(0+(-1)=-1\); \(F(0,1)=-2\), so the "up" term is \(-2-2=-4\); \(F(1,0)=-2\), so the "left" term is \(-2-2=-4\). The maximum of \(\{-1,-4,-4\}\) is \(-1\), so \(F(1,1)=-1\), achieved via the diagonal (mismatch) term.
Explain, in terms of the Smith–Waterman recurrence, why allowing \(F(i,j)=\max\{0,\ldots\}\) rather than letting the score go negative is what makes the algorithm find local rather than global matches.
Solution
Without the option to reset to \(0\), a poorly matching region anywhere in the two sequences would drag the cumulative score increasingly negative, and that negative score would persist and continue to be carried forward into every later cell, penalising even a strong matching region that happens to occur later in the sequence. Allowing the score to reset to \(0\) whenever the cumulative score would otherwise go negative effectively lets the alignment "start fresh" at the best-scoring point, so that only the best-scoring contiguous matching subsequence (not the whole sequence, warts and all) determines the final traceback (Step 5).
Two students align the same pair of protein sequences, one using a substitution matrix calibrated for closely related sequences and one using a matrix calibrated for distantly related sequences, and get two different "optimal" alignments. Explain how this is possible without either computation being in error.
Solution
The dynamic-programming algorithm guarantees finding the true mathematically optimal alignment only relative to whichever scoring scheme it is given (Result, Scope); it makes no claim about which scoring scheme is the biologically correct one to use. Because the two students supplied genuinely different substitution scores (reflecting different assumptions about expected evolutionary divergence), the algorithm can correctly return different optimal alignments for each — both computations can be entirely correct as dynamic-programming solutions to two different scoring problems, exactly the caveat raised in the Hypotheses' third assumption and the Discussion's common misconception.