diff --git a/DESCRIPTION b/DESCRIPTION index 3ade2f07e459effa09d9f891cf6fa343f0fd7f1b..0ceb4d883186dc6c9eaaeb5c32e863bc903c0dc1 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.5.9012 -Date: 2019-09-09 +Version: 1.9.7.9012 +Date: 2019-09-19 Author: Luke Zappia Authors@R: c(person("Luke", "Zappia", role = c("aut", "cre"), @@ -58,7 +58,7 @@ Suggests: scran, mfa, phenopath, - BASiCS, + BASiCS (>= 1.7.10), zinbwave, SparseDC, BiocManager, diff --git a/NEWS.md b/NEWS.md index 5dd3f2e0807fc341c4589a5c7fcad389174444b2..48494f7754a1be38f9319db2035669da1fe75483 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,15 @@ +### Version 1.9.7.9013 (2019-09-19) + +* Merge master into splotch branch + +## Version 1.9.7 (2019-09-19) + +* Rescale when sampling Phi in BASiCSSimulate + +## Version 1.9.6 (2019-09-17) + +* Adjust BASiCSSimulate to match development version of BASiCS + ### Version 1.9.5.9012 (2019-09-09) * Merge master into splotch branch diff --git a/R/BASiCS-simulate.R b/R/BASiCS-simulate.R index f2f40da935393b8f76003990ca10de021a5a5aa0..2b5ab176af7ef79a0a7cbc05b2e796c4e4396313 100644 --- a/R/BASiCS-simulate.R +++ b/R/BASiCS-simulate.R @@ -48,6 +48,10 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { batch.cells <- getParam(params, "batchCells") 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 if (nrow(gene.params) != nGenes) { warning("Number of gene.params not equal to nGenes, ", @@ -59,17 +63,18 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { mu <- gene.params$Mean delta <- gene.params$Delta - if (nSpikes > 0) { + has.spikes <- nSpikes > 0 + spike.mu <- NULL + if (has.spikes) { spike.mu <- getParam(params, "spike.means") if (length(spike.mu) != nSpikes) { warning("Number of spike-in means not equal to nSpikes, ", "spike.means will be sampled.") selected <- sample(length(spike.mu), nSpikes, replace = TRUE) spike.mu <- spike.mu[selected] + + params <- setParam(params, "spike.mu", spike.mu) } - } else { - # Create dummy spike-ins to get around BASiCS_Sim... - spike.mu <- c(10, 10) } cell.params <- getParam(params, "cell.params") @@ -78,6 +83,10 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { "cell.params will be sampled.") selected <- sample(nrow(cell.params), nCells, replace = TRUE) cell.params <- cell.params[selected, ] + + cell.params$Phi <- (cell.params$Phi / sum(cell.params$Phi)) * nCells + + params <- setParam(params, "cell.params", cell.params) } thetas <- getParam(params, "theta") @@ -87,31 +96,32 @@ BASiCSSimulate <- function(params = newBASiCSParams(), verbose = TRUE, ...) { batches <- unlist(batches) if (verbose) {message("Simulating counts with BASiCS...")} - counts.list <- list() - 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] + if (has.spikes) { + BASiCS.sim <- suppressMessages( + BASiCS::BASiCS_Sim(Mu = mu, Mu_spikes = spike.mu, Delta = delta, + Phi = cell.params$Phi, S = cell.params$S, + Theta = thetas, BatchInfo = batches) + ) + 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, spike.mu, delta, phi, s, theta) - ) - batch.counts <- SummarizedExperiment::assay(BASiCS.sim) - counts.list[[batch]] <- batch.counts + 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...")} cell.names <- paste0("Cell", seq_len(nCells)) gene.names <- paste0("Gene", seq_len(nGenes)) - if (nSpikes > 0) { + if (has.spikes) { 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