From e33763c596356eb53c62c5577b403036c7cc255a Mon Sep 17 00:00:00 2001 From: Aaron Lun Date: Tue, 10 Dec 2024 14:47:49 -0800 Subject: [PATCH] Set alpha=Inf as the default for testEmptyDrops. (#118) This is because a finite alpha is not universally safer; despite representing an overdispersed multinomial, using it to compute p-values for a multinomial-distributed vector will actually be anticonservative! This means that an inaccurate alpha cannot be brushed under the carpet. Previously, I was trusting that any finite alpha would be safer than alpha=Inf, but it seems that this is not true. We can't just ignore poor estimates of alpha that we get from the assumed-ambient counts. So, we keep it simple and set the default to alpha=Inf, which is exactly correct for the multinomial case, i.e., independent sampling of molecules into droplet (which should represent most scenarios, hopefully). This also reduces runtime and aligns with the CellRanger defaults. --- DESCRIPTION | 4 ++-- R/emptyDrops.R | 14 +++++++++----- inst/NEWS.Rd | 4 ++++ man/emptyDrops.Rd | 12 ++++++++---- 4 files changed, 23 insertions(+), 11 deletions(-) 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}{