- The matrix with read counts is filtered (to discard taxa not occurring in enough samples) and normalized sample-wise for each row group (body sites in our case). After normalization, abundances of higher-level taxa are computed as the sum of member taxa abundances.
- For each of four measures (Pearson, Spearman, Kullback-Leibler, Bray Curtis), all allowed pair-wise scores are computed (parent-child relationships between taxa are forbidden) and the top- and bottom-ranked 1,000 edges are kept for each measure, resulting in an initial multigraph with 1,000 x 2 x 4 = 8,000 edges.
- The network is recomputed for 1,000 permutations (where taxon profiles are shuffled). The permutations are carried out edge-wise and for Pearson and Spearman are followed by a renormalization of the full matrix. Renormalization can shift the permutation (null) distribution of correlation coefficients away from zero, thus mitigating compositional effects.
- In the next step, the network is recomputed for 1,000 bootstrapped matrices, that is matrices obtained from the original matrix by sampling columns (with equal probability) with replacement. Bootstrapping computes a confidence interval around the edge score.
- Edges with scores not within the limits of the 95% confidence interval defined by the bootstrap distribution are discarded.
- A measure- and edge-specific p-value is then obtained from the Gauss curve defined by the mean and standard deviation of the bootstrap distribution. The p-value is the area under that curve from the mean of the permutation distribution to the left or right tail (the tail depends on whether the measure is a distance or a similarity). In R, the p-value can be computed thus: p-val = pnorm(mean_permut, mean=mean_boot, sd=sd_boot)
- In the article, we pooled variances of the bootstrap and permutation distribution as follows:
pooled_sd = sqrt((var_boot + var_permut)/2)

and computed the p-value as: p-val = pnorm(mean_permut, mean=mean_boot, sd=pooled_sd)

Variance pooling considers the standard deviation of the null distribution, which is otherwise not taken into account. - Since the p-value is computed in a one-sided test, very high p-values correspond to negative relationships (low similarities and high distances). These are converted into low p-values by computing 1 - p-value for all p-values above 0.5.
- Each edge is now supported by a set of measure-specific p-values, which are dependent, since the measures are correlated (that is the ranks they assign to edges are correlated). To merge the p-values, we used Sime's method, which keeps the minimum p-value as the merged p-value of the edge.
- Merged p-values were then multiple-test corrected with Benjamini-Hochberg (BH). In our article, multiple regression was included as a network inference method, which introduces dependencies between merged edge p-values, thus BH was performed with the Yekutieli correction factor. Without the regressions, BH alone is sufficient. Since p-values of edges not part of the initial edge set are not included in the p-value list, p-values resulting from BH are currently over-optimistic. The initial edge list can of course be extended (by setting less restrictive initial thresholds), but at a computational cost that makes network building intractable.
- Edges with p-values above a confidence interval of 0.05 were discarded.
- Additional filter steps, such as keeping only edges supported by at least two methods, are recommended.