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