R/heatmap.2.R

In gplots: Various R Programming Tools for Plotting Data

## $Id: heatmap.2.R 2132 2016-08-12 17:17:24Z warnes $

heatmap.2 <- function (x,

                       ## dendrogram control
                       Rowv = TRUE,
                       Colv=if(symm)"Rowv" else TRUE,
                       distfun = dist,
                       hclustfun = hclust,
                       dendrogram = c("both","row","column","none"),
                       reorderfun = function(d, w) reorder(d, w),
                       symm = FALSE,

                       ## data scaling
                       scale = c("none","row", "column"),
                       na.rm=TRUE,

                       ## image plot
                       revC = identical(Colv, "Rowv"),
                       add.expr,

                       ## mapping data to colors
                       breaks,
                       symbreaks=any(x < 0, na.rm=TRUE) || scale!="none",

                       ## colors
                       col="heat.colors",

                       ## block sepration
                       colsep,
                       rowsep,
                       sepcolor="white",
                       sepwidth=c(0.05,0.05),

                       ## cell labeling
                       cellnote,
                       notecex=1.0,
                       notecol="cyan",
                       na.color=par("bg"),

                       ## level trace
                       trace=c("column","row","both","none"),
                       tracecol="cyan",
                       hline=median(breaks),
                       vline=median(breaks),
                       linecol=tracecol,

                       ## Row/Column Labeling
                       margins = c(5, 5),
                       ColSideColors,
                       RowSideColors,
                       cexRow = 0.2 + 1/log10(nr),
                       cexCol = 0.2 + 1/log10(nc),
                       labRow = NULL,
                       labCol = NULL,
                       srtRow = NULL,
                       srtCol = NULL,
                       adjRow = c(0,NA),
                       adjCol = c(NA,0),
                       offsetRow = 0.5,
                       offsetCol = 0.5,
                       colRow = NULL,
                       colCol = NULL,

                       ## color key + density info
                       key = TRUE,
                       keysize = 1.5,
                       keygrid=255,
                       density.info=c("histogram","density","none"),
                       denscol=tracecol,
                       symkey = any(x < 0, na.rm=TRUE) || symbreaks,
                       densadj = 0.25,
                       key.title = NULL,
                       key.xlab = NULL,
                       key.ylab = NULL,
                       key.xtickfun = NULL,
                       key.ytickfun = NULL,
                       key.par=list(),

                       ## plot labels
                       main = NULL,
                       xlab = NULL,
                       ylab = NULL,

                       ## plot layout
                       lmat = NULL,
                       lhei = NULL,
                       lwid = NULL,

                       ## extras
                       extrafun=NULL,
                       ...
                       )
{
  scale01 <- function(x, low=min(x), high=max(x) )
    {
      x <- (x-low)/(high - low)
      x
    }

  retval <- list()

  scale <- if(symm && missing(scale)) "none" else match.arg(scale)
  dendrogram <- match.arg(dendrogram)
  trace <- match.arg(trace)
  density.info <- match.arg(density.info)

  if(length(col)==1 && is.character(col) )
    col <- get(col, mode="function")

  if(!missing(breaks) && any(duplicated(breaks)) )
      stop("breaks may not contian duplicate values")

  if(!missing(breaks) && (scale!="none"))
    warning("Using scale=\"row\" or scale=\"column\" when breaks are",
            "specified can produce unpredictable results.",
            "Please consider using only one or the other.")

  if ( is.null(Rowv) || is.na(Rowv) )
    Rowv <- FALSE
  if ( is.null(Colv) || is.na(Colv) )
    Colv <- FALSE
  else if( all(Colv=="Rowv") )
    Colv <- Rowv


  if(length(di <- dim(x)) != 2 || !is.numeric(x))
    stop("`x' must be a numeric matrix")

  nr <- di[1]
  nc <- di[2]

  if(nr <= 1 || nc <= 1)
    stop("`x' must have at least 2 rows and 2 columns")

  if(!is.numeric(margins) || length(margins) != 2)
    stop("`margins' must be a numeric vector of length 2")

  if(missing(cellnote))
    cellnote <- matrix("", ncol=ncol(x), nrow=nrow(x))

  if(!inherits(Rowv, "dendrogram")) {
    ## Check if Rowv and dendrogram arguments are consistent
    if (
          (
             ( is.logical(Rowv) && !isTRUE(Rowv) )
              ||
             ( is.null(Rowv) )
          )
          &&
          ( dendrogram %in% c("both","row") )
       )
      {
        warning("Discrepancy: Rowv is FALSE, while dendrogram is `",
                dendrogram, "'. Omitting row dendogram.")

        if (dendrogram=="both")
          dendrogram <- "column"
        else
          dendrogram <- "none"

      }
  }

  if(!inherits(Colv, "dendrogram")) {
    ## Check if Colv and dendrogram arguments are consistent
    if (
         (
           (is.logical(Colv) && !isTRUE(Colv) )
             ||
            (is.null(Colv))
         )
         &&
         ( dendrogram %in% c("both","column")) )
      {
        warning("Discrepancy: Colv is FALSE, while dendrogram is `",
                dendrogram, "'. Omitting column dendogram.")

        if (dendrogram=="both")
          dendrogram <- "row"
        else
          dendrogram <- "none"

      }
  }


  ## by default order by row/col mean
  ## if(is.null(Rowv)) Rowv <- rowMeans(x, na.rm = na.rm)
  ## if(is.null(Colv)) Colv <- colMeans(x, na.rm = na.rm)

  ## get the dendrograms and reordering indices

  ## if( dendrogram %in% c("both","row") )
  ##  { ## dendrogram option is used *only* for display purposes
  if(inherits(Rowv, "dendrogram"))
    {
      ddr <- Rowv ## use Rowv 'as-is', when it is dendrogram
      rowInd <- order.dendrogram(ddr)
      if(length(rowInd)>nr || any(rowInd<1 | rowInd > nr ))
         stop("Rowv dendrogram doesn't match size of x")
      if (length(rowInd) < nr)
          nr <- length(rowInd)
    }
  else if (is.integer(Rowv))
    {
      ## Compute dendrogram and do reordering based on given vector
      distr <- distfun(x)
      hcr <- hclustfun(distr)
      ddr <- as.dendrogram(hcr)
      ddr <- reorderfun(ddr, Rowv)

      rowInd <- order.dendrogram(ddr)
      if(nr != length(rowInd))
        stop("row dendrogram ordering gave index of wrong length")
    }
  else if (isTRUE(Rowv))
    { ## If TRUE, compute dendrogram and do reordering based on rowMeans
      Rowv <- rowMeans(x, na.rm = na.rm)
      distr <- distfun(x)
      hcr <- hclustfun(distr)
      ddr <- as.dendrogram(hcr)
      ddr <- reorderfun(ddr, Rowv)

      rowInd <- order.dendrogram(ddr)
      if(nr != length(rowInd))
        stop("row dendrogram ordering gave index of wrong length")
    }
  else if(!isTRUE(Rowv))
    {
      rowInd <- nr:1
      ddr <- as.dendrogram(hclust(dist(diag(nr))))
    }
  else
    {
      rowInd <- nr:1
      ddr <- as.dendrogram(Rowv)
    }

  if(inherits(Colv, "dendrogram"))
    {
      ddc <- Colv ## use Colv 'as-is', when it is dendrogram
      colInd <- order.dendrogram(ddc)
      if(length(colInd)>nc || any(colInd<1 | colInd > nc ))
         stop("Colv dendrogram doesn't match size of x")
      if (length(colInd) < nc)
          nc <- length(colInd)
    }
  else if(identical(Colv, "Rowv")) {
    if(nr != nc)
      stop('Colv = "Rowv" but nrow(x) != ncol(x)')
    if(exists("ddr"))
      {
        ddc <- ddr
        colInd <- order.dendrogram(ddc)
      } else
    colInd <- rowInd
  } else if(is.integer(Colv))
    {## Compute dendrogram and do reordering based on given vector
      distc <- distfun(if(symm)x else t(x))
      hcc <- hclustfun(distc)
      ddc <- as.dendrogram(hcc)
      ddc <- reorderfun(ddc, Colv)

      colInd <- order.dendrogram(ddc)
      if(nc != length(colInd))
        stop("column dendrogram ordering gave index of wrong length")
    }
  else if (isTRUE(Colv))
    {## If TRUE, compute dendrogram and do reordering based on rowMeans
      Colv <- colMeans(x, na.rm = na.rm)
      distc <- distfun(if(symm)x else t(x))
      hcc <- hclustfun(distc)
      ddc <- as.dendrogram(hcc)
      ddc <- reorderfun(ddc, Colv)

      colInd <- order.dendrogram(ddc)
      if(nc != length(colInd))
        stop("column dendrogram ordering gave index of wrong length")
    }
  else if(!isTRUE(Colv))
    {
      colInd <- 1:nc
      ddc <- as.dendrogram(hclust(dist(diag(nc))))
    }
  else
    {
      colInd <- 1:nc
      ddc <- as.dendrogram(Colv)
    }

  retval$rowInd <- rowInd
  retval$colInd <- colInd
  retval$call <- match.call()


  ## reorder x & cellnote
  x <- x[rowInd, colInd]
  x.unscaled <- x
  cellnote <- cellnote[rowInd, colInd]

  if(is.null(labRow))
    labRow <- if(is.null(rownames(x))) (1:nr)[rowInd] else rownames(x)
  else
    labRow <- labRow[rowInd]

  if(is.null(labCol))
    labCol <- if(is.null(colnames(x))) (1:nc)[colInd] else colnames(x)
  else
    labCol <- labCol[colInd]

  if(!is.null(colRow))
    colRow <- colRow[rowInd]

  if(!is.null(colCol))
    colCol <- colCol[colInd]

  if(scale == "row") {
    retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
    x <- sweep(x, 1, rm)
    retval$rowSDs <-  sx <- apply(x, 1, sd, na.rm = na.rm)
    x <- sweep(x, 1, sx, "/")
  }
  else if(scale == "column") {
    retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
    x <- sweep(x, 2, rm)
    retval$colSDs <-  sx <- apply(x, 2, sd, na.rm = na.rm)
    x <- sweep(x, 2, sx, "/")
  }

  ## Set up breaks and force values outside the range into the endmost bins
  if(missing(breaks) || is.null(breaks) || length(breaks)<1 )
    {
      if( missing(col) ||  is.function(col) )
        breaks <- 16
      else
        breaks <- length(col)+1
    }

  if(length(breaks)==1)
    {
      if(!symbreaks)
        breaks <- seq( min(x, na.rm=na.rm), max(x,na.rm=na.rm), length=breaks)
      else
        {
          extreme <- max(abs(x), na.rm=TRUE)
          breaks <- seq( -extreme, extreme, length=breaks )
        }
    }

  nbr <- length(breaks)
  ncol <- length(breaks)-1

  if(class(col)=="function")
    col <- col(ncol)

  min.breaks <- min(breaks)
  max.breaks <- max(breaks)

  x[x<min.breaks] <- min.breaks
  x[x>max.breaks] <- max.breaks


  ## Calculate the plot layout
  if( missing(lhei) || is.null(lhei) )
    lhei <- c(keysize, 4)

  if( missing(lwid) || is.null(lwid) )
    lwid <- c(keysize, 4)

  if( missing(lmat) || is.null(lmat) )
    {
      lmat <- rbind(4:3, 2:1)

      if(!missing(ColSideColors)) { ## add middle row to layout
        if(!is.character(ColSideColors) || length(ColSideColors) != nc)
          stop("'ColSideColors' must be a character vector of length ncol(x)")
        lmat <- rbind(lmat[1,]+1, c(NA,1), lmat[2,]+1)
        lhei <- c(lhei[1], 0.2, lhei[2])
      }

      if(!missing(RowSideColors)) { ## add middle column to layout
        if(!is.character(RowSideColors) || length(RowSideColors) != nr)
          stop("'RowSideColors' must be a character vector of length nrow(x)")
        lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-1), 1), lmat[,2]+1)
        lwid <- c(lwid[1], 0.2, lwid[2])
      }

      lmat[is.na(lmat)] <- 0
    }

  if(length(lhei) != nrow(lmat))
    stop("lhei must have length = nrow(lmat) = ", nrow(lmat))

  if(length(lwid) != ncol(lmat))
    stop("lwid must have length = ncol(lmat) =", ncol(lmat))

  ## Graphics `output' -----------------------

  op <- par(no.readonly = TRUE)
  on.exit(par(op))
  layout(lmat, widths = lwid, heights = lhei, respect = FALSE)

  plot.index <- 1

  ## draw the side bars
  if(!missing(RowSideColors)) {
    par(mar = c(margins[1],0, 0,0.5))
    image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    plot.index <- plot.index + 1
  }
  if(!missing(ColSideColors)) {
    par(mar = c(0.5,0, 0,margins[2]))
    image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    plot.index <- plot.index + 1
  }
  ## draw the main carpet
  par(mar = c(margins[1], 0, 0, margins[2]))
  #if(scale != "none" || !symm)
  #  {
      x <- t(x)
      cellnote <- t(cellnote)
  #  }
  if(revC)
    { ## x columns reversed
      iy <- nr:1
      if(exists("ddr"))
        ddr <- rev(ddr)
      x <- x[,iy]
      cellnote <- cellnote[,iy]
    }
  else iy <- 1:nr

  ## display the main carpet
  image(1:nc, 1:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
        axes = FALSE, xlab = "", ylab = "", col=col, breaks=breaks,
        ...)
  retval$carpet <- x
  if(exists("ddr"))
    retval$rowDendrogram <- ddr
  if(exists("ddc"))
    retval$colDendrogram <- ddc
  retval$breaks <- breaks
  retval$col <- col

  ## fill 'na' positions with na.color
  if(!invalid(na.color) & any(is.na(x)))
    {
      mmat <- ifelse(is.na(x), 1, NA)
      image(1:nc, 1:nr, mmat, axes = FALSE, xlab = "", ylab = "",
            col=na.color, add=TRUE)
    }

  ## add column labels
  if(is.null(srtCol) && is.null(colCol))
    axis(1,
         1:nc,
         labels= labCol,
         las= 2,
         line= -0.5 + offsetCol,
         tick= 0,
         cex.axis= cexCol,
         hadj=adjCol[1],
         padj=adjCol[2]
         )
  else
    {
      if(is.null(srtCol) || is.numeric(srtCol))
        {
          if(missing(adjCol) || is.null(adjCol))
            adjCol=c(1,NA)

          if(is.null(srtCol))
            srtCol <- 90

          xpd.orig <- par("xpd")
          par(xpd=NA)
          xpos <- axis(1, 1:nc, labels=rep("", nc), las=2, tick=0)
          text(x=xpos,
               y=par("usr")[3] - (1.0 + offsetCol) * strheight("M"),
               labels=labCol,
               ##pos=1,
               adj=adjCol,
               cex=cexCol,
               srt=srtCol,
               col=colCol
               )
          par(xpd=xpd.orig)
        }
      else
        warning("Invalid value for srtCol ignored.")
    }

  ## add row labels
  if(is.null(srtRow) && is.null(colRow))
    {
      axis(4,
           iy,
           labels=labRow,
           las=2,
           line=-0.5+offsetRow,
           tick=0,
           cex.axis=cexRow,
           hadj=adjRow[1],
           padj=adjRow[2]
           )
    }
  else
    {
      if(is.null(srtRow) || is.numeric(srtRow))
        {
          xpd.orig <- par("xpd")
          par(xpd=NA)
          ypos <- axis(4, iy, labels=rep("", nr), las=2, line= -0.5, tick=0)
          text(x=par("usr")[2] + (1.0 + offsetRow) * strwidth("M"),
               y=ypos,
               labels=labRow,
               adj=adjRow,
               cex=cexRow,
               srt=srtRow,
               col=colRow
               )
          par(xpd=xpd.orig)
        }
      else
        warning("Invalid value for srtRow ignored.")
    }



  ## add row and column headings (xlab, ylab)
  if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
  if(!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25)

  ## perform user-specified function
  if (!missing(add.expr))
    eval(substitute(add.expr))

  ## add 'background' colored spaces to visually separate sections
  if(!missing(colsep))
      for(csep in colsep)
        rect(xleft =csep+0.5,               ybottom=0,
             xright=csep+0.5+sepwidth[1],   ytop=ncol(x)+1,
             lty=1, lwd=1, col=sepcolor, border=sepcolor)

  if(!missing(rowsep))
    for(rsep in rowsep)
      rect(xleft =0,          ybottom= (ncol(x)+1-rsep)-0.5,
           xright=nrow(x)+1,  ytop   = (ncol(x)+1-rsep)-0.5 - sepwidth[2],
           lty=1, lwd=1, col=sepcolor, border=sepcolor)


  ## show traces
  min.scale <- min(breaks)
  max.scale <- max(breaks)
  x.scaled  <- scale01(t(x), min.scale, max.scale)

  if(trace %in% c("both","column") )
    {
      retval$vline <- vline
      vline.vals <- scale01(vline, min.scale, max.scale)
      for( i in 1:length(colInd) )
        {
          if(!is.null(vline))
            {
              abline(v=i-0.5 + vline.vals, col=linecol, lty=2)
            }
          xv <- rep(i, nrow(x.scaled)) + x.scaled[,i] - 0.5
          xv <- c(xv[1], xv)
          yv <- 1:length(xv)-0.5
          lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
        }
    }


  if(trace %in% c("both","row") )
    {
      retval$hline <- hline
      hline.vals <- scale01(hline, min.scale, max.scale)
      for( i in 1:length(rowInd) )
        {
          if(!is.null(hline))
            {
              abline(h=i - 0.5 + hline.vals, col=linecol, lty=2)
            }
          yv <- rep(i, ncol(x.scaled)) + x.scaled[i,] - 0.5
          yv <- rev(c(yv[1], yv))
          xv <- length(yv):1-0.5
          lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
        }
    }



  if(!missing(cellnote))
    text(x=c(row(cellnote)),
         y=c(col(cellnote)),
         labels=c(cellnote),
         col=notecol,
         cex=notecex)

  plot.index <- plot.index + 1

  ## increment plot.index and then do
  ##   latout_set( lmat, plot.index )
  ## to set to the correct plot region, instead of
  ## relying on plot.new().

  ## the two dendrograms :
  par(mar = c(margins[1], 0, 0, 0))
  if( dendrogram %in% c("both","row") )
    {
        flag <- try(
            plot.dendrogram(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
            )
        if("try-error" %in% class(flag))
            {
                cond <- attr(flag, "condition")
                if(!is.null(cond) && conditionMessage(cond)=="evaluation nested too deeply: infinite recursion / options(expressions=)?")
                    stop('Row dendrogram too deeply nested, recursion limit exceeded.  Try increasing option("expressions"=...).')
            }
    }
  else
    plot.new()

  par(mar = c(0, 0, if(!is.null(main)) 5 else 0, margins[2]))

  if( dendrogram %in% c("both","column") )
    {
        flag <- try(
            plot.dendrogram(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
            )
        if("try-error" %in% class(flag))
            {
                cond <- attr(flag, "condition")
                if(!is.null(cond) && conditionMessage(cond)=="evaluation nested too deeply: infinite recursion / options(expressions=)?")
                    stop('Column dendrogram too deeply nested, recursion limit exceeded.  Try increasing option("expressions"=...).')
            }
    }
  else
    plot.new()

  ## title
  if(!is.null(main)) title(main, cex.main = 1.5*op[["cex.main"]])

  ## Add the color-key
  if(key)
    {
      mar <- c(5, 4, 2, 1)
      if (!is.null(key.xlab) && is.na(key.xlab))
          mar[1] <- 2
      if (!is.null(key.ylab) && is.na(key.ylab))
          mar[2] <- 2
      if (!is.null(key.title) && is.na(key.title))
          mar[3] <- 1
      par(mar = mar, cex=0.75, mgp=c(2, 1, 0))
      if (length(key.par) > 0)
          do.call(par, key.par)
      tmpbreaks <- breaks

      if(symkey)
        {
          max.raw <- max(abs(c(x,breaks)),na.rm=TRUE)
          min.raw <- -max.raw
          tmpbreaks[1] <- -max(abs(x), na.rm=TRUE)
          tmpbreaks[length(tmpbreaks)] <- max(abs(x), na.rm=TRUE)
        }
      else
        {
          min.raw <- min.breaks
          max.raw <- max.breaks
        }

      z <- seq(min.raw, max.raw, length.out=keygrid)
      image(z=matrix(z, ncol=1),
            #x=matrix(z),
            #y=0,
            col=col,
            breaks=tmpbreaks,
            xlab="",
            xaxt="n",
            yaxt="n")

      if (is.null(key.xtickfun)) {
        usr <- par("usr")
        par(usr=c(usr[1], usr[2], 0, 1))
        lv <- pretty(breaks)
        xv <- scale01(as.numeric(lv), min.raw, max.raw)
        xargs <- list(at=xv, labels=lv)
      } else {
        xargs <- key.xtickfun()
      }


      xargs$side <- 1
      do.call(axis, xargs)
      if (is.null(key.xlab)) {
        if(scale=="row")
            key.xlab <- "Row Z-Score"
        else if(scale=="column")
           key.xlab <- "Column Z-Score"
        else
            key.xlab <- "Value"
      }
      if (!is.na(key.xlab)) {
          mtext(side=1, key.xlab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
      }

      if(density.info=="density")
        {
          dens <- density(x, adjust=densadj, na.rm=TRUE,
                          from=min.scale, to=max.scale)
          omit <- dens$x < min(breaks) | dens$x > max(breaks)
          dens$x <- dens$x[!omit]
          dens$y <- dens$y[!omit]
          dens$x <- scale01(dens$x, min.raw, max.raw)
          lines(dens$x, dens$y / max(dens$y) * 0.95, col=denscol, lwd=1)
          if (is.null(key.ytickfun)) {
              yargs <- list(at=pretty(dens$y)/max(dens$y) * 0.95, labels=pretty(dens$y))
          } else {
              yargs <- key.ytickfun()
          }
          yargs$side <- 2
          do.call(axis, yargs)
          if (is.null(key.title))
              key.title <- "Color Key\nand Density Plot"
          if (!is.na(key.title))
            title(key.title)
          par(cex=0.5)
          if (is.null(key.ylab))
              key.ylab <- "Density"
          if (!is.na(key.ylab))
            mtext(side=2,key.ylab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
        }
      else if(density.info=="histogram")
        {
          h <- hist(x, plot=FALSE, breaks=breaks)
          hx <- scale01(breaks, min.raw, max.raw)
          hy <- c(h$counts, h$counts[length(h$counts)])
          lines(hx, hy/max(hy)*0.95, lwd=1, type="s", col=denscol)
          if (is.null(key.ytickfun)) {
              yargs <- list(at=pretty(hy)/max(hy) * 0.95, labels=pretty(hy))
          } else {
              yargs <- key.ytickfun()
          }
          yargs$side <- 2
          do.call(axis, yargs)
          if (is.null(key.title))
              key.title <- "Color Key\nand Histogram"
          if (!is.na(key.title))
            title(key.title)
          par(cex=0.5)
          if (is.null(key.ylab))
              key.ylab <- "Count"
          if (!is.na(key.ylab))
            mtext(side=2,key.ylab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
        }
      else
          if (is.null(key.title))
              title("Color Key")

      if(trace %in% c("both","column") )
          {
              vline.vals <- scale01(vline, min.raw, max.raw)
              if(!is.null(vline))
                  {
                      abline(v=vline.vals, col=linecol, lty=2)
                  }
          }


      if(trace %in% c("both","row") )
          {
              hline.vals <- scale01(hline, min.raw, max.raw)
              if(!is.null(hline))
                  {
                      abline(v=hline.vals, col=linecol, lty=2)

                  }
          }

    }
  else
    {
      par(mar=c(0, 0, 0, 0))
      plot.new()
    }
  ## Create a table showing how colors match to (transformed) data ranges
  retval$colorTable <- data.frame(
                             low=retval$breaks[-length(retval$breaks)],
                             high=retval$breaks[-1],
                             color=retval$col
                             )

  # Store layout information, suggested by Jenny Drnevich
  retval$layout <- list(lmat = lmat,
                        lhei = lhei,
                        lwid = lwid
                        )


  ## If user has provided an extra function, call it.
  if(!is.null(extrafun))
      extrafun()

  invisible( retval )
}