Network building - suggestion of a pipeline and discussion

Ensemble pipeline

1. 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.
2. 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.
3. 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.
4. 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.
5. Edges with scores not within the limits of the 95% confidence interval defined by the bootstrap distribution are discarded.
6. 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)
7. 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.
8. 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.
9. 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.
10. 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.
11. Edges with p-values above a confidence interval of 0.05 were discarded.
12. Additional filter steps, such as keeping only edges supported by at least two methods, are recommended.

Further considerations