NetCoMi: network construction and comparison for microbiome data in R

1The rapid development of high-throughput sequencing techniques offers new possibilities for investigating the microbiome across different habitats and provides the opportunity to discover relationships between the composition of microbial communities and their environment. Sequencing technologies include shotgun metagenomic sequencing, where whole genomes of microbes in the sample are sequenced, and targeted amplicon sequencing, where microbial abundances are derived from sequencing a specific marker gene (e.g. 16S rRNA, 18S rRNA or ITS). Following McLaren et al., we refer to marker-gene and metagenomic sequencing techniques as MGS techniques.

2MGS techniques allow the quantification of taxon abundances in the sample in form of operational taxonomic units (OTUs), amplicon sequencing variants (ASVs), or metagenomic OTUs (mOTUs). However, due to limitations in the sequencing process, MGS measurements often provide only relative or compositional information with each component expressing the relative frequencies of a taxon in the sample. Furthermore, the observed MGS measurements represent only a sample of the true microbial composition present in the biological material. Accordingly, it is likely that not all taxa occurring in a sample are measured due to technical limitations in library preparation and the sequencing process. The observed number of reads is thus only a noisy measurement reflecting the probability of the corresponding organisms to be present. Moreover, MGS data collections comprise a high amount of zero counts, contain samples with varying sequencing depths (sum of counts per sample), and are usually high-dimensional, i.e. the number of taxa is much larger than the sample size.

3A common exploratory analysis approach for microbiome survey data is the estimation of microbe–microbe association networks, allowing for high-level insights into the global structure of microbial communities. Although a wide range of general-purpose tools is available for visualizing and analysing networks, including the R packages igraph, statnet and network, the special nature of current MGS data requires careful statistical treatment. For instance, statistical network analysis methods ignoring the compositional structure of the data may lead to spurious results (also known as compositional effects). Existing compositionally aware approaches for measuring and estimating microbial associations include correlation estimators, such as SparCC, partial correlation estimators and proportionality.

4Although the analysis of a single microbial association network can provide insights into the general organizational structure of a microbial community, researchers are often more interested in how microbial associations change across different conditions. For instance, it is often desired to find relationships between microbial compositions, their inherent connectivity and an underlying phenotype, e.g. the health status of patients. This task thus requires the quantitative comparison of networks across conditions.

5Current approaches for comparing networks between two conditions can be divided into two types: (i) differential association analysis focusing on differences in the strength of single associations, and (ii) differential network analysis, analysing differences between network metrics and network structure between two conditions. Differential associations can be used as the basis for constructing differential networks, comprising only differentially associated node sets. Existing tools for network comparison either require pre-computed networks or are tailored for specific biological purposes, including protein–protein interactions or gene functions.

6In this paper, we introduce NetCoMi (Network Construction and comparison for Microbiome data), a comprehensive R package, that specifically focuses on microbiome data analysis. NetCoMi integrates previously disjoint microbial network inference and analysis tasks into a single coherent computational workflow and allows the user to construct, analyse and compare microbial association networks in a fast and reproducible manner. The complete NetCoMi framework is shown in Figure 1. NetCoMi provides a wide range of existing methods for data normalization, zero handling, edge filtering, and a selection of association measures, which can be combined in a modular fashion to generate the microbial networks. For network analysis, several local and global network properties are provided, which can be visualized in network plots to enable a descriptive comparison. Quantitative comparison between two different networks is available via the integration of appropriate statistical tests. Furthermore, our package enables (i) the generation of differential microbial networks and (ii) the construction of sample similarity networks (using, e.g. the Bray-Curtis measure), which can serve as a high-level visual summary of the heterogeneity of the microbiome sample collection.

7The process of constructing microbial association networks starts with a matrix containing non-negative read counts originating from a sequencing process. Although NetCoMi can handle any kind of count abundance data as input, our main focus here is on compositional read count data coming from MGS experiments. The total read counts of a sample with taxa are a composition summing up to a constant, the sequencing depth. The sequencing depth differs from sample to sample and is predefined by technical factors leading to sparse data with many zeros. Thus, preprocessing steps (step 1 in Figure 1) are recommended, or even mandatory depending on the association measure (Supplementary Table S1).

8To simplify the graphical interpretation of an association network and the computational processes, it is reasonable to filter out a certain set of taxa as first data preparation step (step 1a in Figure 1). See Table S2 for the options available in NetCoMi.

9The excess number of zeros in the data is a major challenge for analysing microbiome data because parametric as well as non-parametric models may become invalid for data with a large amount of zeros. Moreover, many compositionally aware measures are based on so-called log-ratios. Log-ratios have been proposed by Aitchison as the basis for statistical analyses of compositional data as they are independent of the total sum of counts. More precisely, for two variables and the log-ratio of relative abundances is equal to the log-ratio of the absolute abundances. However, log-ratios cannot be computed if the count matrix contains any zeros, making zero handling necessary (step 1b in Figure 1). Several zero replacement strategies have been proposed. Table 1 gives an overview of the different types of zeros that have been suggested as well as existing approaches for their treatment.

10Normalization techniques are required to make read counts comparable across different samples (step 1c in Figure 1). The normalization approaches included in NetCoMi are summarized in Table 2. A description of these methods is available in. Note that forcing the read counts of each sample to a unique sum (as done with total sum scaling) does not change the compositional structure. Rather, proportions are always compositional, even if the original data are not. Aitchison suggested using the centered log-ratio (clr) transformation to move compositional data from the simplex to real space. Badri et al. have shown that variance-stabilizing transformations (VSTs) as well as the clr-transformation produce very similar Pearson correlation estimates, which are more consistent across different sample sizes than TSS, CSS and COM methods.

11Association estimation is the next step in our workflow (step 2a in Figure 1) to obtain statistical relations between the taxa. Common association measures include correlation, proportionality and conditional dependence. For all three types of association, compositionally aware approaches have been proposed, which are summarized in Table 3. Further information on these measures is available in Supplementary material 1. To ensure wide applicability of NetCoMi, the package comprises also traditional association measures, which are not suitable for application on read count data in their original form.

12Compositionality implies that if the absolute abundance of a single taxon in the sample increases, the perceived relative abundance of all other taxa decreases. Applying traditional correlation measures such as Pearson’s correlation coefficient to compositions can lead to spurious negative correlations, which do not reflect underlying biological relationships. It has been shown that the resulting bias is stronger, the lower the diversity of the data.

13A possible approach to reduce compositional effects is to normalize or transform the data in a compositionally aware manner and apply standard correlation coefficients to the transformed data. A common method is the clr transformation (Table 2), which is a variation of the aforementioned log-ratio approach. Traditional correlation measures provided in NetCoMi are the Pearson correlation, Spearman’s rank correlation and the biweight midcorrelation (bicor). Bicor as part of the R package WGCNA is more robust to outliers than Pearson’s correlation because it is based on the median instead of the mean of observations.

14Popular compositionally aware correlation estimators include SparCC (Sparse Correlations for Compositional data), CCLasso (Correlation inference for Compositional data through Lasso) and CCREPE (Compositionality Corrected by REnormalization and PErmutation), also called ReBoot method. SparCC is one of the first approaches developed for inferring correlations for compositional data and has become a widely used method. However, SparCC has some limitations, namely that the estimated correlation matrix is not necessarily positive definite and may have values outside [-1,1]. Furthermore, the basic algorithm is repeated iteratively to reinforce the sparsity assumption and to account for uncertainties due to random sampling, leading to high computational complexity. CCLasso is also based on log-ratios but aims to avoid the disadvantages of SparCC by using a latent variable model to infer a positive definite correlation matrix directly. The CCREPE approach operates directly on relative read counts, which are permuted and renormalized in order to detect correlations induced by compositionality alone.

15Lovell et al. argue that correlations cannot be inferred from the relative abundances in a compositionally aware manner without any assumptions and propose proportionality as an alternative association measure for compositional data. If two relative abundances are proportional, then their corresponding absolute abundances are proportional as well:. Thus, proportionality is identical for the observed (relative) read counts and the true unobserved counts. Lovell et al. suggest proportionality measures based on the log-ratio variance, which is zero when and are perfectly proportional. This variance, however, lacks a scale that would make the strength of association comparable. For this reason, the proposed proportionality measures and are modifications of the log-ratio variance based on clr transformed data that come with a scale. Due to its analogy to correlations, we included in the NetCoMi package (Table 3 for the formula). is a symmetric measure with values in [-1,1], where 1 corresponds to perfect proportionality.

16Conditional dependence expresses the relation between two variables conditioned on all other variables in the data set. Hence, this is a measure of direct relations between each two taxa, whereas (marginal) correlations cannot differentiate between direct and indirect dependencies. Three estimators of conditional dependence are included in NetCoMi: SPRING (Semi-Parametric Rank-based approach for INference in Graphical model), SPIEC-EASI (Sparse Inverse Covariance Estimation for Ecological Association Inference) and gCoda.

17All three approaches use graphical models to infer the conditional independence structure from the data. For multivariate Gaussian data, the graph structure can be inferred from the non-zero elements of the inverse covariance matrix, where each entry is related to scaled negative partial correlation. Loh and Wainwright relaxed the Gaussianity assumption and established relationships between the inverse covariance matrix and the edges of a graph for discrete data. All three approaches for estimating a conditional dependence graph presume that the assumptions of the data generation process are fulfilled so that the graph structure can be reliably inferred from the count matrix. We also consider the assumptions on the data generation process as satisfied and use conditional dependence and partial correlation equivalently in the following.

18gCoda uses a Majorization-Minimization algorithm to infer from the data based on maximizing a penalized likelihood. In SPIEC-EASI, two approaches are provided to obtain from the observed count data: neighborhood selection and sparse inverse covariance selection (also known as “graphical lasso”). SPRING uses a semi-parametric rank-based correlation estimator which can account for excess zeros in the data and applies neighborhood selection to infer conditional dependencies. All three measures assume a sparsely connected underlying network.

19Using one of the aforementioned association measures, an association matrix with entries expressing the relation between pairs of taxa and is computed. The next step is sparsification and transformation into distances and similarities (step 2b in Figure 1) resulting in an adjacency matrix with entries as numerical representation of the microbial network, where nodes (or vertices) represent the taxa. The different options available in NetCoMi are illustrated in Figure 2 (Supplementary Figure S1 for a more detailed version of this chart).

20Since the estimated associations are generally different from zero, using them directly as adjacency matrix results in a dense network, where all nodes are connected to each other and consequently only weighted network measures are meaningful. Instead, the association matrix is usually sparsified to select edges of interest.

21One possible sparsification strategy consists of defining a cut-off value (or threshold) so that only taxa with an absolute association value above this threshold are connected. This filtering method is available for all types of association. The conditional independence measures SPRING, SPIEC-EASI and gCoda already include a model selection approach, making the filter step unnecessary.

22For correlations, statistical tests are available as alternative sparsification method, allowing only significant associations to be included in the network. Student’s t-test and a bootstrap approach are implemented for identifying correlations significantly different from zero. The resulting P-values need to be adjusted to account for multiple testing. NetCoMi includes all adjustment methods available in the R function p.adjust (stats) and, in addition, two methods for multiple testing under dependence having higher power than the common Benjamini-Yekutieli method: (i) fdrtool for controlling the local false discovery rate, and (ii) the adaptive Benjamini-Hochberg method where the proportion of true null hypotheses is estimated using convex decreasing density estimation as proposed by Langaas et al..

23P-values arising from bootstrapping might be problematic because they can be exactly zero if none of the bootstrap-correlations is more extreme than the observed one. Since P-values of exactly zero cannot be corrected for multiple testing, NetCoMi’s bootstrap p-values are corrected by adding pseudo counts: where is the number of generated samples with a test statistic at least as extreme as the observed one, and is the number of repetitions.

24As illustrated in Figure S1, the sparsified associations are then transformed into dissimilarities/distances, which are needed for computing network measures based on shortest paths (Section 2.4). Following van Dongen and Enright, we included different distance metrics depending on how negative associations should be handled (Figure S1). Available are the options: (i) “unsigned”:, leading to a low distance between strongly associated taxa (positively as well as negatively), (ii) “signed”:, where the distance is highest for strongly negative associated taxa, and (iii) “signedPos”: “signed” distance with setting negative associations to zero. A dissimilarity measure based on the topological overlap matrix (TOM) is also available in NetCoMi.

25The adjacency matrix contains similarities of the form. These similarity values are used for network plot and network metrics based on connection strength (Section 2.4). NetCoMi also offers two options for constructing an unweighted network via generating a binary adjacency matrix from the sparse association matrix where the user can decide whether or not negative associations should be included in the unweighted network (Figure S1).

26Furthermore, we included adjacency matrix constructions via the soft-thresholding approach from the WGCNA package, which is only available for correlations (see blue path in Figure 2). Rather than using hard thresholding (e.g. via cutoff value or statistical tests), Zhang and Horvath suggested raising the estimated correlations to the power of a predefined value greater than 1 to determine the adjacencies (Supplementary material 2.2 for more details). Small correlation values are thus pushed toward zero becoming less important in the network. The resulting similarities are used as edge weights for network plotting and network metrics based on connection strength. Note that this approach leads to a fully connected network. Following Zhang and Horvath, only TOM dissimilarities are available in NetCoMi for computing shortest paths and clusters when soft-thresholding is used.

27We analyse a constructed network by calculating network summary metrics (step 3 in Figure 1), which are amenable to group comparisons. Alternatively, NetCoMi offers the possibility to analyse and visualize single networks.

28Several network statistics require shortest path calculations. A path between two vertices and is a sequence of edges connecting these vertices such that and are each at one end of the sequence and no vertices are repeated. The length of a path is the sum of edge weights, where the weight of an edge is a real non-negative number associated with this edge. The shortest path between two nodes is the path with minimum length. In NetCoMi, edge weights are defined in two ways: (i) for properties based on shortest paths, dissimilarity values are used implying that the path length between two taxa is shorter the higher the association between the taxa is. (ii) For properties based on connection strength, the corresponding similarities are used as edge weights (Section 2.3 and Figure S1 for details on distance and similarity calculation).

29Network centrality measures offer insights into the role of individual taxa within the microbial community. We consider degree, betweenness, closeness and eigenvector centrality (Table 4). Using these measures, so-called hubs (or keystone taxa) can be determined. This analysis step is of high interest to researchers since hub nodes may correspond to taxa with a particularly important role in the microbial community. What is deemed “important”, depends on the scientific context. A selection of definitions for hub nodes is given in Table 4.

30Clustering methods are appropriate to identify functional groups within a microbial community. A cluster (or module) is a group of nodes that are highly connected to one another but have a small number of connections to nodes outside their module. The user of NetCoMi can choose between one of the clustering methods provided by the igraph package and hierarchical clustering (R package hclust from stats package). Both similarities and dissimilarities can be used for clustering (Table 4).

31Global network properties are defined for the whole network and offer an insight into the overall network structure. Typical measures are the average path length, edge and vertex connectivity, modularity, and the clustering coefficient (Table 4).

32Using dissimilarity between samples rather than association measures among taxa leads to networks where nodes represent subjects or samples rather than taxa. Analogous to microbiome ordination plots, these sample networks express how similar microbial compositions between subjects are and thus provide insights into the global heterogeneity of the microbiome sample collection.

33The following dissimilarity measures are available in NetCoMi: Euclidean distance, Bray-Curtis dissimilarity, Kullback-Leibler divergence (KLD), Jeffrey’s divergence, Jensen-Shannon divergence, compositional KLD, and Aitchison distance. Details are available in Table S3. Only Aitchison’s distance and the compositional KLD are suitable for application on MGS data whereas the others may induce compositional effects when applied to raw count data without appropriate transformation (Table 2).

34The workflow for constructing, analysing and comparing sample similarity networks is described in Supplementary material 2.5. For these networks, the same network properties as for microbial networks are available in NetCoMi (Table 4). While the estimated dissimilarity values between two samples and are used for network properties based on shortest paths, the corresponding similarities, calculated by are used for properties based on connection strength and edge weights in the network plot. Accordingly, highly connected nodes are subjects with a microbial composition similar to many other subjects. Furthermore, as in usual cluster analysis, clusters represent subjects with similar bacterial composition but with the advantage, that the solution is visualized in the network plot.

35NetCoMi’s network comparison module (step 5 in Figure 1) focuses on investigating the following questions in a quantitative fashion: (i) Is the overall network structure different between two groups? (ii) Are hub taxa different between the two microbial communities? (iii) Do the microorganisms build different “functional” groups? (iv) Are single pairs of taxa differentially associated among the groups? To meet these objectives, NetCoMi offers several network property comparison modes as well as the estimated associations itself between two groups. These approaches are commonly used in other fields of application and have been adapted to the microbiome context. Each method includes statistical tests for significance.

36To perform network comparison, the count matrix is split into two groups according to the user-defined group indicator vector. All steps of network construction and analysis described in (Section 2) are performed for both subsets separately. The estimated associations, as well as network characteristics, are then compared using the methods described below.

37We next detail NetCoMi’s differential network analysis capabilities. All approaches are applicable for microbial association networks and sample similarity networks, respectively.

38NetCoMi uses permutation tests to assess for each centrality measure and taxon whether the calculated centrality value is significantly different between the two groups. The null hypothesis of these tests is defined as, where and denote the centrality value of taxon in group 1 and 2, respectively. A standard non-parametric permutation procedure is used to generate a sampling distribution of the differences under the null hypothesis (Supplementary material 3.1). The same approach is used to test for significant group differences in the global network characteristics listed in Table 4.

39A set of most central nodes can be defined in two ways: In the first approach, under the assumption that the centrality values are log-normal distributed, the set of most central nodes contains nodes with a centrality value greater than a certain quantile of the fitted log-normal distribution. In the second approach, the empirical quantile can be used directly without any distributional assumption. Both approaches are included in NetCoMi, where the quantile can be freely chosen in each case.

40The Jaccard index (Supplementary material 3.2) can then be used for assessing how different the two sets of most central nodes (regarding a certain centrality measure) are between the groups. This index ranges from zero to one, where a value of one corresponds to two equal sets and zero means that the sets have no members in common. Following Real and Vargas, we included an approach to test whether the observed value of Jaccard’s index is significantly different from that expected at random (Supplementary material 3.2). Note that this approach cannot make a statement on whether the two sets of most central nodes are significantly different.

41NetCoMi offers several network partitioning and clustering algorithms (Table 4). One way to assess the agreement of two partitions is via the Rand index (Supplementary material 3.3). Like Jaccard’s index, the Rand index ranges from 0 to 1, where 1 indicates that the clusters are exactly equal in both groups. The original Rand index is dependent on the number of clusters making it difficult to interpret. Instead, NetCoMi uses the adjusted Rand index. The adjusted values lie in the range, where 1 corresponds to identical clusterings and 0 to the expected value for two random clusterings. Consequently, positive index values imply that two clusterings are more similar and negative values less similar than expected at random. Following, NetCoMi uses a permutation procedure to test whether a calculated value is significantly different from zero. However, this test does not signify whether the clusterings are significantly different between the groups. Details about this approach and its implementation in R are given in Supplementary material 3.3.

42NetCoMi also allows to test differences between the estimated associations themselves rather than network properties. This analysis is referred to as differential association analysis. Table 5 shows the three approaches available in NetCoMi (Supplementary material 4 for details).

43Fisher’s z-test is a common method for comparing two correlation coefficients, assuming normally distributed z values (the transformed estimated correlations) and thus a correct specification of their variance. Novel approaches without normality assumption have been proposed, but are restricted to Pearson correlations. Therefore, we implemented a resampling-based procedure as non-parametric alternative, which is applicable to association measures other than correlation. The Discordant method as the third available method is also based on Fisher’s z values, but groups correlations with a similar magnitude and direction based on mixture models.

44These approaches enable the construction of a differential network, where only differentially associated taxa are connected. More precisely, two taxa are connected if their association is either significantly different between the two groups (to a user-defined significance level) or identified as being different by the Discordant method.

45We use data from the GABRIEL Advanced Surveys (GABRIELA) to illustrate the application of NetCoMi version 1.0.1. GABRIELA is a multi-center study, carried out in rural areas of southern Germany, Switzerland, Austria and Poland, which provides new insights into the causes of the protective effect of exposure to farming environments for the development of asthma, hay fever and atopy. The study comprises the collection of biomaterial and environmental samples including mattress dust from children’s rooms and nasal swabs, for which 16S rRNA amplicon sequencing data are available. After some preprocessing steps (Table S4), a total of bacterial genera remain for mattress dust for a subset of subjects. The nasal data set consists of genera for subjects.

46The NetCoMi functions corresponding to the main steps of our proposed workflow (Figure 1) are netConstruct for network construction (Section 2.3), netAnalyze for network characterization (Section 2.4), and netCompare for network comparison (Section 3). An overview of the exported (and thus usable) NetCoMi functions together with their main arguments is given in Figure 3. In the following, we use the SPRING method as a measure of partial correlation for constructing exemplary microbial association networks and Aitchison’s distance as a dissimilarity measure for constructing sample similarity networks. The construction of a differential network is described in Supplementary material 5.5. We also demonstrate the integration of external networks into NetCoMi’s workflow using BAnOCC-derived microbial association networks (Supplementary material 6). In Supplementary material 7, we investigate the variability of microbial networks among different network construction methods. The results reveal strong differences between the networks based on different association measures as well as different normalization methods whereas the zero replacement methods available in NetCoMi lead to quite similar networks, if the remaining arguments are fixed.

47In principle, the netConstruct function allows the specification of any combination of methods for association estimation, zero treatment and normalization. However, since NetCoMi is mainly designed to handle compositional data, a warning is returned if a chosen combination is not compositionally aware.

48To generate the network shown in Figure 4A, we pass the combined count matrix containing dust samples from Ulm and Munich to netConstruct. Filter parameters are set in such a way that only the most frequent taxa are included in the analyses leading to a read count matrix for our data. Depending on the association measure, a method for zero handling, normalization and sparsification could be chosen (Figure S1 and Table S1). Since these steps are already included in SPRING, they are skipped in our example. If associations have already been estimated in advance, the external association matrix can be passed to netConstruct instead of a count matrix. In this mode, the user only needs to define a sparsification method and a dissimilarity function.

49Besides the available options for transforming the associations into dissimilarities (Supplementary Figure S1), the function also accepts a user-defined dissimilarity function. The choice of dissimilarity function also influences the handling of negative associations. In Figure 4, we use the “signed” distance metric, where strong negative correlations lead to a high dissimilarity (and thus to a low edge weight in the adjacency matrix). Figure S4 shows a network plot where the “unsigned” metric is used.

50netConstruct returns an object of class microNet, which can be directly passed to netAnalyze to compute network properties (Table S5). Applying the plot function to the output of netAnalyze leads to a network visualization, as shown in Figure 4A. In this plot, network characteristics are emphasized in different ways: hubs are highlighted, clusters are marked by different node colors, and node sizes are scaled according to eigenvector centrality. Alternatively, node colors and shapes can represent features such as taxonomic rank. Node sizes can be scaled according to other centrality measures or absolute/relative abundances of the corresponding taxa. Node positions are defined using the Fruchterman-Reingold algorithm. This algorithm provides a force-directed layout aimed at high readability of the network by placing pairs of nodes with a high absolute edge weight close together and those with low edge weight further apart.

51The plot method implemented in NetCoMi includes several options for selecting nodes or edges of interest to facilitate the readability of the network plot without influencing the calculated network measures. In Figure 4B, for instance, we display the 50 bacteria with highest degree, facilitating network interpretability.

52NetCoMi identified four clusters and the following five hub nodes: Aerococcus, Atopostipes, Brachybacterium, Rikenellaceae RC9 gut group and [Eubacterium] coprostanoligenes group. The strongest positive association is between Ezakiella and Peptoniphilus with a partial correlation of 0.55 (in the green cluster in Figure 4A). The strongest negative correlation is -0.09 between Alistipes (in the blue cluster) and Lactobacillus (in the red cluster).

53For network comparison, the combined count matrix is again passed to netConstruct, but this time with an additional binary vector assigning the samples to one of the two centers. This leads to the network plots shown in Figure 5 (Table S7 for the summary of network properties). The layout computed for the Munich network is used for both networks to facilitate the graphical comparison and making differences clearly visible. We observe only slight differences in the estimated associations. Furthermore, both networks show similar clustering and agree on three (out of five) hub nodes: Atopostipes, Brachybacterium and [Eubacterium] coprostanoligenes group.

54The quantitative comparison is done by passing the R object returned from netAnalyze to netCompare. Comparisons of all global measures included in NetCoMi and the five genera with the highest absolute group difference for degree and eigenvector centrality, respectively, are given in Table 6. Supplementary Table S8 extends the output to the 10 genera with the highest absolute group difference, and also includes betweenness and closeness centrality. For none of the four centrality measures, any significant differences are observed, confirming the descriptive analyses. Pseudomonas, Pedobacter and Rikenellaceae RC9 gut group as hub taxa in only one of the networks show high differences in eigenvector centrality. However, the differences are not deemed significant (Table 6). Global network properties (upper part of Table 6) are also not significantly different for.

55Table 7 summarizes Jaccard indices expressing the similarity of sets of most central nodes and the hub nodes among the two centers. They do not imply any group differences, which would be indicated by a small probability P. Similarly, the adjusted Rand index indicates a high similarity of the two clusterings, which is also highlighted in Figure 5.

56We next consider sample similarity networks from mattress dust and nasal swabs of the same subjects. The process of constructing and comparing sample similarity networks is analogous to the one for association networks. Figure 6 shows the two networks using Aitchison’s distance.

57NetCoMi’s quantitative network analysis (Tables S10–S13) reveals strong differences between these networks. The sets of most central nodes are significantly different for all four centrality measures (shown by Jaccard’s index). Hub nodes are also completely different (Jaccard index of zero). We also observe several significantly different global network properties (Table S12). Furthermore, the node’s degree, betweenness and closeness centrality values differ significantly between the groups for several subjects (Table S13), implying that a single subject plays a different role dependent on the investigated microbial habitat.

58Clustering analysis broadly identifies three sample groups for mattress dust and nasal swab samples, respectively. In Figure 6, we highlight the partial overlap between the two clustering solutions via color coding. Clusters that have at least 100 nodes in common are plotted by matching colors. This reveals that the red cluster, seen in both networks, shares similar samples across the two habitats.

59NetCoMi’s plotting functionality also allows to draw nodes in different shapes, corresponding to additional categorical covariate information about the samples. For instance, in Figure 6, the two node shapes correspond to the two different study centers (Ulm and Munich). This feature could potentially highlights confounding of groups of samples and available covariates. For our example here, however, we observe no noticeable pattern in the clusters with respect to study center.

60With NetCoMi we offer an easy-to-use and versatile, integrative R package for the construction, analysis and comparison of microbial networks derived from MGS data. Our package provides a wide variety of compositionally aware association measures, including SparCC, proportionality, SPIEC-EASI and SPRING. The latter method also enables the analysis of recent quantitative microbiome data sets when both amplicon and quantitative cell count or spike-in control data are available. NetCoMi also incorporates standard association measures, thus widening the scope of the package beyond applications to compositional data, and it connects to the popular WGCNA package, enabling principled soft-thresholding of correlations and dissimilarity transformations based on topological overlap. The package includes a dedicated list of methods for handling excess zeros in the count matrix and for data normalization in order to account for the special characteristics of the underlying MGS data prior to association estimation.

61A unique feature of the NetCoMi framework is its ability to perform differential network and differential association analysis in a statistically principled fashion. Differential network analysis allows not only to uncover the global role of a taxon in the overall network structure but also its changing influence under varying conditions. Differential association analysis, on the other hand, can directly assess which associations significantly change across conditions, providing concrete hypotheses for follow-up biological perturbation experiments.

62Similar to phyloseq’s plot_net function, NetCoMi also enables network representation and comparison of the MGS data samples themselves, using popular sample dissimilarity or distance measures, such as the Bray-Curtis dissimilarity and the Aitchison distance. Network analysis of the resulting sample-to-sample or subject-to-subject networks can give insights into the heterogeneity of the collected data. For instance, identified hub subjects are subjects with “representative” bacterial compositions that may comprise archetypical microbial patterns in the studied population. Sample similarity network analysis thus extends standard sample ordination or cluster analysis, making fine-grain structures of the available microbial sample collection visible.

63Whereas methods and tools for the individual analysis steps, such as biological network estimation (e.g. for recent contributions) and (differential) biological network analysis are available, NetCoMi offers a unique and modular R software framework that integrates the complete process of estimating, analysing and comparing microbial networks. From an end user’s perspective with a specific microbiome dataset and scientific question in mind, NetCoMi will thus facilitate both faster development and reproducibility of the microbial network analysis workflow.

64Even though the modular design of NetCoMi allows the user to perform a wide variety of computational workflows, going from the primary data all the way to potentially significant network features, every step of the analysis still warrants careful scientific consideration. For instance, the choice of an association or dissimilarity measure will likely affect all further steps of network analysis and comparison, as we have shown for microbial association networks. However, there is no general consensus in the community about the “right” way to estimate and analyse microbial networks. This is reflected in the heterogeneity of recent simulation studies to assess and compare the performance of compositionally aware as well as traditional association measures. To put some of these studies into context, we give a selection of simulation studies examining the association methods used in NetCoMi in Supplementary Table S17.

65Common shortcomings in current simulation studies include the lack of a universal standard (i) to generate realistic synthetic microbial data with a prescribed ground truth, (ii) to perform comparable model selection and (iii) to report generalizable performance metrics. Comparative studies often concentrate on the performance in edge recovery, for instance, via precision-recall curves or distances between true and estimated associations. However, networks derived from penalized marginal correlations (such as SparCC) and partial correlations (such as SPRING) are statistically difficult to compare, thus requiring special modifications which, in turn, limits cross-study comparisons (e.g. Yoon et al. where SparCC correlations are transformed into SparCC partial correlations). Moreover, the sole focus on edge recovery may obfuscate other aspects of correct model recovery, including the shape of degree distributions or detecting hub nodes. For instance, methods with a similar edge recovery precision may greatly vary regarding their ability to determine hub nodes. Thus, if the correct detection of hub nodes is of major interest to the user, present comparative microbial network studies will give little guidance.

66Finally, we posit that simulation studies accompanying articles that introduce new methods might also be inherently biased. Neutral comparison studies, i.e. studies that are independent of any new method development are rare in our context or do often not include recently published methods. The main impediment regarding neutral comparison studies is, however, the fact that, to date, no large-scale ground truth network of real biological microbial interactions is available. Such biological gold-standard networks, as available in other contexts (e.g. gene-gene or transcription-factor-gene interactions), would greatly facilitate future comparative studies.

67In the absence of a “best method” for microbial network inference and analysis, NetCoMi is intended to give researchers the possibility to apply a consistent and reproducible analysis workflow on their data. Ideally, the selection of the workflow building blocks should be set up once and independent of any hypothesis about the data, thus avoiding the fallacy of starting “fishing” for results that best suit a previously formulated hypothesis. NetCoMi can, however, serve as an ideal tool for principled sensitivity analysis of the inferred results, for instance, by assessing how different normalization and zero handling methods affect the estimated networks, their structural properties and their comparison. Finally, we envision NetCoMi to provide a useful framework for future simulation studies that evaluate and compare the performance of different association measures and network inference tools in a reproducible fashion.

68The current version of NetCoMi is designed to model networks from a single domain of life, e.g. bacteria, fungi or viruses. However, microbes from different domains of life often share the same habitat and likely influence each other. Joint cross-domain network inference already revealed considerable alterations of the overall network structure and network features, compared to their single-domain counterparts. Extending NetCoMi to cross-domain network analysis is thus an important future development goal. Likewise, environmental factors, such as chemical gradients and temperature, as well as batch effects are known to influence microbial abundances and composition and thus bias network estimation. In the current NetCoMi version, we assume that the user has already corrected the microbiome data for these latent influences. However, several inference methods can directly incorporate known or unknown latent factors into network learning. Including or connecting these approaches with NetCoMi will likely increase the robustness and generalizability of future workflows.

69A core feature of NetCoMi is the use of statistical tests at various stages of the computational workflow. For instance, statistical tests can be employed for edge selection in network sparsification. Since statistical power depends on sample size, the sparsified structure of a network will likely depend on the number of available samples. A comprehensive understanding between number of samples, sparsification and network structure is currently elusive. NetCoMi relies on permutation tests for several statistical tests. The lower limit of P-values arising from permutation tests is directly related to the number of available permutations (Section 4). This already required proper adjustment of the calculation for small numbers. For extended simulation studies, NetCoMi’s dependence on permutation tests may prove to be computer intensive. Thus, integration of less demanding alternatives to permutation tests would represent a welcome feature in NetCoMi.

70Finally, despite incorporating a comprehensive list of methods in our R package, we do not claim completeness. This can hardly be achieved in the vibrant field of microbiome research where new methods are constantly developed. We alleviated this shortcoming via NetCoMi’s modular structure which allows certain parts of our workflow to be combined with external methods. For instance, users can input a user-defined association or dissimilarity matrix rather than a data matrix, and then proceed with NetCoMi’s standardized network analysis modules.

71In summary, we believe that NetCoMi is a useful addition to the modern microbiome data analysis toolbox, enabling rapid and reproducible microbial network estimation and comparison and ideally leading to robust hypotheses about the role of microbes in health and disease.

72Current high-throughput sequencing count data carry only relative or compositional information, thus requiring dedicated statistical analysis methods.

73NetCoMi is a comprehensive R package that implements the complete workflow of constructing, analysing, and comparing microbial association networks.

74NetCoMi integrates an extensive list of methods that take into account the special characteristics of marker-gene and metagenomic sequencing data, including methods for zero count handling, normalization and association estimation.

75The package also offers functionality for constructing sample similarity networks as well as differential networks including appropriate methods for identifying differentially associated taxa.

76The data underlying this article will be shared on reasonable request to the corresponding author upon approval by the involved institutional review boards.