diff --git a/DESCRIPTION b/DESCRIPTION index b23c104..dd49a82 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: DropletUtils -Version: 1.27.1 -Date: 2024-12-04 +Version: 1.27.2 +Date: 2024-12-10 Title: Utilities for Handling Single-Cell Droplet Data Authors@R: c( person("Aaron", "Lun", role = "aut"), diff --git a/R/emptyDrops.R b/R/emptyDrops.R index 65cedee..2cce499 100644 --- a/R/emptyDrops.R +++ b/R/emptyDrops.R @@ -75,16 +75,20 @@ #' The latter is a \dQuote{no frills} version that is largely intended for use within other functions. #' #' @section Handling overdispersion: +#' By default, \code{alpha=Inf} which means that the count vector is modelled with a multinomial distribution. +#' This is appropriate when molecules are independently sampled into each droplet. +#' #' If \code{alpha} is set to a positive number, sampling is assumed to follow a Dirichlet-multinomial (DM) distribution. #' The parameter vector of the DM distribution is defined as the estimated ambient profile scaled by \code{alpha}. -#' Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules. +#' Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules (e.g., aggregates, PCR duplication). #' If \code{alpha=NULL}, a maximum likelihood estimate is obtained from the count profiles for all barcodes with totals less than or equal to \code{lower}. -#' If \code{alpha=Inf}, the sampling of molecules is modelled with a multinomial distribution. #' #' Users can check whether the model is suitable by extracting the p-values for all barcodes with \code{test.ambient=TRUE}. #' Under the null hypothesis, the p-values for presumed ambient barcodes (i.e., with total counts less than or equal to \code{lower}) should be uniformly distributed. -#' Skews in the p-value distribution are indicative of an inaccuracy in the model and/or its estimates (of \code{alpha} or the ambient profile). -#' +#' Skews in the p-value distribution are indicative of an inaccuracy in the model. +#' For example, an inaccurate \code{alpha} or ambient profile will manifest in the overenrichment of low p-values. +#' Conversely, very sparse data will often exhibit in a enrichment of p-values at 1 as the Good-Turing probabilities in the ambient profile cannot be zero. +#' #' @section \code{NA} values in the results: #' We assume that barcodes with total UMI counts less than or equal to \code{lower} correspond to empty droplets. #' These are used to estimate the ambient expression profile against which the remaining barcodes are tested. @@ -181,7 +185,7 @@ NULL #' @export #' @rdname emptyDrops -testEmptyDrops <- function(m, lower=100, niters=10000, test.ambient=FALSE, ignore=NULL, alpha=NULL, round=TRUE, by.rank=NULL, known.empty=NULL, BPPARAM=SerialParam()) { +testEmptyDrops <- function(m, lower=100, niters=10000, test.ambient=FALSE, ignore=NULL, alpha=Inf, round=TRUE, by.rank=NULL, known.empty=NULL, BPPARAM=SerialParam()) { ambfun <- function(mat, totals) { assumed.empty <- .get_putative_empty(totals, lower, by.rank, known.empty) astats <- .compute_ambient_stats(mat, totals, assumed.empty) diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 32ce03e..3fd9af4 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -6,6 +6,10 @@ \item Use a more stable algorithm for identifying the knee point in \code{barcodeRanks()}. The new algorithm is based on maximizing the distance from a line between the plateau and the inflection point. Previously, we tried to minimize the signed curvature but this was susceptible to many local minima due to the instability of the empirical second derivative, even after smoothing. + +\item Set \code{alpha=Inf} as the default for \code{testEmptyDrops()}. +This is motivated by the realization that an underestimated \code{alpha} can still yield anticonservative p-values and is not universally safer than \code{alpha=Inf}. +Defaulting \code{alpha=Inf} is preferable as it is at least correct in the expected case of multinomial sampling. }} \section{Version 1.18.0}{\itemize{ diff --git a/man/emptyDrops.Rd b/man/emptyDrops.Rd index 3b8efb6..2320278 100644 --- a/man/emptyDrops.Rd +++ b/man/emptyDrops.Rd @@ -13,7 +13,7 @@ testEmptyDrops( niters = 10000, test.ambient = FALSE, ignore = NULL, - alpha = NULL, + alpha = Inf, round = TRUE, by.rank = NULL, known.empty = NULL, @@ -144,15 +144,19 @@ The latter is a \dQuote{no frills} version that is largely intended for use with \section{Handling overdispersion}{ +By default, \code{alpha=Inf} which means that the count vector is modelled with a multinomial distribution. +This is appropriate when molecules are independently sampled into each droplet. + If \code{alpha} is set to a positive number, sampling is assumed to follow a Dirichlet-multinomial (DM) distribution. The parameter vector of the DM distribution is defined as the estimated ambient profile scaled by \code{alpha}. -Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules. +Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules (e.g., aggregates, PCR duplication). If \code{alpha=NULL}, a maximum likelihood estimate is obtained from the count profiles for all barcodes with totals less than or equal to \code{lower}. -If \code{alpha=Inf}, the sampling of molecules is modelled with a multinomial distribution. Users can check whether the model is suitable by extracting the p-values for all barcodes with \code{test.ambient=TRUE}. Under the null hypothesis, the p-values for presumed ambient barcodes (i.e., with total counts less than or equal to \code{lower}) should be uniformly distributed. -Skews in the p-value distribution are indicative of an inaccuracy in the model and/or its estimates (of \code{alpha} or the ambient profile). +Skews in the p-value distribution are indicative of an inaccuracy in the model. +For example, an inaccurate \code{alpha} or ambient profile will manifest in the overenrichment of low p-values. +Conversely, very sparse data will often exhibit in a enrichment of p-values at 1 as the Good-Turing probabilities in the ambient profile cannot be zero. } \section{\code{NA} values in the results}{