Conduct tip shuffle randomization of phylogenetic meta-analysis

phylometa_rand_fxn_one.R

#####################################################################
# Created by Scott Chamberlain										
# Ecology and Evolutionary Biology Dept., Rice University 			
# Houston, TX 77005, USA												
# myrmecocystus@gmail.com		
# http://schamber.wordpress.com/									
#####################################################################

##### Directions
#Put your phylogeny, in format required by phylometa (see Phylometa manual), in your working directory
#Put your meta-analysis dataset, in format required by phylometa (see Phylometa manual), in your working directory
#Set working directory
#Use below functions
#Beware: only use a moderator variable with up to 4 groups 

# Install packages
install.packages(c("plyr","ggplot2","ape","stringr")) # if not already installed

# Load packages
require(ggplot2); require(ape); require(stringr)

# Call and run functions (used below) in the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY]
source("Z:/Mac/R_stuff/Blog_etc/Phylometa_randomization/phylometa_fxn.r")

####### Fxn to do randomization equivalent of Phylometa analysis
# phy = name of phylogeny file
# dat = name of data file in Phylometa format
# reps = number of randomized trees per real tree
# groups = number of levels in grouping/moderator variable
# Note that trees are written to directory so that Phylometa can call them
# Output is a jpeg of the histograms of randomized results, and vertical line for observed value
# 	Also, p-values for fixed and random-effects are output in the console
PhylogMetaRand <- function(phy, dat, reps, groups) {

	# Observed data
	obsout <- phylometa.process(system(paste("\"new_phyloMeta_1.2b.exe\"", phy, dat), intern=T), groups)
	subobsout <- subset(obsout[4][[1]], Group == "All studies")
	obsfixedeffect <- subobsout[1,3]
	obsrandeffect <- subobsout[2,3]

	# Randomization of data by shuffling tips on tree
  	  	phylo <- read.tree(phy)
  	  	rtree_list <- list()
    	for(i in 1:reps) {
        	temp <- phylo
        	randtips <- sample(temp$tip.label, length(temp$tip.label))
        	temp$tip.label <- randtips
        	rtree_list[[i]] <- temp
    	}
    	
		treenames <- list()
    	for(i in 1:reps){
        	treenames[[i]] <- paste("tree",i,".txt","",sep="")
    	}
    	WriteTrees2(rtree_list)
    
    	fixedrandeffects <- data.frame(matrix(NA, reps, 2))
    	for(i in 1:length(treenames)) {
        	cat(i, "\n")
        	out <- system(paste("\"new_phyloMeta_1.2b.exe\"", treenames[[i]], dat), intern=T)
        	fixrandout <- subset(phylometa.process(out, groups)[4][[1]], Group == "All studies")
        	fixedrandeffects[i,] <- c(fixrandout[1,3], fixrandout[2,3])
	    }
    	names(fixedrandeffects) <- c("fixed", "random")

	# Plot results
	fixed_ <- ggplot(fixedrandeffects, aes(x = fixed)) + 
    	geom_histogram() +
    	geom_vline(xintercept = obsfixedeffect, colour="red", size=1)
	random_ <- ggplot(fixedrandeffects, aes(x = random)) + 
    	geom_histogram() +
	    geom_vline(xintercept = obsrandeffect, colour="red", size=1)
	jpeg(paste(dat, ".jpeg", sep=""))
	arrange_(fixed_, random_)
	dev.off()

	# Calculate randomization P-value
	num_greater_than_fixed <- length(which(fixedrandeffects$fixed >= obsfixedeffect))
	pval_fixed <- num_greater_than_fixed / reps
	num_greater_than_random <- length(which(fixedrandeffects$random >= obsrandeffect))
	pval_random <- num_greater_than_random / reps
	pvals <- c(pval_fixed, pval_random)
	names(pvals) <- c("pval_fixed","pval_random")
	pvals
}

# E.g.
setwd("Z:/Mac/R_stuff/Blog_etc/Phylometa_randomization/")
PhylogMetaRand(phy = "phylogeny.txt", dat = "metadata_2g.txt", reps = 100, groups = 2)