help with a heatmap with deseq - legend?

[vebaev]

Hi,
I have made with deseq a heatmap, but it does not show a legend/colorkey and I need to have it?

I did not see elsewhere how colorkey/legends is made when heatmap is generated from sequencing copy number values, it should not de in log2 values or?

[Wolfgang Huber]

Dear Vebaev

DESeq does not do heatmaps. The programming language R, in which DESeq is implemented, has functions to make heatmaps, and I recommend that you have a look at the manual pages of some of these functions, for instance

## install (do only once)
source('http://www.bioconductor.org/biocLite.R')
biocLite("gplots")

## load (do each time you need it)
library("gplots")
help("heatmap.2")

As to transformation - the log2-transformation is problematic due to its singularity at 0 and its large slope for values close to 0. Doing no transformation is also problematic since then the dynamic range of the colour scale tends to be "wasted" on a sparse set of large values. My recommdation is to have a look at Section 6 of the DESeq vignette, Variance stabilizing transformation, where in Fig.13 you also see an example heatmap (but using the more primitive function heatmap that does not offer a color bar).

Hope this helps,

Wolfgang.

[gringer]

Here's an R function for producing a heatmap based on the vsd data (together with colour legend and dendrogram):

make.heatmap <- function(data.vsd, columns = dim(data.vsd)[2], col = colorRampPalette(blues9)(100), log = FALSE, top = dim(data.vsd)[1], method = "canberra"){
  hm.dist <- dist(t(head(data.vsd[,columns],top)), method = method);
  hm.dend <- as.dendrogram(hclust(hm.dist));
  scaled.mat <- as.matrix(hm.dist);
  if(log){
    scaled.mat <- log(scaled.mat);
  }
  diag(scaled.mat) <- NA;
  scaled.mat <- scaled.mat[order.dendrogram(hm.dend),rev(order.dendrogram(hm.dend))];
  scaled.range <- t(as.matrix(seq(range(scaled.mat, na.rm = TRUE)[1],
                                  range(scaled.mat, na.rm = TRUE)[2], length.out = 100)));

  layout(matrix(1:3,1,3),widths = c(1,2,8));
  par(mar = c(10,3,0.5,0.5));
  image(1,scaled.range[1,],scaled.range, col = col,
        xlab = "", ylab = "", xaxt = "n");
  plot(hm.dend, axes = FALSE, horiz = TRUE, yaxs = "i", leaflab = "none");
  par(mar = c(10,0.5,0.5,10));
  image(1:length(rownames(scaled.mat)),1:length(colnames(scaled.mat)),scaled.mat,
        col = col,
        xlab = "", ylab = "", xaxt = "n", yaxt = "n");
  axis(4, 1:length(colnames(scaled.mat)), labels = sub("_seq1","",colnames(scaled.mat)), las = 2,
       line = -0.5, tick = 0, cex.axis = 0.2 + 1/log10(length(colnames(scaled.mat))));
  axis(1, 1:length(rownames(scaled.mat)), labels = rownames(scaled.mat), las = 2,
       line = -0.5, tick = 0, cex.axis = 0.2 + 1/log10(6));
  invisible(scaled.mat);
}

I recommend top <= 1000 if you're looking at all transcripts.

[vebaev]

Thanks both of you, I will try your suggestions!