@@ -12,7 +12,6 @@ Imports:

     plyr,

     dplyr,

     foreach,

-    gtools,

     ggplot2,

     plotly,

     RColorBrewer,

@@ -1,6 +1,6 @@

 MIT License

-Copyright (c) 2017 Joshua D Campbell

+Copyright (c) 2018 Joshua D Campbell

 Permission is hereby granted, free of charge, to any person obtaining a copy

 of this software and associated documentation files (the "Software"), to deal

@@ -10,27 +10,41 @@ available_models = c("celda_C", "celda_G", "celda_CG")

 #' @param counts A count matrix

 #' @param model Which celda sub-model to run. Options include "celda_C" (cell clustering), "celda_G" (gene clustering), "celda_CG" (gene and cell clustering)

 #' @param sample.label A numeric vector, character vector, or factor indicating the originating sample for each cell (column) in the count matrix. By default, every cell will be assumed to be from an independent sample.

-#' @param nchains The number of chains of Gibbs sampling to run for every combination of K/L parameters. Defaults to 1

+#' @param K Number of desired cell subpopulation clusters. Required for celda_C and celda_CG models.

+#' @param L Number of desired gene clusters. Required for celda_G and celda_CG models.

+#' @param alpha Non-zero concentration parameter for sample Dirichlet distribution (celda_C / celda_CG only)

+#' @param beta Non-zero concentration parameter for gene Dirichlet distribution

+#' @param gamma The Dirichlet distribution parameter for Psi; adds a pseudocount to each gene within each transcriptional state (celda_G / celda_CG only)

+#' @param delta The Dirichlet distribution parameter for Eta; adds a gene pseudocount to the numbers of genes each state (celda_G / celda_CG only)

+#' @param nchains The number of chains of Gibbs sampling to run for every combination of K/L parameters. Defaults to 1 (celda_G / celda_CG only)

 #' @param bestChainsOnly Return only the best chain (by final log-likelihood) per K/L combination.

 #' @param cores The number of cores to use for parallell Gibbs sampling. Defaults to 1.

 #' @param seed The base seed for random number generation

 #' @param verbose Print log messages during celda chain execution

 #' @param logfile_prefix Prefix for log files from worker threads and main process. 

-#' @param ... Model specific parameters

 #' @return Object of class "celda_list", which contains results for all model parameter combinations and summaries of the run parameters

 #' @import foreach

 #' @export

-celda = function(counts, model, sample.label=NULL, nchains=1, bestChainsOnly=TRUE, cores=1, seed=12345, verbose=FALSE, logfile_prefix="Celda", ...) {

+celda = function(counts, model, sample.label=NULL, K=NULL, L=NULL, alpha=1, beta=1, 

+                 delta=1, gamma=1, max.iter=200, z.init=NULL, y.init=NULL,

+                 stop.iter=10, split.on.iter=10, process.counts=FALSE, 

+                 nchains=1, bestChainsOnly=TRUE, cores=1, 

+                 seed=12345, verbose=FALSE, logfile_prefix="Celda") {

+ 

   message(paste(Sys.time(), "Starting celda."))

-  validateArgs(counts, model, sample.label, nchains, cores, seed, ...)

+  params.list = buildParamList(counts, model, sample.label, K, L, alpha, beta, delta,

+                               gamma, max.iter, z.init, y.init, stop.iter, split.on.iter,

+                               process.counts, nchains, cores, seed)

+  

   # Redirect stderr from the worker threads if user asks for verbose

   logfile = paste0(logfile_prefix, "_main_log.txt")

+  params.list$logfile = logfile

   cl = if (verbose) parallel::makeCluster(cores, outfile=logfile) else parallel::makeCluster(cores)

   doParallel::registerDoParallel(cl)

   # Details for each model parameter / chain combination 

-  runs = expand.grid(chain=1:nchains, ...)

+  runs = expand.grid(plyr::compact(list(chain=1:nchains, K=K, L=L)))

   runs$index = as.numeric(rownames(runs))

   if (verbose) print(runs)

@@ -42,20 +56,22 @@ celda = function(counts, model, sample.label=NULL, nchains=1, bestChainsOnly=TRU

   # count matrix was used.

   counts = processCounts(counts)

   count.checksum = digest::digest(counts, algo="md5")

+  params.list$count.checksum = count.checksum

+   

+   res.list = foreach(i = 1:nrow(runs), .export=model, .combine = c, .multicombine=TRUE) %dopar% {

+    chain.params = params.list

+    chain.params$seed = all.seeds[ifelse(i %% nchains == 0, nchains, i %% nchains)]

-  

-  res.list = foreach(i = 1:nrow(runs), .export=model, .combine = c, .multicombine=TRUE) %dopar% {

-    chain.seed = all.seeds[ifelse(i %% nchains == 0, nchains, i %% nchains)]

-    

-    if (verbose) {

+    if (isTRUE(verbose)) {

       ## Generate a unique log file name based on given prefix and parameters

-      logfile = paste0(logfile_prefix, "_", paste(paste(colnames(runs), runs[i,], sep="-"), collapse="_"), "_Seed-", chain.seed, "_log.txt")

-      res = do.call(model, c(list(counts=counts, sample.label=sample.label, count.checksum=count.checksum, seed=chain.seed, logfile=logfile, process.counts=F), c(runs[i,-1])))

+      chain.params$logfile = paste0(logfile_prefix, "_",  paste(paste(colnames(runs), runs[i,], sep="-"), collapse="_"),  "_Seed-", chain.params$seed, "_log.txt")

+      res = do.call(model, chain.params)

     } else {

-      res = suppressMessages(do.call(model, c(list(counts=counts, sample.label=sample.label, count.checksum=count.checksum, seed=chain.seed, logfile=NULL, process.counts=F), c(runs[i,-1]))))

+      chain.params$logfile = NULL

+      res = suppressMessages(do.call(model, chain.params))

     }

     return(list(res))

-  }  

+  }

   parallel::stopCluster(cl)

   celda.res = list(run.params=runs, res.list=res.list, 

                    content.type=model, count.checksum=count.checksum)

@@ -80,10 +96,43 @@ celda = function(counts, model, sample.label=NULL, nchains=1, bestChainsOnly=TRU

 }

+# Build a list of parameters tailored to the specific celda model being run,

+# validating the provided parameters along the way

+buildParamList = function(counts, model, sample.label, K, L, alpha, beta, delta,

+                          gamma, max.iter, z.init, y.init, stop.iter, split.on.iter,

+                          process.counts, nchains, cores, seed) {

+  

+  validateArgs(counts, model, sample.label, nchains, cores, seed, K=K, L=L)

+  

+  params.list = list(counts=counts,

+                     sample.label=sample.label,

+                     max.iter=max.iter,

+                     stop.iter=stop.iter,

+                     split.on.iter=split.on.iter,

+                     process.counts=process.counts)

+  

+  if (model %in% c("celda_C", "celda_CG")) {

+    params.list$alpha = alpha

+    params.list$beta = beta

+    params.list$K = K

+    z.init=z.init

+  } 

+  if (model %in% c("celda_G", "celda_CG")) {

+    params.list$beta = beta

+    params.list$delta = delta

+    params.list$gamma = gamma

+    params.list$L = L

+    y.init=y.init

+  }

+  

+  return(params.list)

+}

+

+

 # Sanity check arguments to celda() to ensure a smooth run.

 # See parameter descriptions from celda() documentation.

 validateArgs = function(counts, model, sample.label, 

-                         nchains, cores, seed, K=NULL, L=NULL, ...) {

+                         nchains, cores, seed, K=NULL, L=NULL) { #, ...) {

   model_args = names(formals(model))

   if ("K" %in% model_args && is.null(K)) {

     stop("Must provide a K parameter when running a celda_C or celda_CG model")

@@ -84,15 +84,13 @@ simulateCells.celda_C = function(model, S=10, C.Range=c(10, 100), N.Range=c(100,

 #' @param z.init Initial values of z. If NULL, z will be randomly sampled. Default NULL.

 #' @param process.counts Whether to cast the counts matrix to integer and round(). Defaults to TRUE.

 #' @param logfile If NULL, messages will be displayed as normal. If set to a file name, messages will be redirected messages to the file. Default NULL.

-#' @param ... additonal parameters

 #' @return An object of class celda_C with clustering results and Gibbs sampling statistics

 #' @export

-celda_C = function(counts, sample.label=NULL, K,

-					alpha=1, beta=1, 

-                   	stop.iter = 10, split.on.iter=10, max.iter=200, 

-                   	count.checksum=NULL, seed=12345,

-                   	z.init = NULL, process.counts=TRUE, logfile=NULL, ...) {

-

+celda_C = function(counts, sample.label=NULL, K, alpha=1, beta=1, 

+                 	 stop.iter = 10, split.on.iter=10, max.iter=200, 

+                 	 count.checksum=NULL, seed=12345,

+                 	 z.init = c(), process.counts=TRUE, logfile=NULL) {

+  

   ## Error checking and variable processing

   if (isTRUE(process.counts)) {

     counts = processCounts(counts)  

@@ -130,6 +128,7 @@ celda_C = function(counts, sample.label=NULL, K,

   iter = 1L

   num.iter.without.improvement = 0L

   do.cell.split = TRUE

+  logMessages(date(), "Max iter:", max.iter, logfile=logfile, append=TRUE)

   while(iter <= max.iter & num.iter.without.improvement <= stop.iter) {

     next.z = cC.calcGibbsProbZ(counts=counts, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.by.C=n.by.C, n.CP=n.CP, z=z, s=s, K=K, nG=nG, nM=nM, alpha=alpha, beta=beta)

@@ -323,13 +323,12 @@ simulateCells.celda_CG = function(model, S=10, C.Range=c(50,100), N.Range=c(500,

 #' @param y.init Initial values of y. If NULL, y will be randomly sampled. Default NULL.

 #' @param logfile The name of the logfile to redirect messages to.

 #' @param count.checksum An MD5 checksum for the provided counts matrix

-#' @param ... Additional parameters

 #' @export

 celda_CG = function(counts, sample.label=NULL, K, L,

-					alpha=1, beta=1, delta=1, gamma=1, 

+                    alpha=1, beta=1, delta=1, gamma=1, 

                     max.iter=200, stop.iter = 10, split.on.iter=10,

                     seed=12345, count.checksum=NULL,

-			        z.init = NULL, y.init = NULL, logfile=NULL, ...) {

+                    z.init = NULL, y.init = NULL, logfile=NULL) {

   ## Error checking and variable processing

   if (isTRUE(processCounts)) {

@@ -198,14 +198,12 @@ cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y,

 #' @param y.init Initial values of y. If NULL, y will be randomly sampled. Default NULL.

 #' @param process.counts Whether to cast the counts matrix to integer and round(). Defaults to TRUE.

 #' @param logfile The name of the logfile to redirect messages to.

-#' @param ...  Additional parameters

 #' @keywords LDA gene clustering gibbs

 #' @export

-celda_G = function(counts, L,

-					beta=1, delta=1, gamma=1,

+celda_G = function(counts, L, beta=1, delta=1, gamma=1,

 					stop.iter=10, split.on.iter=10, max.iter=200,

-                    count.checksum=NULL, seed=12345, 

-                    y.init=NULL, process.counts=TRUE, logfile=NULL, ...) {

+          count.checksum=NULL, seed=12345, 

+          y.init=NULL, process.counts=TRUE, logfile=NULL) {

   ## Error checking and variable processing

   if (isTRUE(process.counts)) {

@@ -4,9 +4,11 @@

 \alias{celda}

 \title{Run the celda Bayesian hierarchical model on a matrix of counts.}

 \usage{

-celda(counts, model, sample.label = NULL, nchains = 1,

-  bestChainsOnly = TRUE, cores = 1, seed = 12345, verbose = FALSE,

-  logfile_prefix = "Celda", ...)

+celda(counts, model, sample.label = NULL, K = NULL, L = NULL, alpha = 1,

+  beta = 1, delta = 1, gamma = 1, max.iter = 200, z.init = NULL,

+  y.init = NULL, stop.iter = 10, split.on.iter = 10,

+  process.counts = FALSE, nchains = 1, bestChainsOnly = TRUE, cores = 1,

+  seed = 12345, verbose = FALSE, logfile_prefix = "Celda")

 }

 \arguments{

 \item{counts}{A count matrix}

@@ -15,7 +17,19 @@ celda(counts, model, sample.label = NULL, nchains = 1,

 \item{sample.label}{A numeric vector, character vector, or factor indicating the originating sample for each cell (column) in the count matrix. By default, every cell will be assumed to be from an independent sample.}

-\item{nchains}{The number of chains of Gibbs sampling to run for every combination of K/L parameters. Defaults to 1}

+\item{K}{Number of desired cell subpopulation clusters. Required for celda_C and celda_CG models.}

+

+\item{L}{Number of desired gene clusters. Required for celda_G and celda_CG models.}

+

+\item{alpha}{Non-zero concentration parameter for sample Dirichlet distribution (celda_C / celda_CG only)}

+

+\item{beta}{Non-zero concentration parameter for gene Dirichlet distribution}

+

+\item{delta}{The Dirichlet distribution parameter for Eta; adds a gene pseudocount to the numbers of genes each state (celda_G / celda_CG only)}

+

+\item{gamma}{The Dirichlet distribution parameter for Psi; adds a pseudocount to each gene within each transcriptional state (celda_G / celda_CG only)}

+

+\item{nchains}{The number of chains of Gibbs sampling to run for every combination of K/L parameters. Defaults to 1 (celda_G / celda_CG only)}

 \item{bestChainsOnly}{Return only the best chain (by final log-likelihood) per K/L combination.}

@@ -26,8 +40,6 @@ celda(counts, model, sample.label = NULL, nchains = 1,

 \item{verbose}{Print log messages during celda chain execution}

 \item{logfile_prefix}{Prefix for log files from worker threads and main process.}

-

-\item{...}{Model specific parameters}

 }

 \value{

 Object of class "celda_list", which contains results for all model parameter combinations and summaries of the run parameters

@@ -6,8 +6,8 @@

 \usage{

 celda_C(counts, sample.label = NULL, K, alpha = 1, beta = 1,

   stop.iter = 10, split.on.iter = 10, max.iter = 200,

-  count.checksum = NULL, seed = 12345, z.init = NULL,

-  process.counts = TRUE, logfile = NULL, ...)

+  count.checksum = NULL, seed = 12345, z.init = c(),

+  process.counts = TRUE, logfile = NULL)

 }

 \arguments{

 \item{counts}{A numeric count matrix}

@@ -35,8 +35,6 @@ celda_C(counts, sample.label = NULL, K, alpha = 1, beta = 1,

 \item{process.counts}{Whether to cast the counts matrix to integer and round(). Defaults to TRUE.}

 \item{logfile}{If NULL, messages will be displayed as normal. If set to a file name, messages will be redirected messages to the file. Default NULL.}

-

-\item{...}{additonal parameters}

 }

 \value{

 An object of class celda_C with clustering results and Gibbs sampling statistics

@@ -7,7 +7,7 @@

 celda_CG(counts, sample.label = NULL, K, L, alpha = 1, beta = 1,

   delta = 1, gamma = 1, max.iter = 200, stop.iter = 10,

   split.on.iter = 10, seed = 12345, count.checksum = NULL,

-  z.init = NULL, y.init = NULL, logfile = NULL, ...)

+  z.init = NULL, y.init = NULL, logfile = NULL)

 }

 \arguments{

 \item{counts}{A numeric count matrix.}

@@ -41,8 +41,6 @@ celda_CG(counts, sample.label = NULL, K, L, alpha = 1, beta = 1,

 \item{y.init}{Initial values of y. If NULL, y will be randomly sampled. Default NULL.}

 \item{logfile}{The name of the logfile to redirect messages to.}

-

-\item{...}{Additional parameters}

 }

 \description{

 celda Cell and Gene Clustering Model

@@ -6,8 +6,7 @@

 \usage{

 celda_G(counts, L, beta = 1, delta = 1, gamma = 1, stop.iter = 10,

   split.on.iter = 10, max.iter = 200, count.checksum = NULL,

-  seed = 12345, y.init = NULL, process.counts = TRUE, logfile = NULL,

-  ...)

+  seed = 12345, y.init = NULL, process.counts = TRUE, logfile = NULL)

 }

 \arguments{

 \item{counts}{A numeric count matrix}

@@ -35,8 +34,6 @@ celda_G(counts, L, beta = 1, delta = 1, gamma = 1, stop.iter = 10,

 \item{process.counts}{Whether to cast the counts matrix to integer and round(). Defaults to TRUE.}

 \item{logfile}{The name of the logfile to redirect messages to.}

-

-\item{...}{Additional parameters}

 }

 \description{

 Provides cluster assignments for all genes in a provided single-cell