Over-Representation Analysis (ORA)

Internal engine that performs Over-Representation Analysis with a hypergeometric test computed directly on phyper. Returns a result table with one row per tested chemical.

Column renaming, multiple-testing correction, chemical-name join, canonical ordering and sort are applied downstream by .format_enrichment_result so the engine stays close to the vocabulary of the test itself.

Usage

ora(
  ChemicalName_GeneSymbols,
  gene_symbols,
  pAdjustMethod = "BH",
  universe = NULL,
  minGSSize = 2,
  maxGSSize = Inf
)

Arguments

ChemicalName_GeneSymbols

A data frame with two columns (term, gene) mapping CTD chemical IDs to HGNC gene symbols. The first two columns are used, whatever their names.

gene_symbols

Character vector of HGNC gene symbols to test for enrichment. Duplicates and NA are removed.

pAdjustMethod

Character. Method for multiple testing correction (default "BH"). Passed to p.adjust.

universe

Vector of background gene identifiers, or NULL (default) to use every gene in ChemicalName_GeneSymbols. Set it to rownames(expr) or to the full tested gene list to restrict the background to measured genes only. Coerced with as.character(), so an integer column of Entrez IDs works as it does for gene_symbols.

The default is a fallback, not a recommendation, and the function says so when it uses it. A hypergeometric test asks how many of \(n\) genes drawn from \(N\) would land in a set. Genes that your experiment could never have detected still sit in \(N\), filling the urn with balls that cannot be drawn, so the draw looks more selective than it was and the p-value comes out too small. The error is anti-conservative: it manufactures significance.

The right universe is every gene that entered your test, not every gene you sequenced and not only the significant ones: a gene filtered out for low expression could not have been selected, so it does not belong in the urn either. On the RNA-seq example bundled with this package the difference is 32 significant chemicals against 19.

minGSSize

Integer. Minimum gene set size after intersection with the background (default 2). The default is chosen for CTD, where the median chemical has 4 target genes: a one-gene set is degenerate, since its p-value equals the ratio of input genes to background whichever gene it contains, so it measures membership rather than enrichment.

maxGSSize

Maximum gene set size after intersection with the background. The default is Inf, that is no upper limit, and that too is a choice made for CTD rather than inherited.

The reasoning mirrors the one for minGSSize, and reaches the opposite conclusion. 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) \approx (M/N)^m\). That floor rises with \(M\), so a large enough set becomes untestable. Measured on CTD, with 28,571 genes in the background and an input list of 169, the largest still-testable set is about 26,600 genes, 93% of the universe. The largest chemical in CTD has 16,536. No chemical is untestable from above, so an upper cut removes sets that could have been declared significant.

What it removes is not marginal. A cut at 500 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: in CTD size tracks how well studied a chemical is, where in GO a large term is one that has stopped meaning anything. Keeping them costs 3% more tests, 8,235 against 7,970.

Set it to a finite value if you have a reason of your own; ctdR does not impose one.

One consequence to be aware of: without a cap the significant results skew towards large sets. This is not the fold enrichment talking, which moves the other way, since \(M\) sits in its denominator. It is that a p-value measures how unlikely an excess is, and the excess is counted in genes. With 154 input genes from a background of 27,444, a fold of 1.5 means 0.4 genes above expectation for a set of 100 (p = 0.43) and 45 genes above for a set of 16,000 (p = 1.4e-15). Sort by p.adjust to rank by evidence, and read foldEnrichment beside it to see how sharp the association is; neither answers the question alone.

Value

A data frame with columns ChemicalID, GeneRatio, BgRatio, pvalue, p.adjust, geneID, Count and foldEnrichment, sorted by pvalue ascending. Returns an empty data frame with the same structure when no gene set can be tested. Emits a message reporting how many chemicals the size filter left untested, since those are absent from the result rather than present with a large p-value.

Details

Earlier versions delegated this step to the clusterProfiler package. That pulled in 59 packages, and a visualization layer this package never used, to reach a single call. The test itself is one line of stats, so it is computed here and that dependency is gone. The p-values are unchanged: the two implementations were compared over 24 configurations of input list, minimum set size and background universe, and the largest absolute difference was exactly 0.

The background is every gene appearing in ChemicalName_GeneSymbols, optionally narrowed by universe. Gene sets are intersected with that background before the size filter, so minGSSize and maxGSSize always refer to the set as actually tested rather than to its nominal size.