Create numeric contrast for DESeq2::results() function

deseq2_numericContrast.R

#' Create numeric contrast for DESeq2::results() function
#' 
#' @param object a DESeqDataSet object.
#' @param group1 logical expression indicating the columns from colData(object) that form part of group1 (numerator).
#' @param group2 logical expression indicating the columns from colData(object) that form part of group2 (denominator).
#' 
#' @examples
#' # example data for a two-factor design with interaction
#' dds <- makeExampleDESeqDataSet()
#' dds$timepoint <- factor(rep(1:3, 4))
#' design(dds) <- ~ condition*timepoint
#' 
#' # compare A vs B for timepoint 2
#' numericContrast(dds,
#'                 group1 = condition == "A" & timepoint == "2",
#'                 group2 = condition == "B" & timepoint == "2")
#' 
#' # compare timepoint 1 vs 2 for condition A
#' numericContrast(dds,
#'                 group1 = condition == "A" & timepoint == "1",
#'                 group2 = condition == "A" & timepoint == "2")
#' 
#' # compare condition A vs B across all timepoints
#' numericContrast(dds,
#'                 group1 = condition == "A",
#'                 group2 = condition == "B")
numericContrast <- function(object, 
                            group1,
                            group2){
  # get logical vector for each group
  group1_condition <- eval(substitute(group1), colData(object), parent.frame())
  group2_condition <- eval(substitute(group2), colData(object), parent.frame())
  
  # extract the model matrix from the object
  modmat <- model.matrix(design(object), colData(object))
  
  # is this approach fine when there are different number of replicates in groups? - I don't think so, see function below instead
  group1_coef <- colMeans(modmat[group1_condition, , drop=FALSE])
  group2_coef <- colMeans(modmat[group2_condition, , drop=FALSE])
  
  # return numeric contrast
  return(group1_coef - group2_coef)
}

#' Create a numeric vector of coefficients from a DESeq model
#' 
#' @param object a DESeqDataSet object.
#' @param subset logical expression indicating the columns from colData(object) that form part of the group.
#' 
#' @examples
#' example data for a two-factor design with interaction
#' dds <- makeExampleDESeqDataSet()
#' dds$timepoint <- factor(rep(1:3, 4))
#' design(dds) <- ~ condition*timepoint
#' 
#' #' compare A vs B for timepoint 2
#' group1 <- getNumericCoef(dds, condition == "A" & timepoint == "2")
#' group2 <- getNumericCoef(dds, condition == "B" & timepoint == "2")
#' results(dds, contrast = group1 - group2)
#' 
#' #' compare timepoint 1 vs 2 for condition A
#' group1 <- getNumericCoef(dds, condition == "A" & timepoint == "1")
#' group2 <- getNumericCoef(dds, condition == "A" & timepoint == "2")
#' results(dds, contrast = group1 - group2)
#' 
#' #' compare condition A vs B across all timepoints
#' group1 <- getNumericCoef(dds, condition == "A")
#' group2 <- getNumericCoef(dds, condition == "B")
#' results(dds, contrast = group1 - group2)
getNumericCoef <- function(object, subset, weighted = FALSE){
  # logical vector of rows
  rows <- eval(substitute(subset), colData(object), parent.frame())
  stopifnot(is.logical(rows))
  
  # get model matrix
  modmat <- model.matrix(design(object), colData(object))
  stopifnot(!is.null(modmat))
  
  # subset matrix
  modmat <- modmat[rows, , drop=FALSE]
  
  # calculate vector
  if(!weighted){
    modmat <- modmat[!duplicated(modmat), , drop=FALSE]
  } 
  
  return(colMeans(modmat))
  
}