R/pca.R

In CytoTree: A Toolkit for Flow And Mass Cytometry Data

#'
#' Calculate principal components in CYT
#'
#' @name runFastPCA
#'
#' @param object a CYT object
#' @param center logical, a logical value indicating whether the variables
#'    should be shifted to be zero centered. Alternately, a vector
#'    of length equal the number of columns of x can be supplied.
#'    The value is passed to scale. See \code{\link[gmodels]{fast.prcomp}}
#' @param scale. logical, a logical value indicating whether the
#'    variables should be scaled to have unit variance before the
#'    analysis takes place. The default is FALSE for consistency
#'    with S, but in general scaling is advisable. Alternatively,
#'    a vector of length equal the number of columns of x can be supplied.
#'    The value is passed to scale. See \code{\link[gmodels]{fast.prcomp}}
#' @param verbose logical. Whether to print calculation progress.
#' @param ... Parameters passing to \code{\link[gmodels]{fast.prcomp}} function
#'
#' @importFrom gmodels fast.prcomp
#' @seealso \code{\link[gmodels]{fast.prcomp}}
#'
#' @export
#' @return A CYT object with PCA
#' @examples
#'
#' cyt.file <- system.file("extdata/cyt.rds", package = "CytoTree")
#' cyt <- readRDS(file = cyt.file)
#' 
#' cyt <- runFastPCA(cyt, verbose = TRUE)
#' 
#'
runFastPCA <- function(object, center = FALSE, scale. = TRUE,
                       verbose = FALSE, ...) {
  # PCA calculation
  if (verbose) message(Sys.time(), " Calculating PCA.")
  if (length(which(object@meta.data$dowsample == 1)) < 10) stop(Sys.time, " Not enough cells, please run processingCluster and choose correct downsampleing.size paramter. ")
  mat <- object@log.data[which(object@meta.data$dowsample == 1), object@markers.idx]
  pca.obj <- fast.prcomp( t(mat), retx = TRUE, center = center, scale. = scale., ...)

  object@pca.sdev <- pca.obj$sdev
  object@pca.value <- pca.obj$rotation
  object@pca.scores <- pca.obj$x

  colnames(object@pca.value) <- paste0("PC_", seq_len(ncol(object@pca.value)))

  if (verbose) message(Sys.time(), " Calculating PCA completed. ")

  return(object)
}