R/nodeResult.R

In fionarhuang/treeclimbR: An algorithm to find optimal signal levels in a tree

#' Extract table with node-level DA/DS results
#'
#' Extract a table with the top-ranked nodes from a DA/DS analysis output
#' (generated by \code{\link{runDA}} or \code{\link{runDS}}).
#'
#' @author Ruizhu Huang, Charlotte Soneson
#' @export
#'
#' @param object The output from \code{\link{runDA}} or \code{\link{runDS}}.
#' @param type Either "DA" (for \strong{object} from \code{\link{runDA}}) or
#'     "DS" (for \strong{object} from \code{\link{runDS}}).
#' @param n An integer indicating the maximum number of entities to return.
#' @param adjust_method A character string specifying the method used to adjust
#'     p-values for multiple testing. See \code{\link[stats]{p.adjust}} for
#'     possible values.
#' @param sort_by A character string specifying the sorting method. This will
#'     be passed to \code{\link[edgeR]{topTags}}.
#'     Possibilities are "PValue" for p-value, "logFC" for absolute log-fold
#'     change or "none" for no sorting.
#' @param p_value A numeric cutoff value for adjusted p-values. This will
#'     be passed to \code{\link[edgeR]{topTags}}. Only entities
#'     with adjusted p-values equal or lower than specified are returned.
#'
#' @returns A data frame with results for all nodes passing the imposed
#'     thresholds. The columns \strong{logFC}, \strong{logCPM},
#'     \strong{PValue}, \strong{FDR}, \strong{F} (or \strong{LR}) are from (the
#'     output table of) \code{\link[edgeR]{topTags}}. The \strong{node} column
#'     stores the node number for each entity. Note: \strong{FDR} is corrected
#'     over all features and nodes when the specified \code{type = "DS"}.
#'
#' @importFrom edgeR topTags
#' @importFrom dplyr arrange slice filter bind_rows desc
#' @importFrom TreeSummarizedExperiment convertNode
#' @importFrom stats p.adjust
#' @importFrom rlang .data
#'
#' @examples
#' suppressPackageStartupMessages({
#'     library(TreeSummarizedExperiment)
#' })
#'
#' lse <- readRDS(system.file("extdata", "da_sim_100_30_18de.rds",
#'                            package = "treeclimbR"))
#' tse <- aggTSE(lse, rowLevel = showNode(tree = rowTree(lse),
#'                                        only.leaf = FALSE))
#' dd <- model.matrix( ~ group, data = colData(tse))
#' out <- runDA(TSE = tse, feature_on_row = TRUE,
#'              assay = "counts", option = "glmQL",
#'              design = dd, contrast = NULL,
#'              normalize = TRUE)
#'
#' ## Top 10 nodes with DA
#' nodeResult(out, n = 10)
#'
nodeResult <- function(object, n = 10, type = c("DA", "DS"),
                       adjust_method = "BH", sort_by = "PValue",
                       p_value = 1) {

    ## Check arguments
    ## -------------------------------------------------------------------------
    .assertVector(x = object, type = "list")
    .assertScalar(x = n, type = "numeric")
    .assertScalar(x = adjust_method, type = "character")
    .assertScalar(x = sort_by, type = "character",
                  validValues = c("PValue", "logFC", "none"))
    .assertScalar(x = p_value, type = "numeric")

    ## Generate topTags table
    ## -------------------------------------------------------------------------
    type <- match.arg(type)
    if (type == "DA") {
        tt <- edgeR::topTags(object = object$edgeR_results, n = n,
                             adjust.method = adjust_method,
                             sort.by = sort_by,
                             p.value = p_value)$table
        ## Add node column
        nod <- TreeSummarizedExperiment::convertNode(
            tree = object$tree, node = rownames(tt))
        ct <- cbind(node = nod, tt)
    } else if (type == "DS") {
        res <- object$edgeR_results
        tt <- lapply(seq_along(res), FUN = function(x) {
            xx <- edgeR::topTags(object = res[[x]], n = Inf,
                                 adjust.method = adjust_method,
                                 sort.by = sort_by,
                                 p.value = 1)$table
            ## Add node column
            nod <- TreeSummarizedExperiment::convertNode(
                tree = object$tree, node = names(res)[x])
            cx <- cbind(xx, node = nod, feature = rownames(xx))
            return(cx)
        })
        ct <- do.call(dplyr::bind_rows, tt)

        ## Estimate FDR across features
        ct$FDR <- stats::p.adjust(p = ct$PValue, method = adjust_method)
        n <- min(nrow(ct), n)
        ## Sort
        if (sort_by == "PValue") {
            ct <- ct |>
                dplyr::arrange(.data$PValue) |>
                dplyr::filter(.data$FDR <= p_value) |>
                dplyr::slice(seq_len(n))
        } else if (sort_by == "logFC") {
            ct <- ct |>
                dplyr::arrange(dplyr::desc(abs(.data$logFC))) |>
                dplyr::filter(.data$FDR <= p_value) |>
                dplyr::slice(seq_len(n))
        } else {
            ## No sorting
            ct <- ct |>
                dplyr::filter(.data$FDR <= p_value) |>
                dplyr::slice(seq_len(n))
        }
    }

    return(ct)
}