R/runDiff.R

In JhuangLab/CytoTree: A Toolkit for Flow And Mass Cytometry Data

#'
#' Calculate differential expression markers
#'
#' @name runDiff
#'
#' @description Calculating differentially expressed markers
#'
#' @param object a CYT object
#' @param branch.id vector. Branch ids use to run differentially expressed markers
#' @param branch.id.2 vector. Branch ids use to run differentially expressed
#'    markers in compare with branch.id
#' @param verbose logic. Whether to print calculation progress.
#'
#' @seealso  \code{bulidTree}
#'
#' @return A CYT object with cluster.id in meta.data
#'
#' @import limma
#' @importFrom stringr str_replace_all fixed
#' @importFrom stats model.matrix
#'
#' @export
#' @return a data.frame with differential expressed markers
#'
#' @examples
#'
#' cyt.file <- system.file("extdata/cyt.rds", package = "CytoTree")
#' cyt <- readRDS(file = cyt.file)
#'
#' DEG.table <- runDiff(cyt)
#'
#' 
#'
#'
#'
runDiff <- function(object, branch.id = NULL, branch.id.2 = NULL, verbose = FALSE) {

  if (verbose) message(Sys.time(), " Calculating differentially expressed markers.")
  if (missing(object)) stop(Sys.time(), " CYT object is missing.")
  if (!"branch.id" %in% colnames(object@meta.data)) stop(Sys.time(), " branch.id is missing, please run buildTree first.")

  all.branch.ids <- unique(object@meta.data$branch.id)

  total.deg.list <- NULL
  branch.contrast <- NULL
  ga <- go <- NULL
  if (length(all.branch.ids) == 1) {
    stop(Sys.time(), " There is only one branch in the tree.")
  } else {
    pdata <- object@meta.data[which(object@meta.data$dowsample == 1), c("cell", "branch.id")]
    edata <- object@log.data[which(object@meta.data$dowsample == 1), object@markers.idx]
    if (is.null(branch.id) & is.null(branch.id.2)) {
      if (verbose) message(Sys.time(), " All branches will be calculated.")
      for (bid in all.branch.ids) {
        pdata$contrast <- "go"
        pdata$contrast[which(pdata$branch.id == bid)] = "ga"
        design <- model.matrix(~ 0 + as.factor(contrast), data = pdata)
        colnames(design) <- stringr::str_replace_all(colnames(design), fixed("as.factor(contrast)"), "")
        fit <- lmFit(t(edata), design)
        contrast <- makeContrasts(ga_go = ga - go,
                                  levels = design)
        fits <- contrasts.fit(fit, contrast)
        ebFit <- eBayes(fits)

        deg_sig_list <- topTable(ebFit, coef = 1, adjust.method = 'fdr', number = Inf)
        deg_sig_list$branch.contrast <- paste0(bid, "_vs_other")
        deg_sig_list$Gene <- rownames(deg_sig_list)
        total.deg.list <- rbind(total.deg.list, deg_sig_list)
      }
    } else if (is.null(branch.id.2)) {
      if (verbose) message(Sys.time(), " Some of branches will be calculated.")
      pdata$contrast <- "go"
      pdata$contrast[pdata$branch.id %in% branch.id] = "ga"
      design <- model.matrix(~ 0 + as.factor(contrast), data = pdata)
      colnames(design) <- str_replace_all(colnames(design), fixed("as.factor(contrast)"), "")
      fit <- lmFit(t(edata), design)
      contrast <- makeContrasts(ga_go = ga - go,
                                levels = design)
      fits <- contrasts.fit(fit, contrast)
      ebFit <- eBayes(fits)

      deg_sig_list <- topTable(ebFit, coef = 1, adjust.method = 'fdr', number = Inf)
      deg_sig_list$branch.contrast <- paste0(paste0(branch.id, collapse = "-"), "_vs_other")
      deg_sig_list$Gene <- rownames(deg_sig_list)
      total.deg.list <- deg_sig_list
    } else {
      if (verbose) message(Sys.time(), " Some of branches will be calculated.")
      pdata$contrast <- "gz"
      pdata$contrast[pdata$branch.id %in% branch.id] = "ga"
      pdata$contrast[pdata$branch.id %in% branch.id.2] = "go"
      design <- model.matrix(~ 0 + as.factor(contrast), data = pdata)
      colnames(design) <- str_replace_all(colnames(design), fixed("as.factor(contrast)"), "")
      fit <- lmFit(t(edata), design)
      contrast <- makeContrasts(ga_go = ga - go,
                                levels = design)
      fits <- contrasts.fit(fit, contrast)
      ebFit <- eBayes(fits)

      deg_sig_list <- topTable(ebFit, coef = 1, adjust.method = 'fdr', number = Inf)
      deg_sig_list$branch.contrast <- paste0(paste0(branch.id, collapse = "-"), "_vs_", paste0(branch.id.2, collapse = "-"))
      deg_sig_list$Gene <- rownames(deg_sig_list)
      total.deg.list <- deg_sig_list
    }
  }
  if (verbose) message(Sys.time(), " Calculating differentially expressed markers completed")
  return(total.deg.list)
}