Analysis of blood microbiome dysbiosis in pulmonary sarcoidosis by decision tree model

Abstract

Pulmonary sarcoidosis is a complex inflammatory disease characterized by granulomas in the lung tissue, leading to breathing difficulties and chest pain. Its etiology remains not fully understood, with factors such as allergies, autoimmune responses and genetics playing a role. This study explores the potential of blood microbiome dysbiosis, defined as an imbalance in the microbial ecosystem, as a missing piece of the puzzle in understanding the etiology of the disease. Our objective was to apply a decision-tree supervised machine learning hierarchical model to distinguish potential patterns of microbiome dysbiosis in blood samples from patients with pulmonary sarcoidosis as compared to healthy age-matched controls. Blood microbiome analysis, being individually-specific and stable, offers a unique perspective. Utilizing 16S rRNA gene amplicon sequencing, we analyzed the blood microbiome composition characterized by non-normally distributed and sparse data. Because of the rarity of the disease in Bulgaria, we studied a relatively small patient group, n = 7. The findings were compared to 21 healthy age-matched controls. Bioinformatics and statistical analysis play a pivotal role in microbiome analysis, especially when discerning associations between taxonomic composition and disorders such as pulmonary sarcoidosis. By analyzing the microbial diversity, we identified alterations in the blood microbiome composition between healthy individuals and those with sarcoidosis, which potentially may trigger the disease. Advanced machine learning techniques provided additional power to the analysis, that might be overlooked by the usual group statistics, confirming the differentiation of the diversity within the studied microbiome.

Keywords:

• Pulmonary sarcoidosis
• blood microbiome dysbiosis
• classification
• decision tree model

Introduction

Pulmonary sarcoidosis, is an inflammatory disease that has long been a subject of medical concern. It is characterized by the formation of granulomas in the lung tissue and lymph nodes a blood microbiome and may lead to breathing problems, chest pain, wheezing, and significant impairment of the patient’s quality of life [1]. The complete evaluation of the etiology of sarcoidosis is essential, and the potential role of pulmonary microbiome dysbiosis [2–4] and blood microbiome dysbiosis [5, 6] is emerging as a significant area of interest.

Dysbiosis refers to the imbalance or maladaptation of the microbial ecosystem within the body, which can lead to various health issues [7]. Among factors such as allergies, autoimmune responses, and genetic predisposition [8], the dysbiosis of the pulmonary and blood microbiome might be the missing piece of the puzzle in unravelling the etiology of pulmonary sarcoidosis.

Blood as a medium is capable of easily translocating microbial species across different body systems, such as the lung-gut axis and gut-brain axis [9, 10]. The composition and abundance of the blood microbiome are highly individually specific [6] and dynamically stable, not influenced dramatically by daily factors, such as diet, as compared to the gut microbiome [11]. This makes it highly probable for the blood microbiome to reflect and be involved in the etiology of systemic diseases [12], such as pulmonary sarcoidosis.

Microbiome analysis is typically based on metagenomic sequencing, either amplicon or whole-genome-based. These methods provide complex microbiome taxonomy descriptions and produce abundance composition datasets, characterized by being non-normally distributed, non-balanced, heterogenous, sparse, stratified and hierarchically organized [13]. The necessary data analytical approaches are a mix of machine learning, bioinformatics and statistics, aiming not only toward the analysis per se but also for a reasonable validation [14, 15].

Pulmonary sarcoidosis is a complex and relatively rare disease in Bulgaria with an incidence of less than 10 cases per 100,000 people [16]. It is difficult to collect a large number of patients, which is also a challenge in data analysis. Thus, throughout the study period, only a small group of eligible patients potentially met the full extent of predefined inclusion/exclusion criteria.

Machine learning methods, with their ability to find meaningful information below the level of noise, i.e. random data variations, can address these challenges [17]. Decision tree methods are particularly interesting in the context of biostatistics upgrade, as they can deal with complex data, are not dependent on data distribution, and can work with smaller sample sizes, among other benefits [14].

In our previous study [18], we applied a more conventional aproach to characterize the blood microbiome dysbiosis in sarcoidosis patients. Here, the aim was to reanalize the patterns of microbiome dysbiosis by applying a decision tree model. The meta-analysis in this study searched for novel taxonomic compositional patterns in blood microbiomes, aiming to uncover potential associations of specific microbial phyla and genera to the etiology of sarcoidosis.

Subjects and methods

Ethical considerations

The study was approved by the Institutional review board/Institutional ethics committee (IRB/IEC) of the National Center of Infectious and Parasitic Diseases at Sofia, Bulgaria (Number IRB00006384; Protocol Number 3/27.07.2020). All patients and healthy volunteers signed informed consent forms.

Study participants and sample collection

The study analyzed the blood microbiome data of patients diagnosed with pulmonary sarcoidosis and 21 healthy age-matched controls [18]. Between November 2020 and November 2022, 7 patients diagnosed with pulmonary sarcoidosis stage two were enlisted in the study. The seven patients (mean age ± SD, 48.33 ± 19.14; 4 males) were diagnosed positively with stage two pulmonary sarcoidosis [18]. None of the patients were receiving antibiotics or other sarcoidosis-related medicinal treatment at the time of the sample collection. The healthy controls included 11 males and 10 females (mean age ± SD, 49.21 ± 17.45). Blood samples of the patients and 21 clinically healthy individuals were collected in 5 mL Vacutainer tubes with K₃EDTA as an anticoagulant (Vacutainer K3E, BD, USA).

DNA extraction and microbiome analysis

DNA extraction: Three mililiters of whole blood were lysed in 10 mL of 0.2 µm filtered and autoclaved dH₂O with DNase I 1 U/mL for 1 h at 37 ⁰C. DNase I was added to eliminate cell-free and human DNA. Microbial cells were collected after centrifugation at 4000 rpm for 20 min (MPW M-diagnostics, cat. 102MD/V56, rotor cat. 12139; Med. Instruments, Warsaw, Poland). The microbial pellet was washed twice with 10 mL dH₂O. Each wash was followed by centrifugation. The cell pellet was resuspended in 1 mL of lysis buffer (500 mmol/L NaCl, 50 mmol/L Tris-HCl, pH 8.0, 50 mmol/L EDTA, and 4% sodium dodecyl sulphate). The cell suspension was subjected three times to vigorous vortexing or optional homogenization with silica/zirconium 0.1/0.3 mm beads (cat. 11079101z and 11079103z, Biospec Products, Bartlesville, OK, USA) on a microtube homogenizer (cat. D1030-E, Benchmark Scientific Inc., Sayreville NJ, USA) for 3 min at 4000 rpm. After the microbial cell lysis, DNA isolation was performed according to the previously described procedure [6,19,20]. The extracted DNA was resuspended in 100 µL of sterile DNA/RNA-free dH2O (cat. W4502, Sigma-Aldrich, Burlington, VT, USA). DNA typical yield was >150 ng/µL with a 260/280 nm ratio > 1.7. 16S rRNA amplicon sequencing was performed in-house. The V4 hypervariable region of the 16S rRNA gene was sequenced with NEXTFLEX® 16S V4 Amplicon-Seq Kit 2.0 for Illumina Platforms, cat. 4203-03, PerkinElmer, Inc. Waltham, MA, USA. Sequencing was performed by Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA, USA). Sequencing data were subjected to quality control measures to ensure accuracy and reliability. The Illumina software MiSeq Reporter version 2.5.1.3 (Illumina Inc., San Diego, CA, USA) and the Illumina Sequence Analysis Viewer version 2.1.8 (Illumina Inc., San Diego, CA, USA) were used for imaging and evaluation of the sequencing run performance. Primary data analysis and quality control (QC), including signal processing and de-multiplexing, were performed by the MiSeq Reporter. The Trimmomatic tool was used to trim low quality bases and sequencing adapters. Remaining short reads below 75 bp were filtered out as well as DNA of human origin.

Data analysis

Data preprocessing. The raw data were preprocessed to remove noise and artifacts Metagenomic taxonomical identification was based on clean data after quality control. Kraken2 taxonomic classification software (v2 2.1.2; SILVA database, version 2022-02-02T162959Z) was used for taxonomic profiling [21, 22]. We used the European server of the Galaxy platform for biostatistics (https://usegalaxy.eu/) [23]. Further, Bayesian Reestimation of Abundance was applied to the taxonomy reports produced by Kraken 2 to compute accurately the microbial abundancy level for each metagenomic sample [21]. The bioinformatic workflow for taxonomic identification is available at https://usegalaxy.eu/u/yhodzhev/w/taxonomic-profiling-and-visualization-of-metagenomic-data. The Operational Taxonomic Units (OTUs) table with genus-level microbiome abundance information was generated using the Pavian package [24]. Sequence clustering was performed with a similarity index of 97% [25, 26]. We set up an abundance threshold of ≥ 10 reads to accept a valid OTU. Microbial genera found in negative controls (water, reagents and air) were excluded from further analysis. OTU tables were subjected to statistical evaluation for microbial composition and raw read abundance for the purposes of the study. In this study, the number of reads per OTU/taxon per condition, being the simplest index average count of OTUs, will be referred to as abundance further in the text.

Statistical analysis

Upon inspection of the overall microbial composition (Figure 1(B,C)), OTU counts from the four predominant phyla and the top ten genera were chosen to characterize the potential association of the microbiome with sarcoidosis.This subsetting was greatly influenced by the pronounced rise in data sparsity observed in subsequent taxa. Descriptive statistics (Mean ± Standard deviation) were performed on taxonomic abundance data for the initial characterization of the sarcoidosis and control groups [27]. For illustration purposes dispersion of the mean was represented by the Standard error of the mean (± SEM, Figure 1(B–E)), which is a conservative estimate of the population dispersion [28]. To assess the abundance data normality, the non-parametric Shapiro–Wilk test was applied [29]. The microbiome abundance data, as anticipated, did not conform to a normal distribution, a conclusion supported by the Shapiro–Wilk test results (results are given in Table 1).

Figure 1. Overall characterization of the core blood microbiome abundance (mean ± SEM) at two levels: (A) abundance difference in sarcoidosis and control groups plotted on a decimal logarithmic scale. (B) Top 4 phyla and (C) top 10 genera. (D) Average divergence between the sarcoidosis and control groups according to individual taxon abundance at phylum and (E) genus level. Asterisks indicate significant differences (Kruskal-Wallis group statistics) as follows: ***p < 0.001; **p < 0.01; *p < 0.05; ns – not significant; SEM – standard error of the mean.

Table 1. Confusion matrix for assessing binary classification model.

		Predicted
		Negative	Positive
Actual	Negative	True negative (TN)	False positive (FP)
	Positive	False negative (FN)	True positive (TP)

Assuming this non-normal distribution, the Kruskal–Wallis test was employed to compare the sarcoidosis and control groups [30]. The Kruskal-Wallis test is an appropriate choice in this context because it is a non-parametric method designed to test for differences between two or more groups without the assumption of normally distributed data. Moreover, the Kruskal-Wallis statistics tests if samples originate from one and the same distribution regardless of the shape of this distribution, even if the central tendency (median) differs making it robust against outliers and skewed distributions [30]. For all statistical calculations SciPy library was used under Python version 3.11.4 [31].

A supervised machine learning analysis was used to enrich information about sarcoidosis related microbiome changes obtained by the simple group comparison [14]. For this purpose, we applied a classification decision tree model, where the predicted outcome is the discrete class to which the data belongs (in our case binary, i.e. sarcoidosis vs. control) on the OTU abundance data at the phylum and genus level [32]. Classification tree algorithms are often utilized to analyze the intricate, non-normally distributed data. We opted for this methodology due to its capacity in managing smaller sample sizes and its agnosticism to data distribution [33]. Machine learning was applied in Python version 3.11.4 using the scikit-learn library [34]. For binary classification, the Gini impurity is a measure of the probability of a randomly chosen case being incorrectly classified: $$\mathit{Gini}=1\text{-}\mathrm{ }\sum _{i=1}^j{p\left(\mathrm{i}\right)}^2$$ where j represents the number of classes in the target variable (i.e. sarcoidosis vs. control); p(i) represents the ratio of Correct/False classified number of samples per node [35].

The Gini impurity is 0 when all values at a node belong to the same class, and it is maximized (0.5) when the cases are randomly split between the classes. The most relevant features (taxa) contributing to the binary classification model were identified using feature importance analysis. Feature selection in the context of binary classification is based on the concept of impurity reduction [35]. Each taxon is evaluated by its ability to split the data in a way that makes the resulting consecutive nodes purer than the parent node (g(i) is closer to 0). The models obtained with the lowest nodes’ impurity (feature selection) at the phylum and genus level were reported in the results section. For a given decision tree, the importance I(f) of a feature f is calculated as follows [35]: $$I\left(f\right)={\displaystyle \sum }_{n\in nodesusingf}{w}_n\Delta i\left(n\right)$$ where $w_n$ is the weighted number of samples reaching node n (the weight being the proportion of samples reaching the node to the total number of samples). $\Delta i(n)$ is the reduction in impurity due to the split at node n. It was calculated as follows: $$\Delta i\left(\text{n}\right)=i\left(\text{n}\right)-p_{left}i\left(n_{left}\right)-p_{right}i\left(n_{right}\right)$$ where i(n) is the impurity of the node n; $p_{\mathit{left}}$ is the proportion of samples going to the left resulting of node n; $i\left(n_{\mathit{left}}\right)$ is the impurity of the left resulting of node n; $p_{\mathit{right}}$ is the proportion of samples going to the right resulting of node n; $i\left(n_{\mathit{right}}\right)$ is the impurity of the right resulting of node n.

The generated binary decision tree models for phylum and genus features were validated using the train-test split approach [35]. The entire dataset was divided into two distinct subsets: training and testing. Specifically, 80% of the samples were allocated to the training set, while the remaining 20% (5 controls vs. 1 sarcoidosis patient) constituted the test set. This partitioning ensured that our model was exposed to a broad spectrum of data during training, while still reserving a substantial portion for validation purposes. Using the designated training set, we trained our binary decision tree classifier. Post-training, the predictive capability of the model was rendered using the test set. The performance of the classification model is assessed by a confusion matrix. For the binary classification model, this matrix consists of four entries (Table 1) [35].

Then, the accuracy of the method can be defined using the terms from a confusion matrix as follows: $$\mathit{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$$

Results

Microbiome analysis

The blood samples of 7 patients with sarcoidosis and 21 healthy individuals were collected and sequenced by a targeted 16S rRNA gene amplicon metagenomic approach. A total of nearly 295 000 raw sequences, including human sequences, were obtained after quality checking and preprocessing. After removing the sequencing reads of human origin, the taxonomic analysis assigned over 45 000 sequences to 1541 bacterial OTUs, including 118 OTUs above the abundance threshold (>10 reads). Figure 1(A) shows that the sarcoidosis group grand mean (855 ± 1179) was marginally lower than that in the controls (1773 ± 1760; (H (1) = 2.79; p > 0.06). OTUs were distributed in 6 bacterial phyla (Figure 1(B)) and 41 genera (Figure 1(C)). Proteobacteria, Actinobacteria, Firmicutes and Bacteroidetes formed a group of 4 most prominent phyla (Figure 1(D)). By applying a classical descriptive statistics approach and using the original data set, we characterized the four most abundant phyla.

Genera Bradyrhizobium (5.14 ± 5.01 vs. 947.14 ±1448.25), Pseudomonas (27.86 ± 66.29 vs. 165.67 ±257.74), Comamonas (1 ± 2.2 vs. 131.52 ± 276.63) and Novosphingobium (0.14 ± 0.38 vs. 66.29 ± 92.05) showed significantly lower abundance in sarcoidosis. Although the abundance of the genus Staphylococcus was also numerically lower in the sarcoidosis group (51.57 ± 52.61 vs. 101.71 ± 142.42), this effect was not significant, suggesting a more consistent distribution across the groups for these genera. Table 2 shows the differences in the distribution of the reads for each phylum between the two groups. Notably, Alcanivorax and Micrococcus genera were completely absent in the sarcoidosis group (Figure 1(E); Table 3).

Table 2. Shapiro–wilk test of normality for the distribution of the OTU reads for sarcoidosis and control group at phylum and genus level.

Feature (Taxon)		Sarcoidosis		Control
		Shapiro-Wilk test	n = 7	Shapiro-Wilk test	n = 21
		W	p	W	p
Philum	Proteobacteria	0.84	0.1005	0.75	0.1005
	Bacteroidetes	0.63	0.0006	0.74	0.0006
	Firmicutes	0.71	0.0052	0.82	0.0052
	Actinobacteria	0.69	0.0029	0.90	0.0029
Genus	Bradyrhizobium (Proteobacteria)	0.88	0.1032	0.67	0.00004
	Pseudomonas (Proteobacteria)	0.49	0.0001	0.65	0.00003
	Comamonas (Proteobacteria)	0.54	0.00001	0.53	0.00021
	Staphylococcus (Firmicutes)	0.91	0.4104	0.74	0.0003
	Streptococcus (Firmicutes)	0.66	0.0014	0.63	0.00002
	Veillonella (Firmicutes)	0.66	0.0012	0.35	0.000001
	Novosphingobium (Proteobacteria)	0.45	0.000014	0.72	0.00006
	Prevotella (Bacteroidetes)	0.57	0.0001	0.43	0.000013
	Micrococcus (Actinobacteria)	Invalid input	data	0.79	0.0005
	Alcanivorax (Proteobacteria)	Invalid input	data	0.82	0.0010

Table 3. Descriptive statistics of sarcoidosis vs. control group at phylum level. Mean ± SD along with the Kruskal–Wallis group statistics results.

Phylum	Sarcoidosis (n = 7)	Control (n = 21)	Kruskal–Wallis test H, df = 1	p value
Proteobacteria	174.57 ± 189.20	1514 ± 1646.69	12.077	0.000511
Firmicutes	448.14 ± 667.35	198.90 ± 229.50	0.341	0.087
Bacteroidetes	224.86 ± 421.02	22.14 ± 32.37	2.935	0.56
Actinobacteria	108.71 ± 152.10	94.57 ± 79.11	0.070	0.79

Abbreviations: df – degrees of freedom, H – Kruscal-Wallis test statistics, p – significancnce level.

Figure 1(E) shows that, among the analyzed genera, only Veillonella (232.71 ± 398.24 vs. 2 ± 6.51) and Prevotella (172.57 ± 358.06 vs. 4.38 ± 11.99) exhibited significantly higher abundance in the sarcoidosis group relative to the control group (Table 4). Numerically, Streptococcus abundance (sarcoidosis − 100.43 ± 172.63 vs. control − 53.95 ± 88.45) followed a similar pattern but did not reach significance (Table 4).

Table 4. Descriptive statistics of sarcoidosis vs. control group at genus level abundance (mean ± SD) along with the Kruskal–Wallis group statistics results.

Genus (Phylum)	Sarcoidosis (n = 7)	Control (n = 21)	Kruskal–Wallis test H, df = 1	p value
Bradyrhizobium (Proteobacteria)	5.14 ± 5.01	947.14 ± 1448.25	8.68	0.0032
Pseudomonas (Proteobacteria)	27.86 ± 66.29	165.67 ± 257.74	9.15	0.0025
Comamonas (Proteobacteria)	1 ± 2.2	131.52 ± 276.63	6.2	0.0128
Staphylococcus (Firmicutes)	51.57 ± 52.61	101.71 ± 142.42	0.04	0.8509
Streptococcus (Firmicutes)	100.43 ± 172.63	53.95 ± 88.45	1.41	0.2352
Veillonella (Firmicutes)	232.71 ± 398.24	2 ± 6.51	18.64	0.00002
Novosphingobium (Proteobacteria)	0.14 ± 0.38	66.29 ± 92.05	8.57	0.0034
Prevotella (Bacteroidetes)	172.57 ± 358.06	4.38 ± 11.99	9.16	0.0025
Micrococcus (Actinobacteria)	0	38.29 ± 47.38	9.58	0.002
Alcanivorax (Proteobacteria)	0	34.76 ± 33.4	12.88	0.0003

Abbreviations: df – degrees of freedom, H – Kruscal-Wallis test statistics, p – significancnse level.

By applying the non-parametric Kruskal–Wallis test, Proteobacteria phylum showed significantly higher abundance in the control vs. sarcoidosis group (Table 3, Figure 1(D)). Numerically, the abundance of Firmicutes, Bacteroidetes and Actinobacteria was higher in the sarcoidosis group relative to the control (Figure 1(C)); however, it did not reach significance (Table 2).

Table 4 presents the results of the Kruskal-Wallis test on various bacterial genera to discern differences in their distributions with the significance of their between-group differences.

Decision tree classification

We built a hierarchical binary decision tree based on the abundance data (Figure 2). OTUs data at phylum and genus levels were subjected to the model. Figure 2(A) demonstrates the branching (bifurcation) into tree construction calculated from the phylum abundance data. The root was based on the phylum Proteobacteria, which splits the samples at an abundance threshold of ≤ 168. The root contained a total of 22 samples: 16 controls and 6 sarcoidosis patients, with respective Gini impurity measured at 0.397. The corresponding left node (Node 1) of the root contained 4 samples, all correctly classified in the sarcoidosis group (Gini = 0). Node 1 did not spit further, i.e. it is an end node. The right Node 2 is based on the phylum Firmicutes abundance with a threshold of 519.5. It contained 18 samples and was classified as a control (16 controls vs. 2 sarcoidoses; Gini = 0.198). The respective right Node 3 is an end node with 22 samples and is classified as a control node (16 control vs. 6 sarcoidosis samples; Gini = 0.397). Node 4 branched out Bacteroidetes phylum at abundance threshold of 519.5 (Gini = 0.444). It was classified as a sarcoidosis node (2 sarcoidoses vs. 1 control sample). It further split into Nodes 5 and 6, which were end nodes with Gini impurity values of 0 (Figure 2(A)).

Figure 2. Decision tree diagrams. The control group is marked in orange; the sarcoidosis group is in blue. (A) Phylum; (B) Genus. Boxes represent nodes of classification with bifurcation. All nodes contain information about the classification process: threshold, number of samples, impurity coefficient and proportion of correctly vs. falsely classified samples. Bifurcation or splitting of the classification is shown with arrows. Nodes that did not bifurcate contained samples from only one class. The intensity of the color shade of each node indicates higher confidence in classification.

In the decision tree classification, the feature importance was calculated for each phylum. Proteobacteria stood out as the most influential genus with a feature importance of 0.593, followed by Firmicutes (0.255) and Bacteroidetes (0.153). Actinobacteria, however, did not show any contribution to the model with a feature importance of 0.

Figure 2(B) shows the decision tree structure at the genus level. It starts with a root node, which was split based on the genus Veillonella at a threshold of 0.5; Gini = 0.397. It consisted of 22 samples and was classified as a control node (16 control vs. 6 sarcoidosis). Following this, the left node (Node 1) of the root, represents instances where Veillonella abundance is less than or equal to 0.5. It became a terminal node containing 4 control samples (Gini = 0). On the contrary, the right node (Node 2), represents conditions where Veillonella exceeds 0.5 and bifurcates further on the basis of the genus Comamonas. The Node 3 split is determined by a threshold of 20.0 and is characterized by a Gini impurity of 0.375. It is classified as sarcoidosis with (2 control vs. 6 sarcoidosis samples). Further, this node bifurcates into two end nodes, Node 4 and Node 5, both indicating a Gini impurity index of 0. Node 4 was classified as sarcoidosis with 6 sarcoidosis samples, whereas Node 5 contained 2 control samples exclusively.

At the genus level, the feature importance of Veillonella was 0.656, followed by Comamonas at 0.344. In contrast, the genera Bradyrhizobium, Pseudomonas, Staphylococcus, Streptococcus, Novosphingobium, Prevotella, Micrococcus and Alcanivorax all expressed lack of contribution, each having a feature importance of 0.

During the validation phase, the decision tree model was tested at the phylum and genus levels. The results from this train-test evaluation were identical. All 6 tested subjects (1 of whom was with sarcoidosis) were classified correctly with values of false positive and false negative classifications of 0.

Discussion

With advancements in 16S rRNA gene metagenomic sequencing techniques, the role of microbiome analysis in various diseases became essential. We performed microbiome analysis to assess the microbial richness and abundance of blood samples in sarcoidosis patients and healthy individuals. For this purpose, the average OTUs count was used as a simple measure of microbiome abundance. The taxonomic compositional features of the OTUs provided insight into which phyla and genera were most influential in the classification process, reflecting their potential relevance in discriminating between sarcoidosis and control conditions.

In this study, a machine learning decision tree algorithm was applied to the blood microbiome data of sarcoidosis patients and healthy controls [18] to reveal any microbial abundance patterns at phylum and genus levels potentially associated with sarcoidosis etiology. A total of 295,000 quality-filtered sequences were assigned to 1541 OTUs, revealing Proteobacteria, Actinobacteria, Firmicutes and Bacteroidetes as dominant phyla in both groups. A stark shift in the abundance of Proteobacteria between sarcoidosis patients and controls (174.57 ± 189.20 vs. 1514 ± 1646.69, Table 2) was observed. This marked difference, potentially suggests its role or presence in the systemic manifestations of sarcoidosis. Other phyla, however, did not exhibit statistically significant differences between the two groups. At the genus level, Veillonella and Prevotella abundance was significantly higher in sarcoidosis than in controls (Figure 1(D), Table 3). In contrast, genera that showed significantly reduced abundance in sarcoidosis as compared to the control samples were Bradyrhizobium, Pseudomonas, Comamonas and Novosphingobium (Figure 1(D), Table 3). Two genera, Alcanivorax and Micrococcus, were absent in the sarcoidosis group.

These results were used as a reference to identify the additional information that machine learning analysis could provide. The decision tree classifier successfully distinguished between the sarcoidosis and control groups with high accuracy on the testing data. The decision tree visualization provided insights into the decision-making process, showing how the abundance of specific phyla and genera was used to classify patients and controls. Proteobacteria (0.593) had the highest feature importance, indicating that it plays a significant role in distinguishing the two groups, which is in concordance with the group statistics. The Actinobacteria importance equaled zero, suggesting it might not be influential in the classification. Interestingly, phylum Firmicutes (0.255, Table 3) had an almost doubled coefficient in comparison to phylum Bacteroidetes (0.153), suggesting that changes in its abundance may be a decisive indicator of the blood microbiome dysbiosis associated with sarcoidosis. The genus Veillonella had the highest feature importance (0.656), indicating that this measure governed the discrimination between the sarcoidosis and control groups. The Comamonas genus contributed substantially to the classification as it represents the second highest feature importance by value (0.344). The remaining genera had a feature importance of 0, indicating that they did not contribute to the decision-making process in the trained model.

Implications to the sarcoidosis etiology

The aforementioned alterations in blood microbiota abundance and composition correspond to those observed in other pulmonary diseases [36], supporting the idea that sarcoidosis may be partially driven by microbiome dysbiosis, which has a complex taxonomical origin. A higher relative abundance of Proteobacteria is linked to various pulmonary conditions like COPD (Chronic Obstructive Pulmonary Disease), asthma and cystic fibrosis [37]. Some species may even be related to infections resulting in granulomas, such as Corynebacterium diphtheriae [36, 37]. Firmicutes, mainly Gram-positive, are abundant in human blood [38] and lungs [37] and have been identified in patients with severe COPD conditions [39].

The Veillonella and Comamonas stood out in their contribution to the decision tree classification model, with significant feature importance values (Figure 2(B)). This observation further supported the dominant role of Proteobacteria at the phylum level shown above. This emphasis on Veillonella is particularly intriguing given its higher abundance in sarcoidosis patients. Veillonela species were previously reported in association with various pulmonary conditions such as asthma and COPD [40–42]. Moreover, еnrichment of lower airway microbiota with oral symbiotic bacteria such as Veillonella could be associated with lung granulomatosis and/or carcinogenesis [43]. The overall decline in the abundance of Comamonas spp. may be linked to localized outgrowth of normally suppressed latent microbial forms, such as C. testosteroni and C. acidovorans. These species were occasionally isolated in the respiratory secretions of patients with cystic fibrosis, suggesting their potential role as opportunistic pathogens [44–46]. Such conditions are frequent complication in stage IV sarcoidosis [47]. The clear bifurcations in the decision tree, based on these genera, further emphasized their potential synergetic contribution in sarcoidosis etiology. However, the direct involvement of Comamonas and Veillonella in granuloma formation or their causative role in pulmonary inflammatory diseases needs further experimental confirmation.

Limitations

This study has some limitations. A key constraint of this study was the possibility of the depletion of microbial DNA by the host DNA. Secondly, in our model we included only the first ten most abundant genera. The count of OTUs is an arbitrary measure of microbial abundance, which may lead to inaccurate comparisons of microbial taxa between different studies. The rigorous selection criteria for patients with stage II sarcoidosis resulted in a smaller patient cohort, but leading to a more homogenous study population. The targeted amplicon sequencing was conducted using a sole kit in a single sequencing run, resulting in comparable contamination outcomes.

Conclusions

This study represents the application of a machine learning decision tree algorithm on blood microbiome data in the context of the ethiology of pulmonary sarcoidosis. By upgrading the tradional methods for microbiome analysis with advanced machine learning techniques, we observed additional structural changes in the blood microbiota that may be associated with sarcoidosis. Our findings suggest a multifaceted involvement and interplay of microorganisms. These results may have implications for understanding the etiology of sarcoidosis and may guide research for further improvements in diagnosis and treatment. Validation with a larger dataset including shotgun sequencing and complementary analyses would be essential to verify these findings and explore their biological relevance.

Acknowledgements

The author expresses his gratitude to his colleagues Stefan Panaiotov, Borislava Tsafarova and Dimitar Vasilev for their valuable suggestions and critical comments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the author (YH), upon reasonable request.

Additional information

Funding

This research was funded by the Bulgarian National Science Fund within National Science Program VIHREN, contract number КP-06-DV/10-21.12.2019, and the European Fund for Regional Development through operational program Science and Education for Smart Growth 2014–2020, grant number BG05M2OP001-1.002-0001-C04.