Matching a single-cell RNA-seq response to perturbation references using perturbMatch
library(perturbMatch)
library(muscData)
library(SingleCellExperiment)
library(scuttle)
library(DESeq2)
library(SummarizedExperiment)
library(singscore)
library(scater)
library(ggplot2)Overview
perturbMatch analyses a differential expression signature regardless
of whether it came from bulk or single-cell RNA-seq. For single-cell
data the recommended route is to aggregate counts into pseudobulk
samples and run a standard differential expression analysis. This
vignette follows that route on the Kang et al. PBMC dataset from
muscData, comparing control and interferon-beta-stimulated
CD14+ monocytes.
The bulk RNA-seq vignette (vignette("perturbMatch"))
introduces the reference database, scoring methods, and plotting
functions in more detail. Here we focus on the single-cell workflow and
two points specific to it: filtering hits that cannot be human
regulators, and validating a candidate regulator at the level of
individual cells.
Load the data
The muscData package distributes the Kang et al. PBMC
study as a SingleCellExperiment. Eight lupus patients each
contributed peripheral blood mononuclear cells that were split into two
aliquots, one left untreated and one stimulated with recombinant
interferon beta for six hours, then profiled by droplet single-cell
RNA-seq. The design therefore pairs a control and a stimulated sample
within every donor, which lets us treat donor as a blocking factor later
on. The object carries raw UMI counts in the counts assay,
and per-cell annotation in colData including the donor of
origin (ind), the stimulation condition
(stim), a coarse cell type label (cell), and a
doublet call (multiplets).
We copy the donor identifier into a tidier donor column
and take a first look at how many cells fall in each cell type and in
each donor by condition combination. The counts confirm the paired
control and stimulated design across the eight donors.
##
## B cells CD14+ Monocytes CD4 T cells CD8 T cells
## 2880 6447 12033 2634
## Dendritic cells FCGR3A+ Monocytes Megakaryocytes NK cells
## 472 1914 346 2330##
## ctrl stim
## Donor101 1027 1363
## Donor1015 3138 2703
## Donor1016 2180 1998
## Donor1039 552 797
## Donor107 645 641
## Donor1244 2250 1755
## Donor1256 2451 2275
## Donor1488 2376 2914Select a cell population
We restrict the analysis to a single, well-defined population, the CD14+ monocytes, keeping only confidently assigned singlets.
sce <- sce[, sce$multiplets == "singlet"]
sce <- sce[, !is.na(sce$cell)]
sce <- sce[, sce$cell == "CD14+ Monocytes"]
sce## class: SingleCellExperiment
## dim: 35635 5385
## metadata(0):
## assays(1): counts
## rownames(35635): MIR1302-10 FAM138A ... MT-ND6 MT-CYB
## rowData names(2): ENSEMBL SYMBOL
## colnames(5385): AAACATACATTTCC-1 AAACATACCAGAAA-1 ... TTTGCATGAGGTTC-1
## TTTGCATGCATGAC-1
## colData names(6): ind stim ... multiplets donor
## reducedDimNames(1): TSNE
## mainExpName: NULL
## altExpNames(0):Aggregate into pseudobulk samples
aggregateAcrossCells() sums counts within each donor and
stimulation condition, giving one pseudobulk profile per biological
replicate.
pb <- aggregateAcrossCells(
sce,
ids = DataFrame(donor = sce$donor, stim = sce$stim)
)
colData(pb)[, c("donor", "stim", "ncells")]## DataFrame with 16 rows and 3 columns
## donor stim ncells
## <character> <factor> <integer>
## 1 Donor101 ctrl 202
## 2 Donor101 stim 253
## 3 Donor1015 ctrl 779
## 4 Donor1015 stim 671
## 5 Donor1016 ctrl 369
## ... ... ... ...
## 12 Donor1244 stim 290
## 13 Donor1256 ctrl 390
## 14 Donor1256 stim 375
## 15 Donor1488 ctrl 295
## 16 Donor1488 stim 348Differential expression
We test for genes responding to interferon stimulation with DESeq2, controlling for donor.
dds <- DESeqDataSet(pb, design = ~ donor + stim)
dds <- DESeq(dds)
res <- results(dds, contrast = c("stim", "stim", "ctrl"))
res |> as.data.frame() |>
dplyr::filter(!is.na(padj)) |>
dplyr::arrange(padj) |>
head(5)## baseMean log2FoldChange lfcSE stat pvalue padj
## ISG20 3031.3359 5.823768 0.2223578 26.19097 3.367411e-151 3.019220e-147
## IL1RN 885.1381 6.692157 0.2656365 25.19291 4.790100e-140 1.431601e-136
## DYNLT1 915.3704 2.759997 0.1095206 25.20070 3.935577e-140 1.431601e-136
## SSB 359.2464 3.429619 0.1400195 24.49386 1.717141e-132 3.848972e-129
## IFITM2 472.9234 4.949523 0.2027138 24.41631 1.147677e-131 2.058014e-128Score candidate regulators
Load the demonstration reference database and map the query onto it.
As in the bulk workflow, prepareQuery() maps gene
identifiers to Entrez and drops genes with a missing adjusted p-value,
which DESeq2 assigns to genes it filtered out before testing.
## Loaded the bundled DEMONSTRATION database (6500 genes x 21 signatures). It is a small illustrative subset of the limma reference and should not be used for real analysis. Use getReferenceDatabase(demo = FALSE) for a full reference database.## Taking gene identifiers from the row names.## Dropped 5854 genes with a missing adjusted p-value (not tested, typically lowly expressed). Set dropLowExpressed = FALSE to keep them.## 'select()' returned 1:many mapping between keys and columns## Mapped 85.0% of query genes to Entrez (7624 / 8966).## 51.3% of mapped genes are present in the reference database.## gene log2FC padj
## 1 26155 -1.1137472 9.713894e-03
## 2 84069 1.8657621 1.391709e-02
## 3 57801 7.6632728 3.677165e-26
## 4 9636 7.7426156 9.789916e-62
## 5 375790 1.8968040 5.284328e-02
## 6 8784 -0.9188896 2.878079e-01The stimulated CD14+ monocytes carry a continuous fold-change
signature, so we score it with cosine.
regulatorResults <- calcEnrichment(
demo_expression_db,
expressionMatrix = query,
method = "cosine"
)
hit_cols <- c("pert_symbol", "pert_type", "entrez_gene_id", "score", "zscore")
topRegulators(regulatorResults, n = 10)[, hit_cols]## pert_symbol pert_type entrez_gene_id score
## GSE182754_E4ORF3_OE_1 E4ORF3 OE NA 0.61481153
## GSE182754_ICP0_OE_1 ICP0 OE NA 0.56136905
## GSE129436_STING1_KO_1 STING1 KO 340061 -0.47116149
## GSE163705_IFI16_KD_1 IFI16 KD 3428 -0.40518807
## GSE213710_FOSL2_OE_1 FOSL2 OE 2355 0.16303436
## GSE263122_VGLL4_OE_1 VGLL4 OE 9686 0.11372127
## GSE114284_IRF1_KO_1 IRF1 KO 3659 -0.09139626
## GSE200312_NR3C1_KD_1 NR3C1 KD 2908 0.06928269
## GSE215822_CCDC88A_KO_2 CCDC88A KO 55704 -0.06873941
## GSE194080_TBX5_KD_1 TBX5 KD 6910 0.06613516
## zscore
## GSE182754_E4ORF3_OE_1 2.4843252
## GSE182754_ICP0_OE_1 2.2607199
## GSE129436_STING1_KO_1 -2.0594238
## GSE163705_IFI16_KD_1 -1.7833887
## GSE213710_FOSL2_OE_1 0.5940737
## GSE263122_VGLL4_OE_1 0.3877461
## GSE114284_IRF1_KO_1 -0.4704728
## GSE200312_NR3C1_KD_1 0.2018135
## GSE215822_CCDC88A_KO_2 -0.3756758
## GSE194080_TBX5_KD_1 0.1886441Filtering hits that cannot be human regulators
The reference collection is assembled from many public perturbation studies, and the sample annotation was largely produced by an LLM-assisted pipeline. That pipeline has erroneously admitted some studies that overexpress viral proteins in human cells to study host-pathogen interactions, and labelled the viral protein as the perturbed target. These signatures make biological sense as top hits for an antiviral interferon response, but the perturbed factor is not a human gene and so cannot be an upstream regulator in the intended sense.
We do not have to spot these by eye. Because the annotation could not
assign a human Entrez ID to a viral protein, the
entrez_gene_id column is NA for exactly these
cases, so we can filter them out explicitly:
scored <- regulatorTable(regulatorResults)
nonhuman <- scored[is.na(scored$entrez_gene_id), ]
nonhuman[, c("pert_symbol", "pert_type", "score", "zscore")]## pert_symbol pert_type score zscore
## GSE182754_ICP0_OE_1 ICP0 OE 0.5613690474 2.26071994
## GSE182754_E4ORF3_OE_1 E4ORF3 OE 0.6148115325 2.48432516
## GSE202869_ORF9b_OE_1 ORF9b OE 0.0005821846 -0.08563181
## GSE202869_S_OE_1 S OE -0.0086188085 -0.12412908Here the two strongest positive hits are the adenoviral protein E4ORF3 and the herpesviral protein ICP0, both known modulators of innate immune signalling. Keeping only the rows with a human Entrez ID leaves the top human candidate regulators:
human_hits <- scored[!is.na(scored$entrez_gene_id), ]
human_hits <- human_hits[order(-abs(human_hits$zscore)), ]
head(human_hits[, c("pert_symbol", "pert_type", "score", "zscore")])## pert_symbol pert_type score zscore
## GSE129436_STING1_KO_1 STING1 KO -0.47116149 -2.0594238
## GSE163705_IFI16_KD_1 IFI16 KD -0.40518807 -1.7833887
## GSE213710_FOSL2_OE_1 FOSL2 OE 0.16303436 0.5940737
## GSE114284_IRF1_KO_1 IRF1 KO -0.09139626 -0.4704728
## GSE263122_VGLL4_OE_1 VGLL4 OE 0.11372127 0.3877461
## GSE215822_CCDC88A_KO_2 CCDC88A KO -0.06873941 -0.3756758The strongest human hits are the DNA-sensing and interferon machinery: STING1 (cGAS-STING pathway), IFI16, and IRF1. Their reference signatures are knockouts or knockdowns whose scores are negative (reversing): removing an antiviral regulator produces the opposite of the interferon response, exactly the relationship expected if that regulator is required for the response.
We can filter at the PerturbMatch object level also:
Visualise the results
Setting inferActivity = TRUE flips the sign of the
knockdown/knockout signatures, so the antiviral regulators move to the
mimicking side: the data are consistent with interferon stimulation
increasing the activity of STING1, IFI16, and
IRF1.
Cell-level scoring with singscore
perturbMatch ranks regulators from the aggregated pseudobulk signature, which collapses every cell in a sample to a single value. This hides any variation between cells. singscore scores each individual cell against a gene signature, so we can bring a candidate regulator back down to the single-cell level and look at how its activity is distributed across the population rather than assuming it is uniform.
The reference contains both a STING1 overexpression and a STING1 knockout signature. They sit on opposite sides of the raw perturbMatch score, the overexpression mimicking the interferon response and the knockout reversing it, but the two are concordant once the direction of each perturbation is taken into account: both say that interferon stimulation raises STING1 activity. We use the knockout here because it is the much stronger, higher-confidence hit.
Because the signature is a knockout, its direction is inverted relative to STING1 activity. A knockout lowers its target’s pathway, so the genes it pushes up mark loss of STING1 signalling and the genes it pushes down mark active signalling. To score STING1 activity we therefore reverse the two sets: the knockout’s down-genes become the activity up-set and its up-genes become the down-set. Because we relabelled nothing, these are gene symbols, matching the row names of the single-cell object.
sting_col <- "GSE129436_STING1_KO_1"
sting_lfc <- as.numeric(assay(demo_expression_db, "log2fc")[, sting_col])
symbols <- rowData(demo_expression_db)$Symbol
ord <- order(sting_lfc)
n <- 100
ko_down <- symbols[head(ord, n)] # most decreased by the STING1 knockout
ko_up <- symbols[tail(ord, n)] # most increased by the STING1 knockout
# Reverse the knockout to score STING1 *activity*: genes that fall when
# STING1 is knocked out are the ones that mark active STING1 signalling.
sting_up <- ko_down
sting_down <- ko_upsingscore ranks all genes within each cell, so we work on log-normalised counts and restrict to reasonably detected genes to keep the matrix compact. The result is one STING1 activity score per cell.
sce <- logNormCounts(sce)
sce <- sce[rowSums(counts(sce) > 0) >= 10, ]
sting_up <- intersect(sting_up, rownames(sce))
sting_down <- intersect(sting_down, rownames(sce))
ranked <- rankGenes(as.matrix(logcounts(sce)))
sce$STING1_activity <- simpleScore(
ranked, upSet = sting_up, downSet = sting_down
)$TotalScoreTo see how that activity is spread across the population, we embed the cells with a UMAP. We first reduce the data with PCA and then run UMAP on the principal components.
set.seed(1)
sce <- runPCA(sce, ncomponents = 20)
sce <- runUMAP(sce, dimred = "PCA")
umap_df <- data.frame(
reducedDim(sce, "UMAP"),
STING1_activity = sce$STING1_activity,
stim = sce$stim
)
colnames(umap_df)[1:2] <- c("UMAP1", "UMAP2")Colouring one point per cell is unreadable with several thousand overlapping cells, so we bin the UMAP into hexagons and colour each hexagon by the mean STING1 activity score of the cells in it. A single-hue scale from near-white to dark blue, with the limits clipped so a few extreme cells cannot dominate it, makes the gradient across the population clear.
p_activity <- ggplot(umap_df, aes(UMAP1, UMAP2, z = STING1_activity)) +
stat_summary_hex(bins = 50, fun = mean) +
scale_fill_gradient(
low = "#EAF3FB", high = "#08519C",
limits = c(0, 0.5), oob = scales::squish,
name = "Mean\nSTING1\nactivity"
) +
labs(title = "STING1 activity") +
theme_bw()
p_condition <- ggplot(umap_df, aes(UMAP1, UMAP2, colour = stim)) +
geom_point(size = 0.35, alpha = 0.6) +
scale_colour_manual(values = c(ctrl = "#4C72B0", stim = "#DD8452")) +
guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))) +
labs(title = "Condition") +
theme_bw()
cowplot::plot_grid(p_activity, p_condition, nrow = 1)
The stimulated and control monocytes separate into two clusters, and the STING1 activity is raised across the stimulated cluster, confirming at single-cell resolution the regulator that perturbMatch flagged from the pseudobulk signature. The binned view also exposes heterogeneity that the pseudobulk value hides: the score is graded rather than uniform, so some regions of the stimulated cluster carry a stronger STING1 programme than others.
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scater_1.41.2 ggplot2_4.0.3
## [3] singscore_1.33.0 DESeq2_1.53.6
## [5] scuttle_1.23.2 muscData_1.27.1
## [7] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
## [9] Biobase_2.73.2 GenomicRanges_1.65.4
## [11] Seqinfo_1.3.2 IRanges_2.47.5
## [13] S4Vectors_0.51.10 MatrixGenerics_1.25.0
## [15] matrixStats_1.5.0 ExperimentHub_3.3.2
## [17] AnnotationHub_4.3.2 BiocFileCache_3.3.0
## [19] dbplyr_2.6.0 BiocGenerics_0.59.12
## [21] generics_0.1.4 perturbMatch_0.99.1
## [23] knitr_1.52 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_2.0.0
## [4] magrittr_2.0.5 ggbeeswarm_0.7.3 farver_2.1.2
## [7] rmarkdown_2.32 vctrs_0.7.3 memoise_2.0.1
## [10] htmltools_0.5.9 S4Arrays_1.13.2 BiocBaseUtils_1.15.1
## [13] curl_8.0.0 BiocNeighbors_2.7.3 Rhdf5lib_2.1.0
## [16] SparseArray_1.13.4 rhdf5_2.57.20 sass_0.4.10
## [19] bslib_0.12.0 plyr_1.8.9 httr2_1.3.0
## [22] cachem_1.1.0 buildtools_1.0.0 lifecycle_1.0.5
## [25] pkgconfig_2.0.3 rsvd_1.0.5 Matrix_1.7-6
## [28] R6_2.6.1 fastmap_1.2.0 digest_0.6.39
## [31] AnnotationDbi_1.75.2 RSpectra_0.16-2 irlba_2.4.1
## [34] RSQLite_3.53.3 org.Hs.eg.db_3.23.1 beachmat_2.29.3
## [37] filelock_1.0.3 labeling_0.4.3 httr_1.4.9
## [40] abind_1.4-8 compiler_4.6.1 bit64_4.8.6
## [43] withr_3.0.3 S7_0.2.2 BiocParallel_1.47.0
## [46] viridis_0.6.5 DBI_1.3.0 HDF5Array_1.41.3
## [49] rappdirs_0.3.4 DelayedArray_0.39.8 tools_4.6.1
## [52] vipor_0.4.7 otel_0.2.0 beeswarm_0.4.0
## [55] glue_1.8.1 h5mread_1.5.3 rhdf5filters_1.25.4
## [58] grid_4.6.1 reshape2_1.4.5 fgsea_1.39.2
## [61] gtable_0.3.6 tidyr_1.3.2 data.table_1.18.6.1
## [64] BiocSingular_1.29.2 ScaledMatrix_1.21.0 XVector_0.53.0
## [67] RcppAnnoy_0.0.23 ggrepel_0.9.8 BiocVersion_3.24.0
## [70] pillar_1.11.1 stringr_1.6.0 limma_3.99.1
## [73] dplyr_1.2.1 lattice_0.23-1 bit_4.6.0
## [76] annotate_1.91.0 tidyselect_1.2.1 locfit_1.5-9.12
## [79] maketools_1.3.2 Biostrings_2.81.9 gridExtra_2.3.1
## [82] edgeR_4.99.6 xfun_0.61 statmod_1.5.2
## [85] stringi_1.8.9 yaml_2.3.12 evaluate_1.0.5
## [88] codetools_0.2-20 tibble_3.3.1 BiocManager_1.30.27
## [91] graph_1.91.0 cli_3.6.6 uwot_0.2.5
## [94] xtable_1.8-8 jquerylib_0.1.4 Rcpp_1.1.2
## [97] png_0.1-9 XML_3.99-0.25 parallel_4.6.1
## [100] blob_1.3.0 viridisLite_0.4.3 GSEABase_1.75.0
## [103] scales_1.4.0 purrr_1.2.2 crayon_1.5.3
## [106] rlang_1.3.0 cowplot_1.2.0 fastmatch_1.1-8
## [109] KEGGREST_1.53.6