Processing nested lists

So perhaps you have all figured this out already, but I was excited to figure out how to finally neatly get all the data frames, lists, vectors, etc. out of a nested list. It is as easy as nesting calls to the apply family of functions, in the case below, using plyr's apply like functions. Take this example:

# Nested lists code, an example
# Make a nested list
mylist <- list()
mylist_ <- list()
for(i in 1:5) {
 for(j in 1:5) {
  mylist[[j]] <- i*j
 }
mylist_[[i]] <- mylist
}
 
# return values from first part of list
laply(mylist_[[1]], identity)
[1] 1 2 3 4 5

 
# return all values
laply(mylist_, function(x) laply(x, identity))
     1  2  3  4  5
[1,] 1  2  3  4  5
[2,] 2  4  6  8 10
[3,] 3  6  9 12 15
[4,] 4  8 12 16 20
[5,] 5 10 15 20 25

 
# perform some function, in this case sqrt of each value
laply(mylist_, function(x) laply(x, function(x) sqrt(x)))

  
        1        2        3        4        5
[1,] 1.000000 1.414214 1.732051 2.000000 2.236068
[2,] 1.414214 2.000000 2.449490 2.828427 3.162278
[3,] 1.732051 2.449490 3.000000 3.464102 3.872983
[4,] 2.000000 2.828427 3.464102 4.000000 4.472136
[5,] 2.236068 3.162278 3.872983 4.472136 5.000000

Created by Pretty R at inside-R.org

Running Phylip's contrast application for trait pairs from R

Here is some code to run Phylip's contrast application from R and get the output within R to easily manipulate yourself. Importantly, the code is written specifically for trait pairs only as the regular expression work in the code specifically grabs data from contast results when only two traits are input. You could easily change the code to do N traits. Note that the p-value calculated for the chi-square statistic is not output from contrast, but is calculated within the function 'PhylipWithinSpContr'. In the code below there are two functions that make a lot of busy work easier: 'WritePhylip' and 'PhylipWithinSpContr'. The first function is nice because the formatting required for data input to Phylip programs is so, well, awkward  - and this function does it for you. The second function runs contrast and retrieves the output data. The example data set I produce in the code below has multiple individuals per species, so that contrasts are calculated taking into account within species variation. Get Phylip's contrast documentation here.

Note that the data input format allows only 10 characters for the species name, so I suggest if your species names are longer than 10 characters use the function abbreviate, for example, to shorten all names to no longer than 10 characters. Also, within the function WritePhylip I concatenate species names and their number of individuals per species leaving plenty of space.

Also, mess around with the options in the "system" call to get what you want. For example, I used "R", "W" and "Y", meaning replace old outfile (R), then turn on within species analyses (W), then accept all options (Y). E..g, if you don't have an old outfile, then you obviously don't need to replace the old file with the "R" command.

(p.s. I have not tried this on a windows machine).

	##############################################
	# Created by Scott Chamberlain
	# Ecology and Evolutionary Biology Dept., Rice University
	# Houston, TX 77005, USA
	# myrmecocystus@gmail.com
	##############################################
	#### Load packages
	require(stringr); require(ape); require(plyr); require(reshape2)
	#### Set directory
	setwd("/Mac/R_stuff/Blog_etc/Phylip_fromR")
	#### Make a phylogeny
	tree <- rcoal(20)
	#### Evolve two traits on the tree
	traits_df <- data.frame(species = rep(tree$tip.label, each = 4),
	xtrait = round(melt(data.frame(replicate(20, rnorm(4, 5, 1))))[,2], 3),
	ytrait = round(melt(data.frame(replicate(20, rnorm(4, 8, 1))))[,2], 3)
	)
	##### Function WritePhylip creates Phylip format data input file
	# x = data frame of 3 rows: species names, first trait, second trait
	WritePhylip <- function(x) {
	tempdat <- x
	names(tempdat)[1] <- "species"
	nsp <- length(unique(tempdat[,1]))
	ntr <- 2
	head <- paste(" ", nsp, ntr, collapse = ' ')
	splist <- list()
	for(i in unique(tempdat[,1])) {
	subdat <- tempdat[tempdat$species == i, ]
	sp <- paste(as.character(subdat[1,1]), length(subdat[,1]), sep = ' ')
	traits <- rbind(sp, cbind(adply(subdat[,-1], 1, paste, collapse = ' ')[,3]))
	splist[[i]] <- traits
	}
	phylipout <- rbind(head, do.call(rbind, splist))
	write.table(phylipout, file = paste(paste(names(tempdat[-1]), collapse="_"), ".txt", sep=""),
	append = F, sep = "\t", col.names = F, row.names = F, quote = F)
	}
	# Write data sets and trees to file
	WritePhylip(traits_df)
	write.tree(tree, "tree.txt")
	#### Run contrast.app, function 'PhylipWithinSpContr'
	PhylipWithinSpContr <- function (data, phylog) {
	exec <- "/Mac/phylip-3.69/exe/contrast"
	system(exec, input = c(data, phylog, "R", "W", "Y"))
	out <- read.table("outfile", skip = 7, sep="\t")
	#### Process data output
	a <- data.frame(sapply(out[1:2,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	b <- data.frame(sapply(out[5:6,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	c <- data.frame(sapply(out[9:10,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	d <- data.frame(sapply(out[13:14,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	e <- data.frame(sapply(out[17:18,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	f <- data.frame(sapply(out[21:22,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	g <- data.frame(sapply(out[26:27,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	h <- data.frame(sapply(out[30:31,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	i <- data.frame(sapply(out[34:35,], function(x) as.numeric(str_split(str_trim(x, "b"), " +")[[1]])))
	logL_withvar_df <- as.numeric(unlist(str_split(str_split(str_trim(out[38:41,], "b")[1], " +")[[1]][6:7], "[,]"))[c(1,3)])
	logL_withoutvar_df <-as.numeric(unlist(str_split(str_split(str_trim(out[38:41,], "b")[2], " +")[[1]][6:7], "[,]"))[c(1,3)])
	chisq_df <- as.numeric(unlist(str_split(str_split(str_trim(out[38:41,], "b")[4], " +")[[1]][4:5], "[,]"))[c(1,3)])
	chisq_p <- pchisq(chisq_df[1], chisq_df[2], lower.tail = F)
	dat <- ldply(list(a, b, c, d, e, f, g, h, i))
	dat <- rbind(dat, logL_withvar_df, logL_withoutvar_df, chisq_df, c(chisq_p, -999))
	names2 <- c(rep(c("VarAIn_VarAest","VarAIn_VarEest","VarAIn_VarAreg","VarAIn_VarAcorr","VarAIn_VarEreg","VarAIn_VarEcorr",
	"VarAOut_VarEest","VarAOut_VarEreg","VarAOut_VarEcorr"), each=2),
	"logL_withvar_df","logL_withoutvar_df","chisq_df","chisq_p")
	dat <- data.frame(names2, dat[,1], dat[,2])
	return(dat)
	}
	#### Run function one trait set
	datout <- PhylipWithinSpContr("xtrait_ytrait.txt", "tree.txt")

view raw phylip_fromR.R hosted with ❤ by GitHub

Here is example output:

> datout
               names2   dat...1.    dat...2.
1      VarAIn_VarAest   0.000110   -0.000017
2      VarAIn_VarAest  -0.000017    0.000155
3      VarAIn_VarEest   0.790783   -0.063097
4      VarAIn_VarEest  -0.063097    0.981216
5      VarAIn_VarAreg   1.000000   -0.107200
6      VarAIn_VarAreg  -0.151800    1.000000
7     VarAIn_VarAcorr   1.000000   -0.127600
8     VarAIn_VarAcorr  -0.127600    1.000000
9      VarAIn_VarEreg   1.000000   -0.064300
10     VarAIn_VarEreg  -0.079800    1.000000
11    VarAIn_VarEcorr   1.000000   -0.071600
12    VarAIn_VarEcorr  -0.071600    1.000000
13    VarAOut_VarEest   0.790734   -0.063104
14    VarAOut_VarEest  -0.063104    0.981169
15    VarAOut_VarEreg   1.000000   -0.064300
16    VarAOut_VarEreg  -0.079800    1.000000
17   VarAOut_VarEcorr   1.000000   -0.071600
18   VarAOut_VarEcorr  -0.071600    1.000000
19    logL_withvar_df -68.779770    6.000000
20 logL_withoutvar_df -68.771450    3.000000
21           chisq_df  -0.016640    3.000000
22            chisq_p   1.000000 -999.000000

Phylometa from R: Randomization via Tip Shuffle

---UPDATE: I am now using code formatting from gist.github, so I replaced the old prettyR code (sorry guys). The github way is much easier and prettier. I hope readers like the change.




I wrote earlier about some code I wrote for running Phylometa (software to do phylogenetic meta-analysis) from R.

I have been concerned about what exactly is the right penalty for including phylogeny in a meta-analysis. E.g.: AIC is calculated from Q in Phylometa, and Q increases with tree size.

So, I wrote some code to shuffle the tips of your tree N number of times, run Phylometa, and extract just the "Phylogenetic MA" part of the output. So, we compare the observed output (without tip shuffling) to the distribution of the tip shuffled output, and we can calculate a P-value from that. The code I wrote simply extracts the pooled effect size for fixed and also random-effects models. But you could change the code to extract whatever you like for the randomization.

I think the point of this code is not to estimate your pooled effects, etc., but may be an alternative way to compare traditional to phylogenetic MA where hopefully simply incorporating a tree is not penalizing the meta-analysis so much that you will always accept the traditional MA as better.

Get the code here, and also below. Get the example tree file and data file, named "phylogeny.txt" and "metadata_2g.txt", respectively below (or use your own data!). You need the file "phylometa_fxn.r" from my website, get here, but just call it using source as seen below.

	#####################################################################
	# 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)

view raw phylometa_rand_fxn_one.R hosted with ❤ by GitHub

As you can see, the observed values fall well within the distribution of values obtained from shuffling tips.  P-values were 0.64 and 0.68 for fixed- and random-effects MA's, respectively. This suggests, to me at least, that the traditional (distribution of tip shuffled analyses, the histograms below) and phylogenetic (red lines) MA's are not different. The way I would use this is as an additional analysis to the actual Phylometa output.

RStudio Beta 2 is Out!

RStudio Beta 2 (v0.93) « RStudio Blog

A new beta version of RStudio is out!

Adjust branch lengths with node ages: comparison of two methods

Here is an approach for comparing two methods of adjusting branch lengths on trees: bladj in the program Phylocom and a fxn written by Gene Hunt at the Smithsonian.

Get the code and example files (tree and node ages) here
Get phylocom here: http://www.phylodiversity.net/phylocom/

Gene Hunt's method has many options you can mess with, including setting tip ages (not available in bladj), setting node ages, and minimum branch length imposed. You will notice that Gene's method may be not the appropriate if you only have extant taxa.

The function AdjBrLens uses as input a newick tree file and a text file of node ages, and uses functions you can simply run by "source" the R file bladjing_twomethods.R file from here.

Note that blad does not like numbers for node names, so you have to put a character in front of a number of just character names for nodes.

# This is where the work happens... 
# Directory below needs to have at least three items:
#  1. phylocom executable for windows or mac
#  2. tree newick file
#  3. node ages file as required by phylocom, see their manual
# Output: trees_out is a list of three trees, the original, bladj, and Gene Hunt's method
# Also, within the function all three trees are written to file as PDFs
setwd("/Mac/R_stuff/Blog_etc/Bladjing") # set working directory
source("bladjing_twomethods.R") # run functions from source file
trees_out <- AdjBrLens("tree.txt", "nodeages.txt")
 
# plot trees of three methods together, 
# with nodes with age estimates labeled
jpeg("threeplots.jpeg", quality=100)
layout(matrix(1:3, 1, 3))
plot(trees_out[[1]])
nodelabels(trees_out[[1]]$node.label, cex = 0.6)
title("original tree")
plot(trees_out[[2]])
nodelabels(trees_out[[2]]$node.label, cex = 0.6)
title("bladj method")
plot(trees_out[[3]])
nodelabels(trees_out[[3]]$node.label, cex = 0.6)
title("gene hunt method, scalePhylo")
dev.off()

Created by Pretty R at inside-R.org


It is sort of hard to see the differences in the branch length changes here, but the individual output trees will reveal the differences better.

Phylometa from R - UDPATE

THIS BLOG HAS MOVED: New version of this post at http://recology.info/2011/04/phylometa-from-r-udpate/

A while back I posted some messy code to run Phylometa from R, especially useful for processing the output data from Phylometa which is not easily done. The code is still quite messy, but it should work now. I have run the code with tens of different data sets and phylogenies so it should work.

I fixed errors when parentheses came up against numbers in the output, and other things. You can use the code for up to 4 levels of your grouping variable. In addition, there are some lines of code to plot the effect sizes with confidence intervals, comparing random and fixed effects models and phylogenetic and traditional models. 

Get the code at my website:

-phylometa_fxn.R

-phylometa_inR_run.R

- Use the first file to do the work, calling the second file using source().
- This new code works with Marc's new version of Phylometa, so please update: http://lajeunesse.myweb.usf.edu/publications.html

Again, please let me know if it doesn't work, if it's worthless, what changes could make it better.

Some notes on tree formatting for Phylometa.
1. Trees cannot have node labels - remove them (e.g., tree$node.label < NULL).
2. Trees cannot have zero length branches. This may seem like a non-problem, but it might be for example if you have resolved polytomies and zero length branches are added to resolve the polytomy.
3. I think you cannot have a branch length on the root branch.