R/cyto_spillover_compute.R

In DillonHammill/CytoExploreR: Interactive Analysis of Cytometry Data

## CYTO_SPILLOVER_COMPUTE ------------------------------------------------------

#' Compute Spillover Matrix
#'
#' \code{cyto_spillover_compute} uses the method described by Bagwell & Adams
#' 1993 to automatically calculate the fluorescent spillover matrix using single
#' stain compensation controls.
#'
#' \code{cyto_spillover_compute} supports spillover matrix calculation for both
#' internal or universal unstained reference populations based on channel
#' selection. \code{cyto_spillover_compute} expects the fluorescent channels of
#' the samples to be pre-transformed. Attempts will be made to transform the
#' data internally (using biexponential transformations) if it looks like the
#' data has not been transformed.
#'
#' \code{cyto_spillover_compute} begins by the user selecting which fluorescent
#' channel is associated with each compensation control from a dropdown menu.
#' The results of these selections are saved to a csv file called
#' "Compensation-Channels.csv" which can be passed to the \code{channel_match}
#' argument on subsequent runs to bypass the channel selection process. In cases
#' where multiple controls are supplied for the same channel, the control with
#' the greatest signal in the designated channel (MedFI) will be used for the
#' calculation.
#'
#' Following channel selection, \code{cyto_spillover_compute} runs through each
#' compensation control and plots the density distribution in the associated
#' channel. If a universal "Unstained" compensation control is supplied, the
#' unstained compensation control will be overlaid onto the plot as a reference
#' for gating. Users can then gate the positive signal for spillover calculation
#' using an interval gate. If no universal unstained compensation control is
#' supplied, users are expected to gate the negative and then the positive
#' signal for each compensation control.
#'
#' The percentage spillover is calculated based on the median fluorescent
#' intensities (MedFI) of the positive populations relative to that of the
#' reference negative population(s). The calculated spillover matrix is returned
#' and written to a named csv file for future use.
#'
#' @param x object of class \code{\link[flowCore:flowSet-class]{flowSet}} or
#'   \code{\link[flowWorkspace:GatingSet-class]{GatingSet}} containing
#'   transformed and gated compensation controls.
#' @param parent name of the population to use for the spillover calculation
#'   when a GatingSet object is supplied, set to the last node of the GatingSet
#'   by default (e.g. "Single Cells"). For greater flexibility, users can
#'   specify a parent population for each control, which will be extracted for
#'   the spillover calculation (e.g. Lymphocytes for CD4 APC or Myeloid Cells
#'   for CD11b FITC). The parent populations for each control can also be
#'   specified in a \code{parent} column in the channel match CSV file or in
#'   \code{cyto_details}.
#' @param axes_trans object of class
#'   \code{\link[flowWorkspace:transformerList]{transformerList}} generated by a
#'   \code{cyto_transformer} which contains the transformer definitions used to
#'   transform the data. Transformer definitions are only required when a
#'   flowSet object is supplied.
#' @param channel_match name of csv file to associate a fluorescent channel with
#'   each of the compensation controls. The \code{channel_match} file should
#'   contain a "name" column with the names of the compensation controls and a
#'   "channel" column to associate a fluorescent channel with each compensation
#'   control. Users need not generate this file by hand as it will be created
#'   following the channel selection process. This information can also be added
#'   directly to the samples using \code{cyto_details_edit}.
#' @param spillover name of the output spillover csv file, set to
#'   \code{"Spillover-Matrix.csv"} by default.
#' @param axes_limits options include \code{"auto"}, \code{"data"} or
#'   \code{"machine"} to use optimised, data or machine limits respectively. Set
#'   to \code{"machine"} by default to use entire axes ranges. Fine control over
#'   axes limits can be obtained by altering the \code{xlim} and \code{ylim}
#'   arguments.
#' @param ... additional arguments passed to \code{\link{cyto_plot}}.
#'
#' @return spillover matrix and write spillover matrix to csv file named in
#'   accordance with \code{spillover}.
#'
#' @examples
#' \dontrun{
#' library(CytoExploreRData)
#'
#' # Bypass directory check for external files
#' options("CytoExploreR_wd_check" = FALSE)
#'
#' # Load in compensation controls
#' fs <- Compensation
#' gs <- GatingSet(Compensation)
#'
#' # Gate using cyto_gate_draw
#' gt <- Compensation_gatingTemplate
#' cyto_gatingTemplate_apply(gs, gt)
#'
#' # Channel match fille
#' cmfile <- system.file("extdata",
#'   "Compensation-Channels.csv",
#'   package = "CytoExploreRData"
#' )
#'
#' # Compute fluorescent spillover matrix
#' spill <- cyto_spillover_compute(cyto_extract(gs, "Single Cells"),
#'   channel_match = cmfile,
#'   spillover = "Example-spillover.csv"
#' )
#'
#' # Compensate samples
#' gs <- cyto_compensate(gs, spill)
#'
#' # Return CytoExploreR_wd_check to default
#' options("CytoExploreR_wd_check" = TRUE)
#' }
#'
#' @importFrom flowCore each_col fsApply sampleNames flowSet Subset
#' @importFrom flowWorkspace pData GatingSet
#' @importFrom methods as is
#' @importFrom utils read.csv write.csv
#' @importFrom stats median
#' @importFrom tools file_ext
#'
#' @seealso \code{\link{cyto_spillover_edit}}
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @references C. B. Bagwell \& E. G. Adams (1993). Fluorescence spectral
#'   overlap compensation for any number of flow cytometry parameters. in:
#'   Annals of the New York Academy of Sciences, 677:167-184.
#'
#' @export
cyto_spillover_compute <- function(x,
                                   parent = NULL,
                                   axes_trans = NULL,
                                   channel_match = NULL,
                                   spillover = NULL,
                                   axes_limits = "machine",
                                   ...) {

  # PREPARE DATA ---------------------------------------------------------------
  
  # COPY
  cyto_copy <- cyto_copy(x)
  
  # CHANNELS 
  channels <- cyto_fluor_channels(cyto_copy)
  
  # TRANSFORMATIONS
  axes_trans <- cyto_transformer_extract(cyto_copy)
  
  # DATA SHOULD BE TRANSFORMED FOR GATING
  if(.all_na(axes_trans) | is.null(axes_trans)){
    axes_trans <- cyto_transformer_biex(cyto_copy,
                                        channels = channels,
                                        type = "biex",
                                        plot = FALSE)
    suppressWarnings(cyto_transform(cyto_copy,
                                    trans = axes_trans,
                                    plot = FALSE))
  }

  # SPILLOVER COMPUTATION ------------------------------------------------------

  # NEGATIVE AND POSITIVE POPULATIONS (TRANSFORMED)
  pops <- .cyto_spillover_pops(cyto_copy,
                               parent = parent,
                               channel_match = channel_match,
                               axes_trans = axes_trans,
                               axes_limits = axes_limits,
                               ...)
  neg_pops <- pops[["negative"]]
  pos_pops <- pops[["positive"]]
  
  # UPDATE DETAILS IN X
  cyto_details(x) <- cyto_details(cyto_copy)
  
  # EXPERIMENT DETAILS
  pd <- cyto_details(pos_pops)
  
  # SAMPLE NAMES
  nms <- cyto_names(pos_pops)
  
  # MEDFI ALL CHANNELS
  neg_stats <- list()
  lapply(seq_along(neg_pops), function(z){
    # CHECK PREVIOUS - PREVENT DUPLICATE COMPUTATION
    if(z > 1 & any(LAPPLY(seq_len(z-1), function(y){
      identical(neg_pops[[z]], neg_pops[[y]])
    }))){
      neg_stats[[z]] <<- neg_stats[[which(LAPPLY(seq_len(z-1), function(y){
        identical(neg_pops[[z]], neg_pops[[y]])
      }))[1]]]
    }else{
      neg_stats[[z]] <<- suppressMessages(
        cyto_stats_compute(neg_pops[[z]],
                           stat = "median",
                           channels = channels,
                           format = "wide",
                           trans = axes_trans)
      )
    }
  })
  neg_stats <- do.call("rbind", neg_stats)
  neg_stats <- neg_stats[, -which(names(neg_stats) %in% colnames(pd))]
  colnames(neg_stats) <- channels
  neg_stats <- neg_stats[, channels]
  neg_stats <- data.matrix(neg_stats)
  rownames(neg_stats) <- nms
  
  # Calculate medFI for all channels in all stained controls
  pos_stats <- suppressMessages(cyto_stats_compute(pos_pops,
                                             stat = "median",
                                             channels = channels,
                                             format = "wide",
                                             trans = axes_trans
  ))
  pos_stats <- pos_stats[, -which(names(pos_stats) %in% colnames(pd))]
  colnames(pos_stats) <- channels
  pos_stats <- pos_stats[, channels]
  pos_stats <- data.matrix(pos_stats)
  rownames(pos_stats) <- nms
  
  # Subtract background fluorescence
  signal <- pos_stats - neg_stats
  signal <- as.matrix(signal)

  # Construct spillover matrix - include values for which there is a control
  spill <- diag(x = 1, 
                nrow = length(channels), 
                ncol = length(channels))
  colnames(spill) <- channels
  rownames(spill) <- channels

  # Normalise each row to stained channel
  lapply(seq(1, nrow(signal), 1), function(z) {
    signal[z, ] <<- signal[z, ] /
      signal[z, match(pd$channel[z], colnames(spill))]
  })
  
  # Insert values into appropriate rows
  rws <- match(pd$channel, rownames(spill))
  spill[rws, ] <- signal

  # No name specified for spillover default to date-Spillover-Matrix.csv
  if (is.null(spillover)) {
    spillover <- paste0(format(Sys.Date(), "%d%m%y"), "-Spillover-Matrix.csv")
  }

  # write spillover matrix to csv file
  if (!is(spillover, "character")) {
    stop("'spillover' should be the name of a csv file.")
  } else {
    spillover <- file_ext_append(spillover, ".csv")
    write.csv(spill, spillover)
  }

  return(spill)
}