Lesson 5.12

Penalizing Insertions and Deletions in Sequence Alignment

6 min read1 code challenge✓ Read

Watch the lecturePenalizing Insertions and Deletions in Sequence Alignment

Affine gap penalties

We have seen that introducing mismatch and indel penalties can produce more biologically adequate global alignments. However, even with this more robust scoring model, the A-domain alignment that we previously constructed (with indel penalty σ = 4) still reveals only 6 of the 8 conserved purple columns corresponding to the non-ribosomal signatures:

STOP and Think Would increasing the indel penalty from σ = 4 to σ = 10 reveal the biologically correct alignment?

In our previously defined linear scoring model, if σ is the penalty for the insertion or deletion of a single symbol, then σ · k is the penalty for the insertion or deletion of an interval of k symbols. This cost model unfortunately results in inadequate scoring for biological sequences. Mutations are often caused by errors in DNA replication that insert or delete an entire interval of k nucleotides as a single event instead of as k independent insertions or deletions. Thus, penalizing such an indel by σ · k represents an excessive penalty. For example, the alignment on the right is more adequate than the alignment on the left, but they would currently receive the same score.

GATCCAG    GATCCAG

GA-C-AG    GA--CAG

A gap is a contiguous sequence of spaces in a row of an alignment. One way to score gaps more appropriately is to define an affine penalty for a gap of length k as σ + ε · (k − 1), where σ is the gap opening penalty, assessed to the first symbol in the gap, and ε is the gap extension penalty, assessed to each additional symbol in the gap. We select ε to be smaller than σ so that the affine penalty for a gap of length k is smaller than the penalty for k independent single-nucleotide indels (σ · k). For example, with affine gap penalties, the alignment on the left above is penalized by 2σ, whereas the alignment on the right is only penalized by σ + ε.

Alignment with Affine Gap Penalties Problem: Construct a highest-scoring global alignment between two strings (with affine gap penalties).

• Input: Two strings, a matrix score, and numbers σ and ε.
• Output: A highest scoring global alignment between these strings, as defined by the scoring matrix score and by the gap opening and extension penalties σ and ε.

STOP and Think How would you modify the alignment graph to solve this problem?

The figure below illustrates how affine gap penalties can be modeled in the alignment graph by introducing a new “long” edge for each gap.

Figure: Representing gaps in the alignment graph on the left as “long” insertion and deletion edges in the alignment graph on the right. For a gap of length k, the weight of the corresponding long edge is equal to σ + ε · (k − 1).

Since we do not know in advance where gaps should be located, we need to add edges accounting for every possible gap. Thus, affine gap penalties can be accommodated by adding all possible vertical and horizontal edges in the alignment graph to represent all possible gaps. Specifically, we add edges connecting (i, j) to both (i + k, j) and (i, j + k) with weights σ + ε · (k − 1) for all possible gap sizes k, as illustrated in the figure below. For two sequences of length n, the number of edges in the resulting alignment graph modeling affine gap penalties increases from O(n²) to O(n³).

STOP and Think Can you design a DAG with just O(n²) edges to solve the Alignment with Affine Gap Penalties Problem?

Figure: Adding edges corresponding to indels for all possible gap sizes adds a large number of edges to the alignment graph.

Building Manhattan on three levels

The trick to decreasing the number of edges in the DAG for the Alignment with Affine Gap Penalties Problem is to increase the number of nodes. To this end, we will build an alignment graph on three levels; for each node (i, j), we will construct three different nodes: (i, j)_lower , (i, j)_middle , and (i, j)_upper . The middle level will contain diagonal edges of weight score(v_i, w_j) representing matches and mismatches. The lower level will have only vertical edges with weight −ε to represent gap extensions in v, and the upper level will have only horizontal edges with weights −ε to represent gap extensions in w (see figure below).

Figure: Building a three-level graph for alignment with affine gap penalties. The lower level corresponds to gap extensions in v, the middle level corresponds to matches and mismatches, and the upper level corresponds to gap extensions in w.

Life in such a three-level city would be difficult because there is currently no way to move between different levels. To address this issue, we add edges responsible for opening and closing gaps. To model gap opening, we connect each node (i, j)_middle to both (i + 1, j)_lower and (i, j + 1)_upper; we then weight these edges with −σ. Closing a gap does not carry a penalty, and so we introduce zero-weight edges connecting nodes (i, j)_lower and (i, j)_upper with the corresponding node (i, j)_middle at the middle level. As a result, a gap of length k starts and ends at the middle level and is charged −σ for the first symbol, −ε for each subsequent symbol, and 0 to close the gap, producing a total penalty of σ + ε · (k − 1), as desired. The figure below illustrates how the path from this lesson's third step traverses the three-level graph.

Exercise BreakProve that the number of edges in the three-level graph described in the figure below is at most 7 · n · m for sequences of length n and m.

The DAG we have just constructed may be complicated, but it uses only O(n · m) edges for sequences of length n and m, and a longest path in this graph still constructs an optimal alignment with affine gap penalties. The three-level alignment graph translates into the system of three recurrence relations shown below. Here, lower_{i, j} , middle_{i, j} , and upper_{i, j} , are the lengths of the longest paths from the source node to (i, j)_lower , (i, j)_middle , and (i, j)_upper , respectively.

The variable lower_{i, j} computes the score of an optimal alignment between the i-prefix of v and the j-prefix of w ending with a deletion (i.e., a vertical edge), whereas the variable upper_{i, j} computes the score of an optimal alignment of these prefixes ending with an insertion (i.e., a horizontal edge), and the variable middle_{i, j} computes the score of an optimal alignment ending with a match or mismatch. The first term in the recurrences for lower_{i, j} and upper_{i, j} corresponds to extending the gap, whereas the second term corresponds to initiating the gap.

Code Challenge Solve the Alignment with Affine Gap Penalties Problem. Solve at Cogniterra or Rosalind.Solve it with the autograderFree to start in the interactive text

• Input: Two amino acid strings v and w (each of length at most 100).
• Output: The maximum alignment score between v and w, followed by an alignment of v and w achieving this maximum score.

STOP and Think Compute an optimal alignment with affine gap penalties for the A-domains considered in the beginning of this lesson and reproduced below. How does varying the gap opening and extension penalties affect the quality of the alignment?

String 1: YAFDLGYTCMFPVLLGGGELHIVQKETYTAPDEIAHYIKEHGITYIKLTPSLFHTIVNTASFAFDANFESLRLIVLGGEKIIPIDVIAFRKMYGHTEFINHYGPTEATIGA

String 2: AFDVSAGDFARALLTGGQLIVCPNEVKMDPASLYAIIKKYDITIFEATPALVIPLMEYIYEQKLDISQLQILIVGSDSCSMEDFKTLVSRFGSTIRIVNSYGVTEACIDS

Exercise Break Design an algorithm for computing optimal local (rather than global) alignment with affine gap penalties.