BRICS RS-98-3 Pedersen et al.: Comparison of Coding DNA BRICS

Basic Research in Computer Science

Comparison of Coding DNA

Christian N. S. Pedersen Rune B. Lyngsø Jotun Hein BRICS Report Series RS-98-3 ISSN 0909-0878 January 1998 Copyright \begin{matrix}c \\ ◯ \end{matrix} 1998, BRICS, Department of Computer Science

University of Aarhus. All rights reserved.

Reproduction of all or part of this work is permitted for educational or research use on condition that this copyright notice is included in any copy.

See back inner page for a list of recent BRICS Report Series publications. Copies may be obtained by contacting:

BRICS Department of Computer Science University of Aarhus Ny Munkegade, building 540 DK–8000 Aarhus C Denmark

Telephone: +45 8942 3360
Telefax: +45 8942 3255
Internet: BRICS@brics.dk

BRICS publications are in general accessible through the World Wide Web and anonymous FTP through these URLs:

http://www.brics.dk
ftp://ftp.brics.dk

This document in subdirectory RS/98/3/

Comparison of coding DNA

Christian N. S. Pedersen^{1}^{∗}  Rune Lyngsø^{1∗}  Jotun Hein^{2∗}

1 BRICS^{†}, Department of Computer Science, University of Aarhus 2 Institute of Biological Sciences, University of Aarhus

Abstract

We discuss a model for the evolutionary distance between two coding DNA sequences which specializes to the DNA/protein model proposed in Hein [3]. We discuss the DNA/protein model in details and present a quadratic time algorithm that computes an optimal alignment of two coding DNA sequences in the model under the assumption of affine gap cost. The algorithm solves a conjecture in [3] and we believe that the constant factor of the running time is sufficiently small to make the algorithm feasible in practice.

1 Introduction

A straightforward model of the evolutionary distance between two coding DNA sequences is to ignore the encoded proteins and compute the distance in some evolutionary model of DNA. We say that such a model is a DNA level model. The evolutionary distance between two sequences in a DNA level model can most often be formulated as a classical alignment problem and be efficiently computed using a dynamic programming approach [7, 9, 10, 11].

It is well known that proteins evolve slower than its coding DNA, so it is usually more reliable to describe the evolutionary distance based on a comparison of the encoded proteins rather than on a comparison of the coding DNA itself. Hence, most often the evolutionary distance between two coding DNA sequences is modeled in terms of amino acid events, such as substitution of a single amino acid and insertion-deletion of consecutive amino acids, necessary to transform the one encoded protein into the other encoded protein. We say that such a model is a protein level model. The evolutionary distance between two coding DNA sequences in a protein level model can most often be formulated as a classical alignment problem of the two encoded proteins. Even though a protein level model is usually more reliable than a DNA level model, it falls short because it postulates that all insertions and deletions on the underlying DNA occur at codon boundaries and because it ignores similarities on the DNA level.

∗E-mail: cstorm@brics.dk, rlyngsoe@brics.dk and jotun@pop.bio.aau.dk

†

Basic Research in Computer Science,

Centre of the Danish National Research Foundation.

In this paper we present a model of the evolutionary distance between two coding DNA sequences in which a nucleotide event is penalized by the change it induces on the DNA as well as on the encoded protein. The model is a natural combination of a DNA level model and a protein level model. The DNA/protein model introduced in Hein [3, 5] is a biological reasonable instance of the general model in which the evolution of coding DNA is idealized to involve only substitution of a single nucleotide and insertion-deletion of a multiple of three nucleotides. Hein [3, 5] presents an O(n^{2}m^{2}) time algorithm for computing the evolutionary distance in the DNA/protein model between two sequences of length n and m. This algorithm assumes certain properties of the cost function. We discuss these properties and present an O(nm) time algorithm that solves the same problem under the assumption of affine gap cost. The practicality of an algorithm not only depends on the asymptotic running time but also on the constant factor hidden by the use of O-notation. To determine the distance between two sequences of length n and m our algorithm computes 400nm table entries. Each computation involves a few additions, table lookups and comparisons. We believe the constant factor is sufficiently small to make the algorithm feasible in practice.

The problem of comparing coding DNA is also discussed by Arvestad [1] and Hua, Jiang and Wu [6]. The models discussed in these papers are inspired by the DNA/protein model in Hein [3, 5] but differs in the interpretation of gap cost. A heuristic algorithm for solving the alignment problem in the DNA/protein model is described In Hein [4]. A related problem of how to compare a coding DNA sequence with a protein has been discussed in [8, 12].

The rest of this paper is organized as follows: In section 2 we introduce and discuss the DNA/protein model. In section 3 we present a simple alignment algorithm. In section 4 we present quadratic time alignment algorithm. In section 5 we discuss improvements and future work.

2 The DNA/protein model

Let a = a_{1}a_{2}a_{3} … a_{3n−2}a_{3n−1}a_{3n} be a coding sequence of DNA of length 3n with a reading frame starting at a_{1}. We introduce the notation a_{1}^{i}a_{2}^{i}a_{3}^{i} to denote the ith codon a_{3i−2}a_{3i−1}a_{3i} and the notation A_{i} to describe the amino acid coded by the ith codon. The amino acid sequence A = A_{1}A_{2} … A_{n} describes the protein coded by a. An evolutionary event e on the DNA that transforms a to a′ will also change the encoded protein from A to A′. Since amino acids are coded by several codons, the proteins A and A′ might be identical. The cost of e should reflect the changes on the DNA as well as the changes on the encoded protein.

cost(a →^{e} a′) = cost_{d}(a →^{e} a′) + cost_{p}(A →^{e} A′)  (1)

We say that cost_{d}(a →^{e} a′) is the cost of e on the DNA level and that cost_{p}(A →^{e} A′) is the cost of e on the protein level. In this paper we assume that the DNA level cost and the protein level cost are combined by addition but other combination functions f : ℝ × ℝ → ℝ could of course also be considered. The cost of a sequence E of evolutionary events e_{1}, e_{2}, …, e_{k}

transforming a^{(0)} to a^{(k)} as a^{(0)} →^{e_{1}} a^{(1)} →^{e_{2}} a^{(2)} →^{e_{3}} ⋯ →^{e_{k}} a^{(k)} is defined as

some function of the costs of each event. In the rest of this paper we will assume that this function is the sum of the costs of each event.

cost(a^{(0)} →^{E} a^{(k)}) = ∑_{i=1}^{k} cost(a^{(i−1)} →^{e_{i}} a^{(i)})  (2)

We define the distance between two coding sequences of DNA a and b according to the parsimony principle as the minimum cost of a sequence of evolutionary events which transforms a to b.

dist(a, b) = min{cost(a →^{E} b) | E is a sequence of events}  (3)

In order to compute dist(a, b) we have to specify the set of allowed evolutionary events and define the cost of each event on the DNA level as well as on the protein level. The choice of evolutionary events and cost function influences both the biological relevance of the distance measure and the computational complexity of computing the distance.

The DNA/protein model introduced in [3] can be described as an instance of the above model where we idealize the evolution of a coding sequence of DNA to involve only substitution of a single nucleotide and insertion or deletion of a multiple of three consecutive nucleotides. The reason why gap lengths are restricted to a multiple of three is because an insertion or deletion of length not divisible by three changes the reading frame. This is called a frame shift and it may change the entire remaining amino acid sequence as illustrated in figure 1. Frame shifts are believed to be rare biological events, so it is not unreasonable to leave them out of the model.

The DNA level cost in the DNA/protein model is defined in the classical way by specifying a substitution cost and a gap cost. The cost of substituting a nucleotide σ with σ′ is c_{d}(σ, σ′) for some metric c_{d} on nucleotides and the cost of inserting or deleting 3k consecutive nucleotides is g_{d}(3k) for some convex^{1} function g_{d} : ℕ → ℝ^{+}. We note that cost_{d}(a →^{e} a′) in the

1A convex function fulfills that f(i + j) ≤ f(i) + f(j). A convex gap cost function implies that an insertion-deletion of a consecutive block of nucleotides is best explained as a single event.

Thr A C G Val G T G Thr A C G Gln C A A Ile A T T ⋯
A C G Thr G A T Gly G A C Asp G C A Ala A A T Asn T ⋯

Figure 1: An insertion or deletion of length not divisible by three changes the reading frame.

DNA/protein model only depends on e (and not a and a′), so the order of events is irrelevant on the DNA level. The protein level cost of a nucleotide event which changes the encoded protein from A to A′ should somehow reflect the difference between protein A and protein A′. Hence, cost_{p}(A →^{e} A′) is defined as the minimum cost of a distance alignment of A and A′ where we allow substitution of a single amino acid and insertion-deletion of consecutive amino acids. The substitution cost is given by a metric c_{p} on amino acids and the gap cost is given by a convex function g_{p} : ℕ → ℝ^{+} fulfilling that g_{p}(0) = 0. We use dist_{p}(A, A′) to denote the minimum cost of such an alignment of A and A′ and note that dist_{p} is commutative.

Except for the restriction on the length of an insertion or a deletion the DNA/protein model allows the traditional set of symbol based nucleotide events. An alignment of two sequences describes a set of substitution or insertion-deletion events necessary to transform the one sequence into the other sequence. The set of events is usually described by a matrix or a path in a graph as illustrated in figure 2.

[ \begin{matrix}T & T & G & C & T & − & − & − & C \\ T & − & − & − & C & A & T & G & C \end{matrix} ]  T T G C T C T C A T G C

Figure 2: An alignment can be described by a matrix or a path in the alignment graph. The above alignment describes three matches and two gaps of combined length six.

The cost of an alignment is the optimal cost of any sequence of the events described by the alignment. Hence, the evolutionary distance in the DNA/protein model between two coding DNA sequences is the cost of an optimal alignment in the model. If the cost of any sequence of events is independent of the order but only depends on the set of events, then an optimal alignment can be computed efficiently using dynamic programming [7, 9, 10, 11]. The protein level cost in the DNA/protein model however depends on the order of events, so we cannot use a classical alignment algorithm to compute an optimal alignment in the DNA/protein model. In the rest of this section we will examine the protein level cost further in order to be able to formulate restrictions that allow an efficient alignment algorithm in the DNA/protein model.

A nucleotide event affects nucleotides in one or more consecutive codons. Since a nucleotide event in the DNA/protein model is assumed not to change the reading frame, then only the amino acids encoded by the affected codons are affected by the nucleotide event. A nucleotide event thus changes protein A = UXV to protein A′ = UX′V where X and X′ are the amino acids affected by nucleotide event. Hein [3] implicitly assumes that dist_{p}(A, A′) is the cost of an alignment of X and X′ describing the minimum number of insertions or deletions. This property is essential to the formulation of alignment algorithms but it inflicts some restrictions on the substitution cost c_{p} and the gap cost g_{p} as summarized in lemma 1.

[ \begin{matrix}A_{1} & A_{2} & ⋯ & A_{i−1} \\ A_{1} & A_{2} & ⋯ & A_{i−1} \end{matrix} \begin{matrix}A_{i} \\ A′_{i} \end{matrix} \begin{matrix}A_{i+1} & ⋯ & A_{n} \\ A_{i+1} & ⋯ & A_{n} \end{matrix} ]

(a) A substitution in the ith codon.

[ \begin{matrix}A_{1} & A_{2} & ⋯ & A_{i−1} & A_{i} \\ A_{1} & A_{2} & ⋯ & A_{i−1} & A_{i} \end{matrix} \begin{matrix}A_{i+1} & ⋯ & A_{i+k} \\ − & ⋯ & − \end{matrix} \begin{matrix}A_{i+k+1} & ⋯ & A_{n} \\ A_{i+k+1} & ⋯ & A_{n} \end{matrix} ]

(b) An insertion-deletion of 3k nucleotides affecting exactly k codons.

[ \begin{matrix}A_{1} & A_{2} & ⋯ & A_{i−1} \\ A_{1} & A_{2} & ⋯ & A_{i−1} \end{matrix} \begin{matrix}A_{i} & ⋯ & A_{j−1} & A_{j} & A_{j+1} & ⋯ & A_{i+k} \\ − & ⋯ & − & υ & − & ⋯ & − \end{matrix} \begin{matrix}A_{i+k+1} & ⋯ & A_{n} \\ A_{i+k+1} & ⋯ & A_{n} \end{matrix} ]

(c) An insertion-deletion of 3k nucleotides affecting k + 1 codons. The remaining amino acid υ is matched with one of the amino acids affected by the deletion.

Figure 3: The protein level cost of a nucleotide event can be determined by considering only the amino acids affected by the event.

Lemma 1 Assume a nucleotide event changes A = UXV to A′ = UX′V. Let n = |A| and k = ||A| − |A′||. If there for any amino acids σ, τ and for all 0 < l ≤ n − k exists 0 ≤ j ≤ k such that c_{p}(σ, τ) + g_{p}(j) + g_{p}(k − j) ≤ g_{p}(l) + g_{p}(l + k), then dist_{p}(A, A′) is the cost of an alignment describing exactly k insertions or deletions. Furthermore dist_{p}(A, A′) only depends on X and X′.

Proof. We will argue that the assumption stated in the lemma implies that dist_{p}(A, A′) is the cost of one of the alignments illustrated in figure 3. These alignments all describe the minimum number of insertions or deletions and only the sub-alignment of X and X′, as illustrated by the shaded parts, contributes to the cost. We split the argumentation depending on the event. Since dist_{p}(A, A′) is equal to dist_{p}(A′, A) then the cost of an insertion transforming A to A′ is equal to the cost of a deletion transforming A′ to A. We thus only consider substitutions and deletions.

A substitution of a nucleotide in the ith codon of A transforms A_{i} to A′_{i}. The alignment in figure 3(a) describes no insertion-deletions and has cost c_{p}(A_{i}, A′_{i}). Any other alignment of A and A′ must describe an equal number of insertions and deletions, so by convexity of g_{p} the cost is at least 2 · g_{p}(l) for some 0 < l ≤ n. The assumption in the lemma implies that c_{p}(A_{i}, A′_{i}) ≤ 2 · g_{p}(l) for any 0 < l ≤ n, so the protein level cost of the substitution is c_{p}(A_{i}, A′_{i}).

A deletion of 3k nucleotides affects k or k + 1 consecutive codons. If the deletion affects exactly k codons then it transforms A = A_{1}A_{2} ⋯ A_{n} to A′ = A_{1}A_{2} ⋯ A_{i}A_{i+k+1}A_{i+k+2} ⋯ A_{n}. Any alignment of A and A′ must describe l insertions and l + k deletions for some 0 ≤ l ≤ n − k, so the cost is at least g_{p}(l) + g_{p}(l + k). The alignment in figure 3(b) describes k deletions and has cost g_{p}(k). The assumption in the lemma and the convexity of g_{p} implies that g_{p}(k) ≤ g_{p}(j) + g_{p}(k − j) ≤ g_{p}(l) + g_{p}(l + k) for all l > 0, so the protein level cost of a deletion affecting k codons is g_{p}(k).

If the deletion affects k + 1 codons, say by deleting the 3k nucleotides a_{3}^{i}a_{1}^{i+1}a_{2}^{i+1}a_{3}^{i+1} ⋯ a_{1}^{i+k}a_{2}^{i+k}, then it transforms A = A_{1}A_{2} ⋯ A_{n} to A′ = A_{1}A_{2} ⋯ A_{i−1}υA_{i+k+1} ⋯ A_{n} where υ is the amino acid coded by a_{1}^{i}a_{2}^{i}a_{3}^{i+k}. We say that υ is the remaining amino acid and a_{1}^{i}a_{2}^{i}a_{3}^{i+k} is the remaining codon. Any alignment of A and A′ describing exactly k deletions must align υ with A_{i+j} for some 0 ≤ j ≤ k, so by convexity of g_{p} the cost is at least g_{p}(j) + c_{p}(A_{i+j}, σ) + g_{p}(k − j). The alignment in figure 3(c) illustrates one of the k + 1 alignments of A and A′ where υ is aligned with an affected amino acids and all non-affected amino acids are aligned. Such an alignment describes exactly k deletions and the cost of the optimal alignment among them has cost

min_{j=0,1,... ,k}{g_{p}(j) + c_{p}(A_{i+j}, σ) + g_{p}(k − j)},  (4)

and is thus optimal for any alignment describing exactly k deletions. Any other alignment of A and A′ must describe l insertions and l+k deletions for some 0 < l ≤ n − k, so the cost is at least g_{p}(l) + g_{p}(l + k). The assumption in the lemma implies that the cost given by expression (4) is less than or equal to g_{p}(l)+g_{p}(l +k), so that the protein level cost of a deletion affecting k + 1 codons is given expression (4). □

The assumption in lemma 1 is sufficient to ensure that we can compute the protein level cost of a nucleotide event efficiently, but the formulation of the lemma it is to general to make the assumption necessary. The following example however suggests when the assumption is necessary. Consider a deletion of three nucleotides that transforms the six amino acids ABEFCD to ABGCD, i.e. X = EF and X′ = G. If we assume that c_{p}(E, G) ≤ c_{p}(F, G) then the cost of the alignment in figure 4 (left) is c_{p}(E, G) + g_{p}(1) while the cost of the alignment in figure 4 (right) is g_{p}(2)+g_{p}(1). If the assumption in lemma 1 does not hold then g_{p}(2)+g_{p}(1) might be less than c_{p}(E, G)+g_{p}(1) because c_{p}(E, G) can be arbitrary large. Hence, the protein level cost of the deletion would not be the cost of an alignment describing the minimum number of insertion-deletions.

[ \begin{matrix}A & B & E & F & C & D \\ A & B & G & − & C & D \end{matrix} ]  [ \begin{matrix}A & B & E & F & − & C & D \\ A & B & − & − & G & C & D \end{matrix} ]

Figure 4: Two alignments of the amino acids ABEFCD and ABGCD.

An often used gap cost is the affine gap cost function g_{p}(k) = α_{p}+β_{p}k for some α_{p}, β_{p} ≥ 0. With affine gap cost the assumption in lemma 1 becomes c_{p}(σ, τ) + α_{p} + β_{p}k ≤ 2 · α_{p} + β_{p}(k + 2l) for any amino acids σ, τ and all lengths 0 < l ≤ n − k. This simplifies to c_{p}(σ, τ) ≤ α_{p} + 2β_{p} for all amino acids σ, τ. This is a biological reasonable assumption since it reflects that insertions and deletions are rare compared to substitutions.

3 A simple alignment algorithm

Let a_{1}a_{2} ⋯ a_{3n} and b_{1}b_{2} ⋯ b_{3m} be two coding sequences of DNA. We want to compute an optimal alignment of a and b in the DNA/protein model. An alignment of a and b describes a set of events between a and b but not the order of the events. The DNA level cost of an alignment is easy to determine as it is independent of the order, so it is just the sum of the cost of all the events described by the alignment. The protein level cost however depends on the order of the events. An obvious way to determine the protein level cost of an alignment is to minimize over all possible sequences of the events described by the alignment. This method is, however, not feasible in practice due to the factorial number of possible sequences one has to consider. If the substitution cost c_{p} and gap cost g_{p} fulfill lemma 1, then the protein level cost of a nucleotide event only depends on the affected codons. We can use this property to decompose the computation of the protein level cost of an alignment into smaller subproblems.

We decompose the alignment into codon alignments. A codon alignment is a minimal part of the alignment that corresponds to a path connecting two nodes (3i′, 3j′) and (3i, 3j) in the alignment graph. The alignment in figure 2 is a codon alignment. We can decompose any alignment uniquely into codon alignments as illustrated in figure 6. The assumption that an insertion or deletion has length a multiple of three implies that there are eleven distinct types of codon alignments as illustrated in figure 5.

[figure: Type 1:; Type 2:; Type 3:; Type 4:; Type 5:; Type 6:; Type 7:; Type 8:; Type 9:; Type 10:; Type 11:]

Figure 5: The eleven types of codon alignments. The length of an insertion or deletion is a multiple of three.

We observe that nucleotide events described by two different codon alignment in the decomposition of an alignment do not affect the same codons. Hence, the protein level cost of the alignment is the sum of the protein level cost of each of the codon alignments in the decomposition. The protein level cost of a codon alignment can be determined by minimizing over the possible sequences of the up to five nucleotide events described by the codon alignment. Hein [3] describes how the decomposition into codon alignments makes it possible to compute the optimal alignment of a and b in time O(n^{2}m^{2}). The algorithm can be summarized as follows.

Let D(i, j) denote the cost of an optimal alignment of a_{1}a_{2} ⋯ a_{3i} and b_{1}b_{2} ⋯ b_{3j}. If i ≤ 0 or j ≤ 0 then we define D(i, j) to be infinity. An optimal alignment of a_{1}a_{2} ⋯ a_{3i} and b_{1}b_{2} ⋯ b_{3j} can be decomposed into codon alignments ca_{1}, ca_{2}, …, ca_{k}. If the last codon alignment ca_{k} is an alignment of a_{3i′+1}a_{3i′+2} ⋯ a_{3i} and b_{3j′+1}b_{3j′+2} ⋯ b_{3j} for some (i′, j′) < (i, j)^{2}, then D(i, j) is equal to D(i′, j′) + cost(ca_{k}). We can thus compute D(i, j) by minimizing the expression D(i′, j′) + cost(ca) over all (i′, j′) < (i, j) and the up to four possible codon alignments ca of a_{3i′+1}a_{3i′+2} ⋯ a_{3i} and b_{3j′+1}b_{3j′+2} ⋯ b_{3j}. If we assume that D(i′, j′) is known for all (i′, j′) < (i, j), then Hein [3] argues that we can compute D(i, j) in time O(ij). By dynamic programming this implies that we can compute D(n, m) in time O(n^{2}m^{2}).

2We say that (i′, j′) < (i, j) iff i′ ≤ i ∧ j′ ≤ j ∧ (i′ ≠ i ∨ j′ ≠ j).

[figure]

Figure 6: An alignment of two sequences decomposed into codon alignments.

4 An improved algorithm

The improved algorithm is similar to the simple algorithm in the sense that we compute the cost of an optimal alignment by minimizing over all possible last codon alignments. We define D^{t}(i, j) to be the cost of an optimal alignment of a_{1}a_{2} ⋯ a_{3i} and b_{1}b_{2} ⋯ b_{3j} under the assumption that the last codon alignment is of type t. If i ≤ 0 or j ≤ 0 then we define D^{t}(i, j) to be infinity. We compute D(i, j) as

D(i, j) = min_{t=1,2,... ,11} D^{t}(i, j).  (5)

In the rest of this paper we assume that the amino acid substitution cost c_{p} and gap cost g_{p} fulfill lemma 1 and that the combined gap cost function g(k) = g_{d}(3k) + g_{p}(k) is affine α + βk for some α, β ≥ 0. These assumptions make it possible to compute D^{t}(i, j) in constant time if D(k, l) has been computed (by the above expression) for all (k, l) < (i, j). This implies that we can compute D(n, m) in time O(nm).

We divide the explanation of how to compute D^{t}(i, j) in constant time according to the number of gaps within a codon (internal gaps) described by a codon alignment of type t. Codon alignments of type 1–3 describe no internal gaps, codon alignments of type 4–7 describe one internal gap and codon alignments of type 8–11 describe two internal gaps. In each case we use the fact that c_{p} and g_{p} fulfill lemma 1 to compute D^{t}(i, j) as the cost of the last codon alignment (of type t) plus the cost of the remaining alignment.

Codon alignments with no internal gaps

A codon alignment of type 1 describes three substitutions. We introduce c_{p}^{∗} : {A, C, G, T}^{3} × {A, C, G, T}^{3} → ℝ such that c_{p}^{∗}(σ_{1}σ_{2}σ_{3}, τ_{1}τ_{2}τ_{3}) is the cost of a codon alignment of type 1 of codon σ_{1}σ_{2}σ_{3} and codon τ_{1}τ_{2}τ_{3}. The cost c_{p}^{∗}(σ_{1}σ_{2}σ_{3}, τ_{1}τ_{2}τ_{3}) is the minimum cost^{3} of a sequence of the three substitutions σ_{1} → τ_{1}, σ_{2} → τ_{2} and σ_{3} → τ_{3}. The cost D^{1}(i, j) is the cost of the last codon alignment of type 1 plus the cost of the remaining alignment.

D^{1}(i, j) = D(i − 1, j − 1) + c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, b_{1}^{j}b_{2}^{j}b_{3}^{j})  (6)

A codon alignment of type 2 or type 3 describes a gap between codons. Since the combined gap cost function is affine we can use the technique introduced in [2] saying that a gap ending in (i, j) is either a continuation of an existing gap ending in (i − 1, j) or (i, j − 1), or a start of a new gap.

D^{2}(i, j) = min{D(i, j − 1) + α + β, D^{2}(i, j − 1) + β}  (7)
D^{3}(3i, 3j) = min{D(i − 1, j) + α + β, D^{3}(i − 1, j) + β}  (8)

Codon alignments with one internal gap

We will describe how to compute D^{6}(i, j). The other three cases where the last codon alignment describes one internal gap are handled similarly. The last codon alignment of type 6 describes three substitutions and one deletion. If the deletion has length k (a deletion of 3k nucleotides), then the last codon alignment is an alignment of a_{1}^{i′}a_{2}^{i′}a_{3}^{i′} ⋯ a_{1}^{i}a_{2}^{i}a_{3}^{i} and b_{1}^{j}b_{2}^{j}b_{3}^{j} where i′ = i − k. This is illustrated in figure 7.

[figure: a_{1}^{i′}; a_{2}^{i′}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; a_{3}^{i}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; (3i,3j); (3i−3k−3,3j−3)]

Figure 7: The last codon alignment of type 6.

The cost D^{6}(i, j) is the cost of the last codon alignment plus the cost of the remaining alignment. The cost of the remaining alignment is D(i − k − 1, j − 1) and the cost of the last codon alignment is the minimum cost of a sequence of the four described events. Any sequence of the four events can be divided into three steps: The substitutions occurring before the deletion, the deletion and the substitutions occurring after the deletion. Figure 8 illustrates the three steps of the evolution of a_{1}^{i′}a_{2}^{i′}a_{3}^{i′} ⋯ a_{1}^{i}a_{2}^{i}a_{3}^{i} to b_{1}^{j}b_{2}^{j}b_{3}^{j}. The nucleotides x_{1}, x_{2} and x_{3} are the result of the up to three substitutions before the deletion. For example, if the substitution a_{1}^{i′} → b_{1}^{j} occurs before the deletion, then x_{1} is b_{1}^{j}, otherwise it is a_{1}^{i′}. We say that x_{1} ∈ {a_{1}^{i′}, b_{1}^{j}}, x_{2} ∈ {a_{2}^{i′}, b_{2}^{j}} and x_{3} ∈ {a_{3}^{i}, b_{3}^{j}} are the status of the three substitutions before the deletion. The status of the substitutions before the deletion is used extensively in the computation of the cost of the last codon alignment.

3We use the term cost to denote the DNA level cost plus the protein level cost.

[figure: a_{1}^{i′}; a_{2}^{i′}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; a_{3}^{i}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; x_{1}; x_{2}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{del}; x_{1}; x_{2}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}]

Figure 8: The evolution of a_{1}^{i′}a_{2}^{i′}a_{3}^{i′} ⋯ a_{1}^{i}a_{2}^{i}a_{3}^{i} to b_{1}^{j}b_{2}^{j}b_{3}^{j} described by the last codon alignment.

To compute the cost of the three substitutions in the last codon alignment we use the status of the substitutions before the deletion to split the occurrence into two steps. We say that the substitutions a_{1}^{i′} → x_{1}, a_{2}^{i′} → x_{2} and a_{3}^{i} → x_{3} occur before the deletion and that the substitutions x_{1} → b_{1}^{j}, x_{2} → b_{2}^{j} and x_{3} → b_{3}^{j} occur after the deletion. Since an identical substitution has cost zero, then the cost of the three substitutions in the last codon alignment is equal to the cost of the six substitutions obtained by splitting the occurrence of the three substitutions into two steps. The substitutions occurring before the deletion change codon a_{1}^{i′}a_{2}^{i′}a_{3}^{i′} to x_{1}x_{2}a_{3}^{i′} and codon a_{1}^{i}a_{2}^{i}a_{3}^{i} to a_{1}^{i}a_{2}^{i}x_{3}. The substitutions occurring after the deletion change codon x_{1}x_{2}x_{3} to b_{1}^{j}b_{2}^{j}b_{3}^{j}. We recall that the cost of changing codon σ_{1}σ_{2}σ_{3} to codon τ_{1}τ_{2}τ_{3} by a sequence of the substitutions σ_{1} → τ_{1}, σ_{2} → τ_{2} and σ_{3} → τ_{3} is c_{p}^{∗}(σ_{1}σ_{2}σ_{3}, τ_{1}τ_{2}τ_{3}), so the cost of the three substitutions is

cost(subs) = c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, x_{1}x_{2}a_{3}^{i′}) +

c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}x_{3}) + c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}). (9)

The cost of the deletion of 3k nucleotides in the last codon alignment is the sum of the DNA level cost g_{d}(3k) and the protein level cost as given by expression (4). By using the combined gap cost function g(k) = g_{d}(3k) + g_{p}(k) = α + βk and our knowledge of the status of the three substitutions before the deletion, i.e. the remaining codon of the deletion, we can formulate this sum as

cost(del) = min \{ \begin{matrix}α + βk + c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3})^{4} \\ 2α + βk + min_{0<l<k} c_{p}(a_{1}^{i−l}a_{2}^{i−l}a_{3}^{i−l}, x_{1}x_{2}x_{3}) \\ α + βk + c_{p}(x_{1}x_{2}a_{3}^{i′}, x_{1}x_{2}x_{3}) \end{matrix} . (10)

The cost of the deletion depends on the deletion length, the remaining codon x_{1}x_{2}x_{3} and a witness. The witness encodes the amino acid aligned with the remaining amino acid (see figure 3(c)). The witness can be the end-codon a_{1}^{i}a_{2}^{i}x_{3}, the start-codon x_{1}x_{2}a_{3}^{i′} or one of the internal codons a_{1}^{i−l}a_{2}^{i−l}a_{3}^{i−l} for some 0 < l < k.

4We use c_{p}(σ_{1}σ_{2}σ_{3}, τ_{1}τ_{2}τ_{3}) as a convenient notation for c_{p}(σ, τ) where σ and τ are the amino acids coded by the codons σ_{1}σ_{2}σ_{3} and τ_{1}τ_{2}τ_{3} respectively.

We compute D^{6}(i, j) by minimizing the cost of the last codon alignment plus the cost of the remaining alignment over all possible last codon alignments of type 6. This is done by minimizing the sum cost(subs)+cost(del)+ D(i − k − 1, j − 1) over all possible combinations of deletion length k and remaining codon x_{1}x_{2}x_{3}. A combination of deletion length k and remaining codon x_{1}x_{2}x_{3} is possible if x_{1} ∈ {a_{1}^{i′}, b_{1}^{j}}, x_{2} ∈ {a_{2}^{i′}, b_{2}^{j}} and x_{3} ∈ {a_{3}^{i}, b_{3}^{j}} where i′ = i − k. The terms c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}x_{3}) and c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}) of cost(subs) do not depend on the deletion length, so we can split the minimization as

D^{6}(i, j) = min_{x_{1}x_{2}x_{3}}{c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}x_{3}) +

c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}) + D_{x_{1}x_{2}x_{3}}^{6}(i, j)} (11)

where

D_{x_{1}x_{2}x_{3}}^{6}(i, j) = min_{0<k<i}{D(i − k − 1, j − 1) +

c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, x_{1}x_{2}a_{3}^{i′}) + cost(del)} (12)

is the minimum cost of the terms that depend on both the deletion length and the remaining codon under the assumption that the remaining codon is x_{1}x_{2}x_{3}. The cost D_{x_{1}x_{2}x_{3}}^{6}(i, j) is defined if there exist a deletion length k such that k and x_{1}x_{2}x_{3} is a possible combination of deletion length and remaining codon. If we expand the term cost(del) we get

D_{x_{1}x_{2}x_{3}}^{6}(i, j) = min_{0<k<i}{len_{x_{1}x_{2}}^{6}(i, j, k) +

min \{ \begin{matrix}c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}) \\ α + min_{0<l<k} c_{p}(a_{1}^{i−l}a_{2}^{i−l}a_{3}^{i−l}, x_{1}x_{2}x_{3}) \\ c_{p}(x_{1}x_{2}a_{3}^{i′}, x_{1}x_{2}x_{3}) \end{matrix} \} (13)

where

len_{x_{1}x_{2}}^{6}(i, j, k) = D(i − k − 1, j − 1) + c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, x_{1}x_{2}a_{3}^{i′}) + α + βk (14)

is the cost of the remaining alignment plus the part of the cost of the last codon alignment that does not depend on the codon a_{1}^{i}a_{2}^{i}a_{3}^{i} and the witness. The cost len_{x_{1}x_{2}}^{6}(i, j, k) is defined if x_{1} ∈ {a_{1}^{i′}, b_{1}^{j}} and x_{2} ∈ {a_{2}^{i′}, b_{2}^{j}} where i′ = i − k.

Since only x_{1} and x_{2} depend on the deletion length (because a_{1}^{i′} and a_{2}^{i′} depend on the deletion length), then there are at most 32 possible remaining codons x_{1}x_{2}x_{3}. If we can compute D_{x_{1}x_{2}x_{3}}^{6}(i, j) in constant time for each of these possible remaining codons, then we can compute D^{6}(i, j) in constant time. To compute D_{x_{1}x_{2}x_{3}}^{6}(i, j) we must determine a combination of witness and deletion length that minimizes the cost. This combination must be one of the four cases illustrated in figure 9, so all we have to do is to compute the minimum cost of these four cases. The cost of case 1–3 is obtained by simplifying expression (13) for a particular witness and deletion length.

[figure: x_{1}; x_{2}; a_{3}^{i′}; z_{1}; z_{2}; z_{3}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; witness]

(a) Case 1

[figure: x_{1}; x_{2}; a_{3}^{i′}; z_{1}; z_{2}; z_{3}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; witness]

(b) Case 2

[figure: x_{1}; x_{2}; z_{3}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; witness]

(c) Case 3

[figure: x_{1}; x_{2}; a_{3}^{i′}; z_{1}; z_{2}; z_{3}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; \underline{not} witness]

(d) Case 4

Figure 9: The four cases in the computation of D_{x_{1}x_{2}x_{3}}^{6}(i, j). We use z_{1}z_{2}z_{3} as notation for a_{1}^{i−1}a_{2}^{i−1}a_{3}^{i−1}.

Case 1: The end-codon is the witness and the deletion length is at least

one. The cost is: min_{0<k<i} len_{x_{1}x_{2}}^{6}(i, j, k) + c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}).

Case 2: The last internal codon is the witness and the deletion length is at

least two. The cost is: min_{1<k<i} len_{x_{1}x_{2}}^{6}(i, j, k)+α+c_{p}(a_{1}^{i−1}a_{2}^{i−1}a_{3}^{i−1}, x_{1}x_{2}x_{3}).

Case 3: The start-codon is the witness and the deletion length is one. The cost is: len_{x_{1}x_{2}}^{6}(i, j, 1) + c_{p}(x_{1}x_{2}a_{3}^{i−1}, x_{1}x_{2}x_{3}).

Case 4: The witness is neither the end-codon nor the last internal codon and the deletion length is at least two. If the witness of D_{x_{1}x_{2}x_{3}}^{6}(i − 1, j) is not the end-codon a_{1}^{i−1}a_{2}^{i−1}x_{3}, then by optimality of D_{x_{1}x_{2}x_{3}}^{6}(i − 1, j) this witness must also be the witness of case 4. If this is the case then the cost of case 4 is D_{x_{1}x_{2}x_{3}}^{6}(i − 1, j) + β.

Case 4 suggests that we can use dynamic programming to keep track of D_{x_{1}x_{2}x_{3}}^{6}(i, j) under the assumption that the end-codon is not the witness, i.e. use dynamic programming to keep track of the minimum cost of case 2– 4. We introduce tables F_{x_{1}x_{2}x_{3}}^{6} corresponding to the 64 combinations of x_{1}x_{2}x_{3}. We maintain that if x_{1}x_{2}x_{3} is a possible remaining codon and the end-codon a_{1}^{i}a_{2}^{i}x_{3} is not the witness of D_{x_{1}x_{2}x_{3}}^{6}(i, j), then F_{x_{1}x_{2}x_{3}}^{6}(i, j) is equal to D_{x_{1}x_{2}x_{3}}^{6}(i, j). If we define F_{x_{1}x_{2}x_{3}}^{6}(0, j) to infinity, then we can compute table entry (i, j) as

F_{x_{1}x_{2}x_{3}}^{6}(i, j) = min \{ \begin{matrix}cost of Case 2 \\ cost of Case 3 \\ F_{x_{1}x_{2}x_{3}}^{6}(i − 1, j) + β \end{matrix}  (15)

In order to compute the cost of case 1 and 2 in constant time we maintain the minimum of len_{x_{1}x_{2}}^{6}(i, j, k) over k by dynamic programming. We introduce tables L_{x_{1}x_{2}}^{6} corresponding to the 16 combinations of x_{1}x_{2} such that L_{x_{1}x_{2}}^{6}(i, j) is equal to min_{0<k<i} len_{x_{1}x_{2}}^{6}(i, j, k). If we define L_{x_{1}x_{2}}^{6}(0, j) to infinity, then we can compute table entry (i, j) as

L_{x_{1}x_{2}}^{6}(i, j) = min \{ \begin{matrix}len_{x_{1}x_{2}}^{6}(i, j, 1) \\ L_{x_{1}x_{2}}^{6}(i − 1, j) + β \end{matrix}  (16)

We are now finally in a position where we can formulate how to compute D_{x_{1}x_{2}x_{3}}^{6}(i, j) in constant time. The cost is the minimum cost of case 1–4. The cost of case 1 is L_{x_{1}x_{2}}^{6}(i, j) + c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}) and the minimum cost of case 2–4 is F_{x_{1}x_{2}x_{3}}^{6}(i, j), so

D_{x_{1}x_{2}x_{3}}^{6}(i, j) = min \{ \begin{matrix}L_{x_{1}x_{2}}^{6}(i, j) + c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}) \\ F_{x_{1}x_{2}x_{3}}^{6}(i, j) \end{matrix}  (17)

To compute D^{6}(i, j) in constant time by expression (11) we must compute D_{x_{1}x_{2}x_{3}}^{6}(i, j) for each of the 32 possible remaining codons. To do this we must compute entry (i, j) in the 16 tables L_{x_{1}x_{2}}^{6} and entry (i, j) in the 64 tables F_{x_{1}x_{2}x_{3}}^{6}. The other three cases where the last codon alignment describes one internal gap (type 4, 5 and 7) are handled similarly. However, if the last codon alignment is of type 4 or 5, then only the first nucleotide x_{1} in the remaining codon depends on the deletion (or insertion) length. This limits the number of possible remaining codons to 16 and implies that only four tables are needed to keep track of min_{0<k<i} len_{x_{1}}^{t}(i, j, k) for t = 4, 5. Hence, to compute D^{t}(i, j) for t = 4, 5, 6, 7, we must compute 2·4+2·16+4·64 = 296 table entries in total.

Codon alignments with two internal gaps

We will describe how to compute D^{8}(i, j). The other three cases where the last codon alignment describes two internal gaps are handled similarly. The last codon alignment of type 8 describes three substitutions and two deletions. If the first deletion has length k′ and the second deletion has length k, then the last codon alignment is an alignment of a_{1}^{i′′}a_{2}^{i′′}a_{3}^{i′′} ⋯ a_{1}^{i′}a_{2}^{i′}a_{3}^{i′} ⋯ a_{1}^{i}a_{2}^{i}a_{3}^{i} and b_{1}^{j}b_{2}^{j}b_{3}^{j} where i′ = i − k and i′′ = i′ − k′. This is illustrated in figure 10.

[figure: a_{1}^{i′′}; a_{2}^{i′′}; a_{3}^{i′′}; a_{1}^{i′}; a_{2}^{i′}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; a_{3}^{i}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; (3i,3j); first deletion; second deletion]

Figure 10: The last codon alignment of type 8

We will compute D^{8}(i, j) in the same way as we computed D^{6}(i, j). We will minimize the cost over all possible combinations of deletion length and remaining codon of the second deletion. This reduces the problem to computing D_{x_{1}x_{2}x_{3}}^{8}(i, j), the cost under the assumption of a certain remaining codon of the second deletion, for each of the 32 possible remaining codons of the second deletion. We will compute D_{x_{1}x_{2}x_{3}}^{8}(i, j) similar to the method described by expression (15), (16) and (17). The method can be used with almost no modifications. All we essentially have to do is to use len_{x_{1}x_{2}}^{8}(i, j, 1) instead of len_{x_{1}x_{2}}^{6}(i, j, 1) in expression (16).

The cost len_{x_{1}x_{2}x_{3}}^{8}(i, j, k) is the part of the cost D^{8}(i, j) that does not depend on the codon a_{1}^{i}a_{2}^{i}a_{3}^{i} and the witness of the second deletion, under the assumption the second deletion has length k and remaining codon x_{1}x_{2}x_{3}. The cost depends on the order of the two deletions in the last codon alignment, so we introduce len_{x_{1}x_{2}x_{3}}^{8′}(i, j, k) and len_{x_{1}x_{2}x_{3}}^{8′′}(i, j, k) to denote the cost when the first deletion occurs before the second deletion and vice versa. We define len_{x_{1}x_{2}x_{3}}^{8}(i, j, k) as min{len_{x_{1}x_{2}x_{3}}^{8′}(i, j, k), len_{x_{1}x_{2}x_{3}}^{8′′}(i, j, k)}. We only have to compute len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1). In the following we will therefore examine the cost D^{8}(i, j) under the assumption that the second deletion has length one and remaining codon x_{1}x_{2}x_{3}. We split the examination depending on the order of the first and second deletion.

Figure 11 illustrates the evolution of the last codon alignment (of type 8) when the second deletion has length one and occurs after the first deletion. The nucleotides y_{1}, y_{2} and y_{3} are the status of the substitutions before the first deletion and the nucleotides x_{1}, x_{2} and x_{3} are the status of the substitutions before the second deletion. We can regard the substitution a_{1}^{i′′} → b_{1}^{j} as occurring in three steps: a^{i′′} → y_{1}, y_{1} → x_{1} and x_{1} → b_{1}^{j} where y_{1} ∈ {a_{1}^{i′′}, x_{1}} and x_{1} ∈ {a_{1}^{i′′}, b_{1}^{j}}. Similarly, we can regard the other two substitutions as occurring as a_{2}^{i′} → y_{2} → x_{2} → b_{2}^{j} and a_{3}^{i} → y_{3} → x_{3} → b_{3}^{j}.

[figure: y_{1}; a_{2}^{i′′}; a_{3}^{i′′}; a_{1}^{i′}; y_{2}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; y_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{1 . del}; y_{1}; y_{2}; a_{1}^{i′}; a_{1}^{i}; a_{2}^{i}; y_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; x_{1}; x_{2}; a_{1}^{i′}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{2 . del}; x_{1}; x_{2}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; b_{1}^{j}; b_{1}^{j}; b_{1}^{j}]

Figure 11: The first deletion occurs before the second deletion and the second deletion has length one.

The cost of the substitutions in the last codon alignment is

cost(subs) = c_{p}^{∗}(a_{1}^{i′′}a_{2}^{i′′}a_{3}^{i′′}, y_{1}a_{2}^{i′′}a_{3}^{i′′}) + c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, a_{1}^{i′}y_{2}a_{3}^{i′}) +

c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}y_{3}) + c_{p}^{∗}(y_{1}y_{2}a_{1}^{i′}, x_{1}x_{2}a_{1}^{i′}) +

c_{p}^{∗}(a_{1}^{i}a_{2}^{i}y_{3}, a_{1}^{i}a_{2}^{i}x_{3}) + c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}). (18)

The deletion length of the first deletion is k′ and the deletion length of the second deletion is one. Recall that i′ = i − 1 and i′′ = i′ − k′. Similar to expression (10) we can formulate the cost of the two deletions as

cost(del_{1}) = min \{ \begin{matrix}α + βk′ + c_{p}(a_{1}^{i′}y_{2}a_{3}^{i′}, y_{1}y_{2}a_{3}^{i′}) \\ 2α + βk′ + min_{0<l<k′} c_{p}(a_{1}^{i′−l}a_{2}^{i′−l}a_{3}^{i′−l}, y_{1}y_{2}a_{3}^{i′}) \\ α + βk′ + c_{p}(y_{1}a_{2}^{i′′}a_{3}^{i′′}, y_{1}y_{2}a_{3}^{i′}) \end{matrix} (19)

cost(del_{2}) = α + β + min \{ \begin{matrix}c_{p}(x_{1}x_{2}a_{3}^{i′}, x_{1}x_{2}x_{3}) \\ c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}) \end{matrix}  (20)

The cost of the remaining alignment is D(i′ − k′ − 1, j − 1), so the minimum cost under the assumption of x_{1}x_{2}x_{3} as the remaining codon of the second deletion is given by the sum cost(subs) + cost(del_{1}) + cost(del_{2}) + D(i′ − k′ − 1, j − 1) minimized over all possible combinations of y_{1}y_{2}y_{3} and k′. The cost len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1) is the part of this cost that does not depend on a_{1}^{i}a_{2}^{i}a_{3}^{i} or the witness of the second deletion. This part includes everything except the terms c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}y_{3}) + c_{p}^{∗}(a_{1}^{i}a_{2}^{i}y_{3}, a_{1}^{i}a_{2}^{i}x_{3}) of cost(subs) and min{c_{p}(x_{1}x_{2}a_{3}^{i′}, x_{1}x_{2}x_{3}), c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3})} of cost(del_{2}). It is easy to verify that the minimum of D(i′ − k′ − 1, j − 1) + c_{p}^{∗}(a_{1}^{i′′}a_{2}^{i′′}a_{3}^{i′′}, y_{1}a_{2}^{i′′}a_{3}^{i′′}) + cost(del_{2}) over the length k′ of the first deletion is D_{y_{1}y_{2}a_{3}^{i′}}^{4} (i′, j), so

len_{x_{1}x_{2}x_{3}}^{8′}(i, j, 1) = α + β + c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}) +

min_{y_{1}y_{2}}{c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, a_{1}^{i′}y_{2}a_{3}^{i′}) + D_{y_{1}y_{2}a_{3}^{i′}}^{4} (i′, j) + c_{p}^{∗}(y_{1}y_{2}a_{3}^{i′}, x_{1}x_{2}a_{3}^{i′})} (21)

where we minimize over y_{1} ∈ {a_{1}^{i′′}, x_{1}} and y_{2} ∈ {a_{2}^{i′}, x_{2}}. The cost len_{x_{1}x_{2}x_{3}}^{8′}(i, j, 1) is defined if x_{1}x_{2}x_{3} allows the second deletion to have length one, i.e. if x_{1} ∈ {a_{1}^{i′′}, b_{1}^{j}}, x_{2} ∈ {a_{2}^{i′}, b_{2}^{j}} and x_{3} ∈ {a_{3}^{i}, b_{3}^{j}}. The nucleotide a_{1}^{i′′} depends on the unknown length of the first deletion, so we must assume that it can be any of the four nucleotides.

Figure 12 illustrates the evolution of the last codon alignment when the second deletion has length one and occurs before the first deletion. The nucleotides z_{1}, x_{2} and x_{3} are the status of the substitutions before the second deletion and the nucleotides y_{1}, y_{2} and y_{3} are the status of the substitutions before the first deletion. The first nucleotide x_{1} in the remaining codon of the second deletion is just a_{1}^{i′}. A detailed description of the cost can be done as above. It would reveal that len_{x_{1}x_{2}x_{3}}^{8′′}(i, j, 1) includes everything except the cost of the substitution a_{3}^{i} → x_{3} → y_{3} and c_{p}(x_{1}x_{2}x_{3}, w_{1}w_{2}w_{3}) where w_{1}w_{2}w_{3} is the witness of the second deletion. It would also reveal the minimum cost of the remaining alignment, the substitution a_{1}^{i′′} → y_{1} and the first deletion is given by D_{y_{1}y_{2}y_{3}}^{4}(i′, j), so

len_{x_{1}x_{2}x_{3}}^{8′′}(i, j, 1) = α + β + c_{p}^{∗}(a_{1}^{i′}a_{2}^{i′}a_{3}^{i′}, x_{1}x_{2}a_{3}^{i′}) +

min_{y_{1}y_{2}y_{3}}{c_{p}^{∗}(a_{1}^{i′}x_{2}x_{3}, a_{1}^{i′}y_{2}y_{3}) + D_{y_{1}y_{2}y_{3}}^{4}(i − 1, j) + c_{p}^{∗}(y_{1}y_{2}y_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j})} (22)

where we minimize over y_{1} ∈ {z_{1}, b_{1}^{j}}, y_{2} ∈ {x_{2}, b_{2}^{j}} and y_{3} ∈ {x_{3}, b_{3}^{j}}. The nucleotide z_{1} depends on the unknown length of the first deletion, so we must assume that z_{1} can be any of the four nucleotides. The cost len_{x_{1}x_{2}x_{3}}^{8′′}(i, j, 1) is defined if x_{1}x_{2}x_{3} allows the second deletion to have length one, i.e. if x_{1} = a_{1}^{i′}, x_{2} ∈ {a_{2}^{i′}, b_{2}^{j}} and x_{3} ∈ {a_{3}^{i}, b_{3}^{j}}.

[figure: z_{1}; a_{2}^{i′′}; a_{3}^{i′′}; x_{1}; x_{2}; a_{3}^{i′}; a_{1}^{i}; a_{2}^{i}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{2 . del}; z_{1}; a_{2}^{i′′}; a_{3}^{i′′}; x_{1}; x_{2}; x_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; y_{1}; a_{2}^{i′′}; a_{3}^{i′′}; x_{1}; y_{2}; y_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{1 . del}; y_{1}; y_{2}; y_{3}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; ⟶^{subs}; b_{1}^{j}; b_{2}^{j}; b_{3}^{j}; b_{1}^{j}; b_{1}^{j}; b_{1}^{j}]

Figure 12: The second deletion occurs before the first deletion and the second deletion has length one.

We are finally in a position where we can formulate how to use the method from the previous section to compute D^{8}(i, j). The cost len_{x_{1}x_{2}x_{3}}^{8}(i, j, k) depends on x_{1}, x_{2} and x_{3}, so instead of 16 tables we need 64 tables to keep track of min_{0<k<i} len_{x_{1}x_{2}x_{3}}^{8}(i, j, k). We still need 64 tables to keep track of the cost under the assumption that the end-codon a_{1}^{i}a_{2}^{i}x_{3} is not the witness (of the second deletion). We compute table entry (i, j) as

L_{x_{1}x_{2}x_{3}}^{8}(i, j) = min \{ \begin{matrix}len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1) \\ L_{x_{1}x_{2}x_{3}}^{8}(i − 1, j) + β \end{matrix}  (23)

F_{x_{1}x_{2}x_{3}}^{8}(i, j) = min \{ \begin{matrix}L_{x_{1}x_{2}x_{3}}^{8}(i − 1, j) + β + α + c_{p}(a_{1}^{i−1}a_{2}^{i−1}a_{3}^{i−1}, x_{1}x_{2}x_{3}) \\ len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1) + c_{p}(x_{1}x_{2}a_{3}^{i−1}, x_{1}x_{2}x_{3}) \\ F_{x_{1}x_{2}x_{3}}^{8}(i − 1, j) + β \end{matrix} (24)

We compute D_{x_{1}x_{2}x_{3}}^{8}(i, j) by using the above tables, and we compute D^{8}(i, j) by minimizing over the 32 possible remaining codons of the second deletion.

D_{x_{1}x_{2}x_{3}}^{8}(i, j) = min \{ \begin{matrix}L_{x_{1}x_{2}x_{3}}^{8}(i, j) + c_{p}(a_{1}^{i}a_{2}^{i}x_{3}, x_{1}x_{2}x_{3}) \\ F_{x_{1}x_{2}x_{3}}^{8}(i, j) \end{matrix}  (25)
D^{8}(i, j) = min_{x_{1}x_{2}x_{3}}{c_{p}^{∗}(a_{1}^{i}a_{2}^{i}a_{3}^{i}, a_{1}^{i}a_{2}^{i}x_{3}) + D_{x_{1}x_{2}x_{3}}^{8}(i, j)}  (26)

This computes D^{8}(i, j) in constant time. The computation require us to compute entry (i, j) in 128 tables. The other three cases where the last codon alignment describes two internal gaps, type 9–11, are handled similarly, so to compute D^{t}(i, j) for t = 8, 9, 10, 11 we compute entry (i, j) in 4·128 = 512 tables. Finally, to compute D(i, j) by expression (5) we compute entry (i, j) in 1 + 11 + 296 + 512 = 820 tables. The 512 entries are explained in this section, the 296 entries was explained in the previous section, the 11 entries correspond to D^{t}(i, j) for t = 1, 2, …, 11 and the last entry corresponds to D(i, j).

5 Improvements and future work

The only real difference between the computation of D^{6}(i, j) and D^{8}(i, j) is between len_{x_{1}x_{2}}^{6}(i, j, 1) and len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1). The similarity stems from the fact that a codon alignment of type 6 and type 8 ends in the same way. By “end in the same way” we mean that the events described on the codon a_{1}^{i}a_{2}^{i}a_{3}^{i} are the same (see figure 5). A codon alignment of type 11 also ends in the same way as a codon alignment of type 6 or 8. The similarity between the computation of D^{t}(i, j) for t = 6, 8, 11 makes it possible to replace the six tables L_{x_{1}x_{2}x_{3}}^{t} and F_{x_{1}x_{2}x_{3}}^{t} for t = 6, 8, 11 by two tables L_{x_{1}x_{2}x_{3}}^{6,8,11} and F_{x_{1}x_{2}x_{3}}^{6,8,11} in which entry (i, j) equals the minimum of entry (i, j) in the tables being replaced. We compute L_{x_{1}x_{2}x_{3}}^{6,8,11}(i, j) and F_{x_{1}x_{2}x_{3}}^{6,8,11}(i, j) similar to expression (23) and (24), except that we use

len_{x_{1}x_{2}x_{3}}^{6,8,11}(i, j, 1) = min_{t=6,8,11} len_{x_{1}x_{2}x_{3}}^{t}(i, j, 1)  (27)

instead of len_{x_{1}x_{2}x_{3}}^{8}(i, j, 1). We use len_{x_{1}x_{2}x_{3}}^{6}(i, j, 1) defined as len_{x_{1}x_{2}}^{6}(i, j, 1)+ c_{p}^{∗}(x_{1}x_{2}x_{3}, b_{1}^{j}b_{2}^{j}b_{3}^{j}) to ensure that len_{x_{1}x_{2}x_{3}}^{t}(i, j, 1) for t = 6, 8, 11 describes the same part of the total cost. The cost D^{6,8,10}(i, j) is the minimum of D^{t}(i, j) over t = 6, 8, 10. To compute D^{6,8,11}(i, j) by expressions similar to expression (25) and (26) we compute entry (i, j) in 1+64+64 = 129 tables. To compute D^{t}(i, j) for t = 6, 8, 10 we compute entry (i, j) in 3+80+128+128 = 339 tables. The combined computations thus saves the computation of 210 table entries. Codon alignments of type 7, 9 and 10 also end in the same way, so we can also combine the computation of D^{t}(i, j) for t = 7, 9, 10. We can thus reduce the number of table entry updates in the computation of D(i, j) to 820 − 2 · (339 − 129) = 400.

We are working on implementing the alignment algorithm in order to compare it to the heuristic alignment algorithm described in [4]. The heuristic algorithm allows frame shifts, so an obvious extension of our exact algorithm would be to allow frame shifts, e.g. to allow insertion-deletions of arbitrary length. This however makes it difficult to split the evaluation of the alignment cost into small independent subproblems (codon alignments) of known size. Another interesting extension would be to annotate the DNA sequence with more information. For example, if the DNA sequence codes in more than one reading frame (overlapping reading frames) then the DNA sequence should be annotated with all the amino acid sequences encoded and the combined cost of a nucleotide event should summarize the cost of changes induced on all the amino acid sequences encoded by the DNA sequence. This extension also makes it difficult to split the evaluation of the alignment cost into small independent subproblems. To implement these extensions efficiently it might be fruitful to investigate reasonable restrictions of the cost functions.

References

[1] Arvestad, L. Aligning coding DNA in the presence of frame-shift errors. In Combinatorial Pattern Matching (1997), vol. 1264 of LNCS, pp. 180–190.

[2] Gotoh, O. An improved algorithm for matching biological sequences. Journal of Molecular Biology 162 (1981), 705–708.

[3] Hein, J. An algorithm combining DNA and protein alignment. Journal of Theoretical Biology 167 (1994), 169–174.

[4] Hein, J., and Støvlbæk, J. Genomic alignment. Journal of Molecular Evolution 38 (1994), 310–316.

[5] Hein, J., and Støvlbæk, J. An algorithm combining DNA and protein alignment. Methods in Enzymology 266 (1996), 402–418.

[6] Hua, Y., Jiang, T., and Wu, B. Aligning DNA sequences to minimize the change in protein. Manuscript, October 1997.

[7] Needleman, S. B., and Wunsch, C. D. A general method applicable to the search for similarities in the amino acid seqeunce of two proteins. Journal of Molecular Biology 48 (1970), 433–443.

[8] Peltola, H., Söderlund, H., and Ukkonen, E. Algorithms for the search of amino acid patterns in nucleic acid sequences. Nuclear Acids Research 14, 1 (1986), 99–107.

[9] Sankoff, D. Matching sequences under deletion/insertion constraints. In Proc. Nat. Acad. of Sci. U.S.A. (1972), vol. 69, pp. 4–6.

[10] Sellers, P. H. On the theory and computation of evolutionary distance. SIAM Journal of Applied Mathematics 26 (1974), 787–793.

[11] Wagner, R. A., and Fisher, M. J. The string to string correction problem. Journal of the ACM 21 (1974), 168–173.

[12] Zhang, Z., Pearson, W. R., and Miller, W. Aligning a DNA sequence with a protein sequence. In Proceedings of the First Annual International Conference on Computational Molecular Biology (1997), ACM, pp. 337–343.

Recent BRICS Report Series Publications

RS-98-3 Christian N. S. Pedersen, Rune B. Lyngsø, and Jotun Hein. Comparison of Coding DNA. January 1998. 20 pp.

RS-98-2 Olivier Danvy. An Extensional Characterization of Lambda-Lifting and Lambda-Dropping. January 1998.

RS-98-1 Olivier Danvy. A Simple Solution to Type Specialization (Extended Abstract). January 1998. 7 pp.

RS-97-53 Olivier Danvy. Online Type-Directed Partial Evaluation. December 1997. 31 pp. Extended version of an article to appear in Third Fuji International Symposium on Functional and Logic Programming, FLOPS ’98 Proceedings (Kyoto, Japan, April 2– 4, 1998).

RS-97-52 Paola Quaglia. On the Finitary Characterization of π-Congruences. December 1997. 59 pp.

RS-97-51 James McKinna and Robert Pollack. Some Lambda Calculus and Type Theory Formalized. December 1997. 43 pp.

RS-97-50 Ivan B. Damgård and Birgit Pfitzmann. Sequential Iteration of Interactive Arguments and an Efficient Zero-Knowledge Argument for NP. December 1997. 19 pp.

RS-97-49 Peter D. Mosses. CASL for ASF+SDF Users. December 1997. 22 pp. Appears in ASF+SDF’97, Proceedings of the 2nd International Workshop on the Theory and Practice of Algebraic Specifications, Electronic Workshops in Computing, http://www.springer.co.uk/ewic/workshops/ASFSDF97. Springer-Verlag, 1997.

RS-97-48 Peter D. Mosses. CoFI: The Common Framework Initiative for Algebraic Specification and Development. December 1997. 24 pp. Appears in Bidoit and Dauchet, editors, Theory and Practice of Software Development. 7th International Joint Conference CAAP/FASE, TAPSOFT ’97 Proceedings, LNCS 1214, 1997, pages 115–137.

RS-97-47 Anders B. Sandholm and Michael I. Schwartzbach. Distributed Safety Controllers for Web Services. December 1997. 20 pp. To appear in European Theory and Practice of Software. 1st Joint Conference FoSSaCS/FASE/ESOP/CC/TACAS, ETAPS ’97 Proceedings, LNCS, 1998.