R/singletGate.R

In flowStats: Statistical methods for the analysis of flow cytometry data

#' Creates a singlet polygon gate using the prediction bands from a robust linear model
#'
#' We construct a singlet gate by applying a robust linear model with
#' \code{\link[MASS]{rlm}}. By default, we model the forward-scatter height
#' (FSC-H)as a function of forward-scatter area (FSC-A). If \code{sidescatter}
#' is given, forward-scatter height is as a function of \code{area} +
#' \code{sidescatter} + \code{sidescatter / area}.
#'
#' Because \code{\link[MASS]{rlm}} relies on iteratively reweighted least
#' squares (IRLS), the runtime to construct a singlet gate is dependent in part
#' on the number of observations in \code{x}. To improve the runtime, we provide
#' an option to subsample randomly a subset of \code{x}. A percentage of
#' observations to subsample can be given in \code{subsample_pct}. By default, no
#' subsampling is applied.
#' 
#' @param x a \code{\link[flowCore:flowFrame-class]{flowFrame}} object
#' @param area character giving the channel name that records the signal
#' intensity as peak area
#' @param height character giving the channel name that records the signal
#' intensity as peak heightchannel name of height
#' @param sidescatter character giving an optional channel name for the
#' sidescatter signal. By default, ignored.
#' @param prediction_level a numeric value between 0 and 1 specifying the level
#' to use for the prediction bands
#' @param subsample_pct a numeric value between 0 and 1 indicating the percentage
#' of observations that should be randomly selected from \code{x} to construct
#' the gate. By default, no subsampling is performed.
#' @param wider_gate logical value. If \code{TRUE}, the prediction bands used to
#' construct the singlet gate use the robust fitted weights, which increase
#' prediction uncertainty, especially for large FSC-A. This leads to wider gates,
#' which are sometimes desired.
#' @param filterId the name for the filter that is returned
#' @param maxit the limit on the number of IWLS iterations passed to \code{\link[MASS]{rlm}}
#' @param ... additional arguments passed to \code{\link[MASS]{rlm}}
#' @return a \code{\link[flowCore]{polygonGate}} object with the singlet gate
#' @rdname gate_singlet
#' @export 
gate_singlet <- function(x, area = "FSC-A", height = "FSC-H", sidescatter = NULL,
                        prediction_level = 0.99, subsample_pct = NULL,
                        wider_gate = FALSE, filterId = "singlet", maxit = 5, ...) {
  flowCore:::checkClass(x, "flowFrame")
  flowCore:::checkClass(area, "character")
  flowCore:::checkClass(height, "character")
  if (!is.null(sidescatter)) {
    flowCore:::checkClass(sidescatter, "character")
  }
  if (length(area) + length(height) != 2) {
    stop("Each of 'area' and 'height' must be 'character' vectors of length 1.")
  }

  x <- data.frame(exprs(x[, c(area, height, sidescatter)]))

  # If specified, subsample from 'x'
  if (!is.null(subsample_pct)) {
    subsample_pct <- as.numeric(subsample_pct)
    if (subsample_pct <= 0 || subsample_pct > 1) {
      warning("The subsampling percentage must be between 0 and 1. ",
              "Setting percentage to default value...", call. = FALSE)
      subsample_pct <- 1
    }
    n <- nrow(x)
    x <- x[sample(x = seq_len(nrow(x)), size = subsample_pct * n), ]
  }
  
  channel_names <- c(area, height)
  area <- make.names(area)
  height <- make.names(height)
  rlm_formula <- paste(make.names(height), make.names(area), sep = " ~ ")
  if (!is.null(sidescatter)) {
    sidescatter <- make.names(sidescatter)
    ssc_ratio <- paste0("I(", sidescatter, " / ", area, ")")
    rlm_formula <- paste(rlm_formula, sidescatter, ssc_ratio, sep = " + ")
  }
  rlm_formula <- as.formula(rlm_formula)

  rlm_fit <- withCallingHandlers(rlm(rlm_formula, data = x, maxit = maxit, ...),
                                 warning = function(w){
                                   if(grepl("failed to converge", conditionMessage(w))){
                                     invokeRestart("muffleWarning")
                                   }
                                 })

  # if (!rlm_fit$converged) {
  #   warning("The IRLS algorithm employed in 'rlm' did not converge.")
  # }

  # Creates polygon gate based on the prediction bands at the minimum and maximum
  # 'area' observation using the trained robust linear model.
  which_min <- which.min(x[[area]])
  which_max <- which.max(x[[area]])
  x_extrema <- x[c(which_min, which_max), ]

  # Prediction weights. By default, these are weighted equally. If 'wider_gate'
  # is set, the weights from 'rlm' are used. The weight for the maximum FSC-A
  # is usually much smaller and yields a more uncertain prediction interval.
  # This leads to a wider gate.
  prediction_weights <- c(1, 1)
  if (wider_gate) {
    prediction_weights <- rlm_fit$w[c(which_min, which_max)]
  }

  predictions <- predict(rlm_fit, x_extrema, interval = "prediction",
                         level = prediction_level, weights = prediction_weights)

  # Create a matrix of the vertices using the prediction bands at the minimum
  # and maximum values of x. The ordering matters. Must go clockwise.
  # Otherwise, the polygon is not convex and makes an X-shape.
  gate_vertices <- rbind(cbind(x_extrema[[area]][1], predictions[1, "lwr"]),
                         cbind(x_extrema[[area]][1], predictions[1, "upr"]),
                         cbind(x_extrema[[area]][2], predictions[2, "upr"]),
                         cbind(x_extrema[[area]][2], predictions[2, "lwr"]))
  colnames(gate_vertices) <- channel_names

  polygonGate(gate_vertices, filterId = filterId)
}

#' @rdname gate_singlet
#' @export 
#' @examples 
#' \dontrun{
#'  # fr is a flowFrame
#'  sg <- gate_singlet(fr, area = "FSC-A", height = "FSC-H")
#'  sg
#'  # plot the gate 
#'  xyplot(`FSC-H` ~ `FSC-A`, fr, filter = sg)
#' }
singletGate <- gate_singlet