Search for gene sets that overlap with genes in a user-supplied list.
searchOverlappingSets(
species,
genes,
counts.only = TRUE,
test.enrichment = TRUE,
config = NULL
)String containing the NCBI taxonomy ID of the species of interest.
Integer vector containing gene indices.
Each gene index refers to a row of the data frame returned by fetchAllGenes.
Duplicate entries are ignored.
Boolean indicating whether to only report the number of overlapping genes for each set.
Boolean indicating whether to compute a hypergeometric p-value for enrichment of genes in each set.
Configuration list, typically created by newConfig.
If NULL, the default configuration is used.
A list containing a overlap data frame and a present integer.
In the overlap data frame, each row represents a set that overlaps with genes.
The data frame contains the following columns:
set, an integer column containing the set index.
This refers to a row of the data frame returned by fetchAllSets.
count, an integer column containing the number of overlaps between the genes in the set and those in genes.
genes, a nested list that contains the entries of genes that overlap with those in the set.
Only reported if counts.only = FALSE.
size column, an integer column containing the size of each set.
Only reported if test.enrichment = TRUE, as it is a by-product of the p-value calculation.
pvalue column, a numeric column containing the hypergeometric p-value for overrepresentation of genes in the set.
Only reported if test.enrichment = TRUE.
The row order is arbitrary.
present specifying the number of genes in genes that are present in at least one set in the Gesel database for species.
present can be used as the number of draws when performing a hypergeometric test for gene set enrichment, instead of length(genes) (see Details).
This ensures that we do not consider genes that are not present in any gene sets in Gesel,
e.g., due to changes in annotation across genome versions or because they are pseudogenes or predicted genes.
Otherwise, unknown genes would inappropriately increase the number of draws and inflate the enrichment p-value.
out <- searchOverlappingSets("9606", 1:10)
overlaps <- out$overlap
head(overlaps)
#> set count size pvalue
#> 1 22 1 82 0.0186933885
#> 2 70 1 1826 0.3487621766
#> 3 368 1 284 0.0634086424
#> 4 451 1 176 0.0397350513
#> 5 584 1 2 0.0004597226
#> 6 681 1 109 0.0247793410
# More details on the overlapping sets.
all.sets <- fetchAllSets("9606")
all.sets[head(overlaps$set),]
#> name description
#> 22 GO:0000049 tRNA binding
#> 70 GO:0000166 nucleotide binding
#> 368 GO:0001525 angiogenesis
#> 451 GO:0001666 response to hypoxia
#> 584 GO:0001869 negative regulation of complement activation, lectin pathway
#> 681 GO:0002020 protease binding
#> size collection number
#> 22 82 1 22
#> 70 1826 1 70
#> 368 284 1 368
#> 451 176 1 451
#> 584 2 1 584
#> 681 109 1 681
# Computing the enrichment p-value manually. We take the upper tail after
# subtracting 1 to ensure that the probability mass of the observed
# number of overlapping genes is included in the p-value.
set.size <- all.sets$size[overlaps$set]
universe <- effectiveNumberOfGenes("9606")
p <- phyper(
q = overlaps$count - 1,
m = set.size,
n = universe - set.size,
k = out$present,
lower.tail=FALSE
)
stopifnot(identical(p, overlaps$pvalue))
# For multiple testing correction, it is necessary to consider all sets
# in the database, as these were implicitly considered during the search
# though only a subset of them are reported by searchOverlappingSets.
fdr <- p.adjust(p, method="BH", n=nrow(all.sets))
summary(fdr <= 0.05)
#> Mode FALSE TRUE
#> logical 1770 34