old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -21,6 +21,7 @@ S3method(simulateCells,celda_G)

 S3method(visualizeModelPerformance,celda_C)

 S3method(visualizeModelPerformance,celda_CG)

 S3method(visualizeModelPerformance,celda_G)

+export(GiniPlot)

 export(available_models)

 export(calculateLoglikFromVariables)

 export(calculatePerplexity)

@@ -37,6 +38,7 @@ export(clusterProbability.celda_G)

 export(compareCountMatrix)

 export(completeClusterHistory)

 export(completeLogLikelihood)

+export(diffExp_MAST)

 export(distinct_colors)

 export(factorizeMatrix)

 export(finalClusterAssignment)

@@ -56,8 +58,10 @@ export(renderInteractiveKLPlot)

 export(renderIterationLikelihoodPlot)

 export(runParams)

 export(seed)

+export(selectModels)

 export(semi_pheatmap)

 export(simulateCells)

+export(stateHeatmap)

 export(topRank)

 export(visualizeModelPerformance)

 import(RColorBrewer)

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -19,7 +19,7 @@ available_models = c("celda_C", "celda_G", "celda_CG")

 #' @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, cores=1, seed=12345, verbose=TRUE, logfile_prefix="Celda", ...) {

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

   message("Starting celda...")

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

old mode 100644

new mode 100755

@@ -204,11 +204,9 @@ celda_C = function(counts, sample.label=NULL, K, alpha=1, beta=1,

     iter = iter + 1    

   }

-  reordered.labels = reorder.label.by.size(z.best, K)

-  z.final.reorder = reordered.labels$new.labels

   names = list(row=rownames(counts), column=colnames(counts), sample=levels(sample.label))

-  result = list(z=z.final.reorder, completeLogLik=ll,  

+  result = list(z=z.best, completeLogLik=ll,  

                 finalLogLik=ll.best, seed=seed, K=K, 

                 sample.label=sample.label, alpha=alpha, 

                 beta=beta, count.checksum=count.checksum, 

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

   class(result) = "celda_C"

+  result = reorder.celdaC(counts = counts, res = result)

+  

   return(result)

 }

@@ -351,6 +351,16 @@ cC.calcLL = function(m.CP.by.S, n.CP.by.G, s, z, K, nS, alpha, beta) {

   return(final)

 }

+reorder.celdaC = function(counts,res){

+  #Reorder K

+  fm <- factorizeMatrix(counts = counts, celda.mod = res)

+  fm.norm <- t(normalizeCounts(t(fm$proportions$gene.states),scale.factor = 1))

+  d <- dist(t(fm.norm),diag = TRUE, upper = TRUE)

+  h <- hclust(d, method = "complete")

+  res <- recodeClusterZ(res,from = h$order,

+                        to = c(1:ncol(fm$counts$gene.states)))

+  return(res)

+}

 #' Generate factorized matrices showing each feature's influence on the celda_C model clustering 

 #' 

old mode 100644

new mode 100755

@@ -93,7 +93,7 @@ calculateLoglikFromVariables.celda_CG = function(counts, s, z, y, K, L, alpha, b

   return(final)

 }

-cCG.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) {

+.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) {

   ## Determine if any TS has 0 genes

   ## Need to remove 0 gene states as this will cause the likelihood to fail

@@ -142,7 +142,7 @@ cCG.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS

   return(final)

 }

-cCG.calcGibbsProbZ = function(m.CP.by.S, n.CP.by.TS, n.CP, L, alpha, beta) {

+.calcGibbsProbZ = function(m.CP.by.S, n.CP.by.TS, n.CP, L, alpha, beta) {

   ## Calculate for "Theta" component

   theta.ll = log(m.CP.by.S + alpha)

@@ -158,7 +158,7 @@ cCG.calcGibbsProbZ = function(m.CP.by.S, n.CP.by.TS, n.CP, L, alpha, beta) {

   return(final)

 }

-cCG.calcGibbsProbY = function(n.CP.by.TS, n.by.TS, nG.by.TS, nG.in.Y, beta, delta, gamma) {

+.calcGibbsProbY = function(n.CP.by.TS, n.by.TS, nG.by.TS, nG.in.Y, beta, delta, gamma) {

   ## Calculate for "Phi" component

   phi.ll = sum(lgamma(n.CP.by.TS + beta))

@@ -176,6 +176,25 @@ cCG.calcGibbsProbY = function(n.CP.by.TS, n.by.TS, nG.by.TS, nG.in.Y, beta, delt

   return(final)

 }

+reorder.celdaCG = function(counts,res){

+  #Reorder K

+  fm <- factorizeMatrix(counts = counts, celda.mod = res)

+  fm.norm <- t(normalizeCounts(t(fm$proportions$population.states),scale.factor = 1))

+  d <- dist(t(fm.norm),diag = TRUE, upper = TRUE)

+  h <- hclust(d, method = "complete")

+  res <- recodeClusterZ(res,from = h$order,

+                        to = c(1:ncol(fm$counts$population.states)))

+  

+  #Reorder L

+  fm <- factorizeMatrix(counts = counts, celda.mod = res)

+  fm.norm <- t(normalizeCounts(t(fm$proportions$population.states),scale.factor = 1))

+  d <- dist((fm.norm),diag = TRUE, upper = TRUE)

+  h <- hclust(d, method = "complete")

+  res <- recodeClusterY(res,from = h$order,

+                        to = c(1:nrow(fm$counts$population.states)))

+  return(res)

+}

+

 #' Simulate cells from the cell/gene clustering generative model

 #' 

 #' @param S The number of samples

@@ -246,8 +265,8 @@ simulateCells.celda_CG = function(model, S=10, C.Range=c(50,100), N.Range=c(500,

   zero.row.idx = which(rowSums(cell.counts) == 0)

   if (length(zero.row.idx > 0)) {

     cell.counts = cell.counts[-zero.row.idx, ]

-  }

-  new$y = new$y[-zero.row.idx]

+    new$y = new$y[-zero.row.idx]

+  } 

   ## Assign gene/cell/sample names 

   rownames(cell.counts) = paste0("Gene_", 1:nrow(cell.counts))

@@ -332,7 +351,7 @@ celda_CG = function(counts, sample.label=NULL, K, L, alpha=1, beta=1,

   nG = nrow(counts)

   nM = ncol(counts)

-  ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

+  ll = .calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

   iter = 1

   continue = TRUE

@@ -358,7 +377,7 @@ celda_CG = function(counts, sample.label=NULL, K, L, alpha=1, beta=1,

         temp.n.CP = n.CP

         temp.n.CP[j] = temp.n.CP[j] + sum(counts[,i])

-        probs[j] = cCG.calcGibbsProbZ(m.CP.by.S=m.CP.by.S[j,s[i]], n.CP.by.TS=temp.n.CP.by.TS, n.CP=temp.n.CP, L=L, alpha=alpha, beta=beta)

+        probs[j] = .calcGibbsProbZ(m.CP.by.S=m.CP.by.S[j,s[i]], n.CP.by.TS=temp.n.CP.by.TS, n.CP=temp.n.CP, L=L, alpha=alpha, beta=beta)

       }  

       ## Sample next state and add back counts

@@ -410,7 +429,7 @@ celda_CG = function(counts, sample.label=NULL, K, L, alpha=1, beta=1,

         temp.nG.by.TS = nG.by.TS + (1 * ADD_PSEUDO)

         temp.nG.by.TS[j] = temp.nG.by.TS[j] + 1

-        probs[j] = cCG.calcGibbsProbY(n.CP.by.TS=temp.n.CP.by.TS,

+        probs[j] = .calcGibbsProbY(n.CP.by.TS=temp.n.CP.by.TS,

 			n.by.TS=temp.n.by.TS,

 			nG.by.TS=temp.nG.by.TS,

 			nG.in.Y=temp.nG.by.TS[j],

@@ -481,7 +500,7 @@ celda_CG = function(counts, sample.label=NULL, K, L, alpha=1, beta=1,

    }

     ## Calculate complete likelihood

-    temp.ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

+    temp.ll = .calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

     if((all(temp.ll > ll)) | iter == 1) {

       z.best = z

       y.best = y

@@ -493,14 +512,12 @@ celda_CG = function(counts, sample.label=NULL, K, L, alpha=1, beta=1,

     iter = iter + 1    

   }

-  ## Peform reordering on final Z and Y assigments:

-  reordered.labels = reorder.labels.by.size.then.counts(counts, z=z.best, 

-                                                        y=y.best, K=K, L=L)

+

   names = list(row=rownames(counts), column=colnames(counts), 

                sample=levels(sample.label))

-  result = list(z=reordered.labels$z, y=reordered.labels$y, completeLogLik=ll, 

+  result = list(z=z.best, y=y.best, completeLogLik=ll, 

                 finalLogLik=ll.best, K=K, L=L, alpha=alpha, 

                 beta=beta, delta=delta, gamma=gamma, seed=seed, 

                 sample.label=sample.label, names=names,

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

   class(result) = "celda_CG" 

+  ## Peform reordering on final Z and Y assigments:

+  result = reorder.celdaCG(counts = counts, res = result)

   return(result)

 }

@@ -637,7 +656,7 @@ clusterProbability.celda_CG = function(counts, celda.mod) {

 	  temp.n.CP = n.CP

 	  temp.n.CP[j] = temp.n.CP[j] + sum(counts[,i])

-	  z.prob[i,j] = cCG.calcGibbsProbZ(m.CP.by.S=m.CP.by.S[j,s[i]], n.CP.by.TS=temp.n.CP.by.TS, n.CP=temp.n.CP, L=L, alpha=alpha, beta=beta)

+	  z.prob[i,j] = .calcGibbsProbZ(m.CP.by.S=m.CP.by.S[j,s[i]], n.CP.by.TS=temp.n.CP.by.TS, n.CP=temp.n.CP, L=L, alpha=alpha, beta=beta)

 	}  

 	m.CP.by.S[z[i],s[i]] = m.CP.by.S[z[i],s[i]] + 1

@@ -667,7 +686,7 @@ clusterProbability.celda_CG = function(counts, celda.mod) {

 	  temp.nG.by.TS = nG.by.TS + (1 * ADD_PSEUDO)

 	  temp.nG.by.TS[j] = temp.nG.by.TS[j] + 1

-	  y.prob[i,j] = cCG.calcGibbsProbY(n.CP.by.TS=temp.n.CP.by.TS,

+	  y.prob[i,j] = .calcGibbsProbY(n.CP.by.TS=temp.n.CP.by.TS,

 		  n.by.TS=temp.n.by.TS,

 		  nG.by.TS=temp.nG.by.TS,

 		  nG.in.Y=temp.nG.by.TS[j],

old mode 100644

new mode 100755

@@ -123,7 +123,16 @@ cG.calcLL = function(n.TS.by.C, n.by.TS, n.by.G, nG.by.TS, nM, nG, L, beta, delt

   return(final)

 }

-

+reorder.celdaG = function(counts,res){

+  #Reorder L

+  fm <- factorizeMatrix(counts = counts, celda.mod = res)

+  fm.norm <- t(normalizeCounts(fm$proportions$gene.states,scale.factor = 1))

+  d <- dist((fm.norm),diag = TRUE, upper = TRUE)

+  h <- hclust(d, method = "complete")

+  res <- recodeClusterY(res,from = h$order,

+                        to = c(1:nrow(fm$counts$cell.states)))

+  return(res)

+}

 # Calculate Log Likelihood For Single Set of Cluster Assignments (Gene Clustering)

@@ -297,16 +306,16 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1, max.iter=25,

     iter <- iter + 1    

   }

-  reordered.labels = reorder.label.by.size(y.best, L)

-  y.final.order = reordered.labels$new.labels

   names = list(row=rownames(counts), column=colnames(counts))  

-  result = list(y=y.final.order, completeLogLik=ll, 

+  result = list(y=y.best, completeLogLik=ll, 

                 finalLogLik=ll.best, L=L, beta=beta, delta=delta, gamma=gamma,

                 count.checksum=count.checksum, seed=seed, names=names)

   class(result) = "celda_G"

+  result = reorder.celdaG(counts = counts, res = result)

+  

   return(result)

 }

@@ -361,8 +370,10 @@ simulateCells.celda_G = function(model, C=100, N.Range=c(500,5000),  G=1000,

   ## Ensure that there are no all-0 rows in the counts matrix, which violates a celda modeling

   ## constraint (columns are guarnteed at least one count):

   zero.row.idx = which(rowSums(cell.counts) == 0)

-  cell.counts = cell.counts[-zero.row.idx, ]

-  y = y[-zero.row.idx]

+  if (length(zero.row.idx > 0)) {

+    cell.counts = cell.counts[-zero.row.idx, ]

+    y = y[-zero.row.idx]

+  }

   rownames(cell.counts) = paste0("Gene_", 1:nrow(cell.counts))

   colnames(cell.counts) = paste0("Cell_", 1:ncol(cell.counts))

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -2,25 +2,26 @@

 # S3 Methods for celda_list objects                                            #

 ################################################################################

-#' Run the celda Bayesian hierarchical model on a matrix of counts.

-#' TODO: - If no chain provided, automatically choose best chain

-#'       - Smarter subsetting of the run.params to DRY this function up

+# TODO: - If no chain provided, automatically choose best chain

+#       - Smarter subsetting of the run.params to DRY this function up

+

+#' Select the best model amongst a list of models generated in a celda run.

 #' 

 #' @param celda.list A celda_list object returned from celda()

-#' @param K The K parameter for the desired model in the results list

-#' @param L The L parameter for the desired model in the results list

-#' @param chain The desired chain for the specified model

-#' @param best Method for choosing best chain automatically. Options are c("perplexity", "loglik"). See documentation for chooseBestModel for details. Overrides chain parameter if provided.

+#' @param K The K parameter for the desired model in the results list. Defaults to NULL.

+#' @param L The L parameter for the desired model in the results list. Defaults to NULL.

+#' @param chain The desired chain for the specified model. Defaults to NULL. Overrides best parameter if provided.

+#' @param best Method for choosing best chain automatically. Options are c("perplexity", "loglik"). See documentation for chooseBestModel for details. Defaults to "loglik."

 #' @return A celda model object matching the provided parameters (of class "celda_C", "celda_G", "celda_CG" accordingly), or NA if one is not found.

 #' @export

-getModel = function(celda.list, K=NULL, L=NULL, chain=1, best=NULL) {

+getModel = function(celda.list, K=NULL, L=NULL, chain=NULL, best="loglik") {

   validateGetModelParams(celda.list, K, L, chain, best)  # Sanity check params

   requested.chain = NA

   run.params = celda.list$run.params

   if (celda.list$content.type == "celda_CG") {

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

+    if (is.null(chain)) {

       matching.chain.idx = run.params[run.params$K == K & run.params$L == L, "index"]

       requested.chain = chooseBestChain(celda.list$res.list[matching.chain.idx], best)

     } else {

@@ -32,7 +33,7 @@ getModel = function(celda.list, K=NULL, L=NULL, chain=1, best=NULL) {

   if (celda.list$content.type == "celda_C") {

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

+    if (is.null(chain)) {

       matching.chain.idx = run.params[run.params$K == K, "index"]

       requested.chain = chooseBestChain(celda.list$res.list[matching.chain.idx], best)

     } else {

@@ -43,7 +44,7 @@ getModel = function(celda.list, K=NULL, L=NULL, chain=1, best=NULL) {

   if (celda.list$content.type == "celda_G") {

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

+    if (is.null(chain)) {

       matching.chain.idx = run.params[run.params$L == L, "index"]

       requested.chain = chooseBestChain(celda.list$res.list[matching.chain.idx], best)

     } else { 

@@ -62,6 +63,30 @@ getModel = function(celda.list, K=NULL, L=NULL, chain=1, best=NULL) {

 }

+#' Select a set of models from a celda_list objects based off of rows in its run.params attribute.

+#' 

+#' @param celda.list A celda_list object returned from celda()

+#' @param run.param.rows Row indices in the celda.list's run params corresponding to the desired models.

+#' @return A celda_list containing celda model objects matching the provided parameters (of class "celda_C", "celda_G", "celda_CG" accordingly), or NA if one is not found.

+#' @export

+selectModels = function(celda.list, run.param.rows) {

+  desired.models = lapply(run.param.rows,

+                          function(row.idx) {

+                            if (!is.numeric(row.idx) | 

+                                row.idx < 0 | 

+                                row.idx > nrow(celda.list$run.params)) {

+                                  stop("Invalid row index provided") 

+                            }

+                            params = celda.list$run.params[row.idx, ]

+                            getModel(celda.list, params$K, params$L, params$chain)

+                          })

+  subsetted.celda.list = list(run.params=celda.list$run.params[run.param.rows, ],

+                              res.list=desired.models, content.type=celda.list$content.type,

+                              count.checksum=celda.list$count.checksum)

+  return(subsetted.celda.list)

+}

+

+

 validateGetModelParams = function(celda.list, K, L, chain, best) {

   if (class(celda.list) != "celda_list") stop("First argument to getModel() should be an object of class 'celda_list'")

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -17,8 +17,8 @@ plotDrGrid <- function(dim1, dim2, matrix, size, xlab, ylab, color_low, color_mi

   colnames(m) <- c(xlab,ylab,"facet",var_label)

   ggplot2::ggplot(m, ggplot2::aes_string(x=xlab, y=ylab)) + ggplot2::geom_point(stat = "identity", size = size, ggplot2::aes_string(color = var_label)) + 

     ggplot2::facet_wrap(~facet) + ggplot2::theme_bw() + ggplot2::scale_colour_gradient2(low = color_low, high = color_high, mid = color_mid, midpoint = (max(m[,4])-min(m[,4]))/2 ,name = gsub("_"," ",var_label)) + 

-    ggplot2::theme(strip.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.spacing = unit(0,"lines"),

-                   panel.background = element_blank(), axis.line = ggplot2::element_line(colour = "black"))

+    ggplot2::theme(strip.background = ggplot2::element_blank(), panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), panel.spacing = unit(0,"lines"),

+                   panel.background = ggplot2::element_blank(), axis.line = ggplot2::element_line(colour = "black"))

 }

 #' Create a scatterplot for each row of a normalized gene expression matrix where x and y axis are from a data dimensionality reduction tool.  

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -30,7 +30,7 @@ An analysis example using celda with RNASeq via vignette('celda-analysis')

 ## New Features and announcements

-The v0.2 release of celda represents a useable implementation of the various celda clustering models.

+The v0.3 release of celda represents a useable implementation of the various celda clustering models.

 Please submit any usability issues or bugs to the issue tracker at https://github.com/compbiomed/celda

 You can discuss celda, or ask the developers usage questions, in our [Google Group.](https://groups.google.com/forum/#!forum/celda-list)

deleted file mode 100644

@@ -1 +0,0 @@

-20892e2fcfc0c39f6f9e6333e5fd3523849044c3

\ No newline at end of file

old mode 100644

new mode 100755

old mode 100644

new mode 100755

@@ -14,9 +14,11 @@ library(celda)

 library(Matrix)

 library(gtools)

 library(ggplot2)

+library(Rtsne)

+library(reshape2)

 ## ------------------------------------------------------------------------

-sim_counts = simulateCells("celda_CG", K=3,L=10)

+sim_counts = simulateCells("celda_CG", K = 3, L = 10)

 ## ------------------------------------------------------------------------

 dim(sim_counts$counts)

@@ -35,11 +37,11 @@ table(sim_counts$z, sim_counts$sample.label)

 ## ---- fig.show='hold', warning = FALSE, message = FALSE------------------

-celdaCG.res = celda_CG(counts = sim_counts$counts, K =3, L = 10, max.iter = 10, sample.label = sim_counts$sample.label)

+celda.res = celda(counts = sim_counts$counts, model = "celda_CG", K = 3, L = 10, max.iter = 10, cores = 1, nchains = 1)

 ## ------------------------------------------------------------------------

-z = celdaCG.res$z

-y = celdaCG.res$y

+z = celda.res$res.list[[1]]$z

+y = celda.res$res.list[[1]]$y

 table(z, sim_counts$z)

 table(y, sim_counts$y)

@@ -47,28 +49,41 @@ table(y, sim_counts$y)

 ## ---- fig.width = 7, fig.height = 7, warning = FALSE, message = FALSE----

 norm.counts <- normalizeCounts(sim_counts$counts)

-render_celda_heatmap(counts=norm.counts,z=z,y=y, cluster.column = FALSE, cluster.row = FALSE)

+renderCeldaHeatmap(counts = norm.counts, z=z, y=y, normalize = NULL, color_scheme = "divergent",cluster_gene = TRUE, cluster_cell = TRUE)

 ## ------------------------------------------------------------------------

-factorized <- factorizeMatrix.celda_CG(counts = sim_counts$counts, celda.obj = celdaCG.res)

+model <- getModel(celda.res, K = 3, L = 10, best = "loglik")

+factorized <- factorizeMatrix(model, sim_counts$count)

-pop.states = t(factorized$proportions$population.states)

+## ------------------------------------------------------------------------

+dim(factorized$proportions$gene.states)

+head(factorized$proportions$gene.states)

+

+## ------------------------------------------------------------------------

+dim(factorized$proportions$cell.states)

+factorized$proportions$cell.states[1:10,1:3]

+

+## ------------------------------------------------------------------------

+pop.states = factorized$proportions$population.states

+dim(pop.states)

 head(pop.states)

 ## ---- fig.width = 7, fig.height = 7--------------------------------------

-render_celda_heatmap(pop.states, col=colorRampPalette(c("white", "blue"))(70), scale.row=F, show_cellnames=T, show_genenames=T, cluster.column=F, cluster.row=F, breaks = NA, z= 1:3, y = 1:10)

+data.pca <- prcomp(t(scale(t(factorized$proportions$cell.states))),scale = F, center = F)

+

+plotDrCluster(dim1 = data.pca$rotation[,1], dim2 = data.pca$rotation[,2], cluster = celda.res$res.list[[1]]$z)

+

+plotDrState(dim1 = data.pca$rotation[,1], dim2 = data.pca$rotation[,2], matrix = factorized$proportions$cell.states, rescale = TRUE)

+

+## ---- fig.width = 7, fig.height = 7--------------------------------------

+renderCeldaHeatmap(pop.states, color_scheme = "sequential", show_cellnames=T, show_genenames=T, breaks = NA, z= 1:ncol(pop.states), y = 1:nrow(pop.states))

 ## ---- fig.width = 7, fig.height = 7--------------------------------------

 rel.states = sweep(pop.states, 1, rowSums(pop.states), "/")

-render_celda_heatmap(rel.states, col=colorRampPalette(c("white","blue"))(70), scale.row=F, show_cellnames=T, show_genenames=T, cluster.column=F, cluster.row=F, breaks = NA, z= 1:3, y = 1:10)

+renderCeldaHeatmap(rel.states, color_scheme = "sequential", show_cellnames=T, show_genenames=T, breaks = NA, z= 1:ncol(rel.states), y = 1:nrow(rel.states))

 ## ------------------------------------------------------------------------

-celda.res.list <- celda(counts = sim_counts$counts,K = 2:4, L = 9:11, nchains = 3, max.iter = 10, model = "celda_CG")

-

-## ---- fig.width = 7, fig.height = 7, message = FALSE---------------------

-visualize.model = visualize_model_performance(celda.res.list,method="perplexity")

-

-visualize.model$K

+celda.res.list <- celda(counts = sim_counts$counts,K = 2:4, L = 9:11, nchains = 3, max.iter = 10, model = "celda_CG", cores = 1)

 ## ------------------------------------------------------------------------

 model = getModel(celda.res.list, K = 3, L = 10, best = "loglik")

@@ -79,37 +94,47 @@ library(celda)

 library(Matrix)

 library(data.table)

 library(pheatmap)

-library(biomaRt)

+library(Rtsne)

 ## ------------------------------------------------------------------------

-load(file="../data/cmatp_data.rda")

-dim(cmatp_data)

-head(rownames(cmatp_data))

-head(colnames(cmatp_data))

+dim(pbmc_data)

+head(rownames(pbmc_data))

+head(colnames(pbmc_data))

 ## ------------------------------------------------------------------------

-cmatp_select <- cmatp_data[rowSums(cmatp_data>4) > 4,]

+pbmc_select <- pbmc_data[rowSums(pbmc_data>4) > 4,]

 ## ---- message = FALSE, verbose = FALSE-----------------------------------

-cmatp.res = celda_CG(cmatp_select, sample.label = rep(1,ncol(cmatp_select)), K = 15, L = 20)

+pbmc.res = celda(pbmc_select, K = 10, L = 20, cores = 1, model = "celda_CG", nchains = 8, max.iter = 100)

+model<-getModel(pbmc.res, K = 10, L = 20, best = "loglik")

 ## ------------------------------------------------------------------------

-factorize.matrix = factorizeMatrix.celda_CG(counts = cmatp_select, celda.obj = cmatp.res)

-top.genes <- topRank(factorize.matrix$proportions$gene.states)

+factorize.matrix = factorizeMatrix(model, counts=pbmc_select)

+top.genes <- topRank(factorize.matrix$proportions$gene.states, n = 25)

 ## ------------------------------------------------------------------------

-top.genes$names$L16

+top.genes$names$L3

 ## ------------------------------------------------------------------------

-top.genes$names$L19

+top.genes$names$L4

 ## ---- fig.width = 7, fig.height = 7--------------------------------------

 top.genes.ix <- unique(unlist(top.genes$index))

-norm.cmatp.counts <- normalizeCounts(cmatp_select)

-render_celda_heatmap(norm.cmatp.counts[top.genes.ix,], z = cmatp.res$z, y = cmatp.res$y[top.genes.ix], cluster.row = FALSE, cluster.column = FALSE)

+norm.pbmc.counts <- normalizeCounts(pbmc_select)

+renderCeldaHeatmap(norm.pbmc.counts[top.genes.ix,], z = model$z, y = model$y[top.genes.ix], normalize = NULL, color_scheme = "divergent")

+

+## ---- message = FALSE----------------------------------------------------

+norm_pbmc <- normalizeCounts(factorize.matrix$counts$cell.states)

+set.seed(123)

+tsne <- Rtsne(t(norm_pbmc), pca = FALSE, max_iter = 2000)

+pbmc_tsne <- tsne$Y

 ## ---- fig.width = 7, fig.height = 7--------------------------------------

-cmatp.states = t(factorize.matrix$proportions$population.states)

-cmatp.rel.states = sweep(cmatp.states, 1, rowSums(cmatp.states), "/")

-render_celda_heatmap(cmatp.rel.states, col=colorRampPalette(c("white","blue","darkgreen","green"))(100), scale.row=F, show_cellnames=T, show_genenames=TRUE, cluster.column=FALSE, cluster.row=FALSE, breaks = NA, z= 1:15, y = 1:20)

+plotDrCluster(dim1 = pbmc_tsne[,1], dim2 = pbmc_tsne[,2], cluster = as.factor(model$z))

+

+plotDrState(dim1 = pbmc_tsne[,1], dim2 = pbmc_tsne[,2], matrix = factorize.matrix$proportions$cell.states, rescale = TRUE)

+

+marker.genes <- c("ENSG00000168685_IL7R","ENSG00000198851_CD3E","ENSG00000105374_NKG7","ENSG00000203747_FCGR3A","ENSG00000090382_LYZ","ENSG00000179639_FCER1A","ENSG00000156738_MS4A1", "ENSG00000163736_PPBP")

+gene.counts <- pbmc_select[marker.genes,]

+plotDrGene(dim1 = pbmc_tsne[,1],dim2 = pbmc_tsne[,2], matrix = gene.counts, rescale = TRUE)

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

 ---

 title: "Analyzing single-cell RNA-seq count data with Celda"

-author: "Josh D. Campbell, Sean Corbett, Masanao Yajima, Zichun Liu, Shiyi Yang, Tianwen Huan, Anahita Bahri, Yusuke Koga"

+author: "Josh D. Campbell, Masanao Yajima, Sean Corbett, Zichun Liu, Shiyi Yang, Tianwen Huan, Anahita Bahri, Zhe Wang, Yusuke Koga, Jiangyuan Liu"

 date: "`r Sys.Date()`"

 output: rmarkdown::html_vignette

 vignette: >

@@ -8,15 +8,15 @@ vignette: >

   %\VignetteEngine{knitr::rmarkdown}

   %\VignetteEncoding{UTF-8}

 ---

-

-# 1 Introduction

-

-"celda" stands for "**CE**llular **L**atent **D**irichlet **A**llocation", which is a suite of Bayesian hierarchical models and supporting functions to perform gene and cell clustering for count data generated by single cell RNA-seq platforms. This algorithm is an extension of the Latent Dirichlet Allocation (LDA) topic modeling framework that has been popular in text mining applications. Celda has advantages over other clustering frameworks:

-

-1) Celda can simultaneously cluster genes into transcriptional states and cells into subpopulations

+  

+## 1 Introduction

+  

+  "celda" stands for "**CE**llular **L**atent **D**irichlet **A**llocation", which is a suite of Bayesian hierarchical models and supporting functions to perform gene and cell clustering for count data generated by single cell RNA-seq platforms. This algorithm is an extension of the Latent Dirichlet Allocation (LDA) topic modeling framework that has been popular in text mining applications. Celda has advantages over other clustering frameworks:

+  

+1) Celda can simultaneously cluster cells into subpopulations and genes into transcriptional states, which are subgroups of genes that are likely to be expressed in similar proportions over different types of cells

 2) Celda uses count-based Dirichlet-multinomial distributions so no additional normalization is required for 3' DGE single cell RNA-seq

 3) These types of models have shown good performance with sparse data.

- 

+

 In this vignette we will demonstrate how to perform cell and gene clustering with simulated and real single cell RNA-seq data using the Bayesian hierarchical models within celda. 

@@ -57,29 +57,31 @@ library(celda)

 library(Matrix)

 library(gtools)

 library(ggplot2)

+library(Rtsne)

+library(reshape2)

 ```

-We use the built-in data generating function `simulateCells.celda_CG` to simulate a single-cell RNA-Seq dataset. `simulateCells.celda_CG` returns a list containing a count matrix where each row is a gene and each column is a cell. The K parameter determines the number of cell clusters while the L parameter determines the number of gene clusters within the dataset.

+We use the built-in data generating function `simulateCells` to simulate a single-cell RNA-Seq dataset. `simulateCells` returns a list containing a count matrix where each row is a gene and each column is a cell. The K parameter determines the number of cell clusters while the L parameter determines the number of gene clusters within the dataset.

 ```{r}

-sim_counts = simulateCells("celda_CG", K=3,L=10)

+sim_counts = simulateCells("celda_CG", K = 3, L = 10)

 ```

-The `counts` variable from `simulateCells.celda_CG` contains the counts matrix.

+The `counts` variable contains the counts matrix.

 The dimensions of counts matrix:

 ```{r}

 dim(sim_counts$counts)

 ```

-The `z` variable from `simulateCells.celda_CG` contains the cluster for each cell.

+The `z` variable contains the cluster for each cell.

 Here is the number of cells in each subpopulation:

 ```{r}

 table(sim_counts$z)

 ```

-The `y` variable from `simulateCells.celda_CG` contains the transcriptional state assignment for each gene. 

+The `y` variable contains the transcriptional state assignment for each gene. 

 Here is the number of genes in each transcriptional state:

 ```{r}

@@ -102,18 +104,18 @@ table(sim_counts$z, sim_counts$sample.label)

 ## 2.1.1 How to Run Celda

-There are currently three models within this package: `celda_C` will cluster cells, `celda_G` will cluster genes, and `celda_CG` will simultaneously cluster cells and genes. Here is the command for `celda_CG`:

+There are currently three models within this package: `celda_C` will cluster cells, `celda_G` will cluster genes, and `celda_CG` will simultaneously cluster cells and genes. The `celda` function will use these models and take in a counts matrix to create a `celda` object that contains the cluster labels for each cell and gene.

 ```{r, fig.show='hold', warning = FALSE, message = FALSE}

-celdaCG.res = celda_CG(counts = sim_counts$counts, K =3, L = 10, max.iter = 10, sample.label = sim_counts$sample.label)

+celda.res = celda(counts = sim_counts$counts, model = "celda_CG", K = 3, L = 10, max.iter = 10, cores = 1, nchains = 1)

 ```

 Here is a comparison between the true cluster labels and the estimated cluster labels.

 ```{r}

-z = celdaCG.res$z

-y = celdaCG.res$y

+z = celda.res$res.list[[1]]$z

+y = celda.res$res.list[[1]]$y

 table(z, sim_counts$z)

 table(y, sim_counts$y)

@@ -124,52 +126,71 @@ We can display the clustering results with a heatmap of the normalized counts:

 ```{r, fig.width = 7, fig.height = 7, warning = FALSE, message = FALSE}

 norm.counts <- normalizeCounts(sim_counts$counts)

-render_celda_heatmap(counts=norm.counts,z=z,y=y, cluster.column = FALSE, cluster.row = FALSE)

+renderCeldaHeatmap(counts = norm.counts, z=z, y=y, normalize = NULL, color_scheme = "divergent",cluster_gene = TRUE, cluster_cell = TRUE)

 ```

 ## 2.1.2 Matrix factorization

-`celda` can also perform matrix factorization to summarize the contribution of each transcriptional state to each cellular subpopulation. 

+`celda` can also perform matrix factorization to summarize the contribution of each transcriptional state to each cellular subpopulation. Matrix factorization decomposes the counts matrix into multiple matrices, using the given `celda` model. The gene.states contains each gene's contribution to the transcriptional state it belongs. The cell.states contains each sample's contribution to each transcriptional state. The population.states contains each transcriptional states' contribution to each of the cell states. Each of these matrices can be viewed as raw counts values or proportional values.

+

+```{r}

+model <- getModel(celda.res, K = 3, L = 10, best = "loglik")

+factorized <- factorizeMatrix(model, sim_counts$count)

+```

+

+```{r}

+dim(factorized$proportions$gene.states)

+head(factorized$proportions$gene.states)

+```

 ```{r}

-factorized <- factorizeMatrix.celda_CG(counts = sim_counts$counts, celda.obj = celdaCG.res)

+dim(factorized$proportions$cell.states)

+factorized$proportions$cell.states[1:10,1:3]

+```

+

+For instance, roughly 30% of the population in cell state K2 is in transcriptional state L1. 

-pop.states = t(factorized$proportions$population.states)

+```{r}

+pop.states = factorized$proportions$population.states

+dim(pop.states)

 head(pop.states)

 ```

-These states can also be visualized with a heatmap:

+`celda` contains several plotting tools for data dimensionality reduction tool outputs. The PCA result for the factorized matrix is plotted by the `plor_dr_cluster` and `plot_dr_state` function based off of the celda clusters and state probabilities for each cell:

 ```{r, fig.width = 7, fig.height = 7}

-render_celda_heatmap(pop.states, col=colorRampPalette(c("white", "blue"))(70), scale.row=F, show_cellnames=T, show_genenames=T, cluster.column=F, cluster.row=F, breaks = NA, z= 1:3, y = 1:10)

+data.pca <- prcomp(t(scale(t(factorized$proportions$cell.states))),scale = F, center = F)

+

+plotDrCluster(dim1 = data.pca$rotation[,1], dim2 = data.pca$rotation[,2], cluster = celda.res$res.list[[1]]$z)

+

+plotDrState(dim1 = data.pca$rotation[,1], dim2 = data.pca$rotation[,2], matrix = factorized$proportions$cell.states, rescale = TRUE)

 ```

-Sometimes it may be helpful to visualize the relative probability of each transcriptional state in each cellular subpopulation. We can do this by normalizing the rows of the probability matrix to sum up to 1:

+The transcriptional and cell states can also be visualized with a heatmap:

+

+```{r, fig.width = 7, fig.height = 7}

+renderCeldaHeatmap(pop.states, color_scheme = "sequential", show_cellnames=T, show_genenames=T, breaks = NA, z= 1:ncol(pop.states), y = 1:nrow(pop.states))

+```

+

+Sometimes it may be helpful to visualize the relative probability of each transcriptional state in each cellular subpopulation. We can do this by z-score normalizing the probability matrix:

 ```{r, fig.width = 7, fig.height = 7}

 rel.states = sweep(pop.states, 1, rowSums(pop.states), "/")

-render_celda_heatmap(rel.states, col=colorRampPalette(c("white","blue"))(70), scale.row=F, show_cellnames=T, show_genenames=T, cluster.column=F, cluster.row=F, breaks = NA, z= 1:3, y = 1:10)

+renderCeldaHeatmap(rel.states, color_scheme = "sequential", show_cellnames=T, show_genenames=T, breaks = NA, z= 1:ncol(rel.states), y = 1:nrow(rel.states))

 ```

 The relative probabilities will tend to reflect the patterns observed in the heatmap as the genes in the heatmaps are Z-score transformed and thus on a relative scale as well.

+

 ## 2.1.3 Identifying the optimal number of states

 The optimal number of K and L will likely not be known a priori. Therefore multiple choices of K and L may need to tested and compared. Additionally, `celda` is sensitive to initial start conditions, so it is good practice to run multiple chains for each combination of K and L to increase the chances of finding the most optimal solution. We have built a wrapper function that will run multiple combinations of K and L with multiple chains for each combination in parallel. 

 ```{r}

-celda.res.list <- celda(counts = sim_counts$counts,K = 2:4, L = 9:11, nchains = 3, max.iter = 10, model = "celda_CG")

-```

-

-Here, `celda` will fit the model for every combination of K between 2 and 4 and L between 9 and 11. Three chains will be run for each combination. The parameter "ncores" can be set to something more than 1 to run multiple chains in parallel.

-

-```{r, fig.width = 7, fig.height = 7, message = FALSE}

-visualize.model = visualize_model_performance(celda.res.list,method="perplexity")

-

-visualize.model$K

+celda.res.list <- celda(counts = sim_counts$counts,K = 2:4, L = 9:11, nchains = 3, max.iter = 10, model = "celda_CG", cores = 1)

 ```

-Perplexity is a measure of goodness-of-fit for the model where lower perplexity is better. The perplexity of each chain can be visualized as a function of K. In this example, K=3 is much better than K=2, but K=4 does not give much improvement over K=3. We would therefore select K=3 in this case. Note: We are still actively testing various methods for assessing model performance. 

+Here, `celda` will fit the model for every combination of K between 2 and 4 and L between 9 and 11. Three chains will be run for each combination. The parameter "cores" can be set to something more than 1 to run multiple chains in parallel.

 We can use the function "getModel" to select a particular model from the list of models once we have selected the optimal K and L.

@@ -179,6 +200,8 @@ model = getModel(celda.res.list, K = 3, L = 10, best = "loglik")

 The chain with the best log likelihood will be returned from all models that were run with K = 3 and L = 10.

+

+

 # 3. PBMC Analysis with celda

 ```{r setup, include=FALSE}

@@ -187,59 +210,72 @@ library(celda)

 library(Matrix)

 library(data.table)

 library(pheatmap)

-library(biomaRt)

+library(Rtsne)

 ```

-We will use the "cmatp_data" dataset, which is a dataset of 3000 Peripheral Blood Mononuclear Cells (PBMC), available from 10X Genomics. The rows are organized by gene names, while the columns are organized by barcodes. For this tutorial, the dataset has been slightly modified so the rownames are comprised of the Ensembl gene ID as well as the gene name. The raw dataset can be found at https://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit.

+We will use the "pbmc_data" dataset, which is a dataset of 3000 Peripheral Blood Mononuclear Cells (PBMC), available from 10X Genomics. The rows are organized by gene names, while the columns are organized by barcodes. For this tutorial, the dataset has been slightly modified so the rownames are comprised of the Ensembl gene ID as well as the gene name. The raw dataset can be found at https://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit.

 We will only include genes with at least 5 counts in a least 5 cells for this analysis to reduce computational time. However, `celda` is capable of analyzing genes with fewer counts. 

 ```{r}

-load(file="../data/cmatp_data.rda")

-dim(cmatp_data)

-head(rownames(cmatp_data))

-head(colnames(cmatp_data))

+dim(pbmc_data)

+head(rownames(pbmc_data))

+head(colnames(pbmc_data))

 ```

 ```{r}

-cmatp_select <- cmatp_data[rowSums(cmatp_data>4) > 4,]

+pbmc_select <- pbmc_data[rowSums(pbmc_data>4) > 4,]

 ```

+

 ### 3.1.1 Analysis of 10X Genomics PBMC data

 ```{r, message = FALSE, verbose = FALSE}

-cmatp.res = celda_CG(cmatp_select, sample.label = rep(1,ncol(cmatp_select)), K = 15, L = 20)

+pbmc.res = celda(pbmc_select, K = 10, L = 20, cores = 1, model = "celda_CG", nchains = 8, max.iter = 100)

+model<-getModel(pbmc.res, K = 10, L = 20, best = "loglik")

 ```

 After matrix factorization, the top genes can be selected using the `topRank` function on the "gene.states" matrix:  

 ```{r}

-factorize.matrix = factorizeMatrix.celda_CG(counts = cmatp_select, celda.obj = cmatp.res)

-top.genes <- topRank(factorize.matrix$proportions$gene.states)

+factorize.matrix = factorizeMatrix(model, counts=pbmc_select)

+top.genes <- topRank(factorize.matrix$proportions$gene.states, n = 25)

 ```

 ```{r}

-top.genes$names$L16

+top.genes$names$L3

 ```

 ```{r}

-top.genes$names$L19

+top.genes$names$L4

 ```

 We can make a heatmap of these top genes as follows: 

 ```{r, fig.width = 7, fig.height = 7}

 top.genes.ix <- unique(unlist(top.genes$index))

-norm.cmatp.counts <- normalizeCounts(cmatp_select)

-render_celda_heatmap(norm.cmatp.counts[top.genes.ix,], z = cmatp.res$z, y = cmatp.res$y[top.genes.ix], cluster.row = FALSE, cluster.column = FALSE)

+norm.pbmc.counts <- normalizeCounts(pbmc_select)

+renderCeldaHeatmap(norm.pbmc.counts[top.genes.ix,], z = model$z, y = model$y[top.genes.ix], normalize = NULL, color_scheme = "divergent")

 ```