R/cyto_stats_compute-methods.R

In DillonHammill/cytoSuite: Compensation, Gating & Visualisation Toolkit for Analysis of Flow Cytometry Data

#' Compute Statistics
#'
#' Calculate and export flow cytometry statistics.
#'
#' @param x object of class \code{\link[flowCore:flowFrame-class]{flowFrame}},
#'   \code{\link[flowCore:flowSet-class]{flowSet}} or
#'   \code{\link[flowWorkspace:GatingSet-class]{GatingSet}}.
#' @param ... additional method-specific arguments.
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @seealso \code{\link{cyto_stats_compute,flowFrame-method}}
#' @seealso \code{\link{cyto_stats_compute,flowSet-method}}
#' @seealso \code{\link{cyto_stats_compute,GatingSet-method}}
#'
#' @export
setGeneric(
  name = "cyto_stats_compute",
  def = function(x, ...) {
    standardGeneric("cyto_stats_compute")
  }
)

#' Compute Statistics - flowFrame Method
#'
#' Calculate and export flow cytometry statistics for a flowFrame.
#'
#' @param x object of class \code{\link[flowCore:flowFrame-class]{flowFrame}}.
#' @param channels names of of channels for which statistic should be
#'   calculated, set to all channels by default.
#' @param trans object of class
#'   \code{\link[flowCore:transformList-class]{transformList}} or
#'   \code{\link[flowWorkspace:transformerList]{transformerList}} generated by
#'   \code{\link[flowCore:logicleTransform]{estimateLogicle}} used to transform
#'   the fluorescent channels of x. The transformation list is required to apply
#'   the inverse transformation such that the required statistics are returned
#'   on the original linear scale.
#' @param stat name of the statistic to calculate, options include
#'   \code{"count"}, \code{"freq"}, \code{"median"}, \code{"mode"},
#'   \code{"mean"}, \code{"geo mean"} or \code{"CV"}.
#' @param density_smooth smoothing parameter passed to
#'   \code{\link[stats:density]{density}} when calculating mode, set to 1.5 by
#'   default.
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @importFrom flowCore exprs inverseLogicleTransform transform transformList
#' @importFrom stats median density sd
#'
#' @seealso \code{\link{cyto_stats_compute,flowSet-method}}
#'
#' @examples
#' library(CytoRSuiteData)
#' 
#' # Load in samples
#' fs <- Activation
#' 
#' # Apply compensation
#' fs <- compensate(fs, fs[[1]]@description$SPILL)
#' 
#' # Transform fluorescent channels
#' trans <- estimateLogicle(fs[[4]], cyto_fluor_channels(fs))
#' fs <- transform(fs, trans)
#' 
#' # Compute statistics
#' cyto_stats_compute(fs[[4]],
#'   channels = c("Alexa Fluor 647-A", "7-AAD-A"),
#'   trans = trans,
#'   stat = "median"
#' )
#' cyto_stats_compute(fs[[1]], stat = "count")
#' @export
setMethod(cyto_stats_compute,
  signature = "flowFrame",
  definition = function(x,
                          channels = NULL,
                          trans = NULL,
                          stat = "median",
                          density_smooth = 1.5) {

    # Check statistic
    stat <- .cyto_stat_check(stat = stat)

    # Assign x to fr
    fr <- x

    # Channels
    if (is.null(channels)) {
      channels <- BiocGenerics::colnames(fr)
    } else {
      channels <- cyto_channel_check(
        x = fr,
        channels = channels,
        plot = FALSE
      )
    }

    # Transformations
    if (is.null(trans) & stat %in%
      c("mean", "median", "mode", "geo mean", "CV", "CVI")) {
      message("'trans' missing - statistics will be returned on current scale.")
      trans <- NULL
    } else if (!is.null(trans)) {
      trans <- cyto_trans_check(trans, inverse = FALSE)
      inv <- cyto_trans_check(trans, inverse = TRUE)
      fr <- transform(fr, inv)
    }

    # Extract data from fr
    fr.exprs <- exprs(fr)

    # Statistics
    if (stat == "count") {
      cnt <- data.frame(count = BiocGenerics::nrow(fr))
      rownames(cnt) <- fr@description$GUID

      return(cnt)
    } else if (stat %in% c("median", "mode", "mean", "CV", "CVI")) {

      # Calculate statistic
      sts <- lapply(channels, function(channel) {
        if (stat == "mean") {
          mean(fr.exprs[, channel])
        } else if (stat == "median") {
          median(fr.exprs[, channel])
        } else if (stat == "mode") {
          d <- density(fr.exprs[, channel], adjust = density_smooth)
          d$x[d$y == max(d$y)]
        } else if (stat == "CV") {
          md <- median(fr.exprs[,channel])
          rSD <- median(abs(fr.exprs[,channel] - md))*1.4826
          rSD/md * 100
        }
      })
      names(sts) <- channels
      sts <- do.call(cbind, sts)
      rownames(sts) <- fr@description$GUID

      return(sts)
    } else if (stat == "geo mean") {

      # Don't transform calculate mean then inverse transform to get geo mean
      sts <- lapply(channels, function(channel) {
        fr.exprs <- exprs(x)[, channel]

        # trans supplied
        if (!is.null(trans)) {

          # Channel is already transformed
          if (channel %in% BiocGenerics::colnames(trans)) {
            inv <- cyto_trans_check(trans, inverse = TRUE)
            sts <- inv@transforms[[channel]]@f(mean(fr.exprs))

            # Channel has not been transformed
          } else if (!channel %in% BiocGenerics::colnames(trans)) {
            sts <- exp(mean(log(fr.exprs)))
          }

          # No trans supplied
        } else {
          sts <- suppressWarnings(exp(mean(log(fr.exprs))))

          if (is.nan(sts)) {
            stop(
              paste0(
                "Supply transformList/transformerList to calculate ",
                stat, "."
              )
            )
          }
        }

        return(sts)
      })
      sts <- do.call(cbind, sts)
      rownames(sts) <- fr@description$GUID
      colnames(sts) <- channels
    }

    return(sts)
  }
)

#' Compute Statistics - flowSet Method
#'
#' Calculate and export flow cytometry statistics for a flowSet.
#'
#' @param x object of class \code{\link[flowCore:flowSet-class]{flowSet}}.
#' @param channels names of of channels for which statistic should be
#'   calculated, set to all channels by default.
#' @param trans object of class
#'   \code{\link[flowCore:transformList-class]{transformList}} or
#'   \code{\link[flowWorkspace:transformerList]{transformerList}} generated by
#'   \code{\link[flowCore:logicleTransform]{estimateLogicle}} used to transform
#'   the fluorescent channels of x. The transformation list is required to apply
#'   the inverse transformation such that the required statistics are returned
#'   on the original linear scale.
#' @param stat name of the statistic to calculate, options include
#'   \code{"count"}, \code{"freq"}, \code{"median"}, \code{"mode"},
#'   \code{"mean"}, \code{"geo mean"} or \code{"CV"}.
#' @param density_smooth smoothing parameter passed to
#'   \code{\link[stats:density]{density}} when calculating mode, set to 1.5 by
#'   default.
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @importFrom flowCore exprs inverseLogicleTransform transform
#' @importFrom flowWorkspace pData
#'
#' @seealso \code{\link{cyto_stats_compute,flowFrame-method}}
#' @seealso \code{\link{cyto_stats_compute,GatingSet-method}}
#'
#' @examples
#' library(CytoRSuiteData)
#' 
#' # Load in samples
#' fs <- Activation
#' 
#' # Apply compensation
#' fs <- compensate(fs, fs[[1]]@description$SPILL)
#' 
#' # Transform fluorescent channels
#' trans <- estimateLogicle(fs[[4]], cyto_fluor_channels(fs))
#' fs <- transform(fs, trans)
#' 
#' # Compute statistics
#' cyto_stats_compute(fs,
#'   channels = c("Alexa Fluor 488-A", "PE-A"),
#'   trans = trans,
#'   stat = "median"
#' )
#' cyto_stats_compute(fs,
#'   trans = trans,
#'   stat = "geo mean"
#' )
#' cyto_stats_compute(fs,
#'   stat = "count"
#' )
#' @export
setMethod(cyto_stats_compute,
  signature = "flowSet",
  definition = function(x,
                          channels = NULL,
                          trans = NULL,
                          stat = "median",
                          density_smooth = 1.5) {

    # Check statistic
    stat <- .cyto_stat_check(stat = stat)

    # Assign x to fs
    fs <- x

    # pData
    pdata <- as.data.frame(pData(fs))

    # cyto_stats_compute
    sts <- fsApply(fs, function(fr) {
      cyto_stats_compute(fr,
        channels = channels,
        trans = trans,
        stat = stat,
        density_smooth = density_smooth
      )
    })

    # Structure of count statistics
    if (stat == "count") {
      sts <- do.call(rbind, sts)
    }

    # Merge with pdata
    pdata[, "name"] <- list(NULL)
    if (class(sts) == "data.frame") {
      sts <- lapply(list(sts), function(x) {
        cbind(pdata, x)
      })
    } else if (class(sts) == "matrix") {
      sts <- lapply(list(sts), function(x) {
        cbind(pdata, x)
      })
    }

    return(sts)
  }
)

#' Compute Statistics - GatingSet Method
#'
#' Calculate and export flow cytometry statistics for a GatingSet.
#'
#' @param x object of class
#'   \code{\link[flowWorkspace:GatingSet-class]{GatingSet}}.
#' @param alias name(s) of the population(s) for which the statistic should be
#'   calculated.
#' @param parent name(s) of the parent population(s) used calculate population
#'   frequencies. The frequency of alias in each parent will be returned as a
#'   percentage.
#' @param channels names of of channels for which statistic should be
#'   calculated, set to all channels by default.
#' @param trans object of class
#'   \code{\link[flowCore:transformList-class]{transformList}} or
#'   \code{\link[flowWorkspace:transformerList]{transformerList}} generated by
#'   \code{\link[flowCore:logicleTransform]{estimateLogicle}} used to transform
#'   the fluorescent channels of x. The transformation list is required to apply
#'   the inverse transformation such that the required statistics are returned
#'   on the original linear scale.
#' @param stat name of the statistic to calculate, options include
#'   \code{"count"}, \code{"freq"}, \code{"median"}, \code{"mode"},
#'   \code{"mean"}, \code{"geo mean"}, \code{"CV"}, or \code{"freq"}.
#' @param density_smooth smoothing parameter passed to
#'   \code{\link[stats:density]{density}} when calculating mode, set to 1.5 by
#'   default.
#' @param save logical indicating whether statistical results should be saved to
#'   a csv file.
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @importFrom flowCore exprs inverseLogicleTransform transform
#' @importFrom flowWorkspace getData sampleNames getTransformations
#' @importFrom utils write.csv
#'
#' @seealso \code{\link{cyto_stats_compute,flowFrame-method}}
#' @seealso \code{\link{cyto_stats_compute,flowSet-method}}
#'
#' @examples
#' library(CytoRSuiteData)
#' 
#' # Load in samples
#' fs <- Activation
#' gs <- GatingSet(fs)
#' 
#' # Apply compensation
#' gs <- compensate(gs, fs[[1]]@description$SPILL)
#' 
#' # Transform fluorescent channels
#' trans <- estimateLogicle(gs[[4]], cyto_fluor_channels(gs))
#' gs <- transform(gs, trans)
#' 
#' # Gate using gate_draw
#' gt <- Activation_gatingTemplate
#' gating(gt, gs)
#' 
#' # Compute statistics
#' cyto_stats_compute(gs,
#'   alias = "T Cells",
#'   channels = c("Alexa Fluor 488-A", "PE-A"),
#'   stat = "median",
#'   save = FALSE
#' )
#' 
#' cyto_stats_compute(gs,
#'   alias = "T Cells",
#'   stat = "geo mean",
#'   save = FALSE
#' )
#' 
#' cyto_stats_compute(gs,
#'   alias = c("CD4 T Cells", "CD8 T Cells"),
#'   stat = "count",
#'   save = FALSE
#' )
#' 
#' cyto_stats_compute(gs,
#'   alias = c("CD4 T Cells", "CD8 T Cells"),
#'   parent = c("Live Cells", "T Cells"),
#'   stat = "freq",
#'   save = FALSE
#' )
#' @export
setMethod(cyto_stats_compute,
  signature = "GatingSet",
  definition = function(x,
                          alias = NULL,
                          parent = NULL,
                          channels = NULL,
                          trans = NULL,
                          stat = "median",
                          density_smooth = 1.5,
                          save = TRUE) {

    # Check statistic
    stat <- .cyto_stat_check(stat = stat)

    # Assign x to gs
    gs <- x

    # transformerList?
    if (!is.null(trans)) {
      trans <- cyto_trans_check(trans, inverse = FALSE)
    }

    # Get trans if not supplied
    if (is.null(trans)) {
      trans <- transformList(
        names(getTransformations(gs[[1]])),
        getTransformations(gs[[1]])
      )
    }

    # Extract population(s)
    fs.lst <- lapply(alias, function(x) getData(gs, x))

    # Population frequencies
    if (stat == "freq") {
      sts <- lapply(fs.lst, function(fs) {
        as.data.frame(cyto_stats_compute(fs,
          channels = channels,
          trans = trans,
          stat = "count",
          density_smooth = density_smooth
        ))
      })
      names(sts) <- alias

      # Get parent counts
      prnts <- lapply(parent, function(x) {
        fs <- getData(gs, x)
        as.data.frame(cyto_stats_compute(fs,
          channels = channels,
          trans = trans,
          stat = "count",
          density_smooth = density_smooth
        ))
      })
      names(prnts) <- parent

      sts <- lapply(sts, function(x) {
        m <- matrix(rep(x[, "count"], length(parent)),
          ncol = length(parent),
          byrow = FALSE
        )
        colnames(m) <- rep("count", length(parent))
        sts <- cbind(x, m)
        return(sts)
      })

      sts <- lapply(sts, function(x) {
        lapply(seq_len(length(parent)), function(y) {
          N <- ncol(x) - length(parent)
          x[, (N + y)] <<- (x[, (N + y)] / prnts[[y]][, "count"]) * 100
          colnames(x)[(N + y)] <<- parent[y]
        })

        return(x)
      })
    } else {
      sts <- lapply(fs.lst, function(fs) {
        sts <- cyto_stats_compute(fs,
          channels = channels,
          trans = trans,
          stat = stat,
          density_smooth = density_smooth
        )
        clnms <- colnames(sts[[1]])
        sts <- as.data.frame(sts)
        colnames(sts) <- clnms
        return(sts)
      })
      names(sts) <- alias
    }

    # Save to csv file
    if (save == TRUE) {
      lapply(seq_along(alias), function(x) {
        write.csv(
          sts[[x]],
          paste(
            format(Sys.Date(), "%d%m%y"),
            "-", alias[x], "-", stat, ".csv",
            sep = ""
          )
        )
      })
    }

    return(sts)
  }
)