A common way to test for differences at higher taxonomic levels (e.g. genus, family, order) is to aggregate counts and then run differential abundance on the collapsed table. That approach is simple, but it has a big flaw: opposing shifts cancel each other out. If one species in a genus goes up and another goes down, aggregation can make it look like nothing changed, even though the biology clearly did.
TaxSEA offers an alternative.
Instead of collapsing your data, you: - Run your usual analysis at the species level - Keep the full ranked list (e.g. log2FC, correlations, model coefficients) - Define taxonomic sets (e.g. all species belonging to the same family) - Test whether those sets are coordinately shifted using enrichment
In other words, you test higher-level taxonomy without throwing away species-level resolution.
This lets you ask questions like: - Are members of this family consistently depleted, even if no single species is significant? - Is this genus showing a directional shift driven by many small effects? - Are higher-level patterns being masked by species-level heterogeneity?
This is particularly powerful in: - antibiotic or perturbation experiments - situations where ecological turnover happens within taxonomic groups
This functionality is available in TaxSEA in the taxon_rank_sets() function.
#### Enrichment testing with taxonomically defined ranks
# --- 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