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

Adjust BASiCSSimulate to match devel BASiCS

parent f9607578
No related branches found
No related tags found
No related merge requests found
Package: splatter Package: splatter
Type: Package Type: Package
Title: Simple Simulation of Single-cell RNA Sequencing Data Title: Simple Simulation of Single-cell RNA Sequencing Data
Version: 1.9.5 Version: 1.9.6
Date: 2019-09-09 Date: 2019-09-17
Author: Luke Zappia Author: Luke Zappia
Authors@R: Authors@R:
c(person("Luke", "Zappia", role = c("aut", "cre"), c(person("Luke", "Zappia", role = c("aut", "cre"),
...@@ -58,7 +58,7 @@ Suggests: ...@@ -58,7 +58,7 @@ Suggests:
scran, scran,
mfa, mfa,
phenopath, phenopath,
BASiCS, BASiCS (>= 1.7.10),
zinbwave, zinbwave,
SparseDC, SparseDC,
BiocManager, BiocManager,
......
## Version 1.9.6 (2019-09-17)
* Adjust BASiCSSimulate to match development version of BASiCS
## Version 1.9.5 (2019-09-09) ## Version 1.9.5 (2019-09-09)
* Fix bug in NAMESPACE * Fix bug in NAMESPACE
......
...@@ -48,6 +48,10 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { ...@@ -48,6 +48,10 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) {
batch.cells <- getParam(params, "batchCells") batch.cells <- getParam(params, "batchCells")
gene.params <- getParam(params, "gene.params") gene.params <- getParam(params, "gene.params")
if (nSpikes == 0 && nBatches == 1) {
stop("If there are no spikes there must be multiple batches")
}
# Sample gene.params if necessary # Sample gene.params if necessary
if (nrow(gene.params) != nGenes) { if (nrow(gene.params) != nGenes) {
warning("Number of gene.params not equal to nGenes, ", warning("Number of gene.params not equal to nGenes, ",
...@@ -59,7 +63,9 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { ...@@ -59,7 +63,9 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) {
mu <- gene.params$Mean mu <- gene.params$Mean
delta <- gene.params$Delta delta <- gene.params$Delta
if (nSpikes > 0) { has.spikes <- nSpikes > 0
spike.mu <- NULL
if (has.spikes) {
spike.mu <- getParam(params, "spike.means") spike.mu <- getParam(params, "spike.means")
if (length(spike.mu) != nSpikes) { if (length(spike.mu) != nSpikes) {
warning("Number of spike-in means not equal to nSpikes, ", warning("Number of spike-in means not equal to nSpikes, ",
...@@ -67,9 +73,6 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { ...@@ -67,9 +73,6 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) {
selected <- sample(length(spike.mu), nSpikes, replace = TRUE) selected <- sample(length(spike.mu), nSpikes, replace = TRUE)
spike.mu <- spike.mu[selected] spike.mu <- spike.mu[selected]
} }
} else {
# Create dummy spike-ins to get around BASiCS_Sim...
spike.mu <- c(10, 10)
} }
cell.params <- getParam(params, "cell.params") cell.params <- getParam(params, "cell.params")
...@@ -87,31 +90,32 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { ...@@ -87,31 +90,32 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) {
batches <- unlist(batches) batches <- unlist(batches)
if (verbose) {message("Simulating counts with BASiCS...")} if (verbose) {message("Simulating counts with BASiCS...")}
counts.list <- list() if (has.spikes) {
for (batch in seq_len(nBatches)) {
batch.cells <- batches == batch
phi <- cell.params[batch.cells, "Phi"]
phi <- (phi / sum(phi)) * sum(batch.cells)
s <- cell.params[batch.cells, "S"]
theta <- thetas[batch]
BASiCS.sim <- suppressMessages( BASiCS.sim <- suppressMessages(
BASiCS::BASiCS_Sim(mu, spike.mu, delta, phi, s, theta) BASiCS::BASiCS_Sim(Mu = mu, Mu_spikes = spike.mu, Delta = delta,
) Phi = cell.params$Phi, S = cell.params$S,
batch.counts <- SummarizedExperiment::assay(BASiCS.sim) Theta = thetas, BatchInfo = batches)
counts.list[[batch]] <- batch.counts )
counts <- SingleCellExperiment::counts(BASiCS.sim)
spike.counts <- SummarizedExperiment::assay(
SingleCellExperiment::altExp(BASiCS.sim)
)
counts <- rbind(counts, spike.counts)
} else {
if (verbose) {message("Simulating without spikes")}
BASiCS.sim <- suppressMessages(
BASiCS::BASiCS_Sim(Mu = mu, Mu_spikes = NULL, Delta = delta,
Phi = NULL, S = cell.params$S,
Theta = thetas, BatchInfo = batches)
)
counts <- SingleCellExperiment::counts(BASiCS.sim)
} }
counts <- do.call(cbind, counts.list)
if (verbose) {message("Creating final dataset...")} if (verbose) {message("Creating final dataset...")}
cell.names <- paste0("Cell", seq_len(nCells)) cell.names <- paste0("Cell", seq_len(nCells))
gene.names <- paste0("Gene", seq_len(nGenes)) gene.names <- paste0("Gene", seq_len(nGenes))
if (nSpikes > 0) { if (has.spikes) {
gene.names <- c(gene.names, paste0("Spike", seq_len(nSpikes))) gene.names <- c(gene.names, paste0("Spike", seq_len(nSpikes)))
} else {
# Remove dummy spikes
counts <- counts[1:(nrow(counts) - 2), ]
spike.mu <- numeric()
} }
rownames(counts) <- gene.names rownames(counts) <- gene.names
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment