Inflammatory bowel disease biomarkers of human gut microbiota selected via different feature selection methods

Dataset and preprocessing steps

This study aims to develop a classification model to aid IBD diagnosis and to discover IBD-associated bio-markers using metagenomics data. In this respect, the raw microbiome DNA sequencing data of 148 IBD patients and 234 control samples were fetched from the MetaHit project (Oudah & Henschel, 2018; ERA000116). Prior to taxonomic analysis, we followed the standard quality analysis, which is proposed by the Human Microbiome Project (HMP) (Consortium et al., 2012; Young, 2011). According to that, firstly the duplicate reads were removed by a modified version of EstimateLibraryComplexity in Picard tool package. The quality trimming was performed using TrimBWAStyle script (Fass, 2010), which is also recommended by the HMP SOP. The human genome contamination removal step was skipped, as the provided data is already free of human contamination. Prior to taxonomic analysis, the reads smaller than 90 bp in length were filtered out.

To estimate the relative abundance of microbial taxa, we use the MetaPhlAn2 tool (Ditzler, Polikar & Rosen, 2015), which is a widely used method in the literature. MetaPhlAn2 firstly assigns reads to microbial clades using a set of clade-specific marker genes; and secondly, it estimates the relative abundance of microbial taxa based on the read coverage. In this study, each DNA sequence is assigned to its microbial species of origin (taxa) using MetaPhlAn2 taxonomic classification tool. The standard relative abundance normalization procedure is followed by dividing the read count of each taxonomic bin by the total number of the reads for a sample. Hence taxonomic abundances values are obtained as real numbers in the range of [0, 1], which sums up to 1 within each sample. Samples containing less than one million reads were discarded. Consequently, the microbial diversity (i.e. which microorganisms exist in what relative abundance) of the gut microbiome for each sample was revealed.

Microbiome sequencing data is categorized into disease states based on the associated metadata. In the original study that provides the metagenomics data (PRJEB2054; PRJEB2054), the healthy controls were selected among the close relatives of IBD patients and the individuals who had antibiotic treatment for at least 4 weeks before fecal sample collection was not included. IBD subjects were in clinical remission for at least 3 months, and had stable maintenance therapy with azathioprine or mesalazine. The patient demographics providing the country, gender, age, body-mass index are available in the original study (PRJEB2054).

At the end of the data preprocessing steps, the dataset included the relative abundance values for 3,302 different taxa for 382 samples. In the lower taxonomic levels (species), there are considerable variations in human gut microbial composition (Eckburg et al., 2005). Hence, species-level information is widely used in metagenomics classification problems. In this respect, the features of kingdom, phylum, class, order, family, genus were removed, and the remaining 1,455 features were considered. When the features of the same species for different strains were combined into one single feature, the final dataset included 1,331 species-level features. As shown in Fig. 1, the final IBD-associated metagenomics dataset is composed of the relative abundance values for 1,331 different species for 148 IBD patients and 234 healthy samples. This dataset is used to develop, train, test, and validate a machine learning model to diagnose IBD.

The methodology of this study can be summarized by the following steps: (i) feature selection to identify the most informative IBD-associated biomarkers; (ii) model construction to classify IBD patients and control samples, followed by the assessment of the generated models using extensive evaluation metrics; (iii) unsupervised learning to detect subgroups of IBD patients and control samples. Figure 2 shows the details of these three steps.

Feature selection

Feature selection is known as an effective preprocessing tool for machine learning problems (Tang, Alelyani & Liu, 2014). However, deciding between the growing numbers of feature selection methods is challenging. Different feature selection methods have their advantages/disadvantages as discussed in different reviews (Chandrashekar & Sahin, 2014; Khalid, Khalil & Nasreen, 2014; Li, Li & Liu, 2017). In a previous work (Manikandan & Abirami, 2021), several state-of-the-art feature selection methods were reviewed in terms of their ability to overcome common problems such as data nonlinearity, correlation and redundancy, noise in the target class, noise in the input features, and having a much higher number of features much higher than the number of samples. Feature selection methods have different usages in the biomedical domain, as presented in Remeseiro & Bolon-Canedo (2019). Moreover, there are different kinds of feature selections, i.e. filter, wrapper, and embedded methods such as (Yousef et al., 2020). Recent feature selection methods make use of the biological knowledge, which is embedded in the machine learning algorithm (Yousef, Sayıcı& Bakir-Gungor, 2021; Yousef, Abdallah & Allmer, 2019; Yousef et al., 2021). Applications of biological domain knowledge based feature selection methods for gene expression data can be found in: Yousef, Kumar & Bakir-Gungor (2021).

In metagenomics studies, the number of predictors (number of taxa or features) is much more than the number of observations (samples). This phenomenon is known as the curse of dimensionality. In this respect, some metagenomics studies focus on the feature selection process rather than focusing on classification. It has been discussed in the literature that although the feature selection process in metagenome-based disease prediction problem is relatively less studied, this area could be as critical as the choice of the classification method (LaPierre et al., 2019). The feature selection process in metagenome-based disease prediction problem could improve the elucidation of the disease mechanisms. Hence, further research in this direction is provoked. To reduce the number of taxa (species or features), in other words, to select informative features, min Redundancy Max Relevance (mRMR) (Ding & Peng, 2005), Lasso (Tibshirani, 1996), Elastic Net (Zou & Hastie, 2005) and iterative sure select algorithm (Duvallet et al., 2017) have been extensively applied previously. Another feature selection method called Fizzy challenges to detect important microbes or functional elements for downstream analysis using classification algorithms (Ditzler et al., 2015). Oudah and Henschel introduced a competing taxonomy-aware feature selection method (Oudah & Henschel, 2018). (Bakir-Gungor et al., 2021) applied CMIM (Brown et al., 2012), FCBF (Fleuret, 2004), mRMR (Ding & Peng, 2005), and Select K best (SKB) (Pedregosa et al., 2011) on Type 2 diabetes associated metagenomics dataset and obtained high-performance metrics. Although those feature selection approaches are well studied in different domains, they are just getting attention in this domain. In a recent review paper (Marcos-Zambrano et al., 2021), some of these methods are reported to achieve good results in human microbiome studies. However, as noted by this review paper, for metagenomics studies, there is no consensus on which feature selection method should be used. In this study, we proposed that multiple feature selection methods should be used and the final features are determined by getting the intersection of the features selected by different methods. Conditional Mutual Information Maximization (CMIM) (Brown et al., 2012), Fast Correlation Based Filter (FCBF) (Fleuret, 2004), Min Redundancy Max Relevance (mRMR) (Ding & Peng, 2005), SKB, Information Gain (IG) and Extreme Gradient Boosting (XGBoost) (Chen & Guestrin, 2016) feature selection algorithms are applied on the metagenomics dataset.

Information Gain (IG) is one of the filter-based feature selection techniques, which leverages the concept of entropy (Gray, 2011; Kent, 1983). More specifically, a weight for each feature is calculated by analyzing the reduction in the class entropy (to which extent the uncertainty in the class prediction drops) when the value of that feature is known. IG is one of the widely used feature selection techniques, and it has applications in different domains (Bolón-Canedo & Alonso-Betanzos, 2019).

mRMR is another mutual information (MI) based method, and it chooses features according to the maximal statistical dependency criteria. Since the implementation of the maximal dependency condition is laborious, mRmR works as an approximation technique to maximize the dependency between the joint distribution of the selected features and the classification variable. In other words, mRMR (Ding & Peng, 2005) attempts to choose the features that have the minimum correlation between themselves (defined as min redundancy); and the maximum correlation with a class to predict (known as max relevance).

CMIM (Brown et al., 2012) first ranks the features according to their conditional entropy and MI with the class to predict; and then selects the feature if it carries additional information. Likewise, FCBF (Fleuret, 2004) ranks the features based on their MI with the class to predict; and then removes the features whose MI is less than a predefined threshold. This method aims to detect relevant features, as well as redundancy among relevant features via calculating feature-class and feature-feature correlations. FCBF picks features in a classifier-independent manner. It chooses the features with high correlation with the target variable, but low correlation with other variables. As a correlation measure, FCBF uses Symmetrical Uncertainty, which harmonizes Shannon entropy and IG. Select K best method scores the features against the class label using a function, and it selects the features having k highest score (Pedregosa et al., 2011). In XGBoost (Chen & Guestrin, 2016) feature selection, the more an attribute is used to make key decisions with decision trees, the higher relative importance it gets. In this study, CMIM, mRMR, FCBF, SKB, IG and XGBoost feature selection methods are implemented in Python 3, using the skfeature and sklearn libraries.

Model construction

Applying supervised learning to the human gut microbiome is useful to detect subsets of microorganisms that are substantially discriminative. Accordingly, one can generate prediction models that can precisely classify unlabeled samples. Recently, Marcos-Zambrano et al. (2021) reviewed machine learning (ML) applications for microbiome studies via analyzing 89 papers. They reported that the most common supervised learning algorithms that were used for microbiome analysis were Random Forest (RF), Support Vector Machines (SVM), Logistic Regression (LR) and k-NN (k nearest neighbor). They concluded that multiple factors need to be taken into account while selecting the ML algorithm (i.e., number of observations, number of features, data type, data quality etc.), and they recommend applying and evaluating more than one method and selecting the one with the highest performance value. Another recent study by Wang & Liu (2020) compared the performances of two ensemble methods (RF, eXtreme Gradient Boosting decision trees (XGBoost)), and two traditional methods (Elastic net and SVM) on 29 human microbiome benchmark datasets. They find that in a few benchmark datasets XGBoost outperforms all other methods; and in the remaining benchmark datasets, XGBoost, RF and Elastic net were comparable. Along this line, in this study, a range of machine learning models were constructed to discriminate IBD samples from controls, using different classification algorithms (RF, Decision Tree, Logitboost, Adaboost, SVM, and stacking ensemble classifiers (i.e., an ensemble of SVM with kNN, an ensemble of the Logitboost with kNN). The working principles of these classification algorithms can be summarized as follows:

Decision trees are the base learners of RF. Each tree is a non-linear model and each model is created with many linear boundaries (Liaw & Wiener, 2002). During the splitting (training) procedure, RF randomly chooses bootstrapped samples from the original data, and randomly selected subsets of features are used to assess the quality of the model. Finally, RF merges several decision trees into a single ensemble model and performs prediction via combining the predictions of individual trees.

In kNN classification, for each sample, the class labels of k training samples nearest to that point (neighbor samples) are investigated, and the majority of those class labels are assigned to the sample. To specify the neighborhood, the most widely used distance metrics are correlation coefficients and Euclidean distance.

SVMs attempt to find out a decision boundary between the classes. This decision boundary helps to assure the maximum feasible distance or margin between the samples that are closest to this decision boundary. The decision boundary is defined via learning from the samples that are closest to the boundary, and hence these samples are called support vectors. In some cases, a linear separation between the classes is not feasible in the original feature space. In such cases, SVM utilizes the kernel trick to determine the decision boundary in a higher-dimensional space (Cortes, 1995).

Boosting was initially proposed (Schapire, 1990) to unite multiple weak classifiers and hence to elevate the classification performance. Freund & Schapire (1997) proposed a more practical and capable boosting algorithm called AdaBoost. AdaBoost, which stands for Adaptive Boosting, is one of the most widely used algorithms with many applications in numerous fields. AdaBoost focuses on the difficult samples by assigning higher weights on them that could not be properly classified with the previous weak classifiers. AdaBoost is capable of reducing the training errors exponentially fast as long as the weak classifiers produce just better results than the random case (Freund & Schapire, 1997). It was reported that AdaBoost had very good generalization capability. Still, like most other classifiers, when dealing with very noisy data, AdaBoost suffers from the overfitting problem (Rätsch, Onoda & Müller, 2001). To overcome this situation, Schapire, Freund et al. (1995) proposed another algorithm called LogitBoost, which could reduce training errors linearly, and hence better generalization. It was designed to decrease the bias and the variance. Logitboost was originally utilized for integrating non-complex classifiers to improve performance in classification. This technique optimizes the multinomial likelihood, which makes it easy to apply in multiclass problems. The derivation of the LogitBoost algorithm can be done by applying logistic regression to the AdaBoost generalized additive model.

Throughout the experiments of this study, 100-fold MCCV is used. In MCCV, some part of the data is randomly chosen (without replacement) as the training set and the rest is used as the test set (Xu & Liang, 2001). This procedure is iterated several times, and hence new training and test sets are randomly generated in each iteration. 90 % is chosen for training and 10 % is chosen for testing. As shown in the workflow Fig. 2, the feature selection methods are applied on the training set. Classification methods are tested on the testing part. The Konstanz Information Miner (KNIME) platform (Berthold et al., 2009) is used for the implementation of the methodology. H20 library in KNIME and Python scikit-learn library (Pedregosa et al., 2011) are also used.

Model performance evaluation

Prediction performances of the generated models were evaluated by using accuracy, F1 Score, and AUC measures. Accuracy is a widely used performance evaluation metric and a reliable measure for balanced datasets. Since there is an uneven class distribution in the dataset of this study, other metrics such as the F1 score and AUC have been utilized to evaluate the performance of the generated models. Precision is defined as the percentage of the predicted cases that are actual cases. On the other hand, recall represents the percentage of actual cases that are correctly identified by the classifier. That is to say that recall determines the rate of falsely predicting healthy. Precision denotes the rate of falsely predicting disease. F1 score is calculated as the harmonic mean of precision and recall. Among different performance metrics, F1 score is a good choice when someone seeks a balance between precision and recall, and when there is an uneven class distribution (a large number of actual negatives). Several classifiers deliver the probability values for their predictions, which can be considered as their confidence values regarding the prediction. The AUC utilizes this information to recapitulate the false prediction rate at different confidence levels. AUC is commonly used as a summary measure of diagnostic accuracy. In real-life examples, there is an overlap between the test results of positive and negative examples. AUC shows how the recall vs. precision relationship changes as the threshold or the cut-off value for identifying a positive is changed. All performance results presented in this study refer to the average of 100-fold MCCV.

Unsupervised learning

One of the goals defined for precision medicine is stratifying patients based on disease subtypes, disease progression, and response to treatments via analyzing molecular profiling data of the patients (Korcsmaros, Schneider & Superti-Furga, 2017). Precision medicine provides big promise to enhance the course of care for IBD patients since it offers the most effective therapy while minimizing the side effects (Weersma et al., 2018). In this respect, studying the relationships between the samples and the relative abundance values of the species can help to detect subgroups of IBD patients; and hence enlighten IBD development mechanisms. To analyze and visualize these relationships, in this study, three unsupervised learning approaches (K-means clustering, principal component analysis (PCA), and hierarchical clustering) were used for the following purposes. Firstly, this study attempts to answer whether there could be any direct relationship between specific species (or a group of species) and IBD subgroups. In this respect, to explore the subgroups of IBD patients and subgroups of healthy samples, K-means clustering is used. K-means (Steinley & Brusco, 2007) is an unsupervised learning algorithm that groups samples by minimizing the distance inside the clusters and maximizing the distance between the clusters. The Euclidean distance metric and Elbow method (Bonaros, 2019) are utilized in this study to find out the optimum number of clusters. In this method, the point where the decline in the error slows down indicates the optimum number of clusters. t-distributed Stochastic Neighbor Embedding (t-SNE) is a nonlinear dimension reduction technique. It is a commonly used graphical approach to guide clustering methods such as K-means in terms of deciding the number of clusters and cluster memberships. In this study, t-SNE is employed for visualizing the identified clusters (subgroups of IBD patients and healthy samples).

Secondly, to observe whether the relative abundance values of the species can induce the formation of two clusters representing IBD patients and controls, IBD-associated metagenomics data were modeled using principal component analysis (PCA). PCA is a dimensionality reduction algorithm that converts a high dimensional space (where each dimension corresponds to a species) to a lower-dimensional space (usually 2D or 3D). In unsupervised machine learning, the class labels (here the diagnosis of IBD, in other words cases and controls) are hidden from the model, leaving the algorithm to impose the most relevant strata.

Thirdly, to better visualize the relationship between the relative abundance values of the selected species and IBD patients/controls, hierarchical clustering was performed using Euclidean distance and Ward variance minimization algorithm as the linkage method. With this analysis, the aim is to reveal the presence of distinct subgroups of IBD patients, corresponding to the ones having complex patterns of microbial species.

Feature selection and classification

Since the dysbiosis between the mucosal immune system and gut microbiota is a hallmark of IBD pathogenesis, the microbiome has the potential to harbor species with biomarker capacity. To detect those species with biomarker capacity, an IBD-associated metagenomics dataset, which includes the relative abundance values for 1,331 different species for 382 samples (following the preprocessing steps in Methodology section), is analyzed. This research effort attempted to remove the irrelevant and redundant features (species) using six different feature selection strategies (FCBF, CMIM, mRMR, SKB, IG and XGBoost). For each feature selection method, the top 100 features are investigated. To evaluate the effects of different classification methods, Decision Tree, Random Forest, LogitBoost, AdaBoost, an ensemble of SVM with kNN, and an ensemble of the Logitboost with kNN are considered. The parameters of c and gamma are optimized for SVM; number of trees is optimized for Random Forest, and n estimators is optimized for Logitboost and Adaboost. By using several metrics as described in the Methods section, the performances of different classifiers are compared using (i) all features (without feature selection); (ii) top 100 features selected using CMIM, mRMR, FCBF, SKB, IG and XGBoost (as presented in Table S1). As shown in the same table, SKB, IG and XGBoost feature selection methods resulted in high accuracy, F1 score, sensitivity, recall values, and high AUC scores for different classifiers throughout the experiments with IBD-associated metagenomics data. It can be observed from the same table that CMIM, FCBF, and MRMR feature selection methods showed signs of poor fitting across all models with low accuracy and high recall. To find the most relevant and informative features, for each feature, the scaled importance values were calculated for three selected feature selection methods (SKB, IG and XGBoost) across different classifiers (presented in Table S2). To eliminate the lowest ranking features among the top 100, the scaled importance value cutoff of 0.5 is applied. As shown in Fig. 3, this procedure generated 23, 57, 96 selected features in SKB, IG, and XGBoost feature selection methods, respectively. 14 of those features were commonly identified in all three-feature selection methods.

As shown in Fig. 4 and Table S1, throughout the experiments using 14 features, it is observed that RF, LogitBoost, and Adaboost classifiers resulted in the top three performance results as compared with other classifiers. Among these three classifiers, using the 14 selected features, the RF classifier generated higher performance results for the above-mentioned promising feature selection methods (SKB, IG, and XGBoost). In RF, the interpretation of the tree model is simple and the model can be easily transformed into a ruleset. Also, as noted in Marcos-Zambrano et al. (2021)s comprehensive review, RF is one of the widely used algorithms in human microbiome studies. As shown in Table S1, even if all 1,331 features are used with Logitboost or with Adaboost, RF with 14 features performs as well as Logitboost and Adaboost, even slightly better. It should be noted that one of the goals of this study is to identify IBD-associated biomarkers of human gut microbiota. Therefore, obtaining satisfactory performance metrics for diagnosing IBD using a small number of selected species is important. The experiments conducted in this study showed that, on Exploration Cohort, the generated RF model using only 14 features resulted in adequate diagnostic accuracy (as shown in Fig. 4 and Table S1). For all these reasons, throughout the rest of the paper, we focused on the results obtained using the RF classifier.

When all 1,331 features are used (without applying feature selection), the generated RF model yielded 0.85 F1 score, 0.92 AUC, and 87 % accuracy on Exploration Cohort, as shown in Fig. 5 and Table S1. By only using the 14 features that are commonly identified in three promising feature selection methods, 0.85 F1-score, 0.93 AUC, and 88% accuracy metrics were obtained using the RF model. Compared with the performance values of the RF model using all features, the model with those selected 14 features performed 1 % higher in terms of accuracy and AUC metrics; 5 % higher in terms of specificity and precision metrics, as shown in Fig. 5 and Table S1. Also, using only those 14 species, the generated RF model yielded the same F1 score (0.85) as the one obtained using all features.

Validation on external data

The performance of the proposed method is also evaluated using an external dataset. A multicenter gut metagenome dataset (Project accession: PRJEB1220) with samples collected from Denmark and Spain, containing Ulcerative colitis (UC) and Crohns disease (CD) cases as well as healthy controls were considered for validation. A random subset of samples of 50 UC patients, 50 CD patients (in total 100 IBD patients) and 100 healthy controls, containing more than 1 Gbp of sequencing reads were subject to taxonomic analysis. MetaPhlAn 3.0 is run using default parameters and the relative abundances of all detected taxonomies are determined. The same preprocessing protocol that is used for the Exploration Cohort (accession: PRJEB2054), and that is presented in the Methods section; is followed for the independent test data (Project accession: PRJEB1220). While the Exploration data included 1331 species, the independent test data included 488 species as features. For the validation data, it is ensured that the species or genus names match, and the same reference database is used. Among the 14 potential taxonomic biomarkers, 10 species (Table S3) are found in the validation dataset. As shown in Table 1, throughout the experiments on the validation data using the selected species as features, it is observed that RF, LogitBoost, and Adaboost classifiers resulted in the top three performance results as compared with other classifiers. Among these three classifiers, RF performed slightly better (0.86 accuracy) than the other two classifiers (0.84, 0.82 accuracy). It is worth noting that the same trend was also observed for the Exploration Cohort, as shown in Fig. 4 and Table S1. Additionally, using 100-fold MCCV and the same set of classifiers, which are utilized in the experimentation with the Exploration Cohort, the performance of the models using 10 selected features are compared against the models using 10 random features. As shown in Table 1, on the validation data, the generated models resulted in higher performance metrics when 10 selected features (species) are used, as compared to the randomly generated 10 features. Especially, in terms of specificity, one can easily observe the sharp decrease to 0.35 when 10 random features are tested, as compared to the obtained specificity value of 0.85 with 10 selected species when the RF classifier is applied on the validation data. The same trend of significant decrease in terms of specificity is observed in the results of all other tested classifiers when 10 random features are tested (Table 1).

As shown in Table 2, while 0.85 ± 0.05 accuracy is obtained on the Exploration Cohort, 0.86 ± 0.07 accuracy is obtained on the Validation Cohort using 10 species, RF classifier and 100-fold MCCV. With the addition of standard deviation values, the results seem quite similar and no significant difference is observed. Moreover, 0.87 ± 0.06 accuracy is obtained using 14 species on the Exploration Cohort, where four of these features (Subdoligranulum_unclassified, Ruminococcus_obeum, Lachnospiraceae_bacterium_1_1_57FAA, Lachnospiraceae_bacterium_2_1_58FAA) are not found in the Validation Cohort. Comparing 0.85 ± 0.05 accuracy value obtained using 10 species with 0.87 ± 0.06 accuracy value obtained using 14 species on the Exploration Cohort, we noticed that 4 additional features are not so necessary for improving the results.

As compared to using all 1,331 features (0.87 ± 0.05 accuracy), those selected 10 features yielded in very similar performance in terms of accuracy and other metrics. On the Validation Cohort, it is also observed that satisfactory performance metrics are obtained using these selected species. That is to say that 10 features are sufficient to accurately classify IBD, compared with using all 1,331 features. For the diagnosis, evaluating the amounts of fewer species is more effective in terms of cost and time. Hence, those 10 features (species) shown in Table S3 are proposed as potential taxonomic biomarkers for IBD. These potential taxonomic biomarkers can be adapted to clinics for facilitating the diagnosis of IBD. If adapted to the clinic, using those 10 features, IBD diagnosis could be performed with 0.85 ± 0.05 accuracy. We thus conclude that the analysis of these 10 potential taxonomic biomarkers in gut microbiota might be used as an auxiliary diagnostic tool for suspected IBD patients.

Feature importance and their correlations

Regarding metagenome-based disease prediction, the most informative features correspond to the microbes that contribute highest to the disease prediction, which strengthens the interpretability of the model (Pasolli et al., 2016; Reiman, Metwally & Dai, 2018; Rahman & Rangwala, 2018; Duvallet et al., 2017). The identification of critical species, the ones having a key role in IBD disease development, can constitute new targets for the development of probiotics to correct microbiota aberrations (Armour et al., 2019; Gueimonde & Collado, 2012). To this end, the feature importance scores of the 14 selected species (as shown in Fig. S1) are investigated. In Fig. S1, while the Y-axis corresponds to the detected species, X-axis corresponds to their relative importance.Some of the identified species were previously reported in the literature as microbiome-associated factors for IBD. Among the top three species that are found as potential IBD biomarkers, Bacteroides xylanisolvens (Ulsemer et al., 2012) is considered as candidate next-generation probiotics, promoting gut health. Bacteroides xylanisolvens is one of the short-chain fatty acid (SCFA)producing bacteria. SCFAs are colonotrophic nutrients and they are immunoregulatory molecules (Puertollano, Kolida & Yaqoob, 2014) that may reduce pro-inflammatory cues within the gut environment. In this study, Bifidobacterium bifidum is identified as the fourth significant species. For the treatment of IBD, a recent review by Jakubczyk, Leszczyńska & Górska (2020) showed the beneficial effect of probiotics including Bifidobacterium bifidum PRL2010 and Bifidobacterium bifidum 231 strains, which are potentially involved in the etiology of IBD. Another research group demonstrated that the supplementation with Bifidobacterium bifidum causes a significant increase in IL-10 level and decrease in IL-1 levels in the colon sections, which confirmed the anti-inflammatory effect (Kumar et al., 2017). This observation verifies the anti-inflammatory effect and hence, supports the modulatory role of Bifidobacterium bifidum in terms of decreasing the inflammation, as well as the clinical symptoms of colitis. With regard to UC, the strains of Bifidobacterium bifidum are considered as promising in sustaining remission (Kato et al., 2004; Kruis et al., 2004). The abundance of Lachnospiraceae bacterium, which is identified in the top 5 important species list, has also been observed previously in IBD cohorts (Nagao-Kitamoto & Kamada, 2017).

To analyze the pairwise correlations of these selected species, SPARCC (Friedman & Alm, 2012), which is able to estimate correlation values from compositional data, is utilized. As noted by Freilich et al. (2018), SPARCC is one of the rigid algorithms that can be used to evaluate correlation in microbiome datasets. SPARCC assumes a sparse data matrix, and the ϕ (Lovell et al., 2015) and ϱ (Erb & Notredame, 2016) metrics (the published versions of which required a non-sparse matrix). These relations were illustrated in Fig. S2 using a heatmap. Between any pair of the identified species, no significant correlation is observed.

Grouping control and IBD samples via K-means clustering algorithm

IBD patients are characterized by genetically and clinically defined subgroups, which have very specific microbial compositions and functions. To this end, analysis of patient microbial metabolomic profiles could be utilized as a predictive clinical tool, which provides a foundation for personalized microbiome-based therapies (Banfi et al., 2021). This research effort also investigated whether some subgroups of IBD patients have a direct relationship with some species. In this respect, a K-means clustering algorithm is used in this study to subgroup samples. As shown in Fig. S3, three subgroups among IBD samples and four subgroups among controls were discovered. t-SNE method was employed for visualizing the identified clusters (subgroups of IBD patients and healthy samples). To this end, based on the relative abundance values of the 14 identified species, two-dimensional t-SNE maps were generated separately for (i) IBD patients subgroups, and (ii) healthy sample subgroups (Fig. 6). Cluster-colored tSNE plots were visually inspected. As shown in Fig. 6, the subgroups of IBD patients and the subgroups of healthy samples were distinct.

Figure 7 illustrates the relative abundance values of the identified species in each of these subgroups. In Fig. 8, The presence of the six selected species was displayed more in detail for each IBD subgroup and controls. It can be concluded from both Figs. 7 and 8 that even though the samples were divided into subgroups, a single species has no direct effect on the development of IBD. Nevertheless, there are a few observations that one can make: (i) Porphyromonas asaccharolytica (shown in Fig. 8B and with grey in Fig. 7) and Peptostreptococcus anaerobius (shown in Fig. 8F and with dark blue in Fig. 7) is observed in all IBD subgroups, and not in the control subgroups; (ii) Bifidobacterium bifidum (shown in Fig. 8A and in pink in Fig. 7) is mainly observed in all IBD subgroups, it is present in very small amounts in the control subgroups; (iii) although Fig. 7 indicates that Eubacterium hallii (shown in orange) is found in all subgroups, the zoomed-in view in Fig. 8C clarifies that the relative abundance values of Eubacterium hallii is higher in IBD subgroups, compared to control subgroups (except for control 3 subgroup including 3 % of the control samples); (iv) a similar observation can be made for Lachnospiraceae bacterium 1_1_57FAA (shown in Fig. 8E and in light pink in Fig. 7) and for Dorea formicigenerans (shown in Fig. 8D and with brown in Fig. 7). It can be concluded from Fig. 7 and Fig. 8 that the relative abundance values of Bifidobacterium bifidum and Porphyromonas asaccharolytica have a role in distinguishing subgroups of IBD patients. Compared to other two subgroups of IBD patients, in the IBD_01 subgroup (including 70% of patients), the average relative abundance values are lower for Bifidobacterium bifidum, and higher for Porphyromonas asaccharolytica.

As a result of the IBD patient subgroup identification steps, the enrichment of Bifidobacterium (Wang et al., 2014) and Peptostreptococcus anaerobius species in IBD patients are in accordance with previous literature(presented in detail in the Discussion section), supporting the stratification analysis performed in this study. In the future, this information could be used to stratify IBD patients more precisely and to offer more effective treatments.

Discrimination of data via principal component analysis (PCA)

To invesigate whether the samples can be divided into two clusters representing control and IBD samples, the first two and three principal components of metagenomic data were obtained using PCA as an unsupervised learning approach. Figure 9 shows that a better separation is observed between control and IBD classes when PCA is performed with the 14 selected features (Figs. 9B, 9D), as compared to using 1,331 species (Figs. 9A, 9C). Figure 9 also points out that compared to the separation using the top two principal components (Figs. 9C and 9D), the third principal component contributes to the separation of controls vs IBD samples (Figs. 9A and 9B). An interactive 3D plot for Figs. 9A, 9B, 9C, and 9D is provided as supplementary material.

However, it is observed that the new feature space, reduced via PCA, does not have a significant contribution in terms of classification performance. It is also important to note that since PCA maps the data into a new feature space, the original feature information is lost during this process. Thus, PCA analysis does not allow for biomarker discovery, as the species information is no longer represented in the new mapped feature space.

Hierarchical clustering of reduced IBD dataset

To better visualize the relationship between the top 14 selected species and the samples, a hierarchical clustering analysis is also conducted. The heatmap in Fig. 10 is obtained using all samples and top-scoring 14 species. Colors represent raw z-scores for samples. While the black color indicates relative abundance values just around the mean, the lighter colors denote the relative abundance values of 1 to 4 standard deviations above the mean. The first column specifies class labels, i.e., IBD patients and healthy samples are shown in red and blue, respectively. The areas that are restricted with red boxes suggest differential relative abundance values for the corresponding species in the associated subgroup. For example, while the red boxes in the 7th and 12th columns (sp9, sp6) indicate excessive levels of Lachnospiraceae bacterium_1_1_57FAA and Ruminococcus obeum in the corresponding IBD subgroups; the red boxes in the 1st, 2nd and 3rd columns (sp2, sp12, sp13) indicate excessive levels of Bacteroides xylanisolvens, Ruminococcus bromii and Subdoligranulum unclassified in the corresponding healthy subgroups. The abundance of Lachnospiraceae bacterium has been observed previously in IBD cohorts (Nagao-Kitamoto & Kamada, 2017). Kang et al. (2010) demonstrated that compared with the healthy subjects, Ruminococcus bromii was observed five times less in the fecal samples of CD patients. For intestinal microorganisms, polysaccharides and dietary cellulose are significant energy sources. R. bromii can degrade some other complex polysaccharides (involving xylan and starch) and the degradation products are known to affect the gut microbial community, as well as human health (Wang & Liu, 2020).