From effef4066bb16fed03e466d1a61448a3382b5386 Mon Sep 17 00:00:00 2001
From: Luke Zappia <lazappi@users.noreply.github.com>
Date: Tue, 16 Jul 2019 17:53:17 +1000
Subject: [PATCH] Add paths parameters

---
 DESCRIPTION               |  2 +-
 NEWS.md                   |  4 +++
 R/AllClasses.R            | 21 ++++++++++--
 R/Params-methods.R        |  7 ++--
 R/SplotchParams-methods.R | 67 ++++++++++++++++++++++++++++++++++-----
 R/params-functions.R      | 14 +++++---
 R/splotch-simulate.R      | 49 ++++++++++++++++++++++++++++
 man/SplotchParams.Rd      |  7 ++++
 8 files changed, 154 insertions(+), 17 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 39cff9a..f3581ec 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: splatter
 Type: Package
 Title: Simple Simulation of Single-cell RNA Sequencing Data
-Version: 1.9.2.9001
+Version: 1.9.2.9002
 Date: 2019-07-16
 Author: Luke Zappia
 Authors@R:
diff --git a/NEWS.md b/NEWS.md
index b1ae280..db32f1f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+### Version 1.9.2.9002 (2019-07-16)
+
+* Add paths parameters
+
 ### Version 1.9.2.9001 (2019-07-16)
 
 * Add mean parameters
diff --git a/R/AllClasses.R b/R/AllClasses.R
index d208885..35fc2af 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -283,6 +283,13 @@ setClass("SplatParams",
 #'             network.}
 #'         }
 #'     }
+#'     \item{\emph{Paths parameters}}{
+#'         \describe{
+#'             \item{\code{[paths.programs]}}{Number of expression programs.}
+#'             \item{\code{[paths.design]}}{data.frame describing path
+#'             structure.}
+#'         }
+#'     }
 #' }
 #'
 #' The parameters not shown in brackets can be estimated from real data using
@@ -300,13 +307,23 @@ setClass("SplotchParams",
                    mean.values = "numeric",
                    network.graph = "ANY",
                    network.nRegs = "numeric",
-                   network.regsSet = "logical"),
+                   network.regsSet = "logical",
+                   paths.nPrograms = "numeric",
+                   paths.design = "data.frame",
+                   paths.means = "list"),
          prototype = prototype(mean.rate = 0.3,
                                mean.shape = 0.6,
                                mean.values = numeric(),
                                network.graph = NULL,
                                network.nRegs = 100,
-                               network.regsSet = FALSE))
+                               network.regsSet = FALSE,
+                               paths.nPrograms = 10,
+                               paths.design = data.frame(
+                                   Path = 1,
+                                   From = 0,
+                                   Steps = 100
+                               ),
+                               paths.means = list()))
 
 #' The LunParams class
 #'
diff --git a/R/Params-methods.R b/R/Params-methods.R
index 4f01bfe..6bb8a00 100644
--- a/R/Params-methods.R
+++ b/R/Params-methods.R
@@ -48,8 +48,11 @@ setMethod("show", "Params", function(object) {
 
     cat("A", crayon::bold("Params"), "object of class",
         crayon::bold(class(object)), "\n")
-    cat("Parameters can be (estimable) or", crayon::blue("[not estimable],"),
-        "'Default' or ", crayon::bold(crayon::green("'NOT DEFAULT'")), "\n\n")
+    cat("Parameters can be (estimable) or",
+        paste0(crayon::blue("[not estimable]"), ","),
+        "'Default' or ", crayon::bold(crayon::green("'NOT DEFAULT'")), "\n")
+    cat(crayon::bgYellow(crayon::black("Secondary")), "parameters are usually",
+        "set during simulation\n\n")
     showPP(object, pp)
     cat(length(slotNames(object)) - 3, "additional parameters", "\n\n")
 })
diff --git a/R/SplotchParams-methods.R b/R/SplotchParams-methods.R
index 9727e82..53ab297 100644
--- a/R/SplotchParams-methods.R
+++ b/R/SplotchParams-methods.R
@@ -26,11 +26,16 @@ setValidity("SplotchParams", function(object) {
                                                       null.ok = TRUE),
                 network.nRegs = checkmate::checkInt(v$network.nRegs,
                                                     lower = 0),
-                network.regsSet = checkmate::checkFlag(v$network.regsSet))
-
-    if (checkmate::testNumeric(v$mean.values, len = 0)) {
-        checks <- c(checks, mean.values = TRUE)
-    } else {
+                network.regsSet = checkmate::checkFlag(v$network.regsSet),
+                paths.nPrograms = checkmate::checkInt(v$paths.nPrograms,
+                                                      lower = 1),
+                paths.design = checkmate::checkDataFrame(v$paths.design,
+                                                         types = "numeric",
+                                                         any.missing = FALSE,
+                                                         min.rows = 1,
+                                                         ncols = 3))
+
+    if (!checkmate::testNumeric(v$mean.values, len = 0)) {
         checks <- c(checks,
                     mean.values = checkmate::checkNumeric(v$mean.values,
                                                           lower = 0,
@@ -39,6 +44,35 @@ setValidity("SplotchParams", function(object) {
                                                           len = v$nGenes))
     }
 
+    if (!all(colnames(v$paths.design) == c("Path", "From", "Steps"))) {
+        checks <- c(checks, paths.design = "Incorrect column names")
+    } else {
+        if (!(0 %in% v$paths.design$From)) {
+            checks <- c(checks, paths.design = paste("origin must be specified",
+                      "in paths.design"))
+        }
+
+        paths.graph <- igraph::graph_from_data_frame(v$paths.design)
+        if (!igraph::is_simple(paths.graph)) {
+            checks <- c(checks, paths.design = "graph is not simple")
+        }
+        if (!igraph::is_connected(paths.graph)) {
+            checks <- c(checks, paths.design = "graph is not connected")
+        }
+        if (!igraph::is_dag(paths.graph)) {
+            checks <- c(checks, paths.design = "graph is not a DAG")
+        }
+    }
+
+    if (!checkmate::testList(v$paths.means, len = 0)) {
+        checks <- c(checks,
+                    paths.means = checkmate::checkList(v$paths.means,
+                                                       types = "matrix",
+                                                       any.missing = FALSE,
+                                                       len = nrow(v$paths.design),
+                                                       names = "unique"))
+    }
+
     if (all(checks == TRUE)) {
         valid <- TRUE
     } else {
@@ -52,14 +86,22 @@ setValidity("SplotchParams", function(object) {
 #' @importFrom methods show
 setMethod("show", "SplotchParams", function(object) {
 
-    pp.top <- list("Mean:" = c("(Rate)"   = "mean.rate",
-                               "(Shape)"  = "mean.shape",
-                               "[Values]" = "mean.values"))
+    pp.top <- list("Mean:" = c("(Rate)"    = "mean.rate",
+                               "(Shape)"   = "mean.shape",
+                               "[Values]*" = "mean.values"))
 
     pp.network <- list("Network:" = c("[Graph]"   = "network.graph",
                                       "[nRegs]"   = "network.nRegs",
                                       "[regsSet]" = "network.regsSet"))
 
+    pp.paths <- list("Paths:" = c("[nPrograms]" = "paths.nPrograms",
+                                  "[Design]"  = "paths.design"))
+
+    paths.means <- getParam(object, "paths.means")
+    if (length(paths.means) == 0) {
+        pp.paths[[1]] <- c(pp.paths[[1]], "[Means]*" = "paths.means")
+    }
+
     callNextMethod()
 
     showPP(object, pp.top)
@@ -82,6 +124,15 @@ setMethod("show", "SplotchParams", function(object) {
                              network.values$`[regsSet]` != FALSE)
         showValues(network.values, network.default)
     }
+
+    showPP(object, pp.paths)
+
+    if (length(paths.means) != 0) {
+        cat(crayon::bgYellow(crayon::bold(crayon::blue("[MEANS]\n"))))
+        cat(crayon::bold(crayon::green(paste(
+            "List of", length(paths.means), "matrices"
+        ))))
+    }
 })
 
 #' @rdname setParam
diff --git a/R/params-functions.R b/R/params-functions.R
index c969ba8..5e7eed9 100644
--- a/R/params-functions.R
+++ b/R/params-functions.R
@@ -122,15 +122,19 @@ showValues <- function(values, not.default) {
     items.per.line <- floor(screen.width / (max.len + 2))
 
     short.names <- names(short.values)
-    short.values <- crayon::col_align(short.values, max.len, "right")
-    short.names <- crayon::col_align(short.names, max.len, "right")
-
     not.est <- !grepl("\\(", short.names)
+    secondary <- grepl("\\*", short.names)
+    short.names <- gsub("\\*", "", short.names)
+
     short.names[not.est] <- crayon::blue(short.names[not.est])
+    short.names[secondary] <- crayon::bgYellow(short.names[secondary])
     short.names[not.default] <- crayon::bold(short.names[not.default])
     short.values[not.default] <- crayon::green(short.values[not.default])
     short.values[not.default] <- crayon::bold(short.values[not.default])
 
+    short.values <- crayon::col_align(short.values, max.len, "right")
+    short.names <- crayon::col_align(short.names, max.len, "right")
+
     names(short.values) <- short.names
 
     values.list <- split(short.values,
@@ -177,6 +181,8 @@ showDFs <- function(dfs, not.default) {
         cat(paste0("\n", name, "\n"))
         cat(msg, "\n")
         print(head(df, n = 4))
-        cat("# ... with", nrow(df) - 4, "more rows\n")
+        if (nrow(df) > 4) {
+            cat("# ... with", nrow(df) - 4, "more rows\n")
+        }
     }
 }
diff --git a/R/splotch-simulate.R b/R/splotch-simulate.R
index 4a7c5eb..e91e4b7 100644
--- a/R/splotch-simulate.R
+++ b/R/splotch-simulate.R
@@ -45,6 +45,7 @@ splotchSimulate <- function(params = newSplotchParams(), verbose = TRUE, ...) {
         if (verbose) {message("Using selected regulators...")}
     }
 
+    # Generate means
     if (length(getParam(params, "mean.values")) == 0) {
         if (verbose) {message("Simulating means...")}
         mean.shape <- getParam(params, "mean.shape")
@@ -55,6 +56,54 @@ splotchSimulate <- function(params = newSplotchParams(), verbose = TRUE, ...) {
         if (verbose) {message("Using defined means...")}
     }
 
+    # Generate paths
+    paths.design <- getParam(params, "paths.design")
+    network.graph <- getParam(params, "network.graph")
+    network.weights <- igraph::as_adjacency_matrix(network.graph,
+                                                   attr = "weight")
+    network.nRegs <- getParam(params, "network.nRegs")
+    network.isReg <- igraph::vertex_attr(network.graph, "IsReg")
+    paths.nPrograms <- getParam(params, "paths.nPrograms")
+    programs.weights <- matrix(rnorm(network.nRegs * paths.nPrograms),
+                               nrow = network.nRegs, ncol = paths.nPrograms)
+    paths.changes <- vector("list", nrow(paths.design))
+    paths.factors <- vector("list", nrow(paths.design))
+    for (path in seq_len(nrow(paths.design))) {
+        nSteps <- paths.design$Steps[path]
+        from <- paths.design$From[path]
+        changes <- matrix(0, nrow = nGenes, ncol = nSteps + 1)
+
+        if (from != 0) {
+            from.changes <- paths.changes[[from]]
+            changes[, 1] <- from.changes[, ncol(from.changes)]
+        }
+
+        for (step in seq_len(nSteps) + 1) {
+            programs.changes <- rnorm(paths.nPrograms)
+            reg.changes <- as.vector(programs.weights %*% programs.changes)
+            changes[network.isReg, step] <- reg.changes
+            change <- as.vector(changes[, step - 1] %*% network.weights)
+            changes[, step] <- changes[, step] + change
+        }
+
+        if (from == 0) {
+            changes <- changes[, 1:nSteps]
+            factors <- matrixStats::rowCumsums(changes)
+        } else {
+            changes <- changes[, 2:nSteps + 1]
+            from.factors <- paths.factors[[from]][, ncol(paths.factors[[from]])]
+            factors <- matrixStats::rowCumsums(changes) + from.factors
+        }
+        paths.changes[[path]] <- changes
+        paths.factors[[path]] <- factors
+    }
+    means.values <- getParam(params, "mean.values")
+    paths.means <- lapply(paths.factors, function(x) {
+        2 ^ x * means.values
+    })
+    names(paths.means) <- paste0("Path", paths.design$Path)
+    params <- setParam(params, "paths.means", paths.means)
+
     # Simulate base gene means
 
     # sim <- SingleCellExperiment(assays = list(counts = counts),
diff --git a/man/SplotchParams.Rd b/man/SplotchParams.Rd
index e03a62b..20c32bd 100644
--- a/man/SplotchParams.Rd
+++ b/man/SplotchParams.Rd
@@ -33,6 +33,13 @@ The Splotch simulation uses the following parameters:
             network.}
         }
     }
+    \item{\emph{Paths parameters}}{
+        \describe{
+            \item{\code{[paths.programs]}}{Number of expression programs.}
+            \item{\code{[paths.design]}}{data.frame describing path
+            structure.}
+        }
+    }
 }
 
 The parameters not shown in brackets can be estimated from real data using
-- 
GitLab