fixed bug which occurs at random probability initiation for a low num…

…ber of cells

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mnem
 Type: Package
 Title: Mixture Nested Effects Models
-Version: 0.99.9
+Version: 0.99.10
 Author: Martin Pirkl
 Maintainer: Martin Pirkl <martin.pirkl@bsse.ethz.ch>
 Description: Mixture Nested Effects Models (mnem) is an extension of Nested Effects Models and allows for the analysis of single cell perturbation data provided by methods like Perturb-Seq (Dixit et al., 2016) or Crop-Seq (Datlinger et al., 2017). In those experiments each of many cells is perturbed by a knock-down of a specific gene, i.e. several cells are perturbed by a knock-down of gene A, several by a knock-down of gene B, ... and so forth. The observed read-out has to be multi-trait and in the case of the Perturb-/Crop-Seq gene are expression profiles for each cell. mnem uses a mixture model to simultaneously cluster the cell population into k clusters and and infer k networks causally linking the perturbed genes for each cluster. The mixture components are inferred via an expectation maximization algorithm.

R/mnems.r

@@ -412,10 +412,7 @@ mnem <- function(D, inference = "em", search = "modules", start = NULL,
                     if (length(probscl) >= s & type %in% "cluster") {
                         probs <- probscl[[s]]
                     } else {
-                        probs <- matrix(log2(sample(c(0,1), k*ncol(data),
-                                                    replace = TRUE,
-                                                    prob = c(0.9, 0.1))), k,
-                                        ncol(data))
+                        probs <- random_probs(k, data)
                     }
                     mw <- apply(getAffinity(probs, affinity = affinity,
                                             norm = TRUE, logtype = logtype,

R/mnems_low.r

@@ -1,4 +1,66 @@
 #' @noRd
+random_probs <- function(k, data, full = FALSE) {
+    probs <- matrix(log2(sample(c(0,1), k*ncol(data),
+                                replace = TRUE,
+                                prob = c(0.9, 0.1))), k,
+                    ncol(data))
+    if (full) {
+        for (i in seq_len(k)) {
+            if (i == 1) { next() }
+            infcells <- which(apply(probs, 2, function(x) {
+                bad <- FALSE
+                if (all(is.infinite(x))) {
+                    bad <- TRUE
+                }
+                return(bad)
+            }))
+            if (i == k) {
+                probs[i, infcells] <- log2(1)
+            } else {
+                probs[i, infcells] <-
+                    log2(sample(c(0,1),
+                                length(infcells),
+                                replace = TRUE,
+                                prob = c(1-1/k, 1/k)))
+            }
+        }
+    }
+    while(any(apply(probs, 1, function(x) {
+        bad <- FALSE
+        if (all(is.infinite(x)) | all(x == 0)) {
+            bad <- TRUE
+        }
+        return(bad)
+    }))) {
+        probs <- matrix(log2(sample(c(0,1), k*ncol(data),
+                                    replace = TRUE,
+                                    prob = c(0.9, 0.1))), k,
+                        ncol(data))
+        if (full) {
+            for (i in seq_len(k)) {
+                if (i == 1) { next() }
+                infcells <- which(apply(probs, 2, function(x) {
+                    bad <- FALSE
+                    if (all(is.infinite(x))) {
+                        bad <- TRUE
+                    }
+                    return(bad)
+                }))
+                if (i == k) {
+                    probs[i, infcells] <- log2(1)
+                } else {
+                    probs[i, infcells] <-
+                        log2(sample(c(0,1),
+                                    length(infcells),
+                                    replace = TRUE,
+                                    prob = c(1-1/k, 1/k)))
+                }
+            }
+        }
+    }
+    return(probs)
+}
+#' @noRd
 sortAdj <- function(res, list = FALSE) {
     resmat <- NULL
     for (i in seq_len(length(res))) {