diff --git a/NAMESPACE b/NAMESPACE
index 7ecb378b092160b094f738ddfa664f6df0af511c..4cfad39e5ef56f527991f6f4713875addfdcd440 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,6 +15,8 @@ S3method(simpleEstimate,SingleCellExperiment)
S3method(simpleEstimate,matrix)
S3method(splatEstimate,SingleCellExperiment)
S3method(splatEstimate,matrix)
+S3method(zinbEstimate,SingleCellExperiment)
+S3method(zinbEstimate,matrix)
export(addGeneLengths)
export(compareSCEs)
export(diffSCEs)
@@ -52,6 +54,7 @@ export(splatSimulateGroups)
export(splatSimulatePaths)
export(splatSimulateSingle)
export(summariseDiff)
+export(zinbEstimate)
export(zinbSimulate)
exportClasses(Lun2Params)
exportClasses(LunParams)
diff --git a/R/ZINBParams-methods.R b/R/ZINBParams-methods.R
index 91b6863787f7b0d54fdd3eb6b2157984206cbbcb..2cddd4da2f4d507dff03f5894ea073989add3679 100644
--- a/R/ZINBParams-methods.R
+++ b/R/ZINBParams-methods.R
@@ -49,7 +49,7 @@ setMethod("setParam", "ZINBParams", function(object, name, value) {
zinbwave::nSamples(value))
}
- object <- cellNextMethod()
+ object <- callNextMethod()
return(object)
})
@@ -90,7 +90,7 @@ setMethod("show", "ZINBParams", function(object) {
cat("Model:", "\n")
cat("ZinbModel with", zinbwave::nFeatures(model), "features,",
zinbwave::nSamples(model), "samples,", zinbwave::nFactors(model),
- "factors and", zinbwave::nParams(model), "parameters", "\n\n")
+ "latent factors and", zinbwave::nParams(model), "parameters", "\n\n")
default <- zinbwave::zinbModel()
for (category in names(pp)) {
diff --git a/R/zinb-estimate.R b/R/zinb-estimate.R
new file mode 100644
index 0000000000000000000000000000000000000000..217283c8e91058b5199ee67462b5437b0ecbce0b
--- /dev/null
+++ b/R/zinb-estimate.R
@@ -0,0 +1,95 @@
+#' Estimate ZINB-WaVE simulation parameters
+#'
+#' Estimate simulation parameters for the ZINB-WaVE simulation from a real
+#' dataset.
+#'
+#' @param counts either a counts matrix or a SingleCellExperiment object
+#' containing count data to estimate parameters from.
+#' @param design.samples design matrix of sample-level covariates.
+#' @param design.genes design matrix of gene-level covariates.
+#' @param common.disp logical. Whether or not a single dispersion for all
+#' features is estimated.
+#' @param iter.init number of iterations to use for initalization.
+#' @param iter.opt number of iterations to use for optimization.
+#' @param stop.opt stopping criterion for optimization.
+#' @param params ZINBParams object to store estimated values in.
+#' @param verbose logical. Whether to print progress messages.
+#' @param BPPARAM A \code{\link[BiocParallel]{BiocParallelParam}} instance
+#' giving the parallel back-end to be used. Default is
+#' \code{\link[BiocParallel]{SerialParam}} which uses a single core.
+#' @param ... additional arguments passes to \code{\link[zinbwave]{zinbFit}}.
+#'
+#' @details
+#' The function is a wrapper around \code{\link[zinbwave]{zinbFit}} that takes
+#' the fitted model and inserts it into a \code{\link{ZINBParams}} object. See
+#' \code{\link{ZINBParams}} for more details on the parameters and
+#' \code{\link[zinbwave]{zinbFit}} for details of the estimation procedure.
+#'
+#' @return ZINBParams object containing the estimated parameters.
+#'
+#' @examples
+#' data("sc_example_counts")
+#' params <- zinbEstimate(sc_example_counts)
+#' params
+#' @importFrom BiocParallel SerialParam
+#' @export
+zinbEstimate <- function(counts, design.samples = NULL, design.genes = NULL,
+ common.disp = TRUE, iter.init = 2, iter.opt = 25,
+ stop.opt = 1e-04, params = newZINBParams(),
+ verbose = TRUE, BPPARAM = SerialParam(), ...) {
+ UseMethod("zinbEstimate")
+}
+
+#' @rdname zinbEstimate
+#' @export
+zinbEstimate.SingleCellExperiment <- function(counts, design.samples = NULL,
+ design.genes = NULL,
+ common.disp = TRUE,
+ iter.init = 2, iter.opt = 25,
+ stop.opt = 1e-04,
+ params = newZINBParams(),
+ verbose = TRUE,
+ BPPARAM = SerialParam(), ...) {
+ counts <- BiocGenerics::counts(counts)
+ zinbEstimate(counts, design.samples, design.genes, common.disp,
+ iter.init, iter.opt, stop.opt, params, verbose, BPPARAM, ...)
+}
+
+#' @rdname zinbEstimate
+#' @export
+zinbEstimate.matrix <- function(counts, design.samples = NULL,
+ design.genes = NULL, common.disp = TRUE,
+ iter.init = 2, iter.opt = 25, stop.opt = 1e-04,
+ params = newZINBParams(), verbose = TRUE,
+ BPPARAM = SerialParam(), ...) {
+
+ checkmate::assertClass(params, "ZINBParams")
+
+ if (verbose) {message("Removing all zero genes...")}
+ counts <- counts[rowSums(counts) > 0, ]
+
+ args.list <- list(Y = counts,
+ commondispersion = common.disp,
+ verbose = verbose,
+ nb.repeat.initialize = iter.init,
+ maxiter.optimize = iter.opt,
+ stop.epsilon.optimize = stop.opt,
+ BPPARAM = BPPARAM)
+
+ if (!is.null(design.samples)) {
+ args.list$X <- design.samples
+ }
+
+ if (!is.null(design.genes)) {
+ args.list$V <- design.genes
+ }
+
+ args.list <- c(args.list, list(...))
+
+ if (verbose) {message("Fitting model...")}
+ model <- do.call(zinbwave::zinbFit, args.list)
+
+ params <- setParams(params, model = model)
+
+ return(params)
+}
diff --git a/man/zinbEstimate.Rd b/man/zinbEstimate.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..84956c9d532364f5793ff5e530328f17abcc2c7c
--- /dev/null
+++ b/man/zinbEstimate.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/zinb-estimate.R
+\name{zinbEstimate}
+\alias{zinbEstimate}
+\alias{zinbEstimate.SingleCellExperiment}
+\alias{zinbEstimate.matrix}
+\title{Estimate ZINB-WaVE simulation parameters}
+\usage{
+zinbEstimate(counts, design.samples = NULL, design.genes = NULL,
+ common.disp = TRUE, iter.init = 2, iter.opt = 25, stop.opt = 1e-04,
+ params = newZINBParams(), verbose = TRUE, BPPARAM = SerialParam(), ...)
+
+\method{zinbEstimate}{SingleCellExperiment}(counts, design.samples = NULL,
+ design.genes = NULL, common.disp = TRUE, iter.init = 2, iter.opt = 25,
+ stop.opt = 1e-04, params = newZINBParams(), verbose = TRUE,
+ BPPARAM = SerialParam(), ...)
+
+\method{zinbEstimate}{matrix}(counts, design.samples = NULL,
+ design.genes = NULL, common.disp = TRUE, iter.init = 2, iter.opt = 25,
+ stop.opt = 1e-04, params = newZINBParams(), verbose = TRUE,
+ BPPARAM = SerialParam(), ...)
+}
+\arguments{
+\item{counts}{either a counts matrix or a SingleCellExperiment object
+containing count data to estimate parameters from.}
+
+\item{design.samples}{design matrix of sample-level covariates.}
+
+\item{design.genes}{design matrix of gene-level covariates.}
+
+\item{common.disp}{logical. Whether or not a single dispersion for all
+features is estimated.}
+
+\item{iter.init}{number of iterations to use for initalization.}
+
+\item{iter.opt}{number of iterations to use for optimization.}
+
+\item{stop.opt}{stopping criterion for optimization.}
+
+\item{params}{ZINBParams object to store estimated values in.}
+
+\item{verbose}{logical. Whether to print progress messages.}
+
+\item{BPPARAM}{A \code{\link[BiocParallel]{BiocParallelParam}} instance
+giving the parallel back-end to be used. Default is
+\code{\link[BiocParallel]{SerialParam}} which uses a single core.}
+
+\item{...}{additional arguments passes to \code{\link[zinbwave]{zinbFit}}.}
+}
+\value{
+ZINBParams object containing the estimated parameters.
+}
+\description{
+Estimate simulation parameters for the ZINB-WaVE simulation from a real
+dataset.
+}
+\details{
+The function is a wrapper around \code{\link[zinbwave]{zinbFit}} that takes
+the fitted model and inserts it into a \code{\link{ZINBParams}} object. See
+\code{\link{ZINBParams}} for more details on the parameters and
+\code{\link[zinbwave]{zinbFit}} for details of the estimation procedure.
+}
+\examples{
+data("sc_example_counts")
+params <- zinbEstimate(sc_example_counts)
+params
+}