R/bdgopt.R

In macs3-project/MACSr: MACS: Model-based Analysis for ChIP-Seq

#' bdgopt
#'
#' Operations on score column of bedGraph file. Note: All regions on
#' the same chromosome in the bedGraph file should be continuous so
#' only bedGraph files from MACS3 are accpetable.
#'
#' @param ifile MACS score in bedGraph. Note: this must be a bedGraph
#'     file covering the ENTIRE genome. REQUIRED
#' @param method Method to modify the score column of bedGraph
#'     file. Available choices are: multiply, add, max, min, or
#'     p2q. 1) multiply, the EXTRAPARAM is required and will be
#'     multiplied to the score column. If you intend to divide the
#'     score column by X, use value of 1/X as EXTRAPARAM. 2) add, the
#'     EXTRAPARAM is required and will be added to the score
#'     column. If you intend to subtract the score column by X, use
#'     value of -X as EXTRAPARAM. 3) max, the EXTRAPARAM is required
#'     and will take the maximum value between score and the
#'     EXTRAPARAM. 4) min, the EXTRAPARAM is required and will take
#'     the minimum value between score and the EXTRAPARAM. 5) p2q,
#'     this will convert p-value scores to q-value scores using
#'     Benjamini-Hochberg process. The EXTRAPARAM is not
#'     required. This method assumes the scores are -log10 p-value
#'     from MACS3. Any other types of score will cause unexpected
#'     errors.", default="p2q"
#' 
#' @param extraparam The extra parameter for METHOD. Check the detail
#'     of -m option.
#' @param outputfile Output filename. Mutually exclusive with
#'     --o-prefix. The number and the order of arguments for --ofile
#'     must be the same as for -m.
#' @param outdir The output directory.
#' @param log Whether to capture logs.
#' @param verbose Set verbose level of runtime message. 0: only show
#'     critical message, 1: show additional warning message, 2: show
#'     process information, 3: show debug messages. DEFAULT:2
#' @return `macsList` object.
#' @export
#' @examples
#' eh <- ExperimentHub::ExperimentHub()
#' CHIP <- eh[["EH4558"]]
#' CTRL <- eh[["EH4563"]]
#' c1 <- callpeak(CHIP, CTRL, gsize = 5.2e7, cutoff_analysis = TRUE,
#'                outdir = tempdir(), name = "callpeak_narrow0",
#'                store_bdg = TRUE)
#' cfile <- grep("treat_pileup.bdg", c1$outputs, value = TRUE)
#' bdgopt(cfile, method = "min", extraparam = 10,
#'        outdir = tempdir(), outputfile = "bdgopt_min.bdg")
bdgopt <- function(ifile,
                   method = c("multiply", "add", "p2q", "max", "min"),
                   extraparam = numeric(),
                   outputfile = character(),
                   outdir = ".", log = TRUE, verbose = 2L){
    method <- match.arg(method)
    ifile <- normalizePath(ifile)
    cl <- basiliskStart(env_macs)
    on.exit(basiliskStop(cl))
    res <- basiliskRun(cl, function(.namespace, outdir){
        opts <- .namespace()$Namespace(ifile = ifile,
                                       method = method,
                                       extraparam = list(extraparam),
                                       ofile = outputfile,
                                       outdir = outdir,
                                       verbose = verbose)
        .bdgopt <- reticulate::import("MACS3.Commands.bdgopt_cmd")
        if(log){
            reticulate::py_capture_output(.bdgopt$run(opts))
        }else{
            .bdgopt$run(opts)
        }
    }, .namespace = .namespace, outdir = outdir)
    if(log){
        message(res)
    }

    ofile <- file.path(outdir, outputfile)
    args <- as.list(match.call())
    macsList(arguments = args, outputs = ofile, log = res)
}