Trait-based phylogenetic inference

Introduction

Another way of thinking about phylogenetic analysis is to focus on the individual letters in a sequence, or on individual traits (or characters) in an organism. Using such an approach allows phylogenetic analysis even if there is no easy way to define a distance between the organisms' features. For example, wings, arms, and flippers are corresponding traits of organisms, but there's no obvious way to define how similar or different they are from each other. Before DNA sequencing became widely used, phylogenetic analyses primarily relied on this trait-based approach. This approach can also be used for DNA or protein sequences, by treating each individual letter as a different trait.

IMPORTANT: The Pevzner and Compeau textbook refers to this analysis as "character-based tree reconstruction". I have chosen to use the word trait in order to avoid the confusion between character, as in a characteristic of an organism; and character, as in a letter in a biological sequence.

We can try to formalize the problem of phylogenetic inference:

Phylogenetic inference: Given a set of organisms (or sequences) and associated traits (nucleotides or amino-acids in the case of sequences), identify the phylogenetic tree that best captures the shared evolutionary history of the organisms.

When reading this formulation, you should feel a bit uneasy. While the input is clear, as is the tree structure we are trying to develop, the objective function is certainly non-rigorous: best captures the shared evolutionary history.

For now, let's try to focus on one formulation of this objective—parsimony. By parsimony, we mean that we want to define the tree that implies the fewest changes between the traits (e.g., an arm evolving into a wing, or an A evolving into a C). We can, thus, add a bit more rigor to the problem described above:

Phylogenetic inference: Given a set of organisms (or sequences) and associated traits (nucleotides or amino-acids in the case of sequences), identify the phylogenetic tree that implies the fewest changes between traits.

If you think about this a little bit, you will note that there's no easy way to measure the number of changes in a tree, if all the traits are represented at the leaves. To address this problem, we redefine the problem once again:

Phylogenetic inference: Given a set of organisms (or sequences) and associated traits (nucleotides or amino-acids in the case of sequences), identify the phylogenetic tree and assignment of traits at internal nodes in the tree, that imply the fewest changes between traits.

Note that we've made the problem a bit more difficult by asking not just for a tree, but also for an assignment of traits at internal nodes. You'll see soon that this does not actually increase the computational complexity of solving the problem.

Before we proceed further, it is important to note that the search space is extremely large. There are $(2n -3)! \over 2^{n-2} (n-2)!$ possible binary trees with n leaves, and for each of them you have an exponential number of possible assignments of traits to internal nodes. If we use DNA letters as traits, given there are $O(n – 1)$ internal nodes in a tree with n leaves, we have $O(4^{n-1})$ different assignments. Thus, the total number of possible solutions is the product of these two numbers:

$4^{n-1} {(2n -3)! \over 2^{n-2} (n-2)!} = 2^n {(2n-3)! \over (n-2)!}$

By now you should be familiar with some strategies for searching for a solution in such a large space. A typical paradigm alternates between two steps:

1. computing the score of a potential solution

2. selecting a new solution that is likely to have a better score

You have seen examples in the Gibbs sampler and k-means clustering. In the context of phylogenetic analysis, the first step is called small parsimony, while the overall goal of reconstructing a tree with the smallest possible parsimony score is referred to as large parsimony.

Solving small parsimony

To briefly recap, we are now focusing on the first step of an algorithm for searching the large space of trees in order to find the phylogenetic tree that best explains the evolutionary history of a group of modern-day organisms. To simplify things, we will assume that the traits we are analyzing are letters in the DNA sequence of a gene that is found in all the organisms of interest. Thus, we can view each trait as an individual character with four possible values: A, C, G, and T. Given a tree, with characters at its leaves, we want to determine the characters assigned to internal nodes such that the parsimony score of the tree is minimized. The parsimony score is simply the number of edges in the tree that have different letters their opposite ends, indicating that an evolutionary change occurred. An example is shown below:

Parsimony score. Edges that represent evolutionary changes are highlighted with the number 1. The total parsimony score of this tree is 4.

To formalize the problem we are trying to solve:

Small parsimony problem: Given a tree and single character labels at its leaves, determine the assignment of characters at internal nodes such that the parsimony score of the tree is minimized.

Before we proceed with a discussion of the algorithm, I would like to mention that usually, the leaves contain more than one character, e.g., the entire sequence of a gene. When computing the parsimony, however, each individual character is independent of the others, thus it is sufficient to solve the problem for a single character. If provided a sequence, the answer is simply the concatenation of the solutions for each separate character, and the total parsimony score is simply the sum of the parsimony scores for each character.

We will attempt to solve the small parsimony problem using dynamic programming. As we did in the case of alignment, during the first stage of the algorithm we will figure out the optimal score, then we will use backtracking to figure out the actual assignment of letters to each internal node of the tree.

To get a better handle on how the algorithm proceeds, let us focus on just one node and its two children. Since we do not know what letter will be assigned to any of the nodes, we will have to keep track of all possible solutions. Specifically, at each node v, we will store an array of length 4 (the number of possible choices) that contains the parsimony score assuming the node is assigned the corresponding letter. More formally, we'll define $s_k(v)$ to be the score of the tree rooted at v if character k is assigned at v. For example, $s_A(v)$ represents the score of the tree rooted at v that has A at the root. See the figure below for an example.

The score at node v. $s_A(v) = 1$ since one edge entails a change from A at the root to C at the leaf, while $s_T(v) = 2$

Thus, at each node in the tree we will store the $s_k(V)$ array:

$s_k(v)$ array at each node in the tree

Since we know the characters at the leaves, and we do not want to allow the algorithm to select any other choice, we will set $s_k(v) = 0$ for the value of k equal to the character assigned to the leaf, and $s_k(v) = \infty$ for all other characters. An example is shown below for a tree that has leaf labels A and C.

Setting the scores at leaves.

All that remains to figure out is the recurrence equation. Here, we'll assume the values are computed for the children, and we want to compute them for the parent. Let us take the example in the figure above. Let us assume we want to fill in the value for $s_A(\text{root})$ . Of course the value depends on the actual choices for the two children, $\text{left}(\text{root})$ and $\text{right}(\text{root})$. If we select A at $\text{left}(\text{root})$ and C at $\text{right}(\text{root})$, we see that $s_A(\text{root}) = 1$ because the right edge (from A to C) implies one evolutionary event. What if we select A at $\text{right}(\text{root})$? Now the right edge has a cost of 0 (from A to A), however the total cost on the right side is $\infty$ because $s_A(\text{right}(v)) = \infty$ and $0 + \infty = \infty$ . Since we are trying to minimize parsimony, the first choice was better. In general:

$s_k(v) = \min_i\{s_i(\text{left}(v)) + \delta_{i,k}\} + min_i\{ s_i(\text{right}(v)) +\delta_{i,k}\}$

where $\delta_{i, k}$ is the cost of going from letter i to letter k: $\delta_{i,k} = (i = k ? 0 : 1)$ . Alternatively written:

$\delta_{i,k} = \left \{\begin{matrix}0 \text{ if } i = k \\ 1 \text{ if } i \ne k \end{matrix} \right .$

In other words, we are minimizing the left and right children of the current node separately, and picking the assignment at each child that minimizes the sum of the transition from parent to child, plus the corresponding value in the child.

Using this recurrence, the tree can be filled in using a post-order traversal. The parsimony score of the entire tree is found at the location $\text{argmin}_k s_k(\text{root})$ , or the smallest value in the s array at the root of the tree.

To figure out the actual letters, we backtrack down the tree, starting with the smallest value in the array located at root. Since the left and right children are independent, we can recurse on the two sub-trees. At each step in the backtracking process, we know the letter assigned to the parent of the node, and we need to determine the letter assigned to the current node. That letter will be the one that led to the parent's letter to be chosen. Specifically, assume we are at node v in the tree and the parent is assigned letter k. We will select letter j at node v as $\text{argmin}_i (s_i(v) + \delta_{i,k})$ . We can then recurse on the two children of v and proceed until we reach the leaves.

As in the case of the alignment algorithm, ties can be broken arbitrarily, yielding one of possibly many solutions that all have the same minimal parsimony score.

The algorithm described in this section is due to David Sankoff and is commonly known as Sankoff's algorithm.

Sankoff's algorithm walk-through

To clarify the algorithm described above, we will follow through a small example comprising four nodes of a larger phylogenetic tree. We assume the top node is the root of this tree, and that some of the nodes have already been filled in during the post-order traversal computation of the parsimony scores. To make it easier to follow the algorithm, we gave each node a name (a – e). A figure showing this partial tree is shown below. The tree can also be represented in the Newick format as:

((c:[3,2,1,2], d:[1,3,2,3])b:[,,,], e:[1,3,2,2])a:[,,,]

The array listed as an annotation for each node represents the score for the letters A, C, G, T, in this specific order. That is $s_A(c) = 3, s_G(d) = 2$, etc., using the notation we used earlier.

Since the values stored at both children of node b are already computed, we have enough information to compute the values that are stored at node b. To compute the value that will be stored for A, we identify the lowest-cost choice on the left and right independently and then add up their values. On the left side (going from b to c) we have the following 4 choices, from which we take the minimum:

$\min\left \{ \begin{matrix} s_A(c) + 0 = 3 \\ s_C(c)+1 =3 \\ s_G(c) + 1 = 2 \\ s_T(c) + 1 = 3 \end{matrix}\right . = 2$

While on the right (going from b to d) we have:

$\min\left \{ \begin{matrix} s_A(d) + 0 = 1 \\ s_C(d)+1 =4 \\ s_G(d) + 1 = 3 \\ s_T(d) + 1 = 4 \end{matrix}\right . = 1$

The value of $s_A(b)$ is the sum of the two values, i.e., 3.

Similarly:

$s_C(b) = \min\left \{ \begin{matrix} s_A(c) + 1 = 4 \\ s_C(c)+0 =2 \\ s_G(c) + 1 = 2 \\ s_T(c) + 1 = 3 \end{matrix}\right . + \min\left \{ \begin{matrix} s_A(d) + 1 = 2 \\ s_C(d)+0 =3 \\ s_G(d) + 1 = 3 \\ s_T(d) + 1 = 4 \end{matrix}\right . = 4$

And so on, yielding the tree:

Now both children of the root (node a) have their values computed, and we can use the same process to fill in the values at the root, yielding the tree:

We can now start the backtracking process. We note that the lowest score at the root (node a) corresponds to the letter A, thus we select that choice, shown as a shaded box in the figure below:

We now have to figure out what letter assignment at nodes b and e leads to the score $s_A(a)$. On the right, we see that the best choice is to select A, since it will yield a total cost of 1. On the left, the best choice is again A, with a total cost of 3, which summed up with the score from the right yields the value 4 represented at node a. Thus, we can make the assignments of letter A for both nodes b and e.

We can now proceed to figure out the letters stored at nodes c and d. At node c, the best choice is the letter G, which yields a total cost of 2 (1 from $s_G(c)$ ) and another 1 from the mutation from A to G). While the cost of going from A to A is 0, $s_A(c) = 3$, thus choosing A at node c would be sub-optimal. At node d, A is the best choice, yielding a total score of 1. This results in the assignments:

In the Newick format, this tree can be now represented as:

((c:G, d:A)b:A, e:A)a:A

I encourage you to draw arbitrary trees and assign arbitrary letters at the leaves, then try to run through the algorithm a couple of times to make sure you understand how it works, and to get the practice necessary to ensure you will be able to solve problems assigned on quizzes and exams.

From parsimony to likelihood

When we discussed sequence alignment, we made a transition from an edit distance formulation, where we tried to minimize the number of edits between two sequences, to an alignment score, that took into account substitutions scores represented in a matrix such as BLOSUM. Both variants could be solved with the same algorithm by simply switching from the minimum to the maximum between a number of choices. The same can be done in the case of phylogenetic tree inference. If we are given a score that defines how likely it is for a character to change into another during an evolutionary event, we can modify the recurrence described above:

$s_k(v) = \max_i \{s_i(\text{left}(v)) + \delta_{i,k}\} + \max_i\{s_i(\text{right}(v))+\delta_{i,k}\}$

In this case the score $\delta_{i, k}$ will be defined by a matrix such as:

	A	C	G	T
A	8	2	3	2
C	2	8	2	3
G	3	2	8	2
T	2	3	2	8

Here the score represents the likelihood a particular letter would evolve into a different one. You see a higher likelihood that a letter stays the same (the numbers on the diagonal being higher) as well as a higher likelihood of a transition (mutations between As and Gs, and between Cs and Ts, which are similar molecules) than a transversion (mutations between As and Ts or Cs; or between Gs and Ts or Cs). As and Gs are chemically similar molecules called purines, while Cs and Ts are molecules called pyrimidines.

Note that the matrix doesn't have to be symmetrical as shown above. The tree has a natural orientation, with evolution occurring from parent to child, thus, we can account for uni-directional evolutionary events, such as the mutation of Cs into Ts by deamination (discussed in the early chapter on patterns in DNA).

	A	C	G	T
A	8	2	3	2
C	2	8	2	5
G	4	2	8	2
T	2	3	2	8

In the matrix above, the values for C-T and G-A (in bold and highlight) are higher than the values for T-C and A-G, accounting for this phenomenon. Sankoff's algorithm can be used without modification, but it's implementation has to be done carefully to ensure the correct directionality is respected.

The scores represented by $\delta_{i, k}$ can also be the actual likelihoods that the two letters have evolved into each other, and can take into account the length of the branch, i.e., the evolutionary time that has passed between two nodes in the tree, which will impact the likelihood that a letter mutates into another. That is, we can define $\delta_{i, k} (t)$ where t is the time that "elapsed" along an edge of the tree. In this case the phylogenetic inference is also referred to as maximum likelihood, as we are trying to find the tree that maximizes the likelihood that the organisms at the leaves are related to each other as described by the tree.

Again, Sankoff's algorithm can be directly applied to solve this problem, except that now the information from the children will be aggregated by multiplication at each node (unlike scores, probabilities get multiplied). In practice, however, most algorithms dealing with probabilities work in logarithmic space to avoid under-flows (after just a few multiplications the numbers become infinitesimally small), thus reverting to an additive formulation for Sankoff's algorithm.

Solving large parsimony

While Sankoff's algorithm provides a way to evaluate (in polynomial time) the quality of a given tree, finding the tree that maximizes the quality is NP-hard. Most algorithms for constructing phylogenetic trees rely on heuristic search through the space of possible trees. A discussion of these algorithms is beyond the scope of the material presented here.

Exercises

1. Using Sankoff's algorithm, calculate the maximum score for the phylogenetic tree shown below, as well as the characters at each internal node. Note that the scoring matrix is not symmetric—the score of evolving A into T is different from the score of evolving T into A. In the scoring matrix, assume that the row represents the parent and the column the child in the tree. For example, if a node is assigned A and it's child C, the score along the branch is 3. If the parent is C and the child is A, the score is only 2. At each internal node in the tree please provide both the four scores corresponding to each individual nucleotide being present at that node, and the inferred ancestral sequence at that node (by circling that letter at the node). If multiple ancestral sequences are possible at each node, please pick one of them (e.g. the answer A/T is not acceptable)

	A	C	G	T
A	5	3	-1	0
C	2	5	0	-2
G	0	1	5	2
T	-1	-1	2	5