Data processing and transformation of standard OTU count data

For this discussion, a table of OTU count data, typical output of 16S rRNA gene sequencing data curation pipelines (e.g., mothur [32], QIIME [33]) are given. The OTU data are stored in a matrix $\mathrm{\text{W}}\in {\mathbb{N}}_0^{n\times p}$ where $w^{(j)}=[w_1^{(j)},w_2^{(j)},\dots ,w_p^{(j)}]$ denotes the p-dimensional row vector of OTU counts from the j^th sample, j = 1, …, n, with total cumulative count $m^{(j)}=\sum _{i=1}^p{w}_i^{(j)}$; ℕ₀ denotes the set of natural numbers {0, 1, 2, …}. As described above, to account for sampling biases, microbiome data is typically transformed by normalizing the raw count data w^(j) with respect to the total count m^(j) of the sample [10]. We thus arrive at vectors of relative abundances or compositions $x^{(j)}=[\frac{w_1^{(j)}}{m^{(j)}},\frac{w_2^{(j)}}{m^{(j)}},\dots ,\frac{w_p^{(j)}}{m^{(j)}}]$ for sample j. Due to this normalization OTU abundances are no longer independent, and the sample space of this p-part composition x^(j) is not the unconstrained Euclidean space but the p-dimensional unit simplex ${𝕊}^p\doteq \{x\mid x_i>0,{\sum }_{i=1}^p{x}_i=1\}$. Thus, OTU compositions from n samples are constrained to lie in the unit simplex, X ∈ 𝕊^n×p. This restriction of the data to the simplex prohibits the application of standard statistical analysis techniques, such as linear regression or empirical covariance estimation. Covariance matrices of compositional data exhibit, for instance, a negative bias due to closure effects.

Major advances in the statistical analysis of compositional data were achieved by Aitchison in the 1980’s [18, 34]. Rather than considering compositions in the simplex, Aitchison proposed log-ratios, $\log [\frac{x_i}{x_j}]$, as a basis for studying compositional data. The simple equivalence $\log [\frac{x_i}{x_j}]=\log [\frac{w_i/m}{w_j/m}]=\log [\frac{w_i}{w_j}]$ implies that statistical inferences drawn from analysis of log-ratios of compositions are equivalent to those that could be drawn from the log-ratios of the unobserved absolute measurements, also termed the basis.

Aitchison also proposed several statistically equivalent log-ratio transformations to remove the unit-sum constraint of compositional data [18]. Here we apply the centered log-ratio (clr) transform: (1) where $g(x)={\left[\underset{i=1}{\prod ^p}{x}_i\right]}^{1/p}$ is the geometric mean of the composition vector. The clr transform is symmetric and isometric with respect to the component parts. The resulting vector z is constrained to a zero sum. The clr transform maps the data from the unit simplex to a p − 1-dimensional Euclidean space, and the corresponding population covariance matrix Γ = Cov[clr(X)] ∈ ℝ^p×p is also singular [18]. The covariance matrix Γ is related to the population covariance of the log-transformed absolute abundances Ω = Cov[logW] via the relationship [34]: (2) where $\mathrm{\text{G}}={\mathrm{\text{I}}}_p-\frac{1}{p}\mathrm{\text{J}}$, I_p is the p-dimensional identity matrix, and J = [j₁, j₂, …, j_i, …, j_p], j_i = [1, 1, …, 1] the p-dimensional all-ones vector. For high-dimensional data, p > > 0, the matrix G is close to the identity matrix, and thus we can assume that a finite sample estimator $\widehat{\Gamma }$ of Γ serves as a good approximation of $\widehat{\Omega }$. This approximation serves as the basis of our network inference scheme. Finally, because real-world OTU data often contain samples with a zero count for low-abundance OTUs, we add a unit pseudo count to the original count data to avoid numerical problems with the clr transform.

Inference of microbial associations from microbial abundance datasets

Our key objective is to learn a network of pairwise taxon-taxon associations (putative interactions) from clr-transformed microbiome compositions Z ∈ ℝ^n×p. We represent the network as an undirected, weighted graph 𝓖 = (V, E), where the vertex set V = {v₁, …, v_p} represents the p taxa (e.g., OTUs) and the edge set E ⊂ V × V the possible associations among them. Our formal approach is to learn a probabilistic graphical model [35] (i) that is consistent with the observed data and (ii) for which the (unknown) graph 𝓖 encodes the conditional dependence structure between the random variables (in our case, the observed taxa). Graphical models over undirected graphs (also known as Markov networks or Markov Random Fields) have a straightforward distributional interpretation when the data are drawn from a probability distribution π(x) that belongs to an exponential family [36, 37]. For example, when the data are drawn from a multivariate normal distribution π(x) = 𝓝(x∣μ, Σ) with mean μ and covariance Σ, the non-zero elements of the off-diagonal entries of the inverse covariance matrix Θ = Σ⁻¹, also termed the precision matrix, defines the adjacency matrix of the graph 𝓖 and thus describes the factorization of the normal distribution into conditionally dependent components [35]. Conversely, if and only if an entry in Θ: Θ_{i, j} = 0, then the two variables are conditionally independent, and there is no edge between v_i and v_j in 𝓖. We seek to estimate the inverse covariance matrix from the data, thereby inferring associations based on conditional independence. This is fundamentally distinct from SparCC and CCREPE (see S1 Table), which essentially estimate pairwise correlations (though other pairwise metrics could be considered for CCREPE). We highlight this key difference in Fig 1. For an intuitive introduction to graphical models in the context of biological networks see Bühlmann et. al, 2014 [38].

Inferring the exact underlying graph structure in the presence of a finite amount of samples is, in general, intractable. However, two types of statistical inference procedures have been useful in high-dimensional statistics due to their provable performance guarantees under assumptions about the sample size n, dimensionality p, underlying graph properties, and the generating distribution [29, 39]. The first approach, neighborhood selection [20, 39], aims at reconstructing the graph on a node-by-node basis where, for each node, a penalized regression problem is solved. The second approach is the penalized maximum likelihood method [22, 23], where the entire graph is reconstructed by solving a global optimization problem, the so-called covariance selection problem [40]. The key advantages of these approaches are that (i) their underlying inference procedures can be formulated as convex (and hence tractable) optimization problems, and (ii) they are applicable even in the underdetermined regime (p > n), provided that certain structural assumptions about the underlying graph hold. One assumption is that the true underlying graph is reasonably sparse, e.g., that the number of taxon-taxon associations scales linearly with the number of measured taxa.

Graphical model inference. The SPIEC-EASI pipeline comprises two types of inference schemes, neighborhood and covariance selection. The neighborhood selection framework, introduced by Meinshausen and Bühlmann [20] and thus often referred as the MB method, tackles graph inference by solving p regularized linear regression problems, leading to local conditional independence structure predictions for each node. Let us denote the i^th column of the data matrix Z by Zⁱ ∈ ℝⁿ and the remaining columns by Z^¬i ∈ ℝ^n×p−1. For each node v_i, we solve the following convex problem: (3) where $\Vert a{\Vert }_1={\sum }_{i=1}^{p-1}\mid a_i\mid$ denotes the L1 norm, and λ ≥ 0 is a scalar tuning parameter. This so-called LASSO problem [41] aims at balancing the least-square fit and the number of necessary predictors (the non-zero components β_j of β) by tuning the λ parameter. We define the local neighborhood of a node v_i as $N_i^{\lambda }=\{j\subset \{1,\dots p\}\backslash i:{\widehat{\beta }}_j^{i,\lambda }\ne 0\}$. The final edge set E of 𝓖 can be defined via the intersection or the union operation of the local neighborhoods. An edge e_{i, j} between node v_i and v_j exists if $j\in N_i^{\lambda }\cap i\in N_j^{\lambda }$ or $j\in N_i^{\lambda }\cup i\in N_j^{\lambda }$. For edges in the set $j\in N_i^{\lambda }\cap i\in N_j^{\lambda }$, the edge weights, e_{i, j} and e_{j, i}, are estimated using the average of the two corresponding β entries. From a theoretical point of view, both edge selection choices are asymptotically consistent under certain technical assumptions [20]. The choice of the λ parameter controls the sparsity of the local neighborhood, which requires tuning [42]. We present our parameter selection strategy at the end of this section.

The second inference approach, (inverse) covariance selection, relies on the following penalized maximum likelihood approach. In the standard Gaussian setting, the related convex optimization problem reads: (4) where PD denotes the set of symmetric positive definite matrices {A: x^T Ax > 0, ∀x ∈ ℝ^p}, $\widehat{\Sigma }$ the empirical covariance estimate, ‖⋅‖₁ the element-wise L1 norm, and λ ≥ 0 a scalar tuning parameter. For λ = 0, the expression is identical to the maximum likelihood estimate of a normal distribution 𝓝(x∣0, Σ). For non-zero λ, the objective function (also referred as the graphical Lasso [22]) encourages sparsity of the underlying precision matrix Θ. The non-zero, off-diagonal entries in Θ define the adjacency matrix of the interaction graph 𝓖 which, similar to MB, depends on the proper choice of the penalty parameter λ. Originally, this estimator was shown to have theoretical guarantees on consistency and recovery only under normality assumptions [43]. However, recent theoretical [29, 44] work shows that distributional assumptions can be considerably relaxed, and the estimator is applicable to a larger class of problems, including inference on discrete (count) data. In addition, nonparametric approaches, such as sparse additive models, can be used to “gaussianize” the data prior to network inference [45]. We thus propose the following estimator for inferring microbial ecological associations. Given clr-transformed OTU data Z ∈ ℝ^n×p, we propose the modified optimization problem: (5) where $\widehat{\Gamma }$ is the empirical covariance estimate of Z, and Ω⁻¹ is the inverse covariance (or precision matrix) of the underlying (unknown) basis. As stated above, $\widehat{\Gamma }$ will be a good approximation for the basis covariance matrix $\widehat{\Omega }$ because p > > 0. The resulting solution is constrained to the set of PD matrices, ensuring that the penalized estimator has full rank p. The non-zero off-diagonal entries of the estimated matrix Ω⁻¹ define the inferred network 𝓖, and their values are the signed edge weights of the graph. To reduce the variance of the estimate, the covariance matrix $\widehat{\Gamma }$ can also be replaced by the empirical correlation matrix $\widehat{R}=D\widehat{\Gamma }D$, where D is a diagonal matrix that contains the inverse of the estimated element-wise standard deviations.

The covariance selection approach has two advantages over the neighborhood selection framework. First, we obtain unique weights associated with each edge in the network. No averaging or subsequent edge selection is necessary. Second, the covariance selection framework provides invertible precision and covariance matrix estimates that can be used in further downstream microbiome analysis tasks, such as regression and discriminant analysis [10].

Model selection. For both neighborhood and covariance selection, the tuning parameter λ ∈ [0, λ_max] controls the sparsity of the final model. Rather than inferring a single graphical model, both methods produce a λ-dependent solution path with the complete and the empty graph as extreme networks. A number of model selection criteria, such as Bayesian Information Criteria [46] and resampling schemes [47], have been used. Here we use a popular model selection scheme known as the Stability Approach to Regularization Selection (StARS) [48]. This method repeatedly takes random subsamples (80% in the standard setting) of the data and estimates the entire graph solution path based on this subsample. For each subsample, the λ-dependent incidence frequencies of individual edges are retained, and a measure of overall edge stability is calculated. StARS selects the λ value at which subsampled non-empty graphs are the least variable (most stable) in terms of edge incidences. For the selected graph, the observed edge frequencies indicate the reproducibility, and likely the predictive power, and are used to rank edges according to confidence.

Theoretical and computational aspects. Learning microbial graphical models with neighborhood or inverse covariance selection schemes has important theoretical and practical advantages over current methods. A wealth of theoretical results are available that characterize conditions for asymptotic and finite sample guarantees for the estimated networks [20, 29, 39, 43, 46]. Under certain model assumptions, the number of samples n necessary to infer the true topology of the graph in the neighborhood selection framework is known to scale as n = O(d³ log(p)), where d is the maximum vertex (or node) degree of the underlying graph (i.e. the maximum size of any local neighborhood). Additional assumptions on the sample covariance matrices reduce the scaling to n = O(d² log(p)) [39]. This implies that graph recovery and precision matrix estimation is indeed possible even in the p > > n regime, and that the underlying graph topology strongly impacts edge recovery. The latter observation means that, even if the number of interactions e is constant, graphs with large hub nodes, perhaps representing keystone species in microbial networks, or, more generally, scale-free graphs with, a few highly connected nodes, will be more difficult to recover than networks with evenly distributed neighborhoods. In addition to these theoretical results, a second advantage is that well-established, efficient, and scalable implementations are available to infer microbial ecological networks from OTU data in practice. Thus, SPIEC-EASI methods will efficiently scale as microbiome dataset dimensions grow (e.g., due to technological advances that increase the number of OTUs detected per sample). The SPIEC-EASI inference engine relies on the R package huge [49], which includes algorithms to solve neighborhood and covariance selection problems [20, 22], as well as the StARS model selection.

Generation of synthetic microbial abundance datasets

Estimating the absolute and comparative performance of network inference schemes from biological data remains a fundamental challenge in biology. In the context of gene regulatory network inference, recent community-wide efforts, such as the DREAM (Dialogue for Reverse Engineering Assessments and Methods) Challenges (http://www.the-dream-project.org/), have considerably advanced our understanding about feasibility, accuracy, and applicability of a large number of developed methods. In the DREAM challenges, both real data from “gold standard” regulatory networks (e.g., networks where the true topology is known from independent experimental evidence) and realistic in-silico data (using, e.g., the GeneNetWeaver pipeline [50]) are included. In the context of microbiome data and microbial ecological networks, neither a gold standard nor a realistic synthetic data generator exist. SPIEC-EASI is accompanied by a set of computational tools that allow the generation of realistic synthetic OTU data. As outlined in Fig 2, real taxa count data serve as input to SPIEC-EASI’s synthetic data generation pipeline. The pipeline enables one to: (i) fit the marginal distributions of the count data to a parametric statistical model and (ii) specify the underlying graphical model architecture (e.g., scale-free).

The NorTA approach. The parametric statistical model and network topology are then combined in the ‘Normal To Anything’ (NorTA) [51] approach to generate synthetic OTU data that resemble real measurements from microbial communities but with user-specified network topologies. NorTA [51] is an approximate technique to generate arbitrary continuous and discrete multivariate distributions, given (1) a target correlation structure R with entries ρ_{i, j} and (2) a target univariate marginal distribution U_i. To achieve this task, NorTA relies on normal-copula functions [51–53]. A n × p matrix of data U is sampled from a normal distribution with zero mean and a p × p correlation matrix R_N. For each marginal U_i, the Normal cumulative distribution function (CDF) is transformed to the target distribution via its inverse CDF. For any target distribution P with CDF Ξ, we can thus generate multivariate correlated data via (6) where U_N ∼ 𝓝(0, R_N) and $\Phi (\mathrm{\text{U}})={\int }_{-\infty }^U\frac{1}{\sqrt{2{\sigma }^2}}{e}^{\frac{-u^2}{2}}du$, the CDF of a univariate normal. In Fig 3a, we illustrate this process for bivariate Poisson and negative binomial data (n = 1000 and correlation ρ_ij = 0.7).

In the original NorTA approach, an element-wise monotone transformation c_U(⋅) with R_N = c_U(R) is applied to account for slight differences in correlation structure between normal and target distribution samples [51]. Here we neglect this transformation step because we observe that the log-transformed data from exponential count distributions, such as the Poisson and Negative binomial, are already close to R, provided that the mean is greater than one, particularly when the counts data are log-transformed (S1 Fig). In practice, SPIEC-EASI relies on routines from base R and VGAM packages [54, 55] to estimate the uniform quantiles of the normal data and to fit the desired CDF with estimated parameters.

Fitting marginal distribution to real OTU data. Prior to fitting marginal distributions to real data, several commonly used pre-proccesing steps are applied. For any given OTU abundance table of size n × K, we first select p non-zero columns. To account for experimental differences in sample sequencing, we then normalize samples to a median sequencing depth by multiplying all counts by the ratio of minimum desirable sampling depth to the total sum of counts in that sample and rounding to the nearest whole number, which is preferable to rarefaction [56]. These filtered and sequencing-depth-normalized data serve as the marginal counts, which are fit to a parametric distribution U_i and used as input to the NorTA approach. The concrete target marginal distribution depends on the actual microbiome dataset. For gut microbiome data (e.g. from HMP or APG), the zero-inflated Negative Binomial (ziNB) distribution is a good choice, as it accounts for both overdispersion [56, 57] and the preponderance of zero-count data points in microbial count datasets. The fitting procedure is done within a maximum likelihood framework. The corresponding optimization problem is solved with the Quasi-Newton methods with box-constraints, as implemented in the optim function in R [54]. In S1 Text, we use quantile-quantile plots to compare ziNB to several other candidate distributions (e.g., lognormal, Poisson, NB) and show that ziNB has superior fit.

Generation of network topologies and correlation matrices. Under normality assumption, the non-zero pattern of the precision matrix corresponds to the adjacency matrix of the underlying undirected graph. We use this property to generate target covariance (correlation) matrices originating from different graph topologies. The pipeline to generate a network structure for simulated data proceeds in three steps: (i) Generate an undirected graph, in the form of an adjacency matrix, with a desired topology and sparsity, (ii) convert the adjacency matrix to a positive-definite precision matrix by assigning positive and negative edge weights and appropriate diagonal entries, and (iii) invert Θ and convert the resulting covariance matrix Σ to a correlation matrix (R = DΣD, where, D is a diagonal matrix with diagonal entries $1/\sqrt{{\sigma }_i}$).

Among many potential graph structures, we focus on three representative network structures: band-like, cluster, and scale-free graphs (see Fig 3b for graphical examples). Maximum network degree strongly impacts network recovery, and thus our choice of network topologies spans a range of maximum degrees (band < cluster < scale-free). In addition, cluster and scale-free lend themselves to hypothetical ecological scenarios. Cluster graphs may be seen as archetypal models for microbial communities that populate different disjoint niches (clusters) and have only few associations across niches. Scale-free graphs, ubiquitous in many other facets of network biology (such as gene regulatory, protein-protein and social networks), serve as a baseline model for a microbial community that comprises (1) a few “keystone” species (hub nodes with many partners) that are essential for coordinating/stabilizing the community and (2) many dependent species that are sparsely connected to each other. The sparsity of the networks is controlled by the number of edges, e < p(p − 1)/2, in the graph. The topologies are generated according to the following algorithms, starting with an empty p × p adjacency matrix:

1. Band: A band-type network consists of a chain of nodes that connect only their nearest neighbors. Let e = e_used + e_available, the number of edges already used and available, respectively. Fill the next available off-diagonal vector with edges if and only if e_available ≥ number of elements in the off-diagonal.
2. Cluster: A cluster network comprises h independent groups of randomly connected nodes. For given p and e we divide the set of nodes into h components of (approximately) identical size and set the number of edges in each component to e_comp = e/h. For each component, we generate a random (Erdös-Renyi) graph for which we randomly assign an edge between two nodes in the cluster with probability $p=\frac{e_{\mathrm{\text{comp}}}}{h(h-1)/2}\cdot$
3. Scale-free: The distribution of degrees, the number of edges per node, in a scale-free network is described by a power law, implying that the central node or nodes (potentially keystone species in an ecological network) have proportionally more connections. We use the standard preferential attachment scheme [58] until p − 1 edges are exhausted.

After generation of these standard adjacency matrices, we randomly remove or add edges until the adjacency matrix has exactly e edges. All schemes generate symmetric adjacency matrices that describes a graph, with entries of 1 if an edge exists and 0 otherwise.

From the adjacency matrices, we generate precision matrices by uniformly sampling non-zero entries Θ_{i, j} ∈ [−Θ_max, −Θ_min]∪[Θ_min, Θ_max], where Θ_min, Θ_max > 0 are model parameters and describe the strength of the conditional dependence among the nodes. To ensure that the precision matrix is positive definite with tunable condition number κ = cond(Θ), we scale the diagonal entries Θ_{i, i} by a constant c using binary search. The precision matrix Θ is then converted to a correlation matrix R to be used as input to the NorTA approach.

Network inference on synthetic microbiome data

Given that no large-scale experimentally validated microbial ecological network exists, we use SPIEC-EASI’s data generator capabilities to synthesize data whose OTU count distributions faithfully resemble microbiome count data. By varying parameters known to influence network recovery (network topology, association strength, sample number) and quantifying performance on resulting networks, we rigorously assess SPIEC-EASI inference relative to state-of-the-art inference methods, SparCC [17] and CCREPE [9], as well as standard Pearson correlation.

Inference of the American Gut network

Thus far we have used the n₁ = 304 first-round AGP samples as a means to construct realistic synthetic microbiome data sets with SPIEC-EASI’s data generation module. In this section, we apply SPIEC-EASI inference methods to construct ecological association networks from the AGP data directly. To do this, we first filter out rare OTUs by selecting only the top 205 OTUs (to match the dimensionality of the synthetic data) in the combined AGP dataset (by frequency of presence) before adding a pseudo-count and total-sum normalization. Although there is no independent means to assess the accuracy of these hypothetical networks, we can assess their reproducibility and consistency. For each method, we first infer a single representative network of taxon-taxon interactions from Round 1 AGP abundance data. For SPIEC-EASI, the StARS model selection approach is used to select the final model network. For SparCC, we use a threshold ρ_t = 0.35 to construct a relevance network from the SparCC-inferred correlation matrix; i.e. an edge between nodes v_i, v_j is present in the SparCC network if ∣ρ_{i, j}∣ > ρ_t [17]. Similarly, we use a q-value cut-off of 10⁻²⁴ to create an interaction network from CCREPE-corrected significance scores of Pearson’s correlation coefficient [9]. For each method, we thus arrive at a reference network that can be considered the hypothetical gold standard. We then use the n₂ = 254 Round 2 AGP samples as an independent test set and learn a new model network from these data alone. We measure consistency between the two network models by computing the Hamming distance between the reference and new network models, i.e., the difference between the upper triangular part of the two adjacency matrices. For the present data, the Hamming distance can vary between p(p − 1)/2 = 20910 (no edges in common) and a minimum of 0 for identical networks. Confidence intervals for Hamming distances can be obtained by combining Round 1 and 2 samples into a unified dataset, repeatedly subsampling these data into two disjoint groups of size n₁ and n₂, and repeating the entire inference procedure.

Fig 6a shows network reproducibility for SPIEC-EASI methods, SparCC, and CCREPE. The S-E(MB) has smallest the Hamming distance, followed by S-E(glasso), SparCC, and CCREPE. In S-E(MB), the edge disagreement is roughly 50 with very small error bars. At the other extreme, CCREPE edge disagreement is 250 edges and highly variable.

These numerical experiments clearly demonstrate that SPIEC-EASI networks are more reproducible than other current methods.

Finally, we use each inference method to construct a candidate American Gut microbiome association network from the unified dataset of size n₁ + n₂ = 558 (Fig 6c). We analyze the differences between the reconstructed networks by quantifying the number of unique and shared predicted edges (Fig 6b). All four methods agree on a core network that consists of 127 edges. These edges are mostly found within OTUs of the same taxonomic group. This phenomenon, termed assortativity, has also been observed in other microbial network studies [11]. Assortativity is one of the most salient features of the AGP-derived networks, and, for all networks, the assortativity coefficients for each network are close to unity (e.g., maximum assortativity, S11 Fig). The SparCC network comprises about twice as many edges as the SPIEC-EASI networks. SparCC infers 147 distinct edges; these additional edges correspond to negative associations between OTUs of Ruminococcaceae (genus Faecalibacterium) and Enterobacteriacae families (various genera) and a dense web of correlations within Enterobacteriacae OTUs. Similarly, CCREPE identified 152 edges uniquely, with many negative edges between Enterobacteriaceae and Lachnospiraceae (genera: Blautia, Roseburia and unknown); additionally, CCREPE uniquely predicted positive edges between the Lachnospiraceae and Ruminococcaceae (genus: Faecalibacterium). Both SPIEC-EASI methods produce relatively sparse networks by comparison. S-E(glasso) infers a total of 271 total edges (with one unique edge), and S-E(MB) infers 206 edges with 25 unique edges. In scale with edge predictions, both CCREPE and SparCC infer networks with large maximum degree (33 and 30, respectively), while the S-E(MB) and S-E(glasso) networks have a maximum degree of sixteen and eight, respectively (S11 Fig). However, even though CCREPE and SparCC predict a similar number of total edges, the global network properties are distinct. CCREPE predicts a higher maximum betweenness centrality and a larger number of nodes in the largest connected component (100).

In summary, analysis of the AGP networks suggests that the SPIEC-EASI inference schemes construct more reproducible taxon-taxon interactions than SparCC and CCREPE and infer considerably sparser model networks than the other two methods. These observations may be explained as follows: SparCC and CCREPE aim to recover correlation networks, which contain both direct edges as well as indirect (e.g., spurious) edges (due to correlation alone). SparCC and CCREPE may recover indirect edges less robustly than direct edges, an explanation that would be consistent with the Hamming distance reproducibility analysis. In addition, all methods’ resulting networks suggest that the topology of the American Gut association network cannot be attributed to a specific network class. Instead, these networks are a mixture of band, scale-free, and cluster network type, and they exhibit high phylogenetic assortativity within highly connected components.