TaxSEA: Taxon Set Enrichment Analysis

📦 Available on Bioconductor
📝 Read the paper in Briefings in Bioinformatics

Note: We found a few minor bugs in the current release. These mostly relate to how the NCBI ID API is contacted. These have been fixed in the development version of TaxSEA. In this latest version there is also a more dynamic integration with BugSigDB. To install this version you will need to have R 4.5 installed and the most recent Bioconductor release. Then run: r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("TaxSEA", version = "devel")

TaxSEA helps microbiome researchers test for enrichment in known microbial signatures, including:

• Metabolite producers
  
• Disease associations
  
• Previously published microbiome signatures

TaxSEA takes as input a vector of genus or species names and a rank. For example log2 fold changes or Spearman's rho. TaxSEA then uses a Kolmogorov-Smirnov test to identify if a particular group of species or genera (i.e. a set of taxa such as butyrate producers) are skewed to one end of the distribution .

Note: Although TaxSEA in principle can be applied to microbiome data from any source, the databases utilized largely cover human associated microbiomes and the human gut microbiome in particular. As such TaxSEA will likely perform best on human gut microbiome data.

Taxon set database

By default TaxSEA utilizes taxon sets generated from five reference databases (gutMGene, GMrepo v2, MiMeDB, mBodyMap, BugSigDB). See below for examples of using custom databases or taxonomically defined taxon sets.

Please cite the appropriate database if using:

• Cheng et al. gutMGene: a comprehensive database for target genes of gut microbes and microbial metabolites Nucleic Acids Res. 2022.
• Dai et al. GMrepo v2: a curated human gut microbiome database with special focus on disease markers and cross-dataset comparison Nucleic Acids Res. 2022.
• Wishart et. al. MiMeDB: the Human Microbial Metabolome Database Nucleic Acids Res. 2023.
• Jin et al. mBodyMap: a curated database for microbes across human body and their associations with health and diseases. Nucleic Acids Res. 2022.
• Geistlinger et al. BugSigDB captures patterns of differential abundance across a broad range of host-associated microbial signatures. Nature Biotech. 2023.

Installation

```r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("TaxSEA")

```

Usage

Quick start

```r library(TaxSEA)

Retrieve taxon sets containing Bifidobacterium longum.

blong.sets <- gettaxonsets(taxon="Bifidobacterium_longum")

Run TaxSEA with test data provided

data(TaxSEAtestdata) taxsearesults <- TaxSEA(taxonranks=TaxSEAtestdata)

Enrichments among metabolite producers from gutMgene and MiMeDB

metabolites.df <- taxsearesults$Metaboliteproducers

Enrichments among health and disease signatures from GMRepoV2 and mBodyMap

disease.df <- taxsearesults$Healthassociations

Enrichments amongh published associations from BugSigDB

bsdb.df <- taxsea_results$BugSigdB

```

Input

All that is required for TaxSEA is a vector in R containing ranks (e.g. log2 fold changes) and names (E.g. species/genus). TaxSEA will not work for ranks higher than species or genus. The input should be for all taxa tested, and not limited to only a pre-defined set (e.g. do not use a threshold for significance or remove any taxa). See example below for format. TaxSEA will lookup and convert taxon names to NCBI taxonomic identifiers. TaxSEA stores a commonly observed identifiers internally and so will only look up whatever is not covered to save time.

Input IDs should be in the format of like one of the following

• Species name. E.g. "Bifidobacterium longum", "Bifidobacterium_longum"
• Genus name. E.g. "Bifidobacterium"
• NCBI ID E.g. 216816

```r

Input IDs with the full taxonomic lineage should be split up. E.g.

x <- "dBacteria.pActinobacteriota.cActinomycetes.oBifidobacteriales.fBifidobacteriaceae.gBifidobacterium" x <- strsplit(x,split="\.")[[1]][6] x <- gsub("g__","",x)

Running this through a vector of IDs may look something like the following

newids <- sapply(as.character(oldids),function(y) {strsplit(x = y,split="\.")[[1]][6]})

newids <- gsub("g","",newids)

Example test data

library(TaxSEA) data("TaxSEAtestdata") head(sample(TaxSEAtestdata),4)

```

Run TaxSEA with test data

```r data("TaxSEAtestdata") taxsearesults <- TaxSEA(taxonranks=TaxSEAtestdata)

Enrichments among metabolite producers from gutMgene and MiMeDB

metabolites.df <- taxsearesults$Metaboliteproducers

Enrichments among health and disease signatures from GMRepoV2 and mBodyMap

disease.df <- taxsearesults$Healthassociations

Enrichments among published associations from BugSigDB

bsdb.df <- taxsea_results$BugSigDB

```

Test data

The test data provided with TaxSEA consists of log2 fold changes comparing between healthy and IBD. The count data for this was downloaded from curatedMetagenomeData and fold changes generated with LinDA.

• Hall et al. A novel Ruminococcus gnavus clade enriched in inflammatory bowel disease patients** Genome Med. 2017 Nov 28;9(1):103.
• Pasolli et al. Accessible, curated metagenomic data through ExperimentHub. Nat Methods. 2017 Oct 31;14(11):1023-1024. doi: 10.1038/nmeth.4468.
• Zhou et al. LinDA: linear models for differential abundance analysis of microbiome compositional data. Genome Biol. 2022 Apr 14;23(1):95.

```r

head(sample(TaxSEAtestdata),3) Bacteroidesthetaiotaomicron BlautiaspCAG257 Ruminococcus bromii 1.908 3.650 -5.038 ```

Output

The output is a list of three data frames providing enrichment results for metabolite produers, health/disease associations, and published signatures from BugSigDB. Each dataframe has 5 columns

• taxonSetName - The name of the taxon set tested
• median_rank - This is simply the median rank across all detected members in the set. This allows you to see the direction of change
• P value - Kolmogorov-Smirnov test P value.
• FDR - P value adjusted for multiple testing.
• TaxonSet - Returns list of taxa in the set to show what is driving the signal

Custom databases

Many users may want to utilise TaxSEA with a custom database. For example for testing if there is a flag in the TaxSEA function "customdb" which expects as input a named list of vectors. This is the same format as the default TaxSEA database. Note: using the customdb flag disables the automatic ID conversion and NCBI API lookup. However we have functionality available via other functions

```r

Perform enrichment analysis using TaxSEA

customtaxsearesults <- TaxSEA(taxonranks = log2foldchanges, customdb = customtaxonsets) customtaxsearesults <- customtaxsearesults$custom_sets ```

Testing for differences in taxonomically defined sets

In addition to taxon sets defined by function or phenotype, users can define sets based on taxonomy. Current methods to test at higher taxonomic levels (e.g., genus or family), involve aggregating counts but with this approach opposing shifts in individual species may cancel each other out, obscuring meaningful biological patterns. For instance, antibiotic treatment may suppress certain species while allowing resistant species within the same genus to expand and occupy the vacant niche, creating an ecological shift that appears as no net change at broader taxonomic levels. Here we utilise data from Chng et al. demonstrating this in a comparsion between Atopic dermatitis and controls.

```r

Applying TaxSEA functionality to taxonomic ranks

This script applies TaxSEA to identify taxonomic enrichment at different taxonomic levels.

Specifically, we analyze enrichment at the family level using metagenomic data

Load required libraries

library(TaxSEA) library(curatedMetagenomicData) library(tidyverse) library(phyloseq) library(MicrobiomeStat) library(dplyr)

Load sample metadata

metadata_all <- sampleMetadata

Filter metadata for the specific study (ChngKR_2016)

metadata <- metadataall %>% filter(studyname == "ChngKR2016") %>% columntorownames('sampleid')

Extract count data using curatedMetagenomicData

cmddata <- curatedMetagenomicData( pattern = "ChngKR2016.relative_abundance", counts = TRUE, dryrun = FALSE )

Convert the extracted data to a count matrix

countsdata <- assay(cmddata[[1]]) countsdata <- countsdata[, rownames(metadata)] # Subset to relevant samples

Filter taxa with at least one sample having counts > 100

countsdata <- countsdata[apply(counts_data > 100, 1, sum) > 0, ]

Extract species names from taxonomic strings

speciesnames <- gsub("s", "", sapply(rownames(countsdata), function(y) strsplit(y, "\|")[[1]][7])) rownames(countsdata) <- speciesnames

Create a taxonomic lineage dataframe

Remove taxonomic prefixes (k, p, c__, etc.) and separate into taxonomic ranks

make data frame of taxon lineages

taxonlineages <- data.frame(Name = speciesnames, Lineage = names(speciesnames)) %>% mutate(Lineage = strremoveall(Lineage, '[kpcofgs]')) %>% separate(col = Lineage, into = c('kingdom', 'phylum', 'class', 'order', 'family', 'genus', 'species'), sep = '\|') %>% mutate(name = Name) %>% removerownames() %>% columntorownames('name')

Perform differential abundance testing using LinDA

metadata$studycondition <- factor(metadata$studycondition, levels = c("control", "AD"))

lindaresults <- linda( feature.dat = countsdata, meta.dat = metadata, formula = '~study_condition', feature.dat.type = 'count', prev.filter = 0.05 )

Extract log2 fold change values for differential taxa

lindaresults <- lindaresults$output$studyconditionAD log2foldchanges <- lindaresults$log2FoldChange names(log2foldchanges) <- rownames(linda_results)

Define the taxonomic rank for enrichment analysis

selectedtaxonlevel <- 'genus' # Modify as needed (e.g., genus, phylum)

Create a named list of species grouped by taxonomic rank

customtaxonsets <- taxonlineages %>% groupby(.data[[selectedtaxonlevel]]) %>% summarise(species = list(species), .groups = "drop") %>% deframe()

Perform enrichment analysis using TaxSEA

customtaxsearesults <- TaxSEA(taxonranks = log2foldchanges, customdb = customtaxonsets) customtaxsearesults <- customtaxsearesults$custom_sets

```

Visualisation of TaxSEA output.

The results above were generated by using TaxSEA on the output of a differential abundance analysis comparing between disease and control. The input was per spcies log2 fold changes between taxa in Inflammatory Bowel disease and control. TaxSEA identified a significant depletion in the producers of certain short chain fatty acids. Using barplots we can show the overall signatures identified as significantly different. We can then highlight the individual species contributing to this signature on a volcano plot.

BugSigDB

The format of BugSigDB is that each publication is entered as a "Study", and within this there is different experiments and signatures. For example one of the signatures may be taxa increased in an experiment, and another signature is taxa that are decreased. Users can find out more by querying the BugSigDB. In recent updates BugSigDB has moved to using the PubMed ID as the study ID (Although not all studies have these yet). See below for an example.

```r library(bugsigdbr) #This package is installable via Bioconductor bsdb <- importBugSigDB() #Import database

E.g. if the BugSigDB identifier you found enriched was bsdb:11/1/1

This is Study 11, Experiment 1, Signature 1

bsdb[bsdb$Study=="Study 11" & bsdb$Experiment=="Experiment 1" & bsdb$Signature=="Signature 1",]

```

TaxSEA database with other enrichment tools

The TaxSEA function by default uses the Kolmogorov Smirnov test and the original idea was inspired by gene set enrichment analyses from RNASeq. Should users wish to use an alternative gene set enrichment analysis tool the database is formatted in such a way that should be possible. See below for an example with fast gene set enrichment analysis (fgsea).

```r library(fgsea) #This package is installable via Bioconductor data(TaxSEA_DB)

Convert input names to NCBI taxon ids

names(TaxSEAtestdata) = getncbitaxonids(names(TaxSEAtestdata)) TaxSEAtestdata = TaxSEAtestdata[!is.na(names(TaxSEAtest_data))]

Run fgsea

fgsearesults <- fgsea(TaxSEAdb, TaxSEAtestdata, minSize=5, maxSize=500) ```