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")
