R/runMikado.R

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

#' Multi-scenario simulator of single-cell epigenomic data
#'
#' Add description here
#'
#' @param path path where the analysis will be saved
#' @param pathdata path where the datasets will be saved
#' @param nsim number of simulated datasets
#'
#' @return Add return here
#'
#' @details Add details here
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#'
#' @rawNamespace import(GenomicRanges, except = c(shift,update))
#' @rawNamespace import(BiocGenerics, except = c(density,mad,weights,residuals,IQR))
#' @import S4Vectors
#'
#' @export
#'
runMikado <- function(path = 'Analysis',pathdata = 'Data',nsim = 100){

    system(paste('mkdir',path))
    path <- normalizePath(file.path(path))
    pathdata <- normalizePath(file.path(pathdata))

    workdir <- getwd()

    # Simulating data

    for(sim in seq_len(nsim)){

        datafile <- file.path(pathdata,paste0('Simulation',sim,'.RData'))
        outfile <- file.path(path,paste0('Simulation',sim,'.RData'))

        sink(file.path(path,paste0('Simulation',sim,'.R')))
        cat('
library(scChIPseqsim)
library(data.table)

load("',datafile,'")

# Running with epigraHMM

Cusanovich2018.epigrahmm <- runCusanovich2018(sc, peaks = TRUE, returnSc = TRUE)
cisTopic.epigrahmm<- runCisTopic(sc = Cusanovich2018.epigrahmm$sc)
SnapATAC.epigrahmm<- runSnapATAC(sc = Cusanovich2018.epigrahmm$sc)
Scasat.epigrahmm<- runScasat(sc = Cusanovich2018.epigrahmm$sc)

dt.epigrahmm <- list()
dt.epigrahmm[["clustering"]] <- rbindlist(list(
            Cusanovich2018.epigrahmm$clustering$metrics,
            cisTopic.epigrahmm$clustering$metrics,
            SnapATAC.epigrahmm$clustering$metrics,
            Scasat.epigrahmm$clustering$metrics))[,Peak := "Pooled"]
dt.epigrahmm[["time"]] <- rbindlist(list(
            data.table(Method = "Cusanovich2018", Time = Cusanovich2018.epigrahmm$time),
            data.table(Method = "cisTopic", Time = cisTopic.epigrahmm$time),
            data.table(Method = "SnapATAC", Time = SnapATAC.epigrahmm$time),
            data.table(Method = "Scasat", Time = Scasat.epigrahmm$time)))[,Peak := "Pooled"]
dt.epigrahmm[["umap"]] <- list(
            "cisTopic" = Cusanovich2018.epigrahmm$umap,
            "Cusanovich2018" = cisTopic.epigrahmm$umap,
            "SnapATAC" = SnapATAC.epigrahmm$umap,
            "Scasat" = Scasat.epigrahmm$umap)

# Running with mikado

Cusanovich2018.mikado <- runCusanovich2018(sc, peaks = "mikado", returnSc = TRUE)
cisTopic.mikado <- runCisTopic(sc = Cusanovich2018.mikado$sc)
SnapATAC.mikado <- runSnapATAC(sc = Cusanovich2018.mikado$sc)
Scasat.mikado <- runScasat(sc = Cusanovich2018.mikado$sc)

dt.mikado <- list()
dt.mikado[["clustering"]] <- rbindlist(list(
            Cusanovich2018.mikado$clustering$metrics,
            cisTopic.mikado$clustering$metrics,
            SnapATAC.mikado$clustering$metrics,
            Scasat.mikado$clustering$metrics))[,Peak := "Differential"]
dt.mikado[["time"]] <- rbindlist(list(
            data.table(Method = "Cusanovich2018", Time = Cusanovich2018.mikado$time),
            data.table(Method = "cisTopic", Time = cisTopic.mikado$time),
            data.table(Method = "SnapATAC", Time = SnapATAC.mikado$time),
            data.table(Method = "Scasat", Time = Scasat.mikado$time)))[,Peak := "Differential"]
dt.mikado[["umap"]] <- list(
            "cisTopic" = Cusanovich2018.mikado$umap,
            "Cusanovich2018" = cisTopic.mikado$umap,
            "SnapATAC" = SnapATAC.mikado$umap,
            "Scasat" = Scasat.mikado$umap)

message("Saving output")

sc.epigrahmm <- Cusanovich2018.epigrahmm$sc
sc.mikado <- Cusanovich2018.mikado$sc

save(sc,sc.epigrahmm,sc.mikado,dt.epigrahmm,dt.mikado,file = "',outfile,'",compress = "xz")
',sep = '')
        sink()

        # Changing directories and submitting job
        setwd(path)
        system(paste0('sbatch -t 4:00:00 R CMD BATCH ',paste0('./Simulation',sim,'.R')))
        setwd(workdir)
    }
}