Finding the MRCA for two or multiple species across a set of trees

A user recently* contacted me with a question about finding the MRCA for two or multiple species across a set of trees in an R "multiPhylo" object. (*Unfortunately, "recently" is running about two weeks behind, so this email was actually sent on 10/15. To that user, I say "sorry!") She nearly had it correct, but there was a small error in how she was using sapply. Here is the correct way, by way of demonstration:

> require(phytools)
Loading required package: phytools
Loading required package: ape
...
> # simulation some random trees
> trees<-rmtree(N=100,n=10)
> sp<-c("t1","t2") # for example
> # now find the MRCA of sp across the set of trees
> ca<-sapply(trees,function(x,sp) findMRCA(x,sp),sp=sp)
> ca
  [1] 11 19 12 12 12 14 15 14 16 13 12 12 11 16 19 14 17 11
[19] 11 11 11 11 11 11 13 18 19 11 15 12 11 19 12 12 11 12
[37] 11 12 13 11 11 12 12 11 18 12 11 11 15 16 11 16 12 12
[55] 11 14 11 11 15 17 11 14 16 19 11 13 18 11 16 13 12 11
[73] 12 14 13 11 13 11 11 14 17 11 12 11 11 11 11 12 11 13
[91] 11 19 19 11 13 12 12 11 11 11
> # now verify (for the first 3 trees)
> par(mfrow=c(3,1),mar=c(1.5,1.5,1.5,1.5))
> plot.phylo(trees[[1]],cex=1.5); nodelabels(cex=1.5)
> title("tree 1",cex.main=2)
> plot.phylo(trees[[2]],cex=1.5); nodelabels(cex=1.5)
> title("tree 2",cex.main=2)
> plot.phylo(trees[[3]],cex=1.5); nodelabels(cex=1.5)
> title("tree 3",cex.main=2)

Cool. It works.

Version of ancThresh posted

I just posted a version of my function for ancestral character estimation under the threshold model (e.g., 1, 2, 3). Direct link to the code is here. I'm sharing this because there has been some interest since I started posting about the method, but I should note that this is a brand new function & method, and the code is "bleeding edge," so should be used only with caution. (I'm also working on a manuscript describing the approach.)

Some more updates to xkcdTree

After discovering an actual xkcd comic featuring a phylogenetic tree, I was determined to incorporate some of the attributes of this comic tree into the function xkcdTree. Specifically, the real xkcd trees had tilted tip labels, a root edge, and an upward orientation (all things that my xkcd trees lacked). All of these features are also popular plotting options which, now that I've figured them out, should be easier to incorporate into other phytools plotting function - such as plotSimmap - so it's not all about plotting cartoon trees!

Tilted labels are easy, we can just use text(...,srt=tilt). However, they do create the problem that, now that the labels are titled, we need a different amount of space on the x & y axes to fit the labels. We can compute these (unsurprisingly) using strwidth and our 6th-grade trig identities, for instance:

dx<-abs(max(strwidth(cw$tip.label,units="figure",cex=fsize, family="xkcd"))*cos(deg2rad(tilt)))
dy<-abs(max(strwidth(cw$tip.label,units="figure",cex=fsize, family="xkcd"))*sin(deg2rad(tilt)))

deg2rad is an internal function for xkcdTree.

To plot a vertical (i.e., upward facing) tree, all we need to do is flip x & y. We can't forget, though, that the space we have calculated for our labels will also need to be transformed - depending on the dimensions we've allocated for our plotting space.

The new version of the function is here, and it will be in the next release of phytools. Let's try it out:

> require(phytools)
Loading required package: phytools
Loading required package: ape
...
> require("extrafont")
Loading required package: extrafont
> tree<-read.tree(text="((mammals,(birds,reptiles)),amphibians);")
> for(i in 1:tree$Nnode+length(tree$tip)) tree<-rotate(tree,node=i)
> tree<-compute.brlen(tree)
> source("xkcdTree.R")
> tree$root.edge<-0.25
> xkcdTree(tree,file="herpetology.pdf",lwd=2,color="black", dim=c(4,4),jitter=0.001,waver=c(0.03,0.03),tilt=-40, right=F)
**** NOTE: use in Windows requires ...
>

Obviously, plotting cartoon trees is a fairly silly exercise - but I have certainly learned a thing or too that will help me as I augment that plotting functions of the phytools package in the future.

Actual xkcd tree

In all my posts about the function for plotting xkcd style trees (e.g., 1, 2, 3, 4, 5, 6), I somehow overlooked the fact that there is an xkcd comic featuring a phylogenetic tree (more specifically, two trees - although the topological correctedness of the second is debatable). Furthermore, it is one of my favorite xkcd comics of all time:



Who knew that I should have been using this as a guide all along!

In the end, though (and with some tweaking of the settings from their default values, and a little post-hoc editing), my function does a pretty good job:

> require(phytools)
Loading required package: phytools
Loading required package: ape
...
> require("extrafont")
Loading required package: extrafont
> tree<-read.tree(text="((mammals,(birds,reptiles)),amphibians);")
> tree<-compute.brlen(tree)
> xkcdTree(tree,file="herpetology.pdf",lwd=2,color="black", dim=c(4,4),jitter=0.001,waver=c(0.03,0.03))
**** NOTE: use in Windows requires ...
>

That's it!

Estimating Pagel's λ for a set of traits in a matrix or data frame

Earlier, I posted about estimating Pagel's λ jointly for a set of traits. A phytools user recently submitted the following comment:

One question I had was how to calculate lambda on multiple traits? I have code that gets me close in that it outputs the lambda values but not the log values or p-value.

This could be done, of course, using a for loop, something like:

> require(phytools); require(geiger)
> # simulate tree & data
> tree<-pbtree(n=100)
> X<-fastBM(lambdaTree(tree,0.7),nsim=10)
> lambda<-matrix(NA,ncol(X),4,dimnames=list(NULL, c("lambda","logL","logL0","P")))
> if(is.data.frame(X)) X<-as.matrix(X)
> for(i in 1:nrow(lambda))
lambda[i,]<-unlist(phylosig(tree,X[,i],method="lambda", test=T))
> lambda
        lambda      logL     logL0            P
[1,] 0.8176128 -176.9442 -199.7361 1.462441e-11
[2,] 0.5321105 -183.5490 -192.2505 3.023510e-05
[3,] 0.7868250 -177.3286 -200.6997 8.096301e-12
[4,] 0.6772734 -182.8100 -198.7038 1.719778e-08
[5,] 0.4893630 -183.6516 -193.9559 5.634189e-06
[6,] 0.7477579 -167.3094 -194.6396 1.432188e-13
[7,] 0.7160706 -182.1853 -207.2614 1.422751e-12
[8,] 0.6029530 -184.5244 -200.5472 1.505981e-08
[9,] 0.6654206 -181.8013 -201.3363 4.088944e-10
[10,] 0.7553931 -182.0996 -200.9483 8.261227e-10

Alternatively, however, we could also do the following using sapply:

> X<-X[tree$tip.label,]
> X<-as.data.frame(X)
> lambda<-t(sapply(X,phylosig,tree=tree,method="lambda", test=T))
[1] "x has no names; assuming x is in the same order as tree$tip.label"
...
> lambda
   lambda    logL      logL0     P          
V1  0.8176128 -176.9442 -199.7361 1.462441e-11
V2  0.5321105 -183.549  -192.2505 3.02351e-05
V3  0.786825  -177.3286 -200.6997 8.096301e-12
V4  0.6772734 -182.81   -198.7038 1.719778e-08
V5  0.489363  -183.6516 -193.9559 5.634189e-06
V6  0.7477579 -167.3094 -194.6396 1.432188e-13
V7  0.7160706 -182.1853 -207.2614 1.422751e-12
V8  0.602953  -184.5244 -200.5472 1.505981e-08
V9  0.6654206 -181.8013 -201.3363 4.088944e-10
V10 0.7553931 -182.0996 -200.9483 8.261227e-10

That's it.

Posterior density for the thresholds during ancestral character estimation using the threshold model

I've been doing more exploration of ancestral character estimation using the threshold model. One thing that I was interested in was how good we can be at estimating the relative positions of the thresholds along our liability axis. Note that since unobserved liabilities have arbitrary scale, only the relative positions of the thresholds have meaning. Specifically, the set of thresholds [1,2,4] with Brownian rate of evolution of the liability of σ²=1.0 has the exact same meaning & interpretation of thresholds [0,2,6] and σ²=2.0.

Nonetheless, since I have simulated liability evolution with a fixed, known σ² (and it is the same σ² that I will use for estimation), it makes sense in this rare case to look at the posterior density for the thresholds on its raw scale. I decided to simulate using a four state (thus 3 threshold) model, in which the four states are "blue", "green", "red", and "yellow"; and the threshold transition points between colors are 0 (blue to green), 1 (green to red) and 4 (red to yellow).

Here was my simulation:

> # load source
> source("ancThresh.R")
> # simulate tree
> tree<-pbtree(n=100,scale=10)
> # simulate liability
> l<-fastBM(tree,internal=TRUE,a=2)
> # translate to discrete character
> x<-rep("blue",length(l))
> x[l>0]<-"green"
> x[l>1]<-"red"
> x[l>4]<-"yellow"
> names(x)<-names(l)
> sequence<-c("blue","green","red","yellow")
> cols<-sequence; names(cols)<-sequence
> res<-ancThresh(tree,x[1:length(tree$tip)],seq=sequence, ngen=1000000,control=list(piecol=cols))
MCMC starting....
gen 1000
gen 2000
...

The posterior sample for the thresholds is in res$par. Unlike in the translation, above, each threshold is named for its upper bound (so, blue = 0; green = 1; red = 4; and yellow = Inf, in this case). Obviously, and again because liabilities are unscaled and uncentered, we fix the lowest of these at a constant value - we choose 0, but this is arbitrary.

Now let's do some computations to plot the posterior density:

> breaks<--0.25:25/5
> breaks
[1] -0.05  0.15  0.35  0.55  0.75  0.95  ...
> green<-hist(res$par[,"green"],breaks=breaks)
> red<-hist(res$par[,"red"],breaks=breaks)
> plot(red$mids,red$density,type="s",ylim=c(0, max(c(red$density,green$density))),xlab="liability", ylab="density",col="red")
> lines(red$mids,green$density,type="s",col="green")

I have also added the generating threshold (colored on each side with the appropriate discrete state) for fun. Not too bad.

It should also be noted that if we look at the ratio of the mean estimates of the two thresholds this is also quite close to the generating ratio of four:

> mean(res$par[,"red"])/mean(res$par[,"green"])
[1] 4.3913

That's it.

Ancestral character estimation under the threshold model, part II

Another quick post about ancestral character estimation under the threshold model from quantitative genetics. For a reminder of the details of this model, please check out my prior posts on the topic, as well as Felsenstein (2012).

The central idea is that this model, in which the value for a discrete character is determined by the evolution of an underlying continuous trait called 'liability,' might be a more appropriate model for many discretely manifested, by complex or polygenic, attributes of organismal phenotype.

I'm estimating under this model by using Bayesian MCMC to sample ancestral and tip liabilities, and the thresholds that result in a change in the discrete character, from their joint posterior probability distribution.

I'm still working on this, but here is an example simulation, analysis, and then some assessments of performance:

> source("ancThresh.R")
> # simulate liabilities
> l<-fastBM(tree,internal=TRUE,sig2=2,a=-1)
> # translate liabilities to threshold character
> x<-rep("blue",length(l))
> x[l>0]<-"green"
> x[l>1]<-"red"
> x[l>6]<-"yellow" > names(x)<-names(l)
> x
    t54      t55      t51 ...
 "blue"   "blue"   "blue" ...
> x11(); res<-ancThresh(tree,x[1:length(tree$tip)],ngen=100000)
**** NOTE: no sequence provided, using alphabetical or numerical order
MCMC starting....
gen 1000
gen 2000
...

> # plot the likelihood profile
> plot(res$par[,"gen"],res$par[,"logLik"],xlab="gen", ylab="logL")
> lines(res$par[,"gen"],res$par[,"logLik"])

> # pull estimates as max PP for each node
> est<-apply(res$ace,1,function(x) names(x)[which(x==max(x))])
> est
    101      102      103 ...
  "red"   "blue"   "blue" ...
> # do they match?
> a<-x[101:199] # get simulated ancestors
> sum(a==est)/length(a)
[1] 0.7878788 # about 79% match
> # what if we restrict to >0.5, 0.7, or 0.9 PP?
> any50<-apply(res$ace,1,function(x) any(x>0.5))
> sum(any50)/tree$Nnode # gives us the % nodes with PP>0.5
[1] 0.969697
> # now the fraction that match, given PP>0.5
> sum(a[any50]==est[any50])/sum(any50)
[1] 0.7916667
> # PP > 0.7
> any70<-apply(res$ace,1,function(x) any(x>0.7))
> sum(any70)/tree$Nnode
[1] 0.5050505
> sum(a[any70]==est[any70])/sum(any70)
[1] 0.96
> # 96% of nodes with PP>0.7 are correct
> # PP > 0.9
> any90<-apply(res$ace,1,function(x) any(x>0.9))
> sum(any90)/tree$Nnode
[1] 0.2929293
> sum(a[any90]==est[any90])/sum(any90)
[1] 1

Cool. I will post the code & more on this very soon.