From fbf8f8ecb1da6e66b4ca91c99a46b20e14ab9e01 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