@@ -6,12 +6,27 @@

 #' across 2 experimental conditions from a real scRNA-seq data set.

 #' 

 #' @param x a \code{\link[SingleCellExperiment]{SingleCellExperiment}}.

-#' @param nc,ns,nk # of cells, samples and clusters to simulate. 

-#'   By default, \code{ns = NULL} will simulated as many samples as 

-#'   available in the reference to avoid duplicated reference samples.

+#' @param ng number of genes to simulate. Importantly, for the library sizes 

+#'   computed by \code{\link{prepSim}} (= \code{exp(x$offset)}) to make sense, 

+#'   the number of simulated genes should match with the number of genes 

+#'   in the reference. To simulate a reduced number of genes, e.g. for 

+#'   testing and development purposes, please set \code{force = TRUE}.

+#' @param nc number of cells to simulate.

+#' @param nk number of clusters to simulate; defaults to the number

+#'   of available reference clusters (\code{nlevels(x$cluster_id)}).  

+#' @param ns number of samples to simulate; defaults to as many as

+#'   available in the reference to avoid duplicated reference samples. 

+#'   Specifically, the number of samples will be set to 

+#'   \code{n = nlevels(x$sample_id)} when \code{dd = FALSE}, 

+#'   \code{n} per group  when \code{dd, paired = TRUE}, and 

+#'   \code{floor(n/2)} per group when \code{dd = TRUE, paired = FALSE}.

+#'   When a larger number samples should be simulated, set \code{force = TRUE}.

 #' @param probs a list of length 3 containing probabilities of a cell belonging

 #'   to each cluster, sample, and group, respectively. List elements must be 

 #'   NULL (equal probabilities) or numeric values in [0, 1] that sum to 1.

+#' @param dd whether or not to simulate differential distributions; if TRUE, 

+#'   two groups are simulated and \code{ns} corresponds to the number of 

+#'   samples per group, else one group with \code{ns} samples is simulated.

 #' @param p_dd numeric vector of length 6.

 #'   Specifies the probability of a gene being

 #'   EE, EP, DE, DP, DM, or DB, respectively.

@@ -39,17 +54,12 @@

 #'   of the branches lengths separating them. 

 #' @param phylo_pars vector of length 2 providing the parameters that control 

 #'   the number of type genes. Passed to an exponential PDF (see details).

-#'   

-#' @param ng # of genes to simulate. Importantly, for the library sizes 

-#'   computed by \code{\link{prepSim}} (= \code{exp(x$offset)}) to make sense, 

-#'   the number of simulated genes should match with the number of genes 

-#'   in the reference. To simulate a reduced number of genes, e.g. for 

-#'   testing and development purposes, please set \code{force = TRUE}.

-#' @param force logical specifying whether to force 

-#'   simulation despite \code{ng != nrow(x)}.

+#' @param force logical specifying whether to force simulation 

+#'   when \code{ng} and/or \code{ns} don't match the number of 

+#'   available reference genes and samples, respectively.

 #'   

 #' @details The simulation of type genes can be performed in 2 ways; 

-#'   (1) via \code{p_type} to simulate independant clusters, OR 

+#'   (1) via \code{p_type} to simulate independent clusters, OR 

 #'   (2) via \code{phylo_tree} to simulate a hierarchical cluster structure.

 #'   

 #'   For (1), a subset of \code{p_type} \% of genes are selected per cluster

@@ -89,11 +99,11 @@

 #'   \item{\code{args}}{a list of the function call's input arguments.}}}}

 #'   

 #' @examples

-#' data(sce)

+#' data(example_sce)

 #' library(SingleCellExperiment)

 #' 

 #' # prep. SCE for simulation

-#' ref <- prepSim(sce)

+#' ref <- prepSim(example_sce)

 #' 

 #' # simulate data

 #' (sim <- simData(ref, nc = 200,

@@ -179,11 +189,17 @@ simData <- function(x,

     if (is.null(x$cluster_id)) {

         x$cluster_id <- factor("foo")

         no_k <- TRUE

-    } else no_k <- FALSE

+    } else {

+        x$cluster_id <- droplevels(factor(x$cluster_id))

+        no_k <- nlevels(x$cluster_id) == 1

+    }

     if (is.null(x$sample_id)) {

         x$sample_id <- factor("foo")

         no_s <- TRUE

-    } else no_s <- FALSE

+    } else {

+        x$sample_id <- droplevels(factor(x$sample_id))

+        no_s <- nlevels(x$sample_id) == 1

+    }

     # store all input arguments to be returned in final output

     args <- c(as.list(environment()))

@@ -192,12 +208,14 @@ simData <- function(x,

     .check_sce(x, req_group = FALSE)

     args_tmp <- .check_args_simData(as.list(environment()))

     nk <- args$nk <- args_tmp$nk

-    ns <- args$ns <- args_tmp$ns

-    

+

     # reference IDs

     nk0 <- length(kids0 <- set_names(levels(x$cluster_id)))

     ns0 <- length(sids0 <- set_names(levels(x$sample_id)))

+    # get number of samples to simulate

+    ns <- .get_ns(ns0, ns, dd, paired, force)

+    

     # simulation IDs

     nk <- length(kids <- set_names(paste0("cluster", seq_len(nk))))

     sids <- set_names(paste0("sample", seq_len(ns)))

@@ -211,8 +229,13 @@ simData <- function(x,

         ref_sids <- cbind(ref_sids, ref_sids)

     } else {

         # draw reference samples at random for each group

-        sidsA <- sample(sids0, ns, ns > ns0)

-        sidsB <- sample(setdiff(sids0, sidsA), ns, ns > ns0)

+        sidsA <- sample(sids0, ns, force && ns > ns0)

+        if (force) {

+            sidsB <- sample(sids0, ns, ns > ns0)

+        } else {

+            sidsB <- setdiff(sids0, sidsA)

+            sidsB <- sample(sidsB, ns)

+        }

         ref_sids <- cbind(sidsA, sidsB)

     }

     dimnames(ref_sids) <- list(sids, gids)

@@ -299,6 +299,7 @@ simData <- function(x,

         bs$cluster_id <- cbind(0, bs$cluster_id)

         names(bs$cluster_id) <- kids0

     }

+    os <- x$offset

     b0 <- bs$beta0

     ds <- rowData(x)$disp

@@ -337,8 +338,8 @@ simData <- function(x,

                 ds_kc <- ds[gs0]

                 lfc_kc <- lfc[[c, k]]

                 bs_ksc <- exp(rowSums(bs_ks[gs0, , drop = FALSE]))

-                ms_g1 <- outer(bs_ksc, exp(x$offset[cs_g1]), "*")

-                ms_g2 <- outer(bs_ksc, exp(x$offset[cs_g2]), "*")

+                ms_g1 <- outer(bs_ksc, exp(os[cs_g1]), "*")

+                ms_g2 <- outer(bs_ksc, exp(os[cs_g2]), "*")

                 re <- .sim(c, cs_g1, cs_g2, ms_g1, ms_g2, ds_kc, lfc_kc, p_ep, p_dp, p_dm)

                 y[gi, unlist(ci)] <- re$cs

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

 #' 

 #' @param x a \code{\link[SingleCellExperiment]{SingleCellExperiment}}.

 #' @param nc,ns,nk # of cells, samples and clusters to simulate. 

+#'   By default, \code{ns = NULL} will simulated as many samples as 

+#'   available in the reference to avoid duplicated reference samples.

 #' @param probs a list of length 3 containing probabilities of a cell belonging

 #'   to each cluster, sample, and group, respectively. List elements must be 

 #'   NULL (equal probabilities) or numeric values in [0, 1] that sum to 1.

@@ -157,12 +159,14 @@

 #' @importFrom S4Vectors split unfactor

 #' @export

-simData <- function(x, nc = 2e3, ns = 3, nk = 3,

-    probs = NULL, p_dd = diag(6)[1, ], paired = FALSE,

+simData <- function(x, 

+    ng = nrow(x), nc = ncol(x), 

+    ns = NULL, nk = NULL, probs = NULL, 

+    dd = TRUE, p_dd = diag(6)[1, ], paired = FALSE,

     p_ep = 0.5, p_dp = 0.3, p_dm = 0.5,

     p_type = 0, lfc = 2, rel_lfc = NULL, 

     phylo_tree = NULL, phylo_pars = c(ifelse(is.null(phylo_tree), 0, 0.1), 3),

-    ng = nrow(x), force = FALSE) {

+    force = FALSE) {

     # throughout this code...

     # k: cluster ID

@@ -171,6 +175,16 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # c: DD category

     # 0: reference

+    # add mock cluster/sample ID if missing

+    if (is.null(x$cluster_id)) {

+        x$cluster_id <- factor("foo")

+        no_k <- TRUE

+    } else no_k <- FALSE

+    if (is.null(x$sample_id)) {

+        x$sample_id <- factor("foo")

+        no_s <- TRUE

+    } else no_s <- FALSE

+    

     # store all input arguments to be returned in final output

     args <- c(as.list(environment()))

@@ -178,6 +192,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     .check_sce(x, req_group = FALSE)

     args_tmp <- .check_args_simData(as.list(environment()))

     nk <- args$nk <- args_tmp$nk

+    ns <- args$ns <- args_tmp$ns

     # reference IDs

     nk0 <- length(kids0 <- set_names(levels(x$cluster_id)))

@@ -190,17 +205,18 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # sample reference clusters & samples

     ref_kids <- setNames(sample(kids0, nk, nk > nk0), kids)

-    if (paired) { 

+    if (!dd || paired) {

         # use same set of reference samples for both groups

         ref_sids <- sample(sids0, ns, ns > ns0)

-        ref_sids <- replicate(length(gids), ref_sids)

+        ref_sids <- cbind(ref_sids, ref_sids)

     } else {

         # draw reference samples at random for each group

-        ref_sids <- replicate(length(gids), 

-            sample(sids0, ns, ns > ns0))

+        sidsA <- sample(sids0, ns, ns > ns0)

+        sidsB <- sample(setdiff(sids0, sidsA), ns, ns > ns0)

+        ref_sids <- cbind(sidsA, sidsB)

     }

     dimnames(ref_sids) <- list(sids, gids)

-    

+

     if (is.null(rel_lfc)) 

         rel_lfc <- rep(1, nk)

     if (is.null(names(rel_lfc))) {

@@ -274,14 +290,17 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

         }), vector("list", length(cats)))

     # compute NB parameters

-    m <- lapply(sids0, function(s) {

-        b <- paste0("beta.", s)

-        b <- exp(rowData(x)[[b]])

-        m <- outer(b, exp(x$offset), "*")

-        dimnames(m) <- dimnames(x); m

-    })

-    d <- rowData(x)$dispersion 

-    names(d) <- rownames(x)

+    bs <- rowData(x)$beta

+    if (!is.null(bs$sample_id)) {

+        bs$sample_id <- cbind(0, bs$sample_id)

+        names(bs$sample_id) <- sids0

+    }

+    if (!is.null(bs$cluster_id)) {

+        bs$cluster_id <- cbind(0, bs$cluster_id)

+        names(bs$cluster_id) <- kids0

+    }

+    b0 <- bs$beta0

+    ds <- rowData(x)$disp

     # initialize list of depth two to store 

     # simulation means in each cluster & group

@@ -292,13 +311,18 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # run simulation -----------------------------------------------------------

     for (k in kids) {

         for (s in sids) {

+            # get output cell indices

+            ci <- cs_idx[[k]][[s]]

+            

             # get reference samples, clusters & cells

             s0 <- ref_sids[s, ]

             k0 <- ref_kids[k]

             cs0 <- cs_by_ks[[k0]][s0]

-            # get output cell indices

-            ci <- cs_idx[[k]][[s]]

+            # get NB parameters

+            bs_ks <- cbind(b0,

+                bs$cluster_id[[k0]],

+                bs$sample_id[[s0[1]]])

             for (c in cats[n_dd[, k] != 0]) {

                 # sample cells to simulate from

@@ -310,12 +334,13 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

                 gi <- gs_idx[[c, k]]

                 # get NB parameters

-                m_g1 <- m[[s0[[1]]]][gs0, cs_g1, drop = FALSE]

-                m_g2 <- m[[s0[[2]]]][gs0, cs_g2, drop = FALSE]

-                d_kc <- d[gs0]

+                ds_kc <- ds[gs0]

                 lfc_kc <- lfc[[c, k]]

+                bs_ksc <- exp(rowSums(bs_ks[gs0, , drop = FALSE]))

+                ms_g1 <- outer(bs_ksc, exp(x$offset[cs_g1]), "*")

+                ms_g2 <- outer(bs_ksc, exp(x$offset[cs_g2]), "*")

-                re <- .sim(c, cs_g1, cs_g2, m_g1, m_g2, d_kc, lfc_kc, p_ep, p_dp, p_dm)

+                re <- .sim(c, cs_g1, cs_g2, ms_g1, ms_g2, ds_kc, lfc_kc, p_ep, p_dp, p_dm)

                 y[gi, unlist(ci)] <- re$cs

                 for (g in gids) sim_mean[[k]][[g]][gi] <- ifelse(

@@ -331,7 +356,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

         category = rep.int(rep(cats, nk), c(n_dd)),

         logFC = unlist(lfc),

         sim_gene = unlist(gs_by_kc),

-        sim_disp = d[unlist(gs_by_kc)]) %>% 

+        sim_disp = ds[unlist(gs_by_kc)]) %>% 

         mutate_at("gene", as.character)

     # add true simulation means

     sim_mean <- sim_mean %>%

@@ -348,14 +373,22 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # construct SCE ------------------------------------------------------------

     # cell metadata including group, sample, cluster IDs

     cd$group_id <- droplevels(cd$group_id)

-    cd$sample_id <- factor(paste(cd$sample_id, cd$group_id, sep = "."))

+    cd$sample_id <- if (dd) {

+        factor(paste(cd$sample_id, cd$group_id, sep = "."))

+    } else factor(cd$sample_id)

     m <- match(levels(cd$sample_id), cd$sample_id)

     gids <- cd$group_id[m]

     o <- order(gids)

     sids <- levels(cd$sample_id)[o]

     cd <- cd %>% 

-        mutate_at("cluster_id", factor, levels = kids) %>% 

-        mutate_at("sample_id", factor, levels = sids) 

+        mutate_at("sample_id", factor, levels = sids) %>% 

+        mutate_at("cluster_id", factor, levels = kids)

+    if (!dd) {

+        cd$group_id <- NULL

+        ref_sids <- ref_sids[, 1]

+    }

+    if (no_s) cd$sample_id <- NULL

+    if (no_k) cd$cluster_id <- NULL

     # gene metadata storing gene classes & specificities

     rd <- DataFrame(class = factor(class, 

         levels = c("state", "shared", "type")))

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

 #' @param p_dd numeric vector of length 6.

 #'   Specifies the probability of a gene being

 #'   EE, EP, DE, DP, DM, or DB, respectively.

-#' @param paired logial specifying whether a paired design should 

+#' @param paired logical specifying whether a paired design should 

 #'   be simulated (both groups use the same set of reference samples) 

 #'   or not (reference samples are drawn at random).

 #' @param p_ep,p_dp,p_dm numeric specifying the proportion of cells

@@ -21,8 +21,8 @@

 #' @param p_type numeric. Probability of EE/EP gene being a type-gene.

 #'   If a gene is of class "type" in a given cluster, a unique mean 

 #'   will be used for that gene in the respective cluster.

-#' @param lfc numeric value to use as mean logFC

-#'   for DE, DP, DM, and DB type of genes.

+#' @param lfc numeric value to use as mean logFC 

+#'   (logarithm base 2) for DE, DP, DM, and DB type of genes.

 #' @param rel_lfc numeric vector of relative logFCs for each cluster. 

 #'   Should be of length \code{nlevels(x$cluster_id)} with 

 #'   \code{levels(x$cluster_id)} as names. 

@@ -32,22 +32,11 @@

 #'   \href{http://evolution.genetics.washington.edu/phylip/newicktree.html}{here}. 

 #'   The distance between the nodes, except for the original branch, will be 

 #'   translated in the number of shared genes between the clusters belonging to 

-#'   these nodes (this relation is controlled with \code{phylo_pars}). The distance 

-#'   between two clusters is defined as the sum of the branches lengths 

-#'   separating them. 

+#'   these nodes (this relation is controlled with \code{phylo_pars}). 

+#'   The distance between two clusters is defined as the sum 

+#'   of the branches lengths separating them. 

 #' @param phylo_pars vector of length 2 providing the parameters that control 

-#'   the number of type genes. Passed to an exponential's PDF:

-#'   \code{N = #genes x gamma1 * e^(-gamma2 x dist)},

-#'   

-#'   \itemize{

-#'   \item  \code{gamma1} is the parameter that controls the percentage of shared 

-#'   genes between the nodes. By default 0.1 if a tree is given, meaning that 

-#'   maximum 10% of the genes can be used as type genes (if gamma2 = 0). 

-#'   However it's advised to tune it depending on the input \code{prep_sce}. 

-#'   \code{gamma2} is the 'penalty' of increasing the distance between clusters 

-#'   (\code{dist}, defined by \code{phylo_tree}), applied on the number of 

-#'   shared genes. Default to 3. 

-#'   }

+#'   the number of type genes. Passed to an exponential PDF (see details).

 #'   

 #' @param ng # of genes to simulate. Importantly, for the library sizes 

 #'   computed by \code{\link{prepSim}} (= \code{exp(x$offset)}) to make sense, 

@@ -58,13 +47,45 @@

 #'   simulation despite \code{ng != nrow(x)}.

 #'   

 #' @details The simulation of type genes can be performed in 2 ways; 

-#'   (1) by defining \code{p_type} and thus simulating independant clusters, OR 

-#'   (2) by defining both \code{phylo_tree} and \code{phylo_pars}, which will 

-#'   simulate a hierarchical structure between the clusters.

+#'   (1) via \code{p_type} to simulate independant clusters, OR 

+#'   (2) via \code{phylo_tree} to simulate a hierarchical cluster structure.

+#'   

+#'   For (1), a subset of \code{p_type} \% of genes are selected per cluster

+#'   to use a different references genes than the remainder of clusters,

+#'   giving rise to cluster-specific NB means for count sampling.

+#'   

+#'   For (2), the number of shared/type genes at each node 

+#'   are given by \code{a*G*e^(-b*d)}, where \itemize{

+#'   \item{\code{a} -- controls the percentage of shared genes between nodes.

+#'   By default, at most 10\% of the genes are reserved as type genes

+#'   (when \code{b} = 0). However, it is advised to tune this parameter 

+#'   depending on the input \code{prep_sce}.}

+#'   \item{\code{b} -- determines how the number of shared genes 

+#'   decreases with increasing distance d between clusters 

+#'   (defined through \code{phylo_tree}).}}

 #'  

 #' @return a \code{\link[SingleCellExperiment]{SingleCellExperiment}}

-#'   containing multiple clusters & samples across 2 groups.

-#' 

+#'   containing multiple clusters & samples across 2 groups 

+#'   as well as the following metadata: \describe{

+#'   \item{cell metadata (\code{colData(.)})}{a \code{DataFrame} containing,

+#'   containing, for each cell, it's cluster, sample, and group ID.}

+#'   \item{gene metadata (\code{rowData(.)})}{a \code{DataFrame} containing, 

+#'   for each gene, it's \code{class} (one of "state", "type", "none") and 

+#'   specificity (\code{specs}; NA for genes of type "state", otherwise 

+#'   a character vector of clusters that share the given gene).}

+#'   \item{experiment metadata (\code{metadata(.)})}{

+#'   \describe{

+#'   \item{\code{experiment_info}}{a \code{data.frame} 

+#'   summarizing the experimental design.}

+#'   \item{\code{n_cells}}{the number of cells for each sample.}

+#'   \item{\code{gene_info}}{a \code{data.frame} containing, for each gene

+#'   in each cluster, it's differential distribution \code{category},

+#'   mean \code{logFC} (NA for genes for categories "ee" and "ep"),

+#'   gene used as reference (\code{sim_gene}), dispersion \code{sim_disp},

+#'   and simulation means for each group \code{sim_mean.A/B}.}

+#'   \item{\code{ref_sids/kidskids}}{the sample/cluster IDs used as reference.}

+#'   \item{\code{args}}{a list of the function call's input arguments.}}}}

+#'   

 #' @examples

 #' data(sce)

 #' library(SingleCellExperiment)

@@ -73,7 +94,7 @@

 #' ref <- prepSim(sce)

 #' 

 #' # simulate data

-#' (sim <- simData(ref, nc = 10,

+#' (sim <- simData(ref, nc = 200,

 #'   p_dd = c(0.9, 0, 0.1, 0, 0, 0),

 #'   ng = 100, force = TRUE,

 #'   probs = list(NULL, NULL, c(1, 0))))

@@ -116,7 +137,7 @@

 #' # view information about shared 'type' genes

 #' table(rowData(sim)$class)

 #' 

-#' @author Helena L Crowell

+#' @author Helena L Crowell & Anthony Sonrel

 #' 

 #' @references 

 #' Crowell, HL, Soneson, C, Germain, P-L, Calini, D, 

@@ -316,10 +337,10 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     sim_mean <- sim_mean %>%

         map(bind_cols) %>% 

         bind_rows(.id = "cluster_id") %>% 

-        mutate_at("cluster_id", factor) %>% 

         mutate(gene = rep(gs, nk))

     gi <- full_join(gi, sim_mean, by = c("gene", "cluster_id")) %>% 

-        rename("sim_mean.A" = "A", "sim_mean.B" = "B")

+        rename("sim_mean.A" = "A", "sim_mean.B" = "B") %>% 

+        mutate_at("cluster_id", factor)

     # reorder

     o <- order(as.numeric(gsub("[a-z]", "", gi$gene)))

     gi <- gi[o, ]; rownames(gi) <- NULL

@@ -332,14 +353,16 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     gids <- cd$group_id[m]

     o <- order(gids)

     sids <- levels(cd$sample_id)[o]

-    ei <- data.frame(sample_id = sids, group_id = gids[o])

-    cd <- cd %>% mutate_at("sample_id", factor, levels = sids)

+    cd <- cd %>% 

+        mutate_at("cluster_id", factor, levels = kids) %>% 

+        mutate_at("sample_id", factor, levels = sids) 

     # gene metadata storing gene classes & specificities

     rd <- DataFrame(class = factor(class, 

         levels = c("state", "shared", "type")))

     rd$specs <- as.list(specs)

     # simulation metadata including used reference samples/cluster, 

     # list of input arguments, and simulated genes' metadata

+    ei <- data.frame(sample_id = sids, group_id = gids[o])

     md <- list(

         experiment_info = ei,

         n_cells = table(cd$sample_id),

@@ -109,7 +109,8 @@

 #' plot(read.dendrogram(text = phylo_tree))

 #' 

 #' # simulate clusters accordingly

-#' sim <- simData(ref, phylo_tree = phylo_tree, 

+#' sim <- simData(ref, 

+#'   phylo_tree = phylo_tree, 

 #'   phylo_pars = c(0.1, 3), 

 #'   ng = 500, force = TRUE)

 #' # view information about shared 'type' genes

@@ -35,11 +35,9 @@

 #'   these nodes (this relation is controlled with \code{phylo_pars}). The distance 

 #'   between two clusters is defined as the sum of the branches lengths 

 #'   separating them. 

-#' @param phylo_pars vector of length 2, defining the parameters that control the 

-#'   number of type genes. It is passed to an adaptation of the 

-#'   exponential's PDF:

-#'   

-#'   \code{N = Ngenes x gamma1 * e^(-gamma2 x dist)} ,

+#' @param phylo_pars vector of length 2 providing the parameters that control 

+#'   the number of type genes. Passed to an exponential's PDF:

+#'   \code{N = #genes x gamma1 * e^(-gamma2 x dist)},

 #'   

 #'   \itemize{

 #'   \item  \code{gamma1} is the parameter that controls the percentage of shared 

@@ -59,11 +57,10 @@

 #' @param force logical specifying whether to force 

 #'   simulation despite \code{ng != nrow(x)}.

 #'   

-#' @details 

-#'  The simulation of type genes can be performed in 2 ways; (1) by defining 

-#'  \code{p_type} and thus simulating independant clusters, OR (2) by defining both 

-#'  \code{phylo_tree} and \code{phylo_pars}, which will simulate a hierarchical structure 

-#'  between the clusters. Only one of the options is allowed. 

+#' @details The simulation of type genes can be performed in 2 ways; 

+#'   (1) by defining \code{p_type} and thus simulating independant clusters, OR 

+#'   (2) by defining both \code{phylo_tree} and \code{phylo_pars}, which will 

+#'   simulate a hierarchical structure between the clusters.

 #'  

 #' @return a \code{\link[SingleCellExperiment]{SingleCellExperiment}}

 #'   containing multiple clusters & samples across 2 groups.

@@ -88,7 +85,7 @@

 #' table(gi$category)

 #' 

 #' # unbalanced sample sizes

-#' sim <- simData(ref, nc = 100,

+#' sim <- simData(ref, nc = 100, ns = 2,

 #'   probs = list(NULL, c(0.25, 0.75), NULL),

 #'   ng = 10, force = TRUE)

 #' table(sim$sample_id)

@@ -142,7 +139,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     probs = NULL, p_dd = diag(6)[1, ], paired = FALSE,

     p_ep = 0.5, p_dp = 0.3, p_dm = 0.5,

     p_type = 0, lfc = 2, rel_lfc = NULL, 

-    phylo_tree = NULL, phylo_pars = c(0, 3),

+    phylo_tree = NULL, phylo_pars = c(ifelse(is.null(phylo_tree), 0, 0.1), 3),

     ng = nrow(x), force = FALSE) {

     # throughout this code...

@@ -152,9 +149,6 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # c: DD category

     # 0: reference

-    # default shared p if phylo tree given

-    if (missing(phylo_pars)) phylo_pars[1] <- ifelse(is.null(phylo_tree), 0, 0.1)

-    

     # store all input arguments to be returned in final output

     args <- c(as.list(environment()))

@@ -35,23 +35,20 @@

 #'   these nodes (this relation is controlled with \code{phylo_pars}). The distance 

 #'   between two clusters is defined as the sum of the branches lengths 

 #'   separating them. 

-#' @param phylo_pars list of length 2, defining the parameters that control the 

-#'   number of shared/ specific type-genes; \itemize{

-#'   \item The first element of the list is a numeric vector of length 2. 

-#'   It defines the number of shared genes as an adaptation of the 

+#' @param phylo_pars vector of length 2, defining the parameters that control the 

+#'   number of type genes. It is passed to an adaptation of the 

 #'   exponential's PDF:

 #'   

 #'   \code{N = Ngenes x gamma1 * e^(-gamma2 x dist)} ,

 #'   

-#'   where \code{gamma1} is the parameter that controls the percentage of shared genes

-#'   between the nodes. By default 0.2, but it's advised to tune it depending 

-#'   on the input \code{prep_sce}. 

-#'   \code{gamma2} is the 'penalty' of increasing 

-#'   the distance between clusters (\code{dist}, defined by \code{phylo_tree}),

-#'   applied on the number of shared genes. Default to -3. 

-#'   \item The second element can be a single numeric or a list of length nk. 

-#'   It is an equivalent of \code{p_type} for each leaf (i.e. cluster) that is applied 

-#'   after the determination of 'shared type genes'. 

+#'   \itemize{

+#'   \item  \code{gamma1} is the parameter that controls the percentage of shared 

+#'   genes between the nodes. By default 0.1 if a tree is given, meaning that 

+#'   maximum 10% of the genes can be used as type genes (if gamma2 = 0). 

+#'   However it's advised to tune it depending on the input \code{prep_sce}. 

+#'   \code{gamma2} is the 'penalty' of increasing the distance between clusters 

+#'   (\code{dist}, defined by \code{phylo_tree}), applied on the number of 

+#'   shared genes. Default to 3. 

 #'   }

 #'   

 #' @param ng # of genes to simulate. Importantly, for the library sizes 

@@ -116,7 +113,7 @@

 #' 

 #' # simulate clusters accordingly

 #' sim <- simData(ref, phylo_tree = phylo_tree, 

-#'   phylo_pars = list(c(0.1, 3), 0.1), 

+#'   phylo_pars = c(0.1, 3), 

 #'   ng = 500, force = TRUE)

 #' # view information about shared 'type' genes

 #' table(rowData(sim)$class)

@@ -145,7 +142,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     probs = NULL, p_dd = diag(6)[1, ], paired = FALSE,

     p_ep = 0.5, p_dp = 0.3, p_dm = 0.5,

     p_type = 0, lfc = 2, rel_lfc = NULL, 

-    phylo_tree = NULL, phylo_pars = list(c(0, 3), 0),

+    phylo_tree = NULL, phylo_pars = c(0, 3),

     ng = nrow(x), force = FALSE) {

     # throughout this code...

@@ -155,6 +152,9 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # c: DD category

     # 0: reference

+    # default shared p if phylo tree given

+    if (missing(phylo_pars)) phylo_pars[1] <- ifelse(is.null(phylo_tree), 0, 0.1)

+    

     # store all input arguments to be returned in final output

     args <- c(as.list(environment()))

@@ -160,40 +160,12 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # check validity of input arguments

     .check_sce(x, req_group = FALSE)

-    .check_args_simData(as.list(environment()))

-    if (!force && ng != nrow(x))

-        stop("Number of simulated genes should match with reference,\n", 

-            "but 'ng != nrow(x)'; please specify 'force = TRUE' if\n", 

-            "simulation should be forced regardlessly (see '?simData').")

-    if (!is.null(phylo_tree) && p_type != 0)

-        stop("Only one of 'p_type' and 'phylo_tree' can be provided.\n", 

-             "Please see the 'Details' section of '?simData'.")

-    if (!length(phylo_pars[[2]]) %in% c(1, nk))

-        stop("The second element of 'phylo_pars' should be correspond\n", 

-             " to the number of clusters ('nk') or of length 1.")

-    if (!is.null(phylo_tree) && phylo_pars[[1]][1] == 0)

-        warning("'phylo_pars[[1]][1]' has been set to 0;\n", 

-                "'phylo_tree' argument will be ignored.")

-    if (!is.null(phylo_tree) && all(phylo_pars[[2]] == 0))

-        warning("'phylo_pars[[2]]' has been set to 0;\n", 

-                "type genes for individual clusters won't be simulated.")

+    args_tmp <- .check_args_simData(as.list(environment()))

+    nk <- args$nk <- args_tmp$nk

     # reference IDs

     nk0 <- length(kids0 <- set_names(levels(x$cluster_id)))

     ns0 <- length(sids0 <- set_names(levels(x$sample_id)))

-

-    # assure number of simulated clusters

-    # matches with specified phylogeny

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

-        kids_phylo <- .get_clusters_from_phylo(phylo_tree)

-        nk_phylo <- length(kids_phylo)

-        ns_phylo <- as.numeric(gsub("[a-z]", "", kids_phylo))

-        if (!all(sort(ns_phylo) == seq_len(nk_phylo)))

-            stop("Some clusters appear to be missing from 'phylo_tree';\n",

-                "please make sure all clusters up to ", 

-                dQuote(kids_phylo[which.max(ns_phylo)]), " are present.")

-        if (nk_phylo != nk) args$nk <- nk <- nk_phylo

-    }

     # simulation IDs

     nk <- length(kids <- set_names(paste0("cluster", seq_len(nk))))

@@ -244,7 +216,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     gs_idx <- .sample_gene_inds(gs, n_dd)

     # for ea. cluster, sample set of genes to simulate from

-    gs_by_k <- setNames(sample(rownames(x), ng, TRUE), gs)

+    gs_by_k <- setNames(sample(rownames(x), ng, ng > nrow(x)), gs)

     gs_by_k <- replicate(nk, gs_by_k)

     colnames(gs_by_k) <- kids

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

 #'   or not (reference samples are drawn at random).

 #' @param p_ep,p_dp,p_dm numeric specifying the proportion of cells

 #'   to be shifted to a different expression state in one group (see details).

-#' @param p_type numeric. Probaility of EE/EP gene being a type-gene.

+#' @param p_type numeric. Probability of EE/EP gene being a type-gene.

 #'   If a gene is of class "type" in a given cluster, a unique mean 

 #'   will be used for that gene in the respective cluster.

 #' @param lfc numeric value to use as mean logFC

@@ -28,19 +28,32 @@

 #'   \code{levels(x$cluster_id)} as names. 

 #'   Defaults to factor of 1 for all clusters.

 #' @param phylo_tree newick tree text representing cluster relations 

-#'   and their relative distance. If a tree is given, the distance between 

-#'   the clusters will be translated in the number of shared genes 

-#'   (this relation is controlled with \code{phylo_pars}). The distance 

+#'   and their relative distance. An explanation of the syntax can be found 

+#'   \href{http://evolution.genetics.washington.edu/phylip/newicktree.html}{here}. 

+#'   The distance between the nodes, except for the original branch, will be 

+#'   translated in the number of shared genes between the clusters belonging to 

+#'   these nodes (this relation is controlled with \code{phylo_pars}). The distance 

 #'   between two clusters is defined as the sum of the branches lengths 

 #'   separating them. 

-#' @param phylo_pars numeric vector of length 2. Defines the number of shared 

-#'   genes as an adaptation of the exponential's PDF:

-#'   N = Ngenes x gamma1 * e^(-gamma2 x dist) ,

-#'   where gamma1 is the parameter that controls the percentage of shared genes. 

-#'   By default, it corresponds to \code{p_type} but it's advised to tune it 

-#'   depending on the input prep_sce. gamma2 is the 'penalty' of increasing 

-#'   the distance between clusters, applied on the number of shared genes. 

-#'   Default to -3. 

+#' @param phylo_pars list of length 2, defining the parameters that control the 

+#'   number of shared/ specific type-genes; \itemize{

+#'   \item The first element of the list is a numeric vector of length 2. 

+#'   It defines the number of shared genes as an adaptation of the 

+#'   exponential's PDF:

+#'   

+#'   \code{N = Ngenes x gamma1 * e^(-gamma2 x dist)} ,

+#'   

+#'   where \code{gamma1} is the parameter that controls the percentage of shared genes

+#'   between the nodes. By default 0.2, but it's advised to tune it depending 

+#'   on the input \code{prep_sce}. 

+#'   \code{gamma2} is the 'penalty' of increasing 

+#'   the distance between clusters (\code{dist}, defined by \code{phylo_tree}),

+#'   applied on the number of shared genes. Default to -3. 

+#'   \item The second element can be a single numeric or a list of length nk. 

+#'   It is an equivalent of \code{p_type} for each leaf (i.e. cluster) that is applied 

+#'   after the determination of 'shared type genes'. 

+#'   }

+#'   

 #' @param ng # of genes to simulate. Importantly, for the library sizes 

 #'   computed by \code{\link{prepSim}} (= \code{exp(x$offset)}) to make sense, 

 #'   the number of simulated genes should match with the number of genes 

@@ -49,6 +62,12 @@

 #' @param force logical specifying whether to force 

 #'   simulation despite \code{ng != nrow(x)}.

 #'   

+#' @details 

+#'  The simulation of type genes can be performed in 2 ways; (1) by defining 

+#'  \code{p_type} and thus simulating independant clusters, OR (2) by defining both 

+#'  \code{phylo_tree} and \code{phylo_pars}, which will simulate a hierarchical structure 

+#'  between the clusters. Only one of the options is allowed. 

+#'  

 #' @return a \code{\link[SingleCellExperiment]{SingleCellExperiment}}

 #'   containing multiple clusters & samples across 2 groups.

 #' 

@@ -96,9 +115,11 @@

 #' plot(read.dendrogram(text = phylo_tree))

 #' 

 #' # simulate clusters accordingly

-#' sim <- simData(sce, phylo_tree = phylo_tree, ng = 500, force = TRUE)

+#' sim <- simData(ref, phylo_tree = phylo_tree, 

+#'   phylo_pars = list(c(0.1, 3), 0.1), 

+#'   ng = 500, force = TRUE)

 #' # view information about shared 'type' genes

-#' table(metadata(sim)$gene_info$shared_class)

+#' table(rowData(sim)$class)

 #' 

 #' @author Helena L Crowell

 #' 

@@ -124,7 +145,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     probs = NULL, p_dd = diag(6)[1, ], paired = FALSE,

     p_ep = 0.5, p_dp = 0.3, p_dm = 0.5,

     p_type = 0, lfc = 2, rel_lfc = NULL, 

-    phylo_tree = NULL, phylo_pars = c(0.2, 3),

+    phylo_tree = NULL, phylo_pars = list(c(0, 3), 0),

     ng = nrow(x), force = FALSE) {

     # throughout this code...

@@ -144,6 +165,18 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

         stop("Number of simulated genes should match with reference,\n", 

             "but 'ng != nrow(x)'; please specify 'force = TRUE' if\n", 

             "simulation should be forced regardlessly (see '?simData').")

+    if (!is.null(phylo_tree) && p_type != 0)

+        stop("Only one of 'p_type' and 'phylo_tree' can be provided.\n", 

+             "Please see the 'Details' section of '?simData'.")

+    if (!length(phylo_pars[[2]]) %in% c(1, nk))

+        stop("The second element of 'phylo_pars' should be correspond\n", 

+             " to the number of clusters ('nk') or of length 1.")

+    if (!is.null(phylo_tree) && phylo_pars[[1]][1] == 0)

+        warning("'phylo_pars[[1]][1]' has been set to 0;\n", 

+                "'phylo_tree' argument will be ignored.")

+    if (!is.null(phylo_tree) && all(phylo_pars[[2]] == 0))

+        warning("'phylo_pars[[2]]' has been set to 0;\n", 

+                "type genes for individual clusters won't be simulated.")

     # reference IDs

     nk0 <- length(kids0 <- set_names(levels(x$cluster_id)))

@@ -221,28 +254,6 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

         gs_by_k <- res$gs_by_k

         class  <- res$class

         specs <- res$specs

-        if (p_type != 0) {

-            # update gene indices 'gs_idx' to avoid imputing

-            # exclusive type genes in shared type genes

-            gs_idx_tmp <- gs_idx

-            for (c in c("ee", "ep"))

-                for (k in kids) {

-                    u <- gs_idx[[c, k]]

-                    gs_idx_tmp[[c, k]] <- u[!u %in% res$used]

-                }

-            # inpute cluster-specific type genes 

-            res <- .impute_type_genes(x, gs_by_k, gs_idx_tmp, p_type)

-            gs_by_k <- res$gs_by_k

-            # update genes classes & specificities

-            is_type <- res$class == "type"

-            # spot checks; any type-gene should not be of class 'shared'

-            # and cluster-specificities should be unassigned at this point

-            stopifnot(

-                all(class[is_type] == "state"),

-                all(is.na(unlist(specs[is_type]))))

-            class[is_type] <- "type"

-            specs[is_type] <- res$specs[is_type]

-        }

     # otherwise, simply impute type-genes w/o phylogeny

     } else if (p_type != 0) {

         res <- .impute_type_genes(x, gs_by_k, gs_idx, p_type)

@@ -159,7 +159,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

             stop("Some clusters appear to be missing from 'phylo_tree';\n",

                 "please make sure all clusters up to ", 

                 dQuote(kids_phylo[which.max(ns_phylo)]), " are present.")

-        if (nk_phylo != nk) nk <- nk_phylo

+        if (nk_phylo != nk) args$nk <- nk <- nk_phylo

     }

     # simulation IDs

@@ -212,7 +212,8 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

     # for ea. cluster, sample set of genes to simulate from

     gs_by_k <- setNames(sample(rownames(x), ng, TRUE), gs)

-    gs_by_k <- set_colnames(replicate(nk, gs_by_k), kids)

+    gs_by_k <- replicate(nk, gs_by_k)

+    colnames(gs_by_k) <- kids

     # when 'phylo_tree' is specified, induce hierarchical cluster structure

     if (!is.null(phylo_tree)) {                                  

@@ -343,7 +344,7 @@ simData <- function(x, nc = 2e3, ns = 3, nk = 3,

         rename("sim_mean.A" = "A", "sim_mean.B" = "B")

     # reorder

     o <- order(as.numeric(gsub("[a-z]", "", gi$gene)))

-    gi <- set_rownames(gi[o, ], NULL)

+    gi <- gi[o, ]; rownames(gi) <- NULL

     # construct SCE ------------------------------------------------------------

     # cell metadata including group, sample, cluster IDs

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

 #' across 2 experimental conditions from a real scRNA-seq data set.

 #' 

 #' @param x a \code{\link[SingleCellExperiment]{SingleCellExperiment}}.

-#' @param ng,nc,ns,nk # of genes, cells, samples and clusters to simulate. 

+#' @param nc,ns,nk # of cells, samples and clusters to simulate. 

 #' @param probs a list of length 3 containing probabilities of a cell belonging

 #'   to each cluster, sample, and group, respectively. List elements must be 

 #'   NULL (equal probabilities) or numeric values in [0, 1] that sum to 1.

@@ -27,13 +27,13 @@

 #'   Should be of length \code{nlevels(x$cluster_id)} with 

 #'   \code{levels(x$cluster_id)} as names. 

 #'   Defaults to factor of 1 for all clusters.

-#' @param cells_phylo newick tree text representing cluster relations 

+#' @param phylo_tree newick tree text representing cluster relations 

 #'   and their relative distance. If a tree is given, the distance between 

 #'   the clusters will be translated in the number of shared genes 

-#'   (this relation is controlled with \code{params_dist}). The distance 

+#'   (this relation is controlled with \code{phylo_pars}). The distance 

 #'   between two clusters is defined as the sum of the branches lengths 

 #'   separating them. 

-#' @param params_dist numeric vector of length 2. Defines the number of shared 

+#' @param phylo_pars numeric vector of length 2. Defines the number of shared 

 #'   genes as an adaptation of the exponential's PDF:

 #'   N = Ngenes x gamma1 * e^(-gamma2 x dist) ,

 #'   where gamma1 is the parameter that controls the percentage of shared genes. 

@@ -41,6 +41,13 @@

 #'   depending on the input prep_sce. gamma2 is the 'penalty' of increasing 

 #'   the distance between clusters, applied on the number of shared genes. 

 #'   Default to -3. 

+#' @param ng # of genes to simulate. Importantly, for the library sizes 

+#'   computed by \code{\link{prepSim}} (= \code{exp(x$offset)}) to make sense, 

+#'   the number of simulated genes should match with the number of genes 

+#'   in the reference. To simulate a reduced number of genes, e.g. for 

+#'   testing and development purposes, please set \code{force = TRUE}.

+#' @param force logical specifying whether to force 

+#'   simulation despite \code{ng != nrow(x)}.

 #'   

 #' @return a \code{\link[SingleCellExperiment]{SingleCellExperiment}}

 #'   containing multiple clusters & samples across 2 groups.

@@ -50,11 +57,13 @@

 #' library(SingleCellExperiment)

 #' 

 #' # prep. SCE for simulation

-#' sce <- prepSim(sce)

+#' ref <- prepSim(sce)

 #' 

 #' # simulate data

-#' (sim <- simData(sce, ng = 100, nc = 10,

-#'   p_dd = c(0.9, 0, 0.1, 0, 0, 0)))

+#' (sim <- simData(ref, nc = 10,

+#'   p_dd = c(0.9, 0, 0.1, 0, 0, 0),

+#'   ng = 100, force = TRUE,

+#'   probs = list(NULL, NULL, c(1, 0))))

 #' 

 #' # simulation metadata

 #' head(gi <- metadata(sim)$gene_info)

@@ -63,30 +72,32 @@

 #' table(gi$category)

 #' 

 #' # unbalanced sample sizes

-#' sim <- simData(sce, ng = 10, nc = 100,

-#'   probs = list(NULL, c(0.25, 0.75), NULL))

+#' sim <- simData(ref, nc = 100,

+#'   probs = list(NULL, c(0.25, 0.75), NULL),