R/findCutpoint.R

In genomaths/MethylIT.utils: MethylIT utility

#' @rdname findCutpoint
#' @title Find a cutoff of divergences of methylation level values
#' @description A function to help on the decision of which is the best cutoff
#'     value for DIMP/DMP predictions. The genome-wide methylation changes that
#'     occurs in any living organism is the result of the superposition of
#'     several stochastic processes: the inherent stochasticity of biological
#'     processes and, particular, ultimately, it derives from the stochasticity
#'     of biochemical processes. On this scenario, there is not way to say with
#'     absolute determinism where a given value of an information divergence is
#'     a true positive value or a true negative value. All what we can do is
#'     the estimation of performance indicators like accuracy, sensitivity,
#'     false positive rate, etc., to evaluate the consequences of our decision
#'     on what we consider a true positive or a true negative. For example, a
#'     difference of methylation levels of 100% observed in a treatment group of
#'     samples in given cytosine position does not means that this difference
#'     will not be observed in some sample from the control group. Without any
#'     doubt about it, such a different can be found in control samples as well.
#'     The fluctuation theorem guaranty such an outcome, which in the current
#'     context is a consequence of the action of second law of thermodynamics on
#'     living organisms.
#' @details Given a numerical vector of cutoff values for the divergences of
#'     methylation level values, or p-values cutoffs, this function search for
#'     the cutoff value that yield the best classification performance for the
#'     specified performance indicator.
#' @param LR A list of GRanges, a GRangesList, a CompressedGRangesList object.
#'     Each GRanges object from the list must have at least two columns: a
#'     column containing the total variation of methylation level (TV,
#'     difference of methylation levels) and a column containing a divergence of
#'     methylation levels (it could be TV or  Hellinger divergence) or a
#'     column with a p-value from where the cutpoint will be found
#'     (see example).
#'
#' @param min.tv Minimum value for the total variation distance (TVD; absolute
#'     value of methylation levels differences, \eqn{TVD = abs(TV)}). Only
#'     sites/ranges k with \eqn{TVD_{k} > min.tv} are analyzed. Defaul
#'     min.tv = 0.25.
#' @param tv.cut A cutoff for the total variation distance to be applied to each
#'     site/range. Sites/ranges k with \eqn{TVD_{k} < tv.cut} are considered
#'     TRUE negatives and sites with \eqn{TVD_{k} > tv.cut} TRUE positives.
#'     Its value must be a number \eqn{0 < tv.cut < 1}. A possible value
#'     for tv.cut would be, e.g., the minimum value of *TV* found in the
#'     treatment group after the potential DMPs are estimated. Default is
#'     \eqn{tv.cut = 0.5}.
#' @param predcuts A numerical vector of possible cutoff values (cutpoints) for
#'     a divergence of methylation levels value or a p-value, according with the
#'     magnitude given in div.col or in pval.col, respectively. For each
#'     cutpoint k the values greater than predcuts[k] are predicted TRUE
#'     (positives), otherwise are predicted FALSE (negatives).
#' @param tv.col Column number where the total variation is located in the
#'     metadata from each GRanges object.
#' @param div.col Column number for divergence of methylation levels used in the
#'     estimation of the cutpoints. Default: NULL. One of the parameter values
#'     div.col or pval.col must be given.
#' @param pval.col Column number for p-value used in the estimation of the
#'     cutpoints. Default: NULL. One of the parameter values div.col or pval.col
#'     must be given.
#' @param stat An integer number indicating the statistic to be used in the
#'     testing. The mapping for statistic names are:
#'     0 = "All" 1 = "Accuracy", 2 = "Sensitivity", 3 = "Specificity",
#'     4 = "Pos Pred Value", 5 = "Neg Pred Value", 6 = "Precision",
#'     7 = "Recall", 8 = "F1",  9 = "Prevalence", 10 = "Detection Rate",
#'     11 = "Detection Prevalence", 12 = "Balanced Accuracy".
#' @param maximize Whether to maximize the performance indicator given in
#'     parameter 'stat'. Default: TRUE.
#' @param num.cores,tasks Paramaters for parallele computation using package
#'     \code{\link[BiocParallel]{BiocParallel-package}}: the number of cores to
#'     use, i.e. at most how many child processes will be run simultaneously
#'     (see \code{\link[BiocParallel]{bplapply}} and the number of tasks per job
#'     (only for Linux OS).
#'
#' @importFrom BiocParallel MulticoreParam bplapply SnowParam
#' @importFrom S4Vectors mcols
#' @importFrom caret confusionMatrix
#' @return A list with the classification repformance results for the best
#'     cutoff value in the ranges of predcuts supplied.
#' @export
#' @author Robersy Sanchez
#' @examples
#' # load simulated data of potential methylated signal
#' data(sim_ps)
#'
#' # Vector of cutoff values
#' cuts = c(2, 5, 10, 15, 18, 20, 21, 22, 25, 27, 30, 35, 40,
#'         45, 50, 55, 60)
#
#' # === To find the cutpoint that maximize the accuracy ===
#' pre.cut.acc = findCutpoint(LR = PS, min.tv = 0.25, tv.cut = 0.5,
#'                             predcuts = cuts, tv.col = 7L, div.col = 9,
#'                             stat = 1, num.cores = 15)

findCutpoint <- function(LR, min.tv=0.25, tv.cut=0.5, predcuts, tv.col,
                       div.col=NULL, pval.col=NULL, stat = 1, maximize = TRUE,
                       num.cores = 1L, tasks=tasks) {
   # statistic names
   # 0 = "All" 1 = "Accuracy", 2 = "Sensitivity", 3 = "Specificity",
   # 4 = "Pos Pred Value", 5 = "Neg Pred Value", 6 = "Precision", 7 = "Recall",
   # 8 = "F1",  9 = "Prevalence", 10 = "Detection Rate",
   # 11 = "Detection Prevalence", 12 = "Balanced Accuracy"

   if (Sys.info()['sysname'] == "Linux") {
       bpparam <- MulticoreParam(workers=num.cores, tasks=tasks)
   } else {
       bpparam <- SnowParam(workers = num.cores, type = "SOCK")
   }

   res <- bplapply(predcuts, function(k, LR)
               classPerform(LR=LR, min.tv=min.tv, tv.cut=tv.cut, cutoff=k,
                           tv.col=tv.col, div.col=div.col, pval.col=pval.col,
                           stat=stat), LR, BPPARAM = bpparam)
   res <- unlist(res)
   if (!all(is.na(res))) {
       if (maximize) cutp = predcuts[which.max(res)]
       else  cutp = predcuts[which.min(res)]
       perf <- classPerform(LR=LR, min.tv=min.tv, tv.cut=tv.cut, cutoff=cutp,
                           tv.col=tv.col, div.col=div.col, pval.col=pval.col,
                           stat=0)

       res = list(optimCutpoint=cutp, Statistic=res, Performance=perf[[1]],
                   FDR=perf[[2]])
   } else res = list(optimCutpoint=NA, Statistic=NA, Performance=NA,
                   FDR=NA)
   return(res)
}