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

---
 DESCRIPTION               |   4 +-
 NAMESPACE                 |   4 ++
 NEWS.md                   |   4 ++
 R/AllClasses.R            |  34 +++++++++---
 R/SplotchParams-methods.R | 108 ++++++++++++++++++++++++++++++++++----
 R/params-functions.R      |   6 ++-
 R/splotch-simulate.R      | 100 +++++++++++++++++++++++++++++++++++
 man/SplatParams.Rd        |   6 +--
 man/SplotchParams.Rd      |  18 ++++++-
 man/newParams.Rd          |   6 ++-
 man/setParam.Rd           |   5 +-
 man/setParams.Rd          |   5 +-
 man/showValues.Rd         |   4 +-
 man/splotchSimulate.Rd    |  29 ++++++++++
 14 files changed, 303 insertions(+), 30 deletions(-)
 create mode 100644 R/splotch-simulate.R
 create mode 100644 man/splotchSimulate.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 8cd1713..39cff9a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: splatter
 Type: Package
 Title: Simple Simulation of Single-cell RNA Sequencing Data
-Version: 1.9.2.9000
-Date: 2019-07-11
+Version: 1.9.2.9001
+Date: 2019-07-16
 Author: Luke Zappia
 Authors@R:
     c(person("Luke", "Zappia", role = c("aut", "cre"),
diff --git a/NAMESPACE b/NAMESPACE
index c9b65bd..afefae4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -47,6 +47,7 @@ export(newSCDDParams)
 export(newSimpleParams)
 export(newSparseDCParams)
 export(newSplatParams)
+export(newSplotchParams)
 export(newZINBParams)
 export(phenoEstimate)
 export(phenoSimulate)
@@ -63,6 +64,7 @@ export(splatSimulate)
 export(splatSimulateGroups)
 export(splatSimulatePaths)
 export(splatSimulateSingle)
+export(splotchSimulate)
 export(summariseDiff)
 export(zinbEstimate)
 export(zinbSimulate)
@@ -75,6 +77,7 @@ exportClasses(SCDDParams)
 exportClasses(SimpleParams)
 exportClasses(SparseDCParams)
 exportClasses(SplatParams)
+exportClasses(SplotchParams)
 exportClasses(ZINBParams)
 importFrom(BiocParallel,SerialParam)
 importFrom(BiocParallel,bplapply)
@@ -125,6 +128,7 @@ importFrom(methods,"slot<-")
 importFrom(methods,as)
 importFrom(methods,callNextMethod)
 importFrom(methods,new)
+importFrom(methods,show)
 importFrom(methods,slot)
 importFrom(methods,slotNames)
 importFrom(methods,validObject)
diff --git a/NEWS.md b/NEWS.md
index d932de8..b1ae280 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+### Version 1.9.2.9001 (2019-07-16)
+
+* Add mean parameters
+
 ### Version 1.9.2.9000 (2019-07-11)
 
 * Add SplotchParams object
diff --git a/R/AllClasses.R b/R/AllClasses.R
index a251d3f..d208885 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -267,8 +267,22 @@ setClass("SplatParams",
 #'     \item{\code{nGenes}}{The number of genes to simulate.}
 #'     \item{\code{nCells}}{The number of cells to simulate.}
 #'     \item{\code{[seed]}}{Seed to use for generating random numbers.}
-#'     \item{\code{[network.graph]}}{Graph containing the gene network.}
-#'     \item{\code{[network.nRegs]}}{Number of regulators in a the network.}
+#'     \item{\emph{Mean parameters}}{
+#'         \describe{
+#'             \item{\code{mean.shape}}{Shape parameter for the mean gamma
+#'             distribution.}
+#'             \item{\code{mean.rate}}{Rate parameter for the mean gamma
+#'             distribution.}
+#'             \item{\code{mean.values}}{Vector of means for each gene.}
+#'         }
+#'     }
+#'     \item{\emph{Network parameters}}{
+#'         \describe{
+#'             \item{\code{[network.graph]}}{Graph containing the gene network.}
+#'             \item{\code{[network.nRegs]}}{Number of regulators in the
+#'             network.}
+#'         }
+#'     }
 #' }
 #'
 #' The parameters not shown in brackets can be estimated from real data using
@@ -281,10 +295,18 @@ setClass("SplatParams",
 #' @exportClass SplotchParams
 setClass("SplotchParams",
          contains = "Params",
-         slots = c(network.graph = "ANY",
-                   network.nRegs = "numeric"),
-         prototype = prototype(network.graph = NULL,
-                               network.nRegs = 100))
+         slots = c(mean.shape = "numeric",
+                   mean.rate = "numeric",
+                   mean.values = "numeric",
+                   network.graph = "ANY",
+                   network.nRegs = "numeric",
+                   network.regsSet = "logical"),
+         prototype = prototype(mean.rate = 0.3,
+                               mean.shape = 0.6,
+                               mean.values = numeric(),
+                               network.graph = NULL,
+                               network.nRegs = 100,
+                               network.regsSet = FALSE))
 
 #' The LunParams class
 #'
diff --git a/R/SplotchParams-methods.R b/R/SplotchParams-methods.R
index b5373f1..9727e82 100644
--- a/R/SplotchParams-methods.R
+++ b/R/SplotchParams-methods.R
@@ -17,13 +17,27 @@ setValidity("SplotchParams", function(object) {
 
     v <- getParams(object, slotNames(object))
 
-    checks <- c(nGenes = checkmate::checkInt(v$nGenes, lower = 1),
+    checks <- c(#nGenes = checkmate::checkInt(v$nGenes, lower = 1),
                 nCells = checkmate::checkInt(v$nCells, lower = 1),
                 seed = checkmate::checkInt(v$seed, lower = 0),
+                mean.rate = checkmate::checkNumber(v$mean.rate, lower = 0),
+                mean.shape = checkmate::checkNumber(v$mean.shape, lower = 0),
                 network.graph = checkmate::checkClass(v$network, "igraph",
                                                       null.ok = TRUE),
                 network.nRegs = checkmate::checkInt(v$network.nRegs,
-                                                    lower = 0))
+                                                    lower = 0),
+                network.regsSet = checkmate::checkFlag(v$network.regsSet))
+
+    if (checkmate::testNumeric(v$mean.values, len = 0)) {
+        checks <- c(checks, mean.values = TRUE)
+    } else {
+        checks <- c(checks,
+                    mean.values = checkmate::checkNumeric(v$mean.values,
+                                                          lower = 0,
+                                                          finite = TRUE,
+                                                          any.missing = FALSE,
+                                                          len = v$nGenes))
+    }
 
     if (all(checks == TRUE)) {
         valid <- TRUE
@@ -35,13 +49,21 @@ setValidity("SplotchParams", function(object) {
     return(valid)
 })
 
+#' @importFrom methods show
 setMethod("show", "SplotchParams", function(object) {
 
-    pp.network <- list("Network:" = c("[Graph]" = "network.graph",
-                                      "[nRegs]" = "network.nRegs"))
+    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"))
 
     callNextMethod()
 
+    showPP(object, pp.top)
+
     network.graph <- getParam(object, "network.graph")
     if (is.null(network.graph)) {
         showPP(object, pp.network)
@@ -53,14 +75,80 @@ setMethod("show", "SplotchParams", function(object) {
             igraph::gsize(network.graph), "edges\n"
         ))))
         show(network.graph)
-        network.nRegs <- getParam(object, "network.nRegs")
-        names(network.nRegs) <- c("[nRegs]")
-        if (network.nRegs == 100) {
-            showValues(network.nRegs, FALSE)
+        network.values <- list("[nRegs]" = getParam(object, "network.nRegs"),
+                               "[regsSet]" = getParam(object,
+                                                      "network.regsSet"))
+        network.default <- c(network.values$`[nRegs]` != 100,
+                             network.values$`[regsSet]` != FALSE)
+        showValues(network.values, network.default)
+    }
+})
+
+#' @rdname setParam
+setMethod("setParam", "SplotchParams", function(object, name, value) {
+    checkmate::assertString(name)
+
+    if (name == "nGenes") {
+        if (!is.null(getParam(object, "network.graph"))) {
+            stop("nGenes can not be changed once network.graph is set")
+        }
+        if (!length(getParam(object, "mean.values")) == 0) {
+            stop("nGenes can not be changed once mean.values is set")
+        }
+    }
+
+    if (name == "mean.values") {
+        if (!is.null(getParam(object, "network.graph"))) {
+            if (length(value) != getParam(object, "nGenes")) {
+                stop("new mean.values does not match the number of genes in ",
+                     "the network")
+            }
         } else {
-            showValues(network.nRegs, TRUE)
+            object <- setParam(object, "nGenes", length(value))
         }
     }
 
-    # showPP(object, pp)
+    if (name == "network.graph") {
+        checkmate::assertClass(value, "igraph")
+        object <- setParamUnchecked(object, "nGenes", igraph::gorder(value))
+        if (!(length(getParam(object, "mean.values")) == 0)) {
+            warning("changing network.graph resets mean.values")
+            object <- setParam(object, "mean.values", numeric())
+        }
+
+        if ("IsReg" %in% igraph::vertex_attr_names(value)) {
+            isReg <- igraph::vertex_attr(value, "IsReg")
+            checkmate::assertLogical(isReg, any.missing = FALSE)
+            object <- setParamUnchecked(object, "network.nRegs", sum(isReg))
+            object <- setParamUnchecked(object, "network.regsSet", TRUE)
+        } else {
+            object <- setParamUnchecked(object, "network.regsSet", FALSE)
+        }
+    }
+
+    if (name == "network.nRegs" && getParam(object, "network.regsSet")) {
+        stop("network.nRegs can not be changed after regulators are set")
+    }
+
+    if (name == "network.regsSet") {
+        stop("network.regsSet can not be set manually")
+    }
+
+    object <- callNextMethod()
+
+    return(object)
+})
+
+#' @rdname setParams
+setMethod("setParams", "SplotchParams", function(object, update = NULL, ...) {
+
+    checkmate::assertList(update, null.ok = TRUE)
+
+    update <- c(update, list(...))
+
+    update <- bringItemsForward(update, c("network.graph"))
+
+    object <- callNextMethod(object, update)
+
+    return(object)
 })
diff --git a/R/params-functions.R b/R/params-functions.R
index 316c15a..c969ba8 100644
--- a/R/params-functions.R
+++ b/R/params-functions.R
@@ -76,8 +76,10 @@ showPP <- function(params, pp) {
         not.default <- sapply(seq_along(values), function(i) {
             !identical(values[i], default.values[i])
         })
-        null.values <- sapply(values, is.null)
-        values[null.values] <- "Not set"
+        empty.values <- sapply(values, function(x) {
+            is.null(x) || length(x) == 0
+        })
+        values[empty.values] <- "Not set"
 
         cat(crayon::bold(category), "\n")
         if (sum(!is.df) > 0) {
diff --git a/R/splotch-simulate.R b/R/splotch-simulate.R
new file mode 100644
index 0000000..4a7c5eb
--- /dev/null
+++ b/R/splotch-simulate.R
@@ -0,0 +1,100 @@
+#' Splotch simulation
+#'
+#' Simulate counts from...
+#'
+#' @param params SplotchParams object containing simulation parameters.
+#' @param verbose logical. Whether to print progress messages
+#' @param ... any additional parameter settings to override what is provided in
+#'        \code{params}.
+#'
+#' @details
+#' Details...
+#'
+#' @return SingleCellExperiment containing simulated counts
+#' @examples
+#' sim <- splotchSimulate()
+#'
+#' @export
+#' @importFrom SingleCellExperiment SingleCellExperiment
+splotchSimulate <- function(params = newSplotchParams(), verbose = TRUE, ...) {
+
+    checkmate::assertClass(params, "SplotchParams")
+    params <- setParams(params, ...)
+
+    # Set random seed
+    seed <- getParam(params, "seed")
+    set.seed(seed)
+
+    # Get the parameters we are going to use
+    nGenes <- getParam(params, "nGenes")
+    network.graph <- getParam(params, "network.graph")
+
+    # Generate network
+    if (is.null(network.graph)) {
+        network.graph <- generateNetwork(nGenes, verbose)
+        params <- setParam(params, "network.graph", network.graph)
+    }
+
+    # Select regulators
+    if (!getParam(params, "network.regsSet")) {
+        network.nRegs <- getParam(params, "network.nRegs")
+        network.graph <- selectRegulators(network.graph, network.nRegs,
+                                          verbose)
+        params <- setParam(params, "network.graph", network.graph)
+    } else {
+        if (verbose) {message("Using selected regulators...")}
+    }
+
+    if (length(getParam(params, "mean.values")) == 0) {
+        if (verbose) {message("Simulating means...")}
+        mean.shape <- getParam(params, "mean.shape")
+        mean.rate <- getParam(params, "mean.rate")
+        mean.values <- rgamma(nGenes, shape = mean.shape, rate = mean.rate)
+        params <- setParam(params, "mean.values", mean.values)
+    } else {
+        if (verbose) {message("Using defined means...")}
+    }
+
+    # Simulate base gene means
+
+    # sim <- SingleCellExperiment(assays = list(counts = counts),
+    #                             rowData = features,
+    #                             colData = cells,
+    #                             metadata = list(params = params))
+    #
+    # return(sim)
+
+    return(params)
+}
+
+generateNetwork <- function(n.nodes, verbose) {
+
+    if (verbose) {message("Generating gene network...")}
+
+    graph.raw <- igraph::sample_forestfire(n.nodes, 0.1)
+    graph.data <- igraph::get.data.frame(graph.raw)
+    graph.data <- graph.data[, c("from", "to")]
+    graph.data$weight <- rnorm(nrow(graph.data))
+    graph <- igraph::graph.data.frame(graph.data)
+    graph <- igraph::set_vertex_attr(graph, "mean",
+                                     value = rnorm(igraph::gorder(graph)))
+
+    return(graph)
+}
+
+selectRegulators <- function(graph, nReg, verbose) {
+
+    if (verbose) {message("Selecting regulators...")}
+
+    out.degree <- igraph::degree(graph, mode = "out")
+    in.degree <- igraph::degree(graph, mode = "in")
+    reg.prob <-  out.degree - in.degree
+    reg.prob <- reg.prob + rnorm(length(reg.prob))
+    reg.prob[reg.prob <= 0] <- 1e-10
+    reg.prob <- reg.prob / sum(reg.prob)
+    reg.nodes <- names(rev(sort(reg.prob))[1:nReg])
+    is.reg <- igraph::V(graph)$name %in% reg.nodes
+    graph <- igraph::set_vertex_attr(graph, "IsReg", value = is.reg)
+
+    return(graph)
+}
diff --git a/man/SplatParams.Rd b/man/SplatParams.Rd
index bb906f9..f2954c5 100644
--- a/man/SplatParams.Rd
+++ b/man/SplatParams.Rd
@@ -6,12 +6,12 @@
 \alias{SplatParams-class}
 \title{The SplatParams class}
 \description{
-S4 class that holds parameters for the Splatter simulation.
+S4 class that holds parameters for the Splat simulation.
 }
 \section{Parameters}{
 
 
-The Splatter simulation requires the following parameters:
+The Splat simulation requires the following parameters:
 
 \describe{
     \item{\code{nGenes}}{The number of genes to simulate.}
@@ -133,7 +133,7 @@ The Splatter simulation requires the following parameters:
 }
 
 The parameters not shown in brackets can be estimated from real data using
-\code{\link{splatEstimate}}. For details of the Splatter simulation
+\code{\link{splatEstimate}}. For details of the Splat simulation
 see \code{\link{splatSimulate}}.
 }
 
diff --git a/man/SplotchParams.Rd b/man/SplotchParams.Rd
index e851d29..e03a62b 100644
--- a/man/SplotchParams.Rd
+++ b/man/SplotchParams.Rd
@@ -17,8 +17,22 @@ The Splotch simulation uses the following parameters:
     \item{\code{nGenes}}{The number of genes to simulate.}
     \item{\code{nCells}}{The number of cells to simulate.}
     \item{\code{[seed]}}{Seed to use for generating random numbers.}
-    \item{\code{[network.graph]}}{Graph containing the gene network.}
-    \item{\code{[network.nRegs]}}{Number of regulators in a the network.}
+    \item{\emph{Mean parameters}}{
+        \describe{
+            \item{\code{mean.shape}}{Shape parameter for the mean gamma
+            distribution.}
+            \item{\code{mean.rate}}{Rate parameter for the mean gamma
+            distribution.}
+            \item{\code{mean.values}}{Vector of means for each gene.}
+        }
+    }
+    \item{\emph{Network parameters}}{
+        \describe{
+            \item{\code{[network.graph]}}{Graph containing the gene network.}
+            \item{\code{[network.nRegs]}}{Number of regulators in the
+            network.}
+        }
+    }
 }
 
 The parameters not shown in brackets can be estimated from real data using
diff --git a/man/newParams.Rd b/man/newParams.Rd
index 5af67ec..6c31e96 100644
--- a/man/newParams.Rd
+++ b/man/newParams.Rd
@@ -2,7 +2,8 @@
 % Please edit documentation in R/AllGenerics.R, R/BASiCSParams-methods.R,
 %   R/Lun2Params-methods.R, R/LunParams-methods.R, R/MFAParams-methods.R,
 %   R/PhenoParams-methods.R, R/SCDDParams-methods.R, R/SimpleParams-methods.R,
-%   R/SparseDCParams-methods.R, R/SplatParams-methods.R, R/ZINBParams-methods.R
+%   R/SparseDCParams-methods.R, R/SplatParams-methods.R,
+%   R/SplotchParams-methods.R, R/ZINBParams-methods.R
 \name{newParams}
 \alias{newParams}
 \alias{newBASiCSParams}
@@ -14,6 +15,7 @@
 \alias{newSimpleParams}
 \alias{newSparseDCParams}
 \alias{newSplatParams}
+\alias{newSplotchParams}
 \alias{newZINBParams}
 \title{New Params}
 \usage{
@@ -35,6 +37,8 @@ newSparseDCParams(...)
 
 newSplatParams(...)
 
+newSplotchParams(...)
+
 newZINBParams(...)
 }
 \arguments{
diff --git a/man/setParam.Rd b/man/setParam.Rd
index e144b40..0aca893 100644
--- a/man/setParam.Rd
+++ b/man/setParam.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/AllGenerics.R, R/BASiCSParams-methods.R,
 %   R/Lun2Params-methods.R, R/LunParams-methods.R, R/Params-methods.R,
 %   R/PhenoParams-methods.R, R/SCDDParams-methods.R, R/SplatParams-methods.R,
-%   R/ZINBParams-methods.R
+%   R/SplotchParams-methods.R, R/ZINBParams-methods.R
 \docType{methods}
 \name{setParam}
 \alias{setParam}
@@ -13,6 +13,7 @@
 \alias{setParam,PhenoParams-method}
 \alias{setParam,SCDDParams-method}
 \alias{setParam,SplatParams-method}
+\alias{setParam,SplotchParams-method}
 \alias{setParam,ZINBParams-method}
 \title{Set a parameter}
 \usage{
@@ -32,6 +33,8 @@ setParam(object, name, value)
 
 \S4method{setParam}{SplatParams}(object, name, value)
 
+\S4method{setParam}{SplotchParams}(object, name, value)
+
 \S4method{setParam}{ZINBParams}(object, name, value)
 }
 \arguments{
diff --git a/man/setParams.Rd b/man/setParams.Rd
index b7265f0..ad217dd 100644
--- a/man/setParams.Rd
+++ b/man/setParams.Rd
@@ -1,11 +1,12 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/AllGenerics.R, R/Params-methods.R,
-%   R/SplatParams-methods.R
+%   R/SplatParams-methods.R, R/SplotchParams-methods.R
 \docType{methods}
 \name{setParams}
 \alias{setParams}
 \alias{setParams,Params-method}
 \alias{setParams,SplatParams-method}
+\alias{setParams,SplotchParams-method}
 \title{Set parameters}
 \usage{
 setParams(object, update = NULL, ...)
@@ -13,6 +14,8 @@ setParams(object, update = NULL, ...)
 \S4method{setParams}{Params}(object, update = NULL, ...)
 
 \S4method{setParams}{SplatParams}(object, update = NULL, ...)
+
+\S4method{setParams}{SplotchParams}(object, update = NULL, ...)
 }
 \arguments{
 \item{object}{Params object to set parameters in.}
diff --git a/man/showValues.Rd b/man/showValues.Rd
index ddd813a..9167cc7 100644
--- a/man/showValues.Rd
+++ b/man/showValues.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/params-functions.R
 \name{showValues}
 \alias{showValues}
-\title{Show vales}
+\title{Show values}
 \usage{
 showValues(values, not.default)
 }
@@ -12,5 +12,5 @@ showValues(values, not.default)
 \item{not.default}{logical vector giving which have changed from the default.}
 }
 \description{
-Function used for pretty printing scale or vector parameters.
+Function used for pretty printing scalar or vector parameters.
 }
diff --git a/man/splotchSimulate.Rd b/man/splotchSimulate.Rd
new file mode 100644
index 0000000..19458d9
--- /dev/null
+++ b/man/splotchSimulate.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/splotch-simulate.R
+\name{splotchSimulate}
+\alias{splotchSimulate}
+\title{Splotch simulation}
+\usage{
+splotchSimulate(params = newSplotchParams(), verbose = TRUE, ...)
+}
+\arguments{
+\item{params}{SplotchParams object containing simulation parameters.}
+
+\item{verbose}{logical. Whether to print progress messages}
+
+\item{...}{any additional parameter settings to override what is provided in
+\code{params}.}
+}
+\value{
+SingleCellExperiment containing simulated counts
+}
+\description{
+Simulate counts from...
+}
+\details{
+Details...
+}
+\examples{
+sim <- splotchSimulate()
+
+}
-- 
GitLab