## 4 Materials and Methods
### Bioinformatic and Statistical Analysis  
Raw sequences were processed using a pipeline that combined PANDAseq [64] and QIIME 2 [65]. Specifically, paired-end sequences were first assembled into a single-end amplicon using PANDAseq, retaining only assembled reads in the range of 350 to 500 nucleotides. An error correction step was then performed using the USEARCH11 'fastq_filter' module [66], which set a maximum error rate of 0.03 to discard low-quality sequences. Sequences were then binned into amplicon sequence variants (ASVs) using the DADA2 pipeline [67], with chimeras being removed at the same time. Taxonomic assignment was performed using the VSEARCH algorithm [68] to align ASVs with the Greengenes database v13.8. Alpha diversity was computed using the 'diversity core-metrics-phylogenetic' plugin from QIIME 2, as well as multiple metrics, including the number of observed ASVs, the Shannon index, and Faith's Phylogenetic Diversity. Beta diversity was estimated via the same plugin using both quantitative and qualitative metrics, such as weighted/unweighted UniFrac, Jaccard, and Euclidean measures, which were used for principal coordinates analysis (pCoA). All statistical analyses were performed using R/rStudio 4.2.2 software. The pwr R package ([[https://CRAN.R-project.org/package=pwr](https://CRAN.R-project.org/package=pwr)]([https://CRAN.R-project.org/package=pwr](https://CRAN.R-project.org/package=pwr)), accessed on 2 June 2022) was implemented to perform a post hoc power analysis. For this procedure, the alpha value was set at 0.05, and the effect size was estimated based on the response rate. The resulting power was 0.88, which could be considered sufficiently robust given the exploratory nature of this study. pCoA plots were generated using the vegan R package v2.6-4 ([[http://www.cran.r-project.org/package=vegan/](http://www.cran.r-project.org/package=vegan/)]([http://www.cran.r-project.org/package=vegan/](http://www.cran.r-project.org/package=vegan/)), accessed on 2 June 2022), and data separation was tested via permutational analysis of variance with pseudo-F ratios (Adonis tests). The variance in beta diversity of a given group, i.e., the distance in beta diversity between samples within the same group, was calculated as the average squared distance to the mean. The Wilcoxon test with continuity correction was used to compare alpha diversity, compositional structure, and metabolite levels between groups. To find multivariate associations among GM profiles, metabolomic data, and patient metadata, we employed multi-block sparse partial least square discriminant analysis (sPLS-DA), as implemented in the DIABLO/mixOmics packages in R [69]. The model developed for this approach was tuned according to the developer's guidelines ([[http://mixomics.org](http://mixomics.org)]([http://mixomics.org](http://mixomics.org)), accessed on 30 June 2022), with repeated cross-validation performed with 10 folds and 100 repeats. The approach allowed us to pick the best number of components for sPLS-DA computation, reducing the overall bit error rate. A similar tuning approach was then implemented, based on the developer's guidance, to select the appropriate features retained from the different layers of information for the DIABLO model. Results were first visualized via plotting the individual values for each data layer (plotIndiv function of the mixOmics package), which represented the contribution of each layer to patient stratification. Next, clustered image maps were generated using cimDiablo to highlight the relationships between the multi-omics panel and the outcome variable. Finally, the network function of mixOmics was used to construct a relevance network that consisted of a set of features linked via a DIABLO association value, which we considered to be a network edge list and plotted as a network plot using Cytoscape [70].