Skip to contents

Chemical Enrichment Analysis Using CTD

Source: R/enrichment_analysis.R
enrichment_CTD.Rd

Identifies chemicals whose known gene targets are significantly enriched in a user-supplied gene list or expression matrix, using data from the Comparative Toxicogenomics Database (CTD).

Four methods are available:

ORA

Over-Representation Analysis (default). Tests whether the overlap between your gene list and each chemical's target genes is larger than expected by chance. Uses a hypergeometric test computed directly on phyper. Input: data frame with column EntrezID (character or numeric Entrez gene IDs) and an optional numeric value column.

GSEA

Gene Set Enrichment Analysis. Uses a ranked gene list (ranked by the numeric column in the input, e.g. p-values or fold changes) to detect chemicals whose targets cluster toward the top or bottom of the ranking. Uses fgsea.

CAMERA

Competitive gene-set test accounting for inter-gene correlation. Uses camera. Input: a numeric expression matrix (genes x samples) or a SummarizedExperiment, plus a design matrix and a contrast.

GSVA

Gene Set Variation Analysis (sample-level scoring). Returns per-sample enrichment scores for each chemical, suitable for downstream clustering or association testing. Uses gsva. Input: a numeric expression matrix (genes x samples) or a SummarizedExperiment, which is also the container the scores are returned in.

Usage

enrichment_CTD(
  x,
  method = c("ORA", "GSEA", "CAMERA", "GSVA"),
  design = NULL,
  contrast = NULL,
  id_type = NULL,
  pAdjustMethod = "BH",
  interaction_types = NULL,
  gene_id_type = c("symbol", "entrez"),
  universe = NULL,
  alpha = NULL,
  alpha_column = NULL,
  assay = NULL,
  ...
)

Arguments

x

The input. Its expected type depends on method:

  • For "ORA" and "GSEA": a data frame with at least two columns, EntrezID (character or numeric Entrez gene IDs) and a numeric value column (e.g. p-value). For ORA this is either the genes you have already selected, in which case say what the background was with universe, or the whole result table with alpha to select from it. For GSEA, an optional column named stat can be added with a signed ranking statistic (e.g. the moderated t-statistic from limma::eBayes()); when present it is used directly for ranking, preserving directionality and avoiding ties. When absent, the second column is transformed via -log10() with a warning. The second column is ignored by ORA.

  • For "CAMERA" and "GSVA": either a numeric expression matrix with genes in rows and samples in columns, or a SummarizedExperiment whose assay holds that matrix (see assay to choose which one). rownames(x) must be either Entrez IDs or HGNC SYMBOLs.

method

Character. Enrichment method: "ORA" (default), "GSEA", "CAMERA", or "GSVA".

design

Design matrix (required when method = "CAMERA").

contrast

Contrast specification for camera (column number, column name, or numeric vector). Required when method = "CAMERA".

id_type

Either "entrez", "symbol", or NULL (default) for auto-detection from rownames(x). Only used when method is "CAMERA" or "GSVA".

pAdjustMethod

Character. Multiple-testing correction applied to the raw p-values, passed to p.adjust. Any value in stats::p.adjust.methods is accepted: "holm", "hochberg", "hommel", "bonferroni", "BH" (Benjamini-Hochberg, the default), "BY", "fdr" (alias for "BH"), or "none". See p.adjust for the meaning of each method. Not used for method = "GSVA" (which returns per-sample scores, not p-values).

interaction_types

Character vector of CTD InteractionActions values to retain when building gene sets, or NULL (default) to use all cached interactions. Values follow the verb\^{}noun convention used by CTD, e.g. "increases\^{}expression", "decreases\^{}expression", "affects\^{}binding". A gene is included in a chemical's set if any of its recorded interaction actions matches one of the specified types. Requires that import_CTD() has been run (the filter is applied to the cached ctd_interactions.rda file). Restricting to expression interactions is recommended for RNA-seq analyses to improve biological specificity.

gene_id_type

Character. Identifier type used in the EnrichedGenes output column: "symbol" (default) returns HGNC gene symbols with Entrez ID as fallback for unmapped genes; "entrez" skips the symbol lookup and returns Entrez IDs directly. Only applies to "ORA" and "GSEA"; ignored by "CAMERA" and "GSVA".

universe

Background gene identifiers for "ORA": every gene that entered your differential test, not only the significant ones and not every gene sequenced. NULL (default) falls back to every gene in the CTD gene sets, and says so, because the package cannot know what your platform measured.

This is the input that decides whether an ORA result means anything. A background wider than what the experiment could detect fills the urn with genes that could never have been drawn, so the overlap looks more selective than it was and p-values come out too small. The error is anti-conservative. On the RNA-seq example bundled with this package, the fallback background returns 32 chemicals at FDR < 0.05 and the correct one returns 19.

ORA is the only method that needs this argument, and the reason is the shape of its input. A bare gene list carries no record of what was measurable, so the background has to be supplied separately. GSEA ranks the whole list you give it, which is already the background; "CAMERA" and "GSVA" intersect the gene sets with rownames(x), so theirs is the set of measured genes by construction. Passing universe to any of those three raises a warning rather than being quietly dropped.

alpha

Significance threshold for "ORA", or NULL (default). Supplying it says that x is the whole result table rather than a list already filtered: the rows whose second column falls below alpha become the genes to test, and every row becomes the background.

Every row of the table becomes the background, the selected genes included: a hypergeometric test draws \(n\) genes from an urn of \(N\), and the drawn ones were in the urn. The background is not the non-significant remainder.

This is the safer way to run ORA, because the background is then derived from the same object as the gene list and cannot disagree with it. Filtering first and passing universe separately asks the caller to reconnect two things that were together a moment earlier, and that reconnection is what goes wrong.

Name the column to threshold with alpha_column: which p-value to judge on is the researcher's decision, not the package's. The function reports the column it used, how many genes passed and how many form the background.

Mutually exclusive with universe, and ignored by the other three methods, which already hold their own background.

alpha_column

Which column alpha applies to: a name, a positive index, or NULL (default) for the second column. Naming it matters on a real result table, where the second column is usually a fold change: limma::topTable() puts logFC there. Pass the whole table and say which p-value to judge on, alpha_column = "padj" or "adj.P.Val", rather than cutting the table down to two columns first.

assay

Which assay to use when x is a SummarizedExperiment: an assay name, a positive index, or NULL (default) for the first assay. Ignored when x is a matrix or a data frame.

...

Additional arguments forwarded to the underlying engine: ora for ORA (universe, minGSSize, maxGSSize; both thresholds are chosen for CTD rather than inherited, minGSSize defaulting to 2 and maxGSSize to Inf, see ora for the measurements behind them), fgseaMultilevel for GSEA (e.g. minSize, maxSize, nproc), camera for CAMERA, gsva for GSVA (e.g. minSize, maxSize).

Value

  • For "ORA", "GSEA", and "CAMERA": a data frame of enrichment results sorted by PValueAdjusted ascending. All three methods share the leading columns ChemicalID, ChemicalName, Method, PValue, PValueAdjusted; method-specific extras follow (see the package vignette for the full per-method schema).

  • For "GSVA": enrichment scores with chemicals (CTD chemical IDs) in rows and samples in columns. The container follows the input: a matrix in returns a numeric matrix, while a SummarizedExperiment in returns a SummarizedExperiment whose assay holds the scores and whose colData is carried over from the input, so sample annotation stays attached to the results.

Details

Before calling this function you must import the CTD data once with import_CTD. If the cached data is not found, the function stops with an informative error message.

Data Licensing Disclaimer

This package does not bundle or redistribute any CTD data. The Comparative Toxicogenomics Database is maintained by NC State University and its data are subject to specific licensing terms. Users must download the data directly from https://ctdbase.org and comply with the CTD Terms of Service (https://ctdbase.org/about/legal.jsp).

See also

import_CTD to import and cache the CTD data; plot_CTD to visualize results.

Examples

# Import the bundled sample data first:
# Examples write to a temporary cache, so running them cannot
# disturb CTD data you have already imported. Set the same option
# yourself to keep an analysis isolated from your main cache.
options(ctdR.cache = tempfile())

sample_file <- system.file(
    "extdata", "CTD_chem_gene_ixns_sample.csv",
    package = "ctdR"
)
import_CTD(sample_file)
#> Reading CTD chemical-gene interactions from: /home/runner/work/_temp/Library/ctdR/extdata/CTD_chem_gene_ixns_sample.csv
#> Filtered to 86 human interactions
#> Mapping genes for 10 chemicals...
#> Warning: 10 ChemicalID(s) appear with more than one ChemicalName in the CTD file; only the first name per ID is retained. Affected IDs: D000082, D001564, D002104, D003907, D004958 ... (and 5 more)
#> CTD data cached successfully in: /tmp/Rtmpt9pODL/file1b1316599d6c
#>   10 chemicals | 17 unique genes | 0 s
#>   CTD release: Mon Jan 01 00:00:00 EST 2024

# ORA / GSEA: prepare a gene list with Entrez IDs and a numeric value
genes <- data.frame(
    EntrezID = c("7124", "3569", "7157", "672", "1956"),
    pvalue = c(0.001, 0.003, 0.01, 0.02, 0.05)
)
ora_results <- enrichment_CTD(genes, method = "ORA")
#> background: all 17 genes in the CTD sets, because no 'universe' was given. If your experiment could only detect some of them, pass those as 'universe': a background wider than what was measurable makes p-values too small.
gsea_results <- enrichment_CTD(genes, method = "GSEA")
#> Warning: GSEA: no 'stat' column found in input. Falling back to -log10(second column) for ranking, which loses directionality and may produce ties at non-significant genes. Add a signed ranking statistic (e.g. the moderated t-statistic from limma::eBayes()) as a column named 'stat' to suppress this warning.
#> Warning: All values in the stats vector are greater than zero and scoreType is "std", maybe you should switch to scoreType = "pos".

# CAMERA / GSVA: expression data plus, for CAMERA, design + contrast.
# Uses the bundled GSE311566 subset, a SummarizedExperiment
# (Dex vs DMSO, female PBMCs; see inst/extdata/README.md).
se <- readRDS(system.file(
    "extdata", "GSE311566_subset.rds", package = "ctdR"
))
d <- model.matrix(~ se$group)
camera_results <- enrichment_CTD(se, method = "CAMERA",
    design = d, contrast = 2)

# GSVA: SummarizedExperiment in, SummarizedExperiment out,
# so the sample annotation stays attached to the scores.
gsva_scores <- enrichment_CTD(se, method = "GSVA")
#> ℹ No assay name provided; using default assay 'logcounts'
#> ℹ GSVA version 2.6.6
#> ℹ Searching for rows with constant values
#> ℹ Calculating GSVA ranks
#> ℹ kcdf='auto' (default)
#> ℹ GSVA dense (classical) algorithm
#> ℹ Row-wise ECDF estimation with Gaussian kernels
#> ℹ Calculating row ECDFs
#> ℹ Calculating column ranks
#> ℹ GSVA dense (classical) algorithm
#> ℹ Calculating GSVA scores for 10 gene sets
#> ✔ Calculations finished
SummarizedExperiment::colData(gsva_scores)$group
#> [1] DMSO DMSO DMSO DMSO Dex  Dex  Dex 
#> Levels: DMSO Dex