Investigation of gut microbiome association with inflammatory bowel disease and depression: a machine learning approach

Database generation

The datasets used for the analyses were retrieved from the Inflammatory Bowel Disease Multi-Omics Database (IBDMDB) (Schirmer et al., 2018), which is part of the Integrative Human Microbiome Project (NIH HMP Working Group et al., 2009). The IBDMDB database contains a wide array of omics data (e.g., 16S and shotgun metagenomic, metatranscriptomic, proteomic and host genomes) of 132 individuals classified by IBD diagnostic in ulcerative colitis, Crohn’s disease and controls. Participants provided bi-weekly stool samples at five hospitals in the United States. Metagenomic and metatranscriptomic data was processed as described in Schirmer et al., 2018 (Abubucker et al., 2012; Truong et al., 2015).

Subject selection

From this dataset, the 70 unique participants who answered an additional self-reported depression and anxiety questionnaire during registration (the answers to which are listed in the HMP2 metadata, column EC to EL) were selected. As the questionnaire model was not specified, only individuals with raw scores over 6 on this test was considered as showing “signs of depression”. To calculate the raw scores, a severity scale was generated, with the following scores: 0, never; 1, rarely; 2, sometimes; 3, often; 4, always. The scores were then summed to give a final total. In the case of individuals undergoing multiple tests, the lower score was used. We selected a low threshold in order to be able to identify putative dysbiotic individuals that were not experiencing severe depression symptoms. All the others were classified as “no sign of depression”. The combination between the test and the IBD diagnosis divided the dataset in six groups: Crohn’s disease with no detectable sign of depression (CD; n=15), Crohn’s disease with signs of depression (CDD; n=20), ulcerative colitis with no sign of depression (UC; n=4), ulcerative colitis with signs of depression (UCD, n=11), signs of depression but no inflammation (nonIBDD; n=7) and the control group: no inflammation/no depression (nonIBD; n=13). As the experimental design of the IBDMDB consisted of a longitudinal study, each subject contributed several times to this study, and all the samples used for this analysis were sequenced by shotgun sequencing as described in Schirmer et al. The resulting datasets for metagenomic and metatranscriptomic consist of 1084 and 566 samples, respectively. The final tables after pre-processing consist of 1486 columns, including Participant ID, data type, diagnostic, sex, mental score, and nested columns on the relative values of the different taxa.

Data analysis

For each of the six groups, abundance matrices of the metagenomic data, metatranscriptomic data, and the combination of metagenomics and metatranscriptomics were used for random forest classification. Each of the datasets was divided randomly into a training set (90% of the individuals) and a validation set (10% of the individuals). Random forest analysis were performed using the library Scikit-learn 0.19.1 (Pedregosa et al., 2011) on the training sets to identify the most important species involved in discriminating the samples without losing predicting power. A 1000-fold cross-validation for the combined dataset, and 500-fold for metagenomic and metatranscriptomic data (see Figure 1), considering one model for each iteration was performed and only the most important species in the construction of this model was retained. Only models with a precision classification >80% were considered, and among the considered models, only species that appeared more in more than one were selected. Afterwards, the validation sets were run with the selected species only to measure the possible loss of predictive capability and computed the area under the receiver operating characteristic (auROC) curve for the prediction of the validation set classes as a performance metric.

Figure 1. Workflow of the k-Fold Cross Validation approach for Random Forest.

First the data gets split into train and validation sets (A). The train dataset will be iterated by the Cross Validation algorithm (B), while the validation set will be spared to test the model trained only with the reduced feature list (C).

Statistical analysis

In order to assess the significance of the differences between the abundances of the selected species, we performed a one-way ANOVA (Scipy 1.0.0, Jones et al., 2001) with a Tukey’s honest significant difference (HSD) post-hoc test. This test makes pair-wise comparisons between the different means to see which classes are different. For clarity, confidence intervals for Tukey’s HSD test can be found in Supplementary Materials (Supplementary Figure 1 and Supplementary Figure 2).

The functional activity of the selected species was retrieved from the HUMAnN metatranscriptomic analyses described above. Only the pathways in which the selected species are involved and those that were different between the groups from the ANOVA test were selected and the correlation between these species was calculated using Spearman’s correlation coefficient. A significance level of 0.05 was applied for all statistical tests.

Species selection and model validation

The random forest cross-validation selection of the most informative species showed a combined list of 24 species, as can be seen in Figure 2. The validation models for DNA, RNA and the combined dataset shows micro-averaged auROC values of 0.96, 0.91 and 0.99, respectively (Supplementary Figure 3–Supplementary Figure 5). This small loss of information suggest a relevant role of the selected species in the interaction of both conditions, while the capability of the model to classify the validation data with with great accuracy shows that our model can generalize its results and it’s not overfitting.

Figure 2. Venn diagram for the species selected for each dataset.

All species exhibited differences in at least one group in a one-way ANOVA (alpha=0.05, Supplementary Table 1), and no significant differences were found between DNA and RNA abundances for these species (Supplementary Table 2). This list of putative species pretends to be a trade-off between the all-relevant and minimal informative approaches. We chose this approach ir order to get as broad of a list as possible while avoiding artifacts related to the longitudinal nature of the dataset.

In order to assess the effect of the small sample size of group UC, the same procedure was made grouping all samples with IBD together. As expected, we see some difference in the species selected. However, the species that showed stronger differences in the previous classification were also the stronger ones, with most of the species overlapping. The interesting exception is Faecalibacterium prausnitzii that was absent.

The non-dysbiotic microbiome

The analyses showed an increase in the number of species from the genus Bacteroides in dysbiotic groups compared with the control (nonIBD) (Figure 3), as has been reported in other dysbiotic samples (Bloom et al., 2011), with the exception of Bacteroides dorei, which is more abundant in nonIBD than in any other group. Aside from Bacteroides dorei, nonIBD samples had a higher abundance of Alistipes shahii and Ruminococcus bromii, while a typical species associated with nonIBD, Faecalibacterium prausnitzii, was significantly decreased in nonIBDD and CD.

Figure 3.

DNA (A) and RNA (B) taxonomic abundances for the selected species. Abundances were quantified by the relative abundances of their sequences, and for each level they should sum to 1 (including unclassified sequences).

Crohn’s disease abundance changes in depression

Both of the Crohn’s disease-related groups (CD and CDD) showed higher abundances of Bacteroides ovatus and Bacteroides uniformis. However, CD samples exhibited higher abundances for several specific species, including Bacteroides xylanisolvens, Parasutterella excrementihominis and Bacteroides fragilis, compared with CDD, but decreased abundance of Faecalibacterium prausnitzii, which did not differ significantly in abundance between nonIBD and CDD groups.

Ulcerative colitis changes in depression

Ulcerative colitis samples had the most distinctive microbiome profile. Several species, including Burkholderiales bacterium 1_1_47, Bacteroides eggerthii and Bacteroides finegoldii were characteristic of this group, and absent in the others, except for B. finegoldii, which was also present in a lower abundance in nonIBD samples. Only UCD samples exhibited an increased abundance of Bacteroides fragilis, Bacteroides vulgatus and Haemophilus pittmaniae, this last species being almost exclusive to the UCD group.

Non-IBD changes in depression

The nonIBDD was the group with the highest number of changes in microbiome diversity when compared with its non-depressed counterpart (Table 1). However, most of those changes followed a similar pattern in other dysbiotic groups.

A notable change was observed in Faecalibacterium prausnitzii, which was present in almost the same abundances in nonIBD, UCD and CDD samples, and a high variability in UC while being significantly lower in CD and nonIBDD (Supplementary Table 3 and Supplementary Table 4). This is particularly interesting, since this species is considered to have anti-inflammatory activity. It seems counterintuitive to find a depleted population of one of the species most associated in the literature with a healthy microbiome compared to an IBD one in a group that doesn’t show any inflammatory process. However, Parabacteroides goldsteinii was increased in nonIBDD and was depleted in all IBD groups in comparison with control samples. The Parabacteroides genre have been associated previously with anti-inflammatory activity (Neff et al., 2016; Schirmer et al., 2016), so the increase in abundance of this bacteria may explain why the nonIBDD microbiome is not associated with inflammation in the gut.

Other than Parabacteroides goldsteinii, nonIBDD samples did not contain other characteristic groups, and, more notably, none of the selected species was specific for depressed or non-depressed phenotypes.

Microbial functional activity

Regarding the functional activity of these species, seven pathways that were more abundant in dysbiotic groups than in nonIBD were identified (Supplementary Figure 1) and were correlated between each other and inversely correlated with most of the others (Supplementary Figure 2 and Supplementary Table 5). Those pathways are folate transformations II, N10-formyl-tetrahydrofolate biosynthesis, de novo L-ornithine biosynthesis, superpathway of pyridoxal 5’phosphate biosynthesis and salvage, phosphopantothenate biosynthesis I, preQ0 biosynthesis and queuosine biosynthesis. Folate (vitamin B9) and pyroxidal 5’-phosphate (vitamin B6) deficiencies have been linked both to depression (Coppen & Bolander-Gouaille, 2005; Hvas et al., 2004; Mitchell et al., 2014), as they are key for the synthesis of several neurotransmitters, and IBD (Pan et al., 2017; Yakut et al., 2010), although this association is not well understood and does not seem to be evidence of causation. Increased levels of L-ornithine derivatives have also been linked to depression (Zheng et al., 2010). However, even if nonIBDD have the highest activity for almost all of these pathways, CD and UC were also significantly increased, while functional activity in CDD was generally lower and non-significant in some pathways. Moreover, UCD did not differ from nonIBD in any of them.

This difference in functional activity again highlights the lack of a concrete pattern of gut microbiome abundance between depressed groups.