R/initializer.R

In plbaldoni/mikado: What the package does (short line)

#' Initializer for mikado
#'
#' Add description here
#'
#' @param object a mikadoDataSet
#' @param control list of control arguments from controlEM()
#'
#' @return
#' Add return here
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#' @references
#' \url{https://github.com/plbaldoni/epigraHMM}
#'
#' @import data.table
#' @importFrom  stats rbinom dbinom glm quasibinomial
#' @rawNamespace import(SummarizedExperiment, except = shift)
#' @rawNamespace import(S4Vectors, except = c(first,second))
#' @importFrom foreach %dopar%
#' @importFrom magrittr %<>%
#'
#' @export
initializer <- function(object,control){
    
    cores = chr = init =  windowSize = j = x = NULL
    
    # Creating control elements
    
    for(i in seq_along(control)){assign(names(control)[i],control[[i]])}
    
    # Registering number of cores
    
    doParallel::registerDoParallel(cores = cores)
    
    # Running the initializers in parallel
    
    if(init == 'HMM'){
        result <- foreach::foreach(i = iterators::iter(seq_len(ncol(object)))) %dopar% {
            result <- tryCatch({S4Vectors::metadata(HMM(object = object[,i],control = control))},
                               error = function(err){list('pi' = rep(NA,2),'gamma' = matrix(NA,2,2),'psi' = c('(Intercept)' = NA,'(Intercept)' = NA),
                                                          'peak' = rep(NA,nrow(object)))})
            return(result)
        }
        # Calculating the euclidean distance on t(sqrt(2*prob)) is equal to the Hellinger distance on t(prob)
        
        peak <- sqrt(2*do.call(cbind,lapply(result,function(x){x[['peak']]})))
    }
    if(init == 'fast'){
        rowm <- do.call(rbind,lapply(seq_len(ncol(object)),function(x){
            data.table(j = x,
                       i = which((scale(data.table::frollmean(assay(object)[,x],n = windowSize, algo = 'fast',fill = 0, align="center"),center = T,scale = F)>0)),
                       x = 1)
        }))
        # Calculating the euclidean distance on t(sqrt(2*prob)) is equal to the Hellinger distance on t(prob)
        
        peak <- sqrt(2*with(rowm,Matrix::sparseMatrix(i = i,j = j,x = x,
                                                    dims = list(nrow(object),ncol(object)),
                                                    dimnames = list(NULL,NULL))))
        result <- lapply(seq_len(ncol(object)), function(x){list('pi' = c(0.999,0.001),
                                                                 'gamma' = matrix(c(0.99,0.01,0.01,0.99),2,2),
                                                                 'psi' = c('(Intercept)' = -4.6,'(Intercept)' = -2.3))})
    }
    
    # Calculating pairwise distances & dendogram
    dist <- stats::as.dist(proxyC::dist(peak,method = 'euclidean',margin = 2))
    # dist[dist==0] <- NA
    clust <- fastcluster::hclust(dist,method = 'ward.D2')
    
    # Returning elements
    
    return(list('dist' = dist,
                'hclust' = clust,
                'par' = lapply(result,function(x){x[c('pi','gamma','psi')]})))
}