Sample FASTQ files were acquired from the European Nucleotide Archive (ENA) browser. The sample metadata including covariates such as BMI, life stage, sex, IBD subtype, disease status and sample type was acquired from the corresponding publication’s supplementary materials or the QIITA microbiome platform (Gonzalez et al., 2018). Samples collected from individuals in North America and with more than 3,000 counts following processing of the raw reads (as described in methods) were retained from each dataset. The dataset ENA accessions and technical information regarding the samples in each dataset are available in Table 1.

The following fifteen studies were included in our dataset:

1. The American Gut cohort is from a large, open platform which collected samples from individuals in the US to identify associations between microbiomes, the environment, and individual’s phenotype (McDonald et al., 2018). We included available samples that did not contain any self-reported diseases in the metadata.

Taxonomy abundance tables were generated from the FASTQ files using QIIME2 (v2020.2) (Bolyen et al., 2019). Only samples from the same study were processed together. Reads were trimmed to remove low quality base pairs (trimming parameters listed in Table 1), chimeras removed, and sequences denoised using Dada2 (Callahan et al., 2016) or Deblur (for GLS and AG due to technical issues of processing these samples with Dada2) (Amir et al., 2017). Closed reference OTU clustering with the Silva 132 99% reference database (Quast et al., 2013; Yilmaz et al., 2014; Glöckner et al., 2017) was performed with the cluster_features_closed_reference function from QIIME2 plugin VSEARCH (v2.7.0) (Rognes et al., 2016) at 99% similarity. The resulting centroid sequences were classified with a Naive Bayes classifier (Wang et al., 2007) at a 99% confidence cut-off. In order to train the Naive Bayes classifier, the sequences of the 16S rRNA hypervariable region sequenced in the respective study (either V3-V4 or V4) were extracted from the Silva 132 99% full length 16S OTU reference with the extract-reads function from the QIIME2 feature-classifier plugin. The extracted reads and the corresponding taxonomy labels were used to train the Naive Bayes classifier with the QIIME2 plugin feature-classifier’s fit-classifier-naive-bayes function. Resulting taxonomic feature tables were collapsed to species (level 7) and genus (level 6) classification for further analysis.

Functional abundance tables were generated using PICRUSt2 (v2.3.0) from the OTU abundance table and representative OTU sequences generated using QIIME2. We generated abundance tables from the six different databases incorporated into PICRUSt2: Clusters of Orthologous Groups of proteins (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KO), Enzyme Commission (EC), Pfam protein domain (PFAM), TIGR protein family (TIGRFAM) and MetaCyc pathways. Each database is independently curated and provides information on different aspects of the functional properties present in the microbiome.

Following taxonomy classification and inference of functional abundance, features present in less than 10% of the samples within each study in the training set were pruned from the dataset. Following pruning of the training set, the test set was subset to those features as well.

When possible, normalization and transformation methods were implemented using python (v3.6.12) and R (v3.6.3) packages with the methods already incorporated. For CLR and ILR transformation, zero values were first replaced with a multiplicative replacement function that replaces zeros with a small positive value (equal to 1/N², where N equals the number of columns) and ensures the sum of the row remains 1 (Martín-Fernández, 2003) prior to transformation with the clr and ilr functions, respectively, from the python package SciKit-Bio (v0.5.2). CLR performs a log transformation of abundance values, which are normalized by the geometric mean of all features. ILR uses a change of coordinate space projection calculation to transform proportional data (or relative abundances) to a new space with an orthonormal basis, in this case the J.J.Egozcue orthonormal basis (Egozcue, 2003).

For TSS normalization, the counts for each feature were divided by the sum of all feature counts in the sample with a custom python function. The method constrains the sample row sum to one, aiming to similarly scale all samples while maintaining biological information of microbial abundances. For ARS normalization, the TSS normalized values were transformed with the sqrt function followed by the arcsin function from the python package Numpy (v1.19.2). The LOG normalization was also applied to the TSS normalized values using the log function from numpy following replacement of all zero values with the multiplicative replacement function.

For VST normalization, we used the varianceStabilizingTransformation function in the R package DESeq2 (v1.26.0). VST aims to factor out the dependence of the variance in the mean abundance of a feature. The method numerically integrates the dispersion relation of the feature mean fitted with a spline, evaluating the transformation for each abundance in the feature. VST normalization was performed by normalizing the training set separately from the test set as the normalization is dependent on all samples present in the dataset.

We explored three methods for batch effect correction: naive zero-centering, an empirical Bayes method, and a negative binomial regression method. The naive zero-centering batch effect correction entails centering the mean of each feature within a batch to zero (Nygaard et al., 2016). We also assessed MMUPHin, an empirical Bayes method designed specifically for zero-inflated microbial abundance data. MMUPHin estimates parameters for the additive and multiplicative batch effects, using normal and inverse gamma distributions, respectively. The estimated parameters are then used to remove the batch effects from the dataset (Ma et al., 2020). MMUPHin was implemented with a custom python script. Lastly, ComBat-seq fits the feature counts to a negative binomial regression model to estimate the batch specific parameters. The batch specific parameters are used to calculate a ‘batch-free’ distribution which the raw counts are mapped to in order to obtain the final corrected data (Zhang et al., 2020). ComBat-seq was implemented with the ComBat_seq function from the R package SVA (v3.38.0) (Parker et al., 2014). We considered a batch as the whole dataset or split a dataset into multiple batches when the metadata indicated different sample preprocessing methods or samples were processed in different locations.

Two variations of MMUPHin and ComBat-seq were implemented to ensure the batch effects were corrected from the training and test sets separately. For variation #1, the test study’s samples were removed to generate the training set. Training set batch effect correction was completed, with the sample type (stool/biopsy) provided as a biological covariate, and the corrected values used for training the ML model. For the test set in variation #1, batch effect correction was performed with the full dataset and then the test study’s samples were collected to form the test set. The corrected test was then used to assess the model’s classification performance (Supplemental Figure S1, variation #1). For variation #2, batch effect correction was completed on the training set prior to model training with both the sample type and disease label (UC/CD/Control) provided as biological covariates. The model’s classification performance was then assessed on samples from the test dataset which were not corrected (Supplemental Figure S1, variation #2). Lastly, feature abundance for some samples following batch effect correction of the OTU dataset with MMUPHin were all zero. To ensure these samples were compatible with the compositional transformation methods, we set all features of these samples to equal 1/N (where n is the number of features) prior to transformation. For all other normalization methods, the feature abundance was not adjusted.

We assessed the ability of the three batch effect correction methods to improve mixing of samples from different batches with the beta-diversity metric Aitchison distance, equivalent to the Euclidean distance between CLR transformed microbiome data (Quinn et al., 2018), and the Local Inverse Simpson Index (LISI) (Polański et al., 2020). The genus abundance dataset was filtered to include features present in 10% of the samples from at least one batch, followed by CLR transformation and batch effect correction applied in the same manner as in our classification pipeline. Dimensionality reduction with principal component analysis (PCA) using the PCA function from the python package SciKit-Learn (v0.22.1) was performed with the resulting clr transformed and batch effect corrected values. The first two principal components were used to generate the scatter plots visualizing the separation of labels corresponding to diagnosis and sample batch. In order to visualize the different labels, a 95% confidence ellipse was added for each disease or batch label in the respective graph.

In order to quantify the mixing or separation of disease and batch labels, the LISI metric was calculated using the first 50 principal components with the compute_lisi function from the python package harmonypy (v0.0.5) (Korsunsky et al., 2019). LISI selects the nearest neighbors of a sample to calculate the inverse Simpson’s index for the diversity of labels surrounding the sample. For the batch integration LISI (iLISI) score, the batch label was provided and for disease LISI (dLISI) score the disease label was provided. The iLISI and the dLISI were calculated for every sample following each batch effect correction method with the first 50 principal components following dimensionality reduction with PCA. The median score for each method was determined and scaled with the overall minimum and maximum scores to a range between 0 and 1 (Tran et al., 2020). A Wilcoxon signed-rank test with Benjamini and Hochberg p-value correction was applied to determine if the values were significantly different. In order to assess the overall effect of the batch effect correction methods, the harmonic mean (also referred to as the F1 score) incorporating both the iLISI and dLISI was calculated as previously described (Lin et al., 2019; Tran et al., 2020).

The deep learning models were built with the python package Tensorflow (v2.2.0). The models were trained for up to 100 epochs with a batch size of 16 and samples shuffled. The best weights were selected using early stopping (EarlyStopping callback) by monitoring the validation loss (5% split of the training set) with a min_delta = 1 × 10⁻³ and patience = 10.

The generalizability of each model, normalization or transformation, and batch effect correction method, were determined through a cross validation strategy which assessed predictive performance on previously unseen batches of samples. As there were 15 studies, we iterated through the full dataset 15 times, generating the training set by removing all samples from a single study to create a separate test set. Feature selection was performed with the training set, followed by normalization and batch effect correction with the respective methods to both the training and test sets. Lastly, the number of non-IBD control and IBD samples in the training set were balanced by randomly subsampling the label with the greater number of samples, while maintaining the proportion of samples from each batch, disease label (UC/CD/Control), and sample type (stool/biopsy).

To measure the performance of the various normalization, transformation, batch effect correction, and model combinations, we calculated two metrics for binary classification: F1 score and Matthews correlation coefficient (MCC) (Wardhani et al., 2019). In the LODO cross-validation, some test sets contained only a single label thereby preventing the calculation of the MCC and F1 Score. In order to be able to include these predictions in the overall performance, the F1 Score and MCC were calculated by combining the test set predictions of all 15 folds. As the number of samples in a study ranged from 23 to 1,279, there was potential for the overall performance metrics to be skewed by the predictions of a single study with a large number of samples. Whereas in a typical 5-fold cross-validation each fold is weighted equally, with the classification performance determined separately for each fold and then an overall average calculated. Therefore, in order to calculate the pipeline’s classification performance with equal weighting to each study, the confusion matrix for each study was generated and normalized by the number of samples in the study. The average proportion of true positives, true negatives, false positives, and false negatives across the 15 studies was used to generate an overall confusion matrix. The overall F1 score and MCC were calculated from the averaged confusion matrix with the following equations.

For the pipeline combinations with the best overall classification performance, the classification accuracy was also determined for each individual dataset. The accuracy was calculated from the normalized confusion matrix with the following equation:

We assessed the difference in classification performance of patients within five different metadata variables (age, BMI, sample type, sex, IBD type), each with two categorical labels coded as 0 or 1. The samples were grouped by the five variables and the classification performance metric within each group calculated. For the logistic regression analysis, the performance metric was input as the dependent variable and the five metadata groups as the independent variables. The MCC score was scaled with the MinMaxScaler from SciKit-Learn to a range from 0 to 1 as required for the logistic function. The logistic function was fit with the Logit function from the statsmodels (v0.11.1) python package.

In order to determine the importance of each taxonomy, we collected the features’ gain value from our second-best pipeline composed of CLR normalized, zero-centered, genus abundance features with an XGBoost Classifier. The gain values were collected from the trained XGBoost classifier in each LODO iteration separately.

Differential taxonomy abundance was performed with Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) (v1.0.5) (Lin and Peddada, 2020). The fold change between control samples and IBD samples (UC and CD) was determined and a Bonferroni multiple comparison correction was applied to the p-values.

In order to assess the cross-batch performance of each pipeline, we implemented a LODO cross validation approach. We collected 16S rRNA gene next generation sequencing data from 15 studies in North America for a total of 7,707 samples, comprising 55% healthy and 45% IBD samples, of which 56% are CD and 44% are UC (Table 2). The mean sequencing depth for each study ranged from 17547.5±7,522.5 to 402901.1±77178.6 reads per sample (Table 1). We included studies that contained IBD and control samples, only IBD samples, or only control samples in order to better recapitulate a diagnostic scenario, where any distribution of IBD and non-IBD samples may be received and processed together for disease classification. As some of the data processing methods share information across samples, it was important to test the pipeline’s performance on datasets with a range of label distributions that may be encountered.

In order to evaluate the different pipelines, we completed 15 cross-validation iterations with the classification model trained on n-1 datasets (all samples from a single dataset were removed) and the model performance assessed on the removed dataset (Figure 1). We evaluated the ability to classify samples from patients with IBD or non-IBD controls using different combinations of three taxonomic feature sets or six functional feature sets, six normalization methods, two transformation methods, six batch effect correction methods, and nine machine learning models (Figure 1). The binary classification performance of each combination of feature set, normalization or transformation, batch effect correction, and machine learning model was assessed with two classification metrics: F1 score and Matthews Correlation Coefficient (MCC) (Wardhani et al., 2019). In order to include datasets with a single label in the overall performance assessment, we calculated the overall classification performance from a confusion matrix comprising the average true positive, true negative, false positive, and false negative proportions across the 15 studies. Performance metrics are reported as median (25% percentile-75% percentile) of all pipelines containing the respective component. We performed a Mann-Whitney U test to determine if the performance was significantly different between different pipeline components.

Machine learning classification models identify decision boundaries within the feature space to separate one datapoint from another. For some ML models (BNB, Linear SVC, LR), these boundaries are linearly constrained, whereas others (RF, KNN, MLP, Radial SVC, XGBoost) can identify more complex, non-linear relationships between features and class. We assessed the generalizability of three linear and five non-linear ML models across the taxonomy and functional feature sets.

The various ML models were sorted by the performance metrics’ median value across all pipelines containing the respective model. Across both performance metrics, the five models with the greatest classification performance were consistently the non-linear models, with the three linear models exhibiting the lowest classification performance (Figure 2A). The two models with the greatest median performance across were XGBoost, with a median F1 Score of 72.5 (64.1–76.4) and MCC of 55.6 (40.6–62.4) and Random Forest, with a median F1 Score of 71.2 (64.1–74.8) and MCC of 53.4 (40.2–59.6).

In addition, some combinations of linear models and functional features performed worse than randomly classifying the samples. For example, COG functional features with a logistic regression model had a median MCC of 8.8. Overall, non-linear models had significantly better classification performance than linear models, with a median F1 score [68.1 (61.4–73.5) vs 60.0 (47.1–69.2), p-value < 0.0001] and MCC [49.1 (35.5–57.7) vs 34.0 (9.6–50.1), p-value < 0.0001] (Figure 2B). Of the non-linear models, XGBoost and Random Forest had significantly higher F1 score and MCC than a MLP, KNN, or radial SVC (Figure 2C).

In order to further assess whether the non-linearity of a model improves classification in the context of microbiome data, we compared linear and non-linear variations of a support vector machine and logistic regression. Comparison of the two variations enables direct analysis of the impact of decision boundary constraints on performance, independent of differences in model architecture (Figure 2D). The non-linear (Radial) version of logistic regression had significantly greater F1 score (64.9 (60.3–71.4) vs 58.5 (39.2–67.9), p-value < 0.0001) and MCC [42.3 (33.8–53.9) vs 31.3 (1.6-47.7), p-value < 0.0001] than the linear (Linear) logistic regression. Additionally, the radial support vector machine had significantly greater F1 score [67.3 (60.4–72.2) vs 60.7 (41.1–66.7), p-value < 0.0001) and MCC (46.6 (33.1–55.2) vs 34.4 (0.2–46.3), p-value < 0.0001] than the linear version. In conclusion, non-linear models provided more accurate IBD classification, likely due to the complex relationships between features and disease labels.

Other ML model architectures, such as convolutional neural networks (CNNs), are commonly used for classification problems with defined structure in the input data, such as image classification. In the context of microbiome data, the CNN MDeep adds structure to OTU features through hierarchical agglomerative clustering of the phylogeny-induced correlation between OTUs (Wang et al., 2021). As MDeep is currently only developed for OTU features, we assessed whether this CNN architecture led to greater classification performance with OTU abundance than our MLP architecture. Comparison of each performance metric across all normalization, transformation, and batch effect correction methods indicated MDeep performance was not significantly different from our MLP model (MDeep F1 Score of 71.6 (67.5–74.0) vs MLP F1 Score of 70.3 (67.9–73.8), p-value > 0.05, and MDeep MCC Score of 54.0 (46.7–58.8) vs MLP F1 Score of 51.9 (47.4–57.9), p-value > 0.05) (Figure 2E).

Taxonomic features (species, genus or OTU) are predominantly used as input for ML models, whereas it is less common to use inferred functional features from PICRUSt2 as input. However, previous studies have identified lower inter-individual variation of the gut microbiome’s inferred functional profile than taxonomy (Davenport et al., 2014; Moustafa et al., 2018), suggesting that functional features may lead to better classification performance and generalizability. We processed the 16S sequencing samples with QIIME2 and PICRUSt2 to obtain taxonomy and functional feature abundance estimates, respectively.

For each ML model, we assessed the performance with taxonomy and functional abundance features in combination with a normalization or transformation and batch effect correction method. The batch effect correction methods were limited to Zero-Centering and no batch effect correction, as they were the only two performed on all datasets. Across all classification performance metrics, the taxonomic features classified IBD samples more effectively than functional features (Figure 3A). Comparison of performance with taxonomy and functional features confirmed the significantly higher F1 score [70.1 (63.8–74.7) vs 59.1 (51.2–64.3), p-value < 0.0001] and MCC [52.3 (41.2–59.5) vs 31.5 (18.1–41.7), p-value < 0.0001] for classification of IBD samples with taxonomic features (Figure 3B). Therefore, ML models using taxonomic features from this dataset lead to better classification of IBD samples than functional features.

We further investigated whether specific taxonomic ranks allow for an improved disease classification. Taxonomic classification consists of seven hierarchical ranks, with kingdom and species at the top and bottom, respectively. Each consecutively lower taxonomy rank provides greater resolution of the gut microbiome’s composition while also increasing data sparsity, which can negatively affect an ML model’s performance (Karlsson and Bostrom, 2014). Previous literature comparing different taxonomy ranks for disease classification indicated that lower ranks, down to genus, improved performance (Bang et al., 2019). We assessed whether the trend for improved classification continued with the species rank and OTUs, despite their increasing sparsity. While no significant performance difference was observed between species and genus ranks, both displayed significantly higher classification performance than OTU features (Figure 3C).

Due to the significantly better performance of non-linear classification models and taxonomic features, our subsequent analysis of normalization, transformations, and batch effect correction methods utilized only taxonomic feature sets and non-linear models.

We assessed normalization and transformation methods which account for different biases commonly observed in next-generation sequencing data: compositionally, heteroskedasticity, and skewness. We selected two transformations designed for compositional data: the isometric log ratio (ILR) and centered log ratio (CLR) (Pawlowsky-Glahn and Egozcue, 2006). We selected two normalization methods which aim to reduce the heteroskedasticity: the arcsine square root (ARS) transformation (Campbell et al., 1970) of the total sum scaling (TSS) values and the variance stabilized transformation (VST) from the R package DESeq2 (Love, 2014). Next, we assessed a log transformation of the TSS values (LOG), which reduces the positive skew commonly seen in the distribution of microbiome data. Lastly, we assessed normalization by TSS alone to remove differences in sequencing depth between samples as well as the effects of not using any normalization.

The compositional transformation methods were the most generalizable across non-linear models (Figure 4A), with median ILR F1 score of 74.3 (71.4–76.9) and MCC of 58.7 (53.6–63.1) and median CLR F1 score 74.2 (71.5–76.9) and MCC 58.5 (53.8–63.1). The compositional transformations were followed by the variance/distribution modifiers ARS (F1 score of 72.5 (69.8–75.8) and MCC of 56.0 (50.8–61.8)) and VST (F1 score 72.0 (65.9–75.0) and MCC 54.8 (44.1–59.9)). Lastly, TSS [F1 score 69.8 (63.8–73.9) and MCC 51.9 (40.5–58.7)] and LOG[F1 score 68.9 (64.1–73.5) and MCC 51.3 (40.4–58.6)] were consistently the lowest performing normalization. Furthermore, the compositional methods led to significantly better F1 score and MCC than the other normalization type (Figure 4B), whereas the variance/distribution modifiers and scaling method were significantly better than no normalization. These results indicate the importance of transformation methods which account for the compositional properties of microbiome data prior to model training.

Various approaches have been proposed to remove technical batch effects from next generation sequencing datasets, of which we selected three relevant to microbiome data. First, the zero-centering method aims to reduce batch effects by centering the mean of each feature within a batch to zero (Nygaard et al., 2016). Second, Meta-analysis Methods with a Uniform Pipeline for Heterogeneity in microbiome studies (MMUPHin) (Ma et al., 2020) implements an empirical Bayes’ approach to estimate and remove batch-specific parameters for each feature. Whereas, ComBat-seq implements a negative binomial regression model to estimate and correct batch effect parameters (Zhang et al., 2020).

We first compared the ability of these three methods to correct batch effects in our dataset. MMUPHin and ComBat-seq were provided the disease and sample type as biological covariates alongside the batch label, whereas zero-centering was blind to biological covariates and only provided the batch label. The performance was evaluated with the beta-diversity Aitchison distance between samples (Figure 5A,B) (Aitchison, 1982) and the Local Inverse Simpson Index (LISI) (Figure 5C) (Polański et al., 2020). LISI evaluates the local neighborhood of a sample with respect to the batch (iLISI) or disease label (dLISI). A higher value indicates the presence of samples with a greater variety of labels in the surrounding neighborhood. A good batch effect correction method will lead to an increased iLISI (improved mixing of batches) and a decreased dLISI (improved separation of disease labels).

Compared to no batch effect correction, MMUPHin improved separation of disease labels (1-dLISI of 0.968 vs 0.951, p < 0.0001), whereas it led to greater separation of batches in our dataset (iLISI of 0.027 vs 0.029, p < 0.0001) (Figure 5A–C). Conversely, ComBat-seq reduced separation of disease labels (1-dLISI of 0.941 vs 0.952, p < 0.0001) and improved mixing of batches (iLISI of 0.032 vs 0.029, p < 0.0001) compared to no batch effect correction (Figure 5A–C). The naive zero-centering method led to the greatest mixing of disease labels (1-dLISI of 0.927 vs 0.951, p < 0.0001) and greatest mixing of batches (iLISI of 0.082 vs 0.029, p < 0.0001) compared to no batch effect correction (Figure 5A–C). For the combined assessment of the dLISI and iLISI scores, we calculated the harmonic mean (F1 score) (Lin et al., 2019; Tran et al., 2020), where a higher mean indicates a better batch effect correction method. The F1 score indicated the performance order of batch effect correction methods was naive zero-centering (0.150), ComBat-seq (0.061), and MMUPHin (0.052), which was lower than no batch effect correction (0.057).

We assessed the effect of three batch effect correction methods on the classification performance of non-linear models. Two aspects of MMUPHin and ComBat-seq were accounted for to ensure a similar scenario of an implemented diagnostic test. First, the disease covariate is commonly provided alongside the batch label, which is unknown for a diagnostic test. Second, information is shared across batches in order to correct the batch effect estimates. For a diagnostic test, new test datasets would be continually received and have to undergo batch effect correction independent of the training dataset. Therefore, we implemented two variations to simulate the diagnostic scenario of obtaining a new dataset with unknown disease labels. The first (#1) method applied batch effect correction to the training and test sets separately providing the sample type (stool/biopsy) as the biological covariate. Whereas the second (#2) method only applied batch effect correction to the training set with the disease and sample type covariates provided (see Methods for detailed description). On the other hand, batches are independent in the zero-centering method allowing it to be implemented without the training and test sets affecting each other.

Of the six different batch effect correction variations, zero-centering and MMUPHin #1 were the most generalizable approach across the non-linear models. Zero-centering led to a significantly greater F1 score of [75.3 (70.3–79.2) vs 72.2 (69.8–75.0), p < 0.01] and MCC [61.6 (55.6–67.0) vs 55.1 (50.8–60.1), p < 0.001] than no batch correction (Figure 6A,B). MMUPHin #1 also led to a significantly greater F1 score of [74.7 (71.5–76.9) vs 72.2 (69.8–75.0), p < 0.05] and MCC [59.4 (53.9–63.1) vs 55.1 (50.8–60.1), p < 0.01] than no batch correction (Figure 6A,B). Whereas, ComBat-seq #1 displayed greater median F1 scores of 71.9 (64.5–73.5) and MCC of 54.5 (41.0–57.4) than no batch reduction, although not significant. Classification performance with MMUPHin #1 and ComBat-seq #1 was more consistent across the different non-linear models than zero-centering, likely due to the poor performance of the radial SVC with Zero-Centering. (Figure 6A). Whereas, MMUPHin #2 and ComBat-seq #2 were the least generalizable (Figure 6A) with significantly lower F1 score of 67.5 (61.5–71.1) (p < 0.0001) and 71.9 (64.5–73.5) (p < 0.05) and MCC of 47.0 (35.2–53.3) (p < 0.0001) and 54.5 (41.0–57.4) (p < 0.05) than no batch effect correction MCC of 54.6 (47.7–59.8), respectively. Therefore, the naive zero-centering MMUPHin #1 method for batch effect correction are the most generalizable approach for IBD classification with non-linear machine learning models.

The samples used to assess the performance of different combinations of normalizations or transformations, batch effect correction, and ML models were drawn from across sample collection methods (i.e., stool and biopsy) and patient demographics (i.e., paediatric and adult samples). While we did not set inclusion criteria for samples based on these differences, previous research has demonstrated distinct differences in microbiome composition between sample types and demographic groups (Walters et al., 2014; Kim et al., 2020; Radjabzadeh et al., 2020). For example, principal coordinate analysis (PCoA) with weighted UniFrac distance of (Durbán et al., 2011) and principal component analysis (PCA) of CLR-transformed taxonomic features indicated paired biopsy and stool samples from the same individual cluster separately (Mas-Lloret et al., 2020).

We compared the model performance for the sample and patient demographics for which we were able to acquire sufficient metadata and have been associated with microbiome alterations: sample type (biopsy vs stool), IBD type (CD vs UC), sex (Female vs Male), BMI (BMI <30 vs BMI >30), and age (Adult vs Pediatric). To assess the performance within each demographic, we included the predictions from taxonomic features (species, genus, OTU) with a compositional transformation method, zero-centering batch effect correction, and a non-linear ML model. Our analysis focused on the MCC performance metric as it is more robust to imbalanced label distribution (Chicco and Jurman, 2020), which occurred when the samples were grouped by the five metadata categories. A logistic regression function was used to assess changes in performance corresponding to each demographic while controlling for the other metadata (Table 3).

IBD classification performance was reduced for biopsy samples compared to stool samples, increased for samples from adult patients compared to pediatric patients, and decreased for samples from patients with BMI less than 30 compared to patients with BMI greater than 30. On the other hand, there was no difference in classification performance for females compared to males or for samples from patients with CD compared to patients with UC (Table 3). The metadata groups with different performance between the two categories coincided with those that are not equally represented in our dataset, highlighting the importance of accounting for different demographic groups in a microbiome based diagnostic test.

Our analysis identified the features, ML models, normalization or transformation methods, and batch effect correction methods which led to the most generalizable performance. In order to determine the best overall combination of features, data processing, and ML model we assessed the top three performing models (Table 4). The top two models consisted of the most generalizable individual components: taxonomic features (genus), non-linear model (XGBoost), compositional transformation (ILR or CLR), and zero-centering to remove batch effects. In addition, the third best model consisted of the lowest performing taxonomic features (OTU) and the lowest performing normalization (LOG), with zero-centering for batch effect correction. Overall, the combination of the most individual most generalizable methods led to the top two pipelines with the best classification performance.

For the three best pipeline combinations (Table 4), we assessed the classification performance on each dataset from the fifteen LODO iterations. Since some studies contained only a single label, either non-IBD, control or IBD, we assessed the performance using the classification accuracy metric (Figure 7). The three pipeline combinations had the lowest performance on a common subset of the fifteen datasets. The low performance datasets were enriched in samples from sample types and patient subgroups that we previously showed negatively affect classification performance (Table 3), such as biopsy samples (MUC, HMP, and GEVERSM), patients with BMI >30 (GLS or QIITA10342), or pediatric patients (GEVERSC). These results further highlight the importance of generating a diverse training set that is representative of the patient demographics the diagnostic model will encounter.

In addition to predicting health diagnoses, machine learning models can be used to identify biomarkers for disease by identifying features important for disease classification. We characterized the feature importance from the second-best overall data processing and ML model pipeline (Table 4). We did not analyze the feature importance of the best-performing model because the ILR normalized values no longer correspond to the input features thereby preventing interpretation of feature importance. For an XGBoost model, the importance corresponds to a feature’s contribution to the model’s decision during training, referred to as the gain value (Chen and Guestrin, 2016). We extracted the features’ gain values from each of the 15 LODO iterations, sorted by the mean of all iterations, and plotted the top fifteen features (Figure 8). In addition, we determined the change in abundance for each taxonomy to assess whether our dataset aligned with previous findings on changes of the microbiome in IBD.

Many of most important taxa are in the short chain fatty acid (SCFA) producing Clostridium XIVa/IV clusters, including bacteria from the Eubacterium, Coprococcus, Lachnospira, and Ruminiclostridium genera (Figure 8A). Aligning with previous studies, these bacteria were decreased, with the exception of Coprococcus 3, in IBD samples vs control samples in our dataset (Figure 8B) (Gevers et al., 2014; Nagao-Kitamoto and Kamada, 2017). Fusobacterium and Veillonellaceae genera, commonly increased in the gut microbiome of IBD patients, were also top contributors to the XGBoost classifier (Figure 8A,B) (Gevers et al., 2014; Glassner et al., 2020). In addition, the Prevotellaceae genus was the second most important feature, with the decreased abundance in IBD samples agreeing with previous studies showing decreased abundance in the gut microbiome of patients with CD and UC (Figure 8A,B) (Chen et al., 2014). Even though some taxa have a low fold change between IBD and non-IBD controls, XGBoost is able to find a cutoff value that contributes to the separation of labels and the final decision. XGBoost classifiers have the best potential for use as a diagnostic test due to their performance as well as their interpretability and utility in identifying disease biomarkers.