From 321c2d0cf42d7fbed71f9988b2ebcb25d071bedf Mon Sep 17 00:00:00 2001 From: Luke Zappia <lazappi@users.noreply.github.com> Date: Thu, 5 Oct 2017 11:31:07 +1100 Subject: [PATCH] Add zinbEstimate --- NAMESPACE | 3 ++ R/ZINBParams-methods.R | 4 +- R/zinb-estimate.R | 95 ++++++++++++++++++++++++++++++++++++++++++ man/zinbEstimate.Rd | 67 +++++++++++++++++++++++++++++ 4 files changed, 167 insertions(+), 2 deletions(-) create mode 100644 R/zinb-estimate.R create mode 100644 man/zinbEstimate.Rd diff --git a/NAMESPACE b/NAMESPACE index 7ecb378..4cfad39 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 91b6863..2cddd4d 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 0000000..217283c --- /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 0000000..84956c9 --- /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 +} -- GitLab