Run the Phylip application 'contrast' from R, and get formatted output.

phylip_fromR.R

##############################################
# Created by Scott Chamberlain
# Ecology and Evolutionary Biology Dept., Rice University
# Houston, TX 77005, USA
# [email protected]
##############################################
#### 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")