---
title: "Introduction to ctdR"
author:
  - name: Luigi Corsaro
    affiliation: ""
    orcid: 0000-0003-1218-230X
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Introduction to ctdR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
# Build this vignette against a throwaway cache. Without it, knitting
# would overwrite whatever CTD data the reader has already imported
# with the ten-chemical sample used for the examples below.
options(ctdR.cache = tempfile())
```

# Overview

**ctdR** identifies chemicals significantly associated
with a set of genes using data from the
[Comparative Toxicogenomics Database (CTD)](https://ctdbase.org).

Four enrichment methods are supported through a single
`enrichment_CTD()` interface; the input shape depends on
the method:

- **ORA** (Over-Representation Analysis) -- gene list +
    hypergeometric test, computed directly on
    `stats::phyper`.
- **GSEA** (Gene Set Enrichment Analysis) -- ranked gene
    list + permutation test. Powered by `fgsea::fgsea`.
- **CAMERA** (competitive gene-set test) -- expression
    matrix or `SummarizedExperiment` + design + contrast.
    Corrects for inter-gene correlation within each
    chemical's gene set. Powered by `limma::camera`.
- **GSVA** (Gene Set Variation Analysis) -- expression
    matrix or `SummarizedExperiment` + per-sample scoring.
    Returns chemical x sample scores in the same container
    it was given. Powered by `GSVA::gsva`.

# Data licensing disclaimer

> **This package does NOT bundle, redistribute, or embed
> any data from the Comparative Toxicogenomics Database.**
> CTD data are created and maintained by NC State
> University and 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).

# Installation

ctdR has been accepted into Bioconductor. For now it is available only in
Bioconductor devel (3.24); after the next Bioconductor release it will be
installable from the standard (release) repository.

## From Bioconductor devel (current installation method)

    if (!require("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install(version = "devel")
    BiocManager::install("ctdR")

## From Bioconductor release (after the next Bioconductor release)

    BiocManager::install("ctdR")

## From GitHub (latest development version)

    # install.packages("devtools")
    devtools::install_github("drake69/ctdR")

# Quick start

The ctdR workflow has three steps:

1. **Download** the CTD data file (once, manually)
2. **Import** the data into ctdR (once)
3. **Analyse** your gene list (as many times as needed)

## Step 1 -- Download CTD data

Download **`CTD_chem_gene_ixns.csv.gz`** from
<https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz>.

Then decompress it:

```bash
gunzip CTD_chem_gene_ixns.csv.gz
```

This produces `CTD_chem_gene_ixns.csv` (several GB
uncompressed).

## Step 2 -- Import into ctdR

In production you would run:

```{r import-real, eval=FALSE}
library(ctdR)
import_CTD("~/Downloads/CTD_chem_gene_ixns.csv")
```

For this vignette we use a small **synthetic** dataset
bundled with the package:

```{r import-sample}
library(ctdR)
sample_file <- system.file(
    "extdata", "CTD_chem_gene_ixns_sample.csv",
    package = "ctdR"
)
import_CTD(sample_file)
```

`import_CTD()` performs the following:

1. Reads the CSV (skipping the 27 CTD header lines).
2. Filters interactions to **Homo sapiens** only
    (OrganismID 9606).
3. Collects Entrez gene IDs for each chemical.
4. Maps Entrez IDs to HGNC gene symbols via
    `org.Hs.eg.db`.
5. Caches the processed data locally via `BiocFileCache`
    (under `tools::R_user_dir("ctdR", "cache")`).

This step takes several minutes with the full CTD file.
You only need to run it once -- or again when you
download a newer CTD release.

### Which CTD release did I use?

CTD is re-released continuously and its downloads are
not versioned in the filename: two copies of
`CTD_chem_gene_ixns.csv` months apart are
indistinguishable from the outside. The release date
lives inside the file, on the `# Report created:` line
of its header, and it is the only thing that identifies
which snapshot a result came from.

`import_CTD()` reads that line, reports it, and keeps
it. Every result then carries it, so an analysis can
state the version of the data behind it without anyone
having to remember:

```{r provenance}
# With no argument: what is in the cache right now.
ctd_provenance()
```

Every result carries the same record, so you can also
ask an object you are about to report on:
`ctd_provenance(ora_results)`.

The record holds the CTD release date, where the file
came from, when you imported it, how much of it was
kept, and the ctdR version. Release date and import
date answer different questions and are both recorded:
the first says which data, the second says when you
took it.

Use `ctd_provenance()` rather than reaching into the
object, because the record is stored wherever the
container provides room for it. Results that carry a
`metadata()` slot keep it there, which is the case for
the `SummarizedExperiment` returned by GSVA and for a
`CTDFile` read with `import()`; the data frames from
ORA, GSEA and CAMERA keep it in an attribute. The
accessor hides the difference.

The record follows subsetting, ordering, `head()` and
the usual **dplyr** verbs. It does not follow `merge()`
or `subset()`, which drop attributes; `ctd_provenance()`
warns rather than quietly reporting nothing, so take
the record before those if you need it after.

## Step 3 -- Run enrichment analysis

### Prepare your gene list

Your input must be a data frame with at least two
columns:

| Column | Description |
|---|---|
| `EntrezID` | Character or numeric Entrez gene IDs |
| *(2nd column)* | Numeric value per gene (e.g. p-value) |

The second column is used for ranking in GSEA and is
ignored in ORA.

```{r prepare-genes}
genes <- data.frame(
    EntrezID = c(
        "7124", "3569", "7157", "672", "1956",
        "4609", "3845", "207", "5290", "3553"
    ),
    pvalue = c(
        0.001, 0.003, 0.005, 0.008, 0.01,
        0.02, 0.03, 0.04, 0.05, 0.06
    )
)
```

### Shared output schema

All three data-frame-returning methods (ORA, GSEA,
CAMERA) share the **same leading columns**. This makes
it trivial to combine results across methods (e.g.
`dplyr::bind_rows()` on a list of result frames).

| Column | Type | Description |
|---|---|---|
| `ChemicalID` | character | CTD chemical identifier |
| `ChemicalName` | character | Human-readable chemical name |
| `Method` | character | One of `"ORA"`, `"GSEA"`, `"CAMERA"` |
| `PValue` | numeric | Raw p-value from the underlying method |
| `PValueAdjusted` | numeric | Multiple-testing-corrected p-value (`pAdjustMethod`) |

Method-specific extras follow these front columns and
differ by method (documented in the sub-sections
below). Rows are sorted by `PValueAdjusted` ascending.

### Over-Representation Analysis (ORA)

ORA tests whether the overlap between your gene list
and each chemical's known gene targets is significantly
larger than expected by chance.

```{r ora}
ora_results <- enrichment_CTD(genes, method = "ORA")
head(ora_results)
```

**ORA method-specific columns** (in addition to the
shared schema above):

| Column | Description |
|---|---|
| `GeneRatio` | Proportion of input genes in the set |
| `BackgroundRatio` | Background ratio |
| `FoldEnrichment` | GeneRatio / BackgroundRatio |
| `EnrichedGenes` | Enriched gene symbols (comma-separated) |
| `Count` | Number of overlapping genes |

### Gene Set Enrichment Analysis (GSEA)

GSEA uses the full ranked gene list to detect chemicals
whose targets cluster toward the top or bottom of the
ranking. For best results, supply a signed ranking
statistic (e.g. the moderated t-statistic from
`limma::eBayes()`) as a column named `stat`. When `stat`
is absent, ctdR falls back to `-log10(pvalue)`, which
loses directionality and may produce ties.

```{r gsea}
# Quick-start: no signed statistic available in this
# toy gene list, so the pvalue fallback is used.
gsea_results <- enrichment_CTD(genes, method = "GSEA")
head(gsea_results)
```

**GSEA method-specific columns** (in addition to the
shared schema above):

| Column | Description |
|---|---|
| `EnrichmentScore` | Enrichment score (`ES` from fgsea) |
| `NormalizedEnrichmentScore` | Normalized enrichment score (`NES` from fgsea): the enrichment score scaled for gene set size. This is the quantity to compare between chemicals |
| `GeneSetSize` | Size of the gene set used |
| `LeadingEdge` | Leading-edge gene subset |
| `EnrichedGenes` | Comma-separated enriched genes |

### End-to-end example with real RNA-seq data (GSE311566)

The hand-coded `genes` data frame above is convenient for the
Quick start but does not exercise the matrix-based methods. This
section ties everything together on a small subset of
[GEO GSE311566](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE311566)
-- human PBMCs treated with **dexamethasone** vs **vehicle (DMSO)**
in female donors. The bundled subset
`inst/extdata/GSE311566_subset.rds` contains 1,500 top-variance
genes (plus the 17 genes referenced by the toy CTD sample),
log2-transformed normalised counts, 7 samples total
(4 DMSO + 3 Dex). See `inst/extdata/README.md` for provenance.

It is stored as a
[`SummarizedExperiment`](https://bioconductor.org/packages/SummarizedExperiment),
the standard Bioconductor container for assay data: the expression
matrix and the sample annotation live in one object, so subsetting or
reordering samples moves both together and cannot silently
desynchronise the group labels from the columns they describe.

```{r load-real}
library(SummarizedExperiment)

se <- readRDS(system.file(
    "extdata", "GSE311566_subset.rds", package = "ctdR"
))
se
table(se$group)
```

The assay is reached with `assay(se)` and the sample annotation with
`colData(se)`, or with the `se$column` shorthand used below.

We compute a deliberately minimal per-gene differential
expression with base R only: a two-sample `t.test` per gene with
BH-adjusted p-values. Production analyses should use `limma`,
`DESeq2`, or `edgeR` for proper count-based modelling; this
ascetic version keeps the example self-contained.

```{r de-real}
expr <- assay(se)
grp  <- se$group

de <- t(apply(expr, 1, function(y) {
    tt <- stats::t.test(y ~ grp)
    c(log2FC = unname(diff(tt$estimate)),
      pvalue = tt$p.value)
}))
de <- as.data.frame(de)
de$padj <- stats::p.adjust(de$pvalue, method = "BH")
de$EntrezID <- rownames(de)
de <- de[, c("EntrezID", "log2FC", "pvalue", "padj")]
head(de[order(de$padj), ])
c(padj_lt_05 = sum(de$padj < 0.05),
  pvalue_lt_05 = sum(de$pvalue < 0.05))
```

#### Real-data GSEA

GSEA uses the full ranked gene list and tends to surface the
expected chemical (Dexamethasone) near the top even with this
small DE. We supply `log2FC` as the `stat` column so that
direction (up/down) is preserved in the ranking.

```{r gsea-real}
gsea_real <- enrichment_CTD(
    data.frame(EntrezID = de$EntrezID,
               pvalue   = de$pvalue,
               stat     = de$log2FC),
    method = "GSEA"
)
head(gsea_real[, c("ChemicalID", "ChemicalName",
                   "PValue", "PValueAdjusted",
                   "NormalizedEnrichmentScore")])
```

#### Real-data ORA

ORA is statistically meaningful only with the **full CTD download**
(~16,000 chemicals × millions of gene interactions, available at
<https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz>). The bundled
`CTD_chem_gene_ixns_sample.csv` covers 10 chemicals × 17 genes — its
sole purpose is API demonstration. With only 7 samples and a universe
this small, ORA is expected to return zero enriched chemicals. Set
`is_subset <- FALSE` in the chunk below to run a fully powered analysis
on the real CTD. All other methods (GSEA, CAMERA, GSVA) work correctly
with the bundled data.

```{r ora-is-subset}
## Set is_subset <- FALSE to run ORA on the full CTD download.
## Requires CTD_chem_gene_ixns.csv.gz from ctdbase.org, decompressed locally.
is_subset <- TRUE
```

```{r ora-real}
sig <- de[de$pvalue < 0.05, c("EntrezID", "pvalue")]

if (is_subset) {
    ora_real <- enrichment_CTD(sig, method = "ORA")
    if (nrow(ora_real)) {
        head(ora_real[, c("ChemicalID", "ChemicalName",
                          "PValue", "PValueAdjusted", "Count")])
    } else {
        message(
            "ORA returned 0 results on the toy CTD sample (expected).\n",
            "Set is_subset <- FALSE and provide CTD_chem_gene_ixns.csv",
            " for a powered analysis."
        )
    }
} else {
    ## Full CTD path -----------------------------------------------------------
    ## 1. Download https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz
    ## 2. gunzip CTD_chem_gene_ixns.csv.gz
    ## 3. Set ctd_path to the decompressed file location
    ctd_path <- "~/Downloads/CTD_chem_gene_ixns.csv"  # adjust to your path
    import_CTD(ctd_path)
    ora_full <- enrichment_CTD(sig, method = "ORA")
    if (nrow(ora_full)) {
        head(ora_full[order(ora_full$PValue),
                      c("ChemicalID", "ChemicalName",
                        "PValue", "PValueAdjusted",
                        "FoldEnrichment", "Count")])
    } else {
        message("No chemicals enriched — check CTD file path and gene list.")
    }
}
```

### CAMERA (competitive test with inter-gene correlation)

CAMERA tests, for each chemical, whether its target
genes show a stronger differential signal than the rest
of the transcriptome, while correcting for inter-gene
correlation within the gene set.

Unlike ORA or GSEA, CAMERA does **not** take a pre-computed
list of differentially expressed genes. Its inputs are the
**full expression data** (all genes by all samples), supplied
either as a matrix or as a `SummarizedExperiment`, together
with a **`design`** matrix and a **`contrast`**:
CAMERA fits the differential model and computes the test
*internally*. The chemical-to-gene mapping (the gene-set
"index", built by `ctdR` from CTD targets intersected with
`rownames(x)`) is handled for you. Note that `design` and
`contrast` are two distinct inputs: `design`
(e.g. `model.matrix(~ se$group)`) describes the *whole*
experimental layout, intercept plus any covariates, while
`contrast` selects *which* coefficient to test (here column 2,
Dex vs DMSO). This separation lets you adjust for batch or
other nuisance covariates while keeping the test focused on a
single comparison.

Passing the `SummarizedExperiment` directly:

```{r camera}
design <- model.matrix(~ se$group)

camera_results <- enrichment_CTD(
    se,
    method   = "CAMERA",
    design   = design,
    contrast = 2  # column 2 of design = Dex vs DMSO
)
head(camera_results)
```

**CAMERA method-specific columns** (in addition to the
shared schema above):

| Column | Description |
|---|---|
| `GeneSetSize` | Gene set size after intersection with `rownames(x)` |
| `Direction` | `"Up"` or `"Down"` |
| `Correlation` | Inter-gene correlation (present only when estimated by `limma::camera`; absent when `inter.gene.cor` is fixed) |

### GSVA (per-sample scoring)

GSVA produces a **per-sample** enrichment score for each
chemical, with chemicals in rows and samples in columns,
suitable for clustering, association tests against
phenotypes, survival analysis, or heatmap visualization.

The result follows the input: give it a matrix and you get a
matrix back, give it a `SummarizedExperiment` and you get one
back, with `colData` carried over. That second form is the
useful one here, because the sample annotation needed to
interpret the scores stays attached to them.

```{r gsva}
gsva_scores <- enrichment_CTD(se, method = "GSVA")
dim(gsva_scores)  # chemicals x samples
head(assay(gsva_scores), 3)

# the treatment groups are still there, next to the scores
colData(gsva_scores)
```

Tune the underlying `GSVA::gsvaParam()` through `...`,
for example to restrict to gene sets of a given size:

```{r gsva-tune}
gsva_strict <- enrichment_CTD(
    se, method = "GSVA",
    minSize = 5, maxSize = 500
)
nrow(gsva_strict)
```

If the object carries several assays, pick one with
`assay = "logcounts"` or by index; the default is the first.

## Visualizing results

`plot_CTD()` creates publication-ready plots from enrichment
results. It auto-detects the method: bar/dot plots of fold
enrichment for ORA/GSEA, bar/dot plots of `-log10(padj)`
coloured by direction for CAMERA, and a sample-level
heatmap of the top-variance chemicals for GSVA.

```{r plot-bar, fig.width=7, fig.height=5}
# ORA bar plot of top enriched chemicals
plot_CTD(ora_results, type = "bar")
```

```{r plot-dot, fig.width=7, fig.height=5}
# ORA dot plot (size = gene count, color = adjusted p-value)
plot_CTD(ora_results, type = "dot", n = 10)
```

```{r plot-camera, fig.width=7, fig.height=5}
# CAMERA bar plot: x-axis = -log10(padj), fill = Direction
plot_CTD(camera_results, type = "bar")
```

```{r plot-gsva, fig.width=7, fig.height=5}
# GSVA heatmap: top-variance chemicals across samples
plot_CTD(gsva_scores)
```

## Adjusting for multiple testing

By default, `enrichment_CTD()` uses the Benjamini-Hochberg
(`"BH"`) method for p-value adjustment. You can change this
via the `pAdjustMethod` parameter:

```{r padj-methods}
# Bonferroni correction (more conservative)
ora_bonf <- enrichment_CTD(genes, method = "ORA",
    pAdjustMethod = "bonferroni")

# No adjustment
ora_raw <- enrichment_CTD(genes, method = "ORA",
    pAdjustMethod = "none")

# Compare the adjusted p-value of the top hit
data.frame(
    method = c("BH (default)", "bonferroni", "none"),
    top_padj = c(min(ora_results$padj),
        min(ora_bonf$padj), min(ora_raw$padj))
)
```

Available methods: `"BH"` (default), `"bonferroni"`,
`"fdr"` (alias for BH), `"none"`.

## Gene set size filters and background universe

Each method applies a gene set size filter and defines a background
universe, but the degree of user control differs:

| Method | Size filter | Configurable | Background universe | Configurable |
|--------|-------------|:---:|---------------------|:---:|
| ORA | `minGSSize` (default 2) / `maxGSSize` (default `Inf`) | ✓ | All genes in TERM2GENE (default) or user-supplied `universe` | ✓ |
| GSEA | `minSize` / `maxSize` | ✓ | No universe concept — uses the full ranked list | — |
| CAMERA | ≥ 2 genes (hard-coded) | ✗ | `rownames(x)` | implicit |
| GSVA | `minSize` / `maxSize` | ✓ | `rownames(x)` | implicit |

ORA, GSEA, and GSVA accept size-filter arguments via `...`, which are
forwarded to the underlying engine. CAMERA's minimum of 2 genes is
hard-coded and cannot be changed.

### The background universe is the choice that matters most

Of everything on this page, this is the argument worth getting right.

A hypergeometric test asks: drawing *n* genes at random from a
background of *N*, how many would land in this chemical's set? If *N*
contains genes your experiment could never have detected, those genes
fill the urn with balls that cannot be drawn. The overlap you observed
then looks more selective than it was, and the p-value comes out too
small. The error is **anti-conservative**: it manufactures significance
rather than hiding it.

`ctdR` defaults to every gene in the CTD gene sets, and says so when it
does. That is a fallback, not a recommendation: the package cannot know
what your platform measured. On the RNA-seq analysis in this package's
tutorial the difference is not subtle:

| Background | *N* | Chemicals significant at FDR < 0.05 |
|---|---:|---:|
| All CTD genes (the default) | 27,444 | **32** |
| Genes that entered the DE test | 24,845 | **19** |

Thirteen of those thirty-two are produced by the background alone.
Rankings move too: benzo(a)pyrene sits fourth on the default background
and tenth on the correct one.

**The simpler way: hand over the whole table.** Rather than filtering
first and then describing the background, pass the result table and say
which p-value to judge on:

```r
enrichment_CTD(de, method = "ORA", alpha = 0.05, alpha_column = "padj")
```

The genes below the threshold become the list to test and **every row**
becomes the background, the selected ones included: the test draws *n*
genes from an urn of *N*, and the drawn genes were in the urn. The
background is everything you studied, not the non-significant
remainder. Both come from one object, so they cannot disagree. `alpha_column` takes a name or an index, and naming it
matters: on a table from `limma::topTable()` the second column is the
log fold change, not a p-value. Which p-value to judge on is your
decision, not the package's, and the function reports the column it
used.

**Or supply the background yourself**, if you have already selected
your genes:

```r
enrichment_CTD(sig_genes, method = "ORA", universe = de$EntrezID)
```

What belongs in `universe` is every gene that entered your differential
test. Not every gene you sequenced, and not only the significant ones. A
gene filtered out for low expression could not have come out
significant, so it does not belong in the urn either; a gene that was
tested and did not reach significance very much does.

The two are mutually exclusive: with `alpha` the background is already
decided.

**Only ORA needs this argument, and it is worth seeing why.** The three
other methods cannot get the background wrong, because their input
carries it. GSEA ranks the whole list you give it, so the list is the
background. CAMERA and GSVA intersect the gene sets with
`rownames(x)`, so their background is the set of genes you measured, by
construction.

ORA is the exception because a bare gene list records only which genes
came out, not which ones could have. That missing information is what
`universe` supplies, and why its absence is the one input error this
package cannot detect for you. Passing it to GSEA, CAMERA or GSVA
raises a warning rather than being silently dropped.

### Why the size thresholds are what they are

In `ctdR` every chemical *is* a gene set, and CTD gene sets are small:
the median chemical has **4** target genes, against 72 for KEGG and 11
for GO. Size thresholds tuned for those collections do not transfer.
Measured on the full CTD chemical-gene file, of 11,067 chemicals:

| `minGSSize` | Chemicals tested | % |
|---:|---:|---:|
| 1 | 10,802 | 97.6 |
| **2** *(default)* | **7,970** | **72.0** |
| 3 | 6,498 | 58.7 |
| 5 | 4,844 | 43.8 |
| 10 | 2,971 | 26.8 |

A `minGSSize` of 10, the conventional value, would test barely a
quarter of them. The rest are not tested and found unremarkable: they
are absent from the output altogether.

One-gene sets are excluded on purpose rather than by accident. For a
set of one gene the hypergeometric p-value is exactly the number of
input genes divided by the background size, the same value whichever
gene the set contains. It measures whether that gene is in your list,
not whether the chemical is enriched. Removing those sets also makes
the remaining ones easier to call significant, since they no longer add
to the multiple-testing burden.

Chemicals removed by the filter are reported rather than dropped
quietly. `ora()` emits a message giving how many chemicals went
untested and whether they fell below or above the thresholds, so a
chemical missing from the results can be told apart from one that was
tested and came back non-significant.

Raise `minGSSize` when you want only well-characterised chemicals, and
lower it to 1 only if you are deliberately looking for single-gene
associations and will read those p-values accordingly.

#### And why there is no upper limit

The same reasoning applied to the other end reaches the opposite
answer, so `maxGSSize` defaults to `Inf`.

A set of *M* genes cannot, even when every one of the *m* input genes
falls inside it, produce a p-value below *C*(*M*,*m*)/*C*(*N*,*m*),
which is roughly (*M*/*N*)^*m*. That floor **rises** with *M*: a large
enough set becomes untestable, just as a small enough one does. The
question is whether CTD contains any such set.

It does not. With the 28,571 genes CTD covers and an input list of 169
differentially expressed genes, the largest still-testable set is about
**26,600** genes, 93% of the universe. The largest chemical in CTD,
benzo(a)pyrene, has 16,536. Every chemical in the database is testable
from above, so an upper cut removes sets that could have been declared
significant.

What a cut removes is also not marginal. At 500 it excludes 265
chemicals, among them benzo(a)pyrene, valproic acid, sodium arsenite,
bisphenol A, aflatoxin B1 and particulate matter. Their sets are large
because the literature on them is large. This is where CTD differs from
GO: there a large term is one that has stopped meaning anything, here a
large set is a chemical somebody has studied for decades. Keeping them
costs 3% more tests, 8,235 against 7,970.

Set `maxGSSize` to a finite value if you have a reason of your own, for
instance to keep the focus on chemicals with sharply defined targets.
`ctdR` does not impose one.

#### What that costs, and how to read around it

Removing the cap has a consequence worth stating rather than
discovering: with no upper limit, the significant results skew towards
large sets. On the analysis in this vignette's tutorial, the median
chemical among the significant hits carries about 4,100 target genes,
against a median of 7 across every chemical tested.

This is worth understanding rather than memorising, because the fold
enrichment appears to say the opposite. Fold is
`(k/n) / (M/N)`: the set size *M* sits in the denominator, so for a
fixed overlap *k* a larger set gives a **lower** fold. Yet larger sets
reach smaller p-values. Both are true, because the two measure
different things.

The p-value does not measure the ratio. It measures how unlikely the
excess is, and the excess is counted in genes, not in proportion. With
154 input genes drawn from a background of 27,444, holding the fold
fixed at 1.5:

| Genes in set | Expected by chance | Observed at fold 1.5 | p |
|---:|---:|---:|---:|
| 100 | 0.6 | 1 | 0.43 |
| 2,000 | 11.2 | 17 | 0.057 |
| 16,000 | 89.8 | 135 | 1.4e-15 |

The same 50% enrichment spans fifteen orders of magnitude. A set of 100
genes at fold 1.5 is 0.4 genes above expectation, which happens
constantly by chance; a set of 16,000 at the same fold is 45 genes
above, which does not. The identical ratio carries incomparable
evidence, and the test is right to say so.

In CTD, set size also measures **how much a chemical has been
studied**: the database curates published literature, so a chemical
measured in many genome-wide experiments accumulates many recorded
genes. It is tempting to read that as a bias, and to suspect that large
sets are padded with genes that respond to everything. Measured, they
are not. The genes in sets above 500 appear in a median of 47 chemicals
each; the genes in sets of 4 or fewer appear in a median of 232. Small
sets are the ones made of the usual suspects, because a chemical
studied once was studied with a targeted assay aimed at a gene somebody
already suspected. Large sets are, per gene, more specific rather than
less.

So finding well-studied chemicals near the top is not an artefact to be
corrected. The evidence behind them is real, and no method can find a
chemical nobody has measured: that is a limit of the available data,
not of the test.

What it does change is what a result means, and this holds for every
enrichment analysis run against a curated database, not only this one.
The question the test answers is not *is this chemical associated with
these genes*, but **is this chemical associated with these genes, as
far as the published literature records**. Between two chemicals of
equal biological relevance, the better-studied one has more opportunity
to overlap your gene list and ranks higher. A chemical absent from the
output is therefore uninvolved only as far as anyone currently knows:
the analysis measures the state of the evidence, and reports the
biology through it.

The practical consequence is how to read the output. Results are sorted
by adjusted p-value, which is the quantity the false discovery rate is
controlled on. Fold enrichment is **not** a sound sort key on its own:
it has no error control, and a two-gene set with both genes hit would
sit at the top of it. Read the two columns together:

| Chemical | Genes in set | Overlap | Fold | Adjusted p |
|---|---:|---:|---:|---:|
| Benzo(a)pyrene | 16,337 | 124 | 1.35 | 4.3e-05 |
| GSK-J4 | 2,331 | 35 | 2.68 | 8.0e-05 |

Both are significant, and the first ranks higher. The second is the
sharper association. Neither column answers the question alone.

`FoldEnrichment`, `Count` and `BackgroundRatio` are in the output for
this reason. Sorting by `FoldEnrichment` **within** the chemicals that
pass your FDR threshold is a reasonable second pass, and gives a
different and equally legitimate reading of the same result.

```{r size-universe}
# ORA: restrict the background to expressed genes. Both size thresholds
# are left at their defaults; maxGSSize is shown here only to make the
# point that capping it is a deliberate choice, not a starting position.
ora_expressed <- enrichment_CTD(
    genes, method = "ORA",
    universe   = c(genes$EntrezID, "7422", "836")  # expressed genes
)

# GSEA: filter gene sets by size
gsea_filtered <- enrichment_CTD(
    genes, method = "GSEA",
    minSize = 3,
    maxSize = 300
)
```

Setting `universe` in ORA to the full list of genes measured in your
experiment (rather than the default of all CTD-annotated genes) is
strongly recommended: it avoids inflated fold-enrichment estimates that
arise when the background is larger than the actual measurement space.

# Choosing among the four methods

| Aspect | ORA | GSEA | CAMERA | GSVA |
|---|---|---|---|---|
| **Input** | Gene list | Ranked gene list | Expression matrix + design + contrast | Expression matrix |
| **Question** | Over-represented? | Cluster at extremes? | Stronger signal than rest of transcriptome (correlation-corrected)? | Per-sample chemical score |
| **Output** | Ranked chemicals (data.frame) | Ranked chemicals (data.frame) | Ranked chemicals (data.frame) | Chemical x sample score matrix |
| **Best for** | DEG lists | Exploratory, ranked | Multi-sample DE experiments, co-regulated gene sets | Patient stratification, downstream tests, heatmaps |

Use **ORA** when you have a well-defined gene list above a
significance threshold; **GSEA** to leverage the full
ranking without a cutoff; **CAMERA** when you have a
proper multi-sample design and want to control the
inflated false-positive rate that ORA/GSEA exhibit on
strongly co-regulated gene sets; **GSVA** when you want a
per-sample score to feed into downstream analyses
(clustering, survival, correlation with phenotypes).

# Interoperability with existing Bioconductor infrastructure

`ctdR`'s contribution is the CTD chemical-gene sets and a unified
interface; the statistical tests themselves are delegated to
established engines (`fgsea`, `limma`, `GSVA`); ORA is a
hypergeometric test computed on `stats::phyper`.
If you would rather run the CTD gene sets through a different engine,
`ctdR` exposes them directly.

## Reading a CTD file the `BiocIO` way

Besides `import_CTD()` (which processes and caches the data), a CTD
file can be read into memory using the standard Bioconductor `import()`
generic via the `CTDFile` class -- so the file behaves like any other
`BiocFile`:

```{r ctdfile, eval=FALSE}
library(BiocIO)
ctd <- import(CTDFile("~/Downloads/CTD_chem_gene_ixns.csv"))
head(ctd)
```

## Handing the CTD gene sets to another engine

`as_genesets_CTD()` returns the chemical gene sets as a named list, in
the exact shape expected by, for example, the `gs` argument of
`EnrichmentBrowser::sbea()`. `ctdR` does **not** depend on
`EnrichmentBrowser` -- you supply it.

```{r asgenesets}
gene_sets <- as_genesets_CTD("entrez")
length(gene_sets)
gene_sets[[1]]
```

```{r sbea, eval=FALSE}
# Illustrative: route the CTD gene sets through EnrichmentBrowser's
# set-based engine. sbea() covers all four of ctdR's methods, GSVA
# included, plus eight more (see EnrichmentBrowser::sbeaMethods()).
library(EnrichmentBrowser)
res <- sbea(method = "camera", se = my_se, gs = as_genesets_CTD("entrez"))
gsRanking(res)
```

The named list is also a one-line conversion away from a
`GeneSetCollection` (`GSEABase::GeneSetCollection()`), for tools that
expect that container.

# Naming and scope: the "CTD" acronym

In `ctdR`, **CTD** always means the *Comparative Toxicogenomics
Database* (<https://ctdbase.org>). The same acronym denotes unrelated
things elsewhere in the R ecosystem: the CRAN package **CTD**
("Connect The Dots") is a graph algorithm for metabolomic and
transcriptomic perturbations, and in **EWCE** a `ctd` object is a
*CellTypeDataset*. `ctdR` is unrelated to both. The `CTDFile` class is
named to avoid this collision -- and to avoid `CTDData`, the class of
the CTD-specific package `CTDquerier`, which was removed from
Bioconductor at release 3.12.

# Updating the CTD data

CTD releases updated data periodically. To update:

1. Download the latest `CTD_chem_gene_ixns.csv.gz` from
    <https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz>.
2. Decompress and re-run `import_CTD()` -- existing
    cache files are overwritten.

```{r reimport, eval=FALSE}
import_CTD("~/Downloads/CTD_chem_gene_ixns.csv")
```

# References

`ctdR` supplies the CTD gene sets and a common interface; the
statistical tests are the work of others, and the methods it uses are
published. Cite the ones your analysis relied on, not the whole list.

**The data.** Every result depends on the CTD release you imported, and
`ctd_provenance()` reports which one that was.

- Davis AP, Wiegers TC, Johnson RJ, Sciaky D, Wiegers J, Mattingly CJ
  (2023). Comparative Toxicogenomics Database (CTD): update 2023.
  *Nucleic Acids Research* 51(D1), D1257-D1262.
  doi:[10.1093/nar/gkac833](https://doi.org/10.1093/nar/gkac833)

**The methods.** One reference per method, plus the implementation where
it differs from the paper that introduced it.

- **GSEA**, the method: Subramanian A, Tamayo P, Mootha VK, *et al.*
  (2005). Gene set enrichment analysis: a knowledge-based approach for
  interpreting genome-wide expression profiles. *PNAS* 102(43),
  15545-15550.
  doi:[10.1073/pnas.0506580102](https://doi.org/10.1073/pnas.0506580102)
- **fgsea**, the implementation `ctdR` calls: Korotkevich G, Sukhov V,
  Budin N, Shpak B, Artyomov MN, Sergushichev A (2021). Fast gene set
  enrichment analysis. *bioRxiv*.
  doi:[10.1101/060012](https://doi.org/10.1101/060012)
- **CAMERA**: Wu D, Smyth GK (2012). Camera: a competitive gene set test
  accounting for inter-gene correlation. *Nucleic Acids Research* 40(17),
  e133. doi:[10.1093/nar/gks461](https://doi.org/10.1093/nar/gks461)
- **limma**, which provides CAMERA: Ritchie ME, Phipson B, Wu D, Hu Y,
  Law CW, Shi W, Smyth GK (2015). limma powers differential expression
  analyses for RNA-sequencing and microarray studies. *Nucleic Acids
  Research* 43(7), e47.
  doi:[10.1093/nar/gkv007](https://doi.org/10.1093/nar/gkv007)
- **GSVA**: Hänzelmann S, Castelo R, Guinney J (2013). GSVA: gene set
  variation analysis for microarray and RNA-seq data. *BMC
  Bioinformatics* 14, 7.
  doi:[10.1186/1471-2105-14-7](https://doi.org/10.1186/1471-2105-14-7)

**Interoperability.** The compatibility section above routes the CTD
gene sets through another engine, which `ctdR` does not depend on:

- Geistlinger L, Csaba G, Zimmer R (2016). Bioconductor's
  EnrichmentBrowser: seamless navigation through combined results of
  set- & network-based enrichment analysis. *BMC Bioinformatics* 17, 45.
  doi:[10.1186/s12859-016-0884-1](https://doi.org/10.1186/s12859-016-0884-1)

**The reasoning behind the size thresholds.** The argument for excluding
one-gene sets, and for not excluding large ones, rests on this:

- Bourgon R, Gentleman R, Huber W (2010). Independent filtering
  increases detection power for high-throughput experiments. *PNAS*
  107(21), 9546-9551.
  doi:[10.1073/pnas.0914005107](https://doi.org/10.1073/pnas.0914005107)

**The example data.** The RNA-seq analysis shown here uses a subset of:

- NCBI Gene Expression Omnibus, accession
  [GSE311566](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE311566):
  Comparing the effects of PFAS and dexamethasone on human PBMCs (2025).

**`ctdR` itself.** `citation("ctdR")` prints the current form, which
includes the version you ran and the archival DOI.

# Session info

```{r session-info}
sessionInfo()
```