Phylogenetic meta-analysis in R using Phylometa
[This article was first published on Recology, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Here is some code to run Phylometa from R. Phylometa is a program that conducts phylogenetic meta-analyses. The great advantage of the approach below is that you can easily run Phylometa from R, and manipulate the output from Phylometa in R. Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Phylometa was created by Marc Lajeunesse at University of South Florida, and is described in his 2009 AmNat paper. Phylometa can be downloaded free here.
Save phylometa_fxn.R to your working directory. Then use the block of code below to call the functions within phylometa_fxn.R. Also, you can download a zip file from my website that has the program, the R code files, and example data files.
The program Phylometa needs to be in the working directory you are calling from.
Let me know what doesn’t work, and what improvements can be made; I’m sure there are many!
########Directions #Place phylometa software to your working directory #Put your phylogeny, in format required by phylometa, in your working directory #Put your meta-analysis dataset, in format required by phylometa, in your working directory #Set working directory #Use below functions #Beware: only use a moderator variable with up to 6 groups ########Install packages install.packages(c("plyr","ggplot2")) library(plyr) library(ggplot2) ########Set the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY] setwd("/Users/Scott/Documents/phylometa") #Call and run functions (used below) in the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY] source("/Users/Scott/Documents/phylometa") ###########################Functions to to a phylogenetic meta-analysis #Define number of groups in moderator variable groups <- 2 ####Run phylometa. Change file names as needed phylometa.run <- system(paste('"phyloMeta_v1-2_beta.exe" phylogeny.txt metadata_2g.txt'),intern=T) ####Process phylometa output #E.g. myoutput <- phylometa.process(phylometa.run,groups) ####Get output from phylometa.run phylometa.output(myoutput) #Prints all five tables phylometa.output.table(myoutput,2) #Prints the table you specify, from 1 to 5, in this example, table 2 is output ################################################### #########Plot effect sizes. These are various ways to look at the data. Go through them to see what they do. Output pdf's are in your working directory #Make table for plotting analysis <- c(rep("fixed",groups+1),rep("random",groups+1)) trad_effsizes <- data.frame(analysis,phylometa.output.table(myoutput,2)) #Tradiational effect size table phylog_effsizes <- data.frame(analysis,phylometa.output.table(myoutput,4)) #Phylogenetic effect size table #The arrange method limits <- aes(ymax = effsize + (CI_high-effsize), ymin = effsize - (effsize-CI_low)) dodge <- position_dodge(width=0.3) plot01 <- ggplot(trad_effsizes,aes(y=effsize,x=analysis,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank(),title="Traditional meta-analysis") + labs(x="Group",y="Effect size") + geom_errorbar(limits, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) plot02 <- ggplot(phylog_effsizes,aes(y=effsize,x=analysis,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank(),title="Phylogenetic meta-analysis") + labs(x="Group",y="Effect size") + geom_errorbar(limits, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) pdf("plots_effsizes_arrange.pdf",width = 8, height = 11) arrange(plot01,plot02,ncol=1) dev.off() #used in the two plotting methods below bothanalyses<-data.frame(tradphy=c(rep("Traditional",(groups*2)+2),rep("Phylogenetic",(groups*2)+2)),fixrand=rep(analysis,2),rbind.fill(phylometa.output.table(myoutput,2),phylometa.output.table(myoutput,4))) #Table of both trad and phylo limits2 <- aes(ymax = effsize + (CI_high-effsize), ymin = effsize - (effsize-CI_low)) dodge <- position_dodge(width=0.3) #The grid/wrap method, version 1 plot03 <- ggplot(bothanalyses,aes(y=effsize,x=tradphy,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank()) + labs(x="Group",y="Effect size") + geom_errorbar(limits2, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) + facet_grid(.~fixrand) pdf("plots_effsizes_wrap1.pdf") plot03 dev.off() #The grid/wrap method, version 2 (excuse the sloppy x-axis labels) plot04 <- ggplot(bothanalyses,aes(y=effsize,x=Group,colour=tradphy)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank()) + labs(x="Group",y="Effect size") + geom_errorbar(limits2, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) + facet_grid(.~fixrand) pdf("plots_effsizes_wrap2.pdf") plot04 dev.off()
Below is an example output figure from the code. This example is from an analysis using 5 groups (i.e., 5 levels in the explanatory variable).
To leave a comment for the author, please follow the link and comment on their blog: Recology.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.