R/runCisTopic.R

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

#' Run cisTopic
#'
#' 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
#'
runCisTopic = function(sc,peaks = TRUE,returnSc = FALSE,clustering = 'all',title = 'cisTopic',
                       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')
    rownames(object) <- with(data.table::as.data.table(SummarizedExperiment::rowRanges(sc$sc)),paste0(seqnames,":",start,"-",end))

    # Starting the clock

    start_time = Sys.time()

    # Running tryCatch

    tryCatch(
        {
            retCatch <- list()

            # Creating cisTopic object

            object <- cisTopic::createcisTopicObject(object)
            object <- cisTopic::addCellMetadata(object, cell.data = colData(sc$sc))

            # Building the model

            object <- cisTopic::runWarpLDAModels(object, topic=c(3,6,9),seed=seed,addModels=FALSE)

            # Selecting the model

            object <- cisTopic::selectModel(object,type='derivative')

            # Dimension reduction

            feature <- cisTopic::modelMatSelection(object, 'cell', 'Probability')

            ret[['feature']] <- feature

            # Running clustering

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

            # 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)
}