A comprehensive list of human autoimmune disease-related case–control studies on gut microbiota was performed in public databases, including NCBI BioProject (https://www.ncbi.nlm.nih.gov/bioproject) and GMrepo (a curated database of human gut metagenomes; https://gmrepo.humangut.info) (Dai et al., 2022; Wu et al., 2020). A total of 1,954 gut microbiota sequences on 10 different autoimmune diseases (Supplementary Table S1, run-level data including 1,043 cases and 911 controls) were collected. AIDs included RA, SpA, MS, Psoriasis, CD, UC, CeD, MG, SLE, and T1D. The specific inclusion and exclusion criteria are as follows: our inclusion criteria included: (1) Study on 16S rRNA amplicon sequencing of gut microbiota with complete disease phenotype metadata. (2) Case–control studies with at least nine valid samples in each case and control group. (3) No antibiotics or probiotic supplements were administered in the past 3 months. The exclusion criteria were as follows: (1) post-treatment follow-up after medication. We divided these 10 diseases into 7 categories, including 3 digestive system diseases, 2 AIDs of the nervous system diseases, 2 musculoskeletal diseases, 1 endocrine system disease, 1 connective tissue disease, and 1 skin disease, according to the NCBI Medical Subject Headings (MeSH, https://meshb.nlm.nih.gov/) database and Human Disease Ontology (DO) database (Schriml et al., 2022).

We extracted all SRA_ID (listed in Supplementary Table S1) from the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra) (Katz et al., 2022) or European Nucleotide Archive (ENA) (Harrison et al., 2021) and obtained related information on samples from the NCBI BioSample database (https://www.ncbi.nlm.nih.gov/biosample). Raw sequencing data (FATSQ files) were then downloaded using the SRA toolkit, and metadata, including age, gender, and country, were collected. Trimmomatic (Bolger et al., 2014) was used to trim the reads and to remove sequencing vectors and low-quality bases. Reads shorter than 50 base pairs were discarded after trimming. To preprocess the sequences, we used QIIME (2023.5) (Bolyen et al., 2019) to demultiplex and quality-filter the data. Representative sequences and their abundances were extracted using the feature table (McDonald et al., 2012a) to generate tables containing amplicon sequence variants (ASVs). Taxonomic assignment of the individual dataset was classified against the Greengenes database (version 13.8123) (McDonald et al., 2012b). Genus-level relative abundance results were retained for subsequent analyses. Subsequently, samples with only two or fewer taxa were excluded from subsequent analyses. Additionally, to minimize the noise caused by low-abundance taxa, we filtered out those with a relative abundance of < 0.001 across all samples.

All statistical analyses were performed using R, version 4.2.2. The ggpubr package (https://github.com/kassambara/ggpubr) was used to perform non-parametric statistical testing between groups and to account for multiple hypothesis testing corrections when necessary. Significant differences in alpha diversity (Shannon index and richness) were determined using the non-parametric Kruskal–Wallis test. Differences in beta diversity (Bray–Curtis distance matrix calculated using the relative abundances of microbial genera) were determined by PERMANOVA using distance matrices (adonis) in the adonis function of the vegan R package with 999 permutations. Principal coordinate analysis (PCoA) based on beta-diversity was used to visualize the clustering of samples based on their genus-level compositional profiles. To adjust these findings for other factors that may affect microbiota, we used microbiota multivariable associations with linear models (MaAsLin) to identify compositional differences while adjusting for age, sex, and country. To account for the sequencing batch effects of all samples treated at different periods, we used the adjust_batch function implemented in the “MMUPHin” R package using project ID as the controlling factor before model development. Detailed information on the effectiveness of the batch effect removal in this study is shown in Supplementary Figure S1.

Binary sub-cohorts were composed of one AID phenotype and its corresponding healthy control, resulting in a total of 10 binary subgroups. The random forest (RF) model was chosen as the binary classifier because its classification performance has been shown to outperform other methods for microbiota data (Pasolli et al., 2016). The RF model was first trained on a randomly selected training set (5-fold stratified cross-validation) and then applied to the withheld test set to assess the final performance. This process was repeated 20 times to obtain a distribution of RF prediction evaluations for the test set, and the mean AUROC value was calculated.

Multi-class models can be implemented in Python 3.8 using standard libraries that are publicly available, including pandas (2.0.3), numpy (1.22.4), scikit-learn (1.3.1), and matplotlib (3.7.3). For each phenotype, the samples were randomly divided into training (70%, n = 1,368) and test (30%, n = 586) sets. Random forest (RF), support vector machine (SVM), K-nearest neighbors (KNN), multilayer perceptron (MLP), and eXtreme Gradient Boosting (XGBoost) were used as classifier models to diagnose different phenotypes based on the taxonomic profiles at the genus level of the gut microbiota. We employed a 5-fold cross-validation and grid search to select the optimal model parameters for the RF, SVM, KNN, and XGBoost models. For the MLP models, we implemented the default settings provided by Scikit-learn. Finally, we evaluated the performance of the five models on the test dataset as the final performance for predicting different diseases. The highly ranked and frequently selected microbiota features were considered predictive signatures for further interpretation.

The area under the receiver operating characteristic curve (AUROC) is a widely used metric for evaluating the performance of classification models and provides a comprehensive assessment of the model's sensitivity and specificity trade-off at different thresholds. The range of AUROC values is typically explained from 0.5 to 1, with a higher value indicating a better ability to distinguish between different classes of samples. The area under the precision-recall curve (AUPRC), as a complementary assessment, considers the trade-offs between precision (or positive predictive value) and recall (or sensitivity) and is more robust for imbalanced datasets. The AUPRC ranges from 0 to 1, with a value of 0 signifying no positive examples identified and a value of 1 indicating perfect identification of all positive examples. In addition, for a more comprehensive evaluation of the model's performance, we employed the F1 score. It is a metric that considers the precision and recall of a model. The F1 score ranges from 0 to 1, with a higher value indicating better overall performance of the model.

We employed the bootstrap method to address the data imbalance and obtain more robust performance evaluations for each model. The bootstrap method involves iteratively resampling the training data, training a new model, and evaluating the model 1,000 times (Wang et al., 2023; Ning et al., 2025). The performance of the model was calculated as the average performance of the individual models developed using the bootstrap method. Bootstrap methods can considerably reduce overfitting in the developed models. To identify the most discriminative features among many bacterial genus features, minimize model complexity, and enhance computational efficiency and interpretability, we generated a learning curve relating the number of features to the model performance.

Considering the potential influence of the three factors, “gender,” “age,” and “geographical location,” on the gut microbiota, we conducted the sensitivity analysis. For the factor “geographical location,” >75% of the samples in this study were collected from the United States (746) and China (724). We evaluated our model using country-based stratification. For factors “gender” and “age,” the Kruskal–Wallis test was conducted to verify the significance of the differences in the abundance of the corresponding genus in the gut microbiota among different age groups and different genders regarding the AID phenotype. Age groups: Juvenile (3–20), Young Adult (21–45), Older Adult (46–65), and Elderly (65+); Gender groups: Female and Male (Supplementary Table S1, “Sensitivity analysis”).

A total of 1,954 gut microbiota sequencing data (1,043 AIDs cases and 911 non-disease controls; sequences based on Illumina platforms) were collected from the NCBI database based on 19 case–control studies. These data could contain up to 127 cases and 247 controls; however, most studies were conducted on a limited number of samples, with median sizes of 55 cases and 48 controls, respectively. The fecal samples for the gut microbiota sequencing data were mainly collected from five continents and 13 countries, most of which were from the United States of America (38.18%), China (32.44%), and Canada (11.36%) (related information was shown in Supplementary Table S1, “Data availability”).

Ecological indices may not be robust indicators for distinguishing disease from healthy status, which results in changes in the structure of the gut microbiota. First, we aimed to determine the differences in the composition of the gut microbiota among the different AIDs. Compared to healthy controls, significant differences in bacterial diversity (Shannon) and richness (number of genera) were observed in AIDs, except for RA. The Shannon index of the digestive system AIDs decreased. Moreover, we found that both indices (Shannon and richness) varied across phenotypes (Figure 1B). The results of microbial diversity among different phenotypes are consistent with a recent study on multi-class disease diagnosis based on gut microbiota and meta-analysis (Su et al., 2022; Gupta et al., 2020). Beta diversity based on Bray–Curtis dissimilarity showed significant differences in gut microbiota composition among individuals with different phenotypes (R = 0.396; F = 36.057; p < 0.001) (Figure 1C). We then used the linear model of MaALin2 after adjusting for age, sex, country, and technical confounders to explore the associations of microbial composition at the genus level with disease phenotypes. We found 192 significant associations between the 11 phenotypes and 62 bacterial taxa at the genus level (FDR < 0.05). Among the 62 genera, > 67.7% were significantly associated with two or more diseases. The genera Haemophilus, Veillonella, Eggerthella, and Rothia were positively associated with AIDs in our results, whereas the opposite was observed with the genera Paraprevotella, SMB53, and Gemmiger (Figure 1D, Supplementary Table S1 “Gut microbiota data”).

The binary classifier of the RF model based on the microbiota could significantly distinguish between healthy and most AIDs (Figure 2A), which indicated that AIDs had different degrees of intestinal flora disturbance compared with the control. To further investigate the discriminatory ability of the gut microbiota in various AIDs and to distinguish between AIDs and healthy controls, we constructed multi-class classifiers.

To select the best multi-class machine learning algorithm, we evaluated five popular algorithms: RF, SVM, KNN, MLP, and XGBoost models. Five models achieved a mean AUROC of all phenotypes of 0.72–0.89 (Figure 2B), suggesting that muti-class disease classification based on the gut microbiota was feasible. Amongst them, the XGBoost multi-class model had an optimal performance and achieved a mean AUROC of 0.89 [interquartile range, IQR (0.87–0.90)], a mean AUPRC of 0.48 [interquartile range, IQR (0.44–0.51)], and a mean F1-score of 0.538 [interquartile range, IQR (0.51–0.57)] for different disease phenotypes in the test set (Figure 2C). Therefore, the XGBoost multi-class classifier was used for further analysis.

The XGBoost model is constructed using a complete training set. The importance rankings of all features at the genus level were obtained by leveraging this model. Subsequently, the features were incorporated in descending order of importance. For each step of adding a feature, the AUROC value of the model was computed, thereby generating a learning curve that depicted the relationship between the number of features and model performance (Figure 2D). The results indicate that the model performance reached a plateau when 77 features were employed. At this stage, a further increase in the number of features did not lead to a significant improvement in performance. Consequently, the first 77 features were selected as the input variables for the final model. The AUROC values for most phenotypes exceed 0.9 (Figure 3A). The macro-AUROC value was 0.9 [IQR (0.88–0.91)], indicating superior performance compared to the binary classifier we trained. This classifier proved to be valuable for effectively distinguishing AIDs based on features derived from the gut microbiota. At the optimal Youden's index threshold, the sensitivities of our XGBoost multi-class model range from 0.66 to 1 at specificities of 0.70 to 0.95 for different diseases with accuracy from 0.74 to 0.94, highlighting good diagnostic performance (Figure 3B). For example, our model achieved an AUROC of 0.95, CD with a sensitivity of 0.94, and a specificity of 0.89 (Figures 3A, B). To better characterize the XGBoost model, we compared its performance on datasets with different split ratios and obtained similar results, which indicated that the model exhibited high stability and good predictive capability without the risk of overfitting (Figure 3C).

Next, we correlated the top 77 bacterial genera that contributed to the model with the different disease phenotypes. Among the 77 genera, 126 significant associations were found between 42 genera and the different disease phenotypes. These 42 genera belonged to the phyla Firmicutes (28 genera), Actinobacteria (6 genera), Proteobacteria (5 genera), Fusobacteria (1 genus), and Bacteroidetes (3 genera) (Figure 3D). Among the selected bacterial genera, significant associations were found between the phenotypes and several genera, except for T1D (genera only). From the perspective of AID phenotypes, CD (23 genera), RA (21 genera), and psoriasis (20 genera) were the three phenotypes with the highest number of related genera (FDR < 0.05). However, in SpA and CeD, fewer related genera were found (n < 5). From the perspective of the gut microbiota, the genera Actinobacteria and Ruminococcaceae II correlated with the largest number of AID phenotypes (six AID phenotypes). The genera Shigella, Clostridium, and Eggerthella also correlated with many AID phenotypes (five AID phenotypes). These genera could serve as shared microbiota features, which may suggest an association with AID phenotypes except for T1D. The genus Dorea, Lachnobacterium, WAL_1855D, Bulleidia, Pseudomonas, and the special genus Prevotella are only associated with one AID phenotype. This may imply that specific changes in the gut microbiota are related to the corresponding AID phenotype. The genera Clostridium, Eggerthella, Haemophilus, Fusobacterium, Subdoligranulum, and Rothia positively correlated with several AID phenotypes. However, Paraprevotella, SMB53, Clostridium II, Gemmiger, and Slackia were negatively correlated with several AID phenotypes. Only Fusobacterium, which belongs to Fusobacteria, was positively correlated with two inflammatory bowel diseases (CD and UC), indicating a potential microbiota feature of inflammatory bowel diseases.

Interestingly, we noted a higher degree of similarity in microbial alterations between diseases within the same system, such as inflammatory bowel diseases (CD and UC). Among these two phenotypes, 11 genera (approximately 50%) were shared between these two phenotypes, exhibiting similar trends in microbial changes. CeD, also an autoimmune disease of the digestive system, presents a completely distinct profile of related genera. Analogous results have also been reported in Psoriasis and SLE. Although shared microbiota features were relatively scarce, eight genera with similar trends in microbial changes were identified. In autoimmune diseases of the nervous system (MS and MG) and musculoskeletal system (RA and SpA), such similarities in microbial changes were not observed. On the other hand, in some AID phenotypes, alterations in the microbiota showed an opposite trend compared to the healthy group. For example, in Psoriasis and MG, ten microbiota features were shared between the two AID phenotypes, including the genera Actinobacteria, Butyricicoccus, Clostridium IV, Blautia, Lactonifactor, Shigella, Anaerostipes, Clostridium I, Parabacteroides, and Ruminococcus II. However, opposite trends were observed. Similar findings were obtained in a comparison of psoriasis and RA. These results suggest that, in different AIDs, the gut microbiota microenvironment may possess completely opposing characteristics.

The data collected for this study were sourced from multiple countries. More than 70% of the study participants were from the United States and China. The model's performance was evaluated using country-based stratification, revealing consistent results (Figure 4). Moreover, the model attained a mean Area Under the Receiver Operating Characteristic Curve (AUROC) of 0.90 (Interquartile Range, IQR: 0.88–0.93) in the United States and 0.91 (0.88–0.93) in China, respectively. Given the potential influence of “gender” and “age” on the gut microbiota, we performed a stratified analysis based on age and gender (Supplementary Table S1, “Sensitivity analyses”). The sensitivity analysis results indicated that “Geographic location,” “Age,” and “Sex” are not the primary factors affecting classification outcomes.

The strengths of this study include the use of gut microbiota data covering a variety of disease phenotypes (10 AIDs), including gut microbiota sequencing data from almost 2,000 participants. However, this study has some limitations that should be acknowledged. First, we aimed to include gut microbiota sequencing data covering the widest possible range of AIDs and the corresponding healthy controls. However, microbiota sequencing data in public databases often lack relevant information, such as host comorbidities, diet, BMI, and treatment/medication conditions. In our study, only sex (70.5%), age (86.9%), and geographical location (100%) were comprehensively available. Previous research has shown that age and geographical location influence gut microbiota composition, and in some diseases, such as systemic lupus erythematosus (SLE), sex also plays a role (Pang et al., 2023; Bradley and Haran, 2024; He et al., 2018). Therefore, we conducted sensitivity analyses based on these three factors (Figure 4 and Supplementary Table S1, “Sensitivity analyses”). Most genera showed significant differences; however, in the ≥65 age group, some genera did not, likely due to the limited sample size rather than age-related effects. Moreover, because of the limitations of publicly available sequencing data (16S limitations), our analyses were restricted to the genus level, and species-level analyses could not be performed. Second, all datasets were retrospectively obtained, and the classifier was not validated in an external prospective cohort. The complex phenotypes of AIDs, prolonged sample collection periods, and limited availability of prospective data in public databases preclude such validation. Future studies should involve collaboration with relevant clinical teams to prospectively collect fecal samples, utilize gut metagenomic sequencing to obtain deep-level microbiota information, assess the classifier's predictive value for disease progression or prognosis, and enhance its clinical applicability. Finally, the data were cross-sectional, limiting our ability to determine the temporal sequence and causality between abnormal gut microbiota and the onset of AIDs.