The exponential growth of High-Throughput Sequencing (HTS) data on the microbial communities presents researchers with an unparalleled opportunity to delve deeper into the association of microorganisms with host phenotype. However, this growth also poses a challenge, as microbial data is complex, sparse, discrete, and prone to zero-inflation. Moreover, current methods for integrating microbiome data and other covariates are severely lacking. Hence, we propose a Bayesian zero-inflated negative binomial (ZINB) regression model that is capable of identifying differentially abundant taxa with distinct phenotypes and quantifying the effects of covariates on these taxa. Our model exhibits excellent performance when tested on simulated data. Upon successfully applying our model to a real multi-ethnic cohort study, we discovered that the prevailing understanding of microbial count data from previous research was overly dogmatic, because only a subset of taxa demonstrated zero inflation in real data. Moreover, we have discovered that dispersion parameters significantly influence the accuracy of model results, and increasing sample size can alleviate this issue. In all, we have presented an innovative integrative Bayesian regression model and a comprehensive pipeline for conducting a multi-ethnic cohort study of children, which facilitates bacterial differential abundance analysis and quantification of microbiome-covariate effects. This approach can be applied to general microbiome studies.