@@ -1,14 +1,16 @@

 Package: celda

 Title: CEllular Latent Dirichlet Allocation

-Version: 0.99.8

+Version: 0.99.9

 Authors@R: c(person("Joshua", "Campbell", email = "camp@bu.edu", role = c("aut", "cre")),

              person("Sean", "Corbett", email = "scorbett@bu.edu", role = c("aut")),

-	     person("Yusuke", "Koga", email="ykoga07@bu.edu", role = c("aut")))

+	     person("Yusuke", "Koga", email="ykoga07@bu.edu", role = c("aut")),

+	     person("Zhe", "Wang", email="zhe@bu.edu", role = c("aut")))

 Description: celda leverages Bayesian hierarchical modeling to cluster genes,

     cells, or both simultaneously from single cell sequencing data.

 Depends:

     R (>= 3.5),

     SummarizedExperiment

+VignetteBuilder: knitr

 Imports:

     plyr,

     foreach,

@@ -31,23 +33,29 @@ Imports:

     GenomicRanges,

     data.table,

     Rcpp,

-    RcppEigen

+    RcppEigen,

+    umap,

+    enrichR,

+    vcd,

+    ggpubr,

+    stringi,

+    SummarizedExperiment

 Suggests:

     testthat,

     knitr,

     roxygen2,

     rmarkdown,

-    vcd,

     corrplot,

     Matrix,

     biomaRt,

     Rtsne,

-    covr

+    covr,

+    M3DExampleData,

+    BiocManager

 LinkingTo: Rcpp, RcppEigen

-VignetteBuilder: knitr

-License: MIT

+License: MIT + file LICENSE

 Encoding: UTF-8

 LazyData: true

-RoxygenNote: 6.1.0

-BugReports: https://github.com/definitelysean/celda/issues

+RoxygenNote: 6.1.1

+BugReports: https://github.com/campbio/celda/issues

 biocViews: SingleCell, GeneExpression, Clustering, Sequencing, Bayesian

@@ -10,6 +10,7 @@ export(celdaHeatmap)

 export(celdaPerplexity)

 export(celdaProbabilityMap)

 export(celdaTsne)

+export(celdaUmap)

 export(celda_C)

 export(celda_CG)

 export(celda_G)

@@ -21,6 +22,7 @@ export(distinct_colors)

 export(factorizeMatrix)

 export(featureModuleLookup)

 export(featureModuleTable)

+export(geneSetEnrich)

 export(logLikelihood)

 export(logLikelihood.celda_C)

 export(logLikelihood.celda_CG)

@@ -59,6 +61,7 @@ export(violinPlot)

 exportMethods(celdaHeatmap)

 exportMethods(celdaProbabilityMap)

 exportMethods(celdaTsne)

+exportMethods(celdaUmap)

 exportMethods(clusterProbability)

 exportMethods(factorizeMatrix)

 exportMethods(featureModuleLookup)

@@ -70,6 +73,7 @@ import(grDevices)

 import(graphics)

 import(grid)

 import(gtable)

+import(matrixStats)

 import(plyr)

 import(scales)

 useDynLib(celda,"_colSumByGroup")

@@ -75,12 +75,12 @@ moduleHeatmap <- function(counts, celda.mod, feature.module = 1, top.cells = 100

   gene_ix = match(rownames(filtered_norm.counts), celda.mod@names$row)

   cell_ix = match(colnames(filtered_norm.counts), celda.mod@names$column)

   z.to.plot = c()

-  if(methods::.hasSlot(celda.mod, "z")){

+  if (is(celda.mod, "celda_C") || is(celda.mod, "celda_CG")){

     cell <- distinct_colors(length(unique(celda.mod@clusters$z)))[sort(unique(celda.mod@clusters$z[cell_ix]))]

     names(cell) <- sort(unique(celda.mod@clusters$z[cell_ix]))

     anno_cell_colors <- list(cell = cell)

     z.to.plot = celda.mod@clusters$z[cell.indices]

-  }else{

+  } else {

     anno_cell_colors <- NULL

   }

   plotHeatmap(

@@ -300,7 +300,14 @@ setGeneric("celdaProbabilityMap",

 #' Embeds cells in two dimensions using tSNE based on celda_CG results.

 #' 

 #' @param counts Integer matrix. Rows represent features and columns represent cells. This matrix should be the same as the one used to generate `celda.mod`.

-#' @param celda.mod Celda model. Options available in `celda::available.models`.

+#' @param celda.mod Celda object of class `celda_CG`. 

+#' @param max.cells Integer. Maximum number of cells to plot. Cells will be randomly subsampled if ncol(counts) > max.cells. Larger numbers of cells requires more memory. Default 25000.

+#' @param min.cluster.size Integer. Do not subsample cell clusters below this threshold. Default 100. 

+#' @param modules Integer vector. Determines which features modules to use for tSNE. If NULL, all modules will be used. Default NULL.

+#' @param perplexity Numeric. Perplexity parameter for tSNE. Default 20.

+#' @param max.iter Integer. Maximum number of iterations in tSNE generation. Default 2500.

+#' @param seed Integer. Passed to `set.seed()`. Default 12345. If NULL, no calls to `set.seed()` are made.

+#' @param ... Additional parameters.

 #' @param ... Additional parameters.

 #' @return Numeric Matrix of dimension `ncol(counts)` x 2, colums representing the "X" and "Y" coordinates in the data's t-SNE represetation.

 #' @examples 

@@ -317,6 +324,30 @@ setGeneric("celdaTsne",

            })

+#' Embeds cells in two dimensions using umap.

+#' 

+#' @param counts Integer matrix. Rows represent features and columns represent cells. This matrix should be the same as the one used to generate `celda.mod`.

+#' @param celda.mod Celda object of class `celda_CG`. 

+#' @param max.cells Integer. Maximum number of cells to plot. Cells will be randomly subsampled if ncol(counts) > max.cells. Larger numbers of cells requires more memory. Default 25000.

+#' @param min.cluster.size Integer. Do not subsample cell clusters below this threshold. Default 100. 

+#' @param modules Integer vector. Determines which features modules to use for tSNE. If NULL, all modules will be used. Default NULL.

+#' @param perplexity Numeric. Perplexity parameter for tSNE. Default 20.

+#' @param max.iter Integer. Maximum number of iterations in tSNE generation. Default 2500.

+#' @param seed Integer. Passed to `set.seed()`. Default 12345. If NULL, no calls to `set.seed()` are made.

+#' @param ... Additional parameters.

+#' @param ... Additional parameters.

+#' @return Numeric Matrix of dimension `ncol(counts)` x 2, colums representing the "X" and "Y" coordinates in the data's t-SNE represetation.

+#' @examples 

+#' tsne.res = celdaTsne(celda.CG.sim$counts, celda.CG.mod)

+#' @export

+setGeneric("celdaUmap",

+           signature = "celda.mod",

+           function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                   modules=NULL, umap.config=umap::umap.defaults) {

+             standardGeneric("celdaUmap")

+           })

+

+

 #' Obtain the gene module of a gene of interest

 #' 

 #' This function will output the corresponding feature module for a specified list of genes from a celda model.

@@ -526,6 +526,7 @@ setMethod("clusterProbability",

 #' @seealso `celda_C()` for clustering cells

 #' @examples

 #' perplexity = perplexity(celda.C.sim$counts, celda.C.mod)

+#' @import matrixStats

 #' @export

 setMethod("perplexity",

           signature(celda.mod = "celda_C"),

@@ -617,55 +618,94 @@ setMethod("celdaTsne",

                      initial.dims=20, modules=NULL, perplexity=20, max.iter=2500, 

                      seed=12345, ...) {

-            counts = processCounts(counts)

-            compareCountMatrix(counts, celda.mod)

-            

-            ## Checking if max.cells and min.cluster.size will work

-            if((max.cells < ncol(counts)) & (max.cells / min.cluster.size < celda.mod@params$K)) {

-              stop(paste0("Cannot distribute ", max.cells, " cells among ", 

-                          celda.mod@params$K, " clusters while maintaining a minumum of ", 

-                          min.cluster.size, 

-                          " cells per cluster. Try increasing 'max.cells' or decreasing 'min.cluster.size'."))

-            }

+

+            prepared.count.info = prepareCountsForDimReduction.celda_C(counts, celda.mod, max.cells,

+                                                                       min.cluster.size, modules)

-            ## Select a subset of cells to sample if greater than 'max.cells'

-            total.cells.to.remove = ncol(counts) - max.cells

-            z.include = rep(TRUE, ncol(counts))

-            if(total.cells.to.remove > 0) {

-          	z.ta = tabulate(celda.mod@clusters$z, celda.mod@params$K)

-          	

-          	## Number of cells that can be sampled from each cluster without 

-          	## going below the minimum threshold

-          	cluster.cells.to.sample = z.ta - min.cluster.size          

-          	cluster.cells.to.sample[cluster.cells.to.sample < 0] = 0

-          	

-          	## Number of cells to sample after exluding smaller clusters

-          	## Rounding can cause number to be off by a few, so ceiling is 

-          	## used with a second round of subtraction

-          	cluster.n.to.sample = ceiling((cluster.cells.to.sample / sum(cluster.cells.to.sample)) * total.cells.to.remove)

-          	diff = sum(cluster.n.to.sample) - total.cells.to.remove 

-          	cluster.n.to.sample[which.max(cluster.n.to.sample)] = cluster.n.to.sample[which.max(cluster.n.to.sample)] - diff

-          

-          	## Perform sampling for each cluster

-          	for(i in which(cluster.n.to.sample > 0)) {

-          	  z.include[sample(which(celda.mod@clusters$z == i), cluster.n.to.sample[i])] = FALSE

-          	}

-            }   

-            cell.ix = which(z.include)

-          

-            norm = t(normalizeCounts(counts[,cell.ix], normalize="proportion", 

-                                     transformation.fun=sqrt))

-            res = calculateTsne(norm, perplexity=perplexity, max.iter=max.iter, 

-                                seed=seed, do.pca=TRUE, 

-                                initial.dims = initial.dims)

+            res = calculateTsne(prepared.count.info$norm, perplexity=perplexity, max.iter=max.iter, 

+                                seed=seed, do.pca=TRUE, initial.dims = initial.dims)

             final = matrix(NA, nrow=ncol(counts), ncol=2)

-            final[cell.ix,] = res

+            final[prepared.count.info$cell.ix,] = res

             rownames(final) = colnames(counts)

             colnames(final) = c("tsne_1", "tsne_2")

             return(final)

           })

+#' @title umap for celda_C

+#' @description Embeds cells in two dimensions using umap based on a `celda_C` model. PCA on the normalized counts is used to reduce the number of features before applying umap. 

+#' 

+#' @param counts Integer matrix. Rows represent features and columns represent cells. This matrix should be the same as the one used to generate `celda.mod`.

+#' @param celda.mod Celda object of class `celda_C`. 

+#' @param max.cells Integer. Maximum number of cells to plot. Cells will be randomly subsampled if ncol(counts) > max.cells. Larger numbers of cells requires more memory. Default 25000.

+#' @param min.cluster.size Integer. Do not subsample cell clusters below this threshold. Default 100. 

+#' @param initial.dims Integer. PCA will be used to reduce the dimentionality of the dataset. The top 'initial.dims' principal components will be used for tSNE. Default 20.

+#' @param perplexity Numeric. Perplexity parameter for tSNE. Default 20.

+#' @param max.iter Integer. Maximum number of iterations in tSNE generation. Default 2500.

+#' @param seed Integer. Passed to `set.seed()`. Default 12345. If NULL, no calls to `set.seed()` are made.

+#' @param ... Additional parameters.

+#' @seealso `celda_C()` for clustering cells and `celdaHeatmap()` for displaying expression

+#' @examples

+#' umap.res = celdaUmap(celda.C.sim$counts, celda.C.mod)

+#' @return A two column matrix of t-SNE coordinates

+#' @export

+setMethod("celdaUmap",

+          signature(celda.mod = "celda_C"),

+          function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                   modules=NULL, umap.config=umap::umap.defaults) {

+            prepared.count.info = prepareCountsForDimReduction.celda_C(counts, celda.mod, max.cells,

+                                                                       min.cluster.size, modules)

+            res = calculateUmap(prepared.count.info$norm, umap.config)

+            final = matrix(NA, nrow=ncol(counts), ncol=2)

+            final[prepared.count.info$cell.ix,] = res

+            rownames(final) = colnames(counts)

+            colnames(final) = c("umap_1", "umap_2")

+            return(final)

+          })

+

+

+prepareCountsForDimReduction.celda_C = function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                                                modules=NULL) {

+  counts = processCounts(counts)

+  compareCountMatrix(counts, celda.mod)

+  

+  ## Checking if max.cells and min.cluster.size will work

+  if((max.cells < ncol(counts)) & (max.cells / min.cluster.size < celda.mod@params$K)) {

+  stop(paste0("Cannot distribute ", max.cells, " cells among ", 

+              celda.mod@params$K, " clusters while maintaining a minumum of ", 

+              min.cluster.size, 

+              " cells per cluster. Try increasing 'max.cells' or decreasing 'min.cluster.size'."))

+  }

+  

+  ## Select a subset of cells to sample if greater than 'max.cells'

+  total.cells.to.remove = ncol(counts) - max.cells

+  z.include = rep(TRUE, ncol(counts))

+  if(total.cells.to.remove > 0) {

+    z.ta = tabulate(celda.mod@clusters$z, celda.mod@params$K)

+    

+    ## Number of cells that can be sampled from each cluster without 

+    ## going below the minimum threshold

+    cluster.cells.to.sample = z.ta - min.cluster.size          

+    cluster.cells.to.sample[cluster.cells.to.sample < 0] = 0

+    

+    ## Number of cells to sample after exluding smaller clusters

+    ## Rounding can cause number to be off by a few, so ceiling is 

+    ## used with a second round of subtraction

+    cluster.n.to.sample = ceiling((cluster.cells.to.sample / sum(cluster.cells.to.sample)) * total.cells.to.remove)

+    diff = sum(cluster.n.to.sample) - total.cells.to.remove 

+    cluster.n.to.sample[which.max(cluster.n.to.sample)] = cluster.n.to.sample[which.max(cluster.n.to.sample)] - diff

+    

+    ## Perform sampling for each cluster

+    for(i in which(cluster.n.to.sample > 0)) {

+      z.include[sample(which(celda.mod@clusters$z == i), cluster.n.to.sample[i])] = FALSE

+    }

+  }   

+  cell.ix = which(z.include)

+  norm = t(normalizeCounts(counts[,cell.ix], normalize="proportion", 

+                         transformation.fun=sqrt))

+  return(list(norm=norm, cell.ix=cell.ix))

+}

+

 #' @title Probability map for a celda_C model

 #' @description Renders probability and relative expression heatmaps to visualize the relationship between cell populations and samples.

@@ -629,6 +629,7 @@ setMethod("clusterProbability",

 #' @seealso `celda_CG()` for clustering features and cells

 #' @examples

 #' perplexity = perplexity(celda.CG.sim$counts, celda.CG.mod)

+#' @import matrixStats

 #' @export

 setMethod("perplexity",

           signature(celda.mod = "celda_CG"),

@@ -738,62 +739,102 @@ setMethod("celdaTsne",

           function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

                      initial.dims=20, modules=NULL, perplexity=20, max.iter=2500, 

                      seed=12345, ...) {

-

-            ## Checking if max.cells and min.cluster.size will work

-            if((max.cells < ncol(counts)) & (max.cells / min.cluster.size < celda.mod@params$K)) {

-              stop(paste0("Cannot distribute ", max.cells, " cells among ",

-                          celda.mod@params$K, " clusters while maintaining a minumum of ", 

-                          min.cluster.size, " cells per cluster. Try increasing 'max.cells' or decreasing 'min.cluster.size'."))

-            }

-            

-            fm = factorizeMatrix(counts=counts, celda.mod=celda.mod, 

-                                 type="counts")    

-            modules.to.use = 1:nrow(fm$counts$cell)

-            if (!is.null(modules)) {

-              if (!all(modules %in% modules.to.use)) {

-                stop("'modules' must be a vector of numbers between 1 and ", 

-                     modules.to.use, ".")

-              }

-              modules.to.use = modules 

-            }

-            

-          

-            ## Select a subset of cells to sample if greater than 'max.cells'

-            total.cells.to.remove = ncol(counts) - max.cells

-            z.include = rep(TRUE, ncol(counts))

-            if(total.cells.to.remove > 0) {

-            	z.ta = tabulate(celda.mod@clusters$z, celda.mod@params$K)

-            	

-            	## Number of cells that can be sampled from each cluster without 

-            	## going below the minimum threshold

-            	cluster.cells.to.sample = z.ta - min.cluster.size          

-            	cluster.cells.to.sample[cluster.cells.to.sample < 0] = 0

-            	

-            	## Number of cells to sample after exluding smaller clusters

-            	## Rounding can cause number to be off by a few, so ceiling is used 

-            	## with a second round of subtraction

-            	cluster.n.to.sample = ceiling((cluster.cells.to.sample / sum(cluster.cells.to.sample)) * total.cells.to.remove)

-            	diff = sum(cluster.n.to.sample) - total.cells.to.remove 

-            	cluster.n.to.sample[which.max(cluster.n.to.sample)] = cluster.n.to.sample[which.max(cluster.n.to.sample)] - diff

-            

-            	## Perform sampling for each cluster

-            	for(i in which(cluster.n.to.sample > 0)) {

-            	  z.include[sample(which(celda.mod@clusters$z == i), cluster.n.to.sample[i])] = FALSE

-            	}

-            }   

-            cell.ix = which(z.include)

-          

-            norm = t(normalizeCounts(fm$counts$cell[modules.to.use,cell.ix], normalize="proportion", transformation.fun=sqrt))

+            prepared.count.info = prepareCountsForDimReduction.celda_CG(counts, celda.mod,

+                                                         max.cells, min.cluster.size,

+                                                         initial.dims, modules, ...)

+            norm = prepared.count.info$norm

             res = calculateTsne(norm, do.pca=FALSE, perplexity=perplexity, max.iter=max.iter, seed=seed)

             final = matrix(NA, nrow=ncol(counts), ncol=2)

-            final[cell.ix,] = res

+            final[prepared.count.info$cell.ix,] = res

             rownames(final) = colnames(counts)

             colnames(final) = c("tsne_1", "tsne_2")

             return(final)

           })

+#' @title umap for celda_CG

+#' @description Embeds cells in two dimensions using umap based on a `celda_CG` model. umap is run on module probabilities to reduce the number of features instead of using PCA. Module probabilities square-root trasformed before applying tSNE. 

+#' 

+#' @param counts Integer matrix. Rows represent features and columns represent cells. This matrix should be the same as the one used to generate `celda.mod`.

+#' @param celda.mod Celda object of class `celda_CG`. 

+#' @param max.cells Integer. Maximum number of cells to plot. Cells will be randomly subsampled if ncol(counts) > max.cells. Larger numbers of cells requires more memory. Default 25000.

+#' @param min.cluster.size Integer. Do not subsample cell clusters below this threshold. Default 100. 

+#' @param modules Integer vector. Determines which features modules to use for tSNE. If NULL, all modules will be used. Default NULL.

+#' @param umap.config Object of class `umap.config`. Configures parameters for umap. Default `umap::umap.defaults`

+#' @seealso `celda_CG()` for clustering features and cells  and `celdaHeatmap()` for displaying expression

+#' @examples

+#' umap.res = celdaUmap(celda.CG.sim$counts, celda.CG.mod)

+#' @return A two column matrix of t-SNE coordinates

+#' @export

+setMethod("celdaUmap",

+          signature(celda.mod = "celda_CG"),

+          function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                   modules=NULL, umap.config=umap::umap.defaults) {

+            prepared.count.info = prepareCountsForDimReduction.celda_CG(counts, celda.mod,

+                                                                        max.cells, min.cluster.size,

+                                                                        initial.dims, modules)

+            umap.res = calculateUmap(prepared.count.info$norm,

+                                     umap.config)

+            final = matrix(NA, nrow=ncol(counts), ncol=2)

+            final[prepared.count.info$cell.ix, ] = umap.res

+            rownames(final) = colnames(counts)

+            colnames(final) = c("umap_1", "umap_2")

+            return(final)

+          })

+

+

+prepareCountsForDimReduction.celda_CG =  function(counts, celda.mod, max.cells=25000, 

+                                                  min.cluster.size=100,

+                                                  initial.dims=20, modules=NULL, ...) {

+   ## Checking if max.cells and min.cluster.size will work

+    if((max.cells < ncol(counts)) & (max.cells / min.cluster.size < celda.mod@params$K)) {

+      stop(paste0("Cannot distribute ", max.cells, " cells among ",

+                  celda.mod@params$K, " clusters while maintaining a minumum of ", 

+                  min.cluster.size, " cells per cluster. Try increasing 'max.cells' or decreasing 'min.cluster.size'."))

+    }

+    

+    fm = factorizeMatrix(counts=counts, celda.mod=celda.mod, 

+                         type="counts")    

+    modules.to.use = 1:nrow(fm$counts$cell)

+    if (!is.null(modules)) {

+      if (!all(modules %in% modules.to.use)) {

+        stop("'modules' must be a vector of numbers between 1 and ", 

+             modules.to.use, ".")

+      }

+      modules.to.use = modules 

+    }

+    

+  

+    ## Select a subset of cells to sample if greater than 'max.cells'

+    total.cells.to.remove = ncol(counts) - max.cells

+    z.include = rep(TRUE, ncol(counts))

+    if(total.cells.to.remove > 0) {

+    	z.ta = tabulate(celda.mod@clusters$z, celda.mod@params$K)

+    	

+    	## Number of cells that can be sampled from each cluster without 

+    	## going below the minimum threshold

+    	cluster.cells.to.sample = z.ta - min.cluster.size          

+    	cluster.cells.to.sample[cluster.cells.to.sample < 0] = 0

+    	

+    	## Number of cells to sample after exluding smaller clusters

+    	## Rounding can cause number to be off by a few, so ceiling is used 

+    	## with a second round of subtraction

+    	cluster.n.to.sample = ceiling((cluster.cells.to.sample / sum(cluster.cells.to.sample)) * total.cells.to.remove)

+    	diff = sum(cluster.n.to.sample) - total.cells.to.remove 

+    	cluster.n.to.sample[which.max(cluster.n.to.sample)] = cluster.n.to.sample[which.max(cluster.n.to.sample)] - diff

+    

+    	## Perform sampling for each cluster

+    	for(i in which(cluster.n.to.sample > 0)) {

+    	  z.include[sample(which(celda.mod@clusters$z == i), cluster.n.to.sample[i])] = FALSE

+    	}

+    }   

+    cell.ix = which(z.include)

+  

+    norm = t(normalizeCounts(fm$counts$cell[modules.to.use,cell.ix], normalize="proportion", transformation.fun=sqrt))

+    return(list(norm=norm, cell.ix=cell.ix))

+}

+

 #' @title Probability map for a celda_CG model

 #' @description Renders probability and relative expression heatmaps to visualize the relationship between features and cell populations or cell populations and samples.

@@ -192,7 +192,7 @@ cG.calcGibbsProbY = function(counts, n.TS.by.C, n.by.TS, nG.by.TS, n.by.G, y, L,

   ix <- sample(1:nG)

   for(i in ix) {

     probs[,i] = cG_CalcGibbsProbY(index=i, counts=counts, nTSbyC=n.TS.by.C, nbyTS=n.by.TS, nGbyTS=nG.by.TS, nbyG=n.by.G, y=y, L=L, nG=nG, lg_beta=lgbeta, lg_gamma=lggamma, lg_delta=lgdelta, delta=delta)

-    if(any(is.na(probs[,i]))) browser()

+    #if(any(is.na(probs[,i]))) browser()

 	## Sample next state and add back counts

 	if(isTRUE(do.sample)) {

 	  prev.y = y[i]

@@ -603,28 +603,10 @@ setMethod("celdaTsne",

                     initial.dims=20, modules=NULL, perplexity=20, max.iter=2500, 

                     seed=12345, ...) {

-            if(max.cells > ncol(counts)) {

-              max.cells = ncol(counts)

-            }

-            fm = factorizeMatrix(counts=counts, celda.mod=celda.mod, 

-                                 type="counts")

-              

-            modules.to.use = 1:nrow(fm$counts$cell)

-            if (!is.null(modules)) {

-          	if (!all(modules %in% modules.to.use)) {

-          	  stop("'modules' must be a vector of numbers between 1 and ", 

-          	       modules.to.use, ".")

-          	}

-          	modules.to.use = modules 

-            }

-             

-            cell.ix = sample(1:ncol(counts), max.cells)

-            norm = t(normalizeCounts(fm$counts$cell[modules.to.use,cell.ix], 

-                                     normalize="proportion", 

-                                     transformation.fun=sqrt))

-          

-            res = calculateTsne(norm, do.pca=FALSE, perplexity=perplexity, 

+            prepared.count.info =  prepareCountsForDimReduction.celda_G(counts, celda.mod, max.cells, 

+                                                                        min.cluster.size, modules) 

+            res = calculateTsne(prepared.count.info$norm, do.pca=FALSE, perplexity=perplexity, 

                                 max.iter=max.iter, seed=seed)

             rownames(res) = colnames(counts)

             colnames(res) = c("tsne_1", "tsne_2")

@@ -632,6 +614,62 @@ setMethod("celdaTsne",

           })

+#' @title umap for celda_G

+#' @description Embeds cells in two dimensions using umap based on a `celda_G` model. umap is run on module probabilities to reduce the number of features instead of using PCA. Module probabilities square-root trasformed before applying tSNE. 

+#' 

+#' @param counts Integer matrix. Rows represent features and columns represent cells. This matrix should be the same as the one used to generate `celda.mod`.

+#' @param celda.mod Celda object of class `celda_CG`. 

+#' @param max.cells Integer. Maximum number of cells to plot. Cells will be randomly subsampled if ncol(counts) > max.cells. Larger numbers of cells requires more memory. Default 25000.

+#' @param min.cluster.size Integer. Do not subsample cell clusters below this threshold. Default 100. 

+#' @param modules Integer vector. Determines which features modules to use for tSNE. If NULL, all modules will be used. Default NULL.

+#' @param umap.config Object of class `umap.config`. Configures parameters for umap. Default `umap::umap.defaults`

+#' @seealso `celda_G()` for clustering features and cells  and `celdaHeatmap()` for displaying expression

+#' @examples

+#' umap.res = celdaUmap(celda.G.sim$counts, celda.G.mod)

+#' @return A two column matrix of t-SNE coordinates

+#' @export

+setMethod("celdaUmap",

+          signature(celda.mod = "celda_G"),

+          function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                   modules=NULL, umap.config=umap::umap.defaults) {

+            prepared.count.info = prepareCountsForDimReduction.celda_G(counts, celda.mod,

+                                                                       max.cells, min.cluster.size,

+                                                                       modules)

+            umap.res = calculateUmap(prepared.count.info$norm, umap.config)

+            final = matrix(NA, nrow=ncol(counts), ncol=2)

+            final[prepared.count.info$cell.ix, ] = umap.res

+            rownames(final) = colnames(counts)

+            colnames(final) = c("umap_1", "umap_2")

+            return(final)

+          })

+

+

+prepareCountsForDimReduction.celda_G = function(counts, celda.mod, max.cells=25000, min.cluster.size=100,

+                                                modules=NULL) {

+  if(max.cells > ncol(counts)) {

+              max.cells = ncol(counts)

+            }

+            

+  fm = factorizeMatrix(counts=counts, celda.mod=celda.mod, 

+                       type="counts")

+    

+  modules.to.use = 1:nrow(fm$counts$cell)

+  if (!is.null(modules)) {

+  if (!all(modules %in% modules.to.use)) {

+    stop("'modules' must be a vector of numbers between 1 and ", 

+         modules.to.use, ".")

+  }

+  modules.to.use = modules 

+  }

+   

+  cell.ix = sample(1:ncol(counts), max.cells)

+  norm = t(normalizeCounts(fm$counts$cell[modules.to.use,cell.ix], 

+                           normalize="proportion", 

+                           transformation.fun=sqrt))

+  return(list(norm=norm, cell.ix=cell.ix))

+}

+

+

 #' @title Lookup the module of a feature

 #' @description Finds the module assignments of given features in a `celda_G()` model

 #'  

@@ -308,7 +308,7 @@ featureModuleTable = function(counts, celda.mod, output.file = NULL){

   if(is.null(output.file)){

     return(res)

   }else{

-    write.table(res, file = output.file, sep = "\t", row.names = F, quote = F)

+    utils::write.table(res, file = output.file, sep = "\t", row.names = FALSE, quote = FALSE)

   }

 }

@@ -334,7 +334,7 @@ violinPlot = function(counts, celda.mod, features){

   colnames(m) = c("Cluster","Feature","Expression")

   color_pal = distinct_colors(length(unique(cluster)))

   p <- ggplot2::ggplot(m, ggplot2::aes(x=Cluster, y=Expression, fill=Cluster)) + 

-    ggplot2::facet_wrap(~Feature) + ggplot2::geom_violin(trim=T, scale = "width") + 

+    ggplot2::facet_wrap(~Feature) + ggplot2::geom_violin(trim=TRUE, scale = "width") + 

     ggplot2::geom_jitter(height = 0, size = 0.1) +

     ggplot2::scale_fill_manual(values = color_pal)

   return(p)

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

 #' @param beta Numeric. Concentration parameter for Phi. Default to be 0.5

 #' @param delta Numeric / Numeric vector. Concentration parameter for Theta. If input as a single numeric value, symmetric values for beta distribution are specified; if input as a vector of lenght 2, the two values will be the shape1 and shape2 paramters of the beta distribution respectively

 #' @param seed Integer. Passed to set.seed(). Default to be 12345. If NULL, no calls to `set.seed()` are made.

+#' @return A list object containing the real expression matrix and contamination

+#' expression matrix as well as other parameters used in the simulation.

 #' @examples 

 #' contamination.sim =  simulateContaminatedMatrix(  K=3,  delta=c(1,9)) 

 #' contamination.sim =  simulateContaminatedMatrix(  K=3,  delta = 1) 

@@ -15,13 +17,13 @@

 simulateContaminatedMatrix = function(C=300, G=100, K=3, N.Range=c(500,1000), beta = 0.5, delta=c(1,2),  seed=12345) {

     if (!is.null(seed)) {

-      set.seed(seed) 

+        set.seed(seed) 

     }

-    if(length(delta)==1) { 

-	    cp.byC = rbeta(n=C, shape1=delta, shape2=delta) 

+    if(length(delta) == 1 ) { 

+        cp.byC = rbeta(n=C, shape1=delta, shape2=delta) 

     } else {

-            cp.byC = rbeta(n=C, shape1=delta[1], shape2=delta[2] ) 

+        cp.byC = rbeta(n=C, shape1=delta[1], shape2=delta[2] ) 

     } 

     z = sample(1:K, size=C, replace=TRUE) 

@@ -32,21 +34,21 @@ simulateContaminatedMatrix = function(C=300, G=100, K=3, N.Range=c(500,1000), be

     }

-    N.byC = sample(min(N.Range):max(N.Range), size=C , replace=TRUE  ) 

+    N.byC = sample(min(N.Range):max(N.Range), size=C, replace=TRUE  ) 

     cN.byC = sapply(1:C, function(i)  rbinom(n=1, size=N.byC[i], p=cp.byC[i] )  ) 

     rN.byC = N.byC - cN.byC 

     phi = rdirichlet(K, rep(beta, G)) 

-    # sample real expressed count matrix 

+    ## sample real expressed count matrix 

     cell.rmat = sapply(1:C, function(i) stats::rmultinom(1, size=rN.byC[i], prob=phi[z[i],] )  )

     rownames(cell.rmat) = paste0("Gene_", 1:G) 

     colnames(cell.rmat) = paste0("Cell_", 1:C) 

-    # sample contamination count matrix 

+    ## sample contamination count matrix 

     n.G.by.K = rowSums(cell.rmat) - colSumByGroup(cell.rmat, group=z, K=K) 

     eta = normalizeCounts(counts=n.G.by.K, normalize="proportion") 

@@ -67,17 +69,17 @@ simulateContaminatedMatrix = function(C=300, G=100, K=3, N.Range=c(500,1000), be

 # eta Numeric matrix. Rows represent features and columns represent cell populations 

 # theta Numeric vector. Proportion of truely expressed transcripts

 decon.calcLL = function(counts, z,  phi, eta, theta){

-  #ll = sum( t(counts) * log( (1-conP )*geneDist[z,] + conP * conDist[z, ] + 1e-20 ) )  # when dist_mat are K x G matrices

-  ll = sum( t(counts) * log( theta * t(phi)[z,] + (1-theta) * t(eta)[z,]  + 1e-20  ))   

-  return(ll)

+    #ll = sum( t(counts) * log( (1-conP )*geneDist[z,] + conP * conDist[z, ] + 1e-20 ) )  # when dist_mat are K x G matrices

+    ll = sum( t(counts) * log( theta * t(phi)[z,] + (1-theta) * t(eta)[z,]  + 1e-20  ))   

+    return(ll)

 }

 # This function calculates the log-likelihood of background distribution decontamination

 #

 # bgDist Numeric matrix. Rows represent feature and columns are the times that the background-distribution has been replicated. 

 bg.calcLL = function(counts, cellDist, bgDist, theta){

-  ll = sum( t(counts) * log( theta * t(cellDist) + (1-theta) * t(bgDist) + 1e-20) ) 

-  return(ll) 

+    ll = sum( t(counts) * log( theta * t(cellDist) + (1-theta) * t(bgDist) + 1e-20) ) 

+    return(ll) 

 }

@@ -88,24 +90,24 @@ bg.calcLL = function(counts, cellDist, bgDist, theta){

 #  theta Numeric vector. Proportion of truely expressed transctripts

 cD.calcEMDecontamination = function(counts, phi, eta, theta,  z, K, beta, delta ) {

-   # Notes: use fix-point iteration to update prior for theta, no need to feed delta anymore

-   log.Pr = log(  t(phi)[z,] + 1e-20) + log(theta + 1e-20 )

-   log.Pc = log(  t(eta)[z,] + 1e-20) + log(1-theta + 1e-20)  

+    ## Notes: use fix-point iteration to update prior for theta, no need to feed delta anymore

+    log.Pr = log(  t(phi)[z,] + 1e-20) + log(theta + 1e-20 )

+    log.Pc = log(  t(eta)[z,] + 1e-20) + log(1-theta + 1e-20)  

-   Pr = exp(log.Pr) / (  exp(log.Pr) + exp(log.Pc) )   

-   Pc = 1 - Pr 

-   delta.v2 = MCMCprecision::fit_dirichlet( matrix( c( Pr, Pc) , ncol = 2 ) )$alpha 

-   

-   est.rmat = t(Pr) * counts 

-   rn.G.by.K = colSumByGroup.numeric(est.rmat, z, K)  

-   cn.G.by.K = rowSums(rn.G.by.K) - rn.G.by.K 

+    Pr = exp(log.Pr) / (  exp(log.Pr) + exp(log.Pc) )   

+    Pc = 1 - Pr 

+    delta.v2 = MCMCprecision::fit_dirichlet( matrix( c( Pr, Pc) , ncol = 2 ) )$alpha 

-   #update parameters 

-   theta = (colSums(est.rmat) + delta.v2[1]) / (colSums(counts) + sum(delta.v2)  )  

-   phi  = normalizeCounts(rn.G.by.K, normalize="proportion",  pseudocount.normalize =beta ) 

-   eta  = normalizeCounts(cn.G.by.K, normalize="proportion",  pseudocount.normalize = beta) 

+    est.rmat = t(Pr) * counts 

+    rn.G.by.K = colSumByGroup.numeric(est.rmat, z, K)  

+    cn.G.by.K = rowSums(rn.G.by.K) - rn.G.by.K 

-   return( list("est.rmat"=est.rmat,  "theta"=theta, "phi"=phi, "eta"=eta, "delta"=delta.v2 ) ) 

+    ## Update parameters 

+    theta = (colSums(est.rmat) + delta.v2[1]) / (colSums(counts) + sum(delta.v2)  )  

+    phi  = normalizeCounts(rn.G.by.K, normalize="proportion",  pseudocount.normalize =beta ) 

+    eta  = normalizeCounts(cn.G.by.K, normalize="proportion",  pseudocount.normalize = beta) 

+

+    return( list("est.rmat"=est.rmat,  "theta"=theta, "phi"=phi, "eta"=eta, "delta"=delta.v2 ) ) 

 }

@@ -113,11 +115,11 @@ cD.calcEMDecontamination = function(counts, phi, eta, theta,  z, K, beta, delta

 # 

 cD.calcEMbgDecontamination = function(counts, cellDist, bgDist, theta, beta, delta){

-  meanN.by.C = apply(counts, 2, mean) 

-  log.Pr = log( t(cellDist) + 1e-20) + log( theta + 1e-20) # + log( t(counts) / meanN.by.C )   # better when without panelty 

-  log.Pc = log( t(bgDist)  + 1e-20) + log( 1-theta + 2e-20) 

+    meanN.by.C = apply(counts, 2, mean) 

+    log.Pr = log( t(cellDist) + 1e-20) + log( theta + 1e-20) # + log( t(counts) / meanN.by.C )   # better when without panelty 

+    log.Pc = log( t(bgDist)  + 1e-20) + log( 1-theta + 2e-20) 

-  Pr = exp( log.Pr) / (exp(log.Pr) + exp( log.Pc) ) 

+    Pr = exp( log.Pr) / (exp(log.Pr) + exp( log.Pc) ) 

   est.rmat = t(Pr) * counts

@@ -125,7 +127,7 @@ cD.calcEMbgDecontamination = function(counts, cellDist, bgDist, theta, beta, del

   theta = (colSums(est.rmat) + delta) / (colSums(counts) + 2*delta) 

   cellDist = normalizeCounts(est.rmat, normalize="proportion", pseudocount.normalize =beta) 

-  return( list("est.rmat"=est.rmat, "theta"=theta, "cellDist"=cellDist) ) 

+    return( list("est.rmat"=est.rmat, "theta"=theta, "cellDist"=cellDist) ) 

 }

@@ -139,48 +141,48 @@ cD.calcEMbgDecontamination = function(counts, cellDist, bgDist, theta, beta, del

 #' @param logfile Character. Messages will be redirected to a file named `logfile`. If NULL, messages will be printed to stdout.  Default NULL

 #' @param verbose Logical. Whether to print log messages. Default TRUE

 #' @param seed Integer. Passed to set.seed(). Default to be 1234567. If NULL, no calls to `set.seed()` are made.

+#' @return A list object which contains the decontaminated count matrix and related parameters

 #' @examples 

 #' decon.c = DecontX( counts = contamination.sim$rmat + contamination.sim$cmat, z=contamination.sim$z, max.iter=3)

 #' decon.bg = DecontX( counts=contamination.sim$rmat + contamination.sim$cmat, max.iter=3 ) 

 #' @export 

-DecontX = function( counts , z=NULL, batch=NULL,  max.iter=200, beta=1e-6, delta=10, logfile=NULL, verbose=TRUE, seed=1234567 ) {

+DecontX = function( counts, z=NULL, batch=NULL,  max.iter=200, beta=1e-6, delta=10, logfile=NULL, verbose=TRUE, seed=1234567 ) {

-  if ( !is.null(batch) ) { 

-      

-      # Set result lists upfront for all cells from different batches 

-      logLikelihood = c() 

-      est.rmat = matrix(NA, ncol = ncol(counts) , nrow = nrow(counts ) , dimnames = list( rownames( counts) , colnames( counts ) )      ) 

-      theta = rep(NA, ncol(counts ) ) 

-      est.conp = rep(NA, ncol(counts) ) 

-

-      batch.index = unique(batch) 

-

-      for ( BATCH in batch.index ) {  

-          counts.BATCH = counts[, batch == BATCH  ]

-          if( !is.null(z) ) { z.BATCH = z[ batch == BATCH ]  }  else { z.BATCH = z } 

-          res.BATCH = DecontXoneBatch( counts = counts.BATCH, z = z.BATCH, batch = BATCH  ,max.iter = max.iter, beta = beta, delta = delta, logfile=logfile, verbose=verbose, seed=seed ) 

-

-          est.rmat[, batch == BATCH ] = res.BATCH$res.list$est.rmat

-          est.conp[ batch == BATCH ] = res.BATCH$res.list$est.conp 

-          theta[  batch == BATCH ] = res.BATCH$res.list$theta 

-

-          if( is.null(logLikelihood)  )   { 

-              logLikelihood = res.BATCH$res.list$logLikelihood  

-          } else { 

-              logLikelihood = addLogLikelihood( logLikelihood, res.BATCH$res.list$logLikelihood ) 

-          }

-      } 

-      

-      run.params = res.BATCH$run.params 

-      method = res.BATCH$method 

-      res.list = list( "logLikelihood"=logLikelihood, "est.rmat"=est.rmat,"est.conp"= est.conp,  "theta"=theta ) 

+    if ( !is.null(batch) ) { 

-      return( list("run.params"=run.params, "res.list"=res.list, "method"=method )  ) 

-  } 

+        ## Set result lists upfront for all cells from different batches 

+        logLikelihood = c() 

+        est.rmat = matrix(NA, ncol = ncol(counts), nrow = nrow(counts ), dimnames = list( rownames( counts), colnames( counts ) )      ) 

+        theta = rep(NA, ncol(counts ) ) 

+        est.conp = rep(NA, ncol(counts) ) 

-  return( DecontXoneBatch( counts = counts, z = z, max.iter = max.iter, beta = beta, delta = delta, logfile=logfile, verbose=verbose, seed=seed )  ) 

+        batch.index = unique(batch) 

+        for ( BATCH in batch.index ) {  

+            counts.BATCH = counts[, batch == BATCH  ]

+            if( !is.null(z) ) { z.BATCH = z[ batch == BATCH ]  }  else { z.BATCH = z } 

+            res.BATCH = DecontXoneBatch( counts = counts.BATCH, z = z.BATCH, batch = BATCH  ,max.iter = max.iter, beta = beta, delta = delta, logfile=logfile, verbose=verbose, seed=seed ) 

+

+            est.rmat[, batch == BATCH ] = res.BATCH$res.list$est.rmat

+            est.conp[ batch == BATCH ] = res.BATCH$res.list$est.conp 

+            theta[  batch == BATCH ] = res.BATCH$res.list$theta 

+

+            if( is.null(logLikelihood)  )   { 

+                logLikelihood = res.BATCH$res.list$logLikelihood  

+            } else { 

+                logLikelihood = addLogLikelihood( logLikelihood, res.BATCH$res.list$logLikelihood ) 

+            }

+        } 

+      

+        run.params = res.BATCH$run.params 

+        method = res.BATCH$method 

+        res.list = list( "logLikelihood"=logLikelihood, "est.rmat"=est.rmat,"est.conp"= est.conp,  "theta"=theta ) 

+

+        return( list("run.params"=run.params, "res.list"=res.list, "method"=method )  ) 

+    } 

+

+    return( DecontXoneBatch( counts = counts, z = z, max.iter = max.iter, beta = beta, delta = delta, logfile=logfile, verbose=verbose, seed=seed )  ) 

 }

@@ -188,185 +190,186 @@ DecontX = function( counts , z=NULL, batch=NULL,  max.iter=200, beta=1e-6, delta

 # This function updates decontamination for one batch  

 #

-DecontXoneBatch = function(counts, z=NULL, batch= NULL,  max.iter=200, beta=1e-6, delta=10, logfile=NULL, verbose=TRUE, seed=1234567 ) {

+DecontXoneBatch = function(counts, z=NULL, batch=NULL, max.iter=200, beta=1e-6, delta=10, logfile=NULL, verbose=TRUE, seed=1234567 ) {

-  checkCounts.decon(counts)

-  checkParameters.decon(proportion.prior=delta, distribution.prior=beta) 

+    checkCounts.decon(counts)

+    checkParameters.decon(proportionPrior=delta, distributionPrior=beta) 

-  nG = nrow(counts)

-  nC = ncol(counts)

-  K =  length(unique(z)) 

+    nG = nrow(counts)

+    nC = ncol(counts)

+    K =  length(unique(z)) 

-  if ( is.null(z) ) {

-      decon.method = "background"

-  } else {

-      decon.method = "clustering"

-      z = processCellLabels(z, num.cells= nC) 

-  }

+    if ( is.null(z) ) {

+        decon.method = "background"

+    } else {

+        decon.method = "clustering"

+        z = processCellLabels(z, num.cells= nC) 

+    }

-  # Set seed for initialization 

-  if (!is.null(seed)) {

-    set.seed(seed)

-  }

+    ## Set seed for initialization 

+    if (!is.null(seed)) {

+        set.seed(seed)

+    }

-  iter = 1L 

-  num.iter.without.improvement = 0L

-  stop.iter = 3L 

+    iter = 1L 

+    num.iter.without.improvement = 0L

+    stop.iter = 3L 

-  logMessages("----------------------------------------------------------------------", logfile=logfile, append=TRUE, verbose=verbose) 

-  logMessages("Start DecontX. Decontamination", logfile=logfile, append=TRUE, verbose=verbose) 

-  if ( !is.null(batch) ) {  logMessages("batch: ",  batch, logfile=logfile, append=TRUE, verbose=verbose)    }

-  logMessages("----------------------------------------------------------------------", logfile=logfile, append=TRUE, verbose=verbose) 

-  start.time = Sys.time()

+    logMessages("----------------------------------------------------------------------", logfile=logfile, append=TRUE, verbose=verbose) 

+    logMessages("Start DecontX. Decontamination", logfile=logfile, append=TRUE, verbose=verbose) 

+    if ( !is.null(batch) ) {  logMessages("batch: ",  batch, logfile=logfile, append=TRUE, verbose=verbose)    }

+    logMessages("----------------------------------------------------------------------", logfile=logfile, append=TRUE, verbose=verbose) 

+    start.time = Sys.time()

-  if( decon.method == "clustering") {

+    if( decon.method == "clustering") {

-  # Initialization

-  theta  = runif(nC, min = 0.1, max = 0.5)  

-  est.rmat = t (t(counts) * theta )       

-  phi =   colSumByGroup.numeric(est.rmat, z, K)

-  eta =   rowSums(phi) - phi 

-  phi = normalizeCounts(phi, normalize="proportion", pseudocount.normalize =beta )

-  eta = normalizeCounts(eta, normalize="proportion", pseudocount.normalize = beta)

-  ll = c()

+        ## Initialization

+        theta  = stats::runif(nC, min = 0.1, max = 0.5)  

+        est.rmat = t (t(counts) * theta )       

+        phi =   colSumByGroup.numeric(est.rmat, z, K)

+        eta =   rowSums(phi) - phi