CytoExploreR: Interactive Analysis of Cytometry Data

Documented in cyto_spillover_spread_compute

## CYTO_SPILLOVER_SPREAD_COMPUTE -----------------------------------------------

#' Compute Spillover Spreading Matrix
#'
#' \code{cyto_spillover_spread_compute} computes the spillover spreading matrix
#' as described by Nguyen et al. (2013) using gated and compensated compensation
#' controls.
#'
#' @param x object of class \code{\link[flowCore:flowSet-class]{flowSet}} or
#'   \code{\link[flowWorkspace:GatingSet-class]{GatingSet}} containing
#'   compensated, transformed and gated compensation controls.
#' @param parent name of the population to use for downstream calculations when
#'   a \code{GatingSet} object is supplied, set to the last gated population by
#'   default (e.g. "Single Cells").
#' @param axes_trans object of class
#'   \code{\link[flowWorkspace:transformerList]{transformerList}} generated by a
#'   \code{cyto_transformer} which contains the transformer definitions to be
#'   used internally to transform the fluorescent channels of the samples as
#'   required.
#' @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.
#' @param compensated logical indicating whether the supplied compensation
#'   controls have been compensated. Compensation controls must be compensated
#'   prior to calculation of spillover spreading matrix. If \code{compensated}
#'   is set to FALSE, the compensation controls will be compensated internally
#'   using the supplied \code{spillover} matrix.
#' @param spillover name of the output spillover matrix csv file to be used
#'   internally to compensate the compensation controls. If no spillover matrix
#'   is supplied \code{spillover_spread_compute}, the spillover matrix attached
#'   to compensation controls will be applied.
#' @param spillover_spread name of the csv file to which the spillover spreading
#'   matrix will be saved, set to \code{"date-Spillover-Spread-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}}.
#'
#' @references Nguyen R, Perfetto S, Mahnke YD, Chattopadhyay P & Roederer M,
#'   (2013). Quantifying spillover spreading for comparing instrument
#'   performance and aiding in multicolor panel design. Cytometry A,
#'   83(3):306-15.
#'
#' @importFrom flowWorkspace pData
#' @importFrom flowCore Subset parameters exprs
#' @importFrom stats quantile na.omit
#' @importFrom grDevices graphics.off dev.new
#' @importFrom methods is
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @seealso \code{\link{cyto_spillover_compute}}
#' @seealso \code{\link{cyto_spillover_edit}}
#'
#' @export
cyto_spillover_spread_compute <- function(x,
                                          parent = NULL,
                                          axes_trans = NULL,
                                          channel_match = NULL,
                                          compensated = FALSE,
                                          spillover = NULL,
                                          spillover_spread = 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)

  # GATINGSET COMPENSATION
  if(is(cyto_copy, "GatingSet")){
    spill <- cyto_spillover_extract(cyto_copy)
    if(length(spill) == 0){
      compensated <- FALSE
    }else{
      compensated <- TRUE
    }
  }
  
  # DATA SHOULD BE COMPENSATED & TRANSFORMED
  if(compensated == FALSE){
    # UNTRANSFORMED DATA
    if(.all_na(axes_trans) | is.null(axes_trans)){
      cyto_copy <- cyto_compensate(cyto_copy, spillover = spillover)
      axes_trans <- cyto_transformer_biex(cyto_copy,
                                          channels = channels,
                                          type = "biex",
                                          plot = FALSE)
      cyto_copy <- suppressWarnings(cyto_transform(cyto_copy,
                                                   trans = axes_trans,
                                                   plot = FALSE))
    # TRANSFORMED DATA
    }else{
      # INVERSE TRANSFORM
      suppressWarnings(cyto_transform(cyto_extract(cyto_copy),
                                      trans = axes_trans,
                                      inverse = TRUE,
                                      plot = FALSE))
      # COMPENSATE
      cyto_compensate(cyto_copy, spillover = spillover)
      # TRANSFORM
      suppressWarnings(cyto_transform(cyto_extract(cyto_copy),
                                      trans = axes_trans,
                                      plot = FALSE))
    }
  }else if(compensated == TRUE){
    # TRANSFORMATIONS
    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 SPREAD 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)
  
  # INVERSE TRANSFORMATIONS
  neg_pops <- cyto_transform(neg_pops,
                             trans = axes_trans,
                             inverse = TRUE,
                             plot = FALSE)
  pos_pops <- cyto_transform(pos_pops,
                             trans = axes_trans,
                             inverse = TRUE,
                             plot = FALSE)
  
  # Calculate spillover spread for each compensation control
  SSM <- lapply(seq_along(neg_pops), function(z){
      
    # Extract NIL and pop
    NIL <- neg_pops[[z]]
    pop <- pos_pops[[z]]
      
    # Channel
    chan <- pd$channel[pd$name == cyto_names(pop)]
      
    # Calculate 50th and 84th percentile for NIL
    NIL_50 <- LAPPLY(channels, function(channel){
      quantile(exprs(NIL)[,channel], 0.5)
    })
    names(NIL_50) <- channels
    
    NIL_84 <- LAPPLY(channels, function(channel){
      quantile(exprs(NIL)[,channel], 0.84)
    })
    names(NIL_84) <- channels
    
    # Combine into a list - NIL_stats
    NIL_stats <- list("50%" = NIL_50, "84%" = NIL_84)
    
    # Calculate 50th and 84th percentile for pop
    pop_50 <- LAPPLY(channels, function(channel){
      quantile(exprs(pop)[,channel], 0.5)
    })
    names(pop_50) <- channels
    
    pop_84 <- LAPPLY(channels, function(channel){
      quantile(exprs(pop)[,channel], 0.84)
    })
    names(pop_84) <- channels
    
    # Combine into a list - pop_stats
    pop_stats <- list("50%" = pop_50, "84%" = pop_84)
    
    # Calculate pop SD and VARIANCE
    pop_SD <- pop_stats[[2]] - pop_stats[[1]]
    pop_VAR <- pop_SD^2
    
    # Calculate NIL SD and VARIANCE
    NIL_SD <- NIL_stats[[2]] - NIL_stats[[1]]
    NIL_VAR <- NIL_SD^2
    
    # Difference in VAR of NIL and pop in other detectors
    VAR_diff <- pop_VAR - NIL_VAR
    
    # Replace negative numbers with 0 - stain has less spread than reference
    if(any(VAR_diff < 0)){
      VAR_diff[VAR_diff < 0] <- 0
    }
    
    # Difference in SD of NIL and pop in other detectors
    SD_diff <- sqrt(VAR_diff)

    # Change in fluorescence in primary detector
    signal <- pop_stats[[1]][match(chan, channels)] - 
      NIL_stats[[1]][match(chan, channels)]

    # Spillover spread values
    return(SD_diff/sqrt(signal))
      
  })
    
  # Make matrix
  SSM <- do.call("rbind", SSM)
    
  # Add rownames (channel associated with each control)
  rownames(SSM) <- pd$channel[match(cyto_names(pos_pops), pd$name)]
  
  # Replace diagonal with NA
  lapply(seq(1, nrow(SSM), 1), function(z) {
    SSM[z, match(rownames(SSM)[z], colnames(SSM))] <<- NA
  })
    
  # Sort SSM by column names
  ind <- na.omit(match(colnames(SSM), rownames(SSM)))
  SSM <- SSM[ind, ]
    
  # Turn off pop-up graphics device (gating)
  graphics.off()
  dev.new()
    
  # Heatmap label size
  ylab_width <- max(nchar(rownames(SSM))) * 0.02
  xlab_width <- max(nchar(colnames(SSM))) * 0.04  
  
  # Plot heatmap
  superheat::superheat(SSM,
                       title = "Spillover Spreading Matrix \n",
                       title.size = 6,
                       title.alignment = "center",
                       left.label.size = ylab_width,
                       left.label.col = "white",
                       left.label.text.alignment = "right",
                       bottom.label.text.angle = 90,
                       bottom.label.size = xlab_width,
                       bottom.label.col = "white",
                       bottom.label.text.alignment = "right",
                       legend.height = 0.2,
                       legend.vspace = 0,
                       row.title = "Fluorochrome \n",
                       row.title.size = 6,
                       heat.na.col = "grey",
                       heat.pal = c("white",
                                    "steelblue",
                                    "steelblue2",
                                    "steelblue3",
                                    "steelblue4",
                                    "purple",
                                    "mediumorchid",
                                    "orchid",
                                    "magenta",
                                    "deeppink2",
                                    "deeppink4",
                                    "red"),
                       heat.pal.values = c(0,
                                           1,
                                           2,
                                           3,
                                           6,
                                           10,
                                           20,
                                           30,
                                           65,
                                           100,
                                           200,
                                           300))
  
  # Default spillover spread file name
  if(is.null(spillover_spread)){
    spillover_spread <- paste0(format(Sys.Date(), "%d%m%y"),
                               "-Spillover-Spread.csv")
  }
  
  # Write to csv file
  if (!is(spillover_spread, "character")) {
    stop("'spillover_spread' should be the name of a csv file.")
  } else {
    spillover_spread <- file_ext_append(spillover_spread, ".csv")
    write.csv(SSM, spillover_spread)
  }

  return(SSM)
}