FastST: an efficient tool for inferring decomposition and directionality of microbial communities

Introduction

The microbiome refers to a collection of microorganisms inhabiting a specific environment, which forms complex and diverse communities of multiple interacting species (Castro-Nallar et al., 2015). Previous studies have demonstrated the profound impact of the microbiome in complex diseases such as diabetes, inflammatory bowel disease, colorectal cancer, and allergy outcomes (Marchesi et al., 2011; Qin et al., 2012; Kostic, Xavier & Gevers, 2014; Fujimura & Lynch, 2015), which indicates the potential of human microbes as biomarkers for disease diagnosis or as therapeutic targets for treatment (Manor et al., 2020). Microbes also play a crucial role in various environments, along with the capacity to maintain a healthy global ecosystem responding to the climate change (Cavicchioli et al., 2019) and participating in essential biogeochemical cycling events, such as carbon and nitrogen fixation (Gougoulias, Clark & Shaw, 2014).

Prediction of the structure and dynamics of the microbial community is essential for understanding microbiome development and ecosystem function (Meroz et al., 2021). Insights into the origins of microbial communities and the movement of microbes across different ecosystems can help unravel the rules that govern microbial ecology. The key, albeit challenging, step in this analysis is to characterize the composition of microbial communities, as they typically comprise multiple source environments, including various contaminants and other microbial communities that have interacted with the sampled habitat (Shenhav et al., 2019). To resolve this, microbial source tracking (MST) approaches have been presented, which aim to quantify the proportion of different microbial samples (sources) in a target microbial community (sink) (Shenhav et al., 2019; Knights et al., 2011; McGhee et al., 2020). Through microbial source tracking (MST), transmissions between different hosts and environments can be tracked and similarities between microbiome samples can be identified (Briscoe, Halperin & Garud, 2023). MST can be applied to the determination of discrete sources of fecal pollution (Unno et al., 2018), tracking bacterial contamination routes of municipal water (Liu et al., 2018) and natural fresh water systems including rivers and lakes (Staley et al., 2018), and disease prevention (Fu & Li, 2014).

The most recent approaches for microbial source tracking includes SourceTracker2 (Knights et al., 2011), FEAST (Shenhav et al., 2019), STENSL (An et al., 2022), and SourceID-NMF (Huang, Cai & Sun, 2024) methods. These methods use species abundance profiles of the sample of interest and potential sources to calculate the percentages of sinks that are attributable to each potential source. Based on a Bayesian framework, SourceTracker2 employs Markov chain Monte Carlo (MCMC) to estimate contamination proportions in metagenomics studies. On the other hand, FEAST determines the fraction of each source environment in a target microbial community by modeling mixture proportions for various source microbial samples in a given sink sample, using expectation maximization for improved computational efficiency. STENSL, derived from the FEAST algorithm, incorporates unsupervised source selection via least-squares optimization with L1-norm regularization. SourceID-NMF employs non-negative matrix factorization (NMF) to identify source contributions to sink communities.

Despite their advancements, current microbial source tracking approaches have limitations that warrant attention for practical applications. Both SourceTracker2 and FEAST are computationally expensive, especially when handling high-dimensional microbiome data, though the latter is considered relatively efficient. Furthermore, these approaches lack the capability to infer directionality, i.e., the requirement to predetermine which is the source/sink prior to analysis. In certain real-world microbiome studies—for example, in hospital sink drains (Laço et al., 2025) and newly opened hospital wards (Lax et al., 2017)—it can be unclear which community is acting as the source and which as the sink, underscoring the critical need for methods to discern the directionality of microbial exchange. Addressing these limitations will enhance the reliability and applicability of microbial source tracking methods.

In the present study, we propose a novel inference framework, named FastST. FastST is able to accurately and efficiently estimate the relative contributions of the sources to the sink of interest, as well as infer the directionality if source and sink relationships are not known a priori. The estimation accuracy and efficiency are achieved by leveraging the generalized least squares (GLS) inference based on a multiple linear regression model. The inference of directionality is further made by selecting the maximized joint likelihood for observing the sink and source data across all enumerated Bayes network models. Through simulations, we demonstrate the outperformance of FastST in terms of mean absolute error and running time, as compared to SourceTracker2, FEAST, STENSL, and SourceID-NMF. In addition, simulations also show that FastST is able to identify the directionality in the source and sink environments.

Methods

Consider a microbiome data of K + 1 observed samples, each consisting of N taxa. We assume that, among the K + 1 samples, one is treated as the sink and the rest as K sources. Denote the observed sink and source data by x = (x₁, …, x_N)^T and y_i = (y_i1, …, y_iN)^T, 1 ≤ i ≤ K, respectively, based on the corresponding random variable model: ${\displaystyle \mathit{X}\sim \mathrm{multinom}\left(C,{\left({\beta }_1,\dots ,{\beta }_N\right)}^T\right)}$ ${\displaystyle {\mathit{Y}}_{\mathbf{i}}\sim \mathrm{multinom}\left(C_i,{\left({\gamma }_{i1},\dots ,{\gamma }_{iN}\right)}^T\right),1\le i\le K}$ where $C={\sum }_{j=1}^N{x}_j,C_i={\sum }_{j=1}^N{y}_{ij}$ are the total taxa counts of the sink and the ith source, and β_j, γ_ij denote the relative abundance of taxa j in the sink and the ith source, respectively. Further assume that, the sink-source relationships are determined by the following mixed-proportion model: (1)${\displaystyle {\beta }_j=\sum _{i=0}^K{\alpha }_i{\gamma }_{ij},1\le j\le N.}$ In this model, an unobserved source Y₀ ∼ multinom(C₀, (γ₀₁, …, γ_0N)^T) is assumed to contribute to the sink with a proportion α₀. Our goal is to unveil the dependence structure between the sink and the sources by estimating the parameters α_i based on the observed data x and y_i, for 1 ≤ i ≤ K.

Generally speaking, the estimation of α_i can be achieved by maximizing the joint likelihood of observing (x, y₁, …, y_K). Since the mixed-proportion model depends on several latent variables, i.e., the taxa abundances β_j, γ_ij in the sink and sources and the variables in the unobserved source, the estimation is usually carried out iteratively. Typical examples include finding (local) maximum likelihood estimates (MLE) through the Expectation–Maximization (EM) algorithm (Shenhav et al., 2019), or finding maximum a posteriori (MAP) estimates by Gibbs sampling (Knights et al., 2011). Although convergence is generally ensured for the MLE or maximum a posteriori (MAP) estimation, from a practical perspective, an obvious disadvantage of these iterative methods lies in their lack of computational efficiency, especially when the number of latent or unobserved variables is high. This greatly affects the practical use of source tracking in high-dimensional microbiome data. To address this issue, we simplify the mixed-proportion model by assuming γ_ij to be known (if not, they can be approximated by the sample estimates y_i/C_i), for 1 ≤ i ≤ K, 1 ≤ j ≤ N, and propose to estimate α_i’s by the generalized least square (GLS) method. This approach is described below.

Multiplying C to both side of Eq. (1), we have ${\displaystyle E\left[X_j\right]=\sum _{i=0}^{K}C{\alpha }_i{\gamma }_{ij},1\le j\le N.}$

For the unobserved source Y₀, since both its taxa abundance (γ₀₁, …, γ_0N)^T and its proportion α₀ are unknown, the estimation of these parameters certainly encounters a nonidentifiability issue. In addition, it is reasonable to treat Y₀ as a nuisance contributor to the sink with a negligible proportion, so we may assume that the product α₀γ_0j is common across 1 ≤ j ≤ N. Letting ${\tilde{\alpha }}_0=C{\alpha }_0{\gamma }_{0j},1\le j\le N$ and ${\tilde{\alpha }}_i=C{\alpha }_i,1\le i\le K$, the above equation describes a multiple linear regression model (2)${\displaystyle X_j={\tilde{\alpha }}_0+\sum _{i=1}^K{\tilde{\alpha }}_i{\gamma }_{ij}+{ϵ}_j,1\le j\le N,}$

where ϵ_j denote the random error term with E[ϵ_j] = 0 and ${\displaystyle Cov\left(ϵ\right)=\left(\begin{array}{cccc}\hfill C{\beta }_1\left(1-{\beta }_1\right)\hfill & \hfill -C{\beta }_1{\beta }_2\hfill & \hfill \cdots \hfill & \hfill -C{\beta }_1{\beta }_N\hfill \\ \hfill -C{\beta }_1{\beta }_2\hfill & \hfill C{\beta }_2\left(1-{\beta }_2\right)\hfill & \hfill \cdots \hfill & \hfill -C{\beta }_2{\beta }_N\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill -C{\beta }_1{\beta }_N\hfill & \hfill -C{\beta }_2{\beta }_N\hfill & \hfill \cdots \hfill & \hfill C{\beta }_N\left(1-{\beta }_N\right)\hfill \end{array}\right).}$

Replacing β_j by its sample estimate x_j/C and denoting this approximated Cov(ϵ) by Σ, we obtain the GLS estimate of $\tilde{\alpha }={\left({\tilde{\alpha }}_0,{\tilde{\alpha }}_1,\dots ,{\tilde{\alpha }}_K\right)}^T$ by ${\displaystyle \widehat{\stackrel{̃}{\alpha }}={\left({\gamma }^T{\Sigma }^{-1}\gamma \right)}^{-1}\left({\gamma }^T{\Sigma }^{-1}x\right),}$

where γ is an N × (K + 1) design matrix with first column being a vector of 1’s. Consequently, the estimate of α = (α₀, α₁, …, α_K)^T is (3)${\displaystyle \hat{\alpha }=\frac{1}{C}{\left({\gamma }^T{\Sigma }^{-1}\gamma \right)}^{-1}\left({\gamma }^T{\Sigma }^{-1}x\right).}$

It is known that, as a multivariate covariance matrix, Cov(ϵ) is positive semidefinite and does not have a unique inverse. To overcome this difficulty and calculate the GLS estimate of $\tilde{\alpha }$, we may either replace Cov(ϵ) with a non-singular matrix by removing any row and corresponding column (and accordingly reduce dimension for γ and x) (May & Johnson, 1998), or alternatively use the Moore–Penrose inverse (pseudoinverse) of Cov(ϵ) (Tanabe & Sagae, 1992).

In practice, it is often not clear which of the K + 1 samples in the microbiome data plays the role of the sink. Denoting the samples by {S₁, …, S_K+1}, we define the kth model, 1 ≤ k ≤ K + 1, by

${\displaystyle \mathrm{Model}k:\mathbf{X}=S_k,{\mathit{y}}_i\in S_{k-},1\le i\le K,}$

where S_k−: = {S₁, …, S_K+1}∖S_k. Our purpose is to determine the correct model from a total of K + 1 candidate models so that the true source and sink relationships are testified. For convenience, let us represent the directionality of Model k by S_k− → S_k. Using a Bayesian network setting, we write the joint likelihood of observing (x, y₁, …, y_K) in Model k as (4)${\displaystyle \begin{array}{cc}\hfill Pr\left(\mathbf{X},{\mathit{y}}_1,\dots ,{\mathit{y}}_K\right)& =Pr\left(\mathbf{X}|{\mathit{y}}_1,\dots ,{\mathit{y}}_K\right)Pr\left({\mathit{y}}_1,\dots ,{\mathit{y}}_K\right)\hfill \\ \hfill & =Pr\left(\mathbf{X}|{\mathit{y}}_1,\dots ,{\mathit{y}}_K\right)\prod _{i=1}^{k}Pr\left({\mathit{y}}_i\right)\hfill \end{array},}$

where the conditional likelihood Pr(x|y₁, …, y_K) is calculated by using the estimated parameters from Eq. (3). We then choose the correct model by maximizing Pr(x, y₁, …, y_K).

We note that, this directionality inference method employs the concept of structure learning for Bayesian networks. Briefly speaking, a Bayesian network is a graphical model that represents a set of random variables and their conditional dependencies via a directed acyclic graph (DAG), and structure learning for Bayesian networks refers to learning the structure of the DAG from data. As shown in Fig. 1, we model the microbiome data by a simple Bayesian network which connects each source y_i, 1 ≤ i ≤ K to the sink x with a directed edge. Our directionality inference is essentially a score-based approach which uses the log-likelihood of the data under the graph structure as a criterion (or score function) to evaluate how well the Bayesian network fits the data. In practice, score-based approaches usually introduce a regularization term (e.g., Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC)) to penalize overfitting of the model and favor simpler models. However, for our directionality inference problem, the regularization term can be omitted safely in the score function because all candidate graphical models are of the same complexity—it is just a matter which variable plays the role of the sink.

Results

Scenario 3: evaluation on real microbiome data

Finally, we evaluated FastST in a fully real-world setting using the Knights et al. dataset, which comprises four 16S rRNA microbiome samples from environmental surfaces—one from a hospital, two from research laboratories, and one from an office building—designated as sinks, along with 180 source samples collected from diverse environments, including human skin, oral cavities, feces, and temperate soils. In this case, the ground-truth source proportions were unknown; therefore, we focused on comparing the patterns of estimated contributions across different tools.

From the result (Fig. 2) for most sink samples, FastST produced source contribution estimates largely consistent with those from FEAST, with only minor discrepancies. Both methods indicated substantial contributions from skin and soil sources for the laboratory and office sinks, while the NICU sample exhibited a broader diversity of contributors. In contrast, FEAST and SourceTracker2 often assigned a dominant proportion of the composition to a single environment type, most frequently skin. STENSL generated patterns similar to FEAST, reflecting its algorithmic derivation from FEAST. Meanwhile, SourceID-NMF showed limited diversity in inferred sources and failed to generate results for two sink samples (Lab 1 and Lab 2), leading to the entire proportion being attributed to unknown sources.

Discussion and Conclusions

In this study, we introduced FastST, a novel microbial source tracking method designed to address the limitations associated with existing approaches such as FEAST, SourceTracker2, STENSL, and SourceID-NMF. FastST predicts the contributions of the sources to the sink, and infers the directionality if the source and sink relationships are not pre-defined. One unique feature of FastST is that it transforms the mixed-proportional sink-source relationships into a multiple linear regression model, which greatly simplifies the inference procedure. The source contributions were predicted by the standard GLS method and the directionality was inferred by selecting the maximized joint likelihood of a Bayes network model. Through simulations with different number of known and major sources, we demonstrated that FastST was able to achieve significant improvements in computational efficiency and accuracy, particularly evident in scenarios involving a large number of known sources. Unlike other methods, FastST exhibited good scalability, meaning that it effectively maintained consistent run-time and accuracy regardless of increasing complexity in terms of data size and dimensionality.

One of the long-standing obstacles to the practical use of microbial source tracking lies in computational efficiency. As shown in the simulation study, when the complexity of microbial communities increases, MST methods usually encounter significant computational difficulties. Due to its Bayesian setup and non-negative matrix factorization, SourceTracker2 and SourceID-NMF are computationally expensive. FEAST and STENSL, while significantly faster than SourceTracker2 at a lower number of known sources, suffered a drastic increase in computation time with a larger dataset. FastST, on the other hand, maintained stable performance with minimal runtime increase, demonstrating suitability for practical applications involving a large number of sources. Additionally, FastST’s capability to accurately infer directionality without predefined source–sink relationships enhances its applicability in real-world scenarios, where distinguishing sources from sinks can be challenging. The high accuracy rates achieved in identifying the correct source–sink relationship underscore its potential for various ecological and public health applications.

It is noteworthy that, our simulation study sets the number of taxa N to be greater than the number of source samples K, which is intuitively reasonable for real applications. When this assumption is not satisfied, for example, there exists a large number of source samples, a potential solution is to conduct a dimension reduction (e.g., PCA) before building a mixed-proportion model for the sink-source relationships. After all, when K is large, both the parameter of the taxa relative abundances γ_ij and the parameter of proportions α_i become difficult to estimate. As a rule of thumb, when K is greater than 100, most of the minor sources have very tiny, negligible contributions to the source, making the inference on α_i worthless.

In conclusion, FastST has a potential to offer a substantial advancement in microbial source tracking by providing rapid and precise estimations and reliably inferring directionality. However, different scenario cases in simulation testing reveals that FastST has robust performance in both accuracy and computational efficiency, and also provide a directionality inference method which hasn’t been firstly presented to the best of our knowledge. Future research should explore extending FastST to incorporate more complex environmental scenarios with real datasets.

Supplemental Information