Hypergeometric test in R

Hypergeometric test

Background

The hypergeometric test is a simple statistical procedure to test for over- or under-representation of items (groups) in a defined global set. For example, if a global set consists of 50 black and 50 white balls, and 20 balls are randomly selected from this set, what is the probability that 17 white balls and only 3 black balls are chosen? The hypergeometric test can answer such questions by using the probability function of the hypergeometric distribution. This distribution is highly related to the binomial distribution which uses a similar definition but with replacement, instead of without.

This document walks through one example where the hypergeometric test is used to determine gene enrichment for pathways of the [subtiwiki pathway categories](http://www.subtiwiki.uni-goettingen.de/v3/category/ for the bacterium Bacillus subtilis.

Libraries and test data

• Install Category package from Bioconductor

BiocManager::install("Category")

• load required libraries

suppressPackageStartupMessages({
  library(Category)
  library(tidyverse)
})

• import category test data (from subtiwiki)

df_subtiwiki <- read_csv("../data/geneCategories-2022-09-02.csv",
  name_repair = function(x) str_replace_all(x, " ", "_"),
  show_col_types = FALSE)

head(df_subtiwiki)

• convert dataframe into a nested list format

subti_categories <- df_subtiwiki %>%
  mutate(category = str_sub(category, 1, 25)) %>%
  select(category, gene_locus) %>%
  group_by(category) %>%
  summarize(gene_locus = list(gene_locus)) %>%
  deframe()

head(subti_categories, n = 3)

$`6S RNA`
[1] "BSU_misc_RNA_32" "BSU_misc_RNA_41"

$`A/P endonucleases`
[1] "BSU_22340" "BSU_25130" "BSU_40880"

$`ABC transporters for the `
 [1] "BSU_01610" "BSU_01610" "BSU_01620" "BSU_01620" "BSU_01630" "BSU_01630" "BSU_03800" "BSU_03800"
 [9] "BSU_03810" "BSU_03810" "BSU_03820" "BSU_03820" "BSU_03830" "BSU_03830" "BSU_07490" "BSU_07490"
[17] "BSU_07500" "BSU_07500" "BSU_07510" "BSU_07510" "BSU_07520" "BSU_07520" "BSU_08440" "BSU_08440"
[25] "BSU_08450" "BSU_08450" "BSU_08460" "BSU_08460" "BSU_10330" "BSU_10330" "BSU_32940" "BSU_32940"
[33] "BSU_33290" "BSU_33290" "BSU_33300" "BSU_33300" "BSU_33310" "BSU_33310" "BSU_33320" "BSU_33320"
[41] "BSU_39610" "BSU_39610"

Prepare input data

• assayed: gene sets that were tested in an experiment, usually all genes/pathways in the universe
• example: we select a subset of 20 subtiwiki categories with each more than 20 genes

assayed <- subti_categories[lapply(subti_categories, length) >= 20][1:20]

• significant: genes within the assayed set that were observed, or fall below a sign. threshold
• example: we select 10% of genes as significant

significant <- lapply(assayed, function(x) sample(x, length(x)/10))

• universe: all unique gene names

universe = assayed %>%
  unlist(use.names = FALSE) %>%
  unique

Run test and inspect results

• run hypergeometric test using hyperg function from Categories

result <- hyperg(
  assayed,
  significant,
  universe,
  representation = "over"
)

result <- mutate(result, across(everything(), unlist))

• expected result: all p-values are very similar (no enrichment of a specific pathway)

head(result)

• what happens if one pathway has more than 10% sign. genes?

significant[[1]] <- sample(assayed[[1]], 20)

• run hypergeometric test again and inspect result
• the first (modified) pathway has lower p-value and higher odds

result2 <- hyperg(
    assayed,
    significant,
    universe,
    representation = "over"
  ) %>% 
  mutate(across(everything(), unlist))

head(result2)

P-value adjustment

• p-values from testing hundreds to thousands of pathways can be highly inflated
• the chance of discovering false positive correlations increases with the number of carried out tests
• to correct for multiple hypothesis testing, p-value adjustment is often advisable
• simply speaking, p-value adjustment works by correcting the observed p-value distribution depending on the number of tests that were carried out
• various methods are available, differing in their conservativeness
• a good option to start with is usually the Benjamini-Hochberg method, implemented in R with stats::p.adjust(p, method = "BH")
• the example below shows how p-values increase (become less “significant”) after adjustment

result2 <- mutate(result2, p_adj = p.adjust(p, method = "BH"))

result2 %>%
  ggplot(aes(x = p, y = p_adj)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = grey(0.5)) +
  geom_point() +
  theme_light() +
  theme(aspect.ratio = 1.0)

Interpretation

The output table from the gene enrichment function can be interpreted as follows:

• p: column with un-adjusted p-values
• odds: intermediate statistic from hypergeometric test – if we picture an urn containing genes from the universe, genes annotated at a given category are white (W) and the rest black (B), then odds ratio is (W_drawn * (B - B_drawn)) / (B_drawn * (W - W_drawn))
• expected: the number of expected genes from the assayed group