R/runSnapATAC.R

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

#' Run SnapATAC
#'
#' Add description here
#'
#' @param sc either a SingleCellExperiment data or a list comprised of elements 'sc' (SingleCellExperiment data) and 'peaks' (called peaks logicals)
#' @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
#'
runSnapATAC = function(sc,peaks = TRUE,returnSc = FALSE,clustering = 'all',title = 'SnapATAC',
                       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(
        {
            # Creating SnapATAC object

            object <- SnapATAC::createSnapFromBmat(Matrix::t(object),
                                                   barcodes = colnames(object),
                                                   bins = SummarizedExperiment::rowRanges(sc$sc))

            # Building the model

            object = SnapATAC::makeBinary(object, mat="bmat")

            object = SnapATAC::runDiffusionMaps(obj = object,input.mat = "bmat")

            # Dimension reduction

            feature <- t(object@smat@dmat)

            ret[['feature']] <- feature

            # Running clustering

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

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

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

    return(ret)
}