Main

Despite substantial progress in experimental techniques and computational methods for culture-independent profiling of the human microbiome, the analysis and interpretation of microbial differential abundance studies remains challenging. A large body of experimental and observational studies on humans and in animal models has reported associations between host-associated microbiomes and the onset, progression and treatment of a variety of diseases, including atherosclerosis1, cardiovascular diseases2, cancers3 and diabetes4. This growing body of published results provides opportunities for synthesis of accumulated knowledge, identification of common patterns across different diseases and exposures and interpretation of new studies by comparison to previous results. However, without a systematic catalog of published differential abundance results, even identical microbial signatures reported in different research fields are unlikely to be noticed. Even within research fields, systematic review in the absence of a catalog or common reporting of differential abundance results are time consuming, static and generally do not summarize all taxa reported.

This situation has parallels to early challenges in the interpretation of differential gene expression analyses5, which have been addressed in the field by gene set enrichment analysis (GSEA). GSEA allows for the comparison of coherent expression patterns among predefined gene signatures that share a biological function or property or that were identified together by a previous study6,7. GSEA is a key tool in gene expression data analysis6, with a wide range of subsequent methods8 to account for correlations between genes9, redundancy of functional annotation10, different types of null hypotheses11 and the application of GSEA for the analysis of genomic regions12, metabolomic data13 and disease phenotypes14.

Analogously, differential microbial abundance analysis can yield lists or ‘signatures’ of microbial clades at multiple taxonomic levels that are associated with a phenotype of interest. The properties shared by these clades are often not obvious but could include common environmental exposures, metabolic or ecological requirements or physiological characteristics. Although nascent attempts to apply concepts of GSEA to the results of microbiome differential abundance analyses exist15,16,17,18, major obstacles have prevented their broad utility and adoption. The most notable obstacle has been the lack of comprehensive databases of signatures designed for enrichment analysis, such as those available for GSEA, including the Gene Ontology19, KEGG20, MSigDB21,22 and GeneSigDB23. Several databases provide important information on microbial physiology and morphology24,25,26,27 but are not designed for enrichment analysis and, by design, exclude the vast majority of experimentally derived microbial signatures associated with cancer, inflammation, diet or other conditions studied in human and other host-associated microbiome research.

This study provides two main contributions to enable high-throughput comparison of published microbial signatures. First, it describes BugSigDB, a database of published microbial signatures of sufficient scale and diversity to capture replicable patterns of differential abundance across a broad spectrum of the host-associated microbiome literature. BugSigDB provides curated published signatures of differentially abundant microbes associated with a wide range of health outcomes, pharmaceutical usage (for example, antibiotics), experiments on animal models, randomized clinical trials and microbial attributes and is built on the technology of Wikipedia to allow community contributions, revisions and review for quality control. Second, we provide a systematic analysis of the results reported by hundreds of published microbiome studies and identify replicated patterns even across 16S amplicon and shotgun sequencing approaches, demonstrating that interpretation of new microbiome studies can be supported by systematic comparison to previously published signatures. Database-wide analysis revealed common patterns of microbe co-occurrence and mutual exclusivity within signatures and identified antibiotic treatment and human immunodeficiency virus (HIV) infection as the experimental conditions with the highest level of consistency in signatures reported by independent studies. ‘Bug set’ enrichment analysis of ten individual colorectal cancer-associated fecal microbiome datasets (N = 663) detected published meta-analysis signatures used as positive controls, supporting direct application of methods adapted from gene set analysis and taxonomy-aware enrichment methods. This analysis also identified commonalities among signatures of other diseases, including elevated frequency of oral pathobionts in the gut, and identified the most common patterns of co-occurrence and mutual exclusion across all conditions or related conditions or unique to one condition. Together, the BugSigDB database and analysis methods described here improve the interpretation of new microbiome studies by systematic comparison to published microbial signatures.

Results

A curated database of published microbial signatures

BugSigDB comprises a comprehensive database of manually curated, host-associated microbial signatures from published microbiome studies of human, mouse and rat (Fig. 1). The database has been simultaneously developed and expanded over the course of more than 4 years with contributions from more than 25 curators trained in-house and currently contains >2,500 microbial signatures extracted from >600 scientific articles (Fig. 1a). The curated papers cover two decades of microbiome research, with the majority of studies being published in the last 5 years (459 of 628 articles, 73.1%; Supplementary Fig. 1). Among these papers are microbiome studies of participants from more than 50 different countries, with more than 50% of the studies originating from China and the United States (201 and 157 studies, respectively; Fig. 1b).

Fig. 1: BugSigDB, a curated database of experimentally derived microbial signatures.
figure 1

ac, BugSigDB is a community-editable collection of published microbiome studies reporting differentially abundant host-associated microbiota (including bacteria, archaea, fungi, protists and viruses) that currently includes >2,500 microbial signatures extracted from >600 scientific articles over the course of more than 4 years (a). These papers report microbiome studies of participants representing different countries and ethnogeographies (b) and microbiome samples from different human body sites (c). d, Number of signatures associated with a specific microbe, with the top 10 most frequently reported microbes annotated. BugSigDB signatures contain taxonomic levels from phylum to strain standardized based on the NCBI Taxonomy29. e, Signature sizes, with more than 50% of the signatures containing five or more microbes. f, Percentage of signatures annotated to major disease categories when classifying the study condition associated with each signature according to the Experimental Factor Ontology32. g, Top 10 most frequently reported microbes and the number of associated signatures stratified by disease category.

Studying microbiome samples from 14 broad body areas comprising more than 60 refined anatomical sites according to the UBERON Anatomy Ontology28 (Supplementary Table 1), the majority of studies in BugSigDB analyzed gut (440 of 628 studies, 70.1%), oral (80 of 628 studies, 12.8%) and vaginal microbiome samples (59 of 628 studies, 9.4%; Fig. 1c). The signatures were generated by both 16S amplicon sequencing (92.5%) and metagenomic shotgun sequencing (MGX; 7.5%; Table 1) and contained taxonomic levels from phylum to strain standardized based on the NCBI Taxonomy29 (Supplementary Fig. 2). BugSigDB is implemented as a semantic MediaWiki30 web interface available at https://bugsigdb.org that supports data entry, semantic validation, quality control and web-based programmatic access to annotations for studies, experiments, signatures and individual taxa (Supplementary Fig. 3 and “Data entry, validation and access” in the Methods).

Table 1 Overview of curated metadata annotations in BugSigDB

Most of the roughly 1,400 unique taxa contained in BugSigDB are reported as differentially abundant in fewer than five signatures (1,009 of 1,370 unique microbes, 73.6%; Fig. 1d). Conversely, genera Streptococcus, Prevotella, Bacteroides and Lactobacillus are each reported as differentially abundant in more than 100 signatures, reflecting the large number of species belonging to these genera and paralleling previous observations of gene signatures with certain genes often being disproportionately reported as differentially expressed31. Signatures contain six microbes, on average, with roughly 55% of the signatures containing five or more microbes (Fig. 1e).

Study conditions associated with each signature are classified according to the Experimental Factor Ontology32 (Fig. 1f), with large proportions associated with different types of cancer (23.1%), anatomical system diseases, such as asthma and endometriosis (22.6%), and metabolic diseases, such as obesity and diabetes (15.3%). Apart from different disease categories, substantial proportions of signatures in BugSigDB are also associated with diet (6%), the use of antibiotics (4.8%) and birth delivery mode (vaginal versus Cesarean section; 4.1%). These distributions approximately reflect the human microbiome literature. Non-human host species remain underrepresented, although the database can support any host organism.

Condition-specific associations were investigated for the ten most-reported taxa (Fig. 1g). For instance, signatures containing the Bacteroides genus are related to metabolic disease (21%) more frequently than all signatures (13%; P = 0.003, χ2 = 8.4, d.f. = 1, two-sided proportion test). Bacteroides is similarly enriched in signatures of antibiotic exposure (11% of Bacteroides-containing signatures being associated with antibiotics exposure versus 5% of all signatures; P = 0.0009, χ2 = 11.2, d.f. = 1). BugSigDB makes such commonalities across groups of related studies straightforward to identify and provides for continuous updates of these associations as the database grows.

Curated metadata and common practices in microbiome research

BugSigDB provides curated metadata at the level of study (study design and automatically generated citation information), experiments within a study that each define one contrast for differential abundance analysis (such as characteristics about study subjects, lab analysis, statistical analysis and α-diversity) and up to two signatures within an experiment, each of which contains one or more taxon (see Table 1 for summary statistics of studies and experiments and “Definition of semantic concepts” in the Methods for details of the design). In this article, we use ‘study group’ to refer to cases in case–control studies, the exposed group in exposure–control studies and whichever group corresponds to the condition of interest with increased microbial relative abundance in other comparisons of two-sample groups.

Signatures are available primarily from observational study designs: case–control (281 studies, 44.7%) and cross-sectional studies (171 studies, 27.2%) were most prevalent, while prospective cohort studies (66, 10.5%), time series/longitudinal studies (51, 8.1%), laboratory studies (29, 4.6%), randomized controlled trials (28, 4.4%) and meta-analyses (6, 0.9%) were also present. Study subject information includes host species (95.9% human of 1,223 experiments), location, condition, body site (Fig. 1), antibiotics exclusion criteria (median exclusion time = 60 d) and sample size in study and control sample groups (median total sample size = 25).

A survey of statistical methods most frequently applied for differential abundance testing in BugSigDB revealed that nonparametric tests, such as Wilcoxon family tests and the related LEfSe software33, were most frequently used, whereas recently suggested tools for differential abundance tests accounting for the compositionality of microbiome data34 were rarely used (Supplementary Results Section 1.1). Furthermore, when stratifying experiments by body site and condition, oral and vaginal samples were frequently reported with increased α-diversity in the study group, as opposed to samples from the gastrointestinal tract, which were frequently found with decreased α-diversity in the study group (Supplementary Tables 24).

Conditions with replicable microbiome changes across studies

BugSigDB facilitates meta-analysis of differential abundance studies and enables the identification of experimental conditions and disease phenotypes where microbiome changes replicate across studies. Focusing on 1,194 signatures derived from human fecal samples in 311 published studies, we computed signature similarity within conditions and assessed whether the resulting similarity exceeds the similarity of randomly sampled signatures (Fig. 2). This simultaneously determines whether signatures of the same phenotype reproduce across studies and whether different phenotypes share similar microbial signatures.

Fig. 2: Analysis of signature similarity in BugSigDB identifies conditions with replicable microbiome changes across studies.
figure 2

a, Signature similarity analysis for conditions from at least two studies for the same body site. Signatures were stratified by direction of abundance change into signatures with decreased abundance (left) and increased abundance (right) in the study group. Conditions are divided into two groups: conditions with highest reproducibility for signatures of decreased abundance (top) and increased abundance (bottom) in the study condition. Shown is the nominal P value (x axis, negative log scale) obtained from testing whether the semantic similarity of signatures for each condition (y axis) exceeds the semantic similarity of randomly sampled signatures (one-sided resampling test). The size of each dot corresponds to the number of studies investigating a condition; HPV, human papillomavirus. be, Example exploration of two top-ranked conditions in more detail: fecal signatures of decreased abundance for individuals infected with HIV (b and c) and fecal signatures of decreased abundance for individuals treated with antibiotics (d and e). Semantic similarity between signatures is shown in b and d. Each node corresponds to a signature. The size of each node is proportional to the number of taxa in a signature. More similar signatures are connected by shorter and thicker edges. Microbial contents of the signatures shown in b and d (y axis) delineating the genera contained in these signatures (x axis) are shown in c and e. Signatures are clustered by semantic similarity, and the relative frequency of each genus across signatures is indicated by a color scale. The display of contained taxa is restricted to genera occurring in at least two signatures (full signatures are shown in Supplementary Figs. 5 and 6). Additional work cited includes refs. 37,38,39,40,41,42,43,44,45,46,47,48,49.

To do this, we applied two alternative approaches for computing similarity between signatures: (1) the more restrictive Jaccard index35 based on pairwise overlaps between signatures harmonized to genus level and (2) the more sensitive semantic similarity36 based on taxonomic distance between signatures of mixed taxonomic levels (“Signature similarity” in the Methods). Hierarchical clustering of signature similarity for both similarity measures was in good agreement, reflecting the dominance of genera reported so far in BugSigDB, but demonstrated better resolution using semantic similarity than the sparse results obtained from the application of Jaccard similarity (Supplementary Fig. 4). The advantages of semantic similarity may grow as taxonomic ranks in BugSigDB become more mixed due to increased reporting of species-level results made possible by shotgun metagenomics.

To assess replication by independent studies of the same condition, we compared semantic similarity between signatures reported for a single condition to the similarity of randomly sampled signatures in repeated simulation and ranked conditions based on the resulting empirical P value (Fig. 2a). Differential abundance signatures of antibiotic treatment and HIV infection were among the most consistent, each investigated by five or more studies in BugSigDB37,38,39,40,41,42,43,44,45,46,47,48,49. Closer inspection of signatures of decreased abundance after antibiotics treatment (semantic similarity = 0.64, one-sided resampling P = 0.0005; Fig. 2d,e and Supplementary Fig. 5) revealed enrichment for genera of fastidious anaerobes that are often short-chain fatty acid producers and displayed frequent loss of Bifidobacterium and Blautia, in agreement with previous reports50. However, signatures of decreased abundance in individuals infected with HIV versus those in uninfected healthy individuals (semantic similarity = 0.68, resampling P = 0.002; Fig. 2b,c and Supplementary Fig. 6) displayed a loss of abundant members of healthy gut microbial communities, as typically observed for diseases associated with gastrointestinal inflammation51, but also resembled a response to antibiotics treatment, a likely side effect of antibiotics often being prescribed for individuals positive for HIV infection to prevent or treat opportunistic and associated infections52. Additional examples of replicable microbiome shifts between studies included similarity (1) among fecal signatures of decreased abundance in individuals with coronavirus disease 2019 (COVID-19), which displayed alterations associated with intensive care units and antibiotics treatment53 (Supplementary Fig. 7), and (2) among stomach signatures of individuals with gastric cancer driven by consistently increased abundance of Streptococcus, Lactobacillus and Prevotella (Supplementary Fig. 8).

Bug set enrichment analysis of colorectal cancer signatures

We integrated BugSigDB signatures with the manually curated metagenomic datasets from curatedMetagenomicData54 to systematically benchmark enrichment methods from the EnrichmentBrowser package7,55 and evaluate whether top-performing gene set enrichment methods can be directly applied to microbiome data. We applied two enrichment methods that have performed well in previous benchmarking of gene set analysis methods7 (overrepresentation analysis (ORA11) and pathway analysis with downweighting of overlapping genes (PADOG31)) to ten colorectal cancer datasets from curatedMetagenomicData. We performed enrichment analysis of all microbiome signatures from BugSigDB simultaneously and used as ‘spike-in’ controls two signatures of the colorectal cancer-associated fecal microbiome derived previously by two independent meta-analyses of individual participant data56,57 from eight of the studies included in this dataset (Fig. 3). The two signatures based on meta-analysis are thus expected to be among the most enriched of all microbiome signatures in the eight training studies for the meta-analyses due to their shared utilization in both datasets and in the remaining two independent studies due to the relevance and large sample size of the colorectal cancer signatures. The two meta-analysis signatures can thus also be considered robust against spurious signals from studies with small sample sizes that were not included in the analysis (Supplementary Results Section 3.2 and Supplementary Tables 5 and 6).

Fig. 3: Bug set enrichment analysis reveals associations of the human microbiome with colorectal cancer and other disease phenotypes.
figure 3

a, ORA of BugSigDB signatures with differentially abundant genera between 662 colorectal cancer samples and 653 control samples pooled from ten published metagenomic datasets. Each dot corresponds to an enriched published signature from BugSigDB (false discovery rate P value (Padj) of <0.05, one-sided Fisher’s exact test, color scale). The size of each dot corresponds to the number of differentially abundant genera in a signature, given as a proportion on the x axis. Asterisks (*) indicate prior meta-analytic signatures from Thomas et al.56 and Wirbel et al.57, which reported differentially abundant species and genera from meta-analysis of eight colorectal cancer datasets. These act as positive controls, as several of the datasets being tested for enrichment were included at the time these signatures were derived. b, Differentially abundant genera (x axis) in the top 15 enriched signatures (y axis) from the overrepresentation analysis in a. The x axis is sorted by occurrence frequency of each genus in decreasing order, and the y axis is sorted by Jaccard similarity between signatures. c, Percentile of ranks (y axis) of both meta-analytic signatures for ORA, PADOG and CBEA relative to all signatures when applied to the ten published metagenomic datasets individually (x axis). Number symbols (#) indicate datasets not included as training sets in the meta-analyses of Thomas et al.56 and Wirbel et al.57. d, Assessment of statistically significant differences of the percentile ranks on the n = 10 independent datasets shown in c between the three methods using a two-sided Wilcoxon’s signed-rank test. Box plots show median (vertical line), interquartile range (box) and ±1.5× interquartile range (whiskers). Additional work cited includes (refs. 97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113).

ORA of 647 BugSigDB signatures yielded 19 signatures enriched in contrasts using 662 colorectal cancer samples and 653 control samples from ten datasets (false discovery rate P value of <0.05, one-sided Fisher’s exact test; Fig. 3a). Only signatures containing five or more genera associated with any condition were included. The two positive-control spike-in signatures from Thomas et al.56 and Wirbel et al.57 were top ranked, as expected. Other enriched signatures included colorectal cancer signatures from Wu et al.58 and Allali et al.59. These are notable because neither study was included in the ten datasets from which the meta-analysis signatures were computed. Additionally, both were based on 16S amplicon sequencing, whereas the meta-analysis signatures were based on shotgun sequencing. This analysis therefore also provides a proof of concept for integrating species-level signatures from shotgun metagenomic data with genus-level signatures derived from 16S amplicon profiles (“Bug set enrichment analysis” in the Methods) and also affords independent replication of the signatures from Thomas et al.56 and Wirbel et al.57.

The presence of 11 enriched signatures (58%) from saliva samples from studies of oral diseases, such as gingivitis, peptic esophagitis and oral carcinoma, is consistent with recent reports that oral-to-gut microbial introgression is a feature of colorectal cancer60 and that periodontal diseases are a well-established risk factor for colorectal cancer61,62. Frequently overlapping genera between the enriched signatures include Fusobacterium, Porphyromonas and Peptostreptococcus, and all displayed strongly increased abundance in individuals with colorectal cancer relative to healthy controls (Fig. 3b).

Although these findings demonstrate the usefulness of ORA as a fast and effective enrichment method for microbiome signatures, the method has known shortcomings in the presence of correlated features63 or an inappropriately large feature universe64. The PADOG method, a top performer in several independent assessments7,65,66, is theoretically superior, as it applies sample permutation to preserve correlations and, by working on the full abundance matrix, does not require thresholding on differential abundance or the definition of a feature universe. In addition, the method downweights frequently overlapping microbes between signatures (such as those displayed in Fig. 1d), leading to increased sensitivity and the identification of more specific signatures for the phenotype under investigation.

We therefore benchmarked PADOG against ORA in each of the ten individual participant shotgun metagenomics datasets of colorectal cancer and compared the rankings of the two spike-in signatures from the meta-analyses of Thomas et al.56 and Wirbel et al.57 (Fig. 3c). On average, PADOG ranked the spike-in signatures better than ORA, although the difference was statistically significant only for the spike-in signature from Thomas et al.56 (P = 0.049, two-sided Wilcoxon signed-rank test; Fig. 3d). This difference was largest for datasets with smaller sample sizes, where the lack in power was more detrimental for ORA than for PADOG (Supplementary Fig. 9).

Despite the apparent effectiveness of established gene set enrichment methods for application to microbiome data, these methods were not developed with microbiome data in mind34,67. Competitive balances for taxonomic enrichment analysis (CBEA18) is a recent taxonomic enrichment method specifically developed for microbiome data that accounts for compositionality via application of an isometric log ratio transformation of relative abundance data for the computation of sample-level enrichment scores. Benchmarked against ORA and PADOG in the colorectal cancer setting (Fig. 3c), CBEA tended to rank the spike-in signature from Thomas et al.56 significantly higher than ORA (P = 0.021, Wilcoxon signed-rank test) but did not display a notable performance gain over PADOG (P = 0.92; Fig. 3d).

Microbe co-occurrence and mutual exclusivity within signatures

BugSigDB enables the exploration of compositional patterns within signatures of differential abundance for different body sites. Focusing on 1,194 signatures of fecal microbiomes from 311 published studies, we analyzed patterns of co-occurrence and mutual exclusivity for individual microbes and groups of microbes (Fig. 4). Inspection of the top 20 genera most frequently reported as differentially abundant in signatures from fecal samples (Fig. 4a) revealed genera predominantly belonging to the phyla Firmicutes (13 genera) and Bacteroidetes (4 genera), in agreement with those being the dominant phyla of the human gut microbiome68,69. Among the 20 were Bacteroides, Prevotella and Ruminococcus, three dominant gut genera that are highly variable in relative abundance70 and have a large effect on gut microbiome clustering71.

Fig. 4: Microbe co-occurrence and mutual exclusivity in BugSigDB reveals recurrent groups of taxa within signatures of differential abundance.
figure 4

a, Microbe–microbe co-occurrence and mutual exclusivity across 1,194 signatures of fecal microbiomes from 311 published studies. ‘Abundance in group 1’ shows the top 20 genera most frequently reported as differentially abundant in the study group of these signatures. Stars indicate microbes that tend to be predominantly unidirectional, that is, reported specifically with increased or decreased abundance in the study group (based on a greater than twofold difference in proportions). ‘Prevalence’ shows the prevalence of these genera in 9,623 healthy adult stool samples pooled from 68 different studies. The top heat map shows co-occurrence between these top 20 genera, that is, in how many signatures these microbes are reported together as differentially abundant with the same direction of abundance change (that is, occurring together in either the signature of increased or decreased abundance). The bottom heat map shows mutual exclusivity between these top 20 genera, that is, in how many signatures these microbes are reported as differentially abundant with opposite direction of abundance change (that is, one microbe in the signature of increased abundance and the other in the signature of decreased abundance or vice versa). b, Spearman correlation between the top 20 genera in healthy adult stool samples as indicated by circles of varying size in a. c, Proportion of signatures where each of the top 20 genera was reported with increased abundance in the study group (y axis) against the prevalence of these genera (x axis) in healthy adult stool samples as shown in a. The correlation coefficient and P value of a two-sided Pearson correlation test are shown at the top. The error bands indicate the 95% confidence interval of the linear regression fit.

When comparing the proportions of signatures where these genera were reported either with increased or decreased abundance in the study sample group (Fig. 4a, Abundance in group), the genera most disproportionately increased in the study group were Enterococcus (64 of 73 signatures, 87.7%), Lactobacillus (73 of 105, 69.5%), Veillonella (50 of 74, 67.6%) and Streptococcus (83 of 128, 64.8%). Pathogenicity of certain Streptococcus and Enterococcus species is well documented72,73, whereas Lactobacillus and Veillonella are typically considered to be commensal gut microbiota of limited pathogenicity74,75. Genera decreased in the study group were Roseburia (77 of 103, 74.8%), a beneficial gut organism with established anti-inflammatory activity76, and Alistipes (61 of 88, 69.3%), for which protective and harmful associations with a range of diseases have been described77.

To elucidate the extent to which frequency of differential abundance in disease phenotypes can be explained by prevalence in the healthy gut microbiome, we contrasted these findings with the prevalence of these genera in 9,623 stool samples from healthy adult participants from 68 different studies in curatedMetagenomicData (Fig. 4a, Prevalence). We observed a strong negative correlation between prevalence, measured as the percentage of control samples in which the genus is observed at non-zero relative abundance, and the proportion of BugSigDB signatures in which the genus is reported with increased abundance in disease (r = −0.84, P = 3 × 10−6, two-sided Pearson correlation test; Fig. 4c). This indicates that across the many diseases and exposures present in BugSigDB, high-prevalence genera tend to be lost in study groups relative to controls, and low-prevalence genera have increased abundance in the study group. This was particularly apparent for Enterococcus, a genus of low prevalence in healthy samples (13%) that was reported almost exclusively with increased abundance in the diseased group (64 of 73 signatures, 88%). The presence of Enterococcus, accompanied by exclusion of other genera of the Clostridia class (below), may therefore be considered as a commonly reported ‘dysbiotic’ signature.

The human gut is a complex ecosystem where microbes compete and cooperate78. To investigate which of these interactions are associated with disease, we next studied patterns of co-occurrence and mutual exclusivity for the top 20 genera most frequently reported as differentially abundant. For each pair of microbes, we counted the number of signatures where both microbes were reported with either the same or opposite direction of change in relative abundance (Fig. 4a, top and bottom heat maps). This resulted in clusters of co-occurrence driven primarily by functional and phylogenetic similarity, with frequent co-occurrence of genera of the phylum Bacteroidetes and blocks of positive associations in the class Clostridia (Fig. 4a, top red heat map, third cluster from top to bottom). However, clusters of mutual exclusivity displayed clear signs of the established Firmicutes–Bacteroidetes gradient in gut microbiomes68 and a strong negative association between Bacteroides and Prevotella within the phylum Bacteroidetes, as previously reported54 and also observed in healthy samples (r = −0.49, P < 2.2 × 10−16, two-sided Spearman correlation test; Fig. 4a, bottom blue heat map, third cluster from top to bottom). Overall, these patterns largely recapitulated the correlation of these genera in healthy adult stool samples (Fig. 4b), arguing against the existence of specific disease-promoting interactions between these genera. Exceptions were patterns of pronounced mutual exclusivity within the class Clostridia in signatures associated with disease, as observed for Clostridium, Enterococcus and Streptococcus, which were not observed in healthy samples.

Shared and exclusive patterns in pooled microbial signatures

BugSigDB provides opportunities for the discovery of microbial biomarkers and reassessment of previous findings across a much larger and more heterogeneous data source than previously possible. To identify similarities between microbial shifts within and across body sites, we aggregated signatures for one body site at a time and within body sites for one condition at a time (“Signature pooling” in the Methods). To account for differences in sample size between studies, we applied a voting approach where each taxon of a pooled signature obtained a weight based on the aggregated sample size of reporting studies and performed hierarchical clustering based on pairwise similarity between the weighted metasignatures (Fig. 5). Clustering of 27 metasignatures, representing body sites studied by at least two studies in BugSigDB and generated from 1,909 individual signatures of either increased or decreased differential abundance in the study group, resulted in two major body site clusters: a cluster primarily composed of oral and nasal sites and a cluster dominated by vaginal and gastrointestinal sites (Fig. 5a). This clustering was largely invariant to the similarity measure used for clustering (Supplementary Fig. 10), confirming the expected dominant effect of host body site of origin and, in particular, the availability of oxygen.

Fig. 5: Pooled genus-level microbial signatures display robust body site specificity and shared and exclusive patterns of gut dysbiosis between disease phenotypes.
figure 5

a, Clustering of 27 metasignatures (y axis) representing host body sites studied by at least two studies in BugSigDB and generated from 1,909 individual signatures. The clustering shows a separation of vaginal and gastrointestinal sites (top) from oral and nasal sites (bottom). A total of 560 different genera (x axis) colored by phylum were observed across the metasignatures. The numbers of genera in each metasignature are shown in the bar plot on the left. The following three panels display representative genera that frequently occur across metasignatures (‘shared’), predominantly occur in either cluster of body sites (‘compartmentalized’) or were reported in only one body site as differentially abundant (‘exclusive’). The size of each dot corresponds to the relative sample size of studies reporting a genus as differentially abundant. b, Clustering of fecal metasignatures drawn from 34 different conditions, each represented by at least two studies in BugSigDB and generated from 504 individual signatures with increased abundance in the study group. The clustering shows a separation of cancer phenotypes (top) from non-neoplastic diseases (bottom). A total of 248 different genera colored by phylum were observed across the fecal metasignatures. The number of genera in each metasignature is shown in the bar plot on the left.

Within body site clustering of fecal metasignatures of 34 different conditions, studied by at least two studies in BugSigDB and generated from 504 signatures of increased relative abundance in the study group, revealed similarities in reported differential abundance patterns between disease phenotypes (Fig. 5b). This included similarities in microbial shifts for (1) HIV infection and different gastrointestinal cancers, both characterized by chronic inflammation of the gastrointestinal tract and microbial signatures that point to shared pathogenic pathways including tryptophan catabolism and butyrate synthesis79, (2) chronic kidney disease and autism, linked through deleterious copy number variants80,81 and pathogenic gut microbiota-derived metabolites produced by species of the Clostridia class82,83, and (3) type 2 diabetes and schizophrenia, consistent with observations that people with schizophrenia are at increased risk of type 2 diabetes and, conversely, that traditional risk factors for type 2 diabetes are common in people with schizophrenia and can affect the gut microbiome, especially obesity, poor diet and a sedentary lifestyle84. A strong enrichment of genera of the Clostridia class drove the similarity between the metasignatures of Hashimoto’s thyroiditis (12 of 14 genera, 85.7%, P = 2.1 × 10−05, two-sided proportion test) and chronic fatigue syndrome (11 of 17, 64.7%, P = 0.003). Chronic fatigue is common in individuals with Hashimoto’s thyroiditis, a disorder that affects the thyroid gland potentially through a gut microbiota-driven thyroid–gut axis85 that may reduce thyroid hormone production and cause extreme fatigue86.

Genera reported as differentially abundant exclusively in one condition, in particular those replicated by independent studies, are candidates for condition-specific biomarkers (Supplementary Fig. 11, prioritized by relative sample size). This included an exclusive abundance increase of the genus Hungatella for colorectal cancer, which is notable given a reported role of Hungatella hathewayi in driving host colonic epithelial cell promoter hypermethylation of tumor suppressor genes in colorectal cancer87. However, an exclusive decrease in abundance was observed for Marvinbryantia in type 2 diabetes, for which a reduction in Marvinbryantia through treatment with glucagon-like peptide-1 receptor agonist drugs, such as liraglutide, has been shown to contribute to treatment success4.

Discussion

We compiled published signatures of microbial differential abundance in the BugSigDB database, assessed applicability of GSEA methods for enrichment analysis of microbial signatures and identified common patterns of co-occurrence and mutual exclusion in differential abundance results across a broad sample of the human microbiome literature. BugSigDB is a Semantic MediaWiki that allows contribution, review and correction by the microbiome research community and is usable through its web interface and through bulk exports compatible with all major GSEA software. BugSigDB supports any taxon and any host species present in the NCBI Taxonomy and therefore non-human hosts and studies of viromes and eukaryotes. BugSigDB has been initially seeded by ~25 trained curators with more than 2,500 manually curated signatures from the figures, tables, main text and supplementary materials of more than 600 primary publications, providing a broadly relevant collection of machine-readable knowledge of microbial differential abundance. Manually curated metadata include study design, geography, health outcomes, host body sites and experimental and statistical methods. Initial analyses of the database (1) identify human diseases with the highest consistency among independently published signatures of differential abundance, (2) demonstrate the capability of established GSEA methods to prioritize colorectal cancer signatures in the analysis of individual participant colorectal cancer datasets and add to evidence of frequent introgression of oral pathobionts into the gut, (3) demonstrate that the prevalence of genera in fecal microbiomes of control populations is strongly correlated with being reported as decreased across diverse study conditions, (4) confirm the broad relevance of the Firmicutes–Bacteroidetes gradient in shaping common patterns of co-occurrence and mutual exclusivity in the published microbiome literature and identify other less dominant patterns and (5) define sample size-adjusted consensus signatures of body sites and conditions that can simplify and clarify future analyses.

There is concern over replicability of human microbiome studies due to numerous sources of variation in complex experimental and computational quantification pipelines88,89. We propose an approach based on semantic similarity that can be used in systematic reviews to evaluate replication of differential abundance signatures of mixed taxonomic levels reported by independent studies, even when different laboratory methods, such as 16S amplicon sequencing versus MGX, were used. This taxonomy-aware framework provides an effective assessment of the replicability of microbiome differential abundance signatures for subsets of the literature, allowing ranking of the relative replicability of microbiome signatures consisting of different taxonomic levels across many disease phenotypes. We identified signatures associated with antibiotics treatment and chronic inflammation of the gastrointestinal tract as having the highest level of consistency or replication in signatures reported by independent studies.

We pooled signatures across host body sites and experimental conditions to expand the analysis of replicability also to the ~45% of signatures containing fewer than five taxa that are too small individually to be effectively compared between studies. Using a voting approach that weights each taxon by sample size of reporting studies, we constructed consensus ‘metasignatures’, revealing shared and specific patterns of gut dysbiosis by disease phenotype. These metasignatures provide a framework for simplifying the interpretation of results from future studies in the context of the published literature, distinguishing specific from generic results and informing enrichment tests by defining a universe of reported abundance changes for a body site or condition of interest.

Although single-species biomarkers are of primary interest for therapeutic interventions, they are not sufficient for capturing complex ecological patterns of co-occurrence and mutual exclusivity and interactions between microbes that may be relevant to health and disease. Inspection of published signatures in BugSigDB is an alternative approach to studying ecological patterns that complements the analysis of co-occurrence and coexclusion in individual participant metagenomic profiles. In agreement with previous results54, but across a much larger corpus of microbiome studies, we confirmed coexclusions not specific to disease, such as the phylum-level Firmicutes–Bacteroides gradient and the genus-level Bacteroides–Prevotella gradient within the Bacteroides phylum. However, patterns of mutual exclusivity for Clostridium, Enterococcus and Streptococcus genera were specific to disease-associated signatures and were not detectable in healthy samples. To distinguish between disease markers and common false positives in biomarker discovery, we compared stool signatures of disease conditions to prevalence in ~10,000 stool specimens from healthy participants. Across the many diseases and exposures present in BugSigDB, prevalent genera in healthy fecal microbiomes tend to be reported with decreased abundance in diseased-associated fecal microbiomes. On the other hand, genera of low prevalence in healthy fecal microbiomes, such as Enterococcus, tend to be reported with increased abundance in the disease group. Genera such as Lactobacillus and Veillonella, which are prevalent in the stools of healthy individuals and are frequently reported as increased in many study conditions, are more likely false positives or at least are not well suited as candidate biomarkers. Future work can use stratification of BugSigDB signatures by experimental, cross-sectional and longitudinal study designs to better infer causality.

BugSigDB enables systematic comparison of microbial signatures from new microbiome studies to previously published results. Although concepts of GSEA are applicable, microbiome data present new challenges, including smaller signature sizes, taxonomic relationships between features and mixed-taxonomy signatures. We therefore benchmarked two gene set enrichment methods (ORA11 and PADOG31) that performed well in previous benchmarking7 and a recent taxonomic enrichment method (CBEA18) developed specifically for microbiome data. Valid application of gene set enrichment methods is limited to analysis at a single taxonomic level, and ORA further requires defining a realistic feature ‘universe’ and significance threshold for differential abundance; nonetheless, all methods performed well in prioritizing signatures of colorectal cancer across numerous colorectal cancer datasets. As PADOG addresses both (1) shortcomings of ORA in the presence of intermicrobe correlation and (2) compositional bias in signature databases with certain taxa occurring more frequently than others, typically a result of technical or biological sampling bias67,90, we recommend PADOG over ORA for the routine application of enrichment analysis to microbial signatures, especially for datasets with smaller sample sizes where a lack in power typically hinders detecting individually differentially abundant features. Recently emerging microbiome-specific enrichment methods such as CBEA18 have the advantage of taking into account the compositional nature of microbiome data; however, they lack the independent benchmarking and implementation of major GSEA approaches. Although we found CBEA to be a sound alternative to ORA and PADOG, we did not observe notable gains of applying CBEA over PADOG in the colorectal cancer benchmark setup, indicating that basic study characteristics, such as sample size and, to a lesser extent, accounting for correlation within microbial signatures, have a larger impact on identifying relevant signatures in practice than mitigating effects of compositionality.

BugSigDB is a large and diverse collection of the currently available literature on microbial differential abundance and thus also has certain limitations inherent to the currently available literature. Because more than 90% of the studies included in the first release of BugSigDB are based on 16S amplicon sequencing, enrichment analyses were performed at the genus level. However, some genera are functionally heterogeneous, such as streptococci, which groups deadly pathogens with common commensals and useful food-fermenting species. Species- and strain-level variations are neglected, although they can contribute to functional differences between individuals that are important in a clinical context91,92. With the availability of more shotgun sequencing studies in the future, it will be possible to perform enrichment analysis at higher taxonomic resolution up to the species or strain level. Furthermore, studies included in BugSigDB are heterogeneous in their design and execution, including antibiotics exclusion time frames ranging from current use to within the previous year. Restoration of baseline microbial composition following antibiotic treatment typically takes around 1 month in children and 1.5 months in adults, although several common species of the gut microbiome might take substantially longer50,93. Inclusion of this and other study information in BugSigDB allows further investigation and sensitivity analysis into potential sources of heterogeneity in the literature. We anticipate that broader adoption of the recently developed STORMS reporting guidelines for human microbiome studies94 will contribute to more efficient extraction of information from the literature for BugSigDB.

Although natural language parsing programs have the potential to complement the manually curated information in BugSigDB, the majority of the curated information is too complex for currently available text mining algorithms. Natural language parsing programs typically extract patterns from unstructured text on a sentence-by-sentence basis95,96, but BugSigDB standardizes microbial signatures and associated experimental, epidemiological and statistical methods from diverse figures, tables, supplements and textual descriptions that often span multiple sentences. This places many key results outside the current capabilities of text mining applications and necessitates manual curation, but improvements in machine learning, using BugSigDB as a gold-standard dataset, may enable more efficient extraction of published microbiome methods and results in the future. Moreover, automated contributions of signatures from differential abundance software via the BugSigDB application programming interface can streamline standardized reporting of results. As the community contributes additional host species and signatures of microbial physiology and morphology, BugSigDB will dynamically expand, leading to complementary insights and improvements to the systematic interpretation of microbiome studies.

Methods

Definition of semantic concepts (study, experiment, signature and taxon)

Taxon

A taxon, or taxonomic unit, is a unit of any rank (that is, kingdom, phylum, class, order, family, genus, species and strain) designating a microbial organism or a group of microbial organisms.

Signature

A microbial signature or set refers here to a simple unordered list of microbial clades (taxa) sharing a common property or response to a study condition.

Experiment

BugSigDB defines experiments as semantic units within studies and records key characteristics about study subjects, lab analysis, statistical analysis and α-diversity. For study subjects, this includes host species, location, condition, body site, antibiotics exclusion and sample size in study and control sample groups. To define the two sample groups that are contrasted for differential abundance, BugSigDB records the diagnostic criteria applied to define the specific condition/phenotype represented in the study group. Recorded lab analysis fields include sequencing type (16S or MGX) and sequencing platform (such as Illumina or Roche454). For 16S rRNA sequencing, the 16S variable region is also recorded. For the statistical analysis, recorded fields include (1) the statistical test or computational tool applied for differential abundance testing, (2) whether multiple testing correction has been applied to adjust for an inflation of false-positive findings, (3) the significance threshold used to render taxa as differentially abundant, (4) confounding factors that have been accounted for by stratification or model adjustment and (5) factors on which study subjects have been matched on in a case–control study, if applicable.

Study

BugSigDB collects and standardizes microbial signatures from published 16S and MGX microbiome studies. Studies are categorized by study design, and each study is associated with a study identifier such as a PubMed ID and/or a DOI, depending on whether studies are indexed in PubMed.

Data entry, validation and access

BugSigDB is implemented as a semantic MediaWiki30 web interface available at https://bugsigdb.org and supports data entry, semantic validation and web-based programmatic access to annotations for studies, experiments, signatures and individual taxa. The semantic Wiki (1) enforces data entry to follow the nomenclature of the NCBI Taxonomy Database29, (2) enforces metadata annotation of signatures to follow established ontologies and controlled vocabulary for body site28, disease condition32 and type of evidence, (3) provides an application programming interface to access all signatures, potentially filtered on taxonomy and metadata attributes, and (4) allows commenting and error reporting on data elements and relationships. The data curation interface provides type-forward autocomplete to assist with valid data entry (including validation against the NCBI Taxonomy29, Experimental Factor Ontology32, UBERON Anatomy Ontology28 and administrator-defined controlled vocabulary for other fields such as statistical test and sequencing methods) to facilitate organization, filtering and comparison of signatures. External contributions from the community, including signatures, annotations and comments, are supported similar to Wikipedia. The quality of contributions is controlled by tagging contributions as verified after review by a trusted editor, a method for flagging suspect entries and the option to exclude from analysis unreviewed contributions or based on elements of study quality, such as sample size, suspected contamination, paper retraction, batch effects, uncontrolled confounding or a combination of these factors. In addition to standard semantic MediaWiki quality control tags on study level, custom methods are available for flagging taxa according to prevalence in frequently investigated host body sites, inclusion in published contamination blacklists114,115 and absence of known association with a host.

Signatures can be searched and browsed by study and experimental attributes and by individual taxa at https://bugsigdb.org. Bulk export of all signatures and associated metadata is available in plain text formats for use in any programming language and software (including .csv and the .gmt standard used by GeneSigDB23 and MSigDB21) of the current database version or as weekly and semiannual snapshots. The companion bugsigdbr R/Bioconductor package (bioconductor.org/packages/bugsigdbr) provides advanced features, such as ontology-based filtering, limitation of taxonomic level, look-up of individual signature and taxon pages and conversion to application-centric formats. The exported files are compatible with most enrichment software and are included by default in our lab’s EnrichmentBrowser R/Bioconductor package55 to facilitate a large number of GSEA methods and visualizations. The BugSigDBStats R package continuously integrates with https://bugsigdb.org and provides weekly updated database statistics in an HTML report page (https://waldronlab.io/BugSigDBStats).

Signature similarity

Signature similarity was computed based on two different measures: (1) Jaccard index based on pairwise overlaps between signatures harmonized to genus level and (2) semantic similarity between signatures of mixed taxonomic levels. Pairwise calculation of Jaccard similarity for genus-level signatures was performed using the calcJaccardSimilarity function of the BugSigDBStats package. Genus-level signatures from BugSigDB were obtained using bugsigdbr’s getSignatures function. Taxonomic clades given at a more specific taxonomic level (species or strain) were transformed by cutting the taxonomic tree at the genus level. The Jaccard index, also known as the Jaccard similarity coefficient, is defined as the size of the intersection divided by the size of the union of two input signatures A and B:

$$J(A,B)=\frac{{\rm{|}}A\cap B{\rm{|}}}{\left|A\cup B\right|}$$
(1)

Note that, by design, 0 ≤ J(A, B) ≤ 1.

Semantic similarity was computed based on Lin’s measure of semantic similarity116, as implemented in the ontologySimilarity package from the ontologyX R package series117. Measures of semantic similarities have been proposed for comparing concepts within a taxonomy118 with numerous applications for biomedical ontologies36. Of note, semantic similarity has conceptual parallels with the computation of UniFrac distance119. Here, individual taxa can be considered nodes of the NCBI Taxonomy when represented as a directed acyclic graph. Computing semantic similarity between two taxa corresponds then to computing a topological similarity such as the shortest path linking the two taxon nodes. More specifically, computing Lin’s measure of semantic similarity between two taxa corresponds to computing the information content of the lowest common ancestor (LCA) of the two taxa116. The more frequently a taxon occurs (that is, the higher its probability of being an ancestor of other taxa in the taxonomy), the lower its information content. If the LCA of two taxa corresponds to a taxon at a higher taxonomic level, these taxa are not very similar, and this is reflected in a low-information content of their LCA. Given pairwise semantic similarities between individual taxa, semantic similarity between two signatures (that is, two taxon sets) is then obtained using a best-match average combination approach120 where each taxon of the first signature is paired with the most similar taxon of the second signature and vice versa.

Bug set enrichment analysis

Metagenomic datasets providing species-level relative abundance for fecal microbiomes of individuals with colorectal cancer and healthy individuals were obtained from curatedMetagenomicData 3.0 (ref. 54). Relative abundance proportions were multiplied by read depth and rounded to the nearest integer before obtaining integer read counts compatible with bulk RNA-sequencing (RNA-seq) tools for differential expression analysis and GSEA. For genus-level analysis, species-level counts were summed across branches using the splitByRanks function from the mia package. Given a recent assessment that reported good performance of bulk RNA-seq tools for microbiome data121, differential abundance analysis was performed following the limma-trend approach122. Read counts were transformed to log (counts per million) using edgeR’s cpm function with a prior count of 3 to dampen the variances of logarithms of low counts123. Genus- and species-level microbial signatures from BugSigDB were obtained using bugsigdbr’s getSignatures function. To keep signatures meaningfully sized, taxonomic clades given at the species or strain level were transformed by cutting the taxonomic tree at the genus level. ORA and PADOG were performed as implemented in the EnrichmentBrowser package55. CBEA was performed as implemented in the CBEA package18.

Taxon co-occurrence

Genus-level signatures from BugSigDB were obtained using the getSignatures function of the bugsigdbr package. Signatures were filtered by body site for fecal samples and stratified by direction of abundance change (increased/decreased). The top 20 most frequently occurring genera in the resulting signatures were reported. Prevalence of these genera was computed as percentage of healthy adult stool samples in which the genus was observed at non-zero relative abundance in metagenomic datasets from curatedMetagenomicData 3.0 (ref. 123). The correlation between prevalence in healthy samples and the proportion of signatures with increased abundance for the top 20 genera was assessed using a Pearson’s correlation test as implemented in the cor.test function of the stats package. Taxon co-occurrence in signatures associated with disease was contrasted against Spearman rank correlation of the top 20 genera in healthy samples using the cor function of the stats package.

Signature pooling

Signatures were pooled for one body site at a time and within body sites for one condition at a time, as implemented in the getMetaSignatures function of the bugsigdbr package. Taxa within a pooled signature were weighted based on the aggregated sample size of the studies that reported this taxon as differentially abundant divided by the total sample size of studies contributing to the pooled signature. Resulting weighted metasignatures were clustered by semantic similarity, where the weights were incorporated into the best-match average combination approach36, as implemented in the weightedBMA function of the BugSigDBStats package. Analysis was restricted to body sites and conditions studied by at least two studies in BugSigDB and containing at least five taxa in the resulting pooled signature. Robustness of the clustering was evaluated by comparing to rank-biased overlap124.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.