R/cutoff.R

In plbaldoni/mixHMM: Differential Peak Calling via Hidden Markov Model With Mixture of Negative Binomial Distributions

#' Determine posterior probability cutoff given an FDR threshold
#'
#' For a given FDR threshold, `cutoff()` outputs a cutoff for the posterior probabilities of windows belonging to a differential state.
#'
#' @param postprob numeric vector of posterior probabilities of windows belonging to differential state
#' @param alpha FDR threshold (default is 0.05)
#'
#' @return Cutoff for posterior probabilities belonging to a differential state
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#' @references \url{https://biostats.bepress.com/jhubiostat/paper115/}
#' @references \url{https://doi.org/10.1093/biostatistics/5.2.155}
#'
#' @examples
#' data(H3K36me3)
#' ChIP = SummarizedExperiment::assay(H3K36me3)
#' offset = ZIMHMM::createOffset(ChIP,method = 'loess',span = 1)
#' group = c(1,1,1,2,2,2,3,3,3)
#' B = 2^(length(unique(group)))-2
#' \dontrun{output = mixHMMConstr(ChIP,Control=NULL,offset,group,B,control = controlPeaks())}
#' \dontrun{cutoff(postprob = output$Prob$PostProb2,alpha = 0.05)}
#'
#' @import data.table
#'
#' @export
cutoff <- function(postprob,alpha = 0.05){
    if(any(postprob<0) | any(postprob>1)){stop('Posterior probabilities outside 0-1')}

    DT <- data.table::data.table(notpp = 1-postprob,Window = 1:length(postprob),key = 'Window')
    return(DT[order(notpp),][,FDR := ((cumsum(notpp)/1:nrow(DT))<alpha)][order(Window),]$FDR)

    # if(sum((1-postprob)*(postprob>0))/sum(postprob>0)<alpha){
    #     return(0)
    # } else{
    #     root <- uniroot(f = function(gamma,alpha,postprob){sum((1-postprob)*(postprob>gamma))/sum(postprob>gamma)-alpha},
    #                     interval = c(0,mean(tail(sort(unique(postprob)),2))),alpha = alpha,postprob = postprob)
    #     return(root$root)
    # }
}