Algorithms for phylogenetics 0b: Trees

Versión en español

Another basic component of phylogenetic analysis' programs are the trees. Trees can be the result of an analysis (as in cladistic analysis' programs) or part of the data used (as in biogeography analysis' programs).

Trees are one of the most well known data structures of computer science. In computers science a tree is a collection of elements (nodes), nodes are related by a "kindship" relation that imposes a hierarchy on the nodes. Kinship is a pairwise relationship: one node is the ancestor and the other the descendant. A node can have an indefinite number of descendants, but, at maximum, only one ancestor.

In general, in phylogenetics, the term node is used only for elements with descendants, and terminal of leafs to elements without descendants (in computer science the used terms are internal node and terminal node, or leaf node. Here I will use the terminology from computer science, and when using node I refer to any element of the tree).

A trees has at maximum a single root node. A root node is an internal node without ancestor. Different from most abstract trees of computer science, the only elements with have the user data are the terminals (from the data matrix). Data stored by internal nodes is inferred through algoritmic means.

There are many ways to represent a tree. For me, the most natural one is using structures and pointers. Here is the basic node

typedef struct tagPHYLONODE PHYLONODE;

struct tagPHYLONODE {
PHYLONODE* anc;
PHYLONODE* left;
PHYLONODE* right;
STATES* recons;
};

This structure is a node from a binary tree, that is, each internal node have two descendants, the pointers left and right, note that they are pointers to PHYLONODE elements. The pointer to ancestor is anc. In this structure, the root node has anc as NULL (i.e. points to no elements), and terminal nodes has left and right as NULL. The array recons store the data (as bitfields).

I usually store the information from a terminal (extracted from data matrix) into an independent structure, for example

typedef struct tagTERMINAL TERMINAL;
struct tagTERMINAL {
char name [32];
STATES* chars;
};

And include a pointer to TERMINAL as part of PHYLONODE

struct tagPHYLONODE {
... //Otros campos aquí
TERMINAL* term;
};

Internal nodes has term as NULL.

A nice property of trees, is that subtrees have the same properties of a tree, then recursive algorithms a very natural way to work with trees. Nevertheless, I always prefer to have an independent structure to store the whole tree

typedef struct tagPHYLOTREE PHYLOTREE;
struct tagPHYLOTREE {
unsigned numNodes;
PHYLONODE* root;
PHYLONODE* nodeArray;
DATAMATRIX* matrix;
};

This structure store a pointer to the root node (root), an array to the nodes of the tree (nodeArray), a pointer to a data matrix (matix) and the number of nodes.

This approach has the advantage of having several independent sets of trees, each tree with his on characterisitcs (id, name, data matrix), that are no properties of any node, but of the whole set of nodes (A typical case: in biogeography, several cladograms form different organism are used).

Different from structures it is possible to store trees as coordinate arrays. This option is used in [1] for the example code, and it is also the way in which TNT macro language uses trees (the language does not support structures). The declarations of arrays is simple

int* anc;
int* left;
int* right;
STATES** chars;

In this case, all arrays are dynamic (they are declared as pointers, but as soon as they are assigned, they can be treated as arrays), and instead of have a pointer, they store the number that is the index of the array. As it is a coordinate array, every element with the same index is the same node. For example anc [17] = 20, assigns the node 20 as the ancestor of node 17. In this case, left [20] or right [20] must be equal to 17. chars [20] has the sates assigned to node 20. If you look at the tree entry in wikipedia (specially the entry on heaps) you will find that is the way used to keep trees.

For me, the bad side of that approximation is that requires an strict bookkeeping (take note that in TNT macro language, TNT automatically updates the information, so that is not a problem!). This is my own experience, I worked with tree reconciliation, and work with someone that does not know working with pointers. So I try a coordinate array approach. But we not know how many nodes would be on the tree after reconciliation, then the code turns messy and extremely difficult to update.

That experience does not show if some alternative is best than the other, just that sometimes the nature of the problem might constraint the way to resolve it, and that some styles of programming are better than others for some programmers ;). In this series I always will use structures.

Just to have some code, and to mark the recursivity of trees, here are a function to store a tree in a parenthetical notation (as in Hennig86, NONA and TNT)

#include < stdio.h >

void PrintTree (FILE* out, PHYLOTREE* tree) {
fprintf (out, "tread \'Simple tree output\'\n");
PrintNode (out, tree - > root)
fprintf (out, ";\np/;");
}

void PrintNode (FILE* out, PHYLONODE* node) {

if (node - > term != NULL) {
fprinf ("%s ", node - > term - > name);
return;
}
fprintf (out, "(");
PrintNode (out, node - > left);
PrintNode (out, node- > right);
fprintf (out, ") ");
}

This functions used standard IO operations. The function PrintTree puts the tree header for Hennig86/NONA/TNT trees. And call PrintNode on the root node. PrintNode checks if the node is a terminal, if this is true, it writes the terminal name, if not, the node is an internal one, so it prints the parenthesis and call recursively PrintNode on each one of its descendants.

I think that some time ago Mike Keesey write a post about trees under a object oriented paradigm. I'm unable to find the post, this one is the most closely alternative :P.

References
[1] Goloboff, P. A. 1997. Self weighted optimization: tree searches and character state reconstruction under implied transformation costs. Cladistics 13: 225-245. doi: 10.1111/j.1096-0031.1997.tb00317.x

Previous post on Algorithms for phylogenetics
0a. Bitfields

Algorithms for phylogenetics 0a: Bitfields

Versión en español

Intro

I want to start a series of post about the algorithms used in phylogenetic analyses. I feel a bit disappointed with the entry on computational phylogenetics (and related subjects) in wikipedia, that are somewhat biased towards model-based methods and "bioinformatics", with a poor representation of algorithms for phylogenetics.

There are two outstanding examples of algorithms for phylogenetics on the web. The book of Wheeler et al. [1], and, in a similar vein, a manuscript by DeLaet [2]. But their presentation of the algorithms is somewhat general, that is be nice for educational purpouse, but in some cases, far away from the actual implementation in computer software.

As the presentation of the algorithms in [1] and [2] is excellent, I want to focus on the implementation. The idea is to go from simpler to more complex algorithms (as far as I can comprehend them xD).

I want to acknowledge Pablo Goloboff. He always have the time to discuss with me about algorithms (and any topic on methodological phylogenetics) :D. Of course, any error in my implementations are my own errors!

As this series progress, I will post the source code on google source or sourceforge, if someone wants to check it, requires an knowledge of C programming language. Thanks to Rubén who showme some tips with the HTML :).

Bitfields

Characters are a basic component of phylogenetic analysis. For simplicity assume that we only have a single character in each terminal. Each character can be stored into a variable, assigning to them 0, 1, 2... as indicated by the character state. This idea is implicit in the original description of Wagner's "groundpland divergence" [3], and in Farris' formalization [4, see 5]. Then operations between characters would be operations of ranges. But it is more simpler to think about characters as a set of states [4, 6]. Then operations between characters would be set operations. This set has a peculiar property: they are perfectly boolean, that is, a taxon has or has not a particular state. This can be translated in two positive consequences: (i) character states can be stored as a bitfield: (ii) character operations are bifield operations.

Computers store numbers as binary numbers. A bitfields use this characteristic to represent a set of positioned bits. This example shows it better:

State    binary representation       "Number" (human)
0        0001                        1
1        0010                        2
2        0100                        4
3        1000                        8
0, 1     0011                        3
1, 3     0101                        5

As can be seen, each state has his on bit position, and the states combination is just the union of its states. Many programming languages include operations to work directly on bits (not, and, or and xor), then the translation from set operations to bit operations is direct.

In this series I will use char to store character states. In C, char has 8 bits, then only 8 states can be stored. I define the bitfield with a typedef, so changing the number of bits on the bitfield (and then the number of states) is easy.

#define BITS_ON_STATE 8
#define ALL_BITS_ON 255
typedef char STATES;

This small example reads the caracters of a taxon, from a file (in), and store it into an array of characters (chars). The number of characters is nc, and the functions getc (in stdio.h) reads a single character from the file, SkipSpaces ignores the blank characters, and isdigit (ctype.h) detects if the character is a number.

for (j = 0; j < nc; ++ j) {
error = SkipSpaces (in);
if (error != NO_ERROR) return error;
c = getc (in);
if (isdigit (c))
  chars [j] = 1 < < (c - '0');
else if ((c == '?') || (c == '-'))
  chars [j] = ALL_BITS_ON;
else return BAD_FORMAT;
}

The operation chars [j] 1 < < (c – '0') substract the ASCII value of '0' (30) from the ASCII value stored in c (a number, so it is form 30 to 39). The resulting value is used to shift the bits of 1. For example if a '0' is readed the operations are:

chars [j] = 1 < < (30 – 30)
chars [j] = 1 < < 0
     [0000 0001 < < 0 = 0000 0001]
chars [j] = 1

If a 5 is readed

chars [j] = 1 < < (35 – 30)
chars [j] = 1 < < 5
      [0000 0001 < < 5 = 0010 0000]
chars [j] = 32

The 1 is shifted five bits to the left.

Note that inapplicables and unknowns are coded with all bits on.

This work nicely with non additive characters. For additive characters, it is possible to “recode it” as several binary characters [7][8], then it is possible to use the character as independent binary characters.

There are more sophisticated usages of bitfields, packing several characters in a single variable (for example, eight 4-bit characters can be stored into a singel 32-bits variable) [8][9]. This packing allows an speed up during searches, because several characters can be examined simultaneously.

References
[1] Wheeler, W. et al. 2006. Dynamic homology and phylogenetic systematics: a unified approach using POY. American Mus. of Natural History, published in cooperation with NASA. online: http://research.amnh.org/scicomp/pdfs/wheeler/Wheeler_etal2006b.pdf
[2] DeLaet, J. 2005. Pseudocode for some tree search algorithms in phylogenetics. Manuscript online: http://www.plantsystematics.org/publications/jdelaet./algora.pdf
[3] Wagner, W.H. 1961. Problems in the classification of ferns. In: Recent advances in botany. Toronto Univ. Press. pp. 841-844.
[4] Farris, J.S. 1970. Methods for computing Wagner trees. Syst. Zool. 19: 83-92.
[5] Goloboff, P.A. et al. 2006. Continuos characters analyzed as such. Cladistics 22: 589-601. doi: 10.1111/j.1096-0031.2006.00122.x
[6] Fitch, W.M. 1971. Toward defining the course of evolution: minimum change for a specific tree topology. Syst. Zool. 20: 406-416.
[7] Farris, J.S. et al. 1970. A numerical approach to phylogenetic systematics. Syst. Zool. 19: 172-189.
[8] Moilanen, A. 1999. Searching for most parsimonious trees with simulated evolutionary optimization. Cladistics 15: 39-50. doi: 10.1111/j.1096-0031.1999.tb00393.x
[9] Goloboff, P.A. 2002. Optimization of polytomies: state and parallel set operations. Mol Phyl. Evol. 22: 269-275. doi: 10.1006/mpev.2001.1049

Core coding

En estos días, algunas de las ideas que tengo de mi proyecto, no están saliendo como yo quería, así que me he dedicado a la programación de algunos detalles computacionales... como el asunto, aunque muy relacionado con computadoras, pero poco con cladística (aparte claro, de ser parte de mi proy xD).. lo deje en mi blog personal :P.

English version:

How can Java (and several other languages) programmers live without pointers?

C--

Bueno, el post esta más relacionado con la programación que con la biogeografía... pero si les interesa, ahí puse algunas de mis experiencias recientes con el C++...

Seguir leyendo

A brief sketch of a program for phylogenetics

The last week is a nice week for me, as I received a grant to begin a graduate scholarship :) As my main objective is coupled with the development of a program, and my own interest in phylogenetics, and in a lesser extent, in software engineering, i will try to show here how I think a phyligenetic-related software would be implemented.

Do not take it for granted! The software projects in which i am involved to date, are smaller and short-term, and then I usually make some (ugly) shortcuts for immediate solutions but with low flexibility to modifications. Here my sketch take into account the possibility of further refinements of each part, so it is far more modular than my small projects.
Every program is composed with two parts. The first, the more messier and the most bored part to write, is the user interface (UI). The second, is the program itself (the core). In small projects you could mix both parts without much pain, but as the program growths leaving both parts mixed makes the program difficult to maintain: you would rewrite several code parts to introduce new input formats.

Ideally you would try to make the core to be totally independent from the UI, but usually it is not possible (unless your program do noting). An effort to construct a well designed communication between the UI and the core, is always well payed in the future.

I start my design with the core, and in every moment thinking how the core would communicate with a black-box UI. In phylogenetics you have two major parts: the matrix, and the tree. Both are composed of more smaller elements. The matrix is a list of labels, usually with some coupled data (taxon-character matrix, a taxon-distribution matrix). The tree is a topology in which each vertex is associated with the matrix labels.

You could think that both the tree and the matrix are the same, after all, the tree is only an specific way to sort the terminals in the matrix (it seems that TreeFitter[1] works in that way). But there is a difference, the terminals in the matrix are unique, but you could have many tress associated with the same matrix, terminal 'a' are independent in two trees, but they could point to the same element in a matrix.
Apart from terminals, that have a fixed content, internal nodes of a tree are coupled with the same sort of data as terminals, but this is a variable content, dependent on the tree topology. So we could implement a single fixed matrix, and one (or more) variable matrix, when the tree is unused, the variable matrix can be discarded.

To communicate the core with the UI, at this moment, I prefer string streams. They could be find in standard library (although with some little spelling differences), in that way the UI always translate between formats to send to the core specific strings, for example the tree in parenthesis format, so the tree construction would not interact with files, or console input, for example.

In that way, the UI is a central of translation from console, file, command dialogs, etc., that can be implemented for different systems (e.g. Win, Linux, command line), but the core remains the same.

[1] Ronquist, F. TreeFitter, program and documentation, available on line at: http://www.ebc.uu.se/systzoo/research/treefitter/treefitter.html