Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("perturbMatch")

## Or install the development version from GitHub
# install.packages("devtools")
devtools::install_github("soulj/perturbMatch", build_vignettes = TRUE)

Gene-identifier mapping depends on AnnotationDbi and org.Hs.eg.db, and the full reference databases are retrieved with ExperimentHub; all three are installed automatically alongside perturbMatch. To reproduce this vignette you will also need the airway example dataset and DESeq2 for the worked example, and GO.db for the regulator GO enrichment:

BiocManager::install(c("airway", "DESeq2", "GO.db"))

Background

Overview

perturbMatch matches a transcriptomics signature against a reference database of single-gene perturbation signatures, asking what transcriptional regulators could be driving the gene expression changes observed in an experiment.

Reference signatures

perturbMatch provides a large collection of single-gene perturbation signatures derived from publicly available RNA-seq datasets. Each signature captures the transcriptional response to overexpression, knockdown, or knockout of a single gene.

The signatures are computed from 4,453 public RNA-seq series in NCBI GEO, covering 3,571 knockdowns, 1,779 knockouts and 1,452 overexpressions over a common set of 39,376 NCBI genes. Series were selected on the basis of detected single gene perturbations and availability of pre-processed count data. The GEO series behind each signature is recorded in its sig_id and accession, so any hit can be traced back to its source experiment.

Assigning a perturbation type, a target gene and case/control sample groups to each series was automated with an LLM-assisted curation pipeline. The model used is recorded per signature in colData(reference)$model. Automated curation is not error-free, so Verifying hits before interpretation shows how to check a result before acting on it.

The reference data are distributed as the perturbMatchData ExperimentHub package and archived at Zenodo under doi:10.5281/zenodo.21792179, released under a CC-BY-4.0 licence.

How it works

perturbMatch compares your experimental gene expression signature against each reference signature using one of several implemented enrichment methods: SignedJaccard, GSEA, cosine or cosineExtreme. These are described in Scoring methods available. Regulators whose perturbation signatures closely match your query are ranked as likely upstream drivers.

Preparing input data

We start from example RNA-seq data, which is analysed with DESeq2 to identify differentially expressed genes. The differential expression results can be passed directly to perturbMatch to be screened against a large collection of potential regulators.

perturbMatch accepts input as a DESeqResults object, a topTable() data frame, or any data frame with gene identifiers and a log2 fold-change column. prepareQuery() (below) converts any of these into the two-column form calcEnrichment() expects, and also maps the gene identifiers to Entrez IDs used by the reference databases.

Correlation-based methods (see Scoring methods available) require the full set of results from a differential expression analysis, i.e. all genes measured in the experiment, not just the statistically significant subset.

The example below uses the airway dataset, which measures the gene expression response to dexamethasone treatment in airway smooth muscle cells. The airway dataset is loaded into the session with data(airway) after library(airway). We perform differential expression analysis on this data with DESeq2 to create an example gene expression response for querying.

library(perturbMatch)
library(SummarizedExperiment)
library(BiocParallel)
library(airway)
library(DESeq2)

# Load example data
data(airway)

# Calculate differential expression using DESeq2
dds_airway <- DESeqDataSet(airway, design = ~ cell + dex) |>
    DESeq() |>
    results(contrast = c("dex", "trt", "untrt"))

We now have a DESeqResults object, dds_airway, containing the log2 fold changes and adjusted p-values for every gene measured in the experiment. This is the input we will pass to perturbMatch.

head(dds_airway)

## log2 fold change (MLE): dex trt vs untrt 
## Wald test p-value: dex trt vs untrt 
## DataFrame with 6 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE      stat      pvalue
##                  <numeric>      <numeric> <numeric> <numeric>   <numeric>
## ENSG00000000003 708.602170     -0.3812540  0.100654 -3.787752 0.000152016
## ENSG00000000005   0.000000             NA        NA        NA          NA
## ENSG00000000419 520.297901      0.2068126  0.112219  1.842943 0.065337292
## ENSG00000000457 237.163037      0.0379204  0.143445  0.264356 0.791505742
## ENSG00000000460  57.932633     -0.0881682  0.287142 -0.307054 0.758801924
## ENSG00000000938   0.318098     -1.3782270  3.499873 -0.393793 0.693733530
##                       padj
##                  <numeric>
## ENSG00000000003 0.00128292
## ENSG00000000005         NA
## ENSG00000000419 0.19646985
## ENSG00000000457 0.91141962
## ENSG00000000460 0.89500478
## ENSG00000000938         NA

Load reference expression data

Select and load the reference gene expression data to compare against your experiment with getReferenceDatabase(). The reference database is returned as a SummarizedExperiment object.

For the purposes of this vignette we use a small bundled demonstration database of 21 signatures (demo = TRUE). The demonstration database should not be used for any real analysis: it is a small subset of the full reference database, chosen to illustrate both this bulk example and the single-cell workflow, and provided for illustration only.

demo_expression_db <- getReferenceDatabase(demo = TRUE)

## 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.

Choosing a full reference database

The full reference databases are retrieved with getReferenceDatabase() (demo = FALSE is the default) from ExperimentHub, where they are registered by the perturbMatchData package. They are between 0.7-3.1 GB, but only need to be downloaded once: ExperimentHub caches them, and metadata() on the returned object records which source they came from. Several builds of the reference are available through the type and batchCorrection arguments:

type = "chrdir" gives signatures summarised with the characteristic-direction method. This is the default.
type = "limma" gives signatures derived from a limma differential expression analysis useful for producing volcano plots of the signatures.

Each type comes in two variants, selected with batchCorrection: the uncorrected signatures (batchCorrection = FALSE, the default) and signatures computed with study or batch effects removed (batchCorrection = TRUE).

# Characteristic-direction reference, without batch-effect correction (default)
ref_chrdir <- getReferenceDatabase(type = "chrdir")

# Characteristic-direction reference, with batch-effect correction
ref_chrdir_bc <- getReferenceDatabase(type = "chrdir", batchCorrection = TRUE)

# limma reference, without batch-effect correction (default)
ref_limma <- getReferenceDatabase(type = "limma")

# limma reference, with batch-effect correction
ref_limma_sva <- getReferenceDatabase(type = "limma", batchCorrection = TRUE)

chrdir references store a single characteristic-direction score, while limma references store both the log2 fold change (log2fc) and a combined pi score (log2FC * -log10(padj), the default), which weights each gene’s fold change by its statistical significance. calcEnrichment()’s refMetric argument selects which of a reference database’s stored metrics to score against (for example refMetric = "log2fc" on a limma reference). The bundled demonstration database used throughout this vignette is a small limma subset, so its default metric is pi.

Data are not loaded into R until required, keeping RAM usage low. The first few rows and columns of the assay can be inspected:

assay(demo_expression_db)[1:3, 1:3]

## <3 x 3> DelayedMatrix object of type "double":
##           GSE200312_NR3C1_KD_1 GSE197354_AFF4_KD_1 GSE125167_ZBTB16_KD_1
## 643837               0.2355057           0.4235750             0.1854626
## 100130417           -0.2795588          -3.6476560            -1.0390347
## 107985728            0.6614488          -1.6465162            -0.5283873

Metadata for each perturbation experiment are stored in the colData slot:

coldata_cols <- c(
    "sig_id", "pert_symbol", "pert_type", "cell_line", "cell_line_efo_label"
)
head(colData(demo_expression_db)[, coldata_cols], 3)

## DataFrame with 3 rows and 5 columns
##                                      sig_id pert_symbol   pert_type
##                                 <character> <character> <character>
## GSE200312_NR3C1_KD_1   GSE200312_NR3C1_KD_1       NR3C1          KD
## GSE197354_AFF4_KD_1     GSE197354_AFF4_KD_1        AFF4          KD
## GSE125167_ZBTB16_KD_1 GSE125167_ZBTB16_KD_1      ZBTB16          KD
##                                    cell_line cell_line_efo_label
##                                  <character>         <character>
## GSE200312_NR3C1_KD_1  Amnion epithelial ce..              amnion
## GSE197354_AFF4_KD_1   mesenchymal stem cells                  NA
## GSE125167_ZBTB16_KD_1 Human Mesenchymal St..                  NA

Any SummarizedExperiment object containing a measure of differential expression can be used as a reference database for perturbMatch, provided colData includes at least the pert_symbol, sig_id, and pert_type columns identifying each experiment. The remaining columns in Table 1 are annotation that the packaged references also provide: provenance of the source experiment (accession, control_samples, case_samples, technology), how the signature was curated (model), the biological context it was generated in (cell_line, cell_line_efo_id, cell_line_efo_label), and per-signature quality information (entrez_gene_id, target_gene_percentile).

Table 1: Reference database colData columns
Column Name	Description
sig_id	Unique identifier for the perturbation signature. Matches the signature matrix column names
pert_symbol	Official gene symbol of the perturbed gene target
pert_id	Internal ID combining the symbol and perturbation type
pert_type	Type of perturbation: `KD` (knockdown), `KO` (knockout) or `OE` (overexpression)
accession	NCBI GEO series accession for the source study
control_samples	Comma-separated list of raw sample IDs used as the control
case_samples	Comma-separated list of raw sample IDs used as the case
technology	Profiling platform used to generate the source data, e.g. `RNA-seq`
cell_line	Cell line or primary cell type the perturbation was performed in, as free text reported by the source study, e.g. `T47D`, `iPSC-derived Macrophages`
model	Curation approach used to annotate the signature: `curated` by hand, or the LLM that extracted the annotation
entrez_gene_id	NCBI Entrez ID of the perturbed gene. `NA` when the target has no human Entrez entry, e.g. a viral gene. Use it to filter out hits that cannot be human upstream regulators
cell_line_efo_id	Experimental Factor Ontology (EFO) identifier matched to `cell_line`, e.g. `EFO:0003043`. `NA` where no confident match was found, so not every signature is annotated
cell_line_efo_label	EFO label for `cell_line_efo_id`, giving a standardised name for the free-text `cell_line` value (`NCI-H1299` for `H1299`). Use it to group or select signatures by cell context
target_gene_percentile	Percentile position of the perturbed gene within its own signature, from most down-regulated (0) to most up-regulated (100). It is a quality check on whether the perturbation worked

Map the query to the reference

The airway DESeq2 result uses Ensembl gene IDs, but the reference database uses Entrez IDs. prepareQuery() maps between the two, reporting how much of the query mapped. Passing reference additionally reports how much of the mapped query is actually present in the reference, a check that the gene IDs line up before scoring.

By default prepareQuery() also drops genes with a missing adjusted p-value. DESeq2 leaves padj as NA for genes it filters out before testing (typically because they are too lowly expressed to be informative), and these carry no reliable signal for signature matching. Set dropLowExpressed = FALSE to keep them.

final_res <- prepareQuery(dds_airway, reference = demo_expression_db)

## Taking gene identifiers from the row names.

## Dropped 15451 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 83.5% of query genes to Entrez (15041 / 18018).

## 38.8% of mapped genes are present in the reference database.

head(final_res)

##        gene    log2FC         padj
## 1       247  9.505975 1.253044e-17
## 2      7704  7.352626 2.255365e-40
## 3 101928819  6.327383 7.206312e-19
## 4      3000  5.885114 1.000057e-04
## 5      2847 -5.325904 2.378691e-04
## 6     79689  5.207158 5.055589e-24

Starting from a table or CSV file

prepareQuery() also accepts any data.frame, so you can read a file with read.csv() (or read.delim() for tab-separated files) to be used as input.

The table needs:

a gene identifier for each row, either in the row names or in a column. Ensembl, Entrez or gene symbols are all fine. prepareQuery() infers the type and maps to Entrez.
a log2 fold-change column named log2FoldChange, logFC, or log2FC.
(optional) an adjusted p-value column named padj or adj.P.Val, which enables the low-expression filtering described above.

Columns are matched by name, so the fold-change and p-value columns can sit anywhere in the table.

Anything prepareQuery() does not recognise such as raw p-values, base means, average expression is ignored, and the returned query always has the same gene/log2FC/padj columns.

To keep this example self-contained we first write the airway result out to a temporary CSV file, deliberately putting the columns in an awkward order:

de_table <- as.data.frame(dds_airway)
de_table <- data.frame(
    baseMean = de_table$baseMean,       # extra column, ignored
    padj     = de_table$padj,           # optional adjusted p-value
    gene_id  = rownames(de_table),      # gene ID, not in the first column
    pvalue   = de_table$pvalue,         # extra column, ignored
    log2FC   = de_table$log2FoldChange  # fold-change column
)

csv_path <- tempfile(fileext = ".csv")
write.csv(de_table, csv_path, row.names = FALSE)

Reading it back in and preparing it then gives the same query as before:

my_de <- read.csv(csv_path)
head(my_de)

##      baseMean        padj         gene_id       pvalue      log2FC
## 1 708.6021697 0.001282918 ENSG00000000003 0.0001520163 -0.38125398
## 2   0.0000000          NA ENSG00000000005           NA          NA
## 3 520.2979006 0.196469846 ENSG00000000419 0.0653372915  0.20681260
## 4 237.1630368 0.911419623 ENSG00000000457 0.7915057416  0.03792043
## 5  57.9326331 0.895004783 ENSG00000000460 0.7588019240 -0.08816818
## 6   0.3180984          NA ENSG00000000938 0.6937335303 -1.37822703

query_from_csv <- prepareQuery(my_de, reference = demo_expression_db)

## Taking gene identifiers from column 'gene_id'.

## Dropped 15451 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 83.5% of query genes to Entrez (15041 / 18018).

## 38.8% of mapped genes are present in the reference database.

head(query_from_csv)

##        gene    log2FC         padj
## 1       247  9.505975 1.253044e-17
## 2      7704  7.352626 2.255365e-40
## 3 101928819  6.327383 7.206312e-19
## 4      3000  5.885114 1.000057e-04
## 5      2847 -5.325904 2.378691e-04
## 6     79689  5.207158 5.055589e-24

Calculating regulator scores

With the query mapped to the reference, we can score every reference signature against it. Here we score the full ranked expression matrix using the cosine method, which is appropriate for continuous fold-change data (see Scoring methods available for alternatives suited to gene-set input).

regulatorResults <- calcEnrichment(
    demo_expression_db,
    expressionMatrix = final_res,
    method = "cosine"
)

Results are stored in a SummarizedExperiment-like PerturbMatch object, with enrichment scores in the enrichment assay and perturbation experiment details in rowData. Printing it gives a quick summary of both:

regulatorResults

## class: PerturbMatch 
## dim: 21 2 
## metadata(9): up_regulated down_regulated ... demo source
## assays(1): enrichment
## rownames(21): GSE200312_NR3C1_KD_1 GSE197354_AFF4_KD_1 ...
##   GSE202869_ORF9b_OE_1 GSE202869_S_OE_1
## rowData names(14): sig_id pert_symbol ... cell_line_efo_label
##   target_gene_percentile
## colnames(2): cosine zscore
## colData names(0):
## NOTE: scored against the bundled DEMONSTRATION database (21 signatures), not a full reference
## up_regulated(250): 247, 7704, 101928819, 3000, 79689 ...
## down_regulated(250): 2847, 148145, 1591, 2903, 26045 ...
## background(5834): 247, 7704, 101928819, 3000, 2847 ...
## parameters: topN = 250, topSig = 250

Perturbation details for each reference signature carry through into rowData (showing a few key columns here):

head(rowData(regulatorResults)[, c("pert_symbol", "pert_type", "cell_line")])

## DataFrame with 6 rows and 3 columns
##                          pert_symbol   pert_type              cell_line
##                          <character> <character>            <character>
## GSE200312_NR3C1_KD_1           NR3C1          KD Amnion epithelial ce..
## GSE197354_AFF4_KD_1             AFF4          KD mesenchymal stem cells
## GSE125167_ZBTB16_KD_1         ZBTB16          KD Human Mesenchymal St..
## GSE179474_SMAD4_KD_1           SMAD4          KD human dental mesench..
## GSE197118_TAPT1_KO_1           TAPT1          KO            fibroblasts
## GSE227512_LINC01638_KD_1   LINC01638          KD            hMSC-tert20

The up-/down-regulated and background gene sets actually used for scoring are recorded on the result and can be retrieved with the corresponding accessors.

upRegulated(regulatorResults) |> head()

## [1] "247"       "7704"      "101928819" "3000"      "79689"     "101928858"

downRegulated(regulatorResults) |> head()

## [1] "2847"   "148145" "1591"   "2903"   "26045"  "6563"

background(regulatorResults) |> head()

## [1] "247"       "7704"      "101928819" "3000"      "2847"      "79689"

View the top hits ranked by absolute score with topRegulators():

cosine_cols <- c("pert_symbol", "pert_type", "score", "zscore")
topRegulators(regulatorResults, n = 5)[, cosine_cols]

##                          pert_symbol pert_type      score    zscore
## GSE197354_AFF4_KD_1             AFF4        KD -0.3804310 -2.121121
## GSE125167_ZBTB16_KD_1         ZBTB16        KD -0.2939681 -1.596064
## GSE227512_LINC01638_KD_1   LINC01638        KD  0.2183534  1.515073
## GSE179474_SMAD4_KD_1           SMAD4        KD -0.2108919 -1.091574
## GSE197118_TAPT1_KO_1           TAPT1        KO -0.2088128 -1.078948

Checking a known target

The airway data measure the response to dexamethasone, a glucocorticoid that acts directly through the glucocorticoid receptor, NR3C1. NR3C1 is therefore the expected upstream regulator, and the small demo database was selected to include a NR3C1 knockdown signature. We rank the reference signatures by absolute score and look up where NR3C1 is ranked:

ranked <- topRegulators(regulatorResults, n = nrow(regulatorResults))
nr3c1_rank <- match("NR3C1", ranked$pert_symbol)
nr3c1_rank

## [1] 6

ranked[nr3c1_rank, cosine_cols]

##                      pert_symbol pert_type      score     zscore
## GSE200312_NR3C1_KD_1       NR3C1        KD -0.1875777 -0.9499949

The NR3C1 knockdown ranks near the top (rank 6 of 21). Its raw score is negative: dexamethasone activates NR3C1 while the reference signature knocks it down, so the two run in opposite directions. The Visualise the results section below turns this into an activity-level statement with inferActivity = TRUE.

A note on zscore and reference set size. zscore standardises each method’s raw score across the reference signatures actually scored in that call, so its scale depends on the size and diversity of that background set. This demo database has only 21 signatures, nowhere near enough for these particular z-scores to be reliable. With the full reference database (thousands of signatures) the same scores become much more informative.

Using an explicit gene set

Instead of a ranked matrix, calcEnrichment() also accepts an explicit query as upRegulated/downRegulated/background character vectors of gene IDs. Use this format when you only have a fixed set of significant genes rather than a full ranked differential expression result. If you do have an expressionMatrix, SignedJaccard and GSEA will derive an up/down gene set from it for you (as cosine/cosineExtreme always require), so there is no need to manually build the sets in this case. Here we derive a significant gene set from the same DESeq2 result used above (adjusted p-value < 0.05), and score it with SignedJaccard:

sig_genes <- final_res[!is.na(final_res$padj) & final_res$padj < 0.05, ]

up_genes <- sig_genes$gene[sig_genes$log2FC > 0]
down_genes <- sig_genes$gene[sig_genes$log2FC < 0]
background_genes <- final_res$gene

results_signedjaccard <- calcEnrichment(
    demo_expression_db,
    upRegulated = up_genes,
    downRegulated = down_genes,
    background = background_genes,
    method = "SignedJaccard"
)

result_cols <- c("pert_symbol", "pert_type", "score", "overlap", "zscore")
topRegulators(results_signedjaccard, n = 5)[, result_cols]

##                          pert_symbol pert_type       score overlap    zscore
## GSE227512_LINC01638_KD_1   LINC01638        KD  0.06209347     357  1.860321
## GSE197354_AFF4_KD_1             AFF4        KD -0.05921616     366 -1.608796
## GSE125167_ZBTB16_KD_1         ZBTB16        KD -0.05416928     380 -1.464469
## GSE194080_TBX5_KD_1             TBX5        KD  0.04537609     336  1.382251
## GSE179474_SMAD4_KD_1           SMAD4        KD -0.04368749     305 -1.164719

Visualise the results

plotSimilarity() displays the top positive (mimicking) and negative (reversing) signatures as a lollipop plot:

plotSimilarity(regulatorResults, topN = 20)

The fill colour shows the sign of the similarity score, orange for signatures that mimic the query (positive score) and blue for those that reverse it (negative score). The point shape shows the direction of the underlying perturbation: up-pointing triangles for gain-of-function overexpression (OE) and down-pointing triangles for loss-of-function knockdown/knockout (KD/KO). Keeping these separate lets you read the perturbation type together with the sign of the score. A reversing (blue) knockdown (down-triangle) signature, for instance, implies the query behaves as though that gene’s activity were increased.

Set inferActivity = TRUE to flip the sign of scores for knockdown/ knockout perturbations, so the plot reflects inferred gene activity rather than raw transcriptional similarity. A knockdown/knockout signature that is anti-correlated with the query (a negative raw score) implies the query behaves as though that gene’s activity had increased, the opposite of what a negative score would mean for an overexpression signature.

For example, AFF4_KD has a raw cosine z-score of about -2.12, placing it on the “Reverse (Negative)” side of the plot above: the knockdown signature runs opposite to the dexamethasone response. Because it is a knockdown, inferActivity = TRUE flips its sign to about 2.12 below, moving it to the “Mimic (Positive)” side: the data are consistent with dexamethasone increasing AFF4 activity, since suppressing AFF4 produces roughly the opposite expression pattern to the treatment.

plotSimilarity(regulatorResults, topN = 20, inferActivity = TRUE)

Inspecting a single reference signature

Before trusting a top hit it is worth looking at the reference signature itself: did the perturbed gene actually move in the source data, and does the overall signature look like a real, well-behaved differential expression result? plotVolcano() draws a volcano plot (log2 fold change against -log10 FDR) for any single column of a reference database that stores per-gene fold changes and p-values (i.e. the limma references and the demo database used here).

Continuing the dexamethasone example, we plot the NR3C1 knockdown signature and highlight the perturbed gene itself with selectedPoints. The reference is indexed by Entrez ID, so we relabel the rows with gene symbols first to make the plot readable.

nr3c1_signature <- demo_expression_db
rownames(nr3c1_signature) <- rowData(nr3c1_signature)$Symbol

plotVolcano(
    nr3c1_signature,
    signature         = "GSE200312_NR3C1_KD_1",
    effectAssay       = "log2fc",
    pvalueAssay       = "padj",
    selectedPoints    = "NR3C1",
    truncateInfPvalue = FALSE
)

NR3C1 itself sits well down the left-hand (down-regulated) side of its own knockdown signature, confirming the intended perturbation took effect.

Where the perturbed gene ranks: `target_gene_percentile`

The volcano plot gives the differential expression of the perturbed gene one signature at a time. colData() records the same information for every reference signature in a single column, target_gene_percentile: the percentile position of the perturbed gene within its own signature, ranked from most down-regulated (0) to most up-regulated (100).

The direction of the perturbation sets the expectation. A knockout or knockdown (KO/KD) suppresses its target, so the perturbed gene should sit near the bottom of its signature and score a low percentile. An overexpression (OE) drives its target up, so it should sit near the top and score a high percentile. A value that contradicts the perturbation type is a potential warning sign that the signature is mislabelled or that the perturbation did not take effect.

pct_cols <- c("pert_symbol", "pert_type", "target_gene_percentile")
pct <- as.data.frame(colData(demo_expression_db)[, pct_cols])
pct <- pct[order(pct$target_gene_percentile), ]
pct

##                          pert_symbol pert_type target_gene_percentile
## GSE200312_NR3C1_KD_1           NR3C1        KD              0.3185293
## GSE87831_NUP153_KD_2          NUP153        KD              0.5174644
## GSE163705_IFI16_KD_1           IFI16        KD              0.5358568
## GSE125167_ZBTB16_KD_1         ZBTB16        KD              0.6448374
## GSE194080_TBX5_KD_1             TBX5        KD              0.9338521
## GSE197118_TAPT1_KO_1           TAPT1        KO              1.0705348
## GSE129436_STING1_KO_1         STING1        KO              1.8564296
## GSE179474_SMAD4_KD_1           SMAD4        KD              4.2386244
## GSE197354_AFF4_KD_1             AFF4        KD              7.9325079
## GSE137155_STAT2_KD_1           STAT2        KD             13.9997828
## GSE165963_ESRRG_KO_1           ESRRG        KO             16.2846348
## GSE215822_CCDC88A_KO_2       CCDC88A        KO             54.1382206
## GSE227512_LINC01638_KD_1   LINC01638        KD             83.7233489
## GSE114284_IRF1_KO_1             IRF1        KO             94.9992435
## GSE213710_FOSL2_OE_1           FOSL2        OE             97.7972816
## GSE263122_VGLL4_OE_1           VGLL4        OE             99.9688807
## GSE165910_STING1_OE_1         STING1        OE            100.0000000
## GSE182754_ICP0_OE_1             ICP0        OE                     NA
## GSE182754_E4ORF3_OE_1         E4ORF3        OE                     NA
## GSE202869_ORF9b_OE_1           ORF9b        OE                     NA
## GSE202869_S_OE_1                   S        OE                     NA

The NR3C1 knockdown from the volcano plot above illustrates the low end:

nr3c1_pct <- pct[pct$pert_symbol == "NR3C1", ]
nr3c1_pct

##                      pert_symbol pert_type target_gene_percentile
## GSE200312_NR3C1_KD_1       NR3C1        KD              0.3185293

NR3C1 sits in the bottom 1% of its own knockdown signature (target_gene_percentile of about 0.3), near the bottom as a KD should be, and matching where it appeared on the left of the volcano plot. At the other extreme, the overexpression signatures (OE) place their targets close to the 100th percentile.

Note the CCDC88A is annotated as a knockout, so its target should sit near the bottom, yet it lands around the 54th percentile, the middle of its own signature:

pct[pct$pert_symbol == "CCDC88A", ]

##                        pert_symbol pert_type target_gene_percentile
## GSE215822_CCDC88A_KO_2     CCDC88A        KO               54.13822

Signatures such as this should be treated with caution as the perturbation may have been weak or incomplete, or the signature may be mislabelled. This is why Verifying hits before interpretation is important before any candidate is followed up.

Filtering on the target percentile

target_gene_percentile is an ordinary colData column, so signatures can be dropped with standard SummarizedExperiment subsetting before scoring. Here we require knockdowns to place their target in the bottom quarter of their own signature, leaving the other perturbation types untouched:

target_pct <- colData(demo_expression_db)$target_gene_percentile
pert_types <- colData(demo_expression_db)$pert_type

drop <- pert_types == "KD" & !is.na(target_pct) & target_pct > 25
filtered_db <- demo_expression_db[, !drop]

c(all = ncol(demo_expression_db), filtered = ncol(filtered_db))

##      all filtered 
##       21       20

filtered_db is a reference database and can be passed straight to calcEnrichment().

Note that a CRISPR knockout (KO) won’t necessarily reduce the target gene expression (only the target protein), so non-low percentile is not by itself evidence that the perturbation failed.

Functional enrichment of the regulator hits

Individual top hits are easier to interpret as a group when you ask what biological processes the mimicking and reversing regulators have in common. regulatorGSEA() does this by ranking the reference regulators by their similarity score and performing gene set enrichment (via fgsea) over Gene Ontology terms.

A regulator gene usually appears as several signatures (different cell lines, knockdown versus overexpression). regulatorGSEA() collapses these to one value per gene by flipping the sign of knockdown/knockout scores so everything is on a common inferred-activity scale (as in plotSimilarity(inferActivity = TRUE)), then aggregating each gene’s signatures (aggregate = "extreme" by default, keeping the single strongest signature, or "mean" to average them).

Gene sets are built from org.Hs.eg.db and readable GO term names are added from GO.db when it is installed.

go <- regulatorGSEA(results_signedjaccard, ontology = "BP", minSize = 1)
head(go[, c("term", "NES", "padj", "size", "leadingEdge")])

##                                                        term       NES      padj
## 1 negative regulation of transcription by RNA polymerase II  1.696715 0.7038893
## 2                                              MAPK cascade  1.121541 0.7038893
## 3                     microtubule cytoskeleton organization -1.129590 0.7038893
## 4                                              ossification  1.272467 0.7038893
## 5                              ovarian follicle development  1.121541 0.7038893
## 6                                  blood vessel development -1.216481 0.7038893
##   size          leadingEdge
## 1    4 ZBTB16, SMAD4, NR3C1
## 2    1                SMAD4
## 3    1              CCDC88A
## 4    3 ZBTB16, SMAD4, TAPT1
## 5    1                SMAD4
## 6    1                 TBX5

plotRegulatorGSEA() draws a dot plot from that table, placing each term by its NES:

plotRegulatorGSEA(go, n = 10, pvalueCutoff = 1)

On this 21-signature demo database the enrichment is underpowered and none of the terms are individually significant, so pvalueCutoff = 1 is used purely to render an illustrative plot and minSize is dropped to 1. On a full reference (thousands of regulators) keep the defaults (minSize = 15, pvalueCutoff = 0.1), where the enrichment becomes meaningful and the terms summarise the biology shared by the top regulators.

Scoring methods available

Table 2: Overview of reference expression scoring methods
Method	Input Type	Description
cosine	Continuous	Classic cosine similarity between the full continuous log2FC vectors of the query and reference.
cosineExtreme	Rank/set-based	Cosine similarity restricted to the top N most extreme genes of the reference and matching query features.
GSEA	Continuous/ranked	Gene Set Enrichment Analysis. It evaluates whether the query gene set or ranked list is enriched at the top or bottom of the reference signature distribution.
SignedJaccard	Set-based	Intersection-over-union (Jaccard index) for gene sets, taking into account the direction (sign) of the log2FC regulation.

Comparing scoring methods

The methods above rank the same reference signatures on different notions of similarity, so they broadly agree but need not identify identical top hits. Having already scored the airway/dexamethasone query with cosine and SignedJaccard above, we complete the comparison by also running GSEA (which, like SignedJaccard, takes an explicit gene set) and cosineExtreme (which, like cosine, takes the ranked matrix):

results_gsea <- calcEnrichment(
    demo_expression_db,
    upRegulated = up_genes,
    downRegulated = down_genes,
    background = background_genes,
    method = "GSEA"
)

results_cosineextreme <- calcEnrichment(
    demo_expression_db,
    expressionMatrix = final_res,
    method = "cosineExtreme"
)

Aligning the four result sets by signature and correlating their z-scores shows how much they overlap. Methods within the same family, the continuous cosine/cosineExtreme pair and the set-based SignedJaccard/GSEA pair, tend to track each other closely, while agreement across the two families is typically weaker, because they measure genuinely different things:

sig_ids <- rownames(regulatorResults)

score_compare <- data.frame(
    cosine        = regulatorTable(regulatorResults)$zscore,
    cosineExtreme = regulatorTable(results_cosineextreme)[sig_ids, "zscore"],
    SignedJaccard = regulatorTable(results_signedjaccard)[sig_ids, "zscore"],
    GSEA          = regulatorTable(results_gsea)[sig_ids, "zscore"]
)

round(cor(score_compare), 2)

##               cosine cosineExtreme SignedJaccard GSEA
## cosine          1.00          0.99          0.98 0.95
## cosineExtreme   0.99          1.00          0.96 0.93
## SignedJaccard   0.98          0.96          1.00 0.97
## GSEA            0.95          0.93          0.97 1.00

On this small demo database the correlations are good, but on a full-size reference (thousands of signatures) cosine and SignedJaccard can produce disconcordant rankings. cosine asks whether the whole fold-change profile of the query is aligned with the reference, so every gene contributes in proportion to its magnitude. SignedJaccard/GSEA first reduce each signature to a discrete up/down gene set and then ask how much those sets overlap in the same direction, so they are driven by which genes cross the threshold rather than by how large their changes are. A regulator can therefore rank highly under one family and modestly under the other depending on whether its signal lives in a broad continuous shift or a concentrated set of strongly moved genes.

In practice, cosine/cosineExtreme are a reasonable default when a full ranked signature is available, while SignedJaccard/GSEA are useful when only a discrete up/down gene set can be defined, or when a set-based interpretation (shared genes, rather than correlated magnitudes) is preferred.

Interpreting the score, zscore, and overlap columns

Every perturbMatch result carries at least two columns beyond the method’s own score: zscore , the number of standard deviations a signature’s score sits from the mean across the reference queried (used for ranking with topRegulators() and for plotSimilarity()) and, for SignedJaccard, overlap, the number of query genes that overlapped the reference gene set at all (in either direction). A large |zscore| backed by a small overlap is weaker evidence than the same |zscore| backed by a large one, since it is based on fewer shared genes:

reliable_hits <- regulatorTable(results_signedjaccard)
reliable_hits <- reliable_hits[reliable_hits$overlap >= 50, ]
reliable_hits <- reliable_hits[order(-abs(reliable_hits$zscore)), ]

head(reliable_hits[, result_cols])

##                          pert_symbol pert_type       score overlap    zscore
## GSE227512_LINC01638_KD_1   LINC01638        KD  0.06209347     357  1.860321
## GSE197354_AFF4_KD_1             AFF4        KD -0.05921616     366 -1.608796
## GSE125167_ZBTB16_KD_1         ZBTB16        KD -0.05416928     380 -1.464469
## GSE194080_TBX5_KD_1             TBX5        KD  0.04537609     336  1.382251
## GSE179474_SMAD4_KD_1           SMAD4        KD -0.04368749     305 -1.164719
## GSE215822_CCDC88A_KO_2       CCDC88A        KO  0.03364153     302  1.046675

Tuning topN

SignedJaccard and GSEA need reference gene sets, not just ranked values, so each reference signature’s continuous scores are first discretised into up-/down-regulated sets. By default this is done with topN: the topN highest- and lowest-ranked reference genes become that signature’s up-/down-regulated sets.

sj_small <- calcEnrichment(
    demo_expression_db, upRegulated = up_genes, downRegulated = down_genes,
    background = background_genes, method = "SignedJaccard", topN = 100
)
sj_large <- calcEnrichment(
    demo_expression_db, upRegulated = up_genes, downRegulated = down_genes,
    background = background_genes, method = "SignedJaccard", topN = 1000
)

NR3C1 is the expected regulator of the dexamethasone target, so we find its rank under each setting to see the effect these settings have:

nr3c1_rank_at <- function(res) {
    ranked <- topRegulators(res, n = nrow(res))
    match("NR3C1", ranked$pert_symbol)
}

c(
    topN_100  = nr3c1_rank_at(sj_small),
    topN_1000 = nr3c1_rank_at(sj_large)
)

##  topN_100 topN_1000 
##         4        12

The NR3C1 knockdown ranks 4 out of 21 with topN = 100 but slips to 12 with topN = 1000: widening the reference gene sets dilutes the overlap that placed it near the top. Comparing a few topN settings this way is a reasonable check that a finding is not strongly dependent on the cutoff chosen.

Parallelisation

Scoring the full reference database (thousands of signatures) is dominated by reading and decompressing parts of the reference database, and this work parallelises well. calcEnrichment() accepts a BiocParallel back-end via its BPPARAM argument. Supplying a multi-core parameter shares the reference signatures across workers. On a typical multi-core machine this runs roughly 3 to 4 times faster than the default.

library(BiocParallel)

regulatorResults <- calcEnrichment(
    ref_limma,
    expressionMatrix = final_res,
    method  = "cosine"
)

Work is split into column wise chunks whose size is controlled by chunkSize. By default (chunkSize = NULL) calcEnrichment() picks a chunk size automatically from the number of reference signatures and the number of workers in BPPARAM. Because each worker holds its own block of the reference in memory, more workers means higher peak RAM usage.

Verifying hits before interpretation

The reference database is assembled from public datasets whose metadata (which samples are cases versus controls, and which gene was perturbed) is majority LLM-annotated. This automated curation makes the collection broad and has been benchmarked against hand curated data, but it is not error-free: mislabelled case/control groups or incorrect perturbation targets are likely to occur.

Before interpreting or following up a top hit, manually verify it. For each candidate regulator, check the original data source (the accession, control_samples, and case_samples columns in rowData()) to confirm that the case/control assignment is correct and that the annotated perturbed gene is the one actually manipulated.

Cell context matters. Each reference signature reflects the perturbation of a gene in a particular cell line or cell type (the cell_line column of rowData()), and a regulator’s downstream program can differ substantially between contexts. A genuine driver may therefore score only modestly if its available reference signatures come from cell types unrelated to your experiment, while signatures generated in a matching context tend to rank more strongly. Consider how closely its cell_line resembles your own system, and potentially give more weight to regulators supported by several signatures across relevant contexts than to those resting on a single signature from an unrelated cell type. Because cell_line is free text taken from the source study, the same cell line can appear under several spellings. Where an ontology match was available, the cell_line_efo_id and cell_line_efo_label columns give a standardised Experimental Factor Ontology name for it, which is the more reliable field to group or filter signatures on.

Saving results to a file

Any of the tables above is an ordinary data.frame, so it can be written out with write.table() for use in Excel or a downstream script. To save the ranked regulators as a tab-separated (TSV) file:

write.table(
    topRegulators(regulatorResults, n = nrow(regulatorResults)),
    file = "perturbmatch_regulators.tsv",
    sep = "\t", quote = FALSE, row.names = FALSE
)

Related packages

perturbMatch screens a query against a curated collection of single-gene genetic perturbation RNA-seq signatures to nominate upstream regulators. Several other Bioconductor packages address related questions (Table 3). In particular, connectivity against drug/compound signatures (from the L1000 platform) is covered by signatureSearch rather than perturbMatch.

Table 3: Related Bioconductor packages
Package	Focus	Contrast with perturbMatch
signatureSearch	Drug/compound connectivity (LINCS L1000, CMap) for drug repurposing	perturbMatch queries genetic perturbations, not chemical ones
decoupleR / dorothea	Transcription-factor activity from curated regulatory networks	Knowledge/regulon-based, whereas perturbMatch is data-driven
viper / RTN	Master-regulator analysis over inferred regulatory networks	Network-inferred regulators, not measured perturbation signatures
gCMAP	General framework for building and querying connectivity maps	perturbMatch includes a curated perturbation reference and scoring methods

Citation

If you use perturbMatch, please cite:

Soul J, Young DA (2026). Automated generation of a gene perturbation transcriptomic atlas using large language models. bioRxiv. doi:10.64898/2026.08.08.743502

citation("perturbMatch")

Session information

This vignette was compiled using:

sessionInfo()

## 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] airway_1.33.2               BiocParallel_1.47.0        
##  [3] scater_1.41.2               ggplot2_4.0.3              
##  [5] singscore_1.33.0            DESeq2_1.53.6              
##  [7] scuttle_1.23.2              muscData_1.27.1            
##  [9] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
## [11] Biobase_2.73.2              GenomicRanges_1.65.4       
## [13] Seqinfo_1.3.2               IRanges_2.47.5             
## [15] S4Vectors_0.51.10           MatrixGenerics_1.25.0      
## [17] matrixStats_1.5.0           ExperimentHub_3.3.2        
## [19] AnnotationHub_4.3.2         BiocFileCache_3.3.0        
## [21] dbplyr_2.6.0                BiocGenerics_0.59.12       
## [23] generics_0.1.4              perturbMatch_0.99.1        
## [25] knitr_1.52                  BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3        sys_3.4.3                
##   [3] jsonlite_2.0.0            magrittr_2.0.5           
##   [5] ggbeeswarm_0.7.3          farver_2.1.2             
##   [7] rmarkdown_2.32            vctrs_0.7.3              
##   [9] DelayedMatrixStats_1.35.0 memoise_2.0.1            
##  [11] htmltools_0.5.9           S4Arrays_1.13.2          
##  [13] BiocBaseUtils_1.15.1      curl_8.0.0               
##  [15] BiocNeighbors_2.7.3       Rhdf5lib_2.1.0           
##  [17] SparseArray_1.13.4        rhdf5_2.57.20            
##  [19] sass_0.4.10               bslib_0.12.0             
##  [21] plyr_1.8.9                httr2_1.3.0              
##  [23] cachem_1.1.0              buildtools_1.0.0         
##  [25] lifecycle_1.0.5           pkgconfig_2.0.3          
##  [27] rsvd_1.0.5                Matrix_1.7-6             
##  [29] R6_2.6.1                  fastmap_1.2.0            
##  [31] digest_0.6.39             AnnotationDbi_1.75.2     
##  [33] RSpectra_0.16-2           irlba_2.4.1              
##  [35] RSQLite_3.53.3            org.Hs.eg.db_3.23.1      
##  [37] beachmat_2.29.3           filelock_1.0.3           
##  [39] labeling_0.4.3            httr_1.4.9               
##  [41] abind_1.4-8               compiler_4.6.1           
##  [43] bit64_4.8.6               withr_3.0.3              
##  [45] S7_0.2.2                  viridis_0.6.5            
##  [47] DBI_1.3.0                 HDF5Array_1.41.3         
##  [49] rappdirs_0.3.4            DelayedArray_0.39.8      
##  [51] tools_4.6.1               vipor_0.4.7              
##  [53] otel_0.2.0                beeswarm_0.4.0           
##  [55] glue_1.8.1                h5mread_1.5.3            
##  [57] rhdf5filters_1.25.4       grid_4.6.1               
##  [59] reshape2_1.4.5            fgsea_1.39.2             
##  [61] gtable_0.3.6              tidyr_1.3.2              
##  [63] data.table_1.18.6.1       BiocSingular_1.29.2      
##  [65] ScaledMatrix_1.21.0       XVector_0.53.0           
##  [67] RcppAnnoy_0.0.23          ggrepel_0.9.8            
##  [69] BiocVersion_3.24.0        pillar_1.11.1            
##  [71] stringr_1.6.0             limma_3.99.1             
##  [73] dplyr_1.2.1               lattice_0.23-1           
##  [75] bit_4.6.0                 annotate_1.91.0          
##  [77] tidyselect_1.2.1          GO.db_3.23.1             
##  [79] locfit_1.5-9.12           maketools_1.3.2          
##  [81] Biostrings_2.81.9         gridExtra_2.3.1          
##  [83] edgeR_4.99.6              xfun_0.61                
##  [85] statmod_1.5.2             stringi_1.8.9            
##  [87] yaml_2.3.12               evaluate_1.0.5           
##  [89] codetools_0.2-20          tibble_3.3.1             
##  [91] BiocManager_1.30.27       graph_1.91.0             
##  [93] cli_3.6.6                 uwot_0.2.5               
##  [95] xtable_1.8-8              jquerylib_0.1.4          
##  [97] Rcpp_1.1.2                png_0.1-9                
##  [99] XML_3.99-0.25             parallel_4.6.1           
## [101] blob_1.3.0                sparseMatrixStats_1.25.0 
## [103] viridisLite_0.4.3         GSEABase_1.75.0          
## [105] scales_1.4.0              purrr_1.2.2              
## [107] crayon_1.5.3              rlang_1.3.0              
## [109] cowplot_1.2.0             fastmatch_1.1-8          
## [111] KEGGREST_1.53.6

Acknowledgements

This package was developed at the University of Liverpool.

Matching a differential expression signature to perturbation references using perturbMatch