General situation
The study involved 64 adult males, comprising 32 cases in the experimental group (age years: 44.44 ± 11.95) and 32 controls in the comparison group with (age years: 40.91 ± 12.28) years. No significant difference in age was observed between the two groups (t= -1.165, P = 0.248). All participants in the case group were farmers (n = 32, 100%), whereas in the control group, occupations comprised farmers (n = 25, 78.1%) and nonagricultural workers (n = 7, 21.9%). Farmers were found to be more susceptible to brucellosis than those in other occupations (X2 = 7.346, P = 0.011), with a higher incidence of the disease among farmers (56.1%). Ethnic distribution analysis revealed predominance of the Hui ethnic group (n = 35, 54.7%), followed by Han (n = 21, 32.8%) and other ethnicities (n = 8, 12.5%). There were no significant differences between the ethnic groups (X2 = 5.697, p = 0.057). Among the 32 male patients with acute brucellosis, occupational exposures were distributed as follows: sheep farming (n = 17, 53.1%), mixed cattle-sheep husbandry (n = 11, 34.4%), cattle ranching (n = 3, 9.4%), and livestock trading (n = 1, 3.1%).The largest subset of patients were sheep farmers, suggesting that working in cattle and sheep farming may enhance the risk of brucellosis infection due to close proximity to these animals. The clinical characteristics of the participants are presented in Table 1. The patients presented with various clinical symptoms, the most common of which included fever, chills, profuse sweating, fatigue, and joint pain.
Serological examination
The RBPT and SAT titers of participants are presented in Table 2. RBPT yielded positive results for the case group. In contrast, the control group exhibited negative results for both the serum RBPT and the SAT. A SAT titer of ≥ 1:100 was considered indicative of a positive result.
Data quality control
Correlation analysis of the QC sample
The metabolome is known for its sensitivity to external interference and its propensity for rapid fluctuations. Consequently, QC is an essential step to statistics by calculating the Pearson correlation coefficients based on the relative quantitative values of metabolites within QC samples. The Pearson correlation coefficients for the quality control samples are plotted in Fig. 1A and B. We found that these coefficients ranged from 0.985 to 1.000. Moreover, the correlation coefficient (R2) is close to 1, indicating an extremely strong correlation among the QC samples, This suggests that the higher the stability throughout the assay process, the superior the quality of the data. The robust data quality demonstrated here provides a solid foundation for subsequent related research.
The Pearson correlation QC samples and PCA score plot of total samples in cationic and anionic models. Both the horizontal and vertical axes correspond to quality control samples, with color intensity gradation reflecting the magnitude of correlation coefficients in this graphical representation, where darker chromatic saturation is proportionally associated with elevated coefficient values. Pearson correlation analysis in cationic mode of QC samples (A). Pearson correlation analysis in anionic mode of QC samples (B). The horizontal coordinate PC1, vertical coordinate PC2 and vertical coordinate PC3 indicate the scores of the first, second and third ranked principal components in the figure, respectively. The scatter dots of different colors indicate the samples of different experimental groupings, and the ellipses are the 95% confidence intervals. Red circle: control sample; green circle: patient sample; blue circle: QC sample. PCA score plot in cationic mode for QC samples (C); PCA score plot in anionic mode for QC samples (D).
PCA analysis of total samples
The peaks derived from the extraction of all experimental samples and QC samples were subjected to PCA. The PCA analysis graph revealed that the QC samples were tightly clustered in Fig. 1C and D, with scores of 7. 80% and 5.40%, respectively. This suggests minimal variability among the QC samples. The distribution of QC samples is clustered together, indicating that the entire method was highly stable, yielding superior data quality and enhanced reliability.
Screening for differential metabolites
PCA
PCA was used to examine the overall distribution trends between the two groups of samples. The PCA scores for the horizontal and vertical axes of the cationic mode are, as below in Fig. 2A, were 7.76% and 6.85%, respectively. PC1 and PC2 represent the first and second principal component scores. The PCA scores for the horizontal and vertical axes of the anionic mode, as shown in Fig. 2B, were 7.98% and 6.76%, respectively. On the score plot, the sample points from different groups are distinctly distributed in separate regions, with a considerable distance between groups. This indicates significant metabolic differences exist among the groups. Within each group, the sample points are more tightly clustered, suggesting a high degree of metabolic similarity within the group.
Multivariate analytical integration of PCA score plots, PLS-DA, and model validation for case-control discrimination in cationic and anionic metabolic models. The horizontal coordinate PC1, vertical coordinate PC2, and vertical coordinate PC3 indicate the scores of the first, second and third ranked principal components in the figure, respectively. The scatter dots of different colors indicate the samples of different experimental groupings, and the ellipses are the 95% confidence intervals. Red circle: control sample; blue circle: patient sample. PCA score plot in cationic mode comparing patient and control samples (A); PCA score plot in anionic mode for the comparison of patient and control samples (B). The scatter plot displays the scores of samples on the first principal component (horizontal axis) and the second principal component (vertical axis). R2Y represents the explanatory power of the model, while Q2Y evaluates the predictive ability of the PLS-DA model. A model is considered well-established when R2Y exceeds Q2Y.PLS-DA model for the cationic mode of the control group with model parameters R2 and Q2 (C); PLS-DA model for the anionic mode of the control group with model parameters R2 and Q2 (D). The sorting validation plot with the horizontal coordinate representing the correlation between Y of the randomized grouping and Y of the original grouping, and the vertical coordinate representing the R2 and Q2 scores. Red dot: control sample; blue dot: patient sample. The sorting validation plot in cationic mode for patient and control samples (E); The sorting validation plot in anionic mode for patient and control samples (F).
PLS-DA score scatterplot and sorting validation plot
The PLS-DA for case-control discrimination in cationic mode is presented in Fig. 2C, the R2Y and Q2Y were 0.94 and 0.71, respectively. The PLS-DA for case-control discrimination in anionic mode is presented in Fig. 2D, the R2Y and Q2Y were 0.93 and 0.64, respectively. The distinct separation of samples between the two groups suggests the potential presence of strong metabolic differences.
The sorting validation plots of the cation and anionic modes are shown in Fig. 2E and F. In both models, the R2 data is greater than the Q2 data (cation: R2 = 0.0,0.80; Q2 = 0.0, -0.40; anion: R2= (0.0,0.78); Q2= (0.0, -0.47) and the intercept between the Q2 regression line and the Y-axis is less than 0. Hence, the models were not overfitted. The models adequately describe the samples and can serve as a prerequisite for the identification of potential biomarker clusters within the model, signifying significant differences in the metabolite profiles of brucellosis between normal individuals and brucellosis patients.
Results of parametric screening for metabolites and differential metabolites
Detection of differential metabolites is presented in Table 3. In the cationic mode, a total of 1,000 metabolites were identified, of which 115 were identified as differential metabolites. Of these, 55 differential metabolites were significantly elevated, and 60 differential metabolites were significantly decreased. In the anionic mode, the analysis detected a total of 556 metabolites, with 74 being screened as differential metabolites. Of these, 38 differential metabolites were significantly elevated, and 36 differential metabolites were significantly decreased.
Metabolite basic information
Status of metabolite classification
The classification results of the metabolites identified in the cationic mode are presented in Fig. 3A, showing nine distinct classes. In the anionic mode, presented in Figs. 3B and 10 classes were detected. Across both modes, the largest group of metabolite species identified was lipids and lipid-like molecules, with organic acids and their derivatives being the next most prevalent.
The classification and KEGG pathway annotation of metabolites in cationic and anionic models. Different colors represent different classification of metabolites, and the size of the color indicates the percentage of each classification of metabolite in the overall metabolites. Classification and proportions of metabolites in the cationic model (A). Classification and proportions of metabolites in the anionic model (B). The horizontal coordinates represent the number of metabolites and the vertical coordinates represent the KEGG pathways annotated to. The figure shows the number of metabolites annotated to each of the secondary classifications under the pathway primary classification. KEGG pathway annotations for the cationic model (C); KEGG Pathway annotations for the anionic model (D).
KEGG pathway annotation classification analysis of metabolite
Metabolic pathways were organized into seven distinct categories within the KEGG pathway database. The identified metabolites were annotated functionally and categorically. The KEGG pathway annotation of metabolites is presented in Fig. 3C and D. The horizontal axis of the figure represents the number of metabolites, while the vertical axis denotes the annotated KEGG pathways. This illustration depicts the count of metabolites annotated within each second-level category under the primary pathway classification. The global and overview maps pathway within metabolism had the highest number of metabolites in both models. Pathway analysis revealed that the metabolic pathway is the most significant pathway that involves the metabolites.
Differential metabolite analysis
Volcanic map of differential metabolites
The volcano map of metabolites in cationic and anionic models were presented in Fig. 4A and B. The volcano plot provides an overview of the distribution of differential metabolites. In the cationic mode, a notably greater number of differential metabolites were downregulated compared to the anionic mode. Additionally, the plot revealed distinct differences in metabolite expression between the case and control groups.
The volcano map, heat map and correlation analysis of metabolites in cationic and anionic models. Horizontal coordinates indicate the fold change of differences in metabolites in different subgroups (log2(fold change)), vertical coordinates indicate the level of significance of differences (-log10(P-value)), and each point in the volcano plot represents a metabolite. Significantly upregulated metabolites are indicated by red dots, and significantly downregulated metabolites by green dots. The size of the dots represents the VIP value. The volcanic maps of differential metabolites in the cationic (A) and anionic models (B). Clustering of samples in the longitudinal direction and clustering of metabolites in the horizontal direction, with shorter clustering branches representing higher similarity. The relationship between the clustering situations of metabolite content between groups can be observed by the horizontal comparison. The heat maps of differential metabolite clustering in the cationic (C) and anionic models (D), respectively. The highest correlation is 1, which is a perfect positive correlation (red), the lowest correlation is -1, which is a perfect negative correlation (blue), and the part with no color indicates a P-value > 0.05. The correlation of the top20 differential metabolites are sorted from smallest to largest P-value value of the cationic (E) and anionic models (F).
Heat map of differential metabolite clustering
Cluster analysis was used to elucidate the metabolic patterns of metabolites across various experimental conditions. Metabolites exhibiting similar metabolic patterns are likely to share functional roles or participate in identical metabolic processes or cellular pathways. The heat map utilizes color intensity within cells to represent the abundance of metabolites in the samples: darker shades indicate higher abundance. The heat map of metabolites in cationic and anionic models were presented in Fig. 4C and D, The cationic mode yielded a greater number of clustered metabolites, whereas fewer metabolites were clustered in the anionic mode.
Analysis of differential metabolites correlation
Metabolites exhibited either synergistic or mutually exclusive relationships with one another. For instance, a positive correlation was observed when a particular type of metabolite consistently changed in the same direction, whereas a negative correlation indicated changes in opposite directions. The objective of the differential metabolite correlation analysis was to assess the consistency of change trends between metabolites. This was achieved by calculating the Pearson correlation coefficient for each pair of metabolites. The correlation analysis of metabolites in cationic and anionic models is presented in Fig. 4E and F, with the P-values sorted from smallest to largest, highlighting the top 20 differential metabolites by correlation. As the linear relationship between two metabolites intensifies, the correlation coefficient approaches 1 for positive correlations and − 1 for negative correlations. Concurrently, a significance statistical test was conducted on the metabolite correlation analysis. In the cationic mode, only one data point exhibited a correlation coefficient close to 1, whereas in the anionic mode, eleven data points showed correlation coefficients near 1.
Identification of differential metabolites following comprehensive database matching
Screening outcomes of serum differential metabolites are presented in Table 4. The differential metabolites initially screened based on the predefined parameters were subsequently validated through comprehensive matching against the mzCloud, mzMasslist, and mzVault databases, leading to the identification of a total of nine differential metabolites. In the cationic mode, five differential metabolites were identified: taurocholic acid (TCA), cis-4-hydroxy-d-proline, inosine, hypoxanthine and L-Iditol. In the anionic mode, the four identified species were taurochenodeoxycholic acid (TCDCA), glycocholic acid (GCA), azelaic acid (AZA), and glycoursodeoxycholic acid (GUDCA). With the exception of cis-4-hydroxy-d-proline, inosine, hypoxanthine, and AZA, all were significantly downregulated.The remaining metabolites were notably upregulated. KEGG pathway enrichment analysis conducted using MetaboAnalyst 6.0 indicated that these differential metabolites were predominantly involved in metabolic pathways including primary bile acid synthesis, purine metabolism, taurine and hypotaurine metabolism, and d-amino acid metabolism.
ROC curve analysis of differential metabolites
The area under the Receiver operating characteristic (ROC) curve, known as the area under the curve (AUC), is utilized to evaluate the sensitivity and specificity of a biomarker in predicting the occurrence of an event. The optimal threshold on the ROC curve determines the sensitivity and specificity for each metabolite. An AUC value of 0.5 indicates that the biomarker is ineffective in predicting the event and possesses no predictive value. As the AUC value exceeds 0.5, its proximity to 1 signifies a higher accuracy in prediction. Generally, an AUC value between 0.5 and 0.7 suggests low prediction accuracy, value between 0.7 and 0.9 indicates a certain prediction accuracy, and values above 0.9 are considered high in prediction accuracy.
In this study, the ROC curve analysis of differential metabolites is presented in Fig. 5, and the results of AUC values, sensitivity and specificity of the differential metabolites are presented in Table 5. The AUC values for the differential metabolites hypoxanthine, l-Iditol, and GUDCA ranged between 0.5 and 0.7 (Fig. 5D and E and Fig. 5I), implying that these biomarkers have low predictive accuracy for the diagnosis of brucellosis. Conversely, the remaining 6 differential metabolites exhibited AUC values between 0.7 and 0.9 (Fig. 5A-C and F-H), suggesting they possess a degree of predictive accuracy for brucellosis diagnosis. The diagnostic performance of the biomarkers was as follows: taurocholic acid demonstrated a sensitivity of 93.75% and specificity of 96.87%; cis-4-hydroxy-d-proline showed 96.87% sensitivity with 90.62% specificity; inosine exhibited 96.87% sensitivity and 81.25% specificity; taurochenodeoxycholic acid achieved both sensitivity and specificity of 96.87%; glycocholic acid presented 96.87% sensitivity with 81.25% specificity; and azelaic acid displayed 96.87% sensitivity accompanied by 93.75% specificity (Table 5).
The ROC curve analysis of differential metabolites. A ROC curve is a curve based on a series of different dichotomous classifications (cut-off values or decision thresholds) with the true positive rate (sensitivity) as the vertical coordinate and the false positive rate (1-specificity) as the horizontal coordinate. The area under the curve (AUC), is used to assess the sensitivity and specificity of the biomarker in predicting the event. The areas under the ROC curve display nine distinct metabolites (A–F).
Box plots of differential metabolites
Box plots of the differential metabolites, derived from comparisons across various groups, provide insights into the changes and variations of these metabolites within different experimental contexts. These box plots are presented in Figs. 6A-I, including TCA (A), cis-4-hydroxy-D-proline (B), inosine (C), hypoxanthine (D), L-Iditol (E), TCDCA (F), GCA (G), AZA (H), and GUDCA (I). The levels of these metabolites exhibited significant variability among the different experimental groups.
The box plots of differential metabolites. Vertical coordinates are quantitative values of metabolites, horizontal coordinates are different subgroups, and the horizontal line in the box (center line) indicates the median of the data. Red boxes: control group; blue boxes: patients group. (A–F) present the box plots for nine different metabolites.
