Just hit ESC. You're welcome.


*A situation where you might have this problem:
I have a treefile with 5000 trees in it.
I import this file into R using the package ape, so the class for this object is multiPhylo[5000].
I enter the command "plot(trees)," like an idiot.
RStudio plots each tree, one at a time, with the text output "Hit <Return> to see next plot:"

I tried "^c" and a bunch of other things. My office mate suggested that I set my coffee cup down on the return key until it stops.

Turns out I just had to hit ESC. Duh!

I'm currently working with a Sanger-generated 10-gene dataset, which includes a few fast-evolving genes. During the initial data exploration phase I used saturation plots to check for potential decay of phylogenetic signal caused by multiple substitutions (saturated nucleotide variation).

There are a lot of ways to visualize saturation. One method, which I believe was originally used by Philippe et al. (1994), is to plot the raw or uncorrected pairwise genetic distances in an alignment against model-corrected genetic distances. If the relationship is approximately linear, then the gene is not saturated; if the line curves or plateaus, there is evidence of saturation. 

Here is an example of an unsaturated gene:

And an example of a saturated gene:

This is a really rough method, which should probably only be used as a preliminary exploration of your data. As far as I know, there is not an established slope value that says definitively, "yes, this gene is saturated." However, I do think it's a useful thing to look at, and it's really easy to do in R. You may want to look into APE's dist.dna command for all of the available models. Here is the R-code I used to make these simple plots:

library(ape)

###Input data: a phylip-format alignment file, converted to a 'DNAbin' object###
dat<-read.dna(file="myData.phy", format = "sequential", as.character=TRUE, skip=0)
dat<-as.DNAbin(dat)

###Convert to genetic distances###
dist<-dist.dna(dat, model="raw")
dist.corrected<-dist.dna(dat, model="TN93")

###Make plot###
plot(dist~dist.corrected, pch=20, col="red", xlab="TrN model distance", ylab="Uncorrected genetic distance", main="Saturation Plot")
abline(0,1, lty=2)
abline(lm(dist~dist.corrected), lwd=3)
lm_coef<-coef(lm(dist~dist.corrected))
text(0.1,0.05,bquote(y == .(lm_coef[2])*x))

Here's a quick way to use R to visualize treespace! It's not necessarily the most rigorous analysis, but I think it's a useful tool for visualizing what's going on between your gene and species trees after an analysis like *BEAST.

Basically I'm computing a matrix of Robinson-Foulds (RF) distances between the trees produced in *BEAST, then plotting these using a Principal Coordinates Analysis (pcoa). Apparently the same analysis can also be implemented in the package TreeSetViz within Mesquite, although I wasn't able to get that to work :( 

I found the tree.dist.matrix function on Dan Warren's github page. Thanks, Dan! Here's the rest of the code I used:

###
library(ape)
library(phangorn)

# import a single nexus file with all of the trees you want to plot. In order to make the colors
# work, each tree name should start with "geneName" and an underscore. Here's an example.
All.trees<-read.nexus(file="/Users/kathryneverson/Desktop/RTreeSpace.trees")

tree.dist.matrix <- function(trees, treenames=names(trees)){
  N <- length(trees)
  
  if(N != length(treenames)){
    stop("Names and tree list must be the same length")
  }
  
  #Create an empty matrix for results
  RF <- matrix(0, N, N)
  
  for(i in 1:(N-1)){
    #print(paste("Tree", i, "of", N))        
    for(j in (i+1):N){
      RFd <- RF.dist(trees[[i]],trees[[j]])
      if(RFd==0) RFd = 0.000000001
      RF[i,j]<-RF[j,i]<-RFd
    }
  }
  
  #Row and column names
  rownames(RF) <- treenames
  colnames(RF) <- treenames
  
  RF
}

All.dist<-tree.dist.matrix(All.trees)
gene <- factor(gsub("_.*", "", rownames(All.dist), perl=TRUE))
levels(gene)
#there should be as many colors in this list as there are genes
cols = c('darkmagenta', 'azure4', 'blue', 'green', 'darkgreen', 'orange2', 'cyan', 'yellow2', 'red', 'black', 'darkred', 'hotpink')[gene]

All.pcoa<-pcoa(All.dist)
biplot(All.pcoa)
plot(All.pcoa$vectors[,1], All.pcoa$vectors[,2], type = "p", xlab = "", ylab = "",
     pch=20, axes = FALSE, main = "Visualization of Tree Space (MDS)", col=cols)
legend('top', legend=levels(gene), cex = 0.5, pch = 20, col =c('darkmagenta', 'azure4', 'blue', 'green', 'darkgreen', 'orange2', 'cyan', 'yellow2', 'red', 'black', 'darkred', 'hotpink'))

###

NOTE: Calculating the RF distances is pretty slow, and I would not recommend trying to include thousands of trees. Here I thinned out my *BEAST tree files and just plotted 100 trees for each gene + the species tree.

I couldn't wait to share a fun new R command I learned last week at the Bodega Applied Phylogenetics Workshop: geoscalePhylo. This command plots a gorgeous chronogram with a colorful geological timescale. I usually format my trees manually in Adobe Illustrator, but this does a really good job.

You'll need two packages: phyloch (if you're reading in a BEAST tree) and strap. Note that phyloch is not available in CRAN and must be downloaded from the author's website. Here's the R code for making this pretty tree:

library(phyloch)
library(strap)
library(OutbreakTools)

t <- read.annotated.nexus("myBeastTree.tre")

t$root.time <- t$root.annotation$height

geoscalePhylo(tree=ladderize(t,right=FALSE), units=c("Period", "Epoch", "Age"), boxes="Epoch",   +cex.tip=0.5, cex.age=0.7, cex.ts=0.7, label.offset=0, x.lim=c(-15,80), lwd=3, width=2)
###

This command is really flexible. "units" allows you to select whether or not to show Periods, Epochs, Ages, Eons, or a user-defined timescale. "boxes" allows you to pick which unit to plot gray boxes on. "width" controls the thickness of the branches. "x.lim" controls the year span. Yes, I know that the labels are squished on the time axis, but you can adjust those.

Compare this to a tree made with the simple plot command:

plot(t, show.tip.label=F)
axisPhylo()

geoscalePhylo is pretty good, right? Special thanks to Tracy Heath for showing us this command.

Hi all, I've run into the following problem a few times now (most recently with the species tree program STAR): Some programs require input trees that are strictly bifurcating. In other words, they are not allowed to have any polytomies/multifurcations/multichotomies. 

I wish my trees were always strictly bifurcating! But sometimes that's just not the case.

Here's a simple workaround in R using the package APE. The crucial command is multi2di, which transforms all polytomies into a series of dichotomies with one or more branches of length zero. Here's my R code:

###Code for a Phylo Object###
library(ape)

myTree<-read.tree("myTree.tre")

is.binary.tree(myTree)
[1] FALSE

myTree.bifurcating<-multi2di(myTree)

is.binary.tree(myTree.bifurcating)
[1] TRUE

write.tree(myTree.bifurcating, file="myTree.bifurcating.tre")

###Code for a Multiphylo Object###
library(ape)

myTrees<-read.tree("myTrees.tre")

is.binary.tree(myTrees[[1]])
[1] FALSE

make.bifurcating<-function(phy){
  for(i in 1:length(phy)){
    phy[[i]]<-multi2di(phy[[i]])
  }
  return(phy)
}

myTrees.bifurcating<-make.bifurcating(myTrees)

is.binary.tree(myTrees.bifurcating[[1]])
[1] TRUE

write.tree(myTree.bifurcating, file="myTrees.bifurcating.trees")

NOTE: I would NOT recommend using poorly resolved trees with multiple polytomies as an input into any program that assumes or depends on the accuracy of gene trees. That said, if you have lots of gene trees with just one or two small polytomies, this method works pretty well.