R/DoTest.R

In YC3/mirNet: microRNA netowrk analyses

#' @title Apply two-sample t-test for each row of input data matrix.
#'
#' @description \code{DoTest} applies two-sample t-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{DoTest}, for example, Normal-0; Cancer-1.
#'
#' @return A table with rows for all genes (ranked by significance) and columns of t-statistic, p-value, adjusted p-value (default to Benjamini–Hochberg procedure).
#'
#' @details This function computes the t-statistic from two-sample comparisons which assumes normality in data. If sample size is too small, it is better to perform a moderated t-test using empirical Bayes method (\code{\link[limma]{eBayes}}).
#'
#' @seealso \code{\link[mirNet]{pairComp}} for applying a moderated t-test for input data matrix, \code{\link[limma]{eBayes}} for Bayes method.
#'
#' @importFrom stats t.test
#' @importFrom stats p.adjust
#'
#' @export DoTest
#'
#' @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
#' DoTest(data = data.m, pheno = class.v)


DoTest <- function(data, pheno){

    tmpfun <- function(data.v, pheno){
        t.o <- t.test(data.v[pheno == 1], data.v[pheno == 0])
        list(t_stat = t.o$stat, p = t.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[, 't_stat'])), index.return = T, decreasing = TRUE)$ix, ]
    
    res.sort
}