R code to generate a phylogeny with associated trait values

The main content of this post is, as the title says, some R code that can be used to generate an image of a phylogeny with representations of trait values of the tips.

But first, some back-story: earlier this year, U of M Plant Biology grad student Derek Nedveck spent about an hour telling me about the virtues of R, the knitr package, open lab notebooks, and replicable science. I enjoyed the conversation but put it aside at the time, mainly since I wasn’t actively working on any particular project (I had just graduated, ok?) and had no knowledge of R.

I’ve learned quite a lot since then, thanks mainly to encouragement from my mentors at XTBG a whole lot of phylogeneticsinR tutorials (and a whole lot more I haven’t linked here — fear not, a more comprehensive list of tutorials I used is on its way). I’ve also been thinking a good bit about transparency in the scientific process, thanks mainly to that conversation with Derek and the awesome folk at rOpenSci. It’s time to put the thoughts into action.

Although I don’t think that this code I am publishing in this post is very valuable as code, I do think that getting into a habit of writing clean code and sharing it with other scientists is important– at least, it is valuable for me personally. As I mentioned above, I taught myself a whole lot of R by using various tutorials that some wonderful people have developed, and will be pleased if I can contribute to this resource over the coming years. I feel like I am  taking an important first step by writing this post: I have learned how to use knitr, become confident about sharing code I wrote (more or less), and have, in a very very very small way, contributed to the community from which I got so much a few weeks ago.

For some background, this is code I wrote for our Ficus leaf trait project at AFEC, based on a small script originally written by Masatoshi Katabuchi (Thanks!! I hope you don’t mind my putting this up here for some reason…). Incidentally, the presentation based on the AFEC project can be accessed here.

I plan on using this code when I have a phylogeny and want to display some relevant traits of the taxa at the tip – here, the traits are continuous, but this can be easily changed for discrete traits. Lots of people probably have code for this sort of thing already, but if feel free to use this if you wish. If you do use it, please contact me – I’d love to know if what I am doing is useful!

Please excuse the crappy formatting in the code as displayed by WordPress. Don’t know how to make it look prettier.

Here goes.


FicusTree <- read.nexus("pruned_tree.nex")
# pruning the tree to remove the outgroup taxa
FicusTree <- drop.tip(FicusTree, c("Antiaropsis_decipiens", "Castilla_elastica", 
    "Poulsenia_armata", "Sparattosyce_dioca"))

#Unfortunately, this code requires that the trait values are ordered according their order in the phylogeny (i.e. if “Species c” is at the top of the phylogeny, followed by “Species F,” then the first two rows in the data table should be “Species c” and “Species F”). I am trying to write some code that doesn't require this

n.sp <- 30  #number of taxa in your phylogeny

# extract the colums we want to work with as vectors:
sla <- TraitVal$sla
# the current version of the code requires this vector to be reversed.
sla <- rev(sla)
toughness <- TraitVal$toughness
toughness <- rev(toughness)
area <- TraitVal$area
area <- rev(area)

n.traits <- 3  #number of traits you want shown

par(mfrow = c(1, 2))
par(mar = c(0, 0, 2, 0))

plot(FicusTree, show.tip.label = T, cex = 0.75)

plot(rep(1:n.traits, each = n.sp), 
rep(1:n.sp, n.traits), xlim = c(0, n.traits + 
    1), axes = F, ylab = "", xlab = "", type = "n")
abline(h = 1:n.sp, col = "gray75")
abline(v = 1:n.traits, col = "gray75")
axis(1, at = 1:n.traits, labels = c("SLA", "Area", "Toughness"), side = 3, col = "white", 
    cex.axis = 1)
points(rep(1:n.traits, each = n.sp), rep(1:n.sp, n.traits), xlim = c(0, n.traits + 
    0.25), pch = 21, col = "darkgreen", cex = (c(sla/100, area/100, toughness/500)))

Here’s the output:

FicusPlot

This zip file contains some dummy files (a phylogeny, data table, and R code). You will be redirected to a FigShare file — figured I’d try that out now too!

Ok! That’s it. Writing this certainly has been a very useful exercise for me; with any luck, it’ll help someone else too.

By the way, I have a separate .html version of this code that knitr actually generated. As I now understand, WordPress can’t handle uploading html files, so if anyone has read this far and has advice on this, please do let me know. I plan to eventually migrate away from WordPress — that move may come sooner than later now.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s