Package {presto}

Title:	Fast Functions for Differential Expression Using Wilcoxon and AUC
Version:	1.1.0
Description:	Scalable implementation of the Wilcoxon rank sum test and the area under the receiver operating characteristic curve (also known as auROC or AUC) statistic. Interfaces to dense and sparse matrices, as well as the genomics analysis frameworks 'Seurat' and 'SingleCellExperiment'.
License:	GPL-3
Encoding:	UTF-8
Depends:	R (≥ 3.5.0)
LinkingTo:	Rcpp, RcppArmadillo
Imports:	Rcpp, data.table, dplyr, methods, tidyr, purrr, tibble, Matrix, rlang, stats, utils
RoxygenNote:	7.3.3
Suggests:	knitr, rmarkdown, testthat, Seurat, SingleCellExperiment, SummarizedExperiment, DelayedArray, broom, BiocStyle, DESeq2
NeedsCompilation:	yes
VignetteBuilder:	knitr
Packaged:	2026-09-20 20:34:12 UTC; ks38
Author:	Ilya Korsunsky [aut], Aparna Nathan [aut], Nghia Millard [aut], Soumya Raychaudhuri [aut], Kamil Slowikowski [cre, ctb], Austin Hartman [ctb]
Maintainer:	Kamil Slowikowski <kslowikowski@gmail.com>
Repository:	CRAN
Date/Publication:	2026-09-30 09:10:02 UTC

presto: Fast differential expression

Description

Scalable implementation of the Wilcoxon rank sum test and the area under the receiver operating characteristic curve (also known as auROC or AUC) statistic. Interfaces to dense and sparse matrices, as well as the genomics analysis frameworks 'Seurat' and 'SingleCellExperiment'.

Author(s)

Maintainer: Kamil Slowikowski kslowikowski@gmail.com (ORCID) [contributor]

Authors:

• Ilya Korsunsky ilya.korsunsky@gmail.com (ORCID)

• Aparna Nathan

• Nghia Millard (ORCID)

• Soumya Raychaudhuri (ORCID)

Other contributors:

• Austin Hartman (ORCID) [contributor]

Pipe operator

Description

Pipe operator

Usage

lhs %>% rhs

Value

return value of rhs function.

Examples

x <- 5 %>% sum(10)

Collapse a single-cell count matrix into pseudobulks

Description

Sums (or averages) the columns of a feature-by-cell count matrix according to one or more cell-metadata columns. The result is a feature-by-pseudobulk matrix where each column pools all cells that share the same combination of metadata values (e.g. all cells from a given donor in a given cluster). The resulting matrix is suitable for bulk-RNA-seq tools such as DESeq2, edgeR, or limma. See the pseudobulk vignette for an end-to-end walkthrough.

Usage

collapse_counts(
  counts_mat,
  meta_data,
  varnames,
  min_cells_per_group = 0,
  keep_n = FALSE,
  how = c("sum", "mean")[1]
)

Arguments

Value

A list with two elements:

• counts_mat - feature-by-pseudobulk numeric matrix.

• meta_data - data.frame with one row per pseudobulk containing the columns named in varnames (and N if keep_n = TRUE).

See Also

pseudobulk_deseq2(), compute_hash()

Examples

m <- matrix(sample.int(8, 100*500, replace=TRUE), nrow=100, ncol=500)
rownames(m) <- paste0("G", 1:100)
colnames(m) <- paste0("C", 1:500)
md1 <- sample(c("a", "b"), 500, replace=TRUE)
md2 <- sample(c("c", "d"), 500, replace=TRUE)
df <- data.frame(md1, md2)
data_collapsed <- collapse_counts(m, df, c("md1", "md2"))
head(data_collapsed$counts_mat)
head(data_collapsed$meta_data)

Compute a unique integer hash for each row of a data.frame

Description

Treats the supplied columns as a composite key: rows that agree on every value in vars_use receive the same integer hash, distinct combinations receive distinct hashes. Used internally by collapse_counts() to identify pseudobulk groups, but exposed for callers who need the same indexing in their own code.

Usage

compute_hash(data_df, vars_use)

Arguments

Value

Integer vector of length nrow(data_df).

See Also

collapse_counts()

Group-wise non-zero counts of a matrix along one axis

Description

For each unique value of the grouping vector y, counts the number of non-zero entries among the corresponding rows (or columns) of X. Used internally by wilcoxauc() to compute the percent-expressed columns (pct_in, pct_out), but exposed as a fast group-wise reduction primitive for both dense and dgCMatrix inputs.

Usage

nnzeroGroups(X, y, MARGIN = 2)

## S3 method for class 'dgCMatrix'
nnzeroGroups(X, y, MARGIN = 2)

## S3 method for class 'matrix'
nnzeroGroups(X, y, MARGIN = 2)

Arguments

Value

Integer matrix of shape ⁠n_groups x n_features⁠, where entry ⁠(g, j)⁠ is the number of observations in group g for which feature j is non-zero.

See Also

sumGroups(), wilcoxauc()

Examples

set.seed(42)
exprs <- matrix(rpois(25 * 150, lambda = 2), nrow = 25,
                dimnames = list(paste0("G", 1:25), NULL))
y <- rep(c("A", "B", "C"), each = 50)
nnz_res <- nnzeroGroups(exprs, y, 1)
nnz_res <- nnzeroGroups(t(exprs), y, 2)

Pseudobulk differential expression with DESeq2

Description

Runs DESeq2::DESeq() on a feature-by-pseudobulk count matrix to identify genes whose expression differs between groups of pseudobulks. Three test designs are supported, selected by mode:

• "one_vs_all" (default) - test each level of the contrast variable against the union of all others. Useful for marker discovery. Dispatched to pseudobulk_one_vs_all().

• "pairwise" - test every ordered pair of levels separately. Useful for high-confidence markers; combine with summarize_dge_pairs() to keep the most conservative pair. Dispatched to pseudobulk_pairwise().

• "within" - split on the first term of dge_formula and test the second term within each split level. Useful for condition / case-vs-control comparisons restricted to one cluster at a time. Dispatched to pseudobulk_within().

Pseudobulk inputs are typically produced by collapse_counts(). Genes with low counts across pseudobulks are filtered out before fitting (controlled by min_counts_per_sample and present_in_min_samples). See the pseudobulk vignette for an end-to-end walkthrough on a real single-cell dataset.

Usage

pseudobulk_deseq2(
  dge_formula,
  meta_data,
  counts_df,
  verbose = TRUE,
  min_counts_per_sample = 10,
  present_in_min_samples = 5,
  collapse_background = TRUE,
  vals_test = NULL,
  mode = c("one_vs_all", "pairwise", "within")[1]
)

Arguments

Value

A long-form data.frame of DESeq2 results. Columns include the group identifier(s) (group, or group1 / group2 in pairwise mode), feature, and the DESeq2 columns baseMean, log2FoldChange, lfcSE, stat, pvalue, padj.

See Also

collapse_counts(), top_markers_dds(), summarize_dge_pairs()

Examples

if (requireNamespace("DESeq2", quietly = TRUE)) {
    ## 40 genes x 300 cells from 2 clusters across 6 donors
    m <- matrix(sample.int(8, 40 * 300, replace = TRUE), nrow = 40)
    rownames(m) <- paste0("G", 1:40)
    colnames(m) <- paste0("C", 1:300)
    meta <- data.frame(
        cluster = sample(c("a", "b"), 300, replace = TRUE),
        donor = sample(paste0("d", 1:6), 300, replace = TRUE)
    )

    ## collapse cells into per-(cluster, donor) pseudobulks
    dc <- collapse_counts(m, meta, c("cluster", "donor"))

    ## test each cluster against the rest
    res <- pseudobulk_deseq2(
        ~cluster,
        dc$meta_data["cluster"],
        dc$counts_mat,
        verbose = FALSE,
        present_in_min_samples = 1,
        collapse_background = FALSE,
        mode = "one_vs_all"
    )
    head(res)
}

Pseudobulk DESeq2: one-vs-all contrasts

Description

For each level of contrast_var, fits a DESeq2 model where that level is the foreground and every other level is pooled as the background. Useful for marker-style differential expression: the return is one set of ⁠(log2FoldChange, padj)⁠ per gene per group, giving a quick view of what's elevated in each group relative to the rest.

Usage

pseudobulk_one_vs_all(
  dge_formula,
  counts_df,
  meta_data,
  contrast_var,
  vals_test,
  collapse_background,
  verbose
)

Arguments

Details

Most users should call pseudobulk_deseq2() with mode = "one_vs_all" rather than this function directly; the wrapper handles formula parsing, gene-count filtering, and dispatch. See the pseudobulk vignette for a worked example.

Value

data.frame of DESeq2 results with columns group, feature, baseMean, log2FoldChange, lfcSE, stat, pvalue, padj. Sorted by stat descending within each group.

See Also

pseudobulk_deseq2(), pseudobulk_pairwise(), pseudobulk_within(), top_markers_dds()

Pseudobulk DESeq2: pairwise contrasts

Description

For each ordered pair of contrast levels in vals_test, fits a DESeq2 model on just those two groups of pseudobulks and returns the foreground-vs-background coefficient. This is more conservative than one-vs-all because a marker has to differentiate the foreground from every other group, not just the average background. Cost grows as O(N^2) in the number of levels, so subset to the levels of interest first.

Usage

pseudobulk_pairwise(
  dge_formula,
  counts_df,
  meta_data,
  contrast_var,
  vals_test,
  verbose,
  min_counts_per_sample,
  present_in_min_samples
)

Arguments

Details

Most users should call pseudobulk_deseq2() with mode = "pairwise" rather than this function directly. The companion summarize_dge_pairs() collapses the directional pairs to a single row per (gene, group).

Value

data.frame of DESeq2 results with columns group1, group2, feature, baseMean, log2FoldChange, lfcSE, stat, pvalue, padj. Each row reports the test where pseudobulks of group1 are the foreground and pseudobulks of group2 are the background.

See Also

pseudobulk_deseq2(), pseudobulk_one_vs_all(), pseudobulk_within(), summarize_dge_pairs()

Pseudobulk DESeq2: within-group contrasts

Description

Splits the pseudobulks by split_var and, within each split level, fits a DESeq2 model whose contrast variable is the second term of dge_formula. Used to test for an effect (e.g. case vs. control) restricted to one cluster at a time, where pooling across clusters would mix biology and batch.

Usage

pseudobulk_within(
  dge_formula,
  counts_df,
  meta_data,
  split_var,
  vals_test,
  verbose,
  min_counts_per_sample,
  present_in_min_samples
)

Arguments

Details

Two-level contrasts (factor or character) yield the standard ⁠<var>_<level2>_vs_<level1>⁠ Wald coefficient. Three or more levels are integer-encoded and the fit returns an ordinal trend. Most users should call pseudobulk_deseq2() with mode = "within" rather than this function directly. See the pseudobulk vignette for a worked example using case-vs-control DGE within each cell cluster.

Value

data.frame of DESeq2 results with columns group (the value of split_var), feature, baseMean, log2FoldChange, lfcSE, stat, pvalue, padj.

See Also

pseudobulk_deseq2(), pseudobulk_one_vs_all(), pseudobulk_pairwise()

Column-wise tied ranks of a matrix

Description

Ranks the entries of each column independently using the average rank for ties, and returns the per-column tie group sizes needed for the Wilcoxon variance correction. Used internally by wilcoxauc() (on a transposed input, so that rows become observations) but exposed as a fast standalone ranking primitive for sparse and dense numeric matrices.

Usage

rank_matrix(X)

## S3 method for class 'dgCMatrix'
rank_matrix(X)

## S3 method for class 'matrix'
rank_matrix(X)

Arguments

Value

List with two elements:

• X_ranked - matrix with the same shape as X containing per-column tied ranks.

• ties - list of integer vectors, one per column, giving the sizes of all tie groups encountered in that column. Used by the Wilcoxon statistic to correct for ties.

See Also

wilcoxauc()

Examples

set.seed(42)
exprs <- matrix(rpois(25 * 150, lambda = 2), nrow = 25,
                dimnames = list(paste0("G", 1:25), NULL))
rank_res <- rank_matrix(exprs)

Group-wise sum of a matrix along one axis

Description

For each unique value of the grouping vector y, sums the corresponding rows (or columns) of X. Used internally by wilcoxauc() and collapse_counts(), but exposed as a fast group-wise reduction primitive that works on both dense matrices and dgCMatrix sparse inputs.

Usage

sumGroups(X, y, MARGIN = 2)

## S3 method for class 'dgCMatrix'
sumGroups(X, y, MARGIN = 2)

## S3 method for class 'matrix'
sumGroups(X, y, MARGIN = 2)

Arguments

Value

Numeric matrix of shape ⁠n_groups x n_features⁠. Row order matches the integer order of factor(y).

See Also

nnzeroGroups(), wilcoxauc()

Examples

set.seed(42)
exprs <- matrix(rpois(25 * 150, lambda = 2), nrow = 25,
                dimnames = list(paste0("G", 1:25), NULL))
y <- rep(c("A", "B", "C"), each = 50)
sumGroups_res <- sumGroups(exprs, y, 1)
sumGroups_res <- sumGroups(t(exprs), y, 2)

Summarize directional pairwise DGE results

Description

Collapses the long-form output of pseudobulk_deseq2(mode = "pairwise") to a single row per ⁠(group, gene)⁠ by selecting one comparison per gene-group pair. With mode = "min" it keeps the worst comparison (most conservative — the gene must beat every other group to have a high statistic). With mode = "max" it keeps the best (most permissive — the gene need only beat one).

Usage

summarize_dge_pairs(dge_res, mode = c("min", "max")[1])

Arguments

Value

data.frame with one row per (group, gene), sorted by stat descending within each group. The group2 column is dropped and group1 is renamed to group.

See Also

pseudobulk_deseq2(), pseudobulk_pairwise()

Top markers per group from wilcoxauc results

Description

Filters and ranks the long-form output of wilcoxauc() to give the most distinguishing features per group. The filter arguments combine multiplicatively, then the top n features per group are kept by descending auc and pivoted into wide form. Counterpart to top_markers_dds() for DESeq2-based pseudobulk results.

Usage

top_markers(
  res,
  n = 10,
  auc_min = 0,
  pval_max = 1,
  padj_max = 1,
  pct_in_min = 0,
  pct_out_max = 100
)

Arguments

Value

tibble in wide form: a rank column (1..n) and one column per group containing the feature name of the top-ranked marker at that rank. Cells are NA for groups with fewer than n features that pass the filters.

See Also

wilcoxauc(), top_markers_dds()

Examples

set.seed(42)
exprs <- matrix(rpois(25 * 150, lambda = 2), nrow = 25,
                dimnames = list(paste0("G", 1:25), NULL))
y <- rep(c("A", "B", "C"), each = 50)

res <- wilcoxauc(exprs, y)

## top 10 markers per group, restricted to nominally significant,
## up-regulated features (auc > 0.5 means in-group > out-of-group).
top_markers(res, n = 10, auc_min = 0.5, pval_max = 0.05)

Top markers per group from pseudobulk DESeq2 results

Description

Filters and ranks the long-form output of pseudobulk_deseq2() to give the most distinguishing features per group. The filtering arguments combine multiplicatively, then the top n features per group are kept by descending Wald statistic and pivoted into wide form. Counterpart to top_markers() for Wilcoxon-based results.

Usage

top_markers_dds(res, n = 10, pval_max = 1, padj_max = 1, lfc_min = 1)

Arguments

Value

tibble in wide form: a rank column (1..n) and one column per group containing the gene identifier of the top-ranked feature at that rank. Cells are NA for groups that have fewer than n features passing the filters.

See Also

pseudobulk_deseq2(), top_markers()

Toy SingleCellExperiment object for examples and tests

Description

Builds a small deterministic SingleCellExperiment (20 genes x 300 cells, two cell types in the cell_type colData column) with the installed SingleCellExperiment version, so no serialized object needs to ship with the package. The same simulated counts as toy_seurat() are stored in the counts assay, with log-normalized values in logcounts. The generator restores the caller's random-number state.

Usage

toy_sce()

Details

Requires the suggested package SingleCellExperiment (Bioconductor).

Value

A SingleCellExperiment with counts and logcounts assays and a cell_type colData column with values "jurkat" and "t293".

See Also

toy_seurat(), wilcoxauc()

Examples

if (requireNamespace("SingleCellExperiment", quietly = TRUE)) {
    object_sce <- toy_sce()
    head(wilcoxauc(object_sce, "cell_type"))
}

Toy Seurat object for examples and tests

Description

Builds a small deterministic Seurat object (20 genes x 300 cells, two cell types in the cell_type metadata column) with the installed Seurat version, so no serialized object needs to ship with the package or be updated when Seurat changes its internal class structure. Counts are Poisson-simulated with a few upregulated marker genes per cell type, and log-normalized into the data layer. The generator restores the caller's random-number state, so calling it does not perturb reproducibility.

Usage

toy_seurat()

Details

Requires the suggested package Seurat.

Value

A Seurat object with counts and data layers and a cell_type metadata column with values "jurkat" and "t293".

See Also

toy_sce(), wilcoxauc()

Examples

if (requireNamespace("Seurat", quietly = TRUE)) {
    object_seurat <- toy_seurat()
    head(wilcoxauc(object_seurat, "cell_type"))
}

Fast Wilcoxon rank-sum test and auROC across groups

Description

For every (feature, group) pair, computes the Wilcoxon rank-sum statistic comparing observations in that group against all other observations, and the area under the ROC curve as a measure of separability. P-values come from the standard Gaussian approximation to the U statistic with a tie correction. Returns one row per (feature, group) with effect-size and percent-expressed columns alongside the test statistics.

Usage

wilcoxauc(X, ...)

## S3 method for class 'seurat'
wilcoxauc(X, ...)

## S3 method for class 'Seurat'
wilcoxauc(
  X,
  group_by = NULL,
  assay = "data",
  groups_use = NULL,
  seurat_assay = "RNA",
  ...
)

## S3 method for class 'SingleCellExperiment'
wilcoxauc(X, group_by = NULL, assay = NULL, groups_use = NULL, ...)

## Default S3 method:
wilcoxauc(
  X,
  y,
  groups_use = NULL,
  verbose = TRUE,
  nthreads = 1,
  transposed = FALSE,
  ...
)

Arguments

Details

Designed to be fast enough to run on whole-genome × hundred-thousand- cell single-cell matrices in seconds. Sparse dgCMatrix inputs are processed without densification. Convenience dispatchers extract the counts matrix and group labels from Seurat and SingleCellExperiment objects. See the getting-started vignette for an end-to-end example on a real dataset.

Value

table with the following columns:

• feature - feature name (e.g. gene name).

• group - group name.

• avgExpr - mean value of feature in group.

• logFC - difference of mean feature values between observations in the group vs out of the group. When the input is log-transformed expression (e.g. Seurat's "data" layer or logcounts), this difference of means is a log fold change. On raw (untransformed) values it is a plain difference of means, not a fold change.

• statistic - Wilcoxon rank sum U statistic.

• auc - area under the receiver operator curve.

• pval - nominal p value.

• padj - Benjamini-Hochberg adjusted p value.

• pct_in - Percent of observations in the group with non-zero feature value.

• pct_out - Percent of observations out of the group with non-zero feature value.

See Also

top_markers() to summarize markers per group; pseudobulk_deseq2() for a count-based pseudobulk alternative.

Examples

## generate a tiny toy dataset
set.seed(42)
exprs <- matrix(rpois(25 * 150, lambda = 2), nrow = 25,
                dimnames = list(paste0("G", 1:25), NULL))
y <- rep(c("A", "B", "C"), each = 50)

## on a dense matrix
head(wilcoxauc(exprs, y))

## restrict the comparison to a subset of groups
head(wilcoxauc(exprs, y, c('A', 'B')))

## on a sparse matrix
exprs_sparse <- as(exprs, 'dgCMatrix')
head(wilcoxauc(exprs_sparse, y))

## on a Seurat object (>= v3)
if (requireNamespace("Seurat", quietly = TRUE) &&
    packageVersion("Seurat") >= "3.0") {
    object_seurat <- toy_seurat()
    head(wilcoxauc(object_seurat, 'cell_type'))
}

## on a SingleCellExperiment object
if (requireNamespace("SingleCellExperiment", quietly = TRUE)) {
    object_sce <- toy_sce()
    head(wilcoxauc(object_sce, 'cell_type'))
}