 Research
 Open access
 Published:
Differential network connectivity analysis for microbiome data adjusted for clinical covariates using jackknife pseudovalues
BMC Bioinformatics volume 25, Article number: 117 (2024)
Abstract
Background
A recent breakthrough in differential network (DN) analysis of microbiome data has been realized with the advent of nextgeneration sequencing technologies. The DN analysis disentangles the microbial coabundance 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 Pseudovalue Information and Estimation for Differential Network Analysis (SOHPIEDNA) that incorporates additional covariates such as continuous age and categorical BMI. SOHPIEDNA is a regression technique adopting jackknife pseudovalues that can be implemented readily for the analysis. We demonstrate through simulations that SOHPIEDNA consistently reaches higher recall and F1score, while maintaining similar precision and accuracy to existing methods (NetCoMi and MDiNE). Lastly, we apply SOHPIEDNA 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 SOHPIEDNA 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
SOHPIEDNA 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.Rproject.org/package=SOHPIE), named SOHPIE (pronounced as Sofie). The source code and user manual can be found at https://github.com/sjahnn/SOHPIEDNA.
Background
The human microbiome is the collective genomes of microbes or microorganisms 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 wellbeing [5] and mental health including depression [6], autism spectrum disorders [7], and human brain diseases [8].
Following the advent of nextgeneration 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 similaritybased 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 coabundances [12]. The abundance matrix or the observed OTU table is used to infer microbial coabundances among taxa through either correlationbased 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 pseudovalues, a component to calculate the biascorrected estimator of leaveoneout jackknife resampling procedure [17]. The pseudovalue technique was first postulated by Andersen and his colleagues [18, 19] in the context of multistate survival models with rightcensored data. Since then, it has been well studied in various disciplines of statistics including the intervalcensored 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 pseudovalues given covariates are shown with the secondorder von Mises expansion [27, 28]. The pseudovalues 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 quasilikelihood generalized linear mixed model [30] and generalized estimating equations framework [23, 31] for pseudovalue regression approach.
Hence, we propose a regression modeling method for differential connectivity (DC) analysis that regresses the jackknife pseudovalues 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 taxontaxon coabundance 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 pseudovalues. 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 Pseudovalue Information and Estimation for Differential Network Analysis (SOHPIEDNA) 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 SOHPIEDNA 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 covariatedependent network, we compared the pseudovalue regression approach with the recent methods available (NetCoMi and MDiNE) that cannot incorporate the additional covariate. Results show that the SOHPIEDNA 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 SOHPIEDNA 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 SOHPIEDNA 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 SOHPIEDNA 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. nonmigraineurs while adjusting for age, sex, exercise frequency, categorical alcohol consumption, oral hygiene behavior, and dog ownership. At the familylevel, 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 highperformance 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 twoweek 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 SOHPIEDNA, a pseudovalue 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 SOHPIEDNA generally maintains higher recall and F1score while maintaining similar precision and accuracy, when compared with the most recent stateoftheart methods: NetCoMi and MDiNE.
In this study, the groupspecific jackknife pseudovalues are calculated. Another way of calculating jackknife pseudovalues is to use the entire sample and introduce the grouplevel 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 WestfallYoung type procedure [40]. The detailed performance of such a modification could be studied elsewhere.
Another issue that we encountered was the incorporation of qvalues, into our procedure. Since our individual tests are not independent, the qvalues 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 SOHPIEDNA 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 realdata application studies with binary groups only.
We analyzed the data from two published studies to showcase the utility of the SOHPIEDNA. Firstly, 6 taxa are found to be significantly DC between migraineurs and nonmigraineurs 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 groupspecific difference of the estimated association matrices between two time points are used for the pseudovalue calculation. As a result, 16 significantly DC taxa are identified between AA and RA after the twoweek diet exchange intervention with the inclusion of covariates. The realworld 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 SOHPIEDNA currently has no feature to address the withinsubject 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 pseudovalues and utilizing a jackknife estimate of the variancecovariance matrix of the pseudovalues 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 pseudovalues. Further studies will be needed to devise an adaptive algorithm to find a threshold value that better classify the jackknife pseudovalues.
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 SPIECEASI [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 reestimation step for the SPIECEASI 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 SOHPIEDNA, a novel pseudovalue regression approach for the DN analysis, which can include additional clinical covariate in the model.
Methods
Compositional correlationbased methods for network estimation
The correlation is a useful proxy measure for identifying coabundances or dependencies among taxa (or OTUs) in a microbial network. The Sparse Correlations for Compositional Data (SparCC) [43] estimates the pairwise correlations of the logratio transformed OTU abundances. Of note, a recent method, namely a Pseudovalue 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 coabundance among taxa is described by a covariance matrix \(T \in \mathbb {R}^{p \times p}\) where the nondiagonal elements \(t_{jk}\) are expressed by
where \(u_j\) and \(u_k\) are the fraction of OTU abundances, \(\sigma _{j}^{2}\) and \(\sigma _{k}^{2}\) are the variances of the logtransformed abundances, and \(\rho _{jk}\) is the correlation of taxa j and k, respectively. Moreover, the variance \(t_{jj}\) is approximated by
where \(j, k \in \{1, \dots p\}\). Then the correlation can be estimated by solving Eqs. 1 and 2:
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.
Pseudovalue 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 nondiagonal entries are either nonzero (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 selfloop (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.
where \(k = 1, \dots , p\).
The jackknife pseudovalues [17] for the \(i^\text {th}\) individual and \(k^\text {th}\) taxon are defined by:
where \(\hat{\theta }_{k(i)}\) is the marginal sum of a taxon calculated based on the reestimated association matrix using the microbiome data eliminating the \(i^\text {th}\) subject.
The computational cost of the reestimation process is dependent on the sample size, as for each taxa k requires n such calculations with the data size of \(n1\). 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 groupspecific 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 groupspecific jackknife pseudovalues. 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 pseudovalues are calculated from \(\tilde{\theta }_{ik} = n_z\hat{\theta }_{k}^z  (n_z1)\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 pseudovalue regression model for the \(i^\text {th}\) individual and \(k^\text {th}\) taxon is
where \(\mu _{i}\) is the kdimensional mean vector of pseudovalue \(\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 Mestimator 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
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 predetermined trimming proportion \(c \in [0.5, 1]\) [50]. For example, one can take \(h = [n(1c)] + 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 tstatistic is defined by \(U_{k} = \hat{\beta }_{k}/SE(\hat{\beta }_{k})\) for \(k = 1, \dots , p\), where \(\hat{\beta }_{k}\) is the leastsquares 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 decisionmaking process, the asymptotically \(\alpha\)level test rejects \(H_{0}\) if \(U_{k} > t_{\alpha /2}\). p values are calculated using a tdistribution 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 qvalue [55] offers a less conservative FDR estimation over the conventional BenjaminiHochberg procedure [56]. The qvalue is estimated from the empirical distribution of the observed p values, and keeps the balance between true positives and false positives [57]. Accordingly, the qvalue is applied to adjust for the multiplicity control in the present paper using fdrtool R package.
Algorithm
The SOHPIEDNA algorithm is described below in Algorithm 1.
Performance evaluations
Construction of adjacency matrices
Generate the scalefree random network (or BarabásiAlbert network) [58] with p nodes using the igraph R package [59]. A network is scalefree if its degree distribution follows a powerlaw 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 nondiagonal 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 spikein 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 covariatedependent taxa. Each index with a value of 1 indicates that the corresponding covariatedependent 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 covariatedependent 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. prespecified 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, F1score, and accuracy. Let \(\Omega ^{z} \in \mathbb {R}^{p \times p}\) be the groupspecific adjacency matrix, where
for \(z = 1, 2\). Next, a nodespecific true connection is calculated
indicating that taxa k has differential connectivity (DC).
In terms of notation, we use \(q_{ks}\) to denote a qvalue [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:
Recall is the fraction of truly DC taxa which are correctly declared to be significant between two comparing groups from the test:
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:
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):
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 zeroinflated, truncated lognormal 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 userspecified percentage of taxa to be “spikedin” 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 spikein association, whereas the other half of the samples have spikein 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 spikein 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 preprocessed 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. gutbrain 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 nondrinking) [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 preprocessed 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 (AfricanAmericans 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 postdietary 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 reestimated 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 pseudovalues as in the previous sections. This additional step is intended to incorporate the temporal change of connectivity of each taxa after dietary interventions.
Availability of data and materials
The original study data of the American Gut Project and the Diet Exchange Study are available in SpiecEasi (https://github.com/zdk123/SpiecEasi) and microbiome (https://bioconductor.org/packages/release/bioc/html/microbiome.html) R packages, respectively. The SOHPIE [69] R package can be downloaded from the CRAN repository https://CRAN.Rproject.org/package=SOHPIE. The source code for SOHPIEDNA is available at (https://github.com/sjahnn/SOHPIEDNA). Please reach out to the authors (Seungjun Ahn, seungjun.ahn@mountsinai.org; Somnath Datta, somnath.datta@ufl.edu) if you have any further inquiries on the data or code.
Abbreviations
 AA:

AfricanAmericans from Pittsburgh area of Pennsylvania
 RA:

Rural South Africans from KwaZulu region
 IBD:

Inflammatory bowel disease
 OTU:

Operational taxonomic units
 DC:

Differentially connected
 DE:

Differential expression
 DN:

Differential network
 MDiNE:

Microbiome differential network estimation
 NetCoMi:

Network construction and comparison for microbiome data
 SOHPIEDNA:

Statistical approach via pseudovalue information and estimation for differential network analysis
 SparCC:

Sparse correlations for compositional data
 SparseDOSSA2:

Sparse data observations for the simulation of synthetic abundances
 CCLasso:

Correlation inference for compositional data through lasso
 SPIECEASI:

Sparse inverse covariance estimation for ecological association inference
 LASSO:

Least absolute shrinkage and selection operator
 LTS:

Least trimmed squares
References
Weinstock G. Genomic approaches to studying the human microbiota. Nature. 2012;489:250–6.
Bhatt AP, Redinbo MR, Bultman SJ. The role of the microbiome in cancer development and therapy. CA Cancer J Clin. 2017;67:326–44.
VujkovicCvijin I, Sortino O, Verheij E, Sklar J, Wit FW, Kootstra NA, et al. HIVassociated gut dysbiosis is independent of sexual practice and correlates with noncommunicable diseases. Nat Commun. 2020;11:2448.
Glassner KL, Abraham BP, Quigley E. The microbiome and inflammatory bowel disease. J Allergy Clin Immunol. 2020;145:16–27.
Lee SH, Yoon SH, Jung Y, Kim N, Min U, Chun J, et al. Emotional wellbeing and gut microbiome profiles by enterotype. Sci Rep. 2020;10:20736.
VallesColomer M, Falony G, Darzi Y, Tigchelaar EF, Wang J, Tito RY, et al. The neuroactive potential of the human gut microbiota in quality of life and depression. Nat Microbiol. 2019;4:623–32.
KrajmalnikBrown R, Lozupone C, Kang DW, Adams JB. Gut bacteria in children with autism spectrum disorders: challenges and promise of studying how a complex community influences a complex disease. Microb Ecol Health Dis. 2015;26:26914.
Mayer EA, Knight R, Mazmanian SK, Cryan JF, Tillisch K. Gut microbes and the brain: paradigm shift in neuroscience. J Neurosci. 2014;34:
Durazzi F, Sala C, Castellani G, Manfreda G, Remondini D, De Cesare A. Comparison between 16S rRNA and shotgun sequencing data for the taxonomic characterization of the gut microbiota. Sci Rep. 2021;11:3030.
Johnson JS, Spakowicz DJ, Hong BY, Petersen LM, Demkowicz P, Chen L, et al. Evaluation of 16S rRNA gene sequencing for species and strainlevel microbiome analysis. Nat Commun. 2019;10:5029.
Reuter JA, Spacek DV, Snyder MP. Highthroughput sequencing technologies. Mol Cell. 2015;58:586–97.
Layeghifard M, Hwang DM, Guttman DS. Disentangling interactions in the microbiome: A network perspective. Trends Microbiol. 2017;25:217–28.
Matchado MS, Lauber M, Reitmeier S, Kacprowski T, Baumbach J, Haller D, et al. Network analysis methods for studying microbial communities: A mini review. Comput Struct Biotechnol J. 2021;9:2687–98.
McGregor K, Labbe A, Greenwood C. MDiNE: a model to estimate differential cooccurrence networks in microbiome studies. Bioinformatics. 2020;36:1840–7.
Peschel S, Muller C, von Mutius E, Boulesteix A, Depner M. NetCoMi: network construction and comparison for microbiome data in R. Brief Bioinform. 2021;22:290.
Lee M, Chang E. Inflammatory bowel diseases (IBD) and the microbiomesearching the crime scene for clues. Gastroenterology. 2021;160:524–37.
Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Philadelphia: Chapman & Hall/CRC; 1993.
Andersen PK, Klein JP, Rosthøj S. Generalised linear models for correlated pseudoobservations, with applications to multistate models. Biometrika. 2003;90:15–27.
Andersen P, Klein J. Regression analysis for multistate models based on a pseudovalue approach, with applications to bone marrow transplantation studies. Scand J Statist. 2007;34:3–16.
Sabathé C, Andersen PK, Helmer C, Gerds TA, JacqminGadda H, Joly P. Regression analysis in an illnessdeath model with intervalcensored data: A pseudovalue approach. Stat Methods Med Res. 2020;29:752–64.
Johansen MN, LundbyeChristensen S, Larsen JM, Parner ET. Regression models for interval censored data using parametric pseudoobservations. BMC Med Res Methodol. 2021;21:36.
Logan BR, Zhang MJ, Klein JP. Marginal models for clustered timetoevent data with competing risks using pseudovalues. Biometrics. 2011;67:1–7.
Ahn KW, Logan BR. Pseudovalue approach for conditional quantile residual lifetime analysis for clustered survival and competing risks data with applications to bone marrow transplant data. Ann Appl Stat. 2016;10:618–37.
Zhao L, Feng D. Deep neural networks for survival analysis using pseudo values. IEEE J Biomed Health Inform. 2020;24:3308–14.
Ginestet PG, Gabriel EE, Sachs MC. Survival stacking with multiple data types using pseudoobservationbasedAUC loss. J Biopharm Stat. 2022. https://doi.org/10.1080/10543406.2022.2041655.
Logan BR, Klein JP, Zhang MJ. Comparing treatments in the presence of crossing survival curves: an application to bone marrow transplantation. Biometrics. 2008;64:733–40.
Graw F, Gerds TA, Schumacher M. On pseudovalues for regression analysis in competing risks models. Lifetime Data Anal. 2009;15:241–55.
Overgaard M, Parner ET, Pedersen J. Asymptotic theory of generalized estimating equations based on jackknife pseudoobservations. Ann Stat. 2017;45:1988–2015.
Klein JP, Gerster M, Andersen PK, Tarima S, Perme MP. SAS and R functions to compute pseudovalues for censored data regression. Comput Methods Programs Biomed. 2008;89:289–300.
Wang Y, Logan B. Testing for center effects on survival and competing risks outcomes using pseudovalue regression. Lifetime Data Anal. 2019;25:206–28.
Ahn K, Mendolia F. Pseudovalue approach for comparing survival medians for dependent data. Stat Med. 2014;33:1531–8.
Zhao S, Shojaie A. Network differential connectivity analysis. Ann Appl Stat. 2022;16:2166–82.
McDonald D, Hyde E, Debelius JW, Morton JT, Gonzalez A, Ackermann G, et al. American gut: an open platform for citizen science microbiome research. mSystems. 2018;3:00031–18.
O’Keefe SJ, Li JV, Lahti L, Ou J, Carbonero F, Mohammed K, et al. Fat, fibre and cancer risk in african americans and rural africans. Nat Commun. 2015;6:6342.
Taur Y, Pamer E. Harnessing microbiota to kill a pathogen: Fixing the microbiota to treat clostridium difficile infections. Nat Med. 2014;20:246–7.
NolanKenney R, Wu F, Hu J, Yang L, Kelly D, Li H, et al. The association between smoking and gut microbiome in bangladesh. Nicotine Tob Res. 2020;22:1339–46.
Chen J, Wang Q, Wang A, Lin Z. Structural and functional characterization of the gut microbiota in elderly women with migraine. Front Cell Infect Microbiol. 2020;9:470.
Nie K, Ma K, Luo W, Shen Z, Yang Z, Xiao M, et al. Roseburia intestinalis: A beneficial gut organism from the discoveries in genus and species. Front Cell Infect Microbiol. 2021;11: 757718.
Hajjar J, Mendoza T, Zhang L, Fu S, PihaPaul SA, Hong DS, et al. Associations between the gut microbiome and fatigue in cancer patients. Sci Rep. 2021;11:5847.
Westfall P, Young SS. ResamplingBased Multiple Testing: Examples and Methods for pValue Adjustment. New York: Wiley; 1993.
Fang H, Huang C, Zhao H, Deng M. CCLasso: correlation inference for compositional data through lasso. Bioinformatics. 2015;31:3172–80.
Kurtz ZD, Muller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol. 2015;11:1004226.
Friedman J, Alm E. Inferring correlation networks from genomic survey data. PLoS Comput Biol. 2012;8:1002687.
Ahn S, Grimes T, Datta S. A pseudovalue regression approach for differential network analysis of coexpression data. BMC Bioinformatics. 2023;24:8.
Ashtiani M, SalehzadehYazdi A, RazaghiMoghadam Z, Hennig H, Wolkenhauer O, Mirzaie M, et al. A systematic survey of centrality measures for proteinprotein interaction networks. BMC Syst Biol. 2018;12:80.
Ozgür A, Vu T, Erkan G, Radev DR. Identifying genedisease associations using centrality on a literature mined geneinteraction network. Bioinformatics. 2008;24:277–85.
Rousseeuw P. Least median of squares regression. J Am Stat Assoc. 1984;79:871–80.
Ahdesmäki M, Lähdesmäki H, Gracey A, Shmulevich L, YliHarja O. Robust regression for periodicity detection in nonuniformly sampled timecourse gene expression data. BMC Bioinformatics. 2007;8:233.
Alfons A, Croux C, Gelper S. Sparse least trimmed squares regression for analyzing highdimensional large data sets. Ann Appl Stat. 2013;7:226–48.
Pison G, Van Aelst S, Willems G. Small sample corrections for LTS and MCD. Metrika. 2002;55:111–23.
Todorov V, Filzmoser P. An objectoriented framework for robust multivariate analysis. J Stat Soft. 2009;32:1–47.
Maechler M, Rousseeuw P, Croux C, Todorov V, Ruckstuhl A, SalibianBarrera M, et al.: Robustbase: Basic Robust Statistics. (2022). R package version 0.950. http://robustbase.rforge.rproject.org/
Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Statist Soc B. 1995;57:289–300.
Benjamini Y. Discovering the false discovery rate. J R Statist Soc B. 2010;72:405–16.
Storey J. A direct approach to false discovery rates. J R Statist Soc B. 2002;64:479–98.
Strimmer K. A unified approach to false discovery rate estimation. BMC Bioinformatics. 2008;9:303.
Storey J, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A. 2003;100:9440–5.
Barabási A, Albert R. Emergence of scaling in random networks. Science. 1999;286:509–12.
Csárdi G, Nepusz T. The igraph software package for complex network research. InterJournal, Complex Systems. 2006;1695:16–27.
Grimes T, Datta S. SeqNet: An R Package for Generating GeneGene Networks and Simulating RNASeq Data. J Stat Softw. 2021;98:10–1863709812.
Ma S, Ren B, Mallick H, Moon YS, Schwager E, Maharjan S, et al. A statistical model for describing and simulating microbial community profiles. PLoS Comput Biol. 2021;17:1008913.
Arzani M, Jahromi SR, Ghorbani Z, Vahabizad F, Martelletti P, Ghaemi A, et al. Gutbrain axis and migraine headache: a comprehensive review. J Headache Pain. 2020;21:1.
Stewart WF, Linet MS, Celentano DD, Van Natta M, Ziegler D, et al. Age and sexspecific incidence rates of migraine with and without visual aura. Am J Epidemiol. 1991;134:1111–20.
Amin FM, Aristeidou S, Baraldi C, CzapinskaCiepiela EK, Ariadni DD, Di Lenola D, et al. The association between migraine and physical exercise. J Headache Pain. 2018;19:83.
Mostofsky E, Bertisch SM, Vgontzas A, Buettner C, Li W, Rueschman M, et al. Prospective cohort study of daily alcoholic beverage intake as a potential trigger of headaches among adults with episodic migraine. Ann Med. 2020;52:386–92.
Leira Y, Ameijeira P, Domínguez C, LópezArias E, ÁvilaGómez P, PérezMato M, et al. Periodontal inflammation is related to increased serum calcitonin generelated peptide levels in patients with chronic migraine. J Periodontol. 2019;90:1088–95.
Koivusilta L, Ojanlatva A. To have or not to have a pet for better health? PLoS One. 2006;1:109.
Lahti L, Shetty S. Microbiome R Package. (2017). Bioconductor. https://doi.org/10.18129/B9.bioc.microbiome
Ahn S, Datta S. SOHPIE: statistical approach via pseudovalue information and estimation for differential network analysis of microbiome data. Bioinformatics. 2024;40(1):btad766. https://doi.org/10.1093/bioinformatics/btad766
Acknowledgements
The authors are exceptionally thankful for the investigators involved with the American Gut Project and the Diet Exchange Study for sharing their data publicly.
Funding
SA was supported by the National Institute on Alcohol Abuse and Alcoholism of the National Institutes of Health under grants number [NIH T32AA025877] during his time at the University of Florida for his doctoral research. Research reported in this publication was supported in part by the National Cancer Institute Cancer Center Support Grant [NIH P30CA19652101] awarded to the Tisch Cancer Institute of the Icahn School of Medicine at Mount Sinai and used the Biostatistics Shared Resource Facility. The content is solely the responsibility of S.A. and does not necessarily represent the official views of the National Institutes of Health.
Author information
Authors and Affiliations
Contributions
SD conceived the use of pseudovalue in the study. SA developed the methodology, performed simulations and data analyses of the study. SA drafted the manuscript and SD provided suggestions when writing the manuscript. All authors have reviewed and edited the manuscript.
Corresponding author
Ethics declarations
Ethical approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests. S.D. is a member of editorial board (Associate Editor), but this did not influence the scientific peerreview process nor did he discuss the submission with other editorial members of the journal.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1
. Comparison of computational time of SOHPIEDNA with that of NetCoMi and MDiNE.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Ahn, S., Datta, S. Differential network connectivity analysis for microbiome data adjusted for clinical covariates using jackknife pseudovalues. BMC Bioinformatics 25, 117 (2024). https://doi.org/10.1186/s12859024056897
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12859024056897