A few days ago, I received the following message from a phytools user.
“I am using the make.simmap function in your R package phytools to get
histories of chromosome number evolution in several taxonomic clades. Some
clades I am interested in have a tight variation in chromosome number…,
and some clades have a considerable variation in chromosome number. In
clades with less variation, the number of states I have in the Q matrix
is small (about eight states), but in clades with considerable variation in
chromosome number, sometimes the Q matrix can have up to 100 states. In
cases where I have a large Q matrix, the make.simmap function takes a
long time to run, and will not generate a stochastic map even after many
days (even with the number of simulations set to 1). It would be much
helpful if you could provide any insight on how I can get the make.simmap
function to work with a large number of states.”
This isn't the first time I've encountered this kind of problem, and what I usually suggest is a model of chromosome evolution in which we assume one rate of increase in chromosome number (the “fission” rate), a second rate decrease in chromosome number (the “fusion” rate), and, lastly, a third rate of chromosome doubling (the whole-genome duplication rate). If we suppose that this is how chromosome number tends to evolve, then we can usually fit this model to our data regardless* of how many levels our character can assume. (*Almost. We need to keep in mind that Q grows with the number of levels of our trait, k, not with the number of parameters in the fitted model. At some point this will make matrix exponentiation of Q computationally prohibitive and this will begin to cause other problems with estimation.)
To see how this might work, I'm going to first imagine (and simulate) a phylogenetic tree under our proposed model. In this case, I'll make my phylogeny 100 units (say, 100 million years) in total depth, and I'll set my fission rate to 0.1 (i.e., one chromosome splitting every 10 million years), my fusion rate to 0.05 (one fusion every 20 million years on average), and my whole genome duplication rate to 0.01 (once every 100 million years). Keep in mind that even though our total tree depth is only 100, we expect more than one whole genome duplication events because the amount of evolutionary time in the tree is the sum of the branch lengths, not the total tree length.
library(phytools)
tree<-pbtree(n=200,scale=100)
plotTree(tree,ftype="off",lwd=1,mar=c(2.6,2.1,1.1,1.1))
axis(1,at=c(0,25,50,75,100))
For this example, I'll set the minimum number of chromosome (haploid number) to 2 and the maximum number to 30. I tested this same code with up to 100 chromosome and it worked fine. (Although it was a bit slower, of course!)
min.n<-2
max.n<-30
levs<-min.n:max.n
k<-length(levs)
Q<-matrix(0,k,k,dimnames=list(levs,levs))
for(i in 1:k) for(j in 1:k){
if(levs[i]==(levs[j]+1)) Q[i,j]<-0.05
if(levs[i]==(levs[j]-1)) Q[i,j]<-0.1
if(levs[j]==(2*levs[i])) Q[i,j]<-Q[i,j]+0.01
}
diag(Q)<--rowSums(Q)
To see what our simulation model looks like, I'm going to first assign
it a class attribute and then graph it using one of phytools built-in
plot methods.
class(Q)<-"Qmatrix"
plot(Q,show.zeros=FALSE,mar=rep(0,4),cex.traits=0.7,
spacer=0.05)
Awesome. This plot shows us exactly the model we had planned to simulate.
Now let's simulate a trait vector, x, using the phytools function
sim.Mk as follows.
x<-sim.Mk(tree,Q,anc=2)
x
## t31 t32 t28 t143 t146 t173 t174 t9 t66 t135 t150 t151 t175 t176 t44 t169 t170 t123 t124 t165 t166 t110 t111 t48 t155 t156 t29 t20 t131 t132 t22 t54 t55
## 30 27 27 12 11 10 10 12 6 9 10 9 7 7 13 12 12 21 20 3 3 9 4 7 4 6 3 4 7 7 8 7 9
## t119 t120 t27 t19 t127 t128 t23 t24 t101 t102 t1 t56 t157 t158 t181 t182 t41 t193 t194 t89 t10 t21 t86 t87 t8 t5 t112 t113 t96 t97 t79 t80 t75
## 4 4 6 18 9 9 5 8 12 10 3 3 3 2 8 8 10 5 6 9 19 19 20 21 21 21 24 23 22 23 8 8 24
## t117 t118 t18 t77 t78 t74 t53 t63 t64 t11 t49 t65 t197 t198 t67 t13 t14 t33 t34 t16 t17 t26 t58 t59 t144 t145 t189 t190 t57 t47 t35 t12 t103
## 11 12 10 11 22 19 15 27 26 30 14 14 17 17 14 18 17 16 17 17 14 19 14 15 11 11 11 11 13 20 9 22 9
## t115 t116 t84 t85 t15 t104 t105 t46 t125 t126 t72 t73 t45 t195 t196 t114 t138 t139 t42 t6 t7 t38 t185 t186 t81 t106 t107 t76 t60 t61 t62 t51 t140
## 11 9 20 12 11 12 15 7 6 6 4 6 5 3 3 5 6 6 8 8 9 10 9 10 11 24 11 23 12 11 13 9 8
## t191 t192 t68 t70 t71 t129 t130 t3 t2 t159 t160 t199 t200 t52 t4 t94 t95 t90 t40 t179 t180 t177 t178 t50 t43 t121 t122 t100 t91 t92 t88 t187 t188
## 8 8 11 11 9 25 24 9 13 9 10 12 12 16 21 9 8 9 9 5 5 5 5 5 6 2 2 4 4 6 5 6 7
## t98 t82 t83 t25 t36 t37 t108 t109 t39 t142 t152 t153 t163 t164 t30 t69 t167 t168 t99 t133 t134 t148 t154 t171 t172 t147 t183 t184 t141 t161 t162 t149 t136
## 6 3 8 2 14 10 11 12 4 4 4 8 6 7 15 15 11 11 11 15 16 17 18 17 17 12 11 10 12 8 9 8 12
## t137 t93
## 13 11
## Levels: 10 11 12 13 14 15 16 17 18 19 2 20 21 22 23 24 25 26 27 3 30 4 5 6 7 8 9
One thing you might notice is that the number of levels can be smaller than the number of possible values of the trait in our case, we can find the observed number of levels by computing…
length(levels(x))
## [1] 27
… and we might find that it is smaller than the number of possible levels (29 in our example). We want to model, however, all 29 levels of our trait – not just those levels that we've observed!
To do that, we'll first convert our data into a 200 (row) by 29 (column)
binary matrix, in which each row/column combination is 1 if the column trait
value is observed for the row species, and 0 otherwise. Fortunately,
phytools::to.matrix does this calculation for us.
levs<-2:30
X<-to.matrix(x,seq=levs)
Now we're ready to build our design matrix for estimation. We'll do this
much as we built our Q matrix for simulation, above, but this time
populating it with only the integer values of 1, 2, or 3 for our
three different rate classes, as follows.
k<-length(levs)
model<-matrix(0,k,k,dimnames=list(levs,levs))
for(i in 1:k) for(j in 1:k){
if(levs[i]==(levs[j]+1)) model[i,j]<-1
if(levs[i]==(levs[j]-1)) model[i,j]<-2
if(levs[j]==(2*levs[i])) model[i,j]<-3
}
We can verify that we've built our model correctly by graphing it.
tmp<-model
class(tmp)<-"Qmatrix"
plot(tmp,show.zeros=FALSE,mar=rep(0,4),cex.traits=0.7,
spacer=0.05)
Now let's proceed to fit our model to our tree and data. This stage is a bit computationally intensive, so we should be prepared to wait a short while.
fitted.model<-fitMk(tree,X,model=model,pi="fitzjohn")
plot(fitted.model,show.zeros=FALSE,mar=rep(0,4),cex.traits=0.7,
spacer=0.03)
We should see that our fitted model is actually a pretty good approximation of the generating process.
Now, I'll generate a single stochastic map. Since I already fit my model
using fitMk, to save time I'll use a fixed value of Q; however, we
could've also let make.simmap both estimate and simulate.
fixedQ<-unclass(as.Qmatrix(fitted.model))
system.time(
smap<-make.simmap(tree,X,Q=fixedQ,
pi="fitzjohn",nsim=1)
)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 2 -0.1206487 0.1108141 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 3 0.0655518 -0.1862005 0.11081411 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 4 0.0000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 5 0.0000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 6 0.0000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000
## 7 0.0000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000
## 8 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 9 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 10 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000
## 11 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000
## 12 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000
## 13 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000
## 14 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141
## 15 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005
## 16 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518
## 17 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 18 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 19 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 20 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 21 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 22 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 23 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 24 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 25 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 26 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 27 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 28 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 29 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 30 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29
## 2 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 3 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 4 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 5 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 6 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 7 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 8 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 9 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 10 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 11 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 12 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 13 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000
## 14 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000
## 15 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 16 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 17 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 18 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 19 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 20 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 21 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 22 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 23 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 24 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000
## 25 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000
## 26 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000
## 27 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000
## 28 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141
## 29 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659
## 30 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518
## 30
## 2 0.00000000
## 3 0.00000000
## 4 0.00000000
## 5 0.00000000
## 6 0.00000000
## 7 0.00000000
## 8 0.00000000
## 9 0.00000000
## 10 0.00000000
## 11 0.00000000
## 12 0.00000000
## 13 0.00000000
## 14 0.00000000
## 15 0.00983463
## 16 0.00000000
## 17 0.00000000
## 18 0.00000000
## 19 0.00000000
## 20 0.00000000
## 21 0.00000000
## 22 0.00000000
## 23 0.00000000
## 24 0.00000000
## 25 0.00000000
## 26 0.00000000
## 27 0.00000000
## 28 0.00000000
## 29 0.11081411
## 30 -0.06555180
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
## 2 3 4 5 6 7 8 9 10 11 12 13
## 2.940788e-01 2.827556e-01 2.225111e-01 1.300080e-01 5.276241e-02 1.461851e-02 2.830492e-03 3.929497e-04 3.918046e-05 2.792319e-06 1.429001e-07 5.331037e-09
## 14 15 16 17 18 19 20 21 22 23 24 25
## 1.480376e-10 3.132127e-12 5.336493e-14 7.184468e-16 7.846039e-18 7.146190e-20 5.580610e-22 3.827361e-24 2.343164e-26 1.289252e-28 6.370658e-31 2.815516e-33
## 26 27 28 29 30
## 1.108554e-35 3.880657e-38 1.207964e-40 3.351626e-43 9.721738e-46
## Done.
## user system elapsed
## 0.93 0.00 0.92
smap
##
## Phylogenetic tree with 200 tips and 199 internal nodes.
##
## Tip labels:
## t31, t32, t28, t143, t146, t173, ...
##
## The tree includes a mapped, 29-state discrete character
## with states:
## 10, 11, 12, 13, 14, 15, ...
##
## Rooted; includes branch lengths.
Now I'll plot my single stochastic character mapped tree as follows.
cols<-setNames(grey.colors(k,0,1,rev=TRUE),levs)
plot(smap,cols,outline=TRUE,lwd=2,ftype="off",ylim=c(-10,Ntip(tree)))
add.color.bar(50,cols,title="chromosome number",lims=c(2,30),
x=0,y=-10,prompt=FALSE,subtitle="")
Obviously, this is just one stochastic map. Let's try to generate 100.
smaps<-make.simmap(tree,X,Q=fixedQ,pi="fitzjohn",nsim=100)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 2 -0.1206487 0.1108141 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 3 0.0655518 -0.1862005 0.11081411 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 4 0.0000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 5 0.0000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 6 0.0000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000
## 7 0.0000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000
## 8 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 9 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 10 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000
## 11 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000 0.00000000 0.0000000
## 12 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141 0.00000000 0.0000000
## 13 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005 0.11081411 0.0000000
## 14 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.18620054 0.1108141
## 15 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1862005
## 16 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518
## 17 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 18 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 19 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 20 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 21 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 22 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 23 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 24 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 25 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 26 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 27 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 28 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 29 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 30 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29
## 2 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 3 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 4 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 5 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 6 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 7 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 8 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 9 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 10 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 11 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 12 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 13 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000 0.00000000 0.0000000
## 14 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00983463 0.0000000
## 15 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 16 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 17 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 18 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 19 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 20 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 21 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 22 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 23 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000
## 24 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000 0.00000000 0.0000000
## 25 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000 0.00000000 0.0000000
## 26 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141 0.00000000 0.0000000
## 27 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659 0.11081411 0.0000000
## 28 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518 -0.17636591 0.1108141
## 29 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.06555180 -0.1763659
## 30 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0655518
## 30
## 2 0.00000000
## 3 0.00000000
## 4 0.00000000
## 5 0.00000000
## 6 0.00000000
## 7 0.00000000
## 8 0.00000000
## 9 0.00000000
## 10 0.00000000
## 11 0.00000000
## 12 0.00000000
## 13 0.00000000
## 14 0.00000000
## 15 0.00983463
## 16 0.00000000
## 17 0.00000000
## 18 0.00000000
## 19 0.00000000
## 20 0.00000000
## 21 0.00000000
## 22 0.00000000
## 23 0.00000000
## 24 0.00000000
## 25 0.00000000
## 26 0.00000000
## 27 0.00000000
## 28 0.00000000
## 29 0.11081411
## 30 -0.06555180
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
## 2 3 4 5 6 7 8 9 10 11 12 13
## 2.940788e-01 2.827556e-01 2.225111e-01 1.300080e-01 5.276241e-02 1.461851e-02 2.830492e-03 3.929497e-04 3.918046e-05 2.792319e-06 1.429001e-07 5.331037e-09
## 14 15 16 17 18 19 20 21 22 23 24 25
## 1.480376e-10 3.132127e-12 5.336493e-14 7.184468e-16 7.846039e-18 7.146190e-20 5.580610e-22 3.827361e-24 2.343164e-26 1.289252e-28 6.370658e-31 2.815516e-33
## 26 27 28 29 30
## 1.108554e-35 3.880657e-38 1.207964e-40 3.351626e-43 9.721738e-46
## Done.
smaps
## 100 phylogenetic trees with mapped discrete characters
I can generate a summary of this set of stochastic maps using another phytools method. Unfortunately, given the large number of states this can be a bit slow too.
obj<-summary(smaps)
plotTree(tree,ftype="off",lwd=1,ylim=c(-10,Ntip(tree)))
add.color.bar(50,cols,title="chromosome number",lims=c(2,30),
x=0,y=-10,prompt=FALSE,subtitle="")
aa<-matrix(0,nrow(obj$ace),k,dimnames=list(rownames(obj$ace),
levs))
for(i in 1:ncol(obj$ace)) aa[,colnames(obj$ace)[i]]<-obj$ace[,i]
nodelabels(pie=aa,piecol=cols,cex=0.5)
Lastly, let's try to plot all our stochastic maps at once using a
visualization similar to a "densityMap" plot. Because this involves
transparency, it takes a while to run!
plotTree(tree,color="black",lwd=4,ftype="off",
ylim=c(-10,Ntip(tree)))
add.color.bar(50,cols,title="chromosome number",lims=c(2,30),
x=0,y=-10,prompt=FALSE,subtitle="")
plotTree(tree,color="white",lwd=2,ftype="off",add=TRUE,
ylim=c(-10,Ntip(tree)))
COLS<-make.transparent(cols,1/length(smap))
for(i in 1:length(smaps)){
plot(smaps[[i]],COLS,lwd=2,ftype="off",add=TRUE,
ylim=c(-10,Ntip(tree)))
}
Overall, that's pretty cool!
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.