Skip to content
Snippets Groups Projects
Commit 7d1fcd34 authored by Luke Zappia's avatar Luke Zappia
Browse files

Add density sampling option for means, lib sizes

parent 84cdda8f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 2
− 2
View file @ 7d1fcd34
Package: splatter
Type: Package
Title: Simple Simulation of Single-cell RNA Sequencing Data
Version: 1.9.4.9009
Date: 2019-08-21
Version: 1.9.4.9010
Date: 2019-08-22
Author: Luke Zappia
Authors@R:
c(person("Luke", "Zappia", role = c("aut", "cre"),
......
NEWS.md
+ 4
− 0
View file @ 7d1fcd34
### Version 1.9.4.9010 (2019-08-22)
* Add density sampling options for means and library sizes
### Version 1.9.4.9009 (2019-08-21)
* Replace library size log-normal with density and rejection sampling
......
R/AllClasses.R
+ 17
− 1
View file @ 7d1fcd34
......@@ -279,7 +279,13 @@ setClass("SplatParams",
#' the expression outlier factor log-normal distribution.}
#' \item{\code{mean.outFacScale}}{Scale (sdlog) parameter for the
#' expression outlier factor log-normal distribution.}
#' \item{\code{mean.values}}{Vector of means for each gene.}
#' \item{\code{mean.dens}}{\code{\link{density}} object describing
#' the log gene mean density.}
#' \item{\code{[mean.method]}}{Method to use for simulating gene
#' means. Either "fit" to sample from a gamma distribution (with
#' expression outliers) or "density" to sample from the provided
#' density object.}
#' \item{\code{[mean.values]}}{Vector of means for each gene.}
#' }
#' }
#' \item{\emph{Network parameters}}{
......@@ -306,6 +312,9 @@ setClass("SplatParams",
#' distribution is used.}
#' \item{\code{lib.dens}}{\code{\link{density}} object describing
#' the library size density.}
#' \item{\code{[lib.method]}}{Method to use for simulating library
#' sizes. Either "fit" to sample from a log-normal distribution or
#' "density" to sample from the provided density object.}
#' }
#' }
#' \item{\emph{Paths parameters}}{
......@@ -331,6 +340,8 @@ setClass("SplotchParams",
mean.outProb = "numeric",
mean.outLoc = "numeric",
mean.outScale = "numeric",
mean.dens = "density",
mean.method = "character",
mean.values = "numeric",
network.graph = "ANY",
network.nRegs = "numeric",
......@@ -341,12 +352,16 @@ setClass("SplotchParams",
lib.loc = "numeric",
lib.scale = "numeric",
lib.dens = "density",
lib.method = "character",
cells.design = "data.frame"),
prototype = prototype(mean.rate = 0.3,
mean.shape = 0.6,
mean.outProb = 0.05,
mean.outLoc = 4,
mean.outScale = 0.5,
mean.dens = density(rgamma(10000, rate = 0.3,
shape = 0.6)),
mean.method = "fit",
mean.values = numeric(),
network.graph = NULL,
network.nRegs = 100,
......@@ -361,6 +376,7 @@ setClass("SplotchParams",
lib.loc = 11,
lib.scale = 0.2,
lib.dens = density(rlnorm(10000, 11, 0.2)),
lib.method = "fit",
cells.design = data.frame(
Path = 1,
Probability = 1,
......
R/SplotchParams-methods.R
+ 15
− 5
View file @ 7d1fcd34
......@@ -25,6 +25,9 @@ setValidity("SplotchParams", function(object) {
mean.outProb = checkNumber(v$mean.outProb, lower = 0, upper = 1),
mean.outLoc = checkNumber(v$mean.outLoc),
mean.outScale = checkNumber(v$mean.outScale, lower = 0),
mean.dens = checkmate::checkClass(v$mean.dens, "density"),
mean.method = checkmate::checkChoice(v$mean.method,
c("fit", "density")),
network.graph = checkmate::checkClass(v$network, "igraph",
null.ok = TRUE),
network.nRegs = checkmate::checkInt(v$network.nRegs,
......@@ -40,6 +43,8 @@ setValidity("SplotchParams", function(object) {
lib.loc = checkmate::checkNumber(v$lib.loc),
lib.scale = checkmate::checkNumber(v$lib.scale, lower = 0),
lib.dens = checkmate::checkClass(v$lib.dens, "density"),
lib.method = checkmate::checkChoice(v$lib.method,
c("fit", "density")),
cells.design = checkmate::checkDataFrame(v$cells.design,
types = "numeric",
any.missing = FALSE,
......@@ -124,6 +129,8 @@ setMethod("show", "SplotchParams", function(object) {
"(Out Prob)" = "mean.outProb",
"(Out Location)" = "mean.outLoc",
"(Out Scale)" = "mean.outScale",
"(Density)" = "mean.density",
"[Method]" = "mean.method",
"[Values]*" = "mean.values"))
pp.network <- list("Network:" = c("[Graph]" = "network.graph",
......@@ -135,7 +142,8 @@ setMethod("show", "SplotchParams", function(object) {
pp.bot <- list("Library size:" = c("(Location)" = "lib.loc",
"(Scale)" = "lib.scale",
"(Density)" = "lib.dens"),
"(Density)" = "lib.dens",
"[Method]" = "lib.method"),
"Cells:" = c("[Design]" = "cells.design"))
paths.means <- getParam(object, "paths.means")
......@@ -205,10 +213,12 @@ setMethod("setParam", "SplotchParams", function(object, name, value) {
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 (getParam(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())
}
object <- setParamUnchecked(object, "nGenes", igraph::gorder(value))
}
if ("IsReg" %in% igraph::vertex_attr_names(value)) {
......
R/splotch-estimate.R
+ 6
− 2
View file @ 7d1fcd34
......@@ -78,7 +78,7 @@ splotchEstMean <- function(norm.counts, params, verbose) {
med <- median(lmeans)
mad <- mad(lmeans)
bound <- med + 1 * mad
bound <- med + 2 * mad
outs <- which(lmeans > bound)
......@@ -88,13 +88,16 @@ splotchEstMean <- function(norm.counts, params, verbose) {
if (length(outs) > 1) {
facs <- means[outs] / median(means)
fit <- selectFit(facs, "lnorm", verbose = verbose)
#fit <- selectFit(facs, "lnorm", verbose = verbose)
fit <- fitdistrplus::fitdist(facs, "lnorm")
params <- setParams(params,
mean.outLoc = unname(fit$estimate["meanlog"]),
mean.outScale = unname(fit$estimate["sdlog"]))
}
params <- setParams(params, mean.dens = density(lmeans))
return(params)
}
......@@ -191,6 +194,7 @@ selectFit <- function(data, distr, weights = NULL, verbose = TRUE) {
}
scores <- fitdistrplus::gofstat(fits)$cvm
# Flatten in case scores is a list
scores.flat <- unlist(scores)
selected <- which(scores.flat == min(scores.flat, na.rm = TRUE))
......
R/splotch-simulate.R
+ 38
− 20
View file @ 7d1fcd34
......@@ -146,20 +146,30 @@ splotchSimGeneMeans <- function(params, verbose) {
if (verbose) {message("Simulating means...")}
nGenes <- getParam(params, "nGenes")
mean.shape <- getParam(params, "mean.shape")
mean.rate <- getParam(params, "mean.rate")
mean.outProb <- getParam(params, "mean.outProb")
mean.outLoc <- getParam(params, "mean.outLoc")
mean.outScale <- getParam(params, "mean.outScale")
mean.values <- rgamma(nGenes, shape = mean.shape, rate = mean.rate)
outlier.facs <- getLNormFactors(nGenes, mean.outProb, 0, mean.outLoc,
mean.outScale)
median.means.gene <- median(mean.values)
outlier.means <- median.means.gene * outlier.facs
is.outlier <- outlier.facs != 1
mean.values[is.outlier] <- outlier.means[is.outlier]
mean.method <- getParam(params, "mean.method")
if (mean.method == "fit") {
if (verbose) {message("Sampling from gamma distribution...")}
mean.shape <- getParam(params, "mean.shape")
mean.rate <- getParam(params, "mean.rate")
mean.outProb <- getParam(params, "mean.outProb")
mean.outLoc <- getParam(params, "mean.outLoc")
mean.outScale <- getParam(params, "mean.outScale")
mean.values <- rgamma(nGenes, shape = mean.shape, rate = mean.rate)
outlier.facs <- getLNormFactors(nGenes, mean.outProb, 0, mean.outLoc,
mean.outScale)
median.means.gene <- median(mean.values)
outlier.means <- median.means.gene * outlier.facs
is.outlier <- outlier.facs != 1
mean.values[is.outlier] <- outlier.means[is.outlier]
} else if (mean.method == "density") {
if (verbose) {message("Sampling from density object...")}
mean.dens <- getParam(params, "mean.dens")
mean.values <- exp(sampleDensity(nGenes, mean.dens, lower = -Inf))
}
params <- setParam(params, "mean.values", mean.values)
......@@ -242,12 +252,20 @@ 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.dens <- getParam(params, "lib.dens")
lib.method <- getParam(params, "lib.method")
# exp.lib.sizes <- rlnorm(nCells, lib.loc, lib.scale)
exp.lib.sizes <- rejectionSample(nCells, lib.dens)
if (lib.method == "fit") {
if (verbose) {message("Sampling from log-normal distribution...")}
lib.loc <- getParam(params, "lib.loc")
lib.scale <- getParam(params, "lib.scale")
exp.lib.sizes <- rlnorm(nCells, lib.loc, lib.scale)
} else if (lib.method == "density") {
if (verbose) {message("Sampling from density object...")}
lib.dens <- getParam(params, "lib.dens")
exp.lib.sizes <- sampleDensity(nCells, lib.dens)
}
colData(sim)$ExpLibSize <- exp.lib.sizes
......@@ -352,7 +370,7 @@ getBetaStepProbs <- function(steps, alpha, beta) {
}
#' @importFrom stats approxfun
rejectionSample <- function(n, dens, lower = 0) {
sampleDensity <- function(n, dens, lower = 0) {
xmin <- min(dens$x)
xmax <- max(dens$x)
......
man/SplotchParams.Rd
+ 10
− 1
View file @ 7d1fcd34
......@@ -29,7 +29,13 @@ The Splotch simulation uses the following parameters:
the expression outlier factor log-normal distribution.}
\item{\code{mean.outFacScale}}{Scale (sdlog) parameter for the
expression outlier factor log-normal distribution.}
\item{\code{mean.values}}{Vector of means for each gene.}
\item{\code{mean.dens}}{\code{\link{density}} object describing
the log gene mean density.}
\item{\code{[mean.method]}}{Method to use for simulating gene
means. Either "fit" to sample from a gamma distribution (with
expression outliers) or "density" to sample from the provided
density object.}
\item{\code{[mean.values]}}{Vector of means for each gene.}
}
}
\item{\emph{Network parameters}}{
......@@ -56,6 +62,9 @@ The Splotch simulation uses the following parameters:
distribution is used.}
\item{\code{lib.dens}}{\code{\link{density}} object describing
the library size density.}
\item{\code{[lib.method]}}{Method to use for simulating library
sizes. Either "fit" to sample from a log-normal distribution or
"density" to sample from the provided density object.}
}
}
\item{\emph{Paths parameters}}{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment