R/runContrasts.R

In DGEobj.utils: Differential Gene Expression (DGE) Analysis Utility Toolkit

#' Build contrast matrix and calculate contrast fits
#'
#' Takes a DGEobj and a named list of contrasts to build. The DGEobj must
#' contain a limma Fit object and associated designMatrix. Returns the DGEobj with
#' contrast fit(s), contrast matrix, and topTable/topTreat dataframes added.
#'
#' The contrastList is a named list.  The values are composed of column names
#' from the designMatrix of the DGEobj.  Each contrast is named to give it a
#' short, recognizable name to be used for display purposes.
#'
#' Example contrastList \cr
#'
#' contrastList = list( \cr
#'    T1 = "treatment1 - control", \cr
#'    T2 = "treatment2 - control" \cr
#' ) \cr
#'
#' where treatment1, treatment2, and control are column names in the designMatrix.
#'
#' The returned DGEobj list contains the new items:
#' \itemize{
#'     \item{"contrastMatrix"} {a matrix}
#'     \item{"Fit.Contrasts"} {a Fit object}
#'     \item{"topTableList"} {a List of dataframes}
#'     \item{"topTreatList"} {a List of dataframes: if enabled}
#' }
#'
#' @param dgeObj A DGEobj object containing a Fit object and design matrix. (Required)
#' @param designMatrixName The name of the design matrix within dgeObj to use for
#'    contrast analysis. (Required)
#' @param contrastList A named list of contrasts. (Required)
#' @param contrastSetName Name for the set of contrasts specified in
#'   contrastList.  Defaults to "<fitName>_cf". Only needed to create 2 or more
#'   contrast sets from the same initial fit.
#' @param runTopTable Runs topTable on the specified contrasts. (Default = TRUE)
#' @param runTopTreat Runs topTreat on the specified contrasts. (Default = FALSE)
#' @param foldChangeThreshold Only applies to topTreat (Default = 1.5)
#' @param runEBayes Runs eBayes after contrast.fit (Default = TRUE)
#' @param robust eBayes robust option. 'statmod' package must be installed to perform required
#' calculations if robust = TRUE (Default = TRUE)
#' @param proportion Proportion of genes expected to be differentially expressed.
#'   (used by eBayes) (Default = 0.01)
#' @param qValue Set TRUE to include Q-values in topTable output. (Default = FALSE)
#' @param IHW Set TRUE to add FDR values from the IHW package. (Default = FALSE)
#' @param verbose Set TRUE to print some information during processing. (Default = FALSE)
#'
#' @return The DGEobj with contrast matrix, fit and topTable/topTreat dataframes added.
#'
#' @examples
#' \dontrun{
#'   # NOTE: Requires the limma and statmod packages
#'
#'    myDGEobj <- readRDS(system.file("exampleObj.RDS", package = "DGEobj"))
#'
#'    # Name the design matrix to be used (see inventory(myDGEobj))
#'    designMatrixName <- "ReplicateGroupDesign"
#'
#'    # Define the named contrasts from design matrix column names
#'    contrastList  <- list(BDL_v_Sham = "ReplicateGroupBDL - ReplicateGroupSham",
#'                          EXT1024_v_BDL = "ReplicateGroupBDL_EXT.1024  - ReplicateGroupBDL",
#'                          Nint_v_BDL = "ReplicateGroupBDL_Nint - ReplicateGroupBDL",
#'                          Sora_v_BDL = "ReplicateGroupBDL_Sora - ReplicateGroupBDL")
#'
#'
#'    myDGEobj <- runContrasts(myDGEobj,
#'                             designMatrixName=designMatrixName,
#'                             contrastList=contrastList,
#'                             contrastSetName = "SecondContrastSet",
#'                             qValue = TRUE,
#'                             IHW = FALSE,
#'                             runTopTable = TRUE,
#'                             runTopTreat = TRUE,
#'                             foldChangeThreshold = 1.25)
#'    DGEobj::inventory(myDGEobj)
#'}
#'
#' @importFrom DGEobj addItem getItem
#' @importFrom assertthat assert_that
#' @importFrom stringr str_c
#'
#' @export
runContrasts <- function(dgeObj,
                         designMatrixName,
                         contrastList,
                         contrastSetName = NULL,
                         runTopTable = TRUE,
                         runTopTreat = FALSE,
                         foldChangeThreshold = 1.5,
                         runEBayes = TRUE,
                         robust = TRUE,
                         proportion = 0.01,
                         qValue = FALSE,
                         IHW = FALSE,
                         verbose = FALSE) {
    assertthat::assert_that(requireNamespace("limma", quietly = TRUE),
                            msg = "limma package is required to create the contrasts.")

    assertthat::assert_that(!missing(dgeObj),
                            !is.null(dgeObj),
                            "DGEobj" %in% class(dgeObj),
                            msg = "dgeObj must be specified and should be of class 'DGEobj'.")
    assertthat::assert_that(!missing(designMatrixName),
                            !is.null(designMatrixName),
                            is.character(designMatrixName),
                            length(designMatrixName) == 1,
                            designMatrixName %in% names(dgeObj),
                            msg = "designMatrixName must be a signular character value and one of dgeobj names.")
    assertthat::assert_that("list" %in% class(contrastList),
                            !missing(contrastList),
                            !is.null(names(contrastList)),
                            msg = "contrastList must specified and must be a named list.")
    assertthat::assert_that(foldChangeThreshold >= 0,
                            msg = "foldChangeThreshold must be greater than or equal to 0.")
    assertthat::assert_that(!(runTopTable == FALSE & runTopTreat == FALSE),
                            msg = "One of runTopTable or runTopTreat must be TRUE.")
    assertthat::assert_that(designMatrixName %in% names(dgeObj),
                            msg = "The specified designMatrixName not found in dgeObj.")

    fitName <- paste(designMatrixName, "_fit", sep = "")
    assertthat::assert_that(fitName %in% names(dgeObj),
                            msg = "The specified fitName object not found in dgeObj.")

    if (any(is.null(runEBayes),
            !is.logical(runEBayes),
            length(runEBayes) != 1)) {
        warning("runEBayes must be a singular logical value. Assigning default value TRUE")
        runEBayes = TRUE
    }

    if (any(is.null(robust),
            !is.logical(robust),
            length(robust) != 1)) {
        warning("robust must be a singular logical value. Assigning default value TRUE")
        robust = TRUE
    }

    if (robust) {
        assertthat::assert_that(requireNamespace("statmod", quietly = TRUE),
                                msg = "'statmod' package is required to run eBayes calculations")
    }

    if (any(is.null(proportion),
            !is.numeric(proportion),
            length(proportion) != 1)) {
        warning("proportion must be a singular numeric value. Assigning default value 0.01")
        proportion = 0.01
    }

    if (any(is.null(qValue),
            !is.logical(qValue),
            length(qValue) != 1)) {
        warning("qValue must be a singular logical value. Assigning default value FALSE")
        qValue = FALSE
    }

    if (any(is.null(IHW),
            !is.logical(IHW),
            length(IHW) != 1)) {
        warning("IHW must be a singular logical value. Assigning default value FALSE")
        IHW = FALSE
    }

    if (any(is.null(verbose),
            !is.logical(verbose),
            length(verbose) != 1)) {
        warning("verbose must be a singular logical value. Assigning default value FALSE")
        verbose = FALSE
    }

    if (any(is.null(contrastSetName),
            !is.character(contrastSetName),
            length(contrastSetName) != 1)) {
        warning("contrastSetName must be a character value. Assigning default value: ", fitName)
            contrastSetName <- fitName
    }

    assertthat::assert_that(!(paste0(contrastSetName, '_cm') %in% names(dgeObj)),
                            msg = "The contrastSetName already exists in dgeObj.")

    funArgs <- match.call()
    do.call("require", list("limma"))

    # Retrieve designMatrix & fit
    designMatrix <- DGEobj::getItem(dgeObj, designMatrixName)
    fit <- DGEobj::getItem(dgeObj, fitName)

    # Run the contrast fit
    ContrastMatrix <- tryCatch({
        do.call("makeContrasts",
                list(contrasts = contrastList,
                     levels    = designMatrix))
    },
    error = function(e) {
        message("Unexpected error: ", e$message, " happened during limma contrast matrix creation")
        return(NULL)
    })

    MyFit.Contrasts <- tryCatch({
        do.call("contrasts.fit",
                list(fit       = fit,
                     contrasts = ContrastMatrix))
    },
    error = function(e) {
        message("Unexpected error: ", e$message, " happened during limma contrast fitting from linear model")
        return(NULL)
    })

    # Run eBayes
    if (runEBayes) {
        if (verbose) {
            .tsmsg(stringr::str_c("running EBayes: proportion = ", proportion))
        }

        MyFit.Contrasts <- tryCatch({
            do.call("eBayes",
                    list(fit        = MyFit.Contrasts,
                         robust     = robust,
                         proportion = proportion))
        },
        error = function(e) {
            message("Unexpected error: ", e$message, " happened during eBayes statistics computation")
            return(NULL)
        })

        MyFit.Contrasts.treat <- tryCatch({
            do.call("treat",
                    list(fit        = MyFit.Contrasts,
                         lfc        = log2(foldChangeThreshold),
                         robust     = robust))
        },
        error = function(e) {
            message("Unexpected error: ", e$message, " happened during eBayes (treat) statistics computation")
            return(NULL)
        })
    }

    # Run topTable on each contrast and add each DF to a list
    if (runTopTable) {
        if (verbose) {
            .tsmsg("Running topTable...")
        }
        # Run topTable via lapply to generate a bunch of contrasts.
        MyCoef <- as.list(1:length(contrastList))
        topTableList <- lapply(MyCoef, function(x) {
            tryCatch({
                do.call("topTable",
                        list(fit     = MyFit.Contrasts,
                             coef    = x,
                             confint = T,
                             number  = Inf,
                             p.value = 1,
                             sort.by = "none"))
            },
            error = function(e) {
                message("Unexpected error: ", e$message, " happened while extracting top-ranked genes from the linear model fit")
                return(NULL)
            })
        })

        # Transfer the contrast names
        names(topTableList) = names(contrastList)

        if (qValue) {
            topTableList <- runQvalue(topTableList)
        }
        if (IHW) {
            IHW_result <- runIHW(topTableList)
            topTableList <- IHW_result[[1]]
        }
    }

    if (runTopTreat) {
        if (verbose) {
            .tsmsg("Running topTreat...")
        }
        # Run topTreat via lapply to generate a bunch of contrasts.
        LFC    <- log2(foldChangeThreshold)
        MyCoef <- as.list(1:length(contrastList))
        topTreatList = lapply(MyCoef, function(x) {
            tryCatch({
                do.call("topTreat",
                        list(fit     = MyFit.Contrasts.treat,
                             coef    = x,
                             confint = T,
                             lfc     = LFC,
                             number  = Inf,
                             p.value = 1,
                             sort.by = "none"))
            },
            error = function(e) {
                message("Unexpected error: ", e$message, " happened while extracting top-ranked genes from the linear model fit")
                return(NULL)
            })
        })

        # Transfer the contrast names
        names(topTreatList) = names(contrastList)
    }

    # Capture the contrast matrix
    dgeObj <- DGEobj::addItem(dgeObj,
                              item = ContrastMatrix,
                              itemName = paste(contrastSetName, "_cm", sep = ""),
                              itemType = "contrastMatrix",
                              funArgs = funArgs,
                              parent = fitName)

    if (runTopTable) {
        # Add the contrast fit
        dgeObj <- DGEobj::addItem(dgeObj,
                                  item = MyFit.Contrasts,
                                  itemName = paste(contrastSetName, "_cf", sep = ""),
                                  itemType = "contrast_fit",
                                  funArgs = funArgs,
                                  parent = fitName)

        # Add the topTable DFs
        listNames <- names(topTableList)
        for (i in 1:length(topTableList)) {
            dgeObj <- DGEobj::addItem(dgeObj,
                                      item = topTableList[[i]],
                                      itemName = listNames[[i]],
                                      itemType = "topTable",
                                      funArgs = funArgs,
                                      parent = paste(contrastSetName, "_cf", sep = ""))
        }
    }

    if (runTopTreat) {
        dgeObj <- DGEobj::addItem(dgeObj,
                                  item = MyFit.Contrasts.treat,
                                  itemName = paste(contrastSetName, "_cft", sep = ""),
                                  itemType = "contrast_fit_treat",
                                  funArgs = funArgs,
                                  parent = fitName)

        listNames <- names(topTreatList)
        for (i in 1:length(topTreatList)) {
            dgeObj <- DGEobj::addItem(dgeObj,
                                      item = topTreatList[[i]],
                                      itemName = paste(listNames[i], "_treat", sep = ""),
                                      itemType = "topTreat",
                                      funArgs = funArgs,
                                      parent = paste(contrastSetName, "_cft", sep = ""))
        }
    }
    dgeObj
}