R/cmapHeatmap.R

In gCMAPWeb: A web interface for gene-set enrichment analyses

##' Function to create an annotated heatmap of gene scores
##' 
##' This function takes a numerical matrix (with samples in columns and genes in rows),
##' preprocesses the data (if desired), determines the optimal height for the
##' heatmap and calls the .ImagePlot function to create the final heatmap figure.
##'
##' @param x numerical matrix with samples in rows and genes in columns.
##' @param col.anno character vector with column annotations to be displayed as (horizontal) annotation bar above the heatmap. If not NULL, must contain one elment for each column of 'x'. 
##' @param row.anno character vector with row annotations to be displayed as (vertical) annotation to the right of the heatmap. If not NULL, must contain one element for each row of 'x'.
##' @param file.name character, the path and filename (without suffix) to save the png file to
##' @param reference.name character, names of the reference cmap, used to construct the html image reference
##' @param url.base character, prefix for the html image reference
##' @param main Character, main title of the plot
##' @param ColorRamp vector of colors used for the heatmap, e.g. generated by a call to colorRampPalette
##' @param col.col named vector with a color for each level of col.anno (e.g. c(up="firebrick", down="blue"))
##' @param row.col named vector with a color for each level of row.anno (e.g. c(correlated="firebrick",anticorrelated="#044381FF"))
##' @param order.by.score logical, should gene scores be reordered independently for each sample ?
##' @param score.cap numerical vector of length two, specifying the limits of the color scale. Scores > max( score.cap) or < min(score.cap) will be set to score.cap. Default: c(-5,5)
##' @param ylab character, y-axis label
##' @param cluster.rows logical, perform hierarchical clustering on significant gene sets ?
##' @return list with two elements, 1. image.html: a character string with html code with the image tag refering to the output png 2. row.order: integer vector with thw row.order obtained after hierarchical clustering or NULL if no clustering was performed
##' @author Thomas Sandmann
cmapHeatmap <- function( x,
                         reference.name,
                         col.anno=NULL,
                         row.anno=NULL, 
                         file.name="heatmap",
                         url.base=NULL,
                         main="Query gene scores",
                         ColorRamp=colorRampPalette(c("#044381FF", "grey95",
                                                      "grey95", "firebrick"))(100),
                         col.col=c(down="black", 
                                   up="grey"), 
                         row.col=c(correlated="#1B9E77", 
                                   anticorrelated="#044381FF",
                                   over="#1B9E77",
                                   under="#044381FF"),
                         order.by.score=TRUE,
                         cluster.rows=TRUE,
                         score.cap=c(-5,5),
                         ylab="Significant datasets"){    
  
  
  ## ensure that the annotation vectors match the dimensions of x
  if( !is.null( row.anno )){
    if(length( row.anno ) != nrow( x) ) {
      stop("Length( row.anno ) != nrow( x )")
    }
  }
  
  if( !is.null( col.anno )){
    if(length( col.anno ) != ncol( x) ) {
      stop("cmapHeatmap: length( col.anno ) != ncol( x )")
    }
  }
  ## remove columns without any valid values
  missing.columns <- apply(x, 2, function(x) all(is.na(x)))
  x<- x[,!missing.columns, drop=FALSE]
  
  if( !is.null( col.anno )){
    col.anno <- col.anno[! missing.columns, drop=FALSE]
  }    
  
  ## cluster & reorder rows
  if( cluster.rows == TRUE && nrow( x ) > 2){
    row.dendrogram <- hclust( dist( x ))
    row.order <- row.dendrogram$order
    branches <- nrow(x) - match( row.names( x )[row.order], 
                                 row.names( x )) + 1
    x <- x[row.order, ]
    if( !is.null( row.anno )){
      row.anno <- row.anno[row.order]
    }
  } else {   
    row.dendrogram <- NULL
    branches <- NULL 
    ## reverse row order
    if( nrow( x ) > 1){
      x <- x[rev( 1:nrow( x)),,drop=FALSE]
      row.anno <- row.anno[rev( 1:nrow( x))]
    }
  }
  
  ## reorder gene scores every row separately according to column-median
  if(order.by.score==TRUE){
    if( is.null( col.anno )){
      col.med <- apply( x, 2, median, na.rm=TRUE)
      x <-x[,order(col.med), drop=FALSE]
    } else {
      stopifnot(length( col.anno ) == ncol( x) )
      if( "down" %in% col.anno){
        x.down <- x[, col.anno =="down", drop=FALSE]
        col.med <- apply( x.down, 2, median, na.rm=TRUE)
        x.down <-x.down[,order(col.med), drop=FALSE]
      } else {
        x.down <- NULL
      } 
      if( "up" %in% col.anno){
        x.up <- x[, col.anno =="up", drop=FALSE]
        col.med <- apply( x.up, 2, median, na.rm=TRUE)
        x.up <-x.up[,order(col.med), drop=FALSE]
      } else {
        x.up <- NULL
      }
      x <- cbind( x.down, x.up)
      col.anno <- c( rep("down", sum( col.anno =="down")), 
                     rep("up", sum( col.anno =="up")))
    }
  }
  
  ## cap maximum and minimum scores
  x[ x < min(score.cap, na.rm=TRUE) ] <- min( score.cap, na.rm=TRUE )
  x[ x > max(score.cap, na.rm=TRUE) ] <- max( score.cap, na.rm=TRUE)
  
  ## dynamically determine plot height and the fraction occupied by the heatmap
  if( nrow( x ) >= 30){
    plot.height=800
    image.fraction=0.8
  } else {
    plot.height=356 + max( 20*nrow( x ), 60)
    image.fraction=max( 20*nrow( x ), 100) / 640
  }
  
  ##----- plot  
  png.file <- paste(file.name, "png", sep=".")
  png(
    png.file,
    height=plot.height,
    width=1000,
    res=200
    )
  .ImagePlot(x,
             ColorRamp=ColorRamp,
             col.col=col.col, 
             row.col=row.col,
             ylab=ylab,
             main=main,
             col.anno=col.anno,
             row.anno=row.anno,
             row.dendrogram=row.dendrogram,
             image.fraction=image.fraction,
             leaflab="none"
             )
  dev.off()
  
  png.file <- paste(
    file.name,
    "large.png",
    sep="."
    )
  png(
    png.file,
    height=plot.height,
    width=1000,
    res=200,
    pointsize=7)
  .ImagePlot(
    x,
    ColorRamp=ColorRamp,
    col.col=col.col, 
    row.col=row.col,
    ylab=ylab,
    main=main,
    col.anno=col.anno,
    row.anno=row.anno,
    row.dendrogram=row.dendrogram,
    image.fraction=image.fraction,
    leaflab="perpendicular")
  dev.off()
  
  image.html <- ifelse(
    is.null( url.base), 
    hwriteImage( file.path(reference.name, "figures", "heatmap.png"),
                link=file.path(reference.name, "figures","heatmap.large.png"),
                table=FALSE, br=TRUE, width=500),
    hwriteImage( file.path(url.base, reference.name, "figures","heatmap.png"), 
                link=file.path(reference.name, "figures","heatmap.large.png"),
                table=FALSE, br=TRUE, width=500) 
    )
  return(
    list(
      image.html = image.html,
      branches = branches
      )
    )
}

##############
##' Function to create an annotated heatmap
##'
##' This function takes a numerical matrix and optional row- and column-annotations
##' and plots an annotated heatmap using the image function.
##' 
##' @param x numerical matrix with samples in rows and genes in columns. 
##' @param ColorRamp vector of colors used for the heatmap, e.g. generated by a call to colorRampPalette 
##' @param col.col named vector with a color for each level of col.anno (e.g. c(up="firebrick", down="blue"))  
##' @param row.col named vector with a color for each level of row.anno (e.g. c(correlated="firebrick",anticorrelated="#044381FF")) 
##' @param main Character, main title of the plot
##' @param col.anno Character vector with column annotations to be displayed as (horizontal) annotation bar above the heatmap. If not NULL, must contain one elment for each column of x
##' @param row.anno character vector with row annotations to be displayed as (vertical) annotation to the right of the heatmap. If not NULL, must contain one element for each row of x
##' @param Either an bbject of calls 'hclust' to add as a row dendrogram or NULL.
##' @param ylab character, y-axis label  
##' @param image.fraction numeric, number between 0 and 1 determining the default fraction of the plot occupied by the heatmap.
##' @return A plot output to the default graphics device
##' @author Thomas Sandmann
.ImagePlot <- function(
  x, 
  ColorRamp=colorRampPalette(c("#044381FF", "grey95",
    "grey95", "firebrick"))(100),
  col.col=c(down="black", 
    up="grey"), 
  row.col=c(correlated="#1B9E77", 
    anticorrelated="#044381FF",
    over="#1B9E77",
    under="#044381FF"),
  main="Query gene scores",
  col.anno=NULL,
  row.anno=NULL,
  row.dendrogram=NULL,
  ylab="Significant datasets",
  image.fraction=0.8,
  leaflab="perpendicular" # either none or perpendicular
  ){
  
  ## suppress y-axis label for small numbers of rows
  if( nrow( x < 10)){
    ylab=""
  }
  
  ColorLevels <- seq(min(x, na.rm=TRUE),
                     max(x, na.rm=TRUE),
                     length=length(ColorRamp)+1)
  
  ## layout the plot
  par(oma=c(1.5,3,5,3))  
  layout(matrix(data=c(1,2,2,2,3,
                       4,5,5,5,6,
                       7,8,9,10,11), ncol=5, nrow=3, 
                byrow=TRUE), 
         heights=c(0.04, image.fraction, 0.2), widths=c(0.1,rep(0.15, 3),0.02)
  ) 
  
  ##--------- Frame 1: empty
  par(mar = c(0.1,0,0,0.05))
  frame()
  ##--------- Frame 2: column annotation bar
  if( !is.null( col.anno )){
    par(mar = c(0.1,0.5,0.5,0.05))
    if( leaflab != "none"){
      par(mar = c(0.1,0.5,1,0.05))
    }
    col.anno.matrix <- matrix( as.integer( factor( col.anno ) ), ncol=1)
    image( col.anno.matrix,
           col=col.col[ as.character( levels( factor( col.anno ))) ],
           axes=FALSE
    )
  } else {
    par(mar = c(0.1,0.5,0.5,0.05))
    frame()
  }
  
  ##--------- Frame 3: empty
  par(mar = c(0,0,0,0))  
  frame()
  ##--------- Frame 4: dendrogram
  par(mar = c(1.4,0,0,0), yaxs="i")  
  if( !is.null( row.dendrogram )){
    if( leaflab != "none"){
      par(mar = c(1.4,0,0,3))
    }
    plot(
      as.dendrogram(
        row.dendrogram
        ),
      horiz=TRUE, 
      axes=FALSE,
      leaflab=leaflab)
  } else {
    frame()
  }
  
  ##--------- Frame 5: heatmap
  par(mar = c(1.4,0.5,0.1,0.05))
  na.grey <- matrix(NA,
                    ncol=ncol(x),
                    nrow=nrow(x)
                    )
  na.grey[ is.na(x)] <- 1
  
  ## first plot NA positions in grey
  image(
    1:ncol(x),
    1:nrow(x),
    t(na.grey),
    col="grey",
    zlim=c(1,1),
    axes=FALSE,
    xlab="",
    ylab=""
    )
  
  ## next overplot real values
  image(
    1:ncol(x),
    1:nrow(x),
    t(x),
    axes=FALSE,
    col=ColorRamp, 
    zlim=c(min(x),max(x)),
    breaks=ColorLevels,
    ylab=ylab,
    add=TRUE)
  box(lwd=1)
  mtext(main,
        side = 3,
        line = 1.2,
        adj=0.6,
        outer = TRUE, 
        font=par()$font.main,
        cex=par()$cex.main)
  mtext(
    ylab,
    side = 4,
    line = 0.1,
    outer = TRUE,
    las=3,
    adj=0.75,
    cex=0.9
    )
  
  ##--------- Frame 6: row annotation bar
  if( !is.null( row.anno )){
    par(mar = c(1.4,0.05,0.1,0.5))
    image(
      matrix( as.integer( factor( row.anno ) ), nrow=1),
      col=row.col[ levels(factor( row.anno ))],
      axes=FALSE
      )
  } else {
    par(
      mar = c(0.5,0.05,0.05,0.5)
      )
    frame()
  }
  
  ##--------- Frame 7: empty
  par(mar = c(1.5,0,0,0))
  frame()
  
  ##--------- Frame 8: color key
  par(mar = c(1.5,1,1.9,1), las=1, cex=0.6)
  if( leaflab != "none"){
    par(mar = c(3.5,1,4,1))
  }
  image(
    ColorLevels,
    1,
    matrix(
      data=ColorLevels,
      nrow=length(ColorLevels),
      ncol=1
      ),
    col=ColorRamp,
    xlab="",ylab="",
    yaxt="n",
    cex=0.5,
    mgp = c(3, 0.5, 0))
  mtext("Gene-level z-score",side = 3,las=1,cex=0.6,line=0.35)  
  
  ##--------- Frame 9: column legend
  if( !is.null( col.anno )){
    par(
      mar = c(0,1,0.5,0),
      las=1,
      cex=1,
      xpd=NA)
    frame()
    legend(
      "right",
      cex=0.6,
      title="Query gene direction",
      lty=c(1,1), lwd=3,
      legend=unique(
        levels(
          factor( col.anno )
          )
        ), 
      col=c(col.col)[as.character(
        unique(levels(
          factor( col.anno )
          )
               )
        )], 
      horiz=FALSE,
      bty="n")
  } else {
    par(
      mar = c(0,1,0.5,0),
      las=1,
      cex=1,
      xpd=NA)
    frame()
  }
  ##--------- Frame 10: row legend
  if( !is.null( row.anno )){
    par(mar = c(0,1,0.5,0), las=1, cex=1, xpd=NA)
    frame()
    legend(
      "right",
      cex=0.6, 
      title="Set direction",
      legend=levels(factor( row.anno )),
      fill=row.col[ levels(factor( row.anno ))], 
      horiz=FALSE,
      bty="n", border=FALSE)    
  } else { 
    par(mar = c(0,1,0.5,1), las=1, cex=1, xpd=NA)
    frame()
  }
}