R/runScasat.R

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

#' Run Scasat
#'
#' 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
#'
runScasat = function(sc,peaks = TRUE,returnSc = FALSE,clustering = 'all',seed = 2020,title = 'Scasat',
                     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 Scasat object

            object <- filterZeros(object)
            object <- Scasat_getJaccardDist(object)

            # Building the model

            # Use the similar number of dimensions as shown in tutorial https://github.com/ManchesterBioinference/Scasat/blob/master/ScAsAT_functions_Buenrostro_All_Bam_Together.ipynb
            fit <- stats::cmdscale(as.dist(object),eig=TRUE, k=15)

            # Dimension reduction

            feature <- as.matrix(t(fit$points))

            ret[['feature']] <- feature

            # Running clustering

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

            # 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 = 'Scasat'),inherits = TRUE)}
    )

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

    return(ret)
}