| Title: | Taxon Set Enrichment Analysis |
|---|---|
| Description: | TaxSEA makes metagenomic data easier to interpret. Rather than reading species one at a time, it tests whether a-priori defined sets of taxa shift together, drawing on public databases (BugSigDB, MiMeDB, GutMGene, mBodyMap, BacDive and GMRepoV2) and on sets collated from the literature. Two approaches are provided. TaxSEA takes taxon names or NCBI IDs with a rank per taxon, such as a fold change or correlation coefficient, and tests each set with a Kolmogorov-Smirnov test. ssTaxSEA scores every set in a single sample as the mean centered log-ratio of its members, giving a per-sample value to use in any downstream model. Both work with any taxonomic profiling technology, including 16S rRNA gene sequencing, shotgun metagenomics and metatranscriptomics. |
| Authors: | Feargal Ryan [aut, cre, fnd] (ORCID: <https://orcid.org/0000-0002-1565-4598>, funding: Supported by NHMRC Investigator Grant) |
| Maintainer: | Feargal Ryan <[email protected]> |
| License: | GPL-3 |
| Version: | 1.5.8 |
| Built: | 2026-10-26 23:27:42 UTC |
| Source: | https://github.com/bioc/TaxSEA |
This function takes a vector of taxon names and returns a vector of NCBI taxonomy IDs by querying the NCBI Entrez API.
get_ncbi_taxon_ids(taxon_names)get_ncbi_taxon_ids(taxon_names)
taxon_names |
A character vector of taxon names |
A character vector of NCBI taxonomy IDs corresponding to the input taxon names
taxon_names <- c("Escherichia coli", "Staphylococcus aureus") taxon_ids <- get_ncbi_taxon_ids(taxon_names)taxon_names <- c("Escherichia coli", "Staphylococcus aureus") taxon_ids <- get_ncbi_taxon_ids(taxon_names)
Retrieve from the TaxSEA database which taxon sets (metabolite producers and disease signatures) contain a taxon of interest.
get_taxon_sets(taxon_to_fetch)get_taxon_sets(taxon_to_fetch)
taxon_to_fetch |
The taxon to search for in the TaxSEA database. |
A character vector containing the names of taxonomic sets where the specified taxon is present.
# Retrieve sets for Bifidobacterium longum get_taxon_sets(taxon_to_fetch = "Bifidobacterium_longum")# Retrieve sets for Bifidobacterium longum get_taxon_sets(taxon_to_fetch = "Bifidobacterium_longum")
A lookup from taxon names to NCBI taxonomy IDs, used to translate the
names supplied to TaxSEA() into the IDs that TaxSEA_db is
built on. Names are present both with spaces and with underscores, and
each ID is also present as its own name so that IDs can be supplied
directly. Species names, former names and IDs of the BacDive_ set
members were added from the BacDive harvest without changing any existing
entry.
data(NCBI_ids)data(NCBI_ids)
A named list where:
Taxon names (e.g. "Bifidobacterium breve", "Bifidobacterium_breve") or NCBI IDs
NCBI taxonomy IDs, as character
NCBI Taxonomy; BacDive (https://bacdive.dsmz.de).
data(NCBI_ids) # Can look up either with or without spaces NCBI_ids["Bifidobacterium_breve"] NCBI_ids["Bifidobacterium breve"]data(NCBI_ids) # Can look up either with or without spaces NCBI_ids["Bifidobacterium_breve"] NCBI_ids["Bifidobacterium breve"]
Computes a per-sample score for each taxon set as the mean centered log-ratio (CLR) of the set's members within that sample.
ssTaxSEA( counts, lookup_missing = FALSE, min_set_size = 5, max_set_size = 300, custom_db = NULL, pseudocount = NULL, scale = c("counts", "proportion"), reference = c("all", "top_k"), top_k = 20 )ssTaxSEA( counts, lookup_missing = FALSE, min_set_size = 5, max_set_size = 300, custom_db = NULL, pseudocount = NULL, scale = c("counts", "proportion"), reference = c("all", "top_k"), top_k = 20 )
counts |
A numeric matrix, data.frame, or
|
lookup_missing |
Logical indicating whether to fetch missing NCBI IDs via the NCBI API. Default is FALSE. |
min_set_size |
Minimum size of taxon sets to include. Default is 5. |
max_set_size |
Maximum size of taxon sets to include. Default is 300. |
custom_db |
A user-provided list of taxon sets. If NULL (default), the built-in TaxSEA database is used (excluding BugSigDB). |
pseudocount |
Numeric value added to every value before the
log transform, to handle zeros. It is on the scale set by
|
scale |
Either |
reference |
Either |
top_k |
Number of most abundant detected taxa used as the
reference when |
For each sample the counts are centered log-ratio transformed, and the score for a taxon set is the mean of the CLR values of the set's members:
The CLR is computed across all taxa supplied in counts,
before any subsetting to set members. This matters: restricting the
matrix to set members first would make the geometric mean that the
CLR divides by depend on the sets being tested, which manufactures
apparent signal in sets that have none.
Because the score is a within-sample quantity, ssTaxSEA does
not need a cohort and is well defined for a single sample.
A numeric matrix of scores with taxon sets as rows and
samples as columns. Higher values indicate that the set's members
are more abundant in that sample, relative to the average taxon
in that same sample (or, with reference = "top_k", relative
to that sample's most abundant taxa).
The standard CLR (reference = "all") centers each sample on
the mean log value of every taxon in the table, zeros included. That
reference moves with richness: when many rare taxa drop out, as after
antibiotics or chemotherapy, the reference falls and every surviving
taxon's score rises, even if its abundance relative to the other
survivors has not changed. Many taxon sets share the same surviving
members, so hundreds of scores can rise together. Adjusting for
richness in a downstream model does not remove this. The reference
also depends on which rows are in the table, so filtering taxa
differently changes every score.
reference = "top_k" centers each sample on the mean log value
of its top_k most abundant detected taxa. Rare taxa are still
scored and still count towards their sets; they just no longer move
everyone else's score. Losing rare taxa, or adding or removing rare
and all-zero rows, leaves the reference unchanged. Zeros never enter
the reference: a sample with fewer than top_k detected taxa
uses all of them. The score then reads as abundance relative to the
sample's dominant taxa.
A pseudocount set from each sample's own data would move with the
rare tail in the same way, which is why pseudocount is a fixed
value. On the count scale a pseudocount of 0.5 reads is a different
proportion in samples of different depth; scale = "proportion"
with a small fixed pseudocount avoids that. The default of 1e-6 is
near the detection limit of deep shotgun profiles; raise it if your
smallest detected proportions are larger. For relative abundance
input, such as MetaPhlAn percentages, use scale = "proportion".
Scores are comparable across samples within a taxon set: a higher score in sample A than sample B means the set's members make up more of sample A's community.
Scores are not comparable across taxon sets. A set of
abundant taxa will score higher than a set of rare taxa in every
sample, regardless of biology, simply because its members are more
abundant. When visualising the matrix, center or scale the rows
first (for example t(scale(t(scores)))), and when comparing
groups, compare a single set's scores between groups rather than
comparing different sets to each other.
No p-values are returned. The score is a descriptive statistic;
to test a hypothesis, compare scores across samples using a test
appropriate to your design (for example wilcox.test
or a linear model on one row of the returned matrix).
TaxSEA for group-level enrichment.
# Toy count matrix: 30 taxa x 8 samples set.seed(42) counts <- matrix(rpois(240, lambda = 10), nrow = 30, ncol = 8) rownames(counts) <- paste0("Taxon_", seq_len(30)) colnames(counts) <- paste0("Sample_", seq_len(8)) # Raise one set's members in the last four samples counts[1:6, 5:8] <- counts[1:6, 5:8] * 5 sets <- list( elevated_set = paste0("Taxon_", 1:6), control_set = paste0("Taxon_", 20:27) ) scores <- ssTaxSEA(counts, custom_db = sets, min_set_size = 3) dim(scores) # sets x samples round(scores, 2) # The planted signal shows up as higher scores in samples 5-8 rowMeans(scores[, 1:4]) - rowMeans(scores[, 5:8]) # A reference that is not moved by changes in rare taxa, for # relative abundance input pct <- sweep(counts, 2, colSums(counts), "/") * 100 ssTaxSEA(pct, custom_db = sets, min_set_size = 3, scale = "proportion", reference = "top_k", top_k = 10)# Toy count matrix: 30 taxa x 8 samples set.seed(42) counts <- matrix(rpois(240, lambda = 10), nrow = 30, ncol = 8) rownames(counts) <- paste0("Taxon_", seq_len(30)) colnames(counts) <- paste0("Sample_", seq_len(8)) # Raise one set's members in the last four samples counts[1:6, 5:8] <- counts[1:6, 5:8] * 5 sets <- list( elevated_set = paste0("Taxon_", 1:6), control_set = paste0("Taxon_", 20:27) ) scores <- ssTaxSEA(counts, custom_db = sets, min_set_size = 3) dim(scores) # sets x samples round(scores, 2) # The planted signal shows up as higher scores in samples 5-8 rowMeans(scores[, 1:4]) - rowMeans(scores[, 5:8]) # A reference that is not moved by changes in rare taxa, for # relative abundance input pct <- sweep(counts, 2, colSums(counts), "/") * 100 ssTaxSEA(pct, custom_db = sets, min_set_size = 3, scale = "proportion", reference = "top_k", top_k = 10)
Groups species by taxonomic ranks and performs TaxSEA enrichment analysis at each rank level. Returns a named list of data frames, one per taxonomic rank (excluding species).
taxon_rank_sets(taxon_ranks, lineage_df, min_set_size = 5, max_set_size = 100)taxon_rank_sets(taxon_ranks, lineage_df, min_set_size = 5, max_set_size = 100)
taxon_ranks |
A named numeric vector of log2 fold changes.
Names should be feature identifiers matching the
|
lineage_df |
Either a data frame or a
Data frame input: Must include a SummarizedExperiment input: Taxonomy is extracted from
|
min_set_size |
Minimum number of species in a set to include in the analysis. Default is 5. |
max_set_size |
Maximum number of species in a set to include in the analysis. Default is 100. |
A named list of data frames, one per taxonomic rank. Each data frame contains columns: taxonSetName, median_rank_of_set_members, PValue, Test_statistic, and FDR. Ranks that produce no valid sets (e.g., due to size filtering) are included as empty data frames with a message.
# --- Example 1: Data frame input --- # Create a lineage data frame (e.g., parsed from curatedMetagenomicData) # The 'species' column must match the names in taxon_ranks. lineage_df <- data.frame( species = c("Cutibacterium_acnes", "Klebsiella_pneumoniae", "Propionibacterium_humerusii", "Moraxella_osloensis", "Enhydrobacter_aerosaccus", "Staphylococcus_capitis", "Staphylococcus_epidermidis", "Staphylococcus_aureus", "Escherichia_coli", "Enterobacter_cloacae", "Pseudomonas_aeruginosa", "Acinetobacter_baumannii", "Lactobacillus_rhamnosus", "Lactobacillus_acidophilus", "Bifidobacterium_longum", "Bifidobacterium_breve"), kingdom = rep("Bacteria", 16), phylum = c("Actinobacteria", "Proteobacteria", "Actinobacteria", "Proteobacteria", "Proteobacteria", "Firmicutes", "Firmicutes", "Firmicutes", "Proteobacteria", "Proteobacteria", "Proteobacteria", "Proteobacteria", "Firmicutes", "Firmicutes", "Actinobacteria", "Actinobacteria"), class = c("Actinobacteria", "Gammaproteobacteria", "Actinobacteria", "Gammaproteobacteria", "Alphaproteobacteria", "Bacilli", "Bacilli", "Bacilli", "Gammaproteobacteria", "Gammaproteobacteria", "Gammaproteobacteria", "Gammaproteobacteria", "Bacilli", "Bacilli", "Actinobacteria", "Actinobacteria"), order = c("Propionibacteriales", "Enterobacterales", "Propionibacteriales", "Pseudomonadales", "Rhodospirillales", "Bacillales", "Bacillales", "Bacillales", "Enterobacterales", "Enterobacterales", "Pseudomonadales", "Pseudomonadales", "Lactobacillales", "Lactobacillales", "Bifidobacteriales", "Bifidobacteriales"), family = c("Propionibacteriaceae", "Enterobacteriaceae", "Propionibacteriaceae", "Moraxellaceae", "Rhodospirillaceae", "Staphylococcaceae", "Staphylococcaceae", "Staphylococcaceae", "Enterobacteriaceae", "Enterobacteriaceae", "Pseudomonadaceae", "Moraxellaceae", "Lactobacillaceae", "Lactobacillaceae", "Bifidobacteriaceae", "Bifidobacteriaceae"), genus = c("Cutibacterium", "Klebsiella", "Cutibacterium", "Moraxella", "Enhydrobacter", "Staphylococcus", "Staphylococcus", "Staphylococcus", "Escherichia", "Enterobacter", "Pseudomonas", "Acinetobacter", "Lactobacillus", "Lactobacillus", "Bifidobacterium", "Bifidobacterium"), stringsAsFactors = FALSE ) set.seed(42) fc <- setNames(rnorm(16), lineage_df$species) results <- taxon_rank_sets(fc, lineage_df, min_set_size = 2) names(results) results$family # --- Example 2: SummarizedExperiment / TreeSummarizedExperiment input --- library(mia) data(GlobalPatterns, package = "mia") tse <- GlobalPatterns # Run differential abundance (e.g., ALDEx2) to get fold changes # aldex_out <- ... (your DA analysis) # fc <- aldex_out$effect # names(fc) <- rownames(aldex_out) # Run taxon rank set enrichment directly from the TSE results <- taxon_rank_sets(fc, tse, min_set_size = 5) names(results) # Kingdom, Phylum, Class, Order, Family, Genus results$Family# --- Example 1: Data frame input --- # Create a lineage data frame (e.g., parsed from curatedMetagenomicData) # The 'species' column must match the names in taxon_ranks. lineage_df <- data.frame( species = c("Cutibacterium_acnes", "Klebsiella_pneumoniae", "Propionibacterium_humerusii", "Moraxella_osloensis", "Enhydrobacter_aerosaccus", "Staphylococcus_capitis", "Staphylococcus_epidermidis", "Staphylococcus_aureus", "Escherichia_coli", "Enterobacter_cloacae", "Pseudomonas_aeruginosa", "Acinetobacter_baumannii", "Lactobacillus_rhamnosus", "Lactobacillus_acidophilus", "Bifidobacterium_longum", "Bifidobacterium_breve"), kingdom = rep("Bacteria", 16), phylum = c("Actinobacteria", "Proteobacteria", "Actinobacteria", "Proteobacteria", "Proteobacteria", "Firmicutes", "Firmicutes", "Firmicutes", "Proteobacteria", "Proteobacteria", "Proteobacteria", "Proteobacteria", "Firmicutes", "Firmicutes", "Actinobacteria", "Actinobacteria"), class = c("Actinobacteria", "Gammaproteobacteria", "Actinobacteria", "Gammaproteobacteria", "Alphaproteobacteria", "Bacilli", "Bacilli", "Bacilli", "Gammaproteobacteria", "Gammaproteobacteria", "Gammaproteobacteria", "Gammaproteobacteria", "Bacilli", "Bacilli", "Actinobacteria", "Actinobacteria"), order = c("Propionibacteriales", "Enterobacterales", "Propionibacteriales", "Pseudomonadales", "Rhodospirillales", "Bacillales", "Bacillales", "Bacillales", "Enterobacterales", "Enterobacterales", "Pseudomonadales", "Pseudomonadales", "Lactobacillales", "Lactobacillales", "Bifidobacteriales", "Bifidobacteriales"), family = c("Propionibacteriaceae", "Enterobacteriaceae", "Propionibacteriaceae", "Moraxellaceae", "Rhodospirillaceae", "Staphylococcaceae", "Staphylococcaceae", "Staphylococcaceae", "Enterobacteriaceae", "Enterobacteriaceae", "Pseudomonadaceae", "Moraxellaceae", "Lactobacillaceae", "Lactobacillaceae", "Bifidobacteriaceae", "Bifidobacteriaceae"), genus = c("Cutibacterium", "Klebsiella", "Cutibacterium", "Moraxella", "Enhydrobacter", "Staphylococcus", "Staphylococcus", "Staphylococcus", "Escherichia", "Enterobacter", "Pseudomonas", "Acinetobacter", "Lactobacillus", "Lactobacillus", "Bifidobacterium", "Bifidobacterium"), stringsAsFactors = FALSE ) set.seed(42) fc <- setNames(rnorm(16), lineage_df$species) results <- taxon_rank_sets(fc, lineage_df, min_set_size = 2) names(results) results$family # --- Example 2: SummarizedExperiment / TreeSummarizedExperiment input --- library(mia) data(GlobalPatterns, package = "mia") tse <- GlobalPatterns # Run differential abundance (e.g., ALDEx2) to get fold changes # aldex_out <- ... (your DA analysis) # fc <- aldex_out$effect # names(fc) <- rownames(aldex_out) # Run taxon rank set enrichment directly from the TSE results <- taxon_rank_sets(fc, tse, min_set_size = 5) names(results) # Kingdom, Phylum, Class, Order, Family, Genus results$Family
Modular TaxSEA implementation supporting enrichment (KS) and
ORA (Fisher).
Provide either taxon_ranks for enrichment or
input_taxa for ORA.
TaxSEA( taxon_ranks = NULL, input_taxa = NULL, mode = NULL, lookup_missing = FALSE, min_set_size = 5, max_set_size = 300, custom_db = NULL, bugsigdb = FALSE )TaxSEA( taxon_ranks = NULL, input_taxa = NULL, mode = NULL, lookup_missing = FALSE, min_set_size = 5, max_set_size = 300, custom_db = NULL, bugsigdb = FALSE )
taxon_ranks |
Named numeric vector of statistics (e.g., log2 fold changes). Required for enrichment. |
input_taxa |
Character vector of taxa to treat as "hits"/selected taxa. Required for ORA. |
mode |
Character. One of |
lookup_missing |
Logical indicating whether to fetch missing NCBI IDs. Default is FALSE. |
min_set_size |
Minimum size of taxon sets to include in the analysis. Default is 5. |
max_set_size |
Maximum size of taxon sets to include in the analysis. Default is 300. |
custom_db |
A user-provided list of taxon sets. If NULL (default), the built-in database is used. |
bugsigdb |
Logical; whether to augment the built-in database with
BugSigDB signatures, which are downloaded (and cached) at run time
via |
A list of data frames with taxon set results.
data("TaxSEA_test_data") res <- TaxSEA(taxon_ranks = TaxSEA_test_data) head(res$All_databases) # Include BugSigDB signatures (downloaded at run time) res_bsdb <- TaxSEA(taxon_ranks = TaxSEA_test_data, bugsigdb = TRUE) head(res_bsdb$BugSigDB) # ORA example (toy): treat taxa with positive values as "hits" hits <- names(TaxSEA_test_data) res_ora <- TaxSEA(input_taxa = hits, mode = "ora") head(res_ora$All_databases)data("TaxSEA_test_data") res <- TaxSEA(taxon_ranks = TaxSEA_test_data) head(res$All_databases) # Include BugSigDB signatures (downloaded at run time) res_bsdb <- TaxSEA(taxon_ranks = TaxSEA_test_data, bugsigdb = TRUE) head(res_bsdb$BugSigDB) # ORA example (toy): treat taxa with positive values as "hits" hits <- names(TaxSEA_test_data) res_ora <- TaxSEA(input_taxa = hits, mode = "ora") head(res_ora$All_databases)
TaxSEA Database A dataset containing taxon sets. Each item in the list is a taxon set, and each member within a taxon set is a taxon.
data(TaxSEA_db)data(TaxSEA_db)
A list of vectors. Each vector contains character strings representing taxa as NCBI taxonomy IDs.
Set names start with their source: MiMeDB_, GutMGene_,
GMRepoV2_, mBodyMap_, BacDive_ and so on.
The BacDive_ sets describe measured phenotypes of cultured strains
(oxygen tolerance, Gram stain, substrate use, enzyme activities, growth
temperature, salt and pH ranges and others), aggregated to species over
every strain BacDive holds. They are built by the bacdive_harvester
pipeline (https://github.com/feargalr/bacdive_harvester); see
inst/scripts/make_BacDive_sets.R. Only sets with at least 3 members
are included. Names read BacDive_<Family>_<value>, optionally
followed by an assay context (e.g. _fermentation,
_reduction) and _negative for species tested and found
negative, e.g. BacDive_Oxygen_facultative_anaerobe,
BacDive_Uses_nitrate_reduction.
Per-set provenance (source family, number of species and strains behind each set) and a record of the build (harvester release, BacDive client version, harvest dates) are shipped with the package:
system.file("extdata", "BacDive_set_provenance.tsv", package = "TaxSEA")
system.file("extdata", "BacDive_build_info.tsv", package = "TaxSEA")
See the README. BacDive sets: Schober et al. (2025) BacDive in 2025: the core database for prokaryotic strain data, via https://bacdive.dsmz.de.
data(TaxSEA_db) all_sets <- names(TaxSEA_db) GABA_producers<-TaxSEA_db[["MiMeDB_producers_of_GABA"]]data(TaxSEA_db) all_sets <- names(TaxSEA_db) GABA_producers<-TaxSEA_db[["MiMeDB_producers_of_GABA"]]
A dataset containing taxon ranks and taxon IDs.
data(TaxSEA_test_data)data(TaxSEA_test_data)
A data frame with two columns:
Character vector representing taxon ranks
Character vector representing taxon IDs
A data frame with columns 'rank' and 'id' representing taxon ranks and taxon IDs, respectively.
See READ ME.
data(TaxSEA_test_data) test_results <- TaxSEA(TaxSEA_test_data)data(TaxSEA_test_data) test_results <- TaxSEA(TaxSEA_test_data)