diff --git a/DESCRIPTION b/DESCRIPTION
index 888d745a0525cac9401d8c44f67fdcda85dc895b..b036793a3ce885f45e7477baa4ccd107a847921f 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.4.9008
-Date: 2019-08-20
+Version: 1.9.4.9009
+Date: 2019-08-21
Author: Luke Zappia
Authors@R:
c(person("Luke", "Zappia", role = c("aut", "cre"),
diff --git a/NAMESPACE b/NAMESPACE
index a2a713a2320e693e1eaa6378f3561d46e76ce219..08a786f5248056391624209b122e8d5908e99d35 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -114,6 +114,7 @@ importFrom(ggplot2,geom_hline)
importFrom(ggplot2,geom_point)
importFrom(ggplot2,geom_smooth)
importFrom(ggplot2,geom_tile)
+importFrom(ggplot2,geom_violin)
importFrom(ggplot2,ggplot)
importFrom(ggplot2,ggtitle)
importFrom(ggplot2,scale_colour_manual)
@@ -136,8 +137,10 @@ importFrom(methods,slot)
importFrom(methods,slotNames)
importFrom(methods,validObject)
importFrom(stats,aggregate)
+importFrom(stats,approxfun)
importFrom(stats,cor)
importFrom(stats,dbeta)
+importFrom(stats,density)
importFrom(stats,dnbinom)
importFrom(stats,ks.test)
importFrom(stats,median)
diff --git a/NEWS.md b/NEWS.md
index 9143752a6ffda540857eb8d74ee4a60a3e328356..2aa1416d811677f2265940608b57d26d644f417e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+### Version 1.9.4.9009 (2019-08-21)
+
+* Replace library size log-normal with density and rejection sampling
+
### Version 1.9.4.9008 (2019-08-20)
* Merge master into splotch branch
diff --git a/R/AllClasses.R b/R/AllClasses.R
index 4124ac3fd8cfc0012ba40b3ba153f362cc675de4..3af99feb6aa582106fa15c63b5b17496ee20d183 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -304,6 +304,8 @@ setClass("SplatParams",
#' \item{\code{lib.scale}}{Scale (sdlog) parameter for the library
#' size log-normal distribution, or sd parameter if a normal
#' distribution is used.}
+#' \item{\code{lib.dens}}{\code{\link{density}} object describing
+#' the library size density.}
#' }
#' }
#' \item{\emph{Paths parameters}}{
@@ -338,6 +340,7 @@ setClass("SplotchParams",
paths.means = "list",
lib.loc = "numeric",
lib.scale = "numeric",
+ lib.dens = "density",
cells.design = "data.frame"),
prototype = prototype(mean.rate = 0.3,
mean.shape = 0.6,
@@ -357,6 +360,7 @@ setClass("SplotchParams",
paths.means = list(),
lib.loc = 11,
lib.scale = 0.2,
+ lib.dens = density(rlnorm(10000, 11, 0.2)),
cells.design = data.frame(
Path = 1,
Probability = 1,
diff --git a/R/SplotchParams-methods.R b/R/SplotchParams-methods.R
index 9b6e736ad000c7f73be0c2b93484eb5ed98127cc..6f4d13b4751f21a8c4eef2f0563feb5234821aea 100644
--- a/R/SplotchParams-methods.R
+++ b/R/SplotchParams-methods.R
@@ -37,8 +37,9 @@ setValidity("SplotchParams", function(object) {
any.missing = FALSE,
min.rows = 1,
ncols = 3),
- lib.loc = checkNumber(v$lib.loc),
- lib.scale = checkNumber(v$lib.scale, lower = 0),
+ lib.loc = checkmate::checkNumber(v$lib.loc),
+ lib.scale = checkmate::checkNumber(v$lib.scale, lower = 0),
+ lib.dens = checkmate::checkClass(v$lib.dens, "density"),
cells.design = checkmate::checkDataFrame(v$cells.design,
types = "numeric",
any.missing = FALSE,
@@ -133,8 +134,9 @@ setMethod("show", "SplotchParams", function(object) {
"[Design]" = "paths.design"))
pp.bot <- list("Library size:" = c("(Location)" = "lib.loc",
- "(Scale)" = "lib.scale"),
- "Cells:" = c("[Design]" = "cells.design"))
+ "(Scale)" = "lib.scale",
+ "(Density)" = "lib.dens"),
+ "Cells:" = c("[Design]" = "cells.design"))
paths.means <- getParam(object, "paths.means")
if (length(paths.means) == 0) {
diff --git a/R/params-functions.R b/R/params-functions.R
index 5e7eed904a0e1ca5c272f17cdb3c7a010a49ad9e..3a05a05a994219f33dfa6783951b8c79a6a0a723 100644
--- a/R/params-functions.R
+++ b/R/params-functions.R
@@ -107,10 +107,19 @@ showValues <- function(values, not.default) {
len = length(values))
short.values <- sapply(values, function(x) {
- if (length(x) > 4) {
- paste0(paste(head(x, n = 4), collapse = ", "), ",...")
+ if (is.list(x)) {
+ classes <- class(x)
+ if (length(classes) == 1 && classes == "list") {
+ paste("List with", length(x), "items")
+ } else {
+ paste("Object of class", paste(classes, collapse = ", "))
+ }
} else {
- paste(x, collapse = ", ")
+ if (length(x) > 4) {
+ paste0(paste(head(x, n = 4), collapse = ", "), ",...")
+ } else {
+ paste(x, collapse = ", ")
+ }
}
})
diff --git a/R/splotch-estimate.R b/R/splotch-estimate.R
index d4bff448926cf4824e4a855ea65aae6979c314bc..523e07b44c2cc1a75b8f9550ee41d08d15e31cff 100644
--- a/R/splotch-estimate.R
+++ b/R/splotch-estimate.R
@@ -98,6 +98,7 @@ splotchEstMean <- function(norm.counts, params, verbose) {
return(params)
}
+#' @importFrom stats density
splotchEstLib <- function(counts, params, verbose) {
if (verbose) {message("Estimating library size parameters...")}
@@ -109,7 +110,10 @@ splotchEstLib <- function(counts, params, verbose) {
lib.loc <- unname(fit$estimate["meanlog"])
lib.scale <- unname(fit$estimate["sdlog"])
- params <- setParams(params, lib.loc = lib.loc, lib.scale = lib.scale)
+ dens <- density(lib.sizes)
+
+ params <- setParams(params, lib.loc = lib.loc, lib.scale = lib.scale,
+ lib.dens = dens)
return(params)
}
diff --git a/R/splotch-simulate.R b/R/splotch-simulate.R
index 3fa1e829eab170ebf04e993679c0d0f3eb6ca44f..1148c61d08f06a05a53ae88efdba85f2a16b7eec 100644
--- a/R/splotch-simulate.R
+++ b/R/splotch-simulate.R
@@ -242,10 +242,13 @@ splotchSimLibSizes <- function(sim, params, verbose) {
if (verbose) {message("Simulating library sizes...")}
nCells <- getParam(params, "nCells")
- lib.loc <- getParam(params, "lib.loc")
- lib.scale <- getParam(params, "lib.scale")
+ # lib.loc <- getParam(params, "lib.loc")
+ # lib.scale <- getParam(params, "lib.scale")
+ lib.dens <- getParam(params, "lib.dens")
+
+ # exp.lib.sizes <- rlnorm(nCells, lib.loc, lib.scale)
+ exp.lib.sizes <- rejectionSample(nCells, lib.dens)
- exp.lib.sizes <- rlnorm(nCells, lib.loc, lib.scale)
colData(sim)$ExpLibSize <- exp.lib.sizes
return(sim)
@@ -347,3 +350,30 @@ getBetaStepProbs <- function(steps, alpha, beta) {
return(probs)
}
+
+#' @importFrom stats approxfun
+rejectionSample <- function(n, dens, lower = 0) {
+
+ xmin <- min(dens$x)
+ xmax <- max(dens$x)
+ ymin <- min(dens$y)
+ ymax <- max(dens$y)
+
+ boundary <- approxfun(dens$x, dens$y)
+
+ values <- c()
+ nsel <- 0
+
+ while(nsel < n) {
+ x <- runif(1e4, xmin, xmax)
+ y <- runif(1e4, ymin, ymax)
+ sel <- y < boundary(x) & x > lower
+
+ nsel <- nsel + sum(sel)
+ values <- c(values, x[sel])
+ }
+
+ values <- values[1:n]
+
+ return(values)
+}
diff --git a/man/SplotchParams.Rd b/man/SplotchParams.Rd
index 904b42ad5703103a0ef480cc44535ce55343a673..2c7861c66e46740d9767144eb84c5bcb21bc8523 100644
--- a/man/SplotchParams.Rd
+++ b/man/SplotchParams.Rd
@@ -54,6 +54,8 @@ The Splotch simulation uses the following parameters:
\item{\code{lib.scale}}{Scale (sdlog) parameter for the library
size log-normal distribution, or sd parameter if a normal
distribution is used.}
+ \item{\code{lib.dens}}{\code{\link{density}} object describing
+ the library size density.}
}
}
\item{\emph{Paths parameters}}{