R/cohorts.R

In BAMMtools: Analysis and Visualization of Macroevolutionary Dynamics on Phylogenetic Trees

##' @title Visualize macroevolutionary cohorts
##'
##' @description Plots the matrix of pairwise correlations in rate regimes
##'     between all tips in a phylogeny.
##'
##' @param x A matrix of pairwise correlations generated by
##'     \code{getCohortMatrix}.
##' @param ephy An object of class \code{bammdata}.
##' @param col A vector of colors passed to the function \code{image}. These
##'     will be used to color the values in \code{x}. See documentation for
##'     \code{image}. If \code{col = 'temperature'}, the color palette from
##'		\code{\link{rich.colors}} from the gplots package will be used. 
##' @param pal The palette to use if \code{use.plot.bammdata=TRUE}. See
##'     options documented in the help file for \code{\link{plot.bammdata}}.  	
##' @param lwd A numeric indicating the width of branches in the phylogeny.  	
##' @param ofs A numeric controlling the offset of the phylogeny from the
##'     matrix plot. Appropriate values will probably be in the interval
##'     [0,0.1].
##' @param use.plot.bammdata Logical. should a phylorate plot be generated?
##' @param useraster A logical indicating whether the function \code{image}
##'     should plot the matrix as a raster.
##' @param LARGE An integer. If trees have more tips than \code{LARGE},
##'     \code{useraster} will be coerced to \code{TRUE}.
##' @param \dots Further arguments passed to \code{plot.bammdata} if
##'     \code{use.plot.bammdata=TRUE}.
##'
##' @details The plotting function creates an image of the \code{BAMM}
##'     correlation matrix between tip lineages of the phylogeny. Each
##'     correlation is the posterior frequency with which a pair of lineages
##'     occurs in the same macroevolutionary rate regime. Correlations are
##'     mapped to a set of colors, with warmer colors corresponding to higher
##'     correlations than cooler colors. The set of colors is specified by the
##'     \code{col} argument and a legend is plotted to guide interpretation of
##'     the color-correlation map. Trees are plotted on the margins of the
##'     matrix image. The correlation between any two tips can be inferred by
##'     finding their intersection within the matrix image.
##'
##'     \strong{IMPORTANT}: the legend DOES NOT apply to the phylorate plots
##'     shown in the margin if \code{use.plot.bammdata=TRUE}.
##'
##' @author Mike Grundler
##'
##' @seealso \code{\link{plot.bammdata}}, \code{\link{getCohortMatrix}},
##'     \code{\link{image}}
##'
##' @examples
##' data(whales, events.whales)
##' ed <- getEventData(whales, events.whales, burnin=0.1, nsamples=500)
##' x <- getCohortMatrix(ed)
##' cohorts(x, ed)
##' cohorts(x, ed, col='temperature')
##' cohorts(x, ed, ofs=0.05, col='temperature')
##' cohorts(x, ed, pal="temperature", col='temperature', use.plot.bammdata=TRUE)
##' # gray scale
##' cohorts(x, ed, col=gray(seq(0.2,0.9,length.out=128)),
##'         use.plot.bammdata=FALSE)
##' @export
cohorts <- function(x, ephy, col, pal, lwd = 1, ofs = 0, use.plot.bammdata = FALSE, useraster = FALSE, LARGE = 500,...) {
	
	if (is.null(dimnames(x) ))
		stop("x must have row and column names");
		
	op <- par(no.readonly = TRUE);
	figs <- matrix(c(0,0.2,0.8,1,
	                 0.2,0.95,0.8+ofs,1,
	                 0,0.2-ofs,0,0.8,
	                 0.2,0.95,0,0.8,
	                 0.95,1,0.25,0.75
	                 ), byrow=TRUE,
	               nrow=5, ncol=4);
	if (dim(x)[1] > LARGE)
		useraster <- TRUE;
	if (missing(pal))
		pal <- "RdYlBu";
	if (missing(col))
		col <- colorRampPalette(get("palettes",.colorEnv)[["RdYlBu"]])(64);
	if (all(col == 'temperature'))
		col <- gplots::rich.colors(64);
	ncolors <- length(col);
	breaks <- quantile(seq(0,1.01,length.out=100),probs=seq(0,1,length.out=ncolors+1));
	
	index <- match(ephy$tip.label, rownames(x));
	x <- x[index, index];
	
	if (use.plot.bammdata) {               
		par(fig = figs[2,], new=FALSE, mar = c(0,0,1,4));
		plot(ephy, pal=pal,lwd=lwd,direction="downwards",...);
		par(fig = figs[3,], new=TRUE, mar = c(5,1,0,0));
		plot(ephy,pal=pal,lwd=lwd,direction="rightwards",...)
		par(fig = figs[4,], new=TRUE, mar = c(5,0,0,4));
		plot(0,0,type="n",axes=FALSE,ann=FALSE,xlim=c(0,1),ylim=c(0,1))
		image(x,axes=FALSE,xlab="",ylab="",col=col,xlim=c(0,1),ylim=c(0,1),breaks=breaks,add=TRUE,useRaster=useraster);
	}
	else {
		phy <- as.phylo.bammdata(ephy);
		bt <- max(ephy$end)
		par(fig = figs[2,], new=FALSE, mar = c(0,0,1,4));
		plot.phylo(phy,edge.width=lwd,direction="downwards",show.tip.label=FALSE,x.lim=c(1,length(phy$tip.label)),y.lim=c(0,bt));
		par(fig = figs[3,], new=TRUE, mar = c(5,1,0,0));
		plot.phylo(phy,edge.width=lwd,direction="rightwards",show.tip.label=FALSE,y.lim=c(1,length(phy$tip.label)),x.lim=c(0,bt));
		par(fig = figs[4,], new=TRUE, mar = c(5,0,0,4));
		gl <- 1:(length(ephy$tip.label)+1);
		plot(0,0,type="n",axes=FALSE,ann=FALSE,xlim=c(1,length(gl)-1),ylim=c(1,length(gl)-1))
		image(gl,gl,x,axes=FALSE,xlab="",ylab="",col=col,xlim=c(1,length(gl)-1),ylim=c(1,length(gl)-1),breaks=breaks,add=TRUE,useRaster=useraster);
	}
	#barLegend(col, quantile(seq(min(x),max(x),length.out=ncolors+1),probs=seq(min(x),max(x),length.out=ncolors+1)),fig=figs[5,],side=2);
	barLegend(col,breaks,fig=figs[5,],side=2);
	par(op);
}