Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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
#' 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
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)
}