Differential network connectivity analysis for microbiome data adjusted for clinical covariates using jackknife pseudo-values

Abstract

Background

A recent breakthrough in differential network (DN) analysis of microbiome data has been realized with the advent of next-generation sequencing technologies. The DN analysis disentangles the microbial co-abundance among taxa by comparing the network properties between two or more graphs under different biological conditions. However, the existing methods to the DN analysis for microbiome data do not adjust for other clinical differences between subjects.

Results

We propose a Statistical Approach via Pseudo-value Information and Estimation for Differential Network Analysis (SOHPIE-DNA) that incorporates additional covariates such as continuous age and categorical BMI. SOHPIE-DNA is a regression technique adopting jackknife pseudo-values that can be implemented readily for the analysis. We demonstrate through simulations that SOHPIE-DNA consistently reaches higher recall and F1-score, while maintaining similar precision and accuracy to existing methods (NetCoMi and MDiNE). Lastly, we apply SOHPIE-DNA on two real datasets from the American Gut Project and the Diet Exchange Study to showcase the utility. The analysis of the Diet Exchange Study is to showcase that SOHPIE-DNA can also be used to incorporate the temporal change of connectivity of taxa with the inclusion of additional covariates. As a result, our method has found taxa that are related to the prevention of intestinal inflammation and severity of fatigue in advanced metastatic cancer patients.

Conclusion

SOHPIE-DNA is the first attempt of introducing the regression framework for the DN analysis in microbiome data. This enables the prediction of characteristics of a connectivity of a network with the presence of additional covariate information in the regression. The R package with a vignette of our methodology is available through the CRAN repository (https://CRAN.R-project.org/package=SOHPIE), named SOHPIE (pronounced as Sofie). The source code and user manual can be found at https://github.com/sjahnn/SOHPIE-DNA.

Background

The human microbiome is the collective genomes of microbes or micro-organisms localized to the various sites of human body [1]. Recent clinical studies have shown that the microbiome has a regulatory role in a wide array of illnesses in humans, such as cancer [2], human immunodeficiency virus [3], and inflammatory bowel disease (IBD) [4]. Moreover, the human microbiome is linked to emotional well-being [5] and mental health including depression [6], autism spectrum disorders [7], and human brain diseases [8].

Following the advent of next-generation sequencing technologies, the taxonomic composition of microbial communities is better characterized by the amplification of small fragments (or amplicon) of the 16S ribosomal RNA (or 16S rRNA) gene. More recently, shotgun metagenomic sequencing has become an alternative for microbial community profiling [9]. Either sequencing platform typically employs similarity-based clustering algorithms to group 16S rRNA sequences into Operational Taxonomic Units (OTU) [10, 11] that are compositional.

The applications of network theory have been successfully utilized to better appraise the complex symbiotic (or dysbiotic) relationship between microbiome and disease states – microbial co-abundances [12]. The abundance matrix or the observed OTU table is used to infer microbial co-abundances among taxa through either correlation-based approaches or probabilistic graphical models.

The differential network (DN) analysis compares the network properties between two or more graphs under different biological conditions, such as degree centrality. Based on the recent review article [13], there are two methods that are newly available to the DN analysis for microbiome data: Microbiome Differential Network Estimation (MDiNE) [14] and Network Construction and comparison for Microbiome data (NetCoMi) [15]. These methods, however, do not assume that the association structure depends on additional binary and continuous covariates.

It has been recognized that the composition of the gut microbiome is central to the pathogenesis of IBD [4, 16]. In addition, the gut microbiome composition in patients with IBD is largely influenced by various factors including the use of antibiotics, diet, and cigarette smoking [4]. In an analogous fashion, it is not unreasonable to speculate that the structure of the microbial networks can also vary depending on these factors. Thereby, there is a need for statistical methods for DN analysis that can include additional predictor variables.

One way to accomplish this goal is to use a regression technique based on pseudo-values, a component to calculate the bias-corrected estimator of leave-one-out jackknife resampling procedure [17]. The pseudo-value technique was first postulated by Andersen and his colleagues [18, 19] in the context of multi-state survival models with right-censored data. Since then, it has been well studied in various disciplines of statistics including the interval-censored data [20, 21], clustered data [22, 23], and machine learning methods [24, 25].

The ultimate benefit of this technique is its straightforward inclusion of additional covariates in the generalized linear model [26]. An asymptotic linearity and consistency of pseudo-values given covariates are shown with the second-order von Mises expansion [27, 28]. The pseudo-values can then be used as the response variable in a regression model with the covariates [29]. Several studies reported that the type I error is well controlled at a nominal level of 0.05 while maintaining a high statistical power under the quasi-likelihood generalized linear mixed model [30] and generalized estimating equations framework [23, 31] for pseudo-value regression approach.

Hence, we propose a regression modeling method for differential connectivity (DC) analysis that regresses the jackknife pseudo-values calculated from a degree centrality of taxa in a microbial network to directly estimate the effects of predictors. The primary focus of the methodological innovation presented in this manuscript is centered on DC analysis, a subset of the broader DN analysis. The findings of DC analysis primarily describe the DC of individual nodes (or taxa) instead of taxon-taxon co-abundance relationships [32]. In this approach, the grouping variable itself could also be included in the regression model along with additional clinical covariates while regressing the pseudo-values. We loosely refer to this as a “multivariable setting”, whereas in “univariable settings” only the grouping variable is utilized in a DN analysis.

In the present study, we introduce Statistical ApprOacH via Pseudo-value Information and Estimation for Differential Network Analysis (SOHPIE-DNA) that can include covariate information in analyzing microbiome data. We firstly demonstrate the plausibility of the proposed method by comparing the model performances with MDiNE and NetCoMi through simulations under multivariable and univariable settings. Of note, the covariate adjustment is the main strength of our proposed method. Therefore, the findings from multivariable simulation setting corroborates the main reason for our methods paper. Furthermore, the SOHPIE-DNA is applied to illustrate its clinical utility by examining real data from the American Gut Project [33] and the Diet Exchange Study [34] to identify DC taxa with presence of covariates. All statistical analyses are performed in R version 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria).

Results

Simulation study

The sample size $n = 20, 50, 200, 500$ are considered for each microbial network with $p = 20, 40$ taxa over 1,000 Monte Carlo replicates. Simulations are repeated to assess the effects of covariates on taxa by changing the effect size, $\delta = 0.05, 0.1, 0.2$, which is described in “Simulated data” section. A new network is generated at each simulation replicate to account for biological variability of the network structure.

The performance metrics provided in “Performance measures” section are computed by comparing the test results with the true network. In the true network setting, a taxon is truly DC between groups if it is connected to at least one different neighbor taxon between groups. Tables 1 and 2 summarize simulation results under the multivariable setting. That is, a continuous covariate is included with the binary group variable in the regression model. To illustrate the utility of the proposed method on covariate-dependent network, we compared the pseudo-value regression approach with the recent methods available (NetCoMi and MDiNE) that cannot incorporate the additional covariate. Results show that the SOHPIE-DNA consistently maintains high recall values in all specifications of taxa, sample sizes, and effect sizes, and outperforms NetCoMi and MDiNE in almost all cases. A higher F1 score of SOHPIE-DNA indicates that the proposed method can achieve a better overall model performance in the presence of additional covariates, compared with the two competing methods. In general, all metrics improve as n increases and/or when the larger effect size is provided ($\delta =0.2)$, as expected. It is worth noting that the MDiNE poses a practical challenge associated with substantially large computational time and costs. For instance, it requires more than 9 days to complete each simulation for $p=40$ and $n=200$ from the University of Florida Research Computing Linux server, HiPerGator 3.0 with 32CPU cores and 4GB of RAM per node, while it takes up to 18 h to execute the same simulation tasks for both the SOHPIE-DNA and NetCoMi with 4CPU cores and 6GB of RAM per node. See Additional file 1: Table S1 for more details.

Table 3 presents results of the univariable setting, where only the binary group variable is included in the model. In other words, only the effect of group was considered when generating random networks. On the whole, a similar pattern is shown in the univariable setting that the SOHPIE-DNA reaches a higher level of recall, compared with NetCoMi and MDiNE. Overall, our method resulted in a higher F1 score when the smaller network is considered. All of the methods suffer from a low precision with a small effect size ($\delta =0.05)$, but eventually improves with a larger effect size ($\delta =0.2)$.

Analysis of the American Gut Project Data

Six out of 138 taxa are found significantly DC between migraineurs vs. non-migraineurs while adjusting for age, sex, exercise frequency, categorical alcohol consumption, oral hygiene behavior, and dog ownership. At the family-level, the DC taxa are members of Ruminococcaceae, Lachnospiraceae, Enterobacteriaceae, Erysipelotrichaceae, and Bacteroidaceae. Of these families, the absence of Lachnospiraceae has been linked to the active or severe Clostridium difficile infection [35]. Erysipelotrichaceae has been associated with dyslipidemic phenotypes and systemic inflammation [36]. Moreover, a recent study [37] reported that the species enriched among migraineurs include Ruminococcus gnavus and Lachnospiraceae bacterium. The computational time for our analysis was about 12 h on the high-performance Linux cluster, HiPerGator 3.0 with 16CPU cores and 4GB of RAM per node.

Analysis of the diet exchange study data

Out of 112 taxa, 16 are predicted to be significantly DC between AA and RA after the two-week dietary exchange intervention while accounting for their age and BMI group. A complete list of DC taxa represent Bacillus, Bacteroides uniformis et rel., Bacteroides vulgatus et rel., Clostridium ramosum et rel., Coprococcus eutactus et rel., Eggerthella lenta et rel., Escherichia coli et rel., Eubacterium hallii et rel., Eubacterium siraeum et rel., Faecalibacterium prausnitzii et rel., Prevotella oralis et rel., Roseburia intestinalis et rel., Ruminococcus gnavus et rel., Staphylococcus, Uncultured Bacteroidetes, and Xanthomonadaceae. Notably, Roseburia intestinalis contributes to the prevention and management of intestinal inflammation and atherosclerosis [38]. Eubacterium hallii has been negatively associated with the fatigue severity scores of patients with advanced metastatic cancer [39]. The analysis took about an hour and 11 min on the HiPerGator 3.0 with 16CPU cores and 4GB of RAM per node.

Discussion

In this manuscript, we introduce the SOHPIE-DNA, a pseudo-value regression approach that determines whether a microbial taxa is significantly DC between groups after adjusting for additional covariates. This study is the first of its kind in the literature to develop a regression modeling for the DN analysis in microbiome data, which includes more than one predictor (e.g., group) in the model and predicts features of connectivity of a network. A simulation study shows that, at least for the scenarios considered, the SOHPIE-DNA generally maintains higher recall and F1-score while maintaining similar precision and accuracy, when compared with the most recent state-of-the-art methods: NetCoMi and MDiNE.

In this study, the group-specific jackknife pseudo-values are calculated. Another way of calculating jackknife pseudo-values is to use the entire sample and introduce the group-level indicator as a covariate into the model. However, in our preliminary simulations, we found that doing it that way led to worse performance.

Albeit not reported, we also looked at the familywise error rate (FWER), as defined to be the probability of at least one false positive and the values were fairly high in some cases. However, our simulation results shown in this paper, still reassure the utility of our proposed method since we generally are not expecting the complete null (where none of the edges to be DC) to hold and the FWER is a stringent measure as generally accepted by many statisticians. In our opinion, the reverse engineering methods such as ours should only a used as a screening tool and any positive discovery should be experimentally validated to alleviate such concerns. Incidentally, if FWER control is deemed to be very important for some situations, our tests could be combined with a Westfall-Young type procedure [40]. The detailed performance of such a modification could be studied elsewhere.

Another issue that we encountered was the incorporation of q-values, into our procedure. Since our individual tests are not independent, the q-values may not have the classical properties. Nevertheless, our tests seem to have reasonable FDR values as can be seen from the empirical results (Tables 1, 2 and 3).

We want to highlight that the SOHPIE-DNA is theoretically feasible to accommodate categorical biological groups, in lieu of binary biological groups. To the best of our knowledge, the use of binary groups has been commonly used for the DN analysis. Further, we have presented our simulation and real-data application studies with binary groups only.

We analyzed the data from two published studies to showcase the utility of the SOHPIE-DNA. Firstly, 6 taxa are found to be significantly DC between migraineurs and non-migraineurs while accounting for covariates using the data from the American Gut Project. A slight modification to the proposed method is grafted for analyzing the Diet Exchange Study data, where the group-specific difference of the estimated association matrices between two time points are used for the pseudo-value calculation. As a result, 16 significantly DC taxa are identified between AA and RA after the two-week diet exchange intervention with the inclusion of covariates. The real-world microbiome data often includes hundreds to thousands of taxa. We recommend that the users should (1) focus on a subset of taxa that are chosen based on experts with biological or clinical knowledge or (2) utilize our method at higher taxonomic levels (such as phylum level).

The latter application demonstrates the capability of assessing the temporal variation in connectivity measures. However, the SOHPIE-DNA currently has no feature to address the within-subject correlation for repeated measurements at different time points. This opens up an avenue for future investigation of longitudinal microbiome studies. One way of handling this is to use a generalized estimating equations (GEE) type approach for the pseudo-values and utilizing a jackknife estimate of the variance-covariance matrix of the pseudo-values at different time points.

Another line of future research direction to extend our work is to consider the idea of variable selection. This will help finding the best prediction model with a subset of phenotypic variables that are more biologically relevant across more heterogeneous study samples.

Additionally, we made an attempt of fitting a model under the generalized linear model for binary outcomes: logistic regression with or without the Firth’s correction, in case of small sample size. It was challenging to appropriately dichotomize the matrices with jackknife pseudo-values. Further studies will be needed to devise an adaptive algorithm to find a threshold value that better classify the jackknife pseudo-values.

As a last remark, it should be emphasized that methods other than SparCC were also considered for network estimation, which includes the CCLasso [41] and SPIEC-EASI [42] with graphical lasso or neighborhood selection algorithms. However, these were not favorable in terms of runtime or due to not being able to run under certain simulation scenarios. For instance, the computational time to complete the re-estimation step for the SPIEC-EASI took more than 200 min for $p = 20$ with $n = 200$ for a single simulation replicate. The CCLasso could not estimate the association matrix with small sample size for a smaller network ($p = 20$ for $n = 20, 40, 60$).

Conclusions

There has been limited research to date that discusses how to adjust for additional covariate information in DN analysis for microbiome data. Herewith, we propose SOHPIE-DNA, a novel pseudo-value regression approach for the DN analysis, which can include additional clinical covariate in the model.

Methods

Compositional correlation-based methods for network estimation

The correlation is a useful proxy measure for identifying co-abundances or dependencies among taxa (or OTUs) in a microbial network. The Sparse Correlations for Compositional Data (SparCC) [43] estimates the pairwise correlations of the log-ratio transformed OTU abundances. Of note, a recent method, namely a Pseudo-value Regression Approach for Network Analysis (PRANA) [44], operates on gene expression data only, which therefore does not use a correlation measure that preserves the compositional profiling.

The co-abundance among taxa is described by a covariance matrix $T \in \mathbb {R}^{p \times p}$ where the non-diagonal elements $t_{jk}$ are expressed by

$$\begin{aligned} t_{{jk}} \equiv & {\text{Var}}\left( {\log \frac{{u_{j} }}{{u_{k} }}} \right) \\ = & {\text{Var}}\left( {\log u_{j} } \right) + {\text{Var}}\left( {\log u_{k} } \right) - 2{\text{Cov}}\left( {\log u_{j} ,\log u_{k} } \right) \\ = & \sigma _{j}^{2} + \sigma _{k}^{2} - 2\rho _{{jk}} \sigma _{j} \sigma _{k} , \\ \end{aligned}$$

(1)

where $u_j$ and $u_k$ are the fraction of OTU abundances, $\sigma _{j}^{2}$ and $\sigma _{k}^{2}$ are the variances of the log-transformed abundances, and $\rho _{jk}$ is the correlation of taxa j and k, respectively. Moreover, the variance $t_{jj}$ is approximated by

$$\begin{aligned} t_{jj} \cong (p-1) \sigma _{j}^{2} + \sum _{k \ne j} \sigma _{k}^{2}, \end{aligned}$$

(2)

where $j, k \in \{1, \dots p\}$. Then the correlation can be estimated by solving Eqs. 1 and 2:

$$\begin{aligned} \hat{\rho }_{jk} = \frac{\hat{\sigma }_{j}^{2} + \hat{\sigma }_{k}^{2} - \hat{t}_{jk}}{2 \hat{\sigma }_{j} \hat{\sigma }_{k}}, \end{aligned}$$

(3)

where $\hat{\sigma }_{j}$, $\hat{\sigma }_{k}$, and $\hat{t}_{jk}$ are the sample estimates of $\sigma _{j}$, $\sigma _{k}$, and $t_{jk}$, respectively.

Furthermore, SparCC takes an iterative approach under the assumption (“sparsity of correlations” as in the original paper) that a small number of strong correlations exists in a true network, which hinders the detection of spurious correlations among taxa.

Besides SparCC, we have attempted to use other compositional correlation measures for our differential network analysis. See the “Discussion” section for further details.

Pseudo-value approach

Consider undirected network estimated from n individuals. It can then be represented by the $p \times p$ association matrix that encodes the pairwise correlations $\hat{\rho }_{jk}$ between a pair of taxa $j,k \in \{1, \dots , p\}$. The association matrix is symmetric ($\hat{\rho }_{jk} =$ $\hat{\rho }_{kj}$) where the non-diagonal entries are either non-zero (i.e., some association between two taxa) or zero (i.e., no association between two taxa). The diagonal entries are all equal to one, because the network is assumed that there is no self-loop (i.e., a node cannot redirect to itself).

The network centrality has been studied to measure the extent of biological or topological importance that a node has in a network [45, 46]. For each taxa k, the network centrality is calculated as the marginal sum of the association matrix.

$$\begin{aligned} \hat{\theta }_{k} = \sum _{j=1}^{p} \hat{\rho }_{jk}, \end{aligned}$$

where $k = 1, \dots , p$.

The jackknife pseudo-values [17] for the $i^\text {th}$ individual and $k^\text {th}$ taxon are defined by:

$$\begin{aligned} \tilde{\theta }_{ik} = n\hat{\theta }_{k} - (n-1)\hat{\theta }_{k(i)}, \end{aligned}$$

(4)

where $\hat{\theta }_{k(i)}$ is the marginal sum of a taxon calculated based on the re-estimated association matrix using the microbiome data eliminating the $i^\text {th}$ subject.

The computational cost of the re-estimation process is dependent on the sample size, as for each taxa k requires n such calculations with the data size of $n-1$. A solution to speed up the processing time is the use of parallel computing such as mclapply function in parallel R package.

Let $Z \in \{1, 2\}$ be a binary group indicator and denote $\mathcal {G}_1 = \{i: Z_{i} = 1 \}$ and $\mathcal {G}_2 = \{i: Z_{i} = 2 \}$. Each group has the same set of p taxa, but group-specific sample size $n_{z} = |\mathcal {G}_z |$ for the two groups $z=1, 2$. Total sample size is $n = \sum _z n_{z}$. The Eq. 4 is used to calculate the group-specific jackknife pseudo-values. That is, for taxon k and group z, we define $\hat{\theta }_{k}^z$ and $\hat{\theta }_{k(i)}^z$, where $i = 1,\ldots , n_{z}$. Then for each $i \in \mathcal {G}_z$, the $k^\text {th}$ taxon jackknife pseudo-values are calculated from $\tilde{\theta }_{ik} = n_z\hat{\theta }_{k}^z - (n_z-1)\hat{\theta }_{k(i)}^z$.

Let ${\textbf {X}} = (X_{1}, \dots , X_{q}$) denote q vector of covariates, such as age at diagnosis, current smoking status, and etc. The pseudo-value regression model for the $i^\text {th}$ individual and $k^\text {th}$ taxon is

$$\begin{aligned} \mu _{i} = E[\tilde{\theta }_{ik} \mid Z_{i}, {\textbf {X}}_{i}] = \alpha _{k} + \beta _{k}Z_{i} + \sum _{m=1}^{q} \gamma _{km} X_{im}, \end{aligned}$$

(5)

where $\mu _{i}$ is the k-dimensional mean vector of pseudo-value $\tilde{\theta }_{ik}$ for the i^th individual, $\alpha _{k}$ is the intercept, $\beta _{k}$ is the regression coefficient for Z, and $\gamma _{k1}, \dots , \gamma _{kq}$ is the set of regression coefficients to be estimated for ${\textbf {X}}$. In our setting, the main parameter of interest is given by $\beta _{k}$, the change in network centrality measure of the $k^{\text {th}}$ taxon between two groups.

The least trimmed squares (LTS), also known as least trimmed sum of squares [47], is then implemented to carry out a robust regression. The main advantages of the LTS estimator over other robust estimators including the M-estimator and least median of squares estimator are its computational efficiency and robustness to outliers in both the response and predictor variables [48, 49].

The LTS estimator is defined by

$$\min _{\alpha _k, \beta _k, \gamma _{k1}, \ldots , \gamma _{kq}} \sum _{i=1}^h r_{(i)} (\alpha _k, \beta _k, \gamma _{k1}, \ldots , \gamma _{kq})^2,$$

where $r_{(i)}$ is the set of ordered absolute values of the residuals sorted in increasing order of absolute value and h may depend on a pre-determined trimming proportion $c \in [0.5, 1]$ [50]. For example, one can take $h = [n(1-c)] + 1$.

Hypothesis testing

We construct the null hypothesis of $H_{0}: \beta _{k} = 0$ against the research hypothesis $H_{1}: \beta _{k} \ne 0$ to test if there is a true difference between groups in the network centrality measure of the $k^\text{th}$ taxon. The t-statistic is defined by $U_{k} = \hat{\beta }_{k}/SE(\hat{\beta }_{k})$ for $k = 1, \dots , p$, where $\hat{\beta }_{k}$ is the least-squares estimator from the robust regression described in the above Eq. 5 and $SE(\hat{\beta }_{k})$ is the standard error of $\hat{\beta }_{k}$. As far as the decision-making process, the asymptotically $\alpha$-level test rejects $H_{0}$ if $|U_{k} |> t_{\alpha /2}$. p values are calculated using a t-distribution as in robustbase R package [51, 52].

Multiple hypothesis testing is a common feature in the DN analysis, and therefore it is crucial to appropriately control the false discovery rate (FDR). The FDR measures the proportion of false discoveries incurred among a set of DC taxa from the test. Most classically, the concept of FDR was pioneered by Benjamini and Hochberg [53], shown to achieve the FDR control, whilst maintaining the adequate statistical power [54]. However, the q-value [55] offers a less conservative FDR estimation over the conventional Benjamini-Hochberg procedure [56]. The q-value is estimated from the empirical distribution of the observed p values, and keeps the balance between true positives and false positives [57]. Accordingly, the q-value is applied to adjust for the multiplicity control in the present paper using fdrtool R package.

Algorithm

The SOHPIE-DNA algorithm is described below in Algorithm 1.

Algorithm 1

SOHPIE-DNA

Performance evaluations

Construction of adjacency matrices

Generate the scale-free random network (or Barabási-Albert network) [58] with p nodes using the igraph R package [59]. A network is scale-free if its degree distribution follows a power-law distribution. In other words, a small portion of “hub” nodes has the highest degree centrality, while most nodes have lower degree centrality.

The two identical $p \times p$ adjacency matrices, where the diagonal entries are 0 and non-diagonal entries are either $\{0, 1\}$, are obtained from this random network. At the end of the data generation phase using SparseDOSSA2 in Simulated Data in Materials subsection, we are able to identify which taxa are spike-in associated with the covariate for each $z = 1, 2$. In order to distinguish networks representative of $z = 1$ (e.g., healthy control) from that of $z = 2$ (e.g., disease group), we keep track of the indices of these covariate-dependent taxa. Each index with a value of 1 indicates that the corresponding covariate-dependent taxon is connected with at least one of the neighboring taxa. We use these indices to perturb the random networks by changing a value from 1 to 0 (i.e. synthetically removing all the connected edges) around the covariate-dependent taxa for each group. This perturbation is further explained and borrowed from the recent paper [60]. The network plots are provided to visually demonstrate the perturbed adjacency matrices (see Figs. 1 and 2). The figures represent two single networks for two particular realizations corresponding to covariate profiles. The effect sizes (i.e. pre-specified proportion of taxa that are associated with the covariate; denoted as $\delta _{1}$ and $\delta _{2}$) control the amount of perturbation. If the effect sizes are different ($\delta _{1} \ne \delta _{2}$), then the covariates are affecting the networks differentially (Fig. 2). See Simulated Data in Materials subsection for further details.

Performance measures

Four performance metrics are adopted to evaluate our proposed method: precision, recall, F1-score, and accuracy. Let $\Omega ^{z} \in \mathbb {R}^{p \times p}$ be the group-specific adjacency matrix, where

$$\begin{aligned} \Omega _{jk}^{z} = {\left\{ \begin{array}{ll} 1 &{} \text {if the two nodes } j \text { and } k \text { are connected}\\ 0 &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$

for $z = 1, 2$. Next, a node-specific true connection is calculated

$$\begin{aligned} \eta _{k} = I\bigg ( \sum _{j=1}^{p} |\Omega _{jk}^{1} - \Omega _{jk}^{2} |> 0 \bigg ), \end{aligned}$$

indicating that taxa k has differential connectivity (DC).

In terms of notation, we use $q_{ks}$ to denote a q-value [55] of taxa k at the simulation replicate s. An error rate control of $\alpha = 0.05$ is used throughout the simulation. In the following, we present the details of each performance metric.

Precision is the fraction of taxa which are declared to be significantly DC from the test that are confirmed as true:

$$\begin{aligned} \text {Precision} = \frac{ \sum _{k = 1}^{p} \eta _{k} \, I(q_{ks}< \alpha ) }{\sum _{k = 1}^{p} I(q_{ks} < \alpha ) }. \end{aligned}$$

Recall is the fraction of truly DC taxa which are correctly declared to be significant between two comparing groups from the test:

$$\begin{aligned} \text {Recall} = \frac{ \sum _{k = 1}^{p} \eta _{k} \, I(q_{ks} < \alpha ) }{\sum _{k = 1}^{p} \eta _{k} }. \end{aligned}$$

The F1 score is the harmonic mean of precision and recall values. A higher F1 score indicates a better overall performance with lower false negative and false positive predictions:

$$\begin{aligned} \text {F1} = 2 \cdot \frac{\text {Precision} \cdot \text {Recall}}{\text {Precision} + \text {Recall}}. \end{aligned}$$

Accuracy is defined as the fraction of total number of taxa that are correctly predicted to be DC. The accuracy ranges from 0 (no correct predictions) to 1 (perfect predictions):

$$\begin{aligned} \text {Accuracy} = \frac{ \sum _{k = 1}^{p} I(\eta _{k}=1) \, I(q_{ks} < \alpha ) + \sum _{k = 1}^{p} I(\eta _{k}=0) \, I(q_{ks} \ge \alpha ) }{\sum _{k = 1}^{p} I(\eta _{k}=1) + \sum _{k = 1}^{p} I(\eta _{k}=0) }. \end{aligned}$$

Materials

Simulated data

The synthetic microbiome dataset are structured with p taxa and n sample size. In the simulation, binary group indicators 1 and 2 are generated from a Bernoulli distribution with equal probabilities and a single continuous covariate $X \sim N(55, 10)$ (e.g., age at diagnosis). We test our proposed method on datasets under two different simulation scenarios: taxa are impacted by the effect of (1) Z and X or (2) Z only, which each corresponds to “multivariable” and “univariable” settings, respectively.

The actual microbial data generation (e.g., OTU counts) given the covariates is described next. In this context, it is perhaps worth mentioning that this part is completely different from generating gene expression data as in PRANA [44]. For each simulation scenario, we generate an OTU table that resembles the dependence structure of covariates Z and/or X on the microbial community (or the network) using the SparseDOSSA2 (Sparse Data Observations for the Simulation of Synthetic Abundances) R package [61]. SparseDOSSA2 adopts a Bayesian Gaussian copula model with zero-inflated, truncated log-normal distributions to capture the marginal distributions of each microbial taxa and to account for the correlation between taxa.

The package has a feature to indicate a user-specified percentage of taxa to be “spiked-in” association with the clinical information (or metadata). This is referred to as the “effect size” of differential abundance $\delta$. To evaluate the effect size of Z under the univariable setting, we generate the data that half of the samples have taxa with no spike-in association, whereas the other half of the samples have spike-in association on $5\%$, $10\%$, or $20\%$ of taxa. The distributions of age in the two groups are different. Therefore, under the multivariable setting, $5\%$, $10\%$, or $20\%$ of taxa have spike-in association with X for each group $z = 1, 2$. In both scenarios, $n_{z} \times p$ matrices for each group $z = 1, 2$ will be available for use.

Application study

The American Gut Project Data A pre-processed OTU table of the human stool microbiome samples from the American Gut Project [33] is available in the SpiecEasi R package, along with the corresponding metadata information. The gut microbiome is involved with the bidirectional relationship between the gastrointestinal system and central nervous system (i.e. gut-brain axis) that impacts on the migraine inflammation [62].

In the analysis, the main variable of interest is a binary variable indicating the migraine headache (yes or no). Age [63], sex [63], exercise frequency ($\ge 3$ days per week or otherwise) [64], and categorical alcohol consumption (heavy, moderate, or non-drinking) [65] are covariates that are included in the multivariable model. Additionally, migraine has been associated with the periodontal inflammation [66] and pet ownership [67], and therefore the oral hygiene behavior such as dental floss frequency ($\ge 3$ times per week or otherwise) and living with a dog (yes or no) were included in the model.

The initial OTU table consists of 138 taxa with 296 subjects. No taxa were removed, however, 28 subjects were excluded due to unidentified sampling body site and missing age or sex information. Hence, 138 taxa and 268 subjects were used for the analysis.

The Diet Exchange Study Data A pre-processed data of the geographical epidemiology study [34] is available in microbiome [68] R package. The aim of the study was to assess the effect of fat and fiber intake of the diet on the composition of the colonic microbiota by switching the diet in study populations with high (African-Americans from Pittsburgh area of Pennsylvania; AA) and low (rural South Africans from KwaZulu region; RA) colon cancer risk for two weeks.

An initial OTU table contains 130 taxa with 38 subjects. After the exclusion of a subject with missing post-dietary intervention data and 18 rare taxa that appear in fewer than 10$\%$ of the samples, 112 taxa with 37 subjects (20 AA and 17 RA) are used for the analysis.

The main predictor variable is binary geographic location (AA or RA). Additional covariates considered in a multivariable model were sex and BMI groups (obese, overweight, or lean).

For each groups separately, we take the difference of the estimated association matrices (as well as the re-estimated association matrices) between two time points, that is, the endoscopy before and after two weeks of dietary change. The differences are then used to calculate the jackknife pseudo-values as in the previous sections. This additional step is intended to incorporate the temporal change of connectivity of each taxa after dietary interventions.

Fig. 1

Network plots visualizing the microbial network ($p = 20$) with a covariate dependence structure that depends on continuous age and binary group information ($\delta _{1} = 0.05$ (left), $\delta _{2} = 0.2$ (right)). This represents the multivariable setting

Fig. 2

Network plots visualizing the microbial network ($p = 20$) without a covariate dependence structure that depends on binary group only ($\delta _{1} = 0$ (left), $\delta _{2} = 0.2$ (right)). This represents the univariable setting

Table 1 The simulation results for the case when the network structure depends on age covariate

Table 2 The simulation results for the case when the network structure depends on age covariate

Table 3 The simulation results for the case when the network structure does not depend on age covariate