R/DGEAR.R

In DGEAR: Differential Gene Expression Analysis with R

#' Differential Gene Expression Analysis with R
#'
#' @description
#' Main orchestration function that runs five statistical tests (Welch t-test,
#' one-way ANOVA, Dunnett's test, Half's modified t-test, and Wilcoxon-Mann-Whitney
#' U-test) on a gene expression matrix, combines their BH-adjusted p-values with
#' Fisher's combined probability method, and identifies differentially expressed
#' genes (DEGs) by majority voting.
#'
#' @param dataframe A numeric matrix or data.frame of gene expression values
#'   (rows = genes, columns = samples).  Raw intensity / count values are
#'   automatically log2-transformed when they appear to be on a linear 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 votting_cutoff Integer. Minimum number of tests (out of 5) that must
#'   independently declare a gene significant for it to be included in the
#'   majority-vote DEG list (default \code{3}).  Must be between 1 and 5.
#' @param annot_df Optional annotation data.frame with columns \code{ID} and
#'   \code{Gene.Symbol} (or \code{Gene.symbol}) that maps probe/row identifiers
#'   to gene symbols.  When \code{NULL} (default) row names of \code{dataframe}
#'   are used directly as gene identifiers.  Typical source: a GEO SOFT family
#'   file parsed with \code{read.delim()}.
#'
#' @details
#' The function internally calls:
#' \itemize{
#'   \item \code{\link{perform_t_test}} — Welch two-sample t-test
#'   \item \code{\link{perform_anova}} — one-way ANOVA
#'   \item \code{\link{perform_dunnett_test}} — Dunnett's test
#'   \item \code{\link{perform_h_test}} — Half's modified t-test
#'   \item \code{\link{perform_wilcox_test}} — Wilcoxon rank-sum test
#' }
#' Each test independently assigns an FDR flag (1 = significant, 0 = not).
#' The five flags are summed per gene; genes whose sum meets or exceeds
#' \code{votting_cutoff} are reported as DEGs (majority voting, Boyer & Moore
#' 1991).  Combined p-values across all five tests are computed with Fisher's
#' method via \code{\link[metapod]{parallelFisher}}.
#'
#' Annotation via \code{annot_df} is entirely optional.  When supplied, the
#' first gene symbol listed for each probe (delimited by \code{ /// }) is used.
#' When absent, row names serve as identifiers, making the function fully
#' self-contained without GEO annotation files.
#'
#' @return A named list with four elements:
#'   \describe{
#'     \item{\code{DEGs}}{Data.frame of gene identifiers that passed majority
#'       voting.}
#'     \item{\code{FDR_Table}}{Wide data.frame with BH-adjusted p-values from
#'       every test, the Fisher-combined FDR, the ensemble voting score, and
#'       log2 fold change for every gene.}
#'     \item{\code{Results_Table}}{Concise data.frame with \code{G_Symbol},
#'       \code{CombineFDR}, \code{log2FC}, and \code{Ensemble} score.}
#'     \item{\code{IndividualTests}}{Named list of the raw output from each
#'       of the five test functions (each containing a \code{Table} and a
#'       \code{DEGs} element).}
#'   }
#'
#' @references
#' Boyer, R.S. and Moore, J.S. (1991). MJRTY — A Fast Majority Vote Algorithm.
#' In \emph{Automated Reasoning: Essays in Honor of Woody Bledsoe}, pp. 105–117.
#' Springer, Dordrecht. \doi{10.1007/978-94-011-3488-0_5}
#'
#' @export
#' @importFrom stats quantile
#' @importFrom metapod parallelFisher
#' @examples
#' library(DGEAR)
#' data("gene_exp_data")
#'
#' ## Basic usage — no annotation file needed
#' result <- DGEAR(dataframe    = gene_exp_data,
#'                 con1         = 1,
#'                 con2         = 10,
#'                 exp1         = 11,
#'                 exp2         = 20,
#'                 alpha        = 0.05,
#'                 votting_cutoff = 2)
#' result$DEGs
#' head(result$FDR_Table)
#'
#' ## With an optional annotation data.frame (GEO SOFT format)
#' ## annot <- read.delim("GSExxxxx_family.soft")
#' ## result <- DGEAR(dataframe = gene_exp_data,
#' ##                 con1 = 1, con2 = 10, exp1 = 11, exp2 = 20,
#' ##                 annot_df = annot)
DGEAR <- function(dataframe, con1, con2, exp1, exp2,
                  alpha = 0.05, votting_cutoff = 3,
                  annot_df = NULL) {
  
  # ---- Validate inputs -------------------------------------------------------
  if (!is.data.frame(dataframe) && !is.matrix(dataframe))
    stop("'dataframe' must be a matrix or data.frame.")
  if (con1 < 1 || con2 > ncol(dataframe) || con1 > con2)
    stop("Invalid control column range (con1:con2).")
  if (exp1 < 1 || exp2 > ncol(dataframe) || exp1 > exp2)
    stop("Invalid experiment column range (exp1:exp2).")
  if (alpha <= 0 || alpha >= 1)
    stop("'alpha' must be strictly between 0 and 1.")
  if (!votting_cutoff %in% 1:5)
    stop("'votting_cutoff' must be an integer between 1 and 5.")
  
  # ---- Run the five statistical tests ----------------------------------------
  message("Running Welch t-test ...")
  t_res  <- perform_t_test(dataframe, con1, con2, exp1, exp2, alpha, annot_df)
  
  message("Running one-way ANOVA ...")
  o_res  <- perform_anova(dataframe, con1, con2, exp1, exp2, alpha, annot_df)
  
  message("Running Dunnett's test ...")
  d_res  <- perform_dunnett_test(dataframe, con1, con2, exp1, exp2, alpha, annot_df)
  
  message("Running Half's t-test ...")
  h_res  <- perform_h_test(dataframe, con1, con2, exp1, exp2, alpha, annot_df)
  
  message("Running Wilcoxon rank-sum test ...")
  u_res  <- perform_wilcox_test(dataframe, con1, con2, exp1, exp2, alpha, annot_df)
  
  # ---- Retrieve shared quantities from the first result ----------------------
  G_Symbol <- t_res$Table$G_Symbol
  log2FC   <- t_res$Table$log2FC
  
  # ---- Ensemble / majority voting --------------------------------------------
  fdr_mat <- data.frame(
    T.FDR = t_res$Table$fdr,
    O.FDR = o_res$Table$fdr,
    D.FDR = d_res$Table$fdr,
    H.FDR = h_res$Table$fdr,
    W.FDR = u_res$Table$fdr
  )
  fdr_mat[is.na(fdr_mat)] <- 0L
  ensemble_score <- rowSums(fdr_mat)
  
  DEGs <- data.frame(
    DEGs = G_Symbol[ensemble_score >= votting_cutoff],
    stringsAsFactors = FALSE
  )
  
  # ---- Combined p-value (Fisher's method via metapod) -----------------------
  P_list <- list(
    P1 = t_res$Table$BH,
    P2 = o_res$Table$BH,
    P3 = d_res$Table$BH,
    P4 = h_res$Table$BH,
    P5 = u_res$Table$BH
  )
  # Replace NA BH values with 1 before combining (conservative)
  P_list <- lapply(P_list, function(p) { p[is.na(p)] <- 1; p })
  CombineFDR <- metapod::parallelFisher(P_list)$p.value
  
  # ---- Summary tables --------------------------------------------------------
  FDR_Table <- data.frame(
    G_Symbol   = G_Symbol,
    log2FC     = log2FC,
    T.FDR      = t_res$Table$BH,
    O.FDR      = o_res$Table$BH,
    D.FDR      = d_res$Table$BH,
    H.FDR      = h_res$Table$BH,
    W.FDR      = u_res$Table$BH,
    CombineFDR = CombineFDR,
    Ensemble   = ensemble_score,
    stringsAsFactors = FALSE
  )
  
  Results_Table <- data.frame(
    G_Symbol   = G_Symbol,
    CombineFDR = CombineFDR,
    log2FC     = log2FC,
    Ensemble   = ensemble_score,
    stringsAsFactors = FALSE
  )
  
  message(sprintf("Done. %d DEGs identified at votting_cutoff = %d.",
                  nrow(DEGs), votting_cutoff))
  
  list(
    DEGs          = DEGs,
    FDR_Table     = FDR_Table,
    Results_Table = Results_Table,
    IndividualTests = list(
      t_test   = t_res,
      anova    = o_res,
      dunnett  = d_res,
      half_t   = h_res,
      wilcoxon = u_res
    )
  )
}