R/runCusanovich2018.R

In plbaldoni/scChIPseqsim: Single-cell ChIP-seq Data Simulator Based on Pseudo-bulk Profiles

#' Run Cusanovich2018
#'
#' Add description here
#'
#' @param sc SingleCellExperiment data
#' @param peaks either TRUE (call peaks), FALSE (do not call peaks), or
#' a vector of logicals of length equal to nrow(sc) indicating the sequence of bins pertaining to peak regions
#' @param returnSc whether or not to return the utilized SingleCellExperiment data and the vector of peaks
#' @param clustering either 'all' (the default), 'none', 'hierarchical', 'louvain', or 'kmeans'
#' @param seed the seed
#' @param legend.title title of the legend
#' @param shape shape of the umap plot
#' @param shape.title title of the shape legend
#' @param title title of the plot
#'
#' @return Add return here
#'
#' @details Add details here
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#'
#' @export
#'
runCusanovich2018 = function(sc,peaks = TRUE,returnSc = FALSE,clustering = 'all',title = 'Cusanovich2018',
                             seed = 2020,legend.title = 'Cell Type',shape = NULL,shape.title = NULL){

    cond = NULL

    # Calling peaks

    if(!is.list(sc)){
        sc <- callPeaks(sc = sc,peaks = peaks)
    }

    # Creating return variable

    ret <- list()
    ret[['true']] <- SummarizedExperiment::colData(sc$sc)$Cluster
    ret[['peaks']] <- sc$peaks

    # Pre processing

    object <- SummarizedExperiment::assay(sc$sc,'counts')

    # Starting the clock

    start_time = Sys.time()

    # Running tryCatch

    tryCatch(
        {
            retCatch <- list()

            # Creating Cusanovich2018 object

            num_pcs <- 20

            object <- filterZeros(object)
            object <- 1*(object>0)

            # Building the model

            SVD = irlba::irlba(Matrix::Diagonal(x=log(1+ncol(object)/ Matrix::rowSums(object))) %*% Matrix::t(Matrix::t(object)/Matrix::colSums(object)),
                               nv=num_pcs,nu=num_pcs)

            # Dimension reduction

            feature <- as.matrix(Matrix::Diagonal(x=SVD$d) %*% Matrix::t(SVD$v))

            ret[['feature']] <- feature

            # Running clustering

            if(!clustering == 'none'){
                ret[['clustering']] <- runClustering(sc = sc$sc,feature = feature,clustering = clustering, method = 'Cusanovich2018')
            }

            # Stopping the clock

            ret[['time']] <- difftime(Sys.time(),start_time,units = 'min')[[1]]

            # Plots

            ret[['umap']] <- plot_umap(x = run_umap(feature), labels = ret[['true']],
                                       title = title,legend.title = legend.title,
                                       shape = shape, shape.title = shape.title)

        },
        error = function(cond){assign('ret',errorCatch(cond,true = SummarizedExperiment::colData(sc$sc)$Cluster,peaks = sc$peaks,method = 'Cusanovich2018'),inherits = TRUE)}
    )

    if(returnSc){
        ret[['sc']] <- sc
    }

    return(ret)
}