Wednesday, November 6, 2019

Simple wrapper function to check a "phylo" object

Recently the following inquiry was sent to the R-sig-phylo email list:

“I am working with the phylo object from the ape package in my own package in which I am manipulating the trees. I would like to check that I have successfully created a valid ape object, but the checkValidPhylo function appears to be solely interactive - it prints out a display, and always returns NULL. Is there a function in the ape package (or can there be?) that would do the checks, but return the results in a way that I can then process in my function? (e.g. return a vector of each of the checks as TRUE/FALSE) (And I don't want anything printed out, since I don't want that output to be printed for users of my function).”

I suggested the following solution (which turned out to be satisfactory).

This simple function captures the printed output of checkValidPhylo, which can find MODERATE, FATAL, or no errors. If no error is detected, the function returns a 0. If at least one FATAL error is detected, it returns a 2, otherwise (if only MODERATE errors are found), it returns a 1:

chk.phylo<-function(x){
    object<-capture.output(checkValidPhylo(x))
    if(length(grep("FATAL",object))>0) return(2)
    else if(length(grep("MODERATE",object))>0) return(1)
    else return(0)
}

Let's try it out:

library(ape)
t1<-rtree(n=10) ## good tree
chk.phylo(t1)
## [1] 0
t2<-t1
t2$Nnode<-9.5 ## tree has non-integer Nnode
## should cause a MODERATE error
chk.phylo(t2)
## [1] 1
t3<-t1
t3$edge<-t3$edge[-4,] ## tree is missing an
## edge, should cause a FATAL error
chk.phylo(t3)
## [1] 2

That's it.

Saturday, November 2, 2019

Showing ancestral state on a phylo.to.map plot

Today I received the following inquiry:

“I am doing phylogeographic analyses of pathogenic bacteria in …, and have found some of the visualization methods helpful for understanding outbreak dynamics. Now here comes the question or maybe a feature request: Is it possible to combine ”plot(summary(simmap.trees)))“ with the ”phylo.to.map(phylo.tree)“ function, to project the ancestral state reconstruction to a world map?”

Indeed, this is already possible because nodelabels( ) & tiplabels( ) from the ape package play perfectly well with phylo.to.map.

One very important word of caution is required, however. Because phylo.to.map( ) does node rotations by default, it will be necessary to match the nodes from the stochastic character mapped trees to the nodes of the rotated phylogeny in the phylo.to.map object.

The easiest way to do that is just to: (1) run phylo.to.map, then; (2) pull out the tree from the resultant "phylo.to.map" object & use that tree for stochastic mapping.

Here's how that might look (using simulated data):

library(phytools)
## here's our tree & data
World.tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  W.ubvct, N.crk, O.tigbryo, A.uxgbavw, S.igzt, M.wuvrb, ...
## 
## Rooted; includes branch lengths.
World
##                      lat        long
## W.ubvct      -58.1670987 -107.067702
## N.crk        -62.4831489 -163.553854
## O.tigbryo    -19.6762651  -65.235665
## A.uxgbavw     -4.4481564  -15.348011
## S.igzt         0.4478175   -7.230006
## M.wuvrb      -24.4885168   76.791949
## I.sfgjcqyphr -27.4762961   90.845604
## P.wtlcrszkfi -26.8759691   56.549824
## Y.kyaltcfq   -23.0423582   71.257658
## N.ojt        -15.6057166   79.638752
## Y.kuvywit    -29.7095111   77.569301
## M.jcfvrbq    -33.2212398   59.108004
## G.pbhsac     -11.1126509   73.307996
## L.tcpsneijl   -2.6186920   60.174359
## A.gprdmcbq    -1.9270332   53.704890
## F.qfpenc      -4.6062239   15.663828
## L.erq         -9.9327773   14.033582
## W.fxeus       75.6131818  -24.214185
## E.qbpilzy     48.7607643 -101.419331
## I.fzuryvaxs   22.6740158  -31.775542
## G.ctas        37.4639540    7.982519
## X.oflhdrym    40.3753670  -15.282336
## Y.baevguzk    35.5968790   -3.151682
## Z.dczjkolt    -6.9465315 -132.877519
## W.uwmjvdal     7.4247476 -132.126099
## A.vnosymup    -0.9018462  -58.811375
x ## discrete character
##      W.ubvct        N.crk    O.tigbryo    A.uxgbavw       S.igzt 
##            b            b            a            a            a 
##      M.wuvrb I.sfgjcqyphr P.wtlcrszkfi   Y.kyaltcfq        N.ojt 
##            b            b            b            b            b 
##    Y.kuvywit    M.jcfvrbq     G.pbhsac  L.tcpsneijl   A.gprdmcbq 
##            b            b            b            b            b 
##     F.qfpenc        L.erq      W.fxeus    E.qbpilzy  I.fzuryvaxs 
##            b            a            a            a            a 
##       G.ctas   X.oflhdrym   Y.baevguzk   Z.dczjkolt   W.uwmjvdal 
##            a            a            a            b            b 
##   A.vnosymup 
##            b 
## Levels: a b
## first phylo.to.map
map.object<-phylo.to.map(World.tree,World,plot=FALSE)
## objective: 226
## objective: 172
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 142
## objective: 142
## now stochastic mapping
mtrees<-make.simmap(map.object$tree,x,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##              a            b
## a -0.007975089  0.007975089
## b  0.007975089 -0.007975089
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5
## Done.
## compute node probabilities
smap.object<-summary(mtrees)

## plot
cols<-setNames(c("blue","red"),levels(x))
plot(map.object,ftype="off",pts=FALSE)
nodelabels(pie=smap.object$ace,
    piecol=cols,cex=0.6)
points(World[,2:1],bg=cols[x[rownames(World)]],
    cex=1.5,pch=21)

plot of chunk unnamed-chunk-1

That's not bad, right?

FYI, the tree & data were simulated as follows:

World.tree<-pbtree(n=26,scale=100)
World.tree$tip.label<-replicate(Ntip(World.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
lat<-fastBM(World.tree,sig2=10,bounds=c(-90,90))
long<-fastBM(World.tree,sig2=80,bounds=c(-180,180))
World<-cbind(lat,long)
Q<-0.01*matrix(c(-1,1,1,-1),2,2,
    dimnames=list(c("a","b"),c("a","b")))
x<-sim.Mk(World.tree,Q)

Wednesday, October 23, 2019

Updated prototype Ishihara color detection test R package & new app

This morning I worked on developing a prototype web application for the Ishira color perception test based on R & the shiny web application framework.

To do this I had to make some small updates to my very basic (so far) R package ishihara.

The first of these involved identifying the pattern-forming circles in the one Ishihara plate that I'd digitized (plate 2, I believe).

Unfortunately, this was not as straightforward as it seems as I myself am abnormally trichromatic so I can't actually see the pattern & had to ask for help.

Here's my helper at work:

The result is now that the pattern can now be visualized by even a colorblind person by applying a transparency filter to the non pattern-forming circles.

This is what I mean:

library(ishihara)
## 
##    *************************************************
##    *                                               *
##    *    This package is just for fun & is not a    *
##    *           medical diagnostic tool.            *
##    *                                               *
##    ***********************************************
data(plate2)
plot(plate2) ## normal

plot of chunk unnamed-chunk-1

plot(plate2,alpha=0.2) ## 80% transparent

plot of chunk unnamed-chunk-1

Users can even create a neat animation where transparency is faded on & off. Try it:

for(i in 50:0){
    plot(plate2,alpha=min(1,i/50))
    Sys.sleep(0.1)
}
for(i in 1:50){
    plot(plate2,alpha=min(1,i/50))
    Sys.sleep(0.1)
}

Check out the app and let me know what you think!