The human microbiome, comprising diverse microorganisms across various body sites, is crucial for understanding health and disease dynamics (Ursell et al., 2012; Wang and LêCao, 2020). Advances in sequencing technologies have deepened our understanding by enabling researchers to not only identify microbial species but also to understand their functional roles within the ecosystem. However, these advancements have also highlighted challenges like batch effects—discrepancies across sample batches that can distort biological insights (Leek and Storey, 2007). That is, batch effects refer to systematic and non-systematic discrepancies or variations that arise during the processing of samples in a study (Lazar et al., 2013; Wang and LêCao, 2020). Systematic batch effects represent the consistent differences across all samples within a batch, and nonsystematic batch effects demonstrate variability that is dependent on the diversity of OTUs present within each individual sample of the same batch (Soneson et al., 2014). Key factors contributing to these effects include variations in sample collection, DNA extraction methods, sequencing protocols, and data analysis techniques. Microbiome data's inherent properties—high zero-inflation, over-dispersion—exacerbate the impact of batch effects. Zero-inflation indicates the presence of many zeros in the data, which highlights the scarcity or absence of certain microbial species in many samples. Over-dispersion arises from individual variability and technical differences in sequencing depth, complicating data analysis.

While numerous methods exist for correcting batch effects in high-throughput data, many assume continuous data and are not directly applicable to microbiome data, which typically consist of count data. Traditional methods, such as the ComBat algorithm (Johnson et al., 2007), typically assume a Gaussian distribution. Based on this, extensions for count data such as RNA-seq have been developed that utilize a negative binomial distribution. This approach is more aligned with the inherent properties of the data and effectively adjusts for consistent batch patterns. However, while this method is applicable for addressing systematic batch effects, it struggles to fully correct for non-systematic batch influences caused by irregular experimental errors or the unique characteristics of individual samples.

Another approach, Meta-analysis Methods with a Uniform Pipeline for Heterogeneity in microbiome studies (MMUPHin) provides a comprehensive solution for managing heterogeneity in microbiome studies (Ma et al., 2022). This method includes joint normalization and meta-analysis techniques specifically designed for the unique characteristics of microbiome data. Its adaptability is especially effective for handling the diverse and often non-parametric nature of microbiome count data. However, MMUPHin assumes the data to be Zero-inflated Gaussian, which is primarily suitable only for certain transformations of relative abundance data, such as taxon counts normalized by each sample's library size. This assumption limits its applicability, indicating a need for more flexible approaches in certain scenarios. Percentile normalization is a method by which data can be converted to a uniform distribution based on percentiles (Gibbons et al., 2018). This transformation helps mitigate the effects of over-dispersion and the high zero count. However, a key limitation of percentile-normalization lies in its potential to oversimplify complex data structures, possibly leading to the loss of meaningful biological variance. This issue can be particularly critical when dealing with highly diverse microbiome samples. Recently, new methods for batch effect correction have been introduced, tailored specifically to design unique characteristics of microbiome data. Among those approaches, Conditional Quantile Regression (ConQuR) (Ling et al., 2021a,b) utilizes a conditional quantile regression approach for batch effect correction, enabling the generation of corrected OTUs. ConQuR independently processes each OTU without assuming a specific distribution, offering flexible handling of data. However, when OTUs across batches exhibit significantly different distributions, each quantile may represent different characteristics, potentially compromising the consistency of the analysis results. While using a reference batch to standardize OTU distributions can address fundamental inter-batch differences, the effectiveness of this approach heavily depends on the extent of these differences and the appropriateness of the chosen reference batch. If the reference batch does not adequately represent the characteristics of other batches, this method may introduce additional distortions.

In our study, we conducted an empirical evaluation of batch effect correction methods, including MMUPHin, Percentile Normalization, ConQuR, ComBat and our proposed model. To demonstrate the effectiveness of our proposed model, we used PERMANOVA R² (Tables 3, 4), Average Silhouette Coefficient for quantification and Principal Coordinates Analysis (PCoA) plots (Figures 1–5) for visualization. A range of dissimilarity metrics commonly used in microbiome research were employed for this evaluation. We compared the values for each method across different dissimilarity metrics (Song et al., 2023): Bray-Curtis, Aitchison, Canberra, and Manhattan. The outcomes of the Kruskal-Wallis test indicate no significant differences between the distributions of batch groups, suggesting that the negative binomial regression model appropriately formed similar distribution pattern by excluding consistent effects across each batch.

Results from the analysis of the Human Immunodeficiency Virus Re-analysis Consortium dataset are presented below. These were analyzed to evaluate the effectiveness of our proposed model in correcting batch effects and maintaining the integrity of the microbiome data. For Bray-Curtis dissimilarity, the proposed model achieved the lowest R² value of 0.0128^*, indicating superior performance in batch effect correction compared to other methods. ConQuR also performed well with a R² value of 0.0149^*, while MMUPHin and ComBat showed moderate performance with R² values of 0.0822^* and 0.0637^*, respectively (Table 4). Regarding Aitchison dissimilarity, the proposed model had a higher R² value of 0.0854^*, indicating less effective correction in this specific metric and ComBat demonstrated better performance with R² values of 0.0564^*. ConQuR showed the highest R² value of 0.0901^*, suggesting the least effective correction in this metric. For Canberra dissimilarity, MMUPHin achieved the lowest R² value of 0.0779^*, closely followed by ComBat with an R² value of 0.0825^*. The proposed model had R² value of 0.1493^*, which was lower than ConQuR's R² value of 0.1628^*, indicating moderate performance. In the case of Manhattan dissimilarity, the proposed model excelled with the lowest R² value of 0.0065^*, demonstrating excellent batch effect correction. ConQuR and ComBat showed moderate performance with R² values of 0.0087^* and 0.0137^*, respectively. This analysis confirms the effectiveness of various methods in addressing batch effects, with significant values (P ≤ 0.05) marked by an asterisk (^*). This calculation was conducted using the “adonis” function of vegan package (2.6–6.1) in R.

Overall, our proposed method consistently achieved lower R² values, indicating a significant reduction in batch effects across various dissimilarity measures. This performance underscores the robustness and superior effectiveness of our method in correcting batch effects in microbiome data. MMUPHin and Percentile Normalization exhibited intermediate performance across most dissimilarity measures, while ConQuR and ComBat showed varying effectiveness depending on the metric used. The effectiveness of batch effect correction varies across different distance metrics, as each metric considers different aspects of the data. Bray-Curtis and Manhattan metrics reflect compositional differences between samples. The Bray-Curtis metric focuses on the relative abundance of shared OTUs between samples, highlighting how similarly the samples are composed in terms of the OTUs they have in common. Similarly, the Manhattan metric considers the absolute differences in abundance, treating all changes equally. Therefore, Bray-Curtis and Manhattan metrics effectively capture overall abundance differences and compositional changes in data. On the other hand, Aitchison and Canberra metrics emphasize ratio-based differences and small values. The Aitchison metric, used for compositional data, is sensitive to variations in the relative abundance of OTUs, considering ratio fluctuations. The Canberra metric assigns greater weight to changes in small values, highlighting differences in proportions with smaller magnitudes. Due to these characteristics, in composite quantile regression models, the effectiveness of adjustments can be diminished when using distance metrics like Aitchison and Canberra, depending on the selected regression coefficients (Ricotta, 2017). This comprehensive evaluation highlights the proposed model's ability to handle batch effects more effectively, ensuring reliable analysis of microbiome data.

The proposed model consistently achieved the lowest R² values across all distance metrics in the HPV dataset, demonstrating optimal batch effect correction performance with values of 0.002^*, 0.0246^*, 0.0264^*, and 0.0029^* for the Bray-Curtis, Aitchison, Canberra, and Manhattan distance metrics, respectively (Table 4). The prominent results of the proposed model with the Bray-Curtis and Manhattan distance metrics can be attributed to the characteristics of the composite quantile regression model described previously.

Our model employs a Negative Binomial approach to model the variability within batch effects, aiming to align batch distributions and exclude generalized effects across batches. We select the reference batch through a methodical process, employing the Kruskal-Wallis test and median absolute deviation to identify batches that consistently reflect the general characteristics of the data. Furthermore, our approach incorporates composite quantile regression which addresses the effects within each batch at the OTU level. This procedure ensures that our model not only corrects for obvious systematic differences but also subtly adjusts for less predictable changes within the microbiome data, leading to a more reliable results of biological data across different batches.

When applied to the HIVRC data, and the HPV study data, this model demonstrated considerable success in reducing batch effects. It showed improved correction results compared to the application of conditional quantile regression independently, consistently across these different datasets. By effectively addressing both systematic and non-systematic batch effects, the proposed method ensures a more reliable and accurate analysis of microbiome data. We also verified that similar results were obtained using an additional dataset (Ewha Women's University Medical Center metagenomic urine control data−16S rRNA sequencing data). Given the similarity of observed patterns, we chose to present results from two representative datasets.

Despite the effectiveness of our proposed method, several limitations should be acknowledged. First, microbiome data often contain a high proportion of zeros that can lead to instability in quantile regression at extreme quantiles (e.g., lower 10% and upper 90%) due to the limited number of non-zero data points. This sparsity can affect the robustness and accuracy of the regression estimates at these extreme quantiles. In practice, the dataset with the lowest zero inflation (31%) showed the best results in correcting batch effect based on quantile estimation and based on this, it was determined that the rate of zero inflation could act as a factor that decreases the stability of the model. Furthermore, the performance of the model may vary depending on the metrics or categories of data used. Different datasets or dissimilarity distance measures may influence the effectiveness of the batch effect correction, indicating that the model's performance is not universally consistent across all possible applications. Lastly, the efficacy of our model may be compromised in research designs where the separation of confounding variables from primary variables is challenging. Addressing batch effects in confounded study designs presents considerable difficulties, and careful experimental design can serve as a critical factor in alleviating batch effects. These limitations highlight the need for further refinement of the model to enhance its robustness and applicability across diverse microbiome datasets. Additionally, the model assumes that the count data follow a negative binomial distribution, which may not always be the case in real-world microbiome data. Also, when various trends or patterns clearly emerge across different quantiles of the data, it may be more appropriate to use regression coefficients specialized for each quantile.