R/heatmap.R

In whtns/varbintools: Sequence Based Copy Number Detection Using Variable Bin Lengths

#!/usr/bin/env Rscript

# Heatmap for all CNV Profiles
# Author: Min Hu
# Last Updated: Nov 3th, 2018



#' Heatmap in Varbin
#'
#' @return
#' @export
#'
#' @examples
varbin_heatmap <- function() {

  set.seed(100)

  for (myStr in commandArgs() )
  {
  	message(" processing command arg:",myStr)
  	if(length(grep ("^-sample=", myStr)) >0)
  	{
  	dataName <- substring(myStr, nchar("-sample=")+1)
  	}
  	if(length(grep("^-outdir=",myStr))>0)
  	{
  	outdir <- substring(myStr, nchar("-outdir=")+1)
  	}
  }
  #data
  final_dir <- paste(outdir,"final_result",sep="/")
  heatmap_dir<-paste(final_dir,"heatmap",sep="/")
  if(!file.exists(heatmap_dir)){
   dir.create(heatmap_dir)
  }

  segfile=list.files(final_dir,pattern="seg.txt")

  segmentInfo  <- read.delim(paste(final_dir,segfile,sep="/"), header= TRUE)
  ##uber input uber.seg.txt, which contains segmented ratio values in columns of each cell, first three columns contain chr, start, absolute position

  ##reads files
  stat_dir<-paste(outdir,"metrics",sep="/")
  stat_file<-list.files(stat_dir,pattern="all_stat_metrics.txt")
  raw_reads<-read.delim(paste(stat_dir,stat_file,sep="/"),header=TRUE,sep="\t")
  rownames(raw_reads)<-raw_reads$Sample.Name
  raw_reads$Sample.Name<-NULL

  ##CNV ratio files

  ratio_file=list.files(final_dir,pattern="ratio.txt")
  cnv_ratio<-read.table(paste(final_dir,ratio_file,sep="/"),header=TRUE)


  ## load library
  #library("devtools")
  n<-ncol(segmentInfo)

  ##get the segment ratio data
  sam <- t(as.matrix(segmentInfo[,4:n]))

  ##assign column names
  colnames(sam) <- as.vector(segmentInfo[,1])

  ##asign sample names as rownames
  rownames(sam) <- colnames(segmentInfo)[4:n]
  rownames(sam) <-gsub("\\.","-",rownames(sam))

  ## filter out the cells with Reads counts less than 1M
  selected_cell<-rownames(raw_reads[raw_reads$TotalReads>1000000,])
  sam<-sam[match(selected_cell,rownames(sam)),]
  dim(sam)

  tryCatch(expr = { library("flowViz")},
           error = function(e) {
             source("https://bioconductor.org/biocLite.R")
             biocLite("flowViz")},
           finally = library("flowViz"))

  #require(IDPmisc)
  #sam<-NaRV.omit(sam)

  # transform seg ratio data for analysis
  mat <-as.matrix(sam) +0.1
  mat <- log2(mat)


  # hierarchycal clustering ############################################


  #library(gplots)

  ###set colobreaks for heatmap
  palettes <- colorRampPalette(c("blue", "white", "red"))(n = 999)
  breaks = c(seq(-2,-0.5,length=50),seq(-0.5,-0.3,length=200),seq(-0.3,0.3,length=300),seq(0.3,0.5,length=400),seq(0.5,3,length=50))

  ##color bars
  chr<- as.numeric(colnames(mat)) %% 2+1
  rbpal <- colorRampPalette(c("grey","black"))
  CHR <- cbind(rbpal(2)[as.numeric(chr)], rbpal(2)[as.numeric(chr)])



  #####heatmap
  #jpeg("/volumes/lab/users/minhu/Projects/2018-09-26-SCNA-kaile/output_1M/final_result/heatmap/heatmap.jpg", height=600, width=800)
  #heatmap.3(as.matrix(mat), dendrogram="none", distfun = dist, hclustfun = function(x) hclust(x, method='ward.D2'),Colv=NA,
  #          ColSideColors=CHR, Rowv=TRUE,notecol="black",col=palettes,breaks=breaks, symm=F,symkey=F,symbreaks=T,trace="none",
  #          cexRow=0.5, plot.row.partition=TRUE)
  #dev.off()

  tryCatch(expr = { library("pheatmap")},
           error = function(e) {
             source("https://bioconductor.org/biocLite.R")
             biocLite("pheatmap")},
           finally = library("pheatmap"))


  tryCatch(expr = { library("RColorBrewer")},
           error = function(e) {
            install.packages("RColorBrewer")},
           finally = library("RColorBrewer"))

  data_col_anno<-as.data.frame(cbind(chr))
  colnames(mat)<-rownames(data_col_anno)

  data_col_anno$chr<-gsub("2","n",data_col_anno$chr)
  data_col_anno$chr<-gsub("1","2n",data_col_anno$chr)

  ann_colors = list(chr = c("n"="#BEBEBE", "2n"= "#000000"))


  my_hclust_cell <- hclust(dist(mat), method = "ward.D2")

  # load package

  tryCatch(expr = { library("dendextend")},
           error = function(e) {
             source("https://bioconductor.org/biocLite.R")
             biocLite("dendextend")},
           finally = library("dendextend"))



  #library(dendextend)
  my_cell_col <- cutree(tree = as.dendrogram(my_hclust_cell), k = 3)
  my_cell_col<-as.data.frame(my_cell_col)
  my_cell_col$my_cell_col<-gsub("1","cluster 1",my_cell_col$my_cell_col)
  my_cell_col$my_cell_col<-gsub("2","cluster 2",my_cell_col$my_cell_col)
  my_cell_col$my_cell_col<-gsub("3","cluster 3",my_cell_col$my_cell_col)
  colnames(my_cell_col) <- "clusters"

  cv <- function(prof) sd(prof) / mean(prof)

  # Input normalized ratio values
  five.step.autocorrelation <- function(prof)
  {
    right <- prof[-(1:5)]
    left <- prof[-((length(prof)-4):length(prof))]
    cor(left, right)
  }


  all<- data.frame()
  for(x in 4:ncol(cnv_ratio)){

    temp<-cnv_ratio[,x]
    temp_cv<-cv(temp)
    all[(x-3),1]<-colnames(cnv_ratio)[x]
    all[(x-3),2]<-temp_cv
    temp_autoco<-five.step.autocorrelation(temp)
    all[(x-3),3]<-temp_autoco
    #all[(x-3),4]<-total_reads_per_cell2[all[(x-3),1],2]

  }
  rownames(all)<-all[,1]
  rownames(all) <-gsub("\\.","-",rownames(all))
  colnames(all)<-c("cell_id","cv","autoco")
  all$cell_id<-NULL

  all<-all[selected_cell,]
  #all<-NaRV.omit(all)


  my_cell_col$Total_Raw_Reads<-raw_reads[match(rownames(my_cell_col),rownames(raw_reads)),]$TotalReads
  my_cell_col$Correlation_of_variation<-all[match(rownames(my_cell_col),rownames(all)),]$cv
  my_cell_col$Autocorrelation_of_variation<-all[match(rownames(my_cell_col),rownames(all)),]$autoco

  my_cell_col$Correlation_of_variation[which(my_cell_col$Correlation_of_variation>1)] =1


  palettes <- colorRampPalette(c("blue", "white", "red"))(n = 999)
  breaks = c(seq(-3,-0.5,length=50),seq(-0.49,-0.3,length=200),seq(-0.29,0.3,length=300),seq(0.31,0.5,length=400),seq(0.51,3,length=50))

  heatmap_file<- paste(heatmap_dir,"heatmap_cluster.jpeg",sep="/")

  jpeg(heatmap_file, height=1000, width=1200)
  pheatmap(mat, annotation_row = my_cell_col, annotation_col = data_col_anno,annotation_colors=ann_colors,annotation_names_col = FALSE, color =palettes,breaks=breaks, clustering_method ="ward.D2",cluster_cols =FALSE,cutree_rows =3,show_colnames=FALSE)

  dev.off()

  heatmap_clusters <-paste(heatmap_dir,"heatmap.cluster.csv",sep="/")
  write.csv(my_cell_col,file=heatmap_clusters)
}

#' Varbin Consensus Complex Heatmap
#'
#' @param seg_df
#' @param seg_df_cells_t
#'
#' @return
#' @export
#'
#' @examples
varbin_consensus_complexheatmap <- function(seg_df, seg_df_cells_t){
  clusters_list <- split(seg_df_cells_t, cl$cluster)

  # we can use apply function to calculate the median
  clusters_consensus <- purrr::map_dfr(clusters_list, function(x) apply(x, 2, median))
  clusters_consensus <- as.data.frame(t(clusters_consensus))

  consensus_row_anno <- ComplexHeatmap::rowAnnotation(clusters = unique(sort(cl$cluster)),
                                      col = list(clusters = row_col))

  chr_lengths <- seg_df %>%
    dplyr::select(abspos, chrom, chrompos) %>%
    dplyr::group_by(chrom) %>%
    dplyr::summarize(n = dplyr::n()) %>%
    pull(n)

  chr_bar <- create_chr_bar(seg_df)

  # heatmap consensus
  ht_consensus <- ComplexHeatmap::Heatmap(as.matrix(log2(clusters_consensus+1e-3)),
                          cluster_columns = FALSE,
                          cluster_rows = FALSE,
                          show_row_names = FALSE,
                          show_column_names = FALSE,
                          use_raster = TRUE,
                          top_annotation = chr_bar,
                          border = TRUE,
                          col = circlize::colorRamp2(breaks = c(-1.5,-0.25,1.5),
                                                     c("dodgerblue3", "white", "firebrick3")),
                          heatmap_legend_param = list(title = "Log2 (Ratio)"))

  ComplexHeatmap::draw(consensus_row_anno + ht_consensus)
}

#' Create Chromosome Bar
#'
#' @param seg_df
#'
#' @return
#' @export
#'
#' @examples
create_chr_bar <- function(seg_df) {
  chr_lengths <- seg_df %>%
    dplyr::select(abspos, chrom, chrompos) %>%
    dplyr::group_by(chrom) %>%
    dplyr::summarize(n = dplyr::n()) %>%
    pull(n)

  chr_binary <- rep(c(2,1), 12)
  chr <- data.frame(chr = rep.int(x = chr_binary, times = chr_lengths))
  # getting lengths for chr numbers annotation
  chr_rl_c <- c(1, cumsum(chr_lengths))
  # creating a data frame to calculate rowMeans
  chr_df <-  data.frame(a = chr_rl_c[1:length(chr_rl_c)-1],b= chr_rl_c[2:length(chr_rl_c)])
  chr_l_means <- round(rowMeans(chr_df))
  chrom.names <- c(1:22,"X", "Y")
  # creating the vector for chr number annotations
  v <- vector(mode = "character",length = cumsum(chr_lengths)[length(chr_lengths)])
  v[chr_l_means] <- chrom.names
  v[is.na(v)] <- ""
  chr_bar <- ComplexHeatmap::HeatmapAnnotation(chr_text = ComplexHeatmap::anno_text(v,
                                                                                    gp = grid::gpar(fontsize = 8)),
                                               df = chr,
                                               show_legend = FALSE,
                                               which = "column",
                                               col = list(chr = c("1" = "grey88", "2" = "black"))
  )
}


#' Varbin Complex Heatmap
#'
#' @param seg_df
#' @param seg_df_cells_t
#'
#' @return
#' @export
#'
#' @examples
varbin_complexheatmap <- function(seg_df, seg_df_cells_t){
  # chromosome annotation top bar
  # getting the vector of chrom lengths

  chr_bar <- create_chr_bar(seg_df)

  heat_row_col <- ComplexHeatmap::rowAnnotation(clusters = sort(cl$cluster),
                                col = list(clusters = row_col))

  # heatmap
  ht <- ComplexHeatmap::Heatmap(as.matrix(log2(seg_df_cells_t+1e-3)),
                cluster_columns = FALSE,
                border = TRUE,
                col = circlize::colorRamp2(breaks = c(-1.5,-0.25,1.5),
                                           c("dodgerblue3", "white", "firebrick3")),
                cluster_rows = FALSE,
                show_row_names = FALSE,
                show_column_names = FALSE,
                use_raster = TRUE,
                top_annotation = chr_bar,
                heatmap_legend_param = list(title = "Log2 (Ratio)"))

  ComplexHeatmap::draw(heat_row_col + ht)
}

#' Heatmap 3
#'
#' @param x
#' @param Rowv
#' @param Colv
#' @param distfun
#' @param hclustfun
#' @param dendrogram
#' @param symm
#' @param scale
#' @param na.rm
#' @param revC
#' @param add.expr
#' @param breaks
#' @param symbreaks
#' @param col
#' @param colsep
#' @param rowsep
#' @param sepcolor
#' @param sepwidth
#' @param cellnote
#' @param notecex
#' @param notecol
#' @param na.color
#' @param trace
#' @param tracecol
#' @param hline
#' @param vline
#' @param linecol
#' @param margins
#' @param ColSideColors
#' @param RowSideColors
#' @param side.height.fraction
#' @param cexRow
#' @param cexCol
#' @param labRow
#' @param labCol
#' @param key
#' @param keysize
#' @param density.info
#' @param denscol
#' @param symkey
#' @param densadj
#' @param main
#' @param xlab
#' @param ylab
#' @param lmat
#' @param lhei
#' @param lwid
#' @param ColSideColorsSize
#' @param RowSideColorsSize
#' @param KeyValueName
#' @param ...
#'
#' @return
#' @export
#'
#' @examples
varbin_heatmap.3 <- 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,
                             symbreaks = max(x < 0, na.rm = TRUE) || scale != "none",
                             col = "heat.colors",
                             colsep,
                             rowsep,
                             sepcolor = "white",
                             sepwidth = c(0.05, 0.05),
                             cellnote,
                             notecex = 1,
                             notecol = "cyan",
                             na.color = par("bg"),
                             trace = c("none", "column","row", "both"),
                             tracecol = "cyan",
                             hline = median(breaks),
                             vline = median(breaks),
                             linecol = tracecol,
                             margins = c(5,5),
                             ColSideColors,
                             RowSideColors,
                             side.height.fraction=0.3,
                             cexRow = 0.2 + 1/log10(nr),
                             cexCol = 0.2 + 1/log10(nc),
                             labRow = NULL,
                             labCol = NULL,
                             key = TRUE,
                             keysize = 1.5,
                             density.info = c("none", "histogram", "density"),
                             denscol = tracecol,
                             symkey = max(x < 0, na.rm = TRUE) || symbreaks,
                             densadj = 0.25,
                             main = NULL,
                             xlab = NULL,
                             ylab = NULL,
                             lmat = NULL,
                             lhei = NULL,
                             lwid = NULL,
                             ColSideColorsSize = 1,
                             RowSideColorsSize = 1,
                             KeyValueName="Value",...){

  invalid <- function (x) {
    if (missing(x) || is.null(x) || length(x) == 0)
      return(TRUE)
    if (is.list(x))
      return(all(sapply(x, invalid)))
    else if (is.vector(x))
      return(all(is.na(x)))
    else return(FALSE)
  }

  x <- as.matrix(x)
  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) && (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 (Colv == "Rowv" && !isTRUE(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
  }
  retval$rowInd <- rowInd
  retval$colInd <- colInd
  retval$call <- match.call()
  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") {
    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, "/")
  }
  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
  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)) {
      #if (!is.matrix(ColSideColors))
      #stop("'ColSideColors' must be a matrix")
      if (!is.character(ColSideColors) || nrow(ColSideColors) != nc)
        stop("'ColSideColors' must be a matrix of nrow(x) rows")
      lmat <- rbind(lmat[1, ] + 1, c(NA, 1), lmat[2, ] + 1)
      #lhei <- c(lhei[1], 0.2, lhei[2])
      lhei=c(lhei[1], side.height.fraction*ColSideColorsSize/2, lhei[2])
    }

    if (!missing(RowSideColors)) {
      #if (!is.matrix(RowSideColors))
      #stop("'RowSideColors' must be a matrix")
      if (!is.character(RowSideColors) || ncol(RowSideColors) != nr)
        stop("'RowSideColors' must be a matrix of ncol(x) columns")
      lmat <- cbind(lmat[, 1] + 1, c(rep(NA, nrow(lmat) - 1), 1), lmat[,2] + 1)
      #lwid <- c(lwid[1], 0.2, lwid[2])
      lwid <- c(lwid[1], side.height.fraction*RowSideColorsSize/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))
  op <- par(no.readonly = TRUE)
  on.exit(par(op))

  layout(lmat, widths = lwid, heights = lhei, respect = FALSE)

  if (!missing(RowSideColors)) {
    if (!is.matrix(RowSideColors)){
      par(mar = c(margins[1], 0, 0, 0.5))
      image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    } else {
      par(mar = c(margins[1], 0, 0, 0.5))
      rsc = t(RowSideColors[,rowInd, drop=F])
      rsc.colors = matrix()
      rsc.names = names(table(rsc))
      rsc.i = 1
      for (rsc.name in rsc.names) {
        rsc.colors[rsc.i] = rsc.name
        rsc[rsc == rsc.name] = rsc.i
        rsc.i = rsc.i + 1
      }
      rsc = matrix(as.numeric(rsc), nrow = dim(rsc)[1])
      image(t(rsc), col = as.vector(rsc.colors), axes = FALSE)
      if (length(rownames(RowSideColors)) > 0) {
        axis(1, 0:(dim(rsc)[2] - 1)/max(1,(dim(rsc)[2] - 1)), rownames(RowSideColors), las = 2, tick = FALSE)
      }
    }
  }

  if (!missing(ColSideColors)) {

    if (!is.matrix(ColSideColors)){
      par(mar = c(0.5, 0, 0, margins[2]))
      image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    } else {
      par(mar = c(0.5, 0, 0, margins[2]))
      csc = ColSideColors[colInd, , drop=F]
      csc.colors = matrix()
      csc.names = names(table(csc))
      csc.i = 1
      for (csc.name in csc.names) {
        csc.colors[csc.i] = csc.name
        csc[csc == csc.name] = csc.i
        csc.i = csc.i + 1
      }
      csc = matrix(as.numeric(csc), nrow = dim(csc)[1])
      image(csc, col = as.vector(csc.colors), axes = FALSE)
      if (length(colnames(ColSideColors)) > 0) {
        axis(2, 0:(dim(csc)[2] - 1)/max(1,(dim(csc)[2] - 1)), colnames(ColSideColors), las = 2, tick = FALSE)
      }
    }
  }

  par(mar = c(margins[1], 0, 0, margins[2]))
  x <- t(x)
  cellnote <- t(cellnote)
  if (revC) {
    iy <- nr:1
    if (exists("ddr"))
      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, ...)
  retval$carpet <- x
  if (exists("ddr"))
    retval$rowDendrogram <- ddr
  if (exists("ddc"))
    retval$colDendrogram <- ddc
  retval$breaks <- breaks
  retval$col <- col
  if (!invalid(na.color) & any(is.na(x))) { # load library(gplots)
    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 = nrow(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")) {
    retval$vline <- vline
    vline.vals <- scale01(vline, min.scale, max.scale)
    for (i in 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 rowInd) {
      if (!is.null(hline)) {
        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 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)
    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(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 = tmpbreaks,
          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, KeyValueName, line = 2)
    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")
      par(cex = 0.5)
      mtext(side = 2, "Density", line = 2)
    }
    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")
      par(cex = 0.5)
      mtext(side = 2, "Count", line = 2)
    }
    else title("Color Key")
  }
  else plot.new()
  retval$colorTable <- data.frame(low = retval$breaks[-length(retval$breaks)],
                                  high = retval$breaks[-1], color = retval$col)
  invisible(retval)
}