R/basesubtract_heatmap.R

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

#' An attaempt to create a basesubtraction option first step at calculating lifetime sparsness
#' @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 basesubtract Whether or not to subtract the baseline firing rate
#' @param baselinesamples Integer indices of the timepoints to use for the
#'   baseline. These are 50ms time points for \code{smSpikes}
#' @param ... Additional parameters passed to \code{\link{heatmap}} function
#' @export
basesubtract_heatmap_cor_dist<-function(cells, odours, col=jet.colors(20),
                                        basesubtract=TRUE, baselinesamples=1:5, ...) {
    # 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 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))

    if (basesubtract==T) {
      allfreqs=baseline_subtract_allfreqs(allfreqs, baselinesamples=baselinesamples)
    }

    # 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
    spcor=cor(allfreqs_odours_matrix,use='pairwise.complete.obs')
    # remove cells with no spikes
    noSpikes=names(which(is.na(spcor[, 1])))
    spcor=spcor[!row.names(spcor)%in%noSpikes, !colnames(spcor)%in%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)
    heatmap(spcor,distfun=function(x) as.dist(1-x),scale='none',symm=T, col=col, ...)
}

#' Subtract baseline spike rate from list of smoothed psth data
#'
#' @param allfreqs A list with one entry per cell, where each entry contains a
#'   matrix of smoothed psth data where rows are timepoints and columns are
#'   odours. This is equivalent to pulling the freq component out of
#'   \code{physplitdata::\link{smSpikes}} object.
#' @inheritParams basesubtract_heatmap_cor_dist
#' @examples
#' cells=c("nm20140714c1", "nm20150305c0", "nm20141128c0", "nm20140901c0")
#' allfreqs=lapply(physplitdata::smSpikes[cells],
#' function(psthsforcell) sapply(psthsforcell,function(x) x$freq))
#' baseline_subtract_allfreqs(allfreqs, baselinesamples=1:5)
#' @export
baseline_subtract_allfreqs <- function (allfreqs, baselinesamples=1:5) {
  cellbaseline=list()
  allfreqs_base=list()

  cells=names(allfreqs)

  for (cell in cells) {
    baseline_psth=allfreqs[[cell]][baselinesamples,]
    cellbaseline[[cell]]=colMeans(baseline_psth)
    allfreqs_base[[cell]]=scale(allfreqs[[cell]],
                                center=cellbaseline[[cell]], scale=FALSE)
  }

  # check that the mean is correctly subtracted by recalculating the mean baseline
  # (should be zero)
  cellbaseline_rebase=list()
  for (cell in cells) {
    cellbaseline_rebase[[cell]]=colMeans(allfreqs_base[[cell]][baselinesamples,])
  }
  stopifnot(all(abs(cellbaseline_rebase[[cell]])<1e-5))

  allfreqs_base
}