sábado, mayo 16, 2009

Algorithms for phylogenetics 0b: Trees


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

viernes, mayo 15, 2009

Algoritmos en filogenia 0b: Árboles


Otro de los componentes básicos de un programa de análisis filogenético son los árboles. Los árboles pueden ser el resultado del análisis (como en los programas de análisis cladístico), o pueden ser parte de los datos usados (como en los programas de análisis biogeográfico).

Los árboles, son además, una de las estructuras de datos más conocidas y usadas en programación. En programación un árbol es una colección de elementos (nodos), los nodos están relacionados entre sí por un "parentesco" que les impone una jerarquía. El parentesco es una relación de pares: un nodo es el ancestro y el otro es el descendiente. Un nodo puede tener un número indefinido de descendientes, más como máximo solo puede tener un ancestro.

En general, en filogenía se llama nodo, únicamente a los elementos que poseen descendientes, y hojas o terminales, a los elementos sin descendientes (en ciencias de computadores se utiliza más nodo interno y nodo terminal o nodo hoja, aquí yo sigo la notación de ciencias de computadores, y utilizo nodos cuando me refiero a todos los elementos del árbol).

Un árbol tiene como máximo un nodo raíz, que es un nodo interno sin ancestro. A diferencia de los árboles más abstractos de ciencias de computadores, los únicos elementos con datos incluidos por el usuario son los terminales (que entran vía la matriz de datos). Los datos guardados por los nodos internos son inferidos mediante algún algoritmo.

Existen varias formas de representar un árbol. Para mi, la más natural es usando estructuras y apuntadores. He aquí un nodo básico
typedef struct tagPHYLONODE PHYLONODE;

struct tagPHYLONODE {
PHYLONODE* anc;
PHYLONODE* left;
PHYLONODE* right;
STATES* recons;
};
Esta estructura es el nodo de un árbol binario, es decir, que cada nodo interno tiene dos descendientes, los apuntadores left, y right, que como pueden notar apuntan a elementos del tipo PHYLONODE. anc es el puntero a el ancestro. En esta estructura el nodo raíz tiene a anc como NULL (es decir que no apunta a ningún elemento), y los nodos terminales tienen left y right como NULL. recons es un arreglo con los caracteres (en forma de campo de bits).

Yo suelo guardar la información del terminal (extraída de la matriz de datos) en una estructura aparte, por ejemplo
typedef struct tagTERMINAL TERMINAL;
struct tagTERMINAL {
char name [32];
STATES* chars;
};
E incluyo un apuntador a TERMINAL como parte de PHYLONODE
struct tagPHYLONODE {
... //Otros campos aquí
TERMINAL* term;
};
Si es un nodo interno, term es NULL.

Una de las ventajas de los árboles es que un subárbol tiene las mismas propiedades de un árbol, por lo que los algoritmos recursivos son muy naturales al trabajarlos con árboles. Aún así, yo prefiero tener siempre una estructura que guarda el árbol
typedef struct tagPHYLOTREE PHYLOTREE;
struct tagPHYLOTREE {
unsigned numNodes;
PHYLONODE* root;
PHYLONODE* nodeArray;
DATAMATRIX* matrix;
};
En esta estructura se guarda un puntero al nodo raíz (root), un array con los nodos del árbol (nodeArray), un apuntador a la matriz de datos (matrix) y el número de nodos.

Esta aproximación tiene la ventaja de que uno bien puede tener varios conjuntos de árboles de forma independiente, por ejemplo, cada árbol con su propio id, o nombre, que no es una característica del nodo raíz, sino del conjunto completo de los nodos.

Aparte de las estructuras es posible guardar los árboles como arrays coordinados. Esta opción, por ejemplo, es usada en [1] para el código que aparece en los ejemplos, y también es la manera de manejar los árboles en el lenguaje de macros de TNT (que no maneja estructuras). La declaración de los arreglos coordinados es más sencilla
int* anc;
int* left;
int* right;
STATES** chars;
En este caso anc, left y right son arreglos dinámicos (por eso están declarados como apuntadores, pero se pueden tratar como arrays), y en vez de contener un apuntador, contienen un numero que es el índice del array. Como es un arreglo coordinado, todos los elementos con el mismo índice son parte del mismo nodo. Así, por ejemplo anc [17] = 20; dice que el ancestro del nodo 17, es el nodo 20. En este caso left [20] o right [20] debe ser igual a 17. chars [20] contiene los estados asignados al nodo 20. Si revisan las entradas sobre árboles en wikipedia (en especial las entradas sobre heaps) notaran que esa es la manera con la que tratan los árboles.

Para mi, lo malo de esa aproximación, es que requiere un control más fuerte sobre el mantenimiento de la información (cabe anotar que en el lenguaje de macros TNT, TNT mantiene todo por nosotros, así que eso no es ningún problema!). En mi propia experiencia, yo trabaje en reconciliación de árboles, y trabajaba con alguien que no sabia manejar apuntadores. Así que trate de hacer la implementación con arreglos coordinados. pero en árboles reconciliados no siempre sabemos cuantos nodos nuevos va a tener el árbol, por lo que mantener el código con los arreglos coordinados se hizo bien problemático.

Lo que muestra esa experiencia, no es que una alternativa es mejor que la otra, solo que a veces, la naturaleza del problema puede complicar la forma de resolverlo, y que además algunos estilos de programación se dan mejor que otros, según quien este haciendo el código ;). De aquí en adelante en otros posts, yo voy a usar estructuras.

Para tener algo de código y resaltar la recursividad de los árboles he aquí una función que guarda un árbol en formato parentética (como Hennig86, NONA y 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, ") ");
}
Estas funciones usas las operaciones IO estándar. La función PrintTree coloca la cabecera de un archivo de árboles para Hennig86/NONA/TNT. Y llama a PrintNode. PrintNode chequea si el nodo actual es un terminal, si es así escribe el nombre del terminal, de lo contrario, como en un nodo interno, imprime los paréntesis y llama recirsivamente a PrintNode en cada uno de sus descendientes.

Me parece que Mike Keesey escribió un post sobre árboles bajo programación orientada a objetos, pero no pude encontrarlo, esta es la alternativa más parecida :P

Referencias
[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


Post anteriores en Algoritmos en filogenia
0a. Campos de bits

miércoles, mayo 06, 2009

TNT macros: Intro


One of the most powerful utilities of TNT is its macro's language. This language allows access to TNT internal variables, and tree and data manipulation. This coupled with the computer power of TNT makes a terrific combination.

Unfortunatelly, maybe because of a lack of an extensive manual, or because there are few papers that explicitly use TNT macros (but they are growing!), then this capability is ignored for several users. In this series of posts, I want to give an introduction of the many possibilities allowed by the usage of macros.

Part of this idea born after I see the book "phylogenetics with R" [1], and the constant exchange with Santi Catalano about the TNT's wiki. I hope to post part of this series on the TNT's wiki, but for the moment, just to get the writing feeling, I post it first on my blog ;).

The style is somewhat orientated to programming--I'm powerfully influenced by Kenighan and Ritchie :)--because I think it is the best way to learn the language.

Notation

The TNT macro language is an interpreted language, this is, each instruction is executed as it is readed, then it is possible to “write” the programs in real time (just like the old BASIC, lisp, or the fashionable python).

This can be very nice, for example, to make some simple mathematical operations directly on TNT.

An example:
> macro = ;
Macro language is ON
> var: dest ;
> set dest 4 + 5 ;
> quote 'dest' ;
9
This simple calculator can introduce us in some particularities of the language.

To activate macros, the command macro = ; must be invoked, macro - ; deactivates the macros.

Each instruction finish with a semi colon (;). Although not necessary, for clarity, an space before the semicolon is a good practice.

The keyword var is used to declare variables.

The keyword set assigns a value to a variable. In the example it is a sum. Set accepts the four basic operations, as well as more complex combinations using parenthesis.

The command quote prints on the screen. In TNT to access the value stored into a variable, the name of the variable must be between simple quotes ('). Note that set assigns a value, so the name is written without quotations, but if the value assigned is stored into a variable, then that values must be in quotations:
set dest 'first' * 'second' ;
This operation assigns to dest the value of first times second.

But more important, is that TNT macros can be written on separate file, and executing like any commando of the program.

If the macro is saved with the ".run" extention, and is on the same working directory of the program (for example c:\bin\tnt) it is possible to call the macro just typing the name of the file. For example, TNT is distributed with the macro stats.run, that calculates consistency and retention index of the trees in memory. To execute it, only is necessary to type:
> stats ;
in the program command line, and the macro runs automatically.

As any file of TNT it can be accessed using proc command (but you lost the parameter list) of the macro. It is better to use the command run (or just write the macro's name). As any command macros can be accessed trough the OS command line.

Exercises

As exercises try to use the different ways to call a macro, as well as accessing values with the commands, think on scripts like this:
macro = ;
var:
numIts
numTrees
itTrees;

set numIts 100 ;
set itTrees 10 ;
set numTrees 'numIts' * 'itTrees' ;

rseed 0;
hold 'numTrees' ;
mult = replic 'numIts' hold 'itTrees' ;
nel * ;
p/ ;
That does noting that you can do more easily by hand, but allows you to get familiarity with the notation.

References
[1] Paradis, E. 2006. Analysis of phylogenetics and evolution with R. Springer

Macros de TNT: Intro


Una de las utilidades más poderosas de TNT es su lenguaje de macros. Este lenguaje permite acceder variables internas de TNT, manipular árboles y datos. Esto en conjunto con el poder de computo de TNT es una excelente combinación.

Infortunadamente, bien sea porque no hay un manual extenso, o porque aún son pocos los que han usado explicitamente macros en sus artículos, esta utilidad suele pasar más o menos desapercibida para los usuarios. En esta serie de posts quiero dar un poco de remedio a esta situación, y dar una especie de introducción a las muchas posibilidades que permite el uso de macros.

Parte de la idea nació a que una vez vi el libro de "R para filogenia" [1], y pues a que he estado en contacto con Santi Catalano, que ha estado trabajando en la wiki de TNT. Espero poder postear después parte de esto a la wiki de TNT, pero por ahora, como para ver como me siento escribiendolo, lo coloco primero en mi blog ;).

El estilo, es bien orientado a la programación--mi influencia más poderosa es Kernighan y Ritchie :)--pues la idea es que se le saque el máximo provecho al lenguaje!

Notación

El lenguaje de macros de TNT es un lenguaje interpretado, esto quiere decir que cada instrucción se ejecuta a medida que se va leyendo, por lo que es posible ir "escribiendo" el programa en tiempo real desde la linea de comando.

Esto puede ser muy útil, por ejemplo, para realizar cálculos matemáticos directamente en TNT.

Por ejemplo:
> macro = ;
Macro language is ON
> var: dest ;
> set dest 4 + 5 ;
> quote 'dest' ;
9
Al escribir esa secuencia uno puede tener una calculadora simple. Esto sirve para introducirnos a las cosas particulares del lenguaje.

Para activar los macros, se usa el comando macro = ; con macro - ; se desactivan los macros.

Al terminar cada instrucción se coloca un punto y coma (;). Aunque no es necesario, por claridad es aconsejable dejar un espacio antes del punto y coma.

La palabra clave var se utiliza para declarar las variables.

Con la palabra clave set se asigna un valor a una variable. En este ejemplo el valor de una suma. Set acepta las cuatro operaciones básicas, así como operaciones más complejas usando paréntesis.

El comando quote imprime en pantalla. En TNT, para acceder al valor guardado en las variables, se coloca el nombre de la variable entre comillas simples. Así en este ejemplo, se imprime lo que esta dentro de dest.

Es importante diferenciar entre el valor guardado en la variable y su nombre. Set asigna a una variable, por ello la variable no va dentro de comillas, pero si lo que es asignado es otra variable, entonces como es un valor, debe estar entre comillas:
set dest 'primero' * 'segundo' ;
Esta instrucción asigna a dest la multiplicación de los valores de primero y segundo.

Sin embargo, la aplicación más importante de los macros de TNT es que pueden guardarse en un archivo y ejecutarse como si se tratara de un comando del programa.

Si los macros se guardan con la extensión ".run", y son colocados en el mismo directorio de trabajo del programa (por ejemplo c:\bin\tnt) es posible llamarlo simplemente usando el nombre del archivo. Por ejemplo, en la instalación de TNT esta incluido stats.run, que calcula el indice de consistencia y retención de los árboles en memoria. Para ejecutarlo, simplemente se escribe
> stats ;
en la linea de comando, y el macro es ejecutado automáticamente.

Como todos los archivos de TNT se pueden acceder usando el comando proc, pero se pierde la lista de parámetros (si los hay) del macro. Por eso es mejor usar el commando run (o simplemente escribir el nombre). Como cualquier comando se puede acceder directamente en la lista de parámetros del programa.

Como ejercicios de esta introducción se puede probar las diferentes maneras de invocar los macros, así como acceder valores desde comandos. Por ejemplo scripts como este:
macro = ;
var:
numIts
numTrees
itTrees;

set numIts 100 ;
set itTrees 10 ;
set numTrees 'numIts' * 'itTrees' ;

rseed 0;
hold 'numTrees' ;
mult = replic 'numIts' hold 'itTrees' ;
nel * ;
p/ ;
Que no hacen nada que no se puede hacer directamente (y segura más fácil!), pero sirven para familiarizarse con el lenguaje ;).

Referencias
[1] Paradis, E. 2006. Analysis of phylogenetics and evolution with R. Springer

viernes, mayo 01, 2009

The TNT wiki

Just few days ago, the TNT wiki was released. The project is led by my friend Santi Catalano, Matt Yoder and Ximo Mengual. They meet at the Hennig Meeting, in Tucumán, with the idea of provide a friendly environment for TNT users.

This is the link: http://tnt.insectmuseum.org/

For the most part, the wiki is an on-line version of the program help. But as it is a wiki, it is possible to provide several examples and tricks.

I happy to help with the pages relevant to implied weights: IW and IW quick tutorial pages. And I hope to write somethings for the scripting pages :).

Pablo is not directly involved with the project. Nonetheless, from time to time he reviewed some content, and of course, he always encourage its development. :).

Also, there is a TNT google group: http://groups.google.com/group/tnt-tree-analysis-using-new-technology, open to questions from any user!

A behemoth analysis, and a wiki, this is a great week for TNT :).

Just in case, this is the TNT is a free software for phylogenetic analysis at any scale. It is sponsored by the Willi Hennig society, and can be downloaded here: http://www.zmuc.dk/public/phylogeny/TNT/

La wiki de TNT

Hace unos días, la wiki de TNT fue lanzada al publico. El proyecto es liderado por mi amigo Santi Catalano, Matt Yoder y Ximo Mengual. Ellos se encontraron en el Hennig meeting, en Tucumán, con la idea de proporcionar una ayuda amigable a los usuarios de TNT.

Este es el enlace: http://tnt.insectmuseum.org/

En su mayor parte, la wiki es una versión en linea de la ayuda del programa. Pero como es una wiki es posible agrandarla con muchos ejemplos y trucos.

Yo colabore con las páginas de pesado implicado: la pagina básica y un tutorial rápido. Y espero escribir algunas cosas para la sección de scripts :).

Pablo no esta directamente asociado con el proyecto. Pero, cada cierto tiempo, el revisa algo del contenido, y por supuesto, siempre a apoyado el desarrollo del sitio :).

Hay además un grupo de google de TNT: http://groups.google.com/group/tnt-tree-analysis-using-new-technology abierto a las preguntas de cualquier usuario!

El análisis de un behemot, y una wiki, esta es una gran semana para TNT :).

Por si acaso, TNT es un software gratuito para análisis filogenético a cualquier escala. Esta patrocinado por la Willi Hennig society, y puede descargarse aquí: http://www.zmuc.dk/public/phylogeny/TNT/