Propensity-score matching for high-MOI Perturb-seq differential expression analysis.
perturbmatch provides cell-level propensity score matching to select balanced control (non-targeting control, NTC) cells for each perturbation in Perturb-seq experiments. In high-MOI screens, cells carry variable numbers of guide RNAs, creating confounding between guide multiplicity and technical covariates (library size, mitochondrial fraction, etc.). Rather than discarding multi-perturbed cells, perturbmatch matches them to appropriate controls within each multiplicity class, preserving statistical power while controlling for confounders.
The core workflow:
- Stratify cells by multiplicity class (singlets, doublets, triplets, etc.) based on the number of detected guides per cell.
- Select target and NTC cells for a given perturbation within each class.
- Match target cells to NTC controls using nearest-neighbor propensity score matching (via MatchIt), balancing on covariates such as log-total counts, mitochondrial fraction, and ribosomal fraction.
- Combine matched cells across multiplicity classes into a single
SingleCellExperimentfor downstream differential expression.
For classes with extreme imbalance (many more NTCs than targets), a weighted-GLM propensity pre-filter removes outlier controls before matching. When the weighted GLM fails, the method falls back to balanced subsampling.
Alternatively, RunPerClassSampling() provides a random selection strategy (no propensity matching) as a drop-in replacement, useful as a baseline comparison or when matching is not needed.
Note: perturbmatch is agnostic to the differential expression method. The package handles cell selection and matching only -- it outputs a labeled SingleCellExperiment that can be passed to any DE framework. In practice, we use single-cell limma-voom as shown in the examples below, but any method that accepts a count matrix and design matrix will work.
# Install from GitHub
devtools::install_github("Genentech/perturbmatch")The package ships a small example SingleCellExperiment (MED19 perturbation from Replogle et al. 2022, ~5300 cells x 3000 genes) so you can run this end-to-end after install:
library(perturbmatch)
library(SingleCellExperiment)
# Load bundled example (Replogle K562 CRISPRi, MED19 target + NTC controls)
sce <- readRDS(system.file("extdata", "example_replogle_MED19.rds", package = "perturbmatch"))
# 1. Select target and NTC cells by multiplicity class
target.by.class <- StratifyCellsByClass(
sce,
get_cells_fn = function(s, nmin, nmax) GetTargetCellsReplogle(s, "MED19", nmin, nmax)
)
ntc.by.class <- StratifyCellsByClass(
sce,
get_cells_fn = function(s, nmin, nmax)
GetNtcCellsReplogle(s, nmin, nmax, exclude.gene = "MED19")
)
# 2a. Run per-class propensity score matching (K=5 NTCs per target cell)
result <- RunPerClassMatching(
sce,
target.by.class = target.by.class,
ntc.by.class = ntc.by.class,
match.K = 5,
match.covariates = c("log.total.counts", "frac.mt", "frac.ribo"),
size.col = "sizeFactor"
)
# 2b. OR use random sampling (no matching) as an alternative
result <- RunPerClassSampling(
sce,
target.by.class = target.by.class,
ntc.by.class = ntc.by.class,
match.K = 5,
size.col = "sizeFactor"
)
# 3. result$sce is a labeled SingleCellExperiment ready for any DE method.
# Cells are labeled "zTarget" or "aNTC" in colData(result$sce)$label.
# Pass counts(result$sce) and a design matrix to your preferred DE framework
# (e.g., limma-voom, edgeR, DESeq2).Adapting to your own data: GetTargetCellsReplogle and GetNtcCellsReplogle require two columns in colData(sce): genestr (comma-separated target gene names per cell, e.g. "MED19", "MED19,PARK7", or "non-targeting" for NTC controls) and nguides.fishash (integer guide count per cell). If your demultiplexing tool outputs guide-level assignments, collapse to gene names and count guides per cell. For fishash-style data with a perturbation indicator matrix, use GetTargetCells()/GetNtcPool() instead (see Example 1 below).
| Function | Description |
|---|---|
GetTargetCells() |
Select target cells by guide name and multiplicity range (fishash-style data) |
GetNtcPool() |
Select NTC control cells, excluding cells co-assigned with the target gene |
GetTargetCellsReplogle() |
Target cell selection for Replogle-style data (comma-delimited genestr column) |
GetNtcCellsReplogle() |
NTC selection for Replogle-style data, with restrictive/permissive modes |
StratifyCellsByClass() |
Convenience wrapper to stratify cells into multiplicity classes (singlets, doublets, ..., 5+) |
| Function | Description |
|---|---|
RunPerClassMatching() |
Top-level function: propensity score matching within each multiplicity class, then combine |
RunPerClassSampling() |
Top-level function: random NTC sampling (no matching). Drop-in replacement for RunPerClassMatching() |
MatchClassCells() |
Per-class matching with weighted-GLM propensity filtering |
SampleClassCells() |
Per-class random sampling. Drop-in replacement for MatchClassCells() |
MatchCellsSce() |
Low-level MatchIt wrapper for nearest-neighbor GLM matching |
These utilities help evaluate DE results but do not perform DE themselves. The DE method is your choice.
| Function | Description |
|---|---|
WrangleLimma() |
Extract a tidy data.frame from a limma fit object (convenience parser) |
ComputeMetrics() |
Compute recall, precision, and cosine similarity vs. a reference result |
rms_zscore() |
Root-mean-square z-score from a toptable |
cosine_sim_tt() |
Cosine similarity between two toptables |
| Function | Description |
|---|---|
SetupSce() |
Label cells as target/NTC, run logNormCounts, prepare matching covariates |
BuildMixedCellsForMoi() |
Build mixed singlet+multiplet cell sets for power analysis |
CheckNullCbind() |
NULL-safe cbind for combining SCE pieces |
This example uses bundled pilot DLD1 Perturb-seq data (ATIC_P1 perturbation, MOI 0.1 and MOI 1.0) to match target and NTC cells across MOI conditions. The example ships a pre-computed indicator matrix so CropQuest is not required.
library(perturbmatch)
library(SingleCellExperiment)
# Load bundled example (DLD1 pilot, ATIC_P1 target, 2 MOI conditions)
obj <- readRDS(system.file("extdata", "example_pilot_ATIC_P1.rds", package = "perturbmatch"))
sce <- obj$sce
ind.mat <- obj$ind.mat
# Split into per-MOI SCE objects
sce.lst <- list(
Low = sce[, colData(sce)$MOI == "MOI_0_1"],
High = sce[, colData(sce)$MOI == "MOI_1_0"]
)
# Subset indicator matrix per MOI
ind.lst <- lapply(sce.lst, function(s) ind.mat[, colnames(s)])
# Stratify cells by multiplicity for each MOI condition
target.cells.lst <- lapply(names(sce.lst), function(moi) {
StratifyCellsByClass(sce.lst[[moi]], function(s, nmin, nmax)
GetTargetCells(s, "ATIC_P1", nmin, nmax, softmatch = FALSE, ind.mat = ind.lst[[moi]]))
})
names(target.cells.lst) <- names(sce.lst)
ntc.pool.lst <- lapply(names(sce.lst), function(moi) {
StratifyCellsByClass(sce.lst[[moi]], function(s, nmin, nmax)
GetNtcPool(s, "ATIC_P1", nmin, nmax, softmatch = TRUE, ind.mat = ind.lst[[moi]]))
})
names(ntc.pool.lst) <- names(sce.lst)
# Reference fit: Low MOI singlets only
ref.result <- RunPerClassMatching(
sce.lst[["Low"]],
target.by.class = list(singlets = target.cells.lst[["Low"]]$singlets),
ntc.by.class = list(singlets = ntc.pool.lst[["Low"]]$singlets),
match.K = 5
)
# High MOI: match singlets + doublets
result <- RunPerClassMatching(
sce.lst[["High"]],
target.by.class = list(singlets = target.cells.lst[["High"]]$singlets,
doublets = target.cells.lst[["High"]]$doublets),
ntc.by.class = list(singlets = ntc.pool.lst[["High"]]$singlets,
doublets = ntc.pool.lst[["High"]]$doublets),
match.K = 5
)For Replogle et al. data where perturbation assignments are stored as a comma-delimited genestr column in colData.
library(perturbmatch)
# Select cells using Replogle-style functions
singlet.targets <- GetTargetCellsReplogle(sce, "CTNNB1", n.min = 1, n.max = 1)
multi.targets <- GetTargetCellsReplogle(sce, "CTNNB1", n.min = 2, n.max = 2)
ntc.singlets <- GetNtcCellsReplogle(sce, n.min = 1, n.max = 1,
restrictive = FALSE, exclude.gene = "CTNNB1")
ntc.multi <- GetNtcCellsReplogle(sce, n.min = 2, n.max = 2,
restrictive = FALSE, exclude.gene = "CTNNB1")
# Match with Replogle-specific settings (sizeFactor instead of total.counts)
result <- RunPerClassMatching(
sce,
target.by.class = list(n1 = singlet.targets, n2 = multi.targets),
ntc.by.class = list(n1 = ntc.singlets, n2 = ntc.multi),
match.K = 5,
match.covariates = c("log.total.counts", "frac.mt", "frac.ribo"),
size.col = "sizeFactor"
)
# Downstream: limma-voom DE and evaluation
tt <- WrangleLimma(fit)
metrics <- ComputeMetrics(tt, tt.ref, dat.ref,
lfc.cutoff = 0.1, lfsr.cutoff = 0.15)RunPerClassSampling() is a drop-in replacement for RunPerClassMatching() that randomly selects K NTCs per target cell without propensity matching. This is useful as a baseline comparison to quantify the benefit of matching, or in settings where covariate balance is already adequate.
# Same interface as RunPerClassMatching -- just swap the function name
result <- RunPerClassSampling(
sce,
target.by.class = target.by.class,
ntc.by.class = ntc.by.class,
match.K = 5
)
# The output is identical in structure: result$sce is a labeled SCE
# with the same "zTarget"/"aNTC" labeling, ready for DE analysis.
# result$match.info records match_method = "random_sample" for each class.The processing scripts support switching between matching and random sampling via a --random_sample flag:
# With propensity score matching (default)
Rscript statistical_power_sublib2.R -guidename CTNNB1_P1 -outfgz out.tsv.gz -outrds out.rds
# With random sampling instead
Rscript statistical_power_sublib2.R -guidename CTNNB1_P1 -outfgz out.tsv.gz -outrds out.rds --random_sampleperturbmatch outputs a labeled SingleCellExperiment with target cells (label == "zTarget") and matched/sampled NTC cells (label == "aNTC"). You then pass this to whatever DE framework you prefer. In the paper we use single-cell limma-voom:
library(limma)
library(edgeR)
# Design matrix: treatment effect + covariates
design <- model.matrix(~ 1 + label + nguides.fishash + log.total.counts + frac.mt + frac.ribo,
data = colData(result$sce))
# Fit limma-voom
dge <- DGEList(counts(result$sce))
dge <- calcNormFactors(dge)
v <- voom(dge, design)
fit <- lmFit(v, design)
fit <- eBayes(fit)
# Wrangle results into a tidy data.frame
tt <- WrangleLimma(fit)But you could equally use edgeR, DESeq2, or any other method -- perturbmatch is upstream of the DE step.
When target cells are rare relative to NTCs (common in high-MOI experiments), the propensity score distribution of controls can be poorly overlapping with targets. MatchClassCells() supports optional propensity score filtering via ps_min and ps_max arguments:
result <- RunPerClassMatching(
sce,
target.by.class = target.by.class,
ntc.by.class = ntc.by.class,
match.K = 5,
ps_min = 0.2,
ps_max = 0.8
)The method uses a weighted GLM where sampling weights equalize effective class sizes, making propensity scores comparable on the [0, 1] scale before applying the fixed thresholds.
- Analysis code: https://github.com/Genentech/perturbmatch_analysis -- figure reproduction scripts and processed data for the paper.
If you use perturbmatch, please cite:
Yeung, J., Tan, J., et al. (2026). Joint analysis of multiply perturbed cells improves statistical power and cost efficiency in Perturb-seq. bioRxiv. https://doi.org/10.64898/2026.07.10.737863
MIT