R/DoWilcoxonTest.R

In YC3/mirNet: microRNA netowrk analyses

#' @title Apply two-sample Wilcoxon tests (‘Mann-Whitney’ test) for each row of input data matrix.
#'
#' @description \code{DoWilcoxonTest} applies two-sample Wilcoxon test for each row of input data matrix. P-values are adjusted using Benjamini–Hochberg procedure.
#'
#' @param data A matrix, the normalized gene/microRNA expression dataset, should be a numeric matrix, with rows referring to genes/microRNAs and columns to samples.
#' @param pheno A vector of sample phenotypes. Sample phenotype in a scientific research could be treatment/control, normal/cancer or smoker/non-smoker. Different phenotypes should each be encoded as 0/1 when inputting to \code{DoWilcoxonTest}, for example, Normal-0; Cancer-1.
#'
#' @return A table with rows for all genes/miRNAs (ranked by significance) and columns of log2 fold change, p-value, adjusted p-value (default to Benjamini–Hochberg procedure).
#'
#' @details This function performs two-sample Wilcoxon tests. It is recommended to use a non-parametric test when data lack normality and t-tools (\code{\link[mirNet]{DoTest}}) can be very misleading.
#'
#' @seealso \code{\link[mirNet]{pairComp}} and \code{\link[limma]{DoTest}}.
#'
#' @importFrom stats wilcox.test
#' @importFrom stats p.adjust
#'
#' @export DoWilcoxonTest
#'
#' @examples
#' # prepare your normalized data matrix
#' data.m <- matrix(rnorm(120), nrow = 20, ncol = 6)
#'
#' # prepare the phenotype info (0-control; 1-treatment)
#' class.v <- c(0, 0, 0, 1, 1, 1)
#'
#' # run function
#' DoWilcoxonTest(data = data.m, pheno = class.v)


DoWilcoxonTest <- function(data, pheno){

    tmpfun <- function(data.v, pheno){

        w.o <- wilcox.test(data.v[pheno == 1], data.v[pheno == 0])

        ave.c <- mean(data.v[pheno == 1])
        ave.n <- mean(data.v[pheno == 0])

        list(logFC = ave.c - ave.n, p = w.o$p.value)
    }

    res.l <- do.call(rbind, apply(data, 1, tmpfun, pheno))

    res <- do.call(rbind, apply(res.l, 1, as.data.frame, stringsAsFactors = FALSE))
    res.adj <- as.matrix(cbind(res, fdr = p.adjust(res[, 'p'], method = 'fdr')))
    res.sort <- res.adj[sort(abs(unlist(res.adj[, 'logFC'])), index.return = T, decreasing = TRUE)$ix, ]
    
    res.sort
}