Hack to fix some or all internal node values in make.simmap (versión español)

Hace unos días alguien me preguntó el siguiente:

“Existe en phytools alguna función que permita asignar un determinado estado a los nodos ancestrales? Es decir, se puede crear un simmap según como uno quiera? O tendría que por ejemplo modificar los mapped.edge directamente? Sólo encontré la función getStates. No existe una que sea setStates?”

La respuesta es que no, hasta ahora no existe ninguna opción en la función make.simmap para especificar ni valores fijos ni probabilidades previas por los nodos interiores en un análisis en que estamos generando mapas estocásticas por la evolución de un rasgo discreto sobre el árbol - con la excepción del raíz global que si podemos ya fijar utilizando el argumento pi.

Pero… si hay una manera en que podemos hacerlo y sea posible que yo pueda adicionar esta funcionalidad como opción en make.simmap al futuro. En este momento, la sola manera de hacerlo es a través de un poquito de codificación básica en R, aka. "scripting."

La idea con esta solución es que voy a adicionar ramas de longitud cero a cada nodo interior por cuáles sabemos sus estados. Como ninguna evolución puede ocurrir a través de una rama de longitud de cero, esta manipulación nos permite de efectivamente fijar los estados por los nodos interiores correspondientes porque tendrán que asumir los valores que pondremos en los tips nuevos de nuestro árbol modificado. Al final, después de correr nuestro análisis make.simmap, podemos remover todas esas ramas adicionadas de nuestros árboles finales.

En este ejemplo, estoy imaginando que tengamos valores que queramos fijar por todos los nodos interiores - pero el lógico, y el algoritmo computacional, deben funcionar igual por casos un que deseemos fijar ciertos pero no todos los nodos.

Empezando con un árbol (tree) y los valores por las especies (en x) y los nodos interiores (en y), continuamos así:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "b" "b" "b" "b" "b" "a" "a" "b" "a" "a" "b" "b" "b" "a" "a" "a" "a" "b" 
##   S   T   U   V   W   X   Y   Z 
## "b" "b" "b" "b" "a" "b" "a" "a"

y

##  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44 
## "a" "a" "b" "b" "b" "b" "a" "a" "b" "b" "b" "b" "a" "b" "a" "a" "b" "b" 
##  45  46  47  48  49  50  51 
## "b" "b" "b" "b" "b" "b" "b"

dotTree(tree,x,colors=setNames(c("blue","red"),c("a","b")))
nodelabels(pie=to.matrix(y,c("a","b")),piecol=c("blue","red"),cex=0.5)

for(i in 1:length(y)){
    if(i==1) new.tree<-tree
    M<-matchNodes(tree,new.tree,method="distances")
    node<-M[which(as.numeric(names(y)[i])==M[,1]),2]
    new.tree<-bind.tip(new.tree,names(y)[i],edge.length=0,where=node)
}

Ya tenemos nuestro árbol modificado. Solo para tener una idea de la manipulación que hemos hecho podemos visualizar el árbol:

plotTree(new.tree)

Y finalmente estamos listos para correr nuestro análisis make.simmap:

trees<-make.simmap(new.tree,c(x,y),nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##           a         b
## a -1.186413  1.186413
## b  1.186413 -1.186413
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5

## Done.

trees<-lapply(trees,drop.tip.simmap,tip=names(y))
class(trees)<-c("multiSimmap","multiPhylo")
trees

## 100 phylogenetic trees with mapped discrete characters

Al final, tenemos un objeto multiSimmap en que hemos fijado los estados de todos o algunos de los nodos interiores y sobre que podemos realizar cualquier análisis subsecuente disponible. Por ejemplo:

obj<-summary(trees)
obj

## 100 trees with a mapped discrete character with states:
##  a, b 
## 
## trees have 18.46 changes between states on average
## 
## changes are of the following types:
##       a,b   b,a
## x->y 7.73 10.73
## 
## mean total time spent in each state is:
##              a         b    total
## raw  6.2708890 9.5060241 15.77691
## prop 0.3974725 0.6025275  1.00000

plot(obj,colors=setNames(c("blue","red"),c("a","b")))

densityMap(trees,lwd=6)

## sorry - this might take a while; please be patient

Chévere!

Los datos y el árbol por este ejemplo fueron simulados utilizando el código que sigue:

tree<-pbtree(n=26,tip.label=LETTERS,scale=2)
Q<-matrix(c(-1,1,1,-1),2,2,dimnames=list(letters[1:2],letters[1:2]))
ttree<-sim.history(tree,Q)
x<-getStates(ttree,"tips")
y<-getStates(ttree,"nodes")