Bakdrive: identifying a minimum set of bacterial species driving interactions across multiple microbial communities

Abstract

Motivation: Interactions among microbes within microbial communities have been shown to play crucial roles in human health. In spite of recent progress, low-level knowledge of bacteria driving microbial interactions within microbiomes remains unknown, limiting our ability to fully decipher and control microbial communities.

Results: We present a novel approach for identifying species driving interactions within microbiomes. Bakdrive infers ecological networks of given metagenomic sequencing samples and identifies minimum sets of driver species (MDS) using control theory. Bakdrive has three key innovations in this space: (i) it leverages inherent information from metagenomic sequencing samples to identify driver species, (ii) it explicitly takes host-specific variation into consideration, and (iii) it does not require a known ecological network. In extensive simulated data, we demonstrate identifying driver species identified from healthy donor samples and introducing them to the disease samples, we can restore the gut microbiome in recurrent Clostridioides difficile (rCDI) infection patients to a healthy state. We also applied Bakdrive to two real datasets, rCDI and Crohn's disease patients, uncovering driver species consistent with previous work. Bakdrive represents a novel approach for capturing microbial interactions.

Availability and implementation: Bakdrive is open-source and available at: https://gitlab.com/treangenlab/bakdrive.

Figure 1.

FMT process simulation and driver species identification. (a) The ecological network G₀ of a metacommunity with 100 species. Red/blue arrows represent positive/negative relationships respectively. The widths of arrows reflect interaction strengths. (b) A subset of species and their corresponding interaction networks G₁ are randomly selected from the pool of 100 species. (c) Simulation of species abundance over time. The thick red line represents the abundance of C. difficile. The FMT process simulation includes the pre-FMT bowel cleansing process, which leads to sudden drops of species abundance in pre-FMT samples. (d) Driver species identification. Directed and weighted bacteria interaction networks are converted to an undirected, unweighted multilayer network where layers are unconnected. Here, the multilayer network consists of three layers: $G_1$($V_1, E_1)$, $G_2$($V_2, E_2)$, $G_3$($V_3, E_3)$. Green nodes represent driver nodes found using the MDSM algorithm. All non-driver nodes directly connect to at least one driver node.

Figure 2.

Efficacy of driver species transplantation with universal ecological networks. (a) Universal. (b) Rewired ecological networks. We simulate the driver species transplantation process with universal/rewired microbial dynamics. In each box plot, the Y-axis shows the recovery degree of C.difficile after the driver species transplantation over 1000 simulated engraftments. The X-axis represents the actual number of species introduced to the diseased samples, while $n$ is the number of driver species identified from a given multilayer network, showing that there is an increase in efficacy as a greater number of the identified driver species are added. The orange boxes represent the recovery degree of C.difficile after colonizing driver species. The green boxes represent the recovery degree of C.difficile after introducing the same number of species, which are randomly selected from the global pool of 100 species. The significant differences of recovery degree by Wilcoxon rank-sum test are labeled by asterisks (*$P$< 0.05).

Figure 3.

Principal coordinate analysis (PCoA) of simulated diseased (red), donor (green), and after driver species transplant (ADT) samples (blue). In each panel, we have 1000 simulations and thus 1000 samples of each type. The panels vary by the number of layers $N$ (input samples) used to calculate the MDS (1, 5, 10, 25, 50, 100). PCoA was done using Bray-Curtis dissimilarity. As $N$ increases the red and blue clusters (pre- and post-transplantation, respectively) become gradually more separated, showing that the MDS affects not just C.difficile abundance but the overall community.

Figure 4.

Experimental setup of rCDI actual FMT versus simulated MDS comparison. Microbiome compositions were taken from actual patient data where samples were available from (a) patient before FMT, (b) donor, and (c) patient after FMT (all shown in blue). 12 patients were considered where eventual recovery was attained (blue). The set of donor samples were used as input to Bakdrive to generate a mock driver species (MDS) community (orange). The actual pre-FMT samples were given a simulated engraftment of the MDS, with the resulting community development simulated (green). Of interest is whether the actual changes following FMT (i.e. (1) to (2)) resemble the simulated changes following MDS engraftment (i.e. (1) to (3)).

Figure 5.

The consistency of species abundance changes between real and simulated samples. (a) Each column represents a patient sample. Each row represents a species that exists in 12 patient samples or belongs to driver species. “*” marks driver species. The driver species abundance changes are labeled as cross, as most driver species do not exist in diseased samples. The size of the bubble represents the percentage of species in patient samples before FMT. Gray shows species abundance changes in opposite directions between simulated and real dataset. Red/blue indicates a given species abundance increases/decreases after FMT in both simulated and real data respectively. White (“None”) indicates that the microbe was not present in the pre-FMT sample and thus no change was measured. In the upper bar plot, agreement quantifies the percentage of species in renormalized diseased samples (3) whose abundances shift in the same direction in both real and simulated data. (b) Recovery degree of pathogens. x-axis represents the target pathogens. Y-axis represents the recovery degree of target pathogens in each sample. We then zoom in the bar plot and only show samples with recovery degree above 0. (c) Principal coordinate analysis (PCoA) of real diseased, after FMT and after-MDST samples. Each sample is represented as a dot. There are a total of 36 samples from 12 patients. Bray-Curtis dissimilarity between samples are computed.