R/T-Test.R

In DGEAR: Differential Gene Expression Analysis with R

#' Welch Two-Sample t-Test for Differential Gene Expression
#'
#' Performs an independent two-sample Welch t-test for every gene (row) in the
#' expression matrix, applies Benjamini-Hochberg (BH) correction, and returns
#' a results table together with the list of significant DEGs.
#'
#' @param dataframe A numeric matrix or data.frame of gene expression values
#'   (rows = genes, columns = samples).  Values are automatically log2-
#'   transformed when they appear to be on a linear / intensity scale.
#' @param con1 Integer. Index of the first control column.
#' @param con2 Integer. Index of the last control column.
#' @param exp1 Integer. Index of the first experiment column.
#' @param exp2 Integer. Index of the last experiment column.
#' @param alpha Numeric significance threshold for BH-adjusted p-values
#'   (default \code{0.05}).
#' @param annot_df Optional annotation data.frame with columns \code{ID} and
#'   \code{Gene.Symbol} (or \code{Gene.symbol}).  When \code{NULL} (default)
#'   row names of \code{dataframe} are used as gene identifiers.
#'
#' @return A named list:
#'   \describe{
#'     \item{Table}{Data.frame with columns \code{G_Symbol}, \code{log2FC},
#'       \code{statistic.t}, \code{p.value}, \code{BH}, \code{fdr}.}
#'     \item{DEGs}{Data.frame of gene identifiers whose BH-adjusted p-value
#'       is \eqn{\le} \code{alpha}.}
#'   }
#'
#' @export
#' @importFrom stats t.test p.adjust na.omit quantile
#' @examples
#' library(DGEAR)
#' data("gene_exp_data")
#' result <- perform_t_test(dataframe = gene_exp_data,
#'                          con1 = 1, con2 = 10,
#'                          exp1 = 11, exp2 = 20)
#' head(result$Table)
#' result$DEGs
perform_t_test <- function(dataframe, con1, con2, exp1, exp2,
                           alpha = 0.05, annot_df = NULL) {
  prep       <- .prepare_data(dataframe, con1, con2, exp1, exp2)
  G_Symbol   <- .resolve_symbols(prep$dataframe, annot_df)
  con        <- prep$con
  exp        <- prep$exp
  log2FC     <- prep$log2FC
  
  t_results <- t(sapply(seq_len(nrow(con)), function(i) {
    res <- t.test(as.numeric(con[i, ]), as.numeric(exp[i, ]))
    c(statistic.t = res$statistic, p.value = res$p.value)
  }))
  
  t_stat         <- data.frame(G_Symbol, log2FC, t_results,
                               stringsAsFactors = FALSE)
  t_stat$BH      <- p.adjust(t_stat$p.value, method = "BH")
  t_stat$fdr     <- as.integer(t_stat$BH <= alpha)
  
  DEGs <- data.frame(DEGs = na.omit(t_stat$G_Symbol[t_stat$fdr == 1L]),
                     stringsAsFactors = FALSE)
  list(Table = t_stat, DEGs = DEGs)
}