Bioconductor Dropletutils

DropletUtils — Comprehensive Skill Guide

Domain: Single-cell / droplet Bioconductor: DropletUtils Paper: Lun ATL et al. (2019). Genome Biology, 20:63.

Utilities for handling droplet-based single-cell RNA-seq data: emptyDrops for ambient RNA removal, barcodeRanks for knee detection, and read10xCounts for CellRanger output loading.

Dependencies & Environment

Package-intrinsic requirements from the Bioconductor landing page — reproduce in any R environment.

• Version: 1.32.0 · Bioconductor: 3.23 · R: ≥ 4.6
• Depends: SingleCellExperiment
• Imports: Matrix, Rcpp, BiocGenerics, S4Vectors, IRanges, GenomicRanges, SummarizedExperiment, BiocParallel, SparseArray, DelayedArray, DelayedMatrixStats, HDF5Array, rhdf5, edgeR, R.utils, dqrng, beachmat, scuttle
• System requirements: C++17, GNU make
• Install: BiocManager::install("DropletUtils")

When to Use

• Load 10x Genomics CellRanger output (barcodes/features/matrix)
• Distinguish real cells from empty droplets (emptyDrops)
• Detect knee/inflection points in barcode-rank plots
• Correct for ambient RNA contamination

Alternatives: CellRanger filtered_feature_bc_matrix, SoupX (ambient RNA)

Do NOT Use When

• Plate-based scRNA-seq (Smart-seq2) — no droplet model needed.
• Pre-filtered CellRanger output — use raw (unfiltered) for emptyDrops.

Data Requirements

Installation

if (!require("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("DropletUtils")
library(DropletUtils)

Workflows

Load 10x Data and Remove Empty Droplets

Read CellRanger output and use emptyDrops to filter real cells

Load unfiltered matrix

library(DropletUtils)
# Point to raw (unfiltered) CellRanger output
sce <- read10xCounts("path/to/raw_feature_bc_matrix/",
                      col.names=TRUE)
dim(sce)  # all barcodes, including empty

Inspect barcode-rank plot

br_out <- barcodeRanks(counts(sce))
plot(br_out$rank, br_out$total + 1, log="xy",
     xlab="Rank", ylab="Total UMI", main="Barcode Rank Plot")
abline(h=metadata(br_out)$knee, col="blue", lty=2)
abline(h=metadata(br_out)$inflection, col="red", lty=2)
legend("bottomleft", c("Knee","Inflection"),
       col=c("blue","red"), lty=2)

Run emptyDrops

set.seed(100)
e_out <- emptyDrops(counts(sce), lower=100, niters=10000,
                    BPPARAM=BiocParallel::SerialParam())

# Filter to real cells at FDR < 0.001
keep <- e_out$FDR <= 0.001
keep[is.na(keep)] <- FALSE
sce_filtered <- sce[, keep]
dim(sce_filtered)
cat("Cells retained:", ncol(sce_filtered), "
")

Downstream QC

library(scater)
# Compute per-cell QC metrics
is_mito <- grepl("^MT-", rownames(sce_filtered))
df <- perCellQCMetrics(sce_filtered, subsets=list(Mito=is_mito))
colData(sce_filtered) <- cbind(colData(sce_filtered), df)

# Remove low-quality cells
qc_out <- quickPerCellQC(df, percent_subsets="subsets_Mito_percent")
sce_filtered <- sce_filtered[, !qc_out$discard]
cat("After QC:", ncol(sce_filtered), "cells
")

Key Functions & Parameters

read10xCounts()()

Load CellRanger output into a SingleCellExperiment

emptyDrops()()

Test each barcode for ambient RNA; returns statistics for classification

barcodeRanks()()

Compute barcode-rank statistics for knee plot visualization

Scientific Assumptions

The method assumes the following about your data and experiment:

• Empty droplets contain only ambient RNA from the solution, not real cells.
• The Monte Carlo test in emptyDrops assumes the ambient profile can be estimated from low-UMI barcodes.
• Cells have significantly higher total UMI counts than empty droplets.

Result Interpretation

• emptyDrops()$FDR: proportion of expected false discoveries; use FDR ≤ 0.001 for strict filtering.
• barcodeRanks()$knee: UMI count at the knee of the rank curve — approximate threshold for cell/empty distinction.
• barcodeRanks()$inflection: lower bound for cell-containing droplets.
• Cells called by emptyDrops at FDR<0.001 can include cells with very low UMI that CellRanger missed.

Best Practices

• Use emptyDrops() on the UNFILTERED CellRanger output (raw_feature_bc_matrix/) for best sensitivity.
• Filter emptyDrops() results: keep <- e.out$FDR <= 0.001 & !is.na(e.out$FDR).
• Check the knee plot (barcodeRanks) before choosing the lower UMI threshold for emptyDrops.
• After emptyDrops, run sce <- sce[, keep] to subset to real cells.
• For standard 10x data, emptyDrops(lower=100, niters=10000) is a good default.

When to Choose This vs Alternatives

Benchmark Evidence

Key studies that have validated or benchmarked this tool:

• Lun ATL et al. (2019). Genome Biol 20:63. PMID:30902100
  emptyDrops recovers more real cells than CellRanger filtering, especially cells with low UMI counts.

Common Errors & Troubleshooting

Error in emptyDrops: all barcodes have total counts <= lower

Cause: lower parameter set too high relative to data

Fix: Reduce lower; inspect barcodeRanks(m) to set appropriate value

Additional Notes from Official Documentation

Extracted from the DropletUtils Bioconductor vignette(s)

Contents

1 Introduction 2 Reading in 10X Genomics data 2.1 From the UMI count matrix 2.2 From the molecule information file 2.1 From the UMI count matrix 2.2 From the molecule information file 3 Downsampling on the reads 4 Computing barcode ranks 5 Detecting empty droplets 6 Demultiplexing hashed libraries 7 Removing swapping effects 7.1 Barcode swapping between samples 7.2 Chimeric reads within cells 7.1 Barcode swapping between samples 7.2 Chimeric reads within cells 8 Session information

1 Introduction

Droplet-based single-cell RNA sequencing (scRNA-seq) technologies allow researchers to obtain transcriptome-wide expression profiles for thousands of cells at once. Briefly, each cell is encapsulated in a droplet in a oil-water emulsion, along with a bead containing reverse transcription primers with a unique barcode sequence. After reverse transcription inside the droplet, each cell’s cDNA is labelled with that barcode (referred to a “cell barcode”). Bursting of the droplets yields a pool of cDNA for library preparation and sequencing. Debarcoding of the sequences can then be performed to obtain the expression profile for each cell. This package implements some general utilities for handling these data after quantification of expression. In particular, we focus on the 10X Genomics

2.1 From the UMI count matrix

The CellRanger pipeline from 10X Genomics will process the raw sequencing data and produce a matrix of UMI counts. Each row of this matrix corresponds to a gene, while each column corresponds to a cell barcode. This is saved in a single directory for each sample, usually named like /outs/filtered_gene_bc_matrices/ 1 1 1 If you use the “filtered” matrix, each column corresponds to a putative cell. If you use the “raw” matrix, all barcodes are loaded, and no distinction is made between cells and empty droplets. . We mock up an example directory below using some simulated data: The matrix.mtx file contains the UMI counts, while the other two files contain the cell barcodes and the gene annotation. We can load this into memory using the read10xCounts function, which re

# To generate the files.
example(write10xCounts, echo=FALSE) 
dir.name <- tmpdir
list.files(dir.name)

## [1] "barcodes.tsv" "genes.tsv"    "matrix.mtx"

2.2 From the molecule information file

CellRanger will also produce a molecule information file ( molecule_info.h5 ) that contains… well, information about the transcript molecules. This includes the UMI sequence 2 2 2 For readers who are unfamiliar with UMIs, they allow reads from different PCR amplicons to be unambiguously assigned to the same original molecule. , the cell barcode sequence, the gene to which it was assigned, and the number of reads covering the molecule. For demonstration purposes, we create an example molecule information file below: We can subsequently load this information into our R session using the read10xMolInfo function: This information can be useful for quality control purposes, especially when the underlying read counts are required, e.g., to investigate sequencing saturation. Note that the func

set.seed(1000)
mol.info.file <- DropletUtils:::simBasicMolInfo(tempfile())
mol.info.file

## [1] "/tmp/Rtmp5yfePp/file1d5dd58cf99eb"

3 Downsampling on the reads

Given multiple batches of very different sequencing depths, it can be beneficial to downsample the deepest batches to match the coverage of the shallowest batches. This avoids differences in technical noise that can drive clustering by batch. The scuttle package provides some utilities to downsample count matrices, but technically speaking, downsampling on the reads is more appropriate as it recapitulates the effect of differences in sequencing depth per cell. This can be achieved by applying the downsampleReads function to the molecule information file containing the read counts: The above code will downsample the reads to 50% of the original coverage across the experiment. However, the function will return a matrix of UMI counts , so the final total count may not actually decrease if th

set.seed(100)
no.sampling <- downsampleReads(mol.info.file, prop=1)
sum(no.sampling)

## [1] 9532

4 Computing barcode ranks

A useful diagnostic for droplet-based data is the barcode rank plot, which shows the (log-)total UMI count for each barcode on the y-axis and the (log-)rank on the x-axis. This is effectively a transposed empirical cumulative density plot with log-transformed axes. It is useful as it allows users to examine the distribution of total counts across barcodes, focusing on those with the largest counts. To demonstrate, let us mock up a count matrix: We compute the statistics using the barcodeRanks function, and then create the plot as shown below. The knee and inflection points on the curve mark the transition between two components of the total count distribution. This is assumed to represent the difference between empty droplets with little RNA and cell-containing droplets with much more RNA

set.seed(0)
my.counts <- DropletUtils:::simCounts()

br.out <- barcodeRanks(my.counts)

# Making a plot.
plot(br.out$rank, br.out$total, log="xy", xlab="Rank", ylab="Total")
o <- order(br.out$rank)
lines(br.out$rank[o], br.out$fitted[o], col="red")

abline(h=metadata(br.out)$knee, col="dodgerblue", lty=2)
abline(h=metadata(br.out)$inflection, col="forestgreen", lty=2)
legend("bottomleft", lty=2, col=c("dodgerblue", "forestgreen"), 
    legend=c("knee", "inflection"))

5 Detecting empty droplets

Empty droplets often contain RNA from the ambient solution, resulting in non-zero counts after debarcoding. The emptyDrops function is designed to distinguish between empty droplets and cells. It does so by testing each barcode’s expression profile for significant deviation from the ambient profile. Given a matrix my.counts containing UMI counts for all barcodes, we call: Droplets with significant deviations from the ambient profile are detected at a specified FDR threshold, e.g., with FDR below 1%. These can be considered to be cell-containing droplets, with a frequency of false positives (i.e., empty droplets) at the specified FDR. Furthermore, droplets with very large counts are automatically retained by setting their p -values to zero. This avoids discarding droplets containing cell

set.seed(100)
e.out <- emptyDrops(my.counts)
e.out

## DataFrame with 11100 rows and 5 columns
##           Total   LogProb    PValue   Limited         FDR
##       <integer> <numeric> <numeric> <logical>   <numeric>
## 1             2        NA        NA        NA          NA
## 2             9        NA        NA        NA          NA
## 3            20        NA        NA        NA          NA
## 4            20        NA        NA        NA          NA
## 5             1        NA        NA        NA          NA
## ...         ...       ...       ...       ...         ...
## 11096       215  -251.745 9.999e-05      TRUE 0.000153287
## 11097       201  -254.403 9.999e-05      TRUE 0.000153287
## 11098       247  -282.276 9.999e-05      TRUE 0.000153287
## 11099       191  -231.535 9.999e-05      TRUE 0.000153287
## 11100       198  -235.842 9.999e-05      TRUE 0.000153287

6 Demultiplexing hashed libraries

Cell hashing experiments can be demultiplexed using the hashedDrops() function on the set of cell-containing barcode libraries. To demonstrate, we will mock up some hash tag oligo (HTO) counts for a population with cells from each of 10 samples. We will also add some doublets and empty droplets for some flavor: Our first task is to identify the barcodes that actually contain cells. If we already did the calling with emptyDrops() , we could just re-use those calls; otherwise we can obtain calls directly from the HTO count matrix, though this requires some fiddling with lower= to match the sequencing depth of the HTO library. Each cell-containing barcode libary is simply assigned to the sample of origin based on its most abundant HTO. The confidence of the assignment is quantified by the lo

set.seed(10000)

# Simulating empty droplets:
nbarcodes <- 1000
nhto <- 10
y <- matrix(rpois(nbarcodes*nhto, 20), nrow=nhto)

# Simulating cells:
ncells <- 100
true.sample <- sample(nhto, ncells, replace=TRUE)
y[cbind(true.sample, seq_len(ncells))] <- 1000

# Simulating doublets:
ndoub <- ncells/10
next.sample <- (true.sample[1:ndoub]  + 1) %% nrow(y)
next.sample[next.sample==0] <- nrow(y)
y[cbind(next.sample, seq_len(ndoub))] <- 500

hto.calls <- emptyDrops(y, lower=500)
has.cell <- hto.calls$FDR <= 0.001
summary(has.cell)

Resources

• Bioconductor page: https://bioconductor.org/packages/release/bioc/html/DropletUtils.html
• Vignette: https://bioconductor.org/packages/release/bioc/html/DropletUtils.html
• GitHub:
• Bioconductor support forum: https://support.bioconductor.org/
• Paper: Lun ATL et al. (2019). Genome Biology, 20:63.

Run this on BioMate

This skill is the knowledge layer — when, why, and how to use dropletutils. To run this analysis on your own data with managed compute, automated QC, and reproducible outputs, use BioMate — free to start.

▶ Open dropletutils on BioMate →