Package 'perturbMatch'

Background genes from a PerturbMatch object

Description

Returns the background gene set recorded in metadata() when the object was built by calcEnrichment. This is the universe of genes against which enrichment was scored.

Usage

background(object)

Arguments

object

A PerturbMatch object.

Value

A character vector of gene identifiers, or NULL if unset.

Examples

pm <- methods::new(
    "PerturbMatch",
    SummarizedExperiment::SummarizedExperiment(
        assays = list(
            enrichment = matrix(0.5, dimnames = list("SIG1", "score"))
        ),
        metadata = list(background = c("GENE1", "GENE2", "GENE3"))))
background(pm)

---

Calculate enrichment

Description

Calculates the enrichment of a regulator expression signature in the query expression signature. This measure is used to prioritise which upstream regulators are driving the observed gene expression changes.

Usage

calcEnrichment(
  object,
  expressionMatrix = NULL,
  upRegulated = NULL,
  downRegulated = NULL,
  background = NULL,
  method = c("SignedJaccard", "GSEA", "cosine", "cosineExtreme"),
  refMetric = NULL,
  threshold = NULL,
  topN = 250L,
  topSig = 250L,
  BPPARAM = BiocParallel::SerialParam(progressbar = TRUE),
  chunkSize = NULL
)

Arguments

object	A SummarizedExperiment reference database with metadata(object)$ref_metric set to one of "chrdir", "log2fc", or "pi".
expressionMatrix	Optional data frame with gene IDs in column 1 and score values in column 2. Required for "cosine" and "cosineExtreme" methods. If supplied, upRegulated, downRegulated, and background are not used.
upRegulated	Character vector of upregulated genes in the experiment. Only used when expressionMatrix is NULL.
downRegulated	Character vector of downregulated genes in the experiment. Only used when expressionMatrix is NULL.
background	Character vector of all genes expressed/measured in the experiment. Only used when expressionMatrix is NULL.
method	Enrichment method to use. One of "SignedJaccard", "GSEA", "cosine", "cosineExtreme". Defaults to "SignedJaccard".
refMetric	Reference metric to score against, one of metadata(object)$ref_metric. When NULL (default), metadata(object)$default_metric is used.
threshold	Optional numeric threshold applied to reference metric values to define reference gene sets (|metric| > threshold). When NULL, topN is used instead. Default: NULL.
topN	Number of top and bottom genes per reference signature to use, ranked by the reference metric. Used only when threshold is NULL. Default: 250.
topSig	Number of top/bottom genes from the query signature to use when deriving gene sets from expressionMatrix. Default: 250.
BPPARAM	A BiocParallelParam object controlling parallelisation. Defaults to BiocParallel::SerialParam(progressbar = TRUE).
chunkSize	Number of reference signatures (columns) per processing block. When NULL (default) a value giving roughly one block per worker in BPPARAM is chosen. Larger values use more memory, and chunkSize >= ncol(object) leaves a single block that cannot be parallelised.

Value

A PerturbMatch object containing the enrichment scores.

Input signature

Two input formats are supported:

• Ranked list: supply expressionMatrix (a data frame of gene IDs and score values). upRegulated, downRegulated, and background are ignored.

• Gene sets: supply upRegulated, downRegulated, and background as character vectors. expressionMatrix should be NULL.

If both are provided, expressionMatrix takes precedence.

Reference metric

Two metadata slots record which metrics a reference supports:

• metadata(object)$ref_metric: the available metrics, each naming a stored assay. limma references hold c("pi", "log2fc") and chrdir references hold "chrdir".

• metadata(object)$default_metric: the metric used when refMetric = NULL.

Examples

# Load the small bundled demo reference database
se <- getReferenceDatabase(demo = TRUE)

# Score a query gene set against every reference signature. Here we just
# use a handful of reference genes to stand in for a real up/down
# signature.
result <- calcEnrichment(
    se,
    upRegulated   = rownames(se)[seq_len(3)],
    downRegulated = rownames(se)[seq(4, 5)],
    background    = rownames(se),
    method        = "SignedJaccard",
    BPPARAM       = BiocParallel::SerialParam(progressbar = FALSE)
)

cols <- c("pert_symbol", "pert_type", "score", "zscore")
topRegulators(result, n = 3)[, cols]

---

Down-regulated genes from a PerturbMatch object

Description

Returns the down-regulated gene set recorded in metadata() when the object was built by calcEnrichment.

Usage

downRegulated(object)

Arguments

object

A PerturbMatch object.

Value

A character vector of gene identifiers, or NULL if unset.

Examples

pm <- methods::new(
    "PerturbMatch",
    SummarizedExperiment::SummarizedExperiment(
        assays = list(
            enrichment = matrix(0.5, dimnames = list("SIG1", "score"))
        ),
        metadata = list(down_regulated = c("GENE3", "GENE4"),
            background = c("GENE3", "GENE4"))))
downRegulated(pm)

---

Enrichment scores from a PerturbMatch object

Description

Returns the matrix stored in the "enrichment" assay.

Usage

enrichment(object)

Arguments

object

A PerturbMatch object.

Value

A matrix of enrichment scores (rows = features, columns = signatures).

Examples

mat <- matrix(runif(3), nrow = 3, dimnames = list(
    c("GENE1", "GENE2", "GENE3"), "score"
))
se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(enrichment = mat)
)
S4Vectors::metadata(se) <- list(
    up_regulated   = "GENE1",
    down_regulated = "GENE3",
    background     = c("GENE1", "GENE2", "GENE3")
)
pm <- new("PerturbMatch", se)
enrichment(pm)

---

Load the perturbMatch reference database

Description

Retrieve a perturbMatch reference expression database as a [SummarizedExperiment::SummarizedExperiment-class] object.

Usage

getReferenceDatabase(
  demo = FALSE,
  type = c("chrdir", "limma"),
  batchCorrection = FALSE,
  source = c("auto", "experimenthub", "zenodo"),
  cacheDir = tools::R_user_dir("perturbMatch", "cache"),
  ask = interactive()
)

Arguments

demo	Logical scalar indicating whether to load the small bundled demonstration dataset instead of a full reference. Defaults to FALSE. The demo ignores type and batchCorrection. See Details for what it is and is not good for.
type	Character scalar specifying which reference database to load. One of "chrdir" (characteristic direction) or "limma" (log2 fold-change and adjusted p-values). Defaults to "chrdir".
batchCorrection	Logical scalar selecting the batch-corrected variant of a reference. Available for both type = "chrdir" and type = "limma": FALSE (default) loads the uncorrected signatures and TRUE loads signatures computed with study/batch effects removed (SVA).
source	Character scalar selecting where the full reference is retrieved from. "experimenthub" takes it from Bioconductor's ExperimentHub and errors if the records are not registered there. "zenodo" downloads it from the Zenodo deposit instead. "auto" (the default) uses ExperimentHub, falling back to Zenodo only when the records are absent from the hub, and announces which of the two it used. The records are registered for Bioconductor 3.24 and later, so an older release reaches the data through the Zenodo fallback. Ignored when demo = TRUE.
cacheDir	Directory holding the BiocFileCache cache of downloaded reference files. Defaults to a per-user cache directory (tools::R_user_dir("perturbMatch", "cache")). Only used by the Zenodo download path, as ExperimentHub manages its own cache.
ask	Logical scalar controlling whether BiocFileCache asks before creating cacheDir. Defaults to interactive(), so a session with a user present is prompted and a scripted one creates the directory. Declining the prompt, or passing TRUE outside an interactive session, puts the cache under tempdir() and so re-downloads the reference in every session.

Details

Each reference is hosted as three files: an HDF5 file of the assay matrices (0.7 GB for chrdir_sva up to 3.1 GB for limma), a serialised colData of the signature annotation, and one gene annotation shared by all four references. They are fetched once, cached, and reassembled into a SummarizedExperiment whose assays are read from disk on demand.

The resources are registered on ExperimentHub, and documented by the perturbMatchData package, which does not need to be installed to use them. The signatures are derived from public NCBI GEO RNA-seq series and released under CC-BY-4.0, doi:10.5281/zenodo.21792179.

The bundled demo mirrors the structure of the full databases so the package can be used offline. It is a 21-signature subset tailored to the vignette examples. The rankings it produces are illustrative only. Any real analysis needs a full reference (demo = FALSE, the default).

Value

A [SummarizedExperiment::SummarizedExperiment-class] object with assay(s) stored as HDF5-backed delayed arrays. The object carries metadata(object)$ref_metric, which tells calcEnrichment how to compute the scoring matrix from the assay(s):

"chrdir"

Single assay named "chrdir", used directly.

"log2fc"

Assay named "log2fc", used directly.

"pi"

Combined $\pi = \log_2FC \times -\log_{10}(padj)$ score stored alongside the "log2fc" and "padj" assays.

ref_metric is fixed when the reference database is built, travels with the data, and does not need to be set or modified by the user. Alongside it are method, batch_correction, default_metric and the build_version and source_doi of the deposit the reference came from. metadata(object)$source records how it was retrieved ("experimenthub", "zenodo" or "bundled"), and metadata(object)$demo whether it is the bundled demo database. All are carried through to the PerturbMatch object returned by calcEnrichment.

Examples

# Load the bundled demonstration database (no download required)
se <- getReferenceDatabase(demo = TRUE)
se
S4Vectors::metadata(se)$ref_metric

# A full reference is a multi-gigabyte one-time download, so it is fetched
# only in an interactive session
if (interactive()) {
    seChrdir <- getReferenceDatabase(type = "chrdir")

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

    # limma reference, without batch-effect correction
    seLimma <- getReferenceDatabase(type = "limma", batchCorrection = FALSE)
}

---

PerturbMatch enrichment results

Description

An S4 class extending SummarizedExperiment that holds the enrichment scores produced by calcEnrichment. Scores are available through enrichment(), or as a tidy table through regulatorTable(). The usual SummarizedExperiment accessors also work.

Details of analysis are stored in metadata() when the object is built. The query gene sets have accessors of their own (upRegulated, downRegulated, background). The scoring parameters are read straight from metadata() and are printed by show(). The reference gene sets come from either threshold or topN.

The class contains an "enrichment" assay of numeric scores, gene sets held as character vectors, and a score_column entry. Objects failing any of these are rejected by validObject.

Arguments

object

A PerturbMatch object returned by calcEnrichment.

Value

• enrichment: matrix of enrichment scores.

• regulatorTable: tidy data.frame of scores plus row metadata.

---

Dot plot of regulator GO enrichment

Description

Draws a dot plot of the GO terms returned by regulatorGSEA. Terms are placed by normalised enrichment score, so those enriched among mimicking regulators (positive NES) sit to the right of zero and those enriched among reversing regulators (negative NES) to the left. Point size shows the GO set size and colour the adjusted p-value.

Usage

plotRegulatorGSEA(
  x,
  n = 15L,
  pvalueCutoff = 0.1,
  wrapWidth = 40,
  fontSize = 12
)

Arguments

x	A data.frame returned by regulatorGSEA.
n	Maximum number of terms to show. Taken by largest |NES| so both directions are represented. Default 15.
pvalueCutoff	Only terms with padj <= pvalueCutoff are shown. Default 0.1.
wrapWidth	Wrap long GO term labels to this many characters per line. NULL disables wrapping. Default 40.
fontSize	Base font size passed to theme_cowplot. Default 12.

Value

A ggplot object.

See Also

regulatorGSEA

Examples

se <- getReferenceDatabase(demo = TRUE)
res <- calcEnrichment(
    se,
    upRegulated   = rownames(se)[seq_len(5)],
    downRegulated = rownames(se)[seq(6, 10)],
    background    = rownames(se),
    BPPARAM       = BiocParallel::SerialParam(progressbar = FALSE)
)
go <- regulatorGSEA(res, minSize = 1)
if (any(go$padj <= 0.5)) plotRegulatorGSEA(go, pvalueCutoff = 0.5)

---

Plot signature similarity scores

Description

Draws the top positive and top negative signatures as a lollipop plot.

Usage

plotSimilarity(object, topN = 10, inferActivity = FALSE)

Arguments

object	A SummarizedExperiment containing scores in its colData.
topN	Integer. Number of top positive and top negative hits to show.
inferActivity	Logical. If TRUE, flip the sign of the score for KD/KO perturbations so it reflects predicted gene activity rather than transcriptional similarity. If FALSE (default), the score keeps its similarity meaning and the perturbation direction is shown by point shape: up-pointing triangles for gain-of-function (OE), down-pointing triangles for loss-of-function (KD/KO).

Value

A ggplot object.

Examples

# A PerturbMatch object as returned by calcEnrichment(): one row per
# reference signature, scores in the 'enrichment' assay and the signature
# annotation in rowData().
enr <- cbind(
    score  = c(0.9, -0.7, -0.4, 0.5, -0.2),
    zscore = c(1.6, -1.1, -0.6, 0.8, -0.3)
)
rownames(enr) <- paste0("SIG", seq_len(5))
se <- SummarizedExperiment::SummarizedExperiment(
    assays  = list(enrichment = enr),
    rowData = S4Vectors::DataFrame(
        pert_symbol = c("GENE_A", "GENE_B", "GENE_C", "GENE_D", "GENE_E"),
        pert_type   = c("OE", "KD", "KO", "OE", "KD")
    )
)
pm <- methods::new("PerturbMatch", se)
plotSimilarity(pm, topN = 3)
plotSimilarity(pm, topN = 3, inferActivity = TRUE)

---

Volcano Plot from a SummarizedExperiment Object

Description

Extracts effect size and p-value assays for a given signature from a SummarizedExperiment and delegates to volcanoPlot.

Usage

plotVolcano(
  object,
  signature,
  effectAssay = "logFC",
  pvalueAssay = "padj",
  numberPoints = 5,
  selectedPoints = NULL,
  logFCThreshold = log2(1.5),
  fdrThreshold = 0.05,
  pointSize = 2,
  labSize = 4,
  truncateInfPvalue = TRUE
)

Arguments

object	A SummarizedExperiment object. Must contain assays named by effectAssay and pvalueAssay, and rownames corresponding to gene identifiers.
signature	A single character string matching one of colnames(object), selecting the column to plot.
effectAssay	Character string giving the name of the assay that holds log2 fold-change values. Default: "logFC".
pvalueAssay	Character string giving the name of the assay that holds adjusted p-values (FDR). Default: "padj".
numberPoints	Integer. Number of top up- and down-regulated genes to label automatically (ranked by p-value, then |logFC|). Set to NULL to suppress automatic labelling. Default: 5.
selectedPoints	Character vector of gene IDs to label regardless of significance. Combined with numberPoints labels when both are supplied. Default: NULL.
logFCThreshold	Numeric. Absolute log2 fold-change cut-off used to classify genes as regulated. Default: log2(1.5).
fdrThreshold	Numeric. Adjusted p-value cut-off. Default: 0.05.
pointSize	Numeric. Size of plotted points passed to geom_point. Default: 2.
labSize	Numeric. Text size for gene labels passed to geom_text_repel. Default: 4.
truncateInfPvalue	Logical. If TRUE, genes with an adjusted p-value of exactly zero (an infinite $-\log_{10}$ FDR) are drawn at the largest finite $-\log_{10}$ FDR in the signature, so they appear as a flat row along the top of the plot. If FALSE, these genes are left at their true infinite value and dropped by the plotting layer. Default: TRUE.

Value

A ggplot object.

See Also

volcanoPlot

Examples

set.seed(1)
n <- 30
se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(
        logFC = matrix(rnorm(n * 2), nrow = n,
            dimnames = list(paste0("Gene", seq_len(n)), c("A", "B"))),
        padj = matrix(runif(n * 2), nrow = n,
            dimnames = list(paste0("Gene", seq_len(n)), c("A", "B")))
    )
)

plotVolcano(se, signature = "A")
plotVolcano(se, signature = "A", selectedPoints = c("Gene1", "Gene2"))

---

Prepare a differential expression result for calcEnrichment

Description

Extracts gene identifiers and log2 fold changes from a differential expression result, maps the identifiers to the Entrez gene IDs used by the perturbMatch reference databases, and reports the fraction of the query that was successfully mapped. The returned data frame can be passed directly to calcEnrichment as expressionMatrix.

Usage

prepareQuery(
  object,
  orgDb = "org.Hs.eg.db",
  keytype = NULL,
  reference = NULL,
  minMappedFraction = 0.5,
  dropLowExpressed = TRUE
)

Arguments

object	A differential expression result. One of: a DESeq2 DESeqResults object. a data.frame (e.g. limma topTable() output, or a plain table read from a CSV/TSV) with gene identifiers in its row names or in a column, and a fold-change column named "log2FoldChange", "logFC" or "log2FC". An adjusted p-value column named "padj" or "adj.P.Val", if present, is carried through. Columns are matched by name, not position, and any further columns are ignored. Only adjusted p-values are recognised, so a "pvalue" or "P.Value" column does not drive dropLowExpressed. Gene identifiers are taken from the row names when those are real names, as in DESeqResults and topTable() output. Positional row names, as produced by data.frame() and read.csv(), are ignored, as are all-digit row names when a column names itself as an identifier. The identifiers are then read from a column named like an identifier ("gene", "gene_id", "symbol", "entrezid", "id", ...) if present, otherwise from the first column that is neither the fold-change nor the adjusted p-value column. Which of the two was used is reported in a message.
orgDb	An AnnotationDbi OrgDb object, or the name of an installed OrgDb package, used for ID mapping. Default "org.Hs.eg.db".
keytype	The type of the query gene identifiers (e.g. "ENSEMBL", "SYMBOL"), passed to mapIds. When NULL (default) it is guessed from the identifiers. Identifiers that are already Entrez IDs are used as-is without mapping.
reference	Optional reference database (a SummarizedExperiment). When supplied, the fraction of the mapped query present in rownames(reference) is also reported.
minMappedFraction	Numeric in $[0, 1]$. A warning is emitted if the fraction of query genes mapped to Entrez falls below this value. Default 0.5.
dropLowExpressed	Logical. When the input carries an adjusted p-value column, tools such as DESeq2 set it to NA for genes that were filtered out before testing, typically because they are too lowly expressed. When TRUE (default) these genes are dropped and the number removed is reported. Has no effect when the input has no adjusted p-value column.

Value

A data.frame with a gene column (Entrez IDs) and a log2FC column, plus a padj column when one was available in object. Rows with a missing fold change or an unmapped identifier are dropped, and identifiers mapping to the same Entrez ID are collapsed to the one with the largest absolute fold change.

See Also

calcEnrichment

Examples

# A limma topTable-style data frame (gene IDs in the row names)
tt <- data.frame(
    logFC = c(2.5, -1.8, 0.3),
    adj.P.Val = c(0.001, 0.02, 0.4),
    row.names = c("100", "1000", "10000")
)
# IDs are already Entrez, so no OrgDb is needed
prepareQuery(tt)

---

GO enrichment of perturbMatch regulator hits by GSEA

Description

Runs gene set enrichment analysis (fgsea) over Gene Ontology terms using the reference regulators ranked by their perturbMatch similarity score. This asks which biological processes are enriched among the regulators whose perturbation signatures most closely mimic the query (positive enrichment) or most strongly reverse it (negative enrichment).

Usage

regulatorGSEA(
  object,
  ontology = c("BP", "MF", "CC"),
  aggregate = c("extreme", "mean"),
  inferActivity = TRUE,
  orgDb = "org.Hs.eg.db",
  minSize = 15L,
  maxSize = 500L,
  ...
)

Arguments

object	A PerturbMatch object from calcEnrichment.
ontology	GO ontology to test: one of "BP" (default), "MF" or "CC".
aggregate	How to collapse the multiple signatures of a regulator into one ranking statistic. "extreme" (default) keeps the single largest-magnitude z-score. "mean" averages the activity z-scores.
inferActivity	Logical. When TRUE (default) the z-score sign is flipped for KD/KO perturbations before aggregating, so the ranking reflects inferred gene activity rather than raw transcriptional similarity.
orgDb	An AnnotationDbi OrgDb object, or the name of an installed OrgDb package, supplying the GO annotation and the symbol/Entrez mapping. Default "org.Hs.eg.db", matching the species of the packaged reference databases.
minSize, maxSize	Minimum and maximum GO gene-set sizes to test, passed to fgsea.
...	Further arguments passed to fgsea.

Details

A reference database usually holds several signatures per regulator gene (different cell lines, and knockdown versus overexpression). These are collapsed to a single ranking statistic per gene before scoring:

1. scores are put on a consistent inferred-activity scale by flipping the sign of the z-score for KD/KO perturbations (as in plotSimilarity(inferActivity = TRUE)), controlled by inferActivity.

2. the per-signature scores of each gene are aggregated with aggregate (the single "extreme" value with the largest magnitude, or the "mean" of the activity z-scores).

The resulting one-value-per-gene ranking is scored against GO gene sets, so the positive end of the ranking (mimicking regulators) and the negative end (reversing regulators) are tested together and separated by the sign of each term's normalised enrichment score (NES).

Gene sets are built from orgDb (using the ancestor-propagated GOALL mapping restricted to ontology), which also supplies the regulator symbol-to-Entrez mapping. Human-readable term names are added from GO.db when it is installed. Otherwise term falls back to the GO accession.

Value

A data.frame, one row per GO term ordered by adjusted p-value, with columns ID (GO accession), term, NES (positive = enriched among mimicking regulators, negative = among reversing regulators), pval, padj, size and leadingEdge (comma-separated regulator gene symbols).

See Also

plotRegulatorGSEA, calcEnrichment

Examples

se <- getReferenceDatabase(demo = TRUE)
res <- calcEnrichment(
    se,
    upRegulated   = rownames(se)[seq_len(5)],
    downRegulated = rownames(se)[seq(6, 10)],
    background    = rownames(se),
    BPPARAM       = BiocParallel::SerialParam(progressbar = FALSE)
)
# minSize is lowered because the demo database is tiny. Keep the default
# on a full reference.
go <- regulatorGSEA(res, minSize = 1)
head(go)

---

Tidy table of regulator scores from a PerturbMatch object

Description

Combines rowData with the enrichment scores into a single data.frame, one row per reference signature.

Usage

regulatorTable(object)

Arguments

object

A PerturbMatch object.

Details

The name deliberately avoids result/results, which would sit one letter from results in a script that uses both.

Value

A data.frame with one row per feature holding every rowData column plus a score column.

Examples

mat <- matrix(runif(3), nrow = 3, dimnames = list(
    c("GENE1", "GENE2", "GENE3"), "score"
))
se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(enrichment = mat)
)
S4Vectors::metadata(se) <- list(
    up_regulated   = "GENE1",
    down_regulated = "GENE3",
    background     = c("GENE1", "GENE2", "GENE3")
)
pm <- new("PerturbMatch", se)
regulatorTable(pm)

---

Show a PerturbMatch object

Description

Prints the inherited SummarizedExperiment summary followed by the gene sets and parameters recorded when the object was built.

Usage

## S4 method for signature 'PerturbMatch'
show(object)

Arguments

object

A PerturbMatch object.

Value

Invisibly returns NULL. Called for its printing side effect.

Examples

pm <- methods::new(
    "PerturbMatch",
    SummarizedExperiment::SummarizedExperiment(
        assays = list(
            enrichment = matrix(0.5, dimnames = list("SIG1", "score"))
        ),
        metadata = list(
            up_regulated   = c("GENE1", "GENE2"),
            down_regulated = "GENE3",
            background     = c("GENE1", "GENE2", "GENE3"),
            threshold      = 0.5,
            topN           = 250L,
            topSig         = 100L
        )
    )
)
show(pm)

---

Top-scoring regulator signatures

Description

Returns the n highest-scoring rows of regulatorTable, ranked by absolute enrichment score.

Usage

topRegulators(object, n = 5, ...)

Arguments

object	A PerturbMatch object.
n	Number of signatures to return, capped at nrow(object). Default 5.
...	Additional arguments passed to methods.

Value

A data.frame with at most n rows, ordered by descending absolute score.

Examples

mat <- matrix(runif(5), nrow = 5, dimnames = list(
    paste0("GENE", seq_len(5)), "score"
))
se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(enrichment = mat)
)
S4Vectors::metadata(se) <- list(
    up_regulated   = "GENE1",
    down_regulated = "GENE5",
    background     = paste0("GENE", seq_len(5))
)
pm <- new("PerturbMatch", se)
topRegulators(pm, n = 3)

---

Up-regulated genes from a PerturbMatch object

Description

Returns the up-regulated gene set recorded in metadata() when the object was built by calcEnrichment.

Usage

upRegulated(object)

Arguments

object

A PerturbMatch object.

Value

A character vector of gene identifiers, or NULL if unset.

Examples

pm <- methods::new(
    "PerturbMatch",
    SummarizedExperiment::SummarizedExperiment(
        assays = list(
            enrichment = matrix(0.5, dimnames = list("SIG1", "score"))
        ),
        metadata = list(up_regulated = c("GENE1", "GENE2"),
            background = c("GENE1", "GENE2"))))
upRegulated(pm)

---

Volcano Plot of Differential Expression Results

Description

Renders a volcano plot of log2 fold-change against $-\log_{10}$(FDR), colouring genes as up-regulated, down-regulated, or unchanged relative to supplied thresholds. The most significant genes at each tail are labelled automatically, and arbitrary genes can be highlighted via selectedPoints.

Usage

volcanoPlot(
  data,
  numberPoints = 5,
  selectedPoints = NULL,
  logFCThreshold = log2(1.5),
  fdrThreshold = 0.05,
  pointSize = 2,
  labSize = 4,
  truncateInfPvalue = TRUE
)

Arguments

data	A data frame with at least three columns, in this order: Gene / feature identifier (character). Log2 fold-change (numeric). Adjusted p-value / FDR (numeric, in the interval $[0, 1]$). Additional columns are ignored.
numberPoints	Integer. Number of top up- and down-regulated genes to label (ranked first by p-value, then by |logFC|). Pass NULL to disable automatic labelling. Default: 5.
selectedPoints	Character vector of gene identifiers to label regardless of their statistics. Merged with automatic labels when numberPoints is also set. Default: NULL.
logFCThreshold	Numeric. Minimum absolute log2 fold-change required for a gene to be called regulated. Default: log2(1.5) ($\approx 0.585$).
fdrThreshold	Numeric. Maximum adjusted p-value for a gene to be called regulated. Shown as a horizontal dashed line. Default: 0.05.
pointSize	Numeric. Passed to geom_point as size. Default: 2.
labSize	Numeric. Passed to geom_text_repel as size. Default: 4.
truncateInfPvalue	Logical. If TRUE, genes with $p = 0$ (an infinite $-\log_{10}$ FDR) are placed at the largest finite $-\log_{10}$ FDR, forming a flat row along the top of the plot. If FALSE, they keep their infinite value and are dropped by the plotting layer. Default: TRUE.

Details

Infinite $-\log_{10}$ p-values (arising from $p = 0$) are, by default, truncated to the largest finite value to keep the plot axis bounded. See truncateInfPvalue to switch this off.

Value

A ggplot object. Colour scale:

firebrick3

Up-regulated (logFC >= logFCThreshold and padj <= fdrThreshold).

dodgerblue3

Down-regulated (logFC <= -logFCThreshold and padj <= fdrThreshold).

black

Unchanged (all other genes).

Examples

set.seed(42)
n <- 500
df <- data.frame(
    gene  = paste0("Gene", seq_len(n)),
    logFC = rnorm(n, sd = 2),
    padj  = rbeta(n, 0.5, 5)
)
volcanoPlot(df)

# Highlight specific genes
volcanoPlot(df, selectedPoints = c("Gene1", "Gene42"))

# Stricter thresholds, no automatic labels
volcanoPlot(df, numberPoints = NULL,
    logFCThreshold = log2(2), fdrThreshold = 0.01)

---