scAmbi: Correcting mapping-ambiguity overdispersion in scRNA-seq

Bootstrap-based overdispersion estimation, correction, and reference comparison

Introduction

Single-cell transcript quantification with Salmon/Alevin produces per-cell bootstrap replicates that capture inferential uncertainty from fragments mapping to multiple transcripts of the same gene^1–3. This uncertainty appears as overdispersion in the estimated counts. If left uncorrected, it inflates variability estimates and distorts downstream analyses including differential expression, clustering, and trajectory inference⁴.

scAmbi estimates per-transcript overdispersion from Alevin bootstraps, computes two additional exploratory components (moment-based and entropy-based), and applies the bootstrap-based estimate by building a corrected assay (RNA_corr) or by supplying model offsets. A per-transcript meanInfRV diagnostic enables distribution-level comparison of inferential uncertainty across transcriptome assemblies quantified with identical reads and parameters.

Relation to edgeR::catchSalmon()

catchSalmon computes per-transcript technical overdispersion from bootstraps pooled across samples, applies EB moderation with limma::squeezeVar(), floors the estimate at 1, and analyzes scaled counts. scAmbi follows the same approach but pools across cells within each sample. The package adds a moment-based component, an annotation prior, and an integrated overdispersion estimate for exploration, and supports offset-based workflows in addition to the corrected assay.

Installation

Install from GitHub:

# install.packages("remotes")
remotes::install_github("sbresnahan/scAmbi")

Or from a local source tarball:

install.packages("scAmbi_0.3.2.tar.gz", repos = NULL, type = "source")

Compatibility with alevin-fry

scAmbi reads Salmon Alevin’s EDS-format output: the transcript-level count matrix (quants_mat.gz) plus per-cell bootstrap replicates (quants_boot_mat.gz, written when Alevin is run with --numCellBootstraps and --dumpFeatures). alevin-fry, the designated successor to Alevin, is not drop-in compatible today: alevin-fry quant -b (--num-bootstraps) resolves UMIs to the gene level and writes Matrix Market output with bootstrap summary statistics (bootstraps_mean.mtx / bootstraps_var.mtx), and full per-replicate bootstrap output is not yet supported in its MTX output. The sufficient statistics scAmbi pools across cells (per-cell bootstrap means and variances) are exactly what alevin-fry’s --summary-stat provides, so support may be added in a future release. In the meantime, fishpond::loadFry() reads alevin-fry output into R/Bioconductor.

Notation and variance model

Transcript-level quantifications exhibit increased technical variance due to read-to-transcript ambiguity. Reads can usually be assigned unambiguously to genes, but they are often compatible with multiple transcripts for the same gene, particularly for genes with many isoforms. This assignment uncertainty disrupts the mean-variance relationship normally observed in RNA-seq data and interferes with standard differential expression methods.

For gene (or transcript) \(g\) in cell \(i\): counts \(y_{gi}\), library size \(N_i\), underlying proportion \(\pi_{gi}\), mean \(\mu_{gi}=N_i\pi_{gi}\). Read-to-transcript ambiguity inflates technical variance beyond Poisson:

\[\mathbb{E}(y_{gi}\mid\pi_{gi})=\mu_{gi}, \qquad \mathrm{Var}(y_{gi}\mid\pi_{gi})=\sigma_g^2\,\mu_{gi}, \; \sigma_g^2>1\]

With biological variation across cells, \(\pi_{gi}\) varies and marginal variance becomes

\[\mathrm{Var}(y_{gi}) = \sigma_g^2\,\mu_{gi} + \phi_g\,\mu_{gi}^2,\]

with \(\phi_g\) the NB-like biological dispersion⁴. The technical overdispersion parameter \(\sigma_g^2\) can be estimated from bootstrap samples and divided out of the transcript counts, yielding scaled counts that follow the standard negative binomial mean-variance relationship and can be analyzed with established differential expression methods.

Integrated overdispersion estimation

scAmbi uses a bootstrap-based overdispersion estimate following Baldoni et al.⁴ for correction. For exploratory purposes, it also computes two additional components per transcript and an integrated estimate by combining the three components on the log scale.

1) Bootstrap-based component (sparse-aware)

Let \(u_{tib}\) be the Alevin bootstrap count for transcript \(t\), cell \(i\), bootstrap \(b = 1..B\). For each cell we traverse the sparse bootstrap columns once, accumulating for each transcript:

• Total sum \(\sum_b u_{tib}\) and sum of squares \(\sum_b u_{tib}^2\) over all \(B\) replicates.
• A seen flag indicating at least one non-zero bootstrap value for that transcript in that cell.

For any transcript marked as seen in a cell, we set \(k = B\) (all replicates, zeros included), compute:

\[\bar{u}_{ti} = \frac{1}{B} \sum_{b=1}^B u_{tib}, \quad s^2_{ti} = \frac{1}{B - 1} \sum_{b=1}^B \left(u_{tib} - \bar{u}_{ti}\right)^2\]

and accumulate across cells the pooled sums \(\sum_i s^2_{ti}\) and \(\sum_i \bar{u}_{ti}\), together with the per-cell inferential relative variance (InfRV, in the sense of Zhu et al.⁵):

\[\mathrm{InfRV}_{ti} = \frac{\max\!\big(s^2_{ti} - \bar{u}_{ti},\, 0\big)}{\bar{u}_{ti} + pc},\]

where \(pc\) is a stabilizing pseudocount (default 1). Subtracting \(\bar{u}_{ti}\) removes the Poisson sampling component, so only excess variance attributable to read-to-transcript ambiguity is measured.

Pooled transcript-level estimate

We pool the sums across cells, mirroring edgeR::catchSalmon(). Pooling weights each cell in proportion to its count; averaging per-cell ratios would weight a 1-count cell the same as a 100-count cell:

\[\mathrm{VMR}_t = \frac{\sum_i s^2_{ti}}{\sum_i \bar{u}_{ti} + pc}.\]

Under the variance model above, \(\mathrm{VMR}_t\) estimates the technical overdispersion \(\sigma^2_t\). A transcript whose quantification is purely Poisson has \(\mathrm{VMR} \approx 1\), and only ambiguity-driven excess variance pushes it above 1.

EB moderation

The pooled VMR is moderated with limma::squeezeVar(), using per-transcript degrees of freedom \(\mathrm{DF}_t \times (B-1)\) (\(\mathrm{DF}_t\) = number of contributing cells) and a prior whose degrees of freedom are estimated from the data, the same strategy as catchSalmon(). If the prior fit fails (e.g., degenerate input), we fall back to fixed-strength shrinkage toward the median with prior df \(d_0 = 10\). The final estimate is floored at 1:

\[\mathrm{OD}^{\mathrm{boot}}_t = \max\!\big(\mathrm{VMR}^{\mathrm{moderated}}_t,\, 1\big) = 1 + \text{pooled Poisson-excess VMR}.\]

Transcripts with \(e_t\) below the expression threshold (min_cells_expr) default to 1.

Notes:

1. We use \(k = B\) for seen transcripts, so implicit zeros contribute to both the mean and the variance.
2. \(\mathrm{OD}^{\mathrm{boot}}_t - 1\) is the pooled analogue of InfRV. The package additionally returns \(\mathrm{meanInfRV}_t = \frac{1}{\mathrm{DF}_t}\sum_i \mathrm{InfRV}_{ti}\), the mean of the per-cell values, as a diagnostic stored in feature_meta. meanInfRV is not used in the correction (see Comparing transcriptome assemblies with meanInfRV).
3. If gene-level OD is desired, aggregate after estimating transcript-level OD (e.g., gene-wise median or counts-weighted mean of transcript ODs).
4. \(\mathrm{OD}^{\mathrm{boot}}_t\) is the default OD estimate for correction.

2) Moment-based component

The moment-based overdispersion estimate captures gene-specific variability by comparing the observed variance to what would be expected under a Poisson model. Sequencing and amplification introduce technical noise that scales with expression level, while biological processes create additional variance through cell-to-cell differences in transcriptional state, cell cycle progression, and stochastic gene expression. Low-abundance genes may also be captured inconsistently across cells due to dropout effects, inflating their apparent variability.

The raw moment-based estimate uses counts across expressing cells only (non-zeros):

\[\widehat{\mathrm{OD}}^{\mathrm{mom}}_g = \frac{\mathrm{Var}(x_g) + \epsilon}{\mathbb{E}(x_g) + \epsilon}\]

where \(\epsilon\) provides numerical stabilization. When rel_eps > 0, both epsilon terms equal rel_eps * mean. When rel_eps = 0 (default), both epsilon terms equal abs_eps (default \(10^{-8}\)). Using only expressing cells avoids bias from technical zeros while capturing the biological variability among cells that detectably express the gene.

Confidence-weighted shrinkage toward the null value of 1 accounts for estimation uncertainty with small sample sizes:

\[\mathrm{OD}^{\mathrm{mom}}_g = \max\!\left(1,\;1+\frac{n_g}{n_g+k}\big(\widehat{\mathrm{OD}}^{\mathrm{mom}}_g-1\big)\right)\]

where \(n_g\) is the number of expressing cells and \(k\) is the shrinkage parameter (default 30). Shrinkage is stronger for genes with fewer expressing cells, reflecting lower confidence in variance estimates from small samples.

This approach may reduce some genuine biological signal along with technical noise, particularly for genes undergoing continuous state transitions or exhibiting meaningful transcriptional heterogeneity. The complexity prior (Section 3) provides gene-specific adjustments that account for expected biological variability in structurally complex genes.

3) Complexity prior (entropy-based)

The complexity prior is motivated by the observation that genes with multiple actively expressed isoforms exhibit increased expression variability due to cell-to-cell differences in splicing regulation and quantification uncertainty when similar isoforms compete for read assignment. The prior combines isoform count with expression evenness to capture this variability.

For each gene \(g\), let \(k_g\) be the number of expressed isoforms (transcripts with mean expression above threshold) and \(E_g\) be the normalized Shannon entropy of isoform expression:

\[E_g = \frac{H_g}{\log(k_g)} = \frac{-\sum_{i} p_i \log(p_i)}{\log(k_g)}\]

where \(p_i\) is the proportion of total gene expression from isoform \(i\). The complexity prior for each transcript is then:

\[\mathrm{OD}^{\mathrm{prior}}_t = 1 + \alpha \cdot \log\big(1 + \max(k_g - 1, 0) \times \max(E_g, 0)\big)\]

where \(\alpha\) is a scaling parameter (default 0.6). This gives higher prior values to transcripts from genes with more expressed isoforms that have relatively even expression levels. Transcripts from single-isoform genes or unmapped transcripts receive a prior of 1.

4) Geometric fusion and adaptive shrinkage

With weights \(w_b=0.8\), \(w_m=0.1\), \(w_p=0.1\):

\[\log \mathrm{OD}^{\mathrm{int}}_t = w_b\log\mathrm{OD}^{\mathrm{boot}}_t + w_m\log\mathrm{OD}^{\mathrm{mom}}_t + w_p\log\mathrm{OD}^{\mathrm{prior}}_t.\]

For additional stability we shrink toward the median of well-observed transcripts (e.g., \(e_t \ge 50\)) with weight \(\alpha_t = \min(e_t/100,1)\):

\[\mathrm{OD}^{\mathrm{int,final}}_t = \alpha_t\,\mathrm{OD}^{\mathrm{int}}_t + (1-\alpha_t)\,\mathrm{median}\big\{\mathrm{OD}^{\mathrm{int}}_{e \ge 50}\big\},\]

then clamp to \([1,100]\).

Applying the correction

\(\mathrm{OD}^{\mathrm{boot}}_t\) is the default OD estimate for correction. The package exposes two interchangeable ways to use the OD estimates.

A) Corrected assay (RNA_corr)

correct_seurat() forms a left-diagonal scaling with entries \(1/\mathrm{OD}_t\):

\[Z = \mathrm{diag}\left(\frac{1}{\mathrm{OD}}\right)\,Y,\]

and stores it as a second assay RNA_corr alongside the raw RNA counts.

The od_source argument selects which estimate drives the correction:

• od_source = "bootstrap" (default): \(\mathrm{OD}^{\mathrm{boot}}_t\), mirroring Baldoni et al. and catchSalmon().
• od_source = "integrated": \(\mathrm{OD}^{\mathrm{int}}_t\), the weighted geometric fusion of all three components (exploratory).
• od_source = "moments" or od_source = "prior": the individual exploratory components.

The chosen vector is stored as scaling_factor in feature_meta and can be retrieved for offset-based workflows via get_od(). Feature-level metadata is stored in both assays with columns:

• OverDisp_integrated, OverDisp_bootstrap, OverDisp_moments, OverDisp_prior
• meanInfRV (bootstrap diagnostic; NA for transcripts below the expression threshold)
• expressing_cells, scaling_factor (the vector selected by od_source), inv_scaling (1 / scaling_factor)

B) Model offsets

When fitting NB-GLMs (e.g., edgeR), you can supply transcript- and sample-specific offsets while leaving the count matrix unmodified.

• Ratio offsets (preferred when RNA_corr exists): use \(\log(\text{corrected}/\text{raw})\) per transcript and sample, e.g., from the ratio of pseudobulked corrected to raw counts.
• OD vector offsets: treat an OD-like vector as a per-transcript multiplicative factor and add \(\log\mathrm{OD}\) to the offset for every sample.
• Arbitrary feature vectors: any column in feature_meta can be used. Interpret as vector_as = "od" for OverDisp_* or scaling_factor, or vector_as = "ratio" for inv_scaling (since RNA_corr = RNA * inv_scaling).

All of these offsets can be constructed from the feature_meta columns retrieved with get_feature() or get_od() and added to the offset slot of your GLM of choice (e.g., edgeR::DGEList(...)$offset).

Comparing transcriptome assemblies with meanInfRV

meanInfRV isolates the mapping-ambiguity-driven component of inferential variance. Its distribution is a useful diagnostic for comparing transcriptome assemblies on identical data, in the spirit of the InfRV diagnostic of Zhu et al.⁵. The comparison is valid only when the same reads, quantifier settings (including --numCellBootstraps), and cell set are used for each assembly. Under those conditions, differences in inferential uncertainty are attributable to the reference.

Recommended practice:

• Compare distributions. Assemblies differ in transcript definition, so per-transcript matching across references is ill-defined. Compare histograms or empirical CDFs of meanInfRV (e.g., medians, upper quantiles, or the fraction of transcripts above a threshold such as 0.2).
• Stratify by expression. Inferential variance is strongly mean-dependent; compare within expression bins (e.g., deciles of mean corrected counts) to avoid confounding by differing expression distributions across assemblies.
• Mind the bootstrap budget. With the typical --numCellBootstraps 20, per-cell variance estimates are noisy. Averaging over cells stabilizes meanInfRV, but small differences (below about 10%) between assemblies should not be over-interpreted.
• 3’-end protocols have weaker signal. In 3’-tag data (e.g., Chromium v3), most ambiguity involves near-identical isoforms sharing 3’ exons, so absolute meanInfRV values are lower and differences between assemblies are compressed relative to 5’-tag protocols, whose reads reach further into the transcript body. The worked example below quantifies both chemistries side by side.

A downward-shifted meanInfRV distribution indicates a reference whose quantifications are less affected by mapping ambiguity for these data.

Example: PBMC 3’ and 5’ libraries against two transcriptome references

Dataset

The worked example uses two public 10x Genomics PBMC libraries from healthy donors:

• 3’ library pbmc_10k_v3 (Chromium 3’ v3; 11,950 cells): 3’-tagged data, where most cDNA reads fall in terminal exons and mapping ambiguity is dominated by isoforms sharing 3’ exons.
• 5’ library 5p_Citrate_CPT (Chromium 5’ v2; 5,355 cells): 5’-tagged data, whose cDNA reads reach further into the transcript body and across splice junctions.

Each library was quantified with Salmon Alevin (--numCellBootstraps 20, --dumpFeatures; library type -l ISR for 3’ and -l ISF for 5’) against two transcriptome references:

• TRAILS, a long-read immune isoform atlas (159,369 isoforms across 17,496 loci, assembled from 29 sorted immune cell subsets)⁶.
• GENCODE v49 (533,740 transcripts)⁷.

Mitochondrial transcripts were excluded from the references before indexing. The TRAILS atlas contains no mitochondrial loci, and the 37 GENCODE MT-gene transcripts (13 mRNAs, 2 rRNAs, 22 tRNAs) were removed from the transcript FASTA, so no mitochondrial rows appear in the quantifications.

Because the same reads, quantifier settings, and bootstrap budget were used for both references, differences in inferential uncertainty are attributable to the reference (see Comparing transcriptome assemblies with meanInfRV). Per-transcript statistics for all four library x reference combinations ship with the package as the pbmc_infrv data object.

Library and quantification statistics

Library	Chemistry	Read pairs	Cells	Reference	Mapping rate	Assigned fragments	UMIs after dedup.
3’ PBMC (pbmc_10k_v3)	10x 3’ v3	638.9M	11,950	TRAILS	41.2%	263.2M	81.7M
				GENCODE v49	48.2%	307.6M	95.4M
5’ PBMC (5p_Citrate_CPT)	10x 5’ v2	122.0M	5,355	TRAILS	54.6%	66.5M	22.8M
				GENCODE v49	62.0%	75.6M	25.4M

Read pairs are implied by assigned fragments / mapping rate. Assigned fragments and mapping rate are from the Salmon logs; UMIs after deduplication are from the Alevin logs. GENCODE v49 gains roughly 7 percentage points of mapping rate over TRAILS on identical reads in both libraries, consistent with its broader transcript catalog. The comparison below asks how the two references differ in the inferential uncertainty of what is quantified.

Reference construction and quantification

Reference construction (TRAILS GTF to transcript FASTA against the hg38/GRCh38 primary assembly; the GENCODE transcript FASTA is used directly after removing MT-gene transcripts):

## TRAILS (the atlas contains no mitochondrial loci)
gffread -w trails_transcripts.fa -g hg38.fa TRAILS.gtf
salmon index -t trails_transcripts.fa -i salmon_idx_trails -p 16

## GENCODE v49: exclude mitochondrial-gene transcripts (GTF seqname chrM)
zcat gencode.v49.annotation.gtf.gz | awk '$1 == "chrM" && $3 == "transcript"' | \
  grep -o 'transcript_id "[^"]*"' | cut -d'"' -f2 | sort -u > mt_tx.txt
seqkit grep -v -f mt_tx.txt gencode.v49.transcripts.fa.gz | \
  gzip > gencode.v49.noMT.transcripts.fa.gz
salmon index -t gencode.v49.noMT.transcripts.fa.gz -i salmon_idx_gencode -p 16 --gencode

Alevin requires a transcript-to-gene map. Passing a transcript-to-self identity table (tx2self.tsv, two identical columns of transcript IDs) keeps quantification at the transcript level, which scAmbi needs.

Quantification (shown for both libraries against TRAILS; the GENCODE runs differ only in the index):

# 3' v3 library (16 bp barcode + 12 bp UMI); ISR orientation
salmon alevin -l ISR \
  -1 pbmc_10k_v3_R1.fastq.gz -2 pbmc_10k_v3_R2.fastq.gz \
  --chromiumV3 -i salmon_idx_trails -p 8 \
  --tgMap tx2self_trails.tsv \
  --numCellBootstraps 20 --dumpFeatures --expectCells 10000 \
  -o quants/3p_trails

# 5' v2 library (16 bp barcode + 10 bp UMI); ISF orientation
salmon alevin -l ISF \
  -1 5p_Citrate_CPT_R1.fastq.gz -2 5p_Citrate_CPT_R2.fastq.gz \
  --chromium -i salmon_idx_trails -p 8 \
  --tgMap tx2self_trails.tsv \
  --numCellBootstraps 20 --dumpFeatures --expectCells 5000 \
  -o quants/5p_trails

Building the pbmc_infrv object

The shipped object is assembled from the four Alevin output directories in three stages (full script: data-raw/pbmc_infrv.R in the source package). None of this runs at vignette build time. The finished object is loaded from package data below.

Stage 1, per-cell collapse. For each library x reference, the EM count matrix (quants_mat.gz) gives mean_count (mean count per cell), expr_cells (number of detecting cells), and pct_expr_cells_threshold (percent of those cells with count >= 10). The bootstrap mean/variance matrices (quants_mean_mat.gz, quants_var_mat.gz) give the per-cell InfRV sums and the pooled VMR accumulators:

library(Matrix); library(eds)

read_eds <- function(f, G, C)
  drop0(as(eds::readEDS(f, numOfGenes = G, numOfOriginalCells = C), "dgCMatrix"))

counts <- read_eds("quants/3p_trails/alevin/quants_mat.gz", G, C)
mean_count <- rowSums(counts) / C
expr_cells <- tabulate(counts@i + 1L, nbins = G)          # detection count
n_ge_10    <- tabulate(counts@i[counts@x >= 10] + 1L, nbins = G)

mean_m <- read_eds("quants/3p_trails/alevin/quants_mean_mat.gz", G, C)
var_m  <- read_eds("quants/3p_trails/alevin/quants_var_mat.gz", G, C)
V <- var_m + 0 * mean_m; M <- mean_m + 0 * var_m          # union support
infrv <- pmax(V@x - M@x, 0) / (M@x + 1)                   # per-cell InfRV

Stage 2, component 1 (OD_boot). The pooled VMR is EB-moderated with limma::squeezeVar() on transcripts detected in >= 10 cells, per-transcript df = expr_cells * (n_boot - 1), floored at 1. This is the same fit od_boot() performs (see Integrated overdispersion estimation, Section 1).

Stage 3, exploratory components 2/3 and the integrated fusion. OD_moments and OD_prior are computed from the EM counts with od_moments() and od_prior() (package defaults). OD_integrated fuses the three components exactly as od_integrated() does: weighted geometric mean (weights 0.8/0.1/0.1), shrinkage toward the median of the expr_cells >= 50 set with weight pmin(expr_cells/100, 1), clamped to [1, 100]:

od_m <- od_moments(counts, min_cells = 10)
od_p <- od_prior(counts, transcript_info = tx2gene, features = tx_ids)

w <- c(bootstrap = 0.8, moments = 0.1, prior = 0.1); w <- w / sum(w)
odi <- exp(w["bootstrap"] * log(pmax(od_boot_vec, 1)) +
           w["moments"]   * log(od_m) +
           w["prior"]     * log(od_p))
hc <- which(expr_cells >= 50)
sw <- pmin(expr_cells / 100, 1)
odi <- pmin(pmax(sw * odi + (1 - sw) * median(odi[hc]), 1), 100)

The pbmc_infrv data object

library(scAmbi)
data("pbmc_infrv")
head(pbmc_infrv)

#>                            transcript_id library reference mean_count
#> 1 028da976-0c67-417e-bb9e-47ad0934e973-1      3p    TRAILS     122.20
#> 2                     ENST00000260379.11      3p    TRAILS      53.69
#> 3   dc4ab983-60ca-4a52-bb2d-5a3512d97312      3p    TRAILS      49.67
#> 4                      ENST00000361575.4      3p    TRAILS      38.68
#> 5                      ENST00000521726.1      3p    TRAILS      37.70
#> 6   5938151a-65fa-4f9b-a8f3-ac185c22f601      3p    TRAILS      37.44
#>   expr_cells pct_expr_cells_threshold meanInfRV OD_boot VMR_raw OD_moments
#> 1      11867                    95.85    5.3310   5.941  5.9410      44.38
#> 2      11652                    92.19    0.1219   1.000  0.9988      24.30
#> 3      11601                    89.78    1.1540   2.040  2.0400      24.14
#> 4      11465                    87.86    0.0961   1.000  0.9533      18.02
#> 5      11539                    92.06    0.1264   1.004  1.0040      13.47
#> 6      11607                    67.15    0.9965   2.154  2.1540      34.33
#>   OD_prior OD_integrated         gene_id gene_symbol
#> 1    2.323         6.613  chr11:65499000        <NA>
#> 2    1.224         1.404 ENSG00000137818       RPLP1
#> 3    1.609         2.550 ENSG00000133112        TPT1
#> 4    1.310         1.372 ENSG00000198918       RPL39
#> 5    1.034         1.306 ENSG00000156482       RPL30
#> 6    2.653         2.901 ENSG00000167526       RPL13
#>                     gene_biotype
#> 1                           <NA>
#> 2                 protein_coding
#> 3                 protein_coding
#> 4                 protein_coding
#> 5 protein_coding_CDS_not_defined
#> 6                 protein_coding

Each row is one transcript in one library x reference combination. The object is unfiltered (every quantified transcript; references were built MT-free):

• transcript_id, gene_id, gene_symbol, gene_biotype: feature identifiers. TRAILS transcript IDs are long-read assembly UUIDs; gene IDs are version-stripped Ensembl gene IDs shared across both references.
• library, reference: 3p/5p and TRAILS/GENCODE.
• mean_count: mean Alevin EM count per cell.
• expr_cells: number of cells in which the transcript was detected.
• pct_expr_cells_threshold: percent of expr_cells with EM count >= 10 (NA when expr_cells = 0).
• meanInfRV: mean per-cell inferential relative variance (diagnostic; not used in the correction).
• OD_boot: EB-moderated pooled bootstrap overdispersion (component 1, floored at 1); the default correction vector. NA outside the fit set (expr_cells < 10).
• VMR_raw: unmoderated pooled variance-to-mean ratio.
• OD_moments, OD_prior, OD_integrated: exploratory components 2 and 3 and their fusion (see Comparing the overdispersion components).

Reference-comparison figure

The figure below compares the two references on identical reads, per library. Transcripts first pass a labelKeep-style abundance filter adapted from fishpond::labelKeep() (kept: expr_cells > 10 and pct_expr_cells_threshold > 0%, i.e., detected in more than 10 cells with at least one cell clearing 10 UMIs). All abundance plots use log10(meanInfRV + 0.01). The pseudocount is required because the median meanInfRV is exactly 0 in every combination.

suppressPackageStartupMessages({
  library(dplyr); library(tidyr); library(ggplot2)
  library(scales); library(patchwork); library(ggpointdensity)
})

PC <- 0.01        # pseudocount for the log10 transform
LK_MINN <- 10     # keep: expr_cells strictly greater than this
LK_PCT  <- 0      # keep: pct_expr_cells_threshold strictly greater than this

cols <- c(TRAILS = "royalblue", GENCODE = "orange")
lib_labels <- c(`3p` = "3' PBMC (10x v3)", `5p` = "5' PBMC (10x v2)")

label_keep <- function(expr_cells, pct_expr_cells_threshold,
                       minN = LK_MINN, minPct = LK_PCT) {
  expr_cells > minN & pct_expr_cells_threshold > minPct
}

## snapshot the pre-filter table: the filter-step panels need every
## transcript that entered the filter, including those later removed
prefilter_tx <- pbmc_infrv

pbmc_infrv_filtered <- pbmc_infrv %>%
  filter(label_keep(expr_cells, pct_expr_cells_threshold))

## sequential filter funnel: step 1 removes transcripts failing
## expr_cells > LK_MINN; step 2 removes the survivors failing
## pct_expr_cells_threshold > LK_PCT (never-expressed transcripts,
## pct = NA, fail both)
step_data <- prefilter_tx %>%
  mutate(fail_expr = !(expr_cells > LK_MINN),
         fail_pct  = is.na(pct_expr_cells_threshold) |
           pct_expr_cells_threshold <= LK_PCT) %>%
  group_by(library, reference) %>%
  summarise(total = n(),
            step1 = sum(fail_expr),
            step2 = sum(!fail_expr & fail_pct),
            kept  = sum(!fail_expr & !fail_pct),
            .groups = "drop")

step_data %>%
  mutate(across(c(step1, step2, kept),
                ~ sprintf("%s (%.1f%%)", format(.x, big.mark = ","),
                          100 * .x / total))) %>%
  as.data.frame()

#>   library reference  total           step1           step2         kept
#> 1      3p    TRAILS 159368    3,961 (2.5%) 151,024 (94.8%) 4,383 (2.8%)
#> 2      3p   GENCODE 517001 236,374 (45.7%) 276,624 (53.5%) 4,003 (0.8%)
#> 3      5p    TRAILS 159368    7,499 (4.7%) 149,331 (93.7%) 2,538 (1.6%)
#> 4      5p   GENCODE 517001 272,417 (52.7%) 242,286 (46.9%) 2,298 (0.4%)

plot_all <- pbmc_infrv_filtered %>%
  filter(!is.na(meanInfRV)) %>%
  mutate(log10_InfRV = log10(meanInfRV + PC))

step_long <- step_data %>%
  select(library, reference, total, step1, step2) %>%
  pivot_longer(c(step1, step2), names_to = "step", values_to = "n") %>%
  mutate(pct = 100 * n / total,
         step = factor(step, levels = c("step1", "step2"),
                       labels = c(sprintf("1. expr_cells > %g", LK_MINN),
                                  sprintf("2. pct_expr_cells_threshold\n> %g%%", LK_PCT))))

## gene-level aggregation: expression-weighted mean of transcript meanInfRV,
## restricted to genes with at least one filtered transcript in BOTH references
gene_level <- plot_all %>%
  filter(!is.na(gene_id)) %>%
  group_by(library, reference, gene_id) %>%
  summarise(wInfRV = weighted.mean(meanInfRV, w = mean_count), .groups = "drop") %>%
  pivot_wider(names_from = reference, values_from = wInfRV) %>%
  filter(!is.na(TRAILS), !is.na(GENCODE)) %>%
  mutate(log10_TRAILS = log10(TRAILS + PC),
         log10_GENCODE = log10(GENCODE + PC))

frac_below <- gene_level %>%
  group_by(library) %>%
  summarise(n_genes = n(),
            frac_TRAILS_lt_GENCODE = mean(TRAILS < GENCODE),
            .groups = "drop")

base_theme <- theme_bw(base_size = 10) +
  theme(strip.text = element_text(face = "bold", size = 10),
        legend.position = "top",
        legend.title = element_blank(),
        plot.tag = element_text(face = "bold", size = 12))

make_filterbar <- function(lib) {
  d <- filter(step_long, library == lib)
  kept <- filter(step_data, library == lib) %>%
    mutate(txt = sprintf("%s %s (%.1f%%)", format(kept, big.mark = ","),
                         reference, 100 * kept / total))
  ggplot(d, aes(x = step, y = pct, fill = reference)) +
    geom_col(position = position_dodge(width = 0.75), width = 0.65) +
    geom_text(aes(label = sprintf("%.1f%%", pct), group = reference),
              position = position_dodge(width = 0.75),
              vjust = -0.4, size = 2.8) +
    scale_fill_manual(values = cols) +
    scale_y_continuous(limits = c(0, 100),
                       expand = expansion(mult = c(0, 0.08))) +
    labs(x = "filter step (applied sequentially)",
         y = "% of transcripts filtered",
         title = lib_labels[[lib]],
         subtitle = sprintf("kept after both steps: %s",
                            paste(kept$txt, collapse = " / "))) +
    base_theme
}

make_ecdf <- function(lib) {
  ggplot(filter(plot_all, library == lib),
         aes(x = log10_InfRV, color = reference)) +
    stat_ecdf(linewidth = 0.8) +
    scale_color_manual(values = cols) +
    labs(x = bquote(log[10]~"(meanInfRV + 0.01)"),
         y = "ECDF", title = lib_labels[[lib]]) +
    base_theme
}

make_scatter <- function(lib) {
  d  <- filter(gene_level, library == lib)
  fb <- frac_below[frac_below$library == lib, ]
  ggplot(d, aes(x = log10_GENCODE, y = log10_TRAILS)) +
    geom_pointdensity(size = 0.8, alpha = 0.9) +
    scale_color_viridis_c(
      option = "viridis",
      guide = guide_colorbar(barwidth = 8, barheight = 0.6,
                             title = "local density",
                             title.position = "top")) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
    labs(x = bquote("GENCODE v49 gene-level"~log[10]~"(meanInfRV + 0.01)"),
         y = bquote("TRAILS gene-level"~log[10]~"(meanInfRV + 0.01)"),
         title = lib_labels[[lib]],
         subtitle = sprintf(paste0(
           "%s shared genes; %.1f%% below diagonal (TRAILS < GENCODE)\n",
           "gene-level InfRV = mean_count-weighted mean of transcript meanInfRV"),
           format(fb$n_genes, big.mark = ","),
           100 * fb$frac_TRAILS_lt_GENCODE)) +
    base_theme +
    theme(plot.subtitle = element_text(size = 8))
}

fig_comparison <- (make_filterbar("3p") | make_filterbar("5p")) /
  (make_ecdf("3p") | make_ecdf("5p")) /
  (make_scatter("3p") | make_scatter("5p")) +
  plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(face = "bold", size = 12))

fig_comparison

Interpretation

Every pattern in the figure traces back to how each catalog was built and to where the reads land. TRAILS is a long-read atlas, so it carries many low-abundance novel isoforms that are detected in more than 10 cells yet never clear 10 UMIs in any single cell; its removals are therefore dominated by the second filter step. GENCODE is a comprehensive annotation, so a large fraction of its transcripts are never or barely expressed in PBMCs at all and are removed at the first step. The 5’ library carries substantially more inferential uncertainty than the 3’ library under either reference because 5’ reads reach into transcript bodies and across splice junctions, where isoforms diverge, whereas 3’ reads concentrate in terminal exons shared across isoforms. The reference ordering is chemistry-dependent for the same reason: in 3’ data, where most fragments map to shared terminal exons, the two references are nearly indistinguishable, while in 5’ data the experimentally resolved junction structures in TRAILS give fragments a unique target more often than GENCODE’s larger, partly redundant catalog, shifting the TRAILS distribution left and placing the expression-weighted gene-level InfRV below the diagonal for a clear majority of shared genes. Because the assemblies define different transcripts, all comparisons here are distributional or gene-level, and differences below about 10% should not be over-interpreted given the 20-bootstrap budget.

Comparing the overdispersion components (exploratory)

The shipped object carries all three overdispersion components plus their integrated fusion. Components 2 (OD_moments), 3 (OD_prior), and OD_integrated are exploratory. Component 1 (OD_boot) remains the default correction vector. The scatter plots below use the labelKeep-filtered transcripts (the analysis set).

od_filtered <- pbmc_infrv_filtered %>%
  mutate(combo = paste0(lib_labels[as.character(library)], " / ", reference))

pair_scatter <- function(xvar, yvar) {
  rho <- od_filtered %>%
    group_by(combo) %>%
    summarise(r = cor(.data[[xvar]], .data[[yvar]], method = "spearman"),
              .groups = "drop")
  d <- od_filtered %>%
    left_join(rho, by = "combo") %>%
    mutate(panel = sprintf("%s\nSpearman rho = %.2f", combo, r))
  ggplot(d, aes(x = .data[[xvar]], y = .data[[yvar]])) +
    geom_pointdensity(size = 0.5, alpha = 0.8) +
    scale_color_viridis_c(option = "viridis", guide = "none") +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
    scale_x_log10() + scale_y_log10() +
    facet_wrap(~ panel, nrow = 1) +
    labs(x = paste0(xvar, " (log10)"), y = paste0(yvar, " (log10)")) +
    theme_bw(base_size = 9) +
    theme(strip.text = element_text(size = 7.5))
}

pair_scatter("OD_boot", "OD_moments") /
  pair_scatter("OD_boot", "OD_integrated") /
  pair_scatter("OD_moments", "OD_prior") +
  plot_annotation(tag_levels = "A")

The three components measure different quantities, and the differences between them are mechanistic. OD_moments sits at 1 for the large majority of transcripts because the moment estimator on nonzero counts is dominated by ordinary sampling noise at these expression levels: after shrinkage, most transcripts have variance at or below the Poisson mean among expressing cells, and only highly variable transcripts form a long right tail. OD_prior is >= 1 everywhere by construction and identical for all transcripts of a gene because it encodes the potential for mapping ambiguity implied by the annotation; most expressed transcripts belong to multi-isoform genes, so the prior assigns them baseline complexity independently of the data. OD_boot is the only component that directly measures realized mapping-ambiguity variance, which is why it remains the default correction vector. OD_integrated is anchored by the bootstrap component, pulled upward where the prior indicates complex loci and shrunk toward the high-confidence median for sparsely detected transcripts, and is best viewed as a sensitivity analysis around the bootstrap estimate.

From uncertainty estimates to count correction

The object above carries the diagnostic layer. The correction workflow operates on Alevin output directories directly. With the PBMC quants above, the full pipeline is:

## per-library overdispersion from bootstraps
od <- od_boot(
  alevin_dir = "quants/5p_trails/alevin",
  n_boot = 20, min_cells_expr = 10
)

## read counts and build a corrected Seurat assay (RNA_corr)
s <- read_alevin(
  sample_id = "5p_trails",
  base_dir  = "quants",
  n_boot    = 20
)
seu <- correct_seurat(
  sample_id = "5p_trails",
  counts = s$counts, od = s$od, feats = s$feats, cells = s$cells,
  od_source = "bootstrap"   # default; "integrated", "moments", "prior" also possible
)

## optional: sanity-check the corrected assay
diagnose_correction(list(`5p_trails` = seu))

Further reading

• The function reference (help(package = "scAmbi")) documents the full API, including od_moments(), od_prior(), od_integrated(), correct_seurat(), diagnose_correction(), get_od(), and get_feature().
• The source package includes data-raw/pbmc_infrv.R, the full reproducible script for building the shipped data object from the public FASTQ files and reference annotations.

Session information

sessionInfo()

#> R version 4.4.3 (2025-02-28)
#> Platform: x86_64-conda-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS/LAPACK: /opt/conda/lib/libopenblasp-r0.3.34.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggpointdensity_0.2.1 patchwork_1.3.2      scales_1.4.0        
#> [4] ggplot2_4.0.3        tidyr_1.3.2          dplyr_1.2.1         
#> [7] scAmbi_0.3.2        
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3     jsonlite_2.0.0         magrittr_2.0.5        
#>   [4] spatstat.utils_3.2-4   farver_2.1.2           rmarkdown_2.31        
#>   [7] ragg_1.5.2             vctrs_0.7.3            ROCR_1.0-12           
#>  [10] spatstat.explore_3.8-2 base64enc_0.1-6        htmltools_0.5.9       
#>  [13] sass_0.4.10            sctransform_0.4.3      parallelly_1.48.0     
#>  [16] KernSmooth_2.23-26     bslib_0.12.0           htmlwidgets_1.6.4     
#>  [19] ica_1.0-3              plyr_1.8.9             plotly_4.12.1         
#>  [22] zoo_1.8-15             cachem_1.1.0           uuid_1.2-2            
#>  [25] igraph_2.3.3           mime_0.13              lifecycle_1.0.5       
#>  [28] pkgconfig_2.0.3        Matrix_1.7-6           R6_2.6.1              
#>  [31] fastmap_1.2.0          fitdistrplus_1.2-6     future_1.75.0         
#>  [34] shiny_1.14.0           digest_0.6.39          Seurat_5.5.1          
#>  [37] tensor_1.5.1           RSpectra_0.16-2        irlba_2.3.7           
#>  [40] textshaping_1.0.1      labeling_0.4.3         progressr_1.0.0       
#>  [43] spatstat.sparse_3.2-0  httr_1.4.8             polyclip_1.10-7       
#>  [46] abind_1.4-8            compiler_4.4.3         withr_3.0.3           
#>  [49] S7_0.2.2               fastDummies_1.7.6      MASS_7.3-66           
#>  [52] tools_4.4.3            lmtest_0.9-40          otel_0.2.0            
#>  [55] httpuv_1.6.17          future.apply_1.20.2    goftest_1.2-3         
#>  [58] glue_1.8.1             nlme_3.1-170           promises_1.5.0        
#>  [61] grid_4.4.3             pbdZMQ_0.3-14          Rtsne_0.17            
#>  [64] cluster_2.1.8.3        reshape2_1.4.5         generics_0.1.4        
#>  [67] gtable_0.3.6           spatstat.data_3.1-9    data.table_1.18.4     
#>  [70] sp_2.2-3               spatstat.geom_3.8-2    RcppAnnoy_0.0.23      
#>  [73] ggrepel_0.9.8          RANN_2.6.2             pillar_1.11.1         
#>  [76] stringr_1.6.0          spam_2.11-4            IRdisplay_1.1         
#>  [79] RcppHNSW_0.7.0         later_1.4.8            splines_4.4.3         
#>  [82] lattice_0.23-1         survival_3.8-9         deldir_2.0-4          
#>  [85] tidyselect_1.2.1       miniUI_0.1.2           pbapply_1.7-4         
#>  [88] knitr_1.51             gridExtra_2.3.1        scattermore_1.2       
#>  [91] xfun_0.60              matrixStats_1.5.0      eds_1.8.0             
#>  [94] stringi_1.8.7          yaml_2.3.12            evaluate_1.0.5        
#>  [97] codetools_0.2-20       tibble_3.3.1           cli_3.6.6             
#> [100] uwot_0.2.5             IRkernel_1.3.2         xtable_1.8-8          
#> [103] reticulate_1.46.0      systemfonts_1.3.2      repr_1.1.7            
#> [106] jquerylib_0.1.4        Rcpp_1.1.2             globals_0.19.1        
#> [109] spatstat.random_3.5-1  png_0.1-9              spatstat.univar_3.2-0 
#> [112] parallel_4.4.3         dotCall64_1.2          listenv_1.0.0         
#> [115] viridisLite_0.4.3      ggridges_0.5.7         SeuratObject_5.4.0    
#> [118] purrr_1.2.2            crayon_1.5.3           rlang_1.3.0           
#> [121] cowplot_1.2.0

References

1.

Srivastava, A., Malik, L., Sarkar, H. & Patro, R. A Bayesian framework for inter-cellular information sharing improves dscRNA-seq quantification. Bioinformatics 36, i292–i299 (2020).

2.

Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14, 417–419 (2017).

3.

Srivastava, A., Malik, L., Smith, T., Sudbery, I. & Patro, R. Alevin efficiently estimates accurate gene abundances from dscRNA-seq data. Genome Biology 20, 65 (2019).

4.

Baldoni, P. L. et al. Dividing out quantification uncertainty allows efficient assessment of differential transcript expression with edgeR. Nucleic Acids Research 52, e13–e13 (2024).

5.

Zhu, A., Srivastava, A., Ibrahim, J. G., Patro, R. & Love, M. I. Nonparametric expression analysis using inferential replicate counts. Nucleic Acids Research 47, e105 (2019).

6.

Inamo, J. et al. Long-read sequencing for 29 immune cell subsets reveals disease-linked isoforms. Nature Communications 15, 4285 (2024).

7.

Frankish, A. et al. GENCODE: Reference annotation for the human and mouse genomes in 2023. Nucleic Acids Research 51, D942–D949 (2023).