We used PubMed to systematically search for studies that analyzed 16S rRNA amplicon sequence microbiome data from the vaginal microbiome in individuals with EC and control groups (see Supplementary Section 1 for search keywords). To reduce variability, we included only those studies that collected specimens via vaginal swabs or uterine lavages, used 16S rRNA gene sequencing (any region), and provided publicly available raw sequence data with pathology labels. We excluded studies that used metagenomics, non-amplicon markers, or sampled endometrial tissue. For studies with longitudinal sampling or multiple anatomical sites, only baseline (cervico) vaginal samples were included. Atypical endometrial hyperplasia was grouped with EC, because they are likely to arise together and are clinically treated similarly. All other pathologies, including simple hyperplasia, were deemed benign. Sequence data and metadata were retrieved from the Short Read Archive or materials provided by study authors.

Our goal was to identify a bioinformatics pipeline that consistently reproduced and replicated key microbiome metrics across various studies, including 1) Alpha diversity, measuring diversity within samples, 2) Beta diversity, to measure clustering of microbial communities based on their composition, and 3) differentially abundant taxa associated with disease status (Douglas and Langille, 2021). We define reproducibility as obtaining similar results using the original authors’ analysis pipeline on their dataset and replicability as achieving comparable results when applying that same pipeline to a different dataset (Hejblum et al., 2020). We extracted bioinformatics workflows from published manuscripts.

Each dataset was reanalyzed using both the original pipeline, the pipelines from other published studies, as well as the DADA2 pipeline, a high-resolution denoiser suitable for identifying rare taxa, particularly relevant in Lactobacillus-dominated vaginal communities (Callahan et al., 2016; Nearing et al., 2018). Samples were excluded based on pipeline-specific quality thresholds or incomplete metadata (Supplementary Section 2). Processing steps for the DADA2 pipeline included, adapter trimming using BBDuk¹, quality filtering using Phred scores and trimming the 3’ ends of reads where the average quality dropped below 20 (see Supplementary Section 3 for read lengths maintained). Forward and reverse reads were denoised and paired-end reads were merged. In studies where there was minimum to no read overlap (<50% samples merging with >12bp overlap; details in Supplementary Section 3), only forward reads were used as previously done by (Sekaran et al., 2023; Abdill et al., 2025). After filtering out bimeras, the Amplicon Sequence Variants (ASVs) abundance table was normalized by the total number of reads sequenced in each sample. Identified ASVs were assigned species-level taxonomic information using the naïve Bayes classifier implemented in QIIME2 based on 1) the Green Genes (GG) database (v13.8), 2) the Genome Taxonomic Database (GTDB) (vbac120) both with uniform taxonomic distribution, and 3) the GTDB database with an expected species distribution using q2-clawback (DeSantis et al., 2006; Fettweis et al., 2012; Bokulich et al., 2018, Bokulich et al., 2022; Bolyen et al., 2019; Kaehler et al., 2019; Rinke et al., 2021).

To merge taxonomies from GG and GTDB, we used superstring matching and RESCRIPt to generate a consensus taxonomy based on the last common ancestor (Robeson et al., 2021; Bokulich et al., 2022). A phylogenetic tree was built by aligning ASVs with MAFFT, processed with FastTree, and midpoint-rooted using the phangorn package (v2.11.1) in R (Price et al., 2010; Katoh and Standley, 2013). Supplementary Section 4 outlines the bioinformatics pipelines implemented.

For each dataset, we compared the performance of the previously implemented pipelines with DADA2. We used the Shannon index to measure alpha diversity. For beta diversity, we tested the marginal significance of available participant characteristics and disease status using the PERMANOVA test (did not implement multiple testing correction), along with distance measures reported in individual studies (such as Bray-Curtis, calculated using the relative abundances of taxonomic features, and (weighted or unweighted) UniFrac, which considers the phylogeny of the taxonomic features as well) (Bray and Curtis, 1957; Lozupone et al., 2011; Anderson, 2017). For the DADA2 pipeline, we used the Bray-Curtis, Jaccard (Real and Vargas, 1996), UniFrac, and Jensen-Shannon Divergence (Fuglede and Topsoe, 2004) metrics. The weighted UniFrac metric was calculated using the phylogenetic tree generated from ASVs. Both alpha and beta diversity were measured at the ASV level. To identify differentially abundant taxa in EC participants, we aggregated the abundance tables to the species level (or the genus taxonomic level, if unassigned at the species level) and used ANCOM-BC with multiple testing adjustment using the Holm method (Lin and Peddada, 2020). All microbiome metrics evaluated were adjusted for potential confounding patient factors, including age, BMI, and ethnicity. The pipeline that consistently demonstrated trends reported in literature for alpha/beta diversity and associated EC taxa was selected for downstream predictive modeling.

We compared several data integration strategies to develop a predictive vaginal microbiome signature for EC (Figure 1). One study, Antonio et al (Walther-António et al., 2016). was set aside for validation. This was selected because it was not too large nor too small and had balanced number of EC and benign diagnoses. The remaining studies were used for model training. Our baseline approach used early integration, where datasets were concatenated into a single high-dimensional matrix (Picard et al., 2021). In the batch-corrected approach, we used ComBat Johnson et al., 2007) to adjust for batch effects and assessed clustering by batch versus disease status using the PERMANOVA test (Anderson, 2017).

We also evaluate a late integration strategy, where separate classifiers were trained on individual datasets and combined to generate ensemble predictions (Mienye and Sun, 2022). This allowed the inclusion of patient-level characteristics when available and permitted the use of different algorithms for different datasets (Mienye and Sun, 2022).

Under late integration, we trained four models; (1) a microbiome-only model, (2) a model using patient characteristics (age, BMI, and ethnicity) only, (3) a model combining patient characteristics and microbiome data, (4) a model using patient characteristics and vaginal pH as a microbiome biomarker, since vaginal pH is a consequence of the vaginal microbiome composition and associated with EC risk (Walther-António et al., 2016).

All models were built with CLR transformed, genus-level count data and adhering to the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis + Artificial Intelligence guidelines (Collins et al., 2024) (checklist in Supplementary Section 5). Taxa present in fewer than 5% of samples were filtered out. For each integration strategy, we evaluated the following algorithms: random forest algorithm (randomForest v4.7-1.1), gradient tree boosting algorithm (xgboost v1.7.7.1), and neural network algorithm (nnet v7.3-19) using the tidymodels (v1.2.0) package. We also evaluated ensembles comprising three base models. Hyperparameter tuning was performed using grid search (details in Supplementary Section 6) with a five-fold stratified cross-validation, repeated three times, and optimized for F1 score at a 0.5 classification threshold. Final ensemble models were selected based on out-of-fold performance and tested on a held-out test. To demonstrate the generalizability of the model that performs the best on the held-out dataset, we implemented a leave-one-study-out (LOSO) validation approach, retraining the models each time to only include the covariates available in the held-out dataset.

We reported sensitivity, specificity, negative predictive value (NPV), positive predictive value (PPV), and area under the receiver operating characteristics (AUROC) using the yardstick (v1.3.1) package. For the LOSO validation, we report pooled metrics. Exact 95% confidence intervals were calculated using epiR (v2.0.80). To correct for class imbalance, we used SMOTENC from the themis package (v1.0.3). For datasets with partially missing data, Multivariate Imputation by Chained Equations (mice v3.16.0) was used to impute missing values, conditioned on other participant data and microbiome profile. All analyses were performed in R (v4.3.3) and RStudio (v2023.06.2). Executable code and pipeline parameters are available on Github².

The systematic literature review yielded 248 articles, of which 11 used 16S rRNA amplicon sequencing to assess vaginal microbiome profiles. We excluded 5 studies that relied on invasive tissue sampling and 1 study that did not provide access to sequencing data, resulting in 5 eligible datasets for analysis.

The first study, published by the Mayo Clinic in 2016, included 22 participants and investigated the vaginal microbiome composition and its putative role in EC (Walther-António et al., 2016). The study reported significant age differences between EC patients and controls with benign conditions (mean age 62 vs 47 years). A follow-up study by the same group in 2019 expanded the cohort to 149 participants and evaluated participant characteristics, including menopausal status, BMI, vaginal pH, and age, and their associations with microbial composition (Walsh et al., 2019). This second cohort demonstrated differences in EC patients compared to those with benign diagnoses. Notably, EC patients were older, had higher BMI and vaginal pH (Table 1).

Tsementzi et al., 2020 analyzed vaginal microbiome composition in EC and cervical cancer (n=38) compared to participants with benign conditions. EC patients had significantly higher vaginal pH (> 4.5) compared to controls. Gressel et al., 2021, examined postmenopausal women (n=27), profiling cervicovaginal microbiota to identify taxa associated with EC. Most recently, Chao et al., 2022 used lavage-based vaginal sampling to compare microbial taxa associated with EC and benign diagnoses (n=31). Patient characteristics (age, BMI, ethnicity, and vaginal pH) were available for all cohorts except Chao et al. (age only) and Gressel et al. (no patient metadata). The available demographic and clinical characteristics of the participants are summarized in Table 1.

Biological specimens across all studies were collected by physicians; all cohorts were recruited in the United States of America except for Chao et al., which enrolled participants in China. Genomic DNA extraction protocol varied: the Mayo Clinic studies used the MoBio PowerSoil Kit, Tsementzi et al. used the DNeasy PowerSoil Kit, and Gressel and Chao et al. used the QIAamp DNA Mini Kit. All studies sequenced the 16S rRNA gene on the Illumina MiSeq platform, targeting various hypervariable regions, with V4 being the common region across all datasets. Study-specific details, including inclusion/exclusion criteria, sampling and storage protocols, DNA extraction, primer design, and controls, are provided in Supplementary Section 7.

Bioinformatic pipelines also varied across studies: Antonio et al. and Walsh et al. employed IM-TORNADO (Jeraldo et al., 2014), an in-house pipeline that concatenates paired end-reads and processes them using Mothur and USEARCH to generate Operational Taxonomic Units (OTUs) (Schloss et al., 2009; Edgar, 2010). Chao et al. used UNOISE to denoise sequences and produce ASVs (Edgar and Flyvbjerg, 2015). Tsementzi et al. Gressel et al. used QIIME2 with VSEARCH and USEARCH, respectively, to generate OTUs (Edgar, 2010; Rognes et al., 2016; Bolyen et al., 2019). We reprocessed all the datasets using those four pipelines and added a pipeline based on DADA2 (a total of 25 data/pipeline combinations).

All included studies originally reported increased vaginal microbial (alpha) diversity in EC patients compared to controls. We were able to reproduce this trend across datasets; however, results varied depending on the bioinformatics pipeline used (Figure 2A). The expected trend was not observed in the Tsementzi cohort when processed using the IM-TORNADO pipeline, nor in the Walsh cohort processed with the Tsementzi pipeline. Similarly, alpha diversity was not recapitulated in the Tsementzi and Chao datasets using the Gressel pipeline or in the Antonio dataset using the Chao pipeline. Among all datasets, the Gressel cohort demonstrated the most consistent trends across pipelines and the best quality profile (see Supplementary Section 3). However, their pipeline, which does not filter for chimeric reads, systematically produced inflated alpha diversity values in other datasets.

In contrast, processing with DADA2 consistently replicated the expected alpha diversity trend across all datasets (with no significant differences observed), suggesting greater robustness to technical variability.

Beta diversity appeared to be less sensitive to preprocessing differences (Figure 2B). Biological trends remained largely consistent across pipelines and metrics. We note that disease status accounted for less than 16% of the variance in all cohorts, regardless of the pipeline or distance metric (see Supplementary Section 8). By contrast, when available, participant-level characteristics, such as ethnicity, age, and BMI, explained a greater proportion of variance.

Using ANCOM-BC, we adjusted for available covariates such as age, BMI, and ethnicity to identify differential abundant taxonomies associated with disease status (Figure 2C). Across all datasets, various species of Peptoniphilus were consistently associated with disease status. Specifically, we observed enrichment of Peptoniphilus urinimassilliensis in the Walsh dataset, Peptoniphilus coxii in the Antonio and Chao datasets, and Peptoniphilus sp000478985 and sp900099555 in the Tsementzi and Gressel datasets, respectively. Notably, Peptoniphilus coxii and sp900099555 were also associated with BMI and age, respectively, in the Antonio and Walsh dataset. In addition to Peptoniphilus, species from the genera Prevotella, Streptococcus, and Blautia were differentially abundant across four of the five datasets, with the exception of Walsh et al., where these taxa were instead associated with BMI or ethnicity, both known EC risk factors.

Differentially abundant taxa reported in the original studies could not be reproduced using the published bioinformatics pipelines but were consistently recovered using the DADA2 pipeline. The DADA2 pipeline consistently outperformed other methods in identifying consistent ecological trends and was therefore selected as the primary preprocessing pipeline for model development.

Four studies (Chao, Gressel, Tsementzi, and Walsh) were used for model training and one for performance evaluation (Antonio). Various data integration and modeling strategies were assessed, as outlined in Figure 1. Performance metrics for all ensemble models are detailed in Table 2 and Supplementary Section 9.

Early integration models, which pooled data prior to modeling, showed limited ability to correctly identify benign cases. Batch correction using ComBat removed study-specific microbiome clustering (adonis PERMANOVA) while preserving variance due to disease status (see Supplementary Section 10). However, models trained on batch-corrected or uncorrected data performed similarly in identifying EC cases, both achieving perfect sensitivity, but very poor specificity (0; 95% CI; [0, 0.31]), which resulted in undefined NPV.

Late integration models using only microbiome data achieved moderate performance in a held-out test set (sensitivity: 0.67 (95% CI; 0.35-0.90) and NPV: 0.64 (95% CI; 0.31-0.89)). Predicting with participant metadata where available: age, BMI, and ethnicity for Tsementzi and Walsh; age only for Chao—improved specificity 0.9 (95% CI; 0.56-1) but reduced sensitivity to 0.67 (95% CI; 0.35-0.90) and NPV to 0.69 (95% CI; 0.39-0.91). Including vaginal pH, a downstream biomarker of microbiome shifts, with other patient characteristics, did not enhance predictive value (Table 2).

The highest-performing model was an ensemble approach that integrated both microbial and host characteristics where available. Applied to the held-out test set, this model achieved perfect sensitivity of 1.0 (95% CI; 0.74-1) and NPV of 1.0 (95% CI; 0.59-1), with a specificity of 0.7 (95% CI; 0.35-0.93) and AUROC of 0.93 (95% CI; 0.71-0.93). Feature importance analysis across ensemble frameworks (Supplementary Section 11) identified Lactobacillus, Prevotella,

Peptoniphilus, Porphyromonas, Peptostreptococcus, and Streptococcus among the top 10 predictors, consistent with ANCOM-BC findings.

The LOSO validation of the ensemble model resulted in a pooled AUROC of 0.7 (95% CI; 0.6-0.7), ranging from 0.7 to 1.0 across individual studies, except for Tsementzi et al., which had the lowest data quality and an AUROC of 0.5 (Supplementary Section 12). The validation on the Chao et al. dataset, achieving perfect discrimination (AUROC of 1.0; 95% CI: [1-1]) and an NPV of 1.0 (95% CI: [0.9-1]), while the discrimation on the Walsh et al. and Gressel et al. datasets was moderate, with AUROC values of 0.8 (95% CI: [0.7-0.8]) and 0.7 (95% CI: 0.4-1), and NPVs of 0.8 (95% CI: 0.7-0.9) and 0.6 (95% CI: 0.4-0.8), respectively.