R/CorrelationHeatmap.R

In DanteR: Proteomics data analysis tool

# Written by Ashoka D. Polpitiya
# for the Department of Energy (PNNL, Richland, WA)
# Copyright 2007, Battelle Memorial Institute
# E-mail: ashoka.polpitiya@pnl.gov
# Website: http://omics.pnl.gov/software
# -------------------------------------------------------------------------
#
# Licensed under the Apache License, Version 2.0; you may not use this file except
# in compliance with the License.  You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
#
# R Plotting functions used in DAnTE
# -------------------------------------------------------------------------


CorrelationHeatmap.dialog <- list(
  label = "Plot Correlation Heatmap",
  x.dataframeItem="T_Data", label = "Data source", 
    #signal = c("default", "get.colnames", "Columns"),
    signal = c("default", "get.dataset.factors", "order.key"),
    signal = c("default", "get.dataset.factors", "color.key"),
  tooltip = "A numeric matrix or data frame",
  method.choiceItem = c(value="pearson", "kendall", "spearman"), label = "Correlation Method", tooltip = "If 'method' is 'kendall' or 'spearman', Kendall's\ntau or Spearman's rho statistic is used to estimate a rank-based\nmeasure of association.  These are more robust and have been\nrecommended if the data do not necessarily come from a bivariate\nnormal distribution.",
  use.choiceItem = c("all.obs", "complete.obs", "na.or.complete", value = "pairwise.complete.obs"), item.labels = c("Everything", "Complete Only", "NA Or Complete", "Pairwise Complete"), label = "Missing Value Behavior", tooltip = "Method for computing covariances in the presence of missing values - Pairwise Complete is usually used",
  BREAK=T,
  orderByFactor.trueFalseItem=F, label = "Order By Factor?", 
    tooltip = "Arrange columns of data by some factor in linked Column Metadata",
    signal = c("default", "toggle.sensitive", "order.key"),
    order.key.choiceItem=NULL, label="Choose Factor", indent=10,
  colorByFactor.trueFalseItem=F, label = "Color Sidebars?",
    tooltip = "Add a sidebar indicating factors in linked Column Metadata",
    signal = c("default", "toggle.sensitive", "color.key"),
    color.key.choiceItem=NULL, label="Choose Factor", indent=10,
  cMap.choiceItem = c("BlackBody", value="GreenRed" , "Heat", "BlueWhiteRed"), label="Color Map",
  tooltip = "Which Color Scheme to use",
  labelscale.rangeItem = c(value=1, from=0, to=1, by=0.1), label="Label Scale"
) 
  
CorrelationHeatmap <- function(x, 
                            file="deleteme.png",
                            cMap="BlackBody",
    method = "pearson",
    use = "pairwise.complete.obs",
                            corRangeMin=0, corRangeMax=1,
                            bkground="white",
                            customColors=c("#FF0000","#00FF00", "#0000FF"),
                            labelscale=1,
                            stamp=NULL,
                            orderByFactor = FALSE,
                            order.key = NULL,
                            colorByFactor = FALSE,
                            color.key = NULL,
                            ...)
{
  #png(filename=file,width=1152,height=864,pointsize=12,bg=bkground,
  #          res=600)
  ##require(Cairo)
  ##CairoPNG(filename=file,width=IMGwidth,height=IMGheight,pointsize=FNTsize,bg=bkground)
  #par(oma=c(3.4, 2, 2, 2), mar=c(25,25,4,1))
  if(orderByFactor && length(order.key)==1)
    x <- OrderByFactor(x, order.key)

  corRange <- c(corRangeMin, corRangeMax)
    # Factor colors
  cbf <- NULL
  if(colorByFactor && length(color.key)==1)
    cbf <- ColorByFactor(x, color.key)
  
  Xcor <- cor(x,method=method,use=use)
  Xcor <- abs(Xcor)
  corRange <- range(Xcor, na.rm=T)

  ##### Correction for correlation range selection #########
  Xcor[Xcor > corRange[2]] <- corRange[2]
  Xcor[Xcor < corRange[1]] <- corRange[1]
  Xcor[1,1] <- corRange[1] # to make sure the lower and upper limits are ..
  Xcor[2,2] <- corRange[2] # represented in the image. These will be
                           # replaced in heatmap.corr with 1
  ##########################################################

  cmap <- colorRampPalette(c("black", "red", "orange","yellow","lightyellow"),
                    space="rgb")(20)
  cmap <- switch (cMap,
        BlackBody = colorRampPalette(
                    c("black", "red", "orange","yellow","lightyellow"),
                    space="rgb")(20),
        GreenRed = colorRampPalette( c("green", "black", "red"),
                   space="rgb")(20),
        Heat = colorRampPalette(c("red", "orange","yellow","lightyellow"),
                    space="rgb")(20),
        BlueWhiteRed =  colorRampPalette( c("blue", "white", "red"),
                   space="rgb")(20),
        Custom = colorRampPalette(customColors, space="rgb")(20)
        )

    library(gplots)
    best.width <- function(v) 2 + max(0.3*nchar(v)) 
 
    margin=c(best.width(colnames(Xcor)),5)

    if(!is.null(cbf))
      heatmap.corr(Xcor, Rowv=FALSE,Colv=FALSE,trace="none", dendrogram="none",
           col=cmap,cexCol=labelscale,cexRow=labelscale,margin=margin,
           RowSideColors = cbf$color, ColSideColors = cbf$color, 
           showLegend = TRUE, Legend = cbf$legend, ...)
    else
      heatmap.corr(Xcor, Rowv=FALSE,Colv=FALSE,trace="none", dendrogram="none",
           col=cmap,cexCol=labelscale,cexRow=labelscale,margin=margin, ...)
            #cellnote=signif(Xcor,digits=2),notecex=.7,notecol=1,...)
    text(1,1,"abc",col="white") # ????
      if (!is.null(stamp))
            mtext(paste("DAnTE : ", format(Sys.time(), "%m-%d-%Y %I:%M%p"),
                  " (", stamp, ")", sep=""),col=1,cex=.6,line=3, side=1, adj=1)
  return(recordPlot())
}
                     
#--------------------------------------------------------
###############################################################
# This is the heatmap.2 command from gplots package
# modified to accomadate for the correlation range selection
###############################################################
heatmap.corr <- function (x, Rowv = TRUE, Colv = if (symm) "Rowv" else TRUE,
    distfun = dist, hclustfun = hclust, dendrogram = c("both",
        "row", "column", "none"), symm = FALSE, scale = c("none",
        "row", "column"), na.rm = TRUE, revC = identical(Colv,
        "Rowv"), add.expr, breaks, col = "heat.colors", colsep,
    rowsep, sepcolor = "white", sepwidth = c(0.05, 0.05), cellnote,
    notecex = 1, notecol = "cyan", na.color = par("bg"), trace = c("column",
        "row", "both", "none"), tracecol = "cyan", hline = median(breaks),
    vline = median(breaks), linecol = tracecol, margins = c(5,
        5), ColSideColors, RowSideColors, cexRow = 0.2 + 1/log10(nr),
    cexCol = 0.2 + 1/log10(nc), labRow = NULL, labCol = NULL,
    key = TRUE, keysize = 1.2, density.info = c("histogram",
        "density", "none"), denscol = tracecol, symkey = min(x <
        0, na.rm = TRUE), densadj = 0.25, main = NULL, xlab = NULL,
    ylab = NULL, showLegend = FALSE, Legend = NULL,
    ...)
{
    scale01 <- function(x, low = min(x), high = max(x)) {
        x <- (x - low)/(high - low)
        x
    }
    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 (!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 ((Colv == "Rowv") && (!isTRUE(Rowv) || is.null(Rowv)))
        Colv <- FALSE
    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")) {
        if (((!isTRUE(Rowv)) || (is.null(Rowv))) && (dendrogram %in%
            c("both", "row"))) {
            if (is.logical(Colv) && (Colv))
                dendrogram <- "column"
            else dedrogram <- "none"
            warning("Discrepancy: Rowv is FALSE, while dendrogram is `",
                dendrogram, "'. Omitting row dendogram.")
        }
    }
    if (!inherits(Colv, "dendrogram")) {
        if (((!isTRUE(Colv)) || (is.null(Colv))) && (dendrogram %in%
            c("both", "column"))) {
            if (is.logical(Rowv) && (Rowv))
                dendrogram <- "row"
            else dendrogram <- "none"
            warning("Discrepancy: Colv is FALSE, while dendrogram is `",
                dendrogram, "'. Omitting column dendogram.")
        }
    }
    if (inherits(Rowv, "dendrogram")) {
        ddr <- Rowv
        rowInd <- order.dendrogram(ddr)
    }
    else if (is.integer(Rowv)) {
        hcr <- hclustfun(distfun(x))
        ddr <- as.dendrogram(hcr)
        ddr <- reorder(ddr, Rowv)
        rowInd <- order.dendrogram(ddr)
        if (nr != length(rowInd))
            stop("row dendrogram ordering gave index of wrong length")
    }
    else if (isTRUE(Rowv)) {
        Rowv <- rowMeans(x, na.rm = na.rm)
        hcr <- hclustfun(distfun(x))
        ddr <- as.dendrogram(hcr)
        ddr <- reorder(ddr, Rowv)
        rowInd <- order.dendrogram(ddr)
        if (nr != length(rowInd))
            stop("row dendrogram ordering gave index of wrong length")
    }
    else {
        rowInd <- nr:1
    }
    if (inherits(Colv, "dendrogram")) {
        ddc <- Colv
        colInd <- order.dendrogram(ddc)
    }
    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)) {
        hcc <- hclustfun(distfun(if (symm)
            x
        else t(x)))
        ddc <- as.dendrogram(hcc)
        ddc <- reorder(ddc, Colv)
        colInd <- order.dendrogram(ddc)
        if (nc != length(colInd))
            stop("column dendrogram ordering gave index of wrong length")
    }
    else if (isTRUE(Colv)) {
        Colv <- colMeans(x, na.rm = na.rm)
        hcc <- hclustfun(distfun(if (symm)
            x
        else t(x)))
        ddc <- as.dendrogram(hcc)
        ddc <- reorder(ddc, Colv)
        colInd <- order.dendrogram(ddc)
        if (nc != length(colInd))
            stop("column dendrogram ordering gave index of wrong length")
    }
    else {
        colInd <- 1:nc
    }
    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 (scale == "row") {
        x <- sweep(x, 1, rowMeans(x, na.rm = na.rm))
        sx <- apply(x, 1, sd, na.rm = na.rm)
        x <- sweep(x, 1, sx, "/")
    }
    else if (scale == "column") {
        x <- sweep(x, 2, colMeans(x, na.rm = na.rm))
        sx <- apply(x, 2, sd, na.rm = na.rm)
        x <- sweep(x, 2, sx, "/")
    }
    if (missing(breaks) || is.null(breaks) || length(breaks) <
        1)
        if (missing(col))
            breaks <- 16
        else breaks <- length(col) + 1
    if (length(breaks) == 1) {
        breaks <- seq(min(x, na.rm = na.rm), max(x, na.rm = na.rm),
            length = breaks)
    }
    nbr <- length(breaks)
    ncol <- length(breaks) - 1
    if (class(col) == "function")
        col <- col(ncol)
    else if (is.character(col) && length(col) == 1)
        col <- do.call(col, list(ncol))
    min.breaks <- min(breaks)
    max.breaks <- max(breaks)
    x[] <- ifelse(x < min.breaks, min.breaks, x)
    x[] <- ifelse(x > max.breaks, max.breaks, x)
    lmat <- rbind(4:3, 2:1)
    lhei <- lwid <- c(keysize, 4)
    if (!missing(ColSideColors)) {
        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)) {
        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
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
    if (!missing(RowSideColors)) {
        par(mar = c(margins[1], 0, 0, 0.5))
        image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    }
    if (!missing(ColSideColors)) {
        par(mar = c(0.5, 0, 0, margins[2]))
        image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    }     
    par(mar = c(margins[1], 0, 0, margins[2]))
    if (!symm || scale != "none") {
        x <- t(x)
        cellnote <- t(cellnote)
    }
    if (revC) {
        iy <- nr:1
        ddr <- rev(ddr)
        x <- x[, iy]
        cellnote <- cellnote[, iy]
    }
    else iy <- 1:nr
    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, ...)
    ############Modification###################
    x.tmp <- matrix(NA, dim(x)[1], dim(x)[2])
    x.tmp[1,dim(x)[2]] <- 1 # this was replaced by min of corr range
    x.tmp[2,dim(x)[2]-1] <- 1 # this was replaced by max of corr range
    image(1:nc, 1:nr, x.tmp, axes = FALSE, xlab = "", ylab = "",
            col = col[length(col)], add = TRUE)
    ###########################################
    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)
    }
    axis(1, 1:nc, labels = labCol, las = 2, line = -0.5, tick = 0,
        cex.axis = cexCol)
    if (!is.null(xlab))
        mtext(xlab, side = 1, line = margins[1] - 1.25)
    axis(4, iy, labels = labRow, las = 2, line = -0.5, tick = 0,
        cex.axis = cexRow)
    if (!is.null(ylab))
        mtext(ylab, side = 4, line = margins[2] - 1.25)
    if (!missing(add.expr))
        eval(substitute(add.expr))
    if (!missing(colsep))
        for (csep in colsep) rect(xleft = csep + 0.5, ybottom = rep(0,
            length(csep)), xright = csep + 0.5 + sepwidth[1],
            ytop = rep(ncol(x) + 1, csep), 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 = ncol(x) + 1, ytop = (ncol(x) +
            1 - rsep) - 0.5 - sepwidth[2], lty = 1, lwd = 1,
            col = sepcolor, border = sepcolor)
    min.scale <- min(breaks)
    max.scale <- max(breaks)
    x.scaled <- scale01(t(x), min.scale, max.scale)
    if (trace %in% c("both", "column")) {
        for (i in colInd) {
            if (!is.null(vline)) {
                vline.vals <- scale01(vline, min.scale, max.scale)
                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")) {
        for (i in rowInd) {
            if (!is.null(hline)) {
                hline.vals <- scale01(hline, min.scale, max.scale)
                abline(h = i + hline, 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)
    par(mar = c(margins[1], 0, 0, 0))
    if (dendrogram %in% c("both", "row")) {
        plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
    }
    else plot.new()
    par(mar = c(0, 0, if (!is.null(main)) 5 else 0, margins[2]))
    if (dendrogram %in% c("both", "column")) {
        plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
    } else if (showLegend && !is.null(Legend)) {

      legend("topleft",names(Legend),col=Legend,pch=19,bg='transparent')
      plot.new()
    }  else {
      plot.new()
    }
    
    if (!is.null(main))
        title(main, cex.main = 1.5 * op[["cex.main"]])
    if(key){
        par(mar = c(5, 4, 2, 1), cex = 0.75)
        if (symkey) {
            max.raw <- max(abs(x), na.rm = TRUE)
            min.raw <- -max.raw
        }
        else {
            min.raw <- min(x, na.rm = TRUE)
            max.raw <- max(x, na.rm = TRUE)
        }
        z <- seq(min.raw, max.raw, length = length(col))
        image(z = matrix(z, ncol = 1), col = col, breaks = breaks,
            xaxt = "n", yaxt = "n")
        par(usr = c(0, 1, 0, 1))
        lv <- pretty(breaks)
        xv <- scale01(as.numeric(lv), min.raw, max.raw)
        axis(1, at = xv, labels = lv)
        if (scale == "row")
            mtext(side = 1, "Row Z-Score", line = 2)
        else if (scale == "column")
            mtext(side = 1, "Column Z-Score", line = 2)
        else mtext(side = 1, "Value", line = 2, cex=0.7)
        if (density.info == "density") {
            dens <- density(x, adjust = densadj, na.rm = TRUE)
            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)
            axis(2, at = pretty(dens$y)/max(dens$y) * 0.95, pretty(dens$y))
            title("Color Key\nand Density Plot", cex.main=.7)
            par(cex = 0.4)
            mtext(side = 2, "Density", line = 2, cex=0.7)
        }
        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)
            axis(2, at = pretty(hy)/max(hy) * 0.95, pretty(hy))
            title("Color Key\nand Histogram", cex.main=.9)
            par(cex = 0.5)
            mtext(side = 2, "Count", line = 2, cex=0.7)
        }
        else title("Color Key", cex.main=.9)
    }
    else plot.new()
    invisible(list(rowInd = rowInd, colInd = colInd))
}