@@ -61,22 +61,22 @@ simulateObservedMatrix = function(C=300, G=100, K=3, N.Range=c(500,1000), beta =

 # This function calculates the log-likelihood

 # 

-# omat Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells

+# counts Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells

 # z Integer vector. Cell population labels

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

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

 # theta Numeric vector. Proportion of truely expressed transcripts

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

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

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

+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)

 }

 # 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(omat, cellDist, bgDist, theta){

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

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

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

   return(ll) 

 }

@@ -86,7 +86,7 @@ bg.calcLL = function(omat, cellDist, bgDist, theta){

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

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

 #  theta Numeric vector. Proportion of truely expressed transctripts

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

+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 )

@@ -96,12 +96,12 @@ cD.calcEMDecontamination = function(omat, phi, eta, theta,  z, K, beta, delta )

    Pc = 1 - Pr 

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

-   est.rmat = t(Pr) * omat 

+   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 

    #update parameters 

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

+   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) 

@@ -111,18 +111,18 @@ cD.calcEMDecontamination = function(omat, phi, eta, theta,  z, K, beta, delta )

 # This function updates decontamination using background distribution 

 # 

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

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

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

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

+  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) ) 

-  est.rmat = t(Pr) * omat

+  est.rmat = t(Pr) * counts

   # Update paramters 

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

+  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) ) 

@@ -130,7 +130,7 @@ cD.calcEMbgDecontamination = function(omat, cellDist, bgDist, theta, beta, delta

 #' This function updates decontamination on dataset with multiple batches 

-#' @param omat Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells. 

+#' @param counts Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells. 

 #' @param z Integer vector. Cell population labels. Default NULL 

 #' @param batch Integer vector. Cell batch labels. Default NULL

 #' @param max.iter Integer. Maximum iterations of EM algorithm. Default to be 200 

@@ -140,26 +140,26 @@ cD.calcEMbgDecontamination = function(omat, cellDist, bgDist, theta, beta, delta

 #' @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.

 #' @examples 

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

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

+#' 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( omat , 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(omat) , nrow = nrow(omat ) , dimnames = list( rownames( omat) , colnames( omat ) )      ) 

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

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

+      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 ) {  

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

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

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

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

+          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 

@@ -179,7 +179,7 @@ DecontX = function( omat , z=NULL, batch=NULL,  max.iter=200, beta=1e-6, delta=1

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

   } 

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

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

 }

@@ -188,13 +188,13 @@ DecontX = function( omat , z=NULL, batch=NULL,  max.iter=200, beta=1e-6, delta=1

 # This function updates decontamination for one batch  

 #

-DecontXoneBatch = function(omat, 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(omat)

+  checkCounts.decon(counts)

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

-  nG = nrow(omat)

-  nC = ncol(omat)

+  nG = nrow(counts)

+  nC = ncol(counts)

   K =  length(unique(z)) 

   if ( is.null(z) ) {

@@ -225,7 +225,7 @@ DecontXoneBatch = function(omat, z=NULL, batch= NULL,  max.iter=200, beta=1e-6,

   # Initialization

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

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

+  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 )

@@ -233,14 +233,14 @@ DecontXoneBatch = function(omat, z=NULL, batch= NULL,  max.iter=200, beta=1e-6,

   ll = c()

-  ll.round =  decon.calcLL(omat=omat, z=z, phi = phi, eta = eta, theta=theta ) 

+  ll.round =  decon.calcLL(counts=counts, z=z, phi = phi, eta = eta, theta=theta ) 

   # EM updates

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

-    next.decon = cD.calcEMDecontamination(omat=omat, phi=phi, eta=eta, theta=theta, z=z, K=K,  beta=beta, delta=delta) 

+    next.decon = cD.calcEMDecontamination(counts=counts, phi=phi, eta=eta, theta=theta, z=z, K=K,  beta=beta, delta=delta) 

     theta = next.decon$theta 

     phi = next.decon$phi

@@ -248,7 +248,7 @@ DecontXoneBatch = function(omat, z=NULL, batch= NULL,  max.iter=200, beta=1e-6,

     delta = next.decon$delta

     # Calculate log-likelihood

-    ll.temp = decon.calcLL(omat=omat, z=z, phi = phi, eta = eta, theta=theta )

+    ll.temp = decon.calcLL(counts=counts, z=z, phi = phi, eta = eta, theta=theta )

     ll  = c(ll, ll.temp)

     ll.round = c( ll.round, round( ll.temp, 2) )  

@@ -267,27 +267,27 @@ DecontXoneBatch = function(omat, z=NULL, batch= NULL,  max.iter=200, beta=1e-6,

   # Initialization

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

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

-  bgDist = rowSums( omat ) / sum( omat) 

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

+  bgDist = rowSums( counts ) / sum( counts) 

   bgDist = matrix( rep( bgDist, nC), ncol=nC) 

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

   ll =c()

-  ll.round = bg.calcLL(omat=omat, cellDist=cellDist, bgDist=bgDist, theta=theta)  

+  ll.round = bg.calcLL(counts=counts, cellDist=cellDist, bgDist=bgDist, theta=theta)  

   # EM updates

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

-    next.decon = cD.calcEMbgDecontamination(omat=omat, cellDist=cellDist, bgDist=bgDist, theta=theta, beta=beta, delta=delta)  

+    next.decon = cD.calcEMbgDecontamination(counts=counts, cellDist=cellDist, bgDist=bgDist, theta=theta, beta=beta, delta=delta)  

     theta = next.decon$theta

     cellDist = next.decon$cellDist

     # Calculate log-likelihood

-    ll.temp = bg.calcLL(omat=omat, cellDist=cellDist, bgDist=bgDist, theta=theta) 

+    ll.temp = bg.calcLL(counts=counts, cellDist=cellDist, bgDist=bgDist, theta=theta) 

     ll  = c(ll, ll.temp) 

     ll.round = c( ll.round, round( ll.temp, 2) )

@@ -301,7 +301,7 @@ DecontXoneBatch = function(omat, z=NULL, batch= NULL,  max.iter=200, beta=1e-6,

   }

-  res.conp  = 1- colSums(next.decon$est.rmat) / colSums(omat)

+  res.conp  = 1- colSums(next.decon$est.rmat) / colSums(counts)

   end.time = Sys.time() 

@@ -334,9 +334,9 @@ checkParameters.decon = function(proportion.prior, distribution.prior) {

 }

 # Make sure provided rmat is the right type

-checkCounts.decon = function(omat) {

-  if ( sum(is.na(omat)) >0   ) {

-    stop("Missing value in 'omat' matrix.") 

+checkCounts.decon = function(counts) {

+  if ( sum(is.na(counts)) >0   ) {

+    stop("Missing value in 'counts' matrix.") 

   }

 }

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

 Highly expressed genes from various cells clusters will be expressed at low levels in other clusters in droplet-based systems due to contamination. DecontX will decompose an observed count matrix into a decontaminated expression matrix and a contamination matrix. The only other parameter needed is a vector of cell cluster labels. 

 ```

-new.counts = DecontX( omat = counts, z = cell.label) 

+new.counts = DecontX( counts = counts, z = cell.label) 

 # Decontaminated matrix: new.counts$res.list$est.rmat

 # Percentage of contamination per cell: new.counts$res.list$est.conp

 ```

@@ -4,11 +4,12 @@

 \alias{DecontX}

 \title{This function updates decontamination on dataset with multiple batches}

 \usage{

-DecontX(omat, z = NULL, batch = NULL, max.iter = 200, beta = 1e-06,

-  delta = 10, logfile = NULL, verbose = TRUE, seed = 1234567)

+DecontX(counts, z = NULL, batch = NULL, max.iter = 200,

+  beta = 1e-06, delta = 10, logfile = NULL, verbose = TRUE,

+  seed = 1234567)

 }

 \arguments{

-\item{omat}{Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells.}

+\item{counts}{Numeric/Integer matrix. Observed count matrix, rows represent features and columns represent cells.}

 \item{z}{Integer vector. Cell population labels. Default NULL}

@@ -30,6 +31,6 @@ DecontX(omat, z = NULL, batch = NULL, max.iter = 200, beta = 1e-06,

 This function updates decontamination on dataset with multiple batches

 }

 \examples{

-decon.c = DecontX( omat = contamination.sim$rmat + contamination.sim$cmat, z=contamination.sim$z, max.iter=3)

-decon.bg = DecontX( omat=contamination.sim$rmat + contamination.sim$cmat, max.iter=3 ) 

+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 ) 

 }

@@ -41,16 +41,16 @@

   test_that( desc = "Testing DecontXoneBatch", {

     expect_equal( model_DecontXoneBatch$res.list$est.conp  , 1 - colSums(model_DecontXoneBatch$res.list$est.rmat) /  colSums( Decon.sim$rmat + Decon.sim$cmat) )

     expect_equal( model_DecontXoneBatch$res.list$theta, (model_DecontXoneBatch$run.params$delta + colSums(model_DecontXoneBatch$res.list$est.rmat) ) / ( 2*model_DecontXoneBatch$run.params$delta + colSums(Decon.sim$cmat + Decon.sim$rmat) ) )

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, beta=-1), "'beta' should be a single positive value.")

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, beta=c(1,1) ), "'beta' should be a single positive value.")

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, delta=-1), "'delta' should be a single positive value.")

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, delta=c(1,1) ), "'delta' should be a single positive value.")

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=c(Decon.sim$z, 1) ), "'z' must be of the same length as the number of cells in the 'counts' matrix.")   

-    expect_error( DecontXoneBatch(omat=Decon.sim$rmat+Decon.sim$cmat, z=rep(1, ncol(Decon.sim$rmat)) ), "'z' must have at least 2 different values.") 

-

-    omat.NA = Decon.sim$rmat + Decon.sim$cmat

-    omat.NA[1,1] = NA

-    expect_error( DecontXoneBatch(omat=omat.NA, z=Decon.sim$z), "Missing value in 'omat' matrix.") 

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, beta=-1), "'beta' should be a single positive value.")

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, beta=c(1,1) ), "'beta' should be a single positive value.")

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, delta=-1), "'delta' should be a single positive value.")

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=Decon.sim$z, delta=c(1,1) ), "'delta' should be a single positive value.")

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=c(Decon.sim$z, 1) ), "'z' must be of the same length as the number of cells in the 'counts' matrix.")   

+    expect_error( DecontXoneBatch(counts=Decon.sim$rmat+Decon.sim$cmat, z=rep(1, ncol(Decon.sim$rmat)) ), "'z' must have at least 2 different values.") 

+

+    counts.NA = Decon.sim$rmat + Decon.sim$cmat

+    counts.NA[1,1] = NA

+    expect_error( DecontXoneBatch(counts=counts.NA, z=Decon.sim$z), "Missing value in 'counts' matrix.") 

     } )

   test_that( desc = " Testing DecontXoneBatch using background distribution", {

@@ -61,18 +61,18 @@

   # logLikelihood

   test_that( desc = "Testing logLikelihood.DecontXoneBatch", {

     z.process = processCellLabels(Decon.sim$z, num.cells=ncol(Decon.sim$rmat) )

-    expect_equal( decon.calcLL(omat=Decon.sim$cmat+Decon.sim$rmat, z=z.process  ,  theta=model_DecontXoneBatch$res.list$theta, eta=model_DecontXoneBatch$res.list$est.ConDist, phi=model_DecontXoneBatch$res.list$est.GeneDist ), model_DecontXoneBatch$res.list$logLikelihood[ model_DecontXoneBatch$run.params$iteration  ] )

+    expect_equal( decon.calcLL(counts=Decon.sim$cmat+Decon.sim$rmat, z=z.process  ,  theta=model_DecontXoneBatch$res.list$theta, eta=model_DecontXoneBatch$res.list$est.ConDist, phi=model_DecontXoneBatch$res.list$est.GeneDist ), model_DecontXoneBatch$res.list$logLikelihood[ model_DecontXoneBatch$run.params$iteration  ] )

     cellDist.model.bg = normalizeCounts( model_DecontXoneBatchbg$res.list$est.rmat, normalize="proportion", pseudocount.normalize= model_DecontXoneBatchbg$run.params$beta) 

     bgDist.model.bg = rowSums( Decon.sim$rmat+ Decon.sim$cmat) / sum( Decon.sim$N.by.C) 

     bgDist.model.bg = matrix( rep(bgDist.model.bg, length(Decon.sim$N.by.C)   ), ncol= length(Decon.sim$N.by.C)  )

-    expect_equal( bg.calcLL( omat=Decon.sim$cmat+Decon.sim$rmat, theta=model_DecontXoneBatchbg$res.list$theta, cellDist= cellDist.model.bg, bgDist= bgDist.model.bg), model_DecontXoneBatchbg$res.list$logLikelihood[ model_DecontXoneBatchbg$run.params$iteration  ] )

+    expect_equal( bg.calcLL( counts=Decon.sim$cmat+Decon.sim$rmat, theta=model_DecontXoneBatchbg$res.list$theta, cellDist= cellDist.model.bg, bgDist= bgDist.model.bg), model_DecontXoneBatchbg$res.list$logLikelihood[ model_DecontXoneBatchbg$run.params$iteration  ] )

     } )

   # decontamination EM updates

   test_that( desc = "Testing decontamination EM updates", {

     z.process = processCellLabels(Decon.sim$z, num.cells=ncol(Decon.sim$rmat) )

-    expect_equal( cD.calcEMDecontamination( omat=Decon.sim$cmat+Decon.sim$rmat, z=z.process, K=length(unique(Decon.sim$z)), theta=model_DecontXoneBatch.iter1$res.list$theta, phi=model_DecontXoneBatch.iter1$res.list$est.GeneDist, eta=model_DecontXoneBatch.iter1$res.list$est.ConDist, beta=model_DecontXoneBatch.iter1$run.params$beta, delta=model_DecontXoneBatch.iter1$run.params$delta)$theta,   model_DecontXoneBatch$res.list$theta )

+    expect_equal( cD.calcEMDecontamination( counts=Decon.sim$cmat+Decon.sim$rmat, z=z.process, K=length(unique(Decon.sim$z)), theta=model_DecontXoneBatch.iter1$res.list$theta, phi=model_DecontXoneBatch.iter1$res.list$est.GeneDist, eta=model_DecontXoneBatch.iter1$res.list$est.ConDist, beta=model_DecontXoneBatch.iter1$run.params$beta, delta=model_DecontXoneBatch.iter1$run.params$delta)$theta,   model_DecontXoneBatch$res.list$theta )

     } )