R/benjaminiHochbergCorrection.R

In nicohuttmann/htmnanalysis: Analysis of Omics Data

#' Does Benjamini-Hochberg correction
#'
#' @param data data containing p-value column
#'
#' @return
#' @export
#'
#'
benjaminiHochbergCorrection <- function(data) {

  if(ok("Should p-values be corrected by Benjamini-Hochberg correction? ")) {

    bhc <- getNumber("False discovery rate? (default = 0.05) ")

    data <- data[order(data[,"p.value"]), ]
    data <- cbind(data, `(i/m)Q` = bhc * (1:nrow(data)) / nrow(data), significant = 0)

    for(i in nrow(data):1) {

      if(data[i, "(i/m)Q"] > data[i, "p.value"]) {
        data[1:i, "significant"] <- 1
        break()
      }

    }

  }

  data

}