Followers

Sunday, April 23, 2017

ggtree function in R

As a student in the Organismic and Evolutionary Biology program, I finally have the occasion to compare trait data with phylogenetic data. Over the past couple of weeks, I have been searching for a way to put these things together graphically.

To make the tree of plant species, I have used the "phylomatic" function in R package "brranching" by Scott Chamberlain.

To associate data with the tree, I used Guangchuang Yu's ggtree package. This makes faceted plots where one facet is the tree and another facet (or more) show the data.

First I will offer my thanks to these gurus who make things possible for us R mortals.

Here was the script I came up with:

library(ggplot2)
library(brranching)

#Make the tree:
taxa.plus<- c("Aesculus carnea", "Antirrhinum majus",
              "Brassica napus", "Catalpa speciosa", "Citrus sinensis",
              "Cucurbita pepo", "Dicentra eximia", "Digitalis purpurea",
              "Echium vulgare","Eupatorium perfoliatum", "Fragaria ananassa",
              "Geranium maculatum", "Helianthus annuus",
              "Impatiens capensis", "Kalmia latifolia",
              "Linaria vulgaris", "Lobelia siphilitica",
              "Lythrum salicaria", "Malus domestica",
              "Monarda_didyma",
              "Penstemon digitalis", "Persea americana",
              "Prunus dulcis", "Rhododendron_prinophyllum",
              "Silene vulgaris",
              "Solanum carolinense", "Solidago canadensis",
              "Tilia_cordata", "Thymus vulgaris",
              "Trifolium pratense", "Vaccinium corymbosum",
              "Verbascum thapsis")
treeget<-phylomatic(taxa = taxa.plus, get="GET")
class(treeget)<-"phylo"
plot(treeget)



#Generate data for the graph:
library(plyr)
library(dplyr)
df1<-data_frame(Species=treeget$tip.label,
                Height=runif(n=length(treeget$tip.label), min=10, max=100),
                Chem=rnorm(n=length(treeget$tip.label), mean=5, sd=3))


#Plot in ggtree
library(ggtree)
Tree<-ggtree(as.phylo(treeget)) + geom_tiplab()
#change axis limits to view tip labels
ggtree(as.phylo(treeget)) + geom_tiplab() +
  xlim_tree(30)

#Adjust limits on tree
Tree.l<-Tree + geom_tiplab() + xlim(0,22)

#Creat facet plot:
fp<-facet_plot(Tree, panel = "Chem", data = df1,
           geom = geom_point, aes(x=Chem),
           stat= 'identity') + ggtree::xlim_expand(xlim=c(-2,2), panel="Chem")
fp



#adjust axis limits:
fp.limited<-fp + xlim_tree(40) + xlim_expand(xlim=c(-5, 25), panel = "Chem") +
  theme_tree2()
fp.limited




#Add geom_pointrange:
df1$error<-rnorm(n=length(treeget$tip.label), mean=2, sd=0.4)
fp.error<-facet_plot (fp.limited, geom=geom_pointrangeh, data=df1,
                      aes(x=Chem, xmax = Chem + error, xmin = Chem - error),
                      panel = "Chem")
fp.error





#Try with geom_errobar instead
fp.error.bar<-facet_plot (fp.limited, geom=geom_errorbarh, data=df1,
                      aes(x=Chem, xmax = Chem + error, xmin = Chem - error),
                      panel = "Chem")
fp.error.bar
#Warning: Ignoring unknown aesthetics: x


#However, you get error without specification of "x"
fp.error.bar2<-facet_plot (fp.limited, geom=geom_errorbarh, data=df1,
                          aes(#x=Chem,
                              xmax = Chem + error, xmin = Chem - error),
                          panel = "Chem")
fp.error.bar2
#Error in FUN(X[[i]], ...) : object 'x' not found





Of course I haven't used any real data, but you could!
Good luck plotting!
Evan