Newer
Older
#' Simple simulation
#'
#' Simulate counts from a simple negative binomial distribution without
#' simulated library sizes, differential expression etc.
#'
#' @param params splatParams 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
#' Uses the following parameters: \code{nCells}, \code{nGenes},
#' \code{mean.shape}, \code{mean.rate}, \code{bcv.common}.
#'
#' Gene means are simulated from a gamma distribution with
#' \code{shape = mean.shape} and \code{rate = mean.rate}. Counts are then
#' simulated from a negative binomial distribution with \code{mu = means} and
#' \code{size = 1 / bcv.common}.
#'
#' Parameters are set in the tiered manner described in \code{\link{splat}}.
#'
#' @return SCESet containing simulated counts
#' @examples
#' sim <- simSimple()
#' @export
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
simSimple <- function(params = defaultParams(), verbose = TRUE, ...) {
if (verbose) {message("Getting parameters...")}
params <- setParams(params, ...)
params <- mergeParams(params, defaultParams())
params <- expandParams(params)
# Set random seed
seed <- getParams(params, "seed")
set.seed(seed)
# Get the parameters we are going to use
nCells <- getParams(params, "nCells")
nGenes <- getParams(params, "nGenes")
mean.shape <- getParams(params, "mean.shape")
mean.rate <- getParams(params, "mean.rate")
bcv.common <- getParams(params, "bcv.common")
if (verbose) {message("Simulating means...")}
means <- rgamma(nGenes, shape = mean.shape, rate = mean.rate)
if (verbose) {message("Simulating counts...")}
counts <- matrix(rnbinom(nGenes * nCells, mu = means,
size = 1 / bcv.common),
nrow = nGenes, ncol = nCells)
if (verbose) {message("Creating final SCESet...")}
cell.names <- paste0("Cell", 1:nCells)
gene.names <- paste0("Gene", 1:nGenes)
rownames(counts) <- gene.names
colnames(counts) <- cell.names
phenos <- new("AnnotatedDataFrame", data = data.frame(Cell = cell.names))
rownames(phenos) <- cell.names
features <- new("AnnotatedDataFrame",
data = data.frame(Gene = gene.names, GeneMean = means))
rownames(features) <- gene.names
sim <- newSCESet(countData = counts, phenoData = phenos,
featureData = features)
return(sim)
}