R/heatmap.R

In sfrechter/physplit.analysis: Packaged and versioned Analysis functions for Shahar's cells

#' heatmap for set of cells and odours based on correlation distance
#'
#' @param cells Character specifying cells to plot OR a data.frame containing
#'   information about those cells
#' @param odours Character vector specifying odours to plot
#' @param col Character vector of colour levels generated by a colour palette
#'   function such as \code{\link{jet.colors}}. Defaults to
#'   \code{jet.colors(20)}.
#' @param labCol character vectors with column labels to use; defaults to cells.
#' @param heatmapfun Which heatmap function to use. There are many besides the
#'   default \code{heatmap}. See examples.
#' @inheritParams stats::heatmap
#' @param ... Additional parameters passed to \code{\link{heatmap}} function
#' @export
#' @importFrom stats heatmap as.dist cor
#' @examples
#' # Find LHNs with Anatomy.type 4 and ephys class 16,17,23 for which we have
#' # spiking response data
#' physplit.c161723=subset(PhySplitDB, cell%in% names(Spikes) &
#'   Anatomy.type==4 & class%in%c(16,17,23))
#' heatmap_cor_dist(physplit.c161723)
#'
#' # repeat but with class information
#' heatmap_cor_dist(physplit.c161723, labRow=physplit.c161723$class,
#'   labCol=NA, RowSideColors=rainbow(3)[factor(physplit.c161723$class)])
#'
#' # Same but set limit for the palette
#' heatmap_cor_dist(physplit.c161723, labRow=physplit.c161723$class,
#'   labCol=NA, RowSideColors=rainbow(3)[factor(physplit.c161723$class)],
#'   zlim=c(-1,1))
#'
#' ## use heatmap.2
#' \dontrun{
#' library(gplots)
#' heatmap_cor_dist(physplit.c161723, heatmapfun=heatmap.2)
#' }
heatmap_cor_dist<-function(cells, odours, col=jet.colors(20), labRow=NULL,
                           labCol=NULL, ColSideColors, RowSideColors,
                           heatmapfun=heatmap,...) {
  spcor=spike_cor_dist(cells, odours)
  # remove cells with no spikes
  noSpikes=names(which(is.na(spcor[, 1])))
  if(length(noSpikes)){
    warning("There are no spikes for these odours for these cells:\n",
            paste(noSpikes, collapse=' '))
    spcor=spcor[!row.names(spcor)%in%noSpikes, !colnames(spcor)%in%noSpikes]
    # now we also need to fix the row/col labels/colours since they will be out of step
    # after we drop some cells
    labRow=labRow[cells!=noSpikes]
    labCol=labCol[cells!=noSpikes]

    if(!missing(ColSideColors))
      ColSideColors=ColSideColors[cells!=noSpikes]
    if(!missing(RowSideColors))
      RowSideColors=RowSideColors[cells!=noSpikes]
  }

  # If there are no odours with spikes in common between 2 cells
  # then the correlation score will still be NA
  # set it to 0 instead so that we can at least use those data in clustering
  spcor[is.na(spcor)]=0

  # The Heatmap!
  # dendrogram is based on distance of 1-correlation score, but the
  # colours in the heatmap are still the correlation scores (ie hot is highly correlated)
  heatmapfun(spcor,distfun=function(x) as.dist(1-x),scale='none',symm=T, col=col,
          labRow=labRow, labCol=labCol, ColSideColors=ColSideColors,
          RowSideColors=RowSideColors, ...)
}

#' @rdname heatmap_cor_dist
#' @export
#' @description \code{spike_cor_dist} calculate the cross-correlation between
#'   the spike responses of a set of cells. It is called by
#'   \code{heatmap_cor_dist} but you may wish to use this data directly.
spike_cor_dist <- function(cells, odours) {
  # First get the database info we need.
  # End up with a data.frame in the order of cells
  if(is.data.frame(cells)) {
    physplit=cells
    cells=physplit$cell
  } else {
    physplit=physplitdata::PhySplitDB[match(cells, physplitdata::PhySplitDB$cell),]
  }
  rownames(physplit)=physplit$cell
  if(length(missing_cells<-setdiff(physplit$cell, names(physplitdata::smSpikes)))) {
    stop("Error: some cells without spike data: ", paste(missing_cells, collapse = " "))
  }

  # if odours not specified, then use all odours
  if(missing(odours)) odours = unique(unlist(sapply(physplitdata::smSpikes,names)))

  # now collect responses that we need
  allfreqs=lapply(physplitdata::smSpikes[cells],
                  function(psthsforcell) sapply(psthsforcell,function(x) x$freq))
  # pad those frequencies with columns of NAs for missing odours
  # also reorder odours into the order given by odours
  # and use sapply => matrix output at the end
  allfreqs_odours_matrix=sapply(allfreqs,addnacols,odours)

  # Calculate correlation distance between responses for all our cells
  cor(allfreqs_odours_matrix,use='pairwise.complete.obs')
}


#' heatmap for set of cells on nblast anatomy distance
#' @param x A score matrix calculated by
#'   \code{\link[nat.nblast]{nblast_allbyall}}
#' @inheritParams stats::heatmap
#' @inheritParams heatmap_cor_dist
#' @importFrom nat.nblast sub_dist_mat
#' @export
#' @importFrom stats hclust
heatmap_anatomy <- function(x, col=jet.colors(20), labRow = NULL,
                            labCol = NULL, ColSideColors, RowSideColors,
                            heatmapfun=heatmap, ...) {
  heatmapfun(1-sub_dist_mat(scoremat = x), distfun = function(x) as.dist(1-x),
          hclustfun = function(x, ...) hclust(x, method='ward.D',...), scale = "none",
          symm = T, col = col, labRow = labRow, labCol = labCol,
          ColSideColors = ColSideColors, RowSideColors = RowSideColors,
          ...)
}