Machine learning-assisted genomic profiling to identify differences between Bacillus Calmette-Guérin (BCG) vaccine strains and non-BCG wild-type Mycobacterium bovis

Strain selection and genome data

In this study, we systematically searched the NCBI database using the keyword “Mycobacterium tuberculosis variant bovis” and identified 76 high-quality genome assemblies at the scaffold level or higher. Among these, 28 genomes were explicitly designated as BCG derivatives in their nomenclature and were included in subsequent analyses. Additionally, four genomes lacking BCG designation in their names but annotated with BCG-related information in the “Modifier” field were classified as taxonomically ambiguous samples. Furthermore, 44 confirmed non-BCG wild-type mutant strains were identified for comparative analysis (Fig. 1). In the end, the accession numbers and information for the 72 clearly taxonomically classified genomes used can be found in Excel File S1).

Computational infrastructure

All computational analyses were executed on a high-performance computing cluster utilizing two-socket Intel Xeon Gold 6230 processors (2.1 GHz base/3.9 GHz Turbo) with 40 physical cores (80 hyperthreaded logical processors) per node, Non-Uniform Memory Access (NUMA)-optimized memory allocation across two domains, and 1 TB of 2933 MHz DDR4 ECC RAM in 6-channel configuration. Workflow orchestration was implemented through Snakemake v7.32 with Docker 24.0 containers (CentOS 8.3 base), where the TensorFlow 2.15.0 framework exploited CPU-level parallelism via Intel MKL-DNN acceleration.

Data preprocessing

The preprocessing pipeline focused on high-quality genome assemblies (≤200 contigs/genome) directly retrieved from NCBI, beginning with assembly integrity validation through (1) cross-referencing NCBI metadata for completeness flags and contiguity metrics (N50 ≥50 kb), (2) taxonomic confirmation via genome-wide average nucleotide identity (ANI ≥98% against M. bovis reference genomes). Raw read processing was omitted as analysis utilized pre-assembled genomes bearing NCBI’s International Nucleotide Sequence Database Collaboration (INSDC) quality certification (contamination scores <1%, assembly gaps <5%).

Bitscore values calculation of core orthologous genes

We annotated the 72 selected genome sequences using Prokka (Seemann, 2014), generating structural annotation files in General Feature Format (GFF) format and protein sequence files in FAA format. Subsequently, a pan-genome comparative analysis was performed with Roary (Page et al., 2015) (default parameters) based on the 72 GFF files, yielding a corresponding gene presence/absence matrix (gene_presence_absence.csv).

For further investigation, the DeltaBS approach was employed (Gardner BinfLab, 2021). This method utilizes profile hidden Markov models (HMMs) to integrate sequence diversity of homologous genes across organisms and capture natural variation patterns. Specifically, the profile HMMs of Gammaproteobacterial proteins (gproNOG.hmm) were retrieved from the eggNOG database (http://eggnog5.embl.de/) as a reference set of conserved bacterial orthologous groups (OGs). Although taxonomically labeled for Gammaproteobacteria, these OGs represent a subset of bacterial OGs (bacNOG) conserved across phyla including Actinobacteria. Annotated protein sequences from Prokka were aligned to their corresponding HMM profiles using the hmmsearch tool in the HMMER3.0 package (http://hmmer.org), generating bitscore values for each sequence. The DeltaBS method does not require phylum-matched HMMs as it uses these profiles as universal references to measure gene conservation. Core orthologous genes of each strain were then extracted using the parse_hmmsearch.pl and parse_bitscores.pl scripts from DeltaBS (Gardner BinfLab, 2021) to construct the DBS metric table (bitscores.tsv). This table served as the input for training and constructing a random forest classifier.

Random forest classifier constructing and training

Given that the 72 genomes collected in this study represent unambiguously classified BCG vaccine strains and nonBCG wild-type Mycobacterium bovis, we use all the 28 complete BCG representative genomes and 44 non-BCG wild-type mutant strains genomes to develop a random forest classifier model for distinguishing between the two groups. Training the model on a set of 3,595 orthologous genes enabled effective differentiation between BCG vaccine strains and non-BCG wild-type Mycobacterium bovis. Model performance was evaluated using out-of-bag (OOB) accuracy as the primary metric. The random forest classifier was implemented using the randomForest and caret packages in R. To optimize model parameters, we systematically evaluated tree configurations. The number of trees (ntree) was set to 10,000, as the out-of-bag (OOB) error rate stabilized beyond this threshold. For mtry (the number of genes randomly sampled as candidate predictors at each node), a grid search was performed across values {1, n/10, n/5, n/3, n/2, n}, where n represents the total number of predictor genes (3,595). mtry = n/10 was selected to minimize collinearity among predictors, yielding the lowest OOB error (0.2).

Model refinement was achieved through iterative feature selection:

• (1)

  Initial model (Model 1): All 3,595 genes were included.

• (2)

  Sparsity pruning: Predictors with variable importance (VI) ≤ 0 were iteratively removed.

• (3)

  Iterative optimization: The model was retrained using the reduced gene set, followed by exclusion of the lowest 50% VI-ranked genes.

This process was repeated until perfect OOB accuracy was attained.

To assess the null hypothesis (no association between predictors and strain type), 1,000 permuted datasets were generated by randomizing labels. The empirical p-value was calculated as the proportion of permuted models achieving accuracy equal to or greater than the original model. Top discriminative genes from the final model were functionally categorized via COG and KEGG databases using COG (https://github.com/transcript/COG) (Westreich, 2019) and KAAS (KEGG Automatic Annotation Server https://www.genome.jp/tools/kaas/) separately (Moriya et al., 2007).

Convolutional neural network based genomic classifier construction based on top predicted genes

The convolutional neural network (CNN)-based genomic classifier was implemented using TensorFlow/Keras for binary classification of BCG strains. Input genomic features (47 numerical columns, 1D vectors) underwent preprocessing: missing values were replaced by feature-wise medians, followed by standardization via StandardScaler. The architecture comprised two cascaded convolutional blocks—each with Conv1D layers (32 and 64 filters, kernel size = 3, ReLU activation), MaxPooling1D (pool size = 2), and dropout regularization (rates = 0.3 and 0.5)—followed by a flattening layer and two dense layers (32 ReLU units and 1 sigmoid unit). Training utilized the Adam optimizer (learning rate = 0.001) with binary cross-entropy loss over 30 epochs (batch size = 32), while addressing class imbalance through inverse-frequency class weights computed by compute_class_weight(). Model robustness was rigorously evaluated through stratified five-fold cross-validation to preserve label distribution across all partitions. Performance metrics (accuracy, precision, recall, F1, area under curve-receiver operating characteristic (AUC-ROC)) and composite scores (weighted averages with default equal weights) were quantified across all validation folds, with final results reported as mean ± standard deviation. The AUC-ROC analysis was performed using roc_curve() from scikit-learn, with mean ROC curves and standard deviation bands generated through interpolation at 100 evenly spaced false positive rate thresholds. Following comprehensive cross-validation, the final model was trained on the complete dataset and serialized as an H5 file for deployment.

Selection and evaluation method

The random forest model was selected for its transparent decision-making process, capacity to manage complex genomic datasets, and resilience to correlated genetic features, while its gene-ranking capability facilitated biological interpretation. The CNN framework was implemented to identify nucleotide-level signatures within linear genomic data, utilizing layered filters to detect localized sequence motifs. To mitigate model overfitting and address uneven class distribution, we incorporated randomized node deactivation (30–50% dropout rates) and adjusted sample weighting. These approaches were adopted based on their proven synergy in microbial studies—Random Forest for pinpointing key biomarkers and CNN for decoding intricate sequence hierarchies. Performance validation employed robust stratified k-fold cross-validation (five folds) with comprehensive metrics. The random forest’s parameters were fine-tuned using inbuilt out-of-bag error estimates, followed by rigorous verification through iterative gene elimination and label permutation tests (1,000 shuffled datasets). For the CNN, performance was systematically evaluated through stratified k-fold cross-validation using multi-metric analysis—classification accuracy, precision-recall balance, F1 consistency, and ROC curve profiling. Biological plausibility of the identified gene sets was confirmed via functional enrichment analysis against the COG and KEGG databases.

Training parameter optimization

To optimize the random forest models, we systematically optimized two key parameters: First, for the number of trees (ntree), we tested values from 1 to 10,000 (specifically 1, 10, 50, 250, 500, 1,000, 1,500, 2,000, 5,000, and 10,000), recording the out-of-bag (OOB) error rate and feature sparsity (proportion of zero-importance features) at each value. We observed that the OOB error rate decreased with increasing tree count and stabilized around 1,000–2,000 trees (Fig. S1), leading us to select ntree =10,000 for final models to ensure robustness through variance reduction without overfitting. Second, for features per split (mtry), we evaluated values including 1, approximately 10% of total features, 20% (1/5), 33% (1/3), 50% (1/2), and 100% of features. The minimal OOB error occurred at approximately 10% of features (Fig. S2), which simultaneously promoted feature sparsity and enhanced identification of biologically informative features. This dual optimization strategy balanced model performance with interpretability while ensuring computational efficiency.

Classification of BCG vaccine strains and non-BCG wild-type Mycobacterium bovis based on informative genes

We constructed a random forest classifier to differentiate strains based on their known labels (BCG vaccine or non-BCG wild-type). This approach identified a set of interpretable predictor genes associated with adaptation to each environmental context. Bitscore values of orthologous genes were used as input for training the random forest model, and its performance was assessed using OOB accuracy. To improve model performance, we implemented iterative feature selection: initially, all 1,052 orthologous genes meeting the selection criteria were used for training. Genes with variable importance (VI) equal to zero were pruned after the initial model training, and subsequent rounds of retraining were performed using only the top 50% of predictor genes until perfect OOB accuracy (100% accuracy) was achieved (Fig. 2A).

In the first round of model building, only 1,510 genes exhibited variable importance (VI) values that were significantly higher than those of the remaining genes (Fig. 2B). In contrast, 2,085 orthologous groups had a VI of zero, meaning these genes did not contribute to improving model accuracy and were excluded from further feature selection. This iterative process led to the development of Model 7, which achieved perfect classification accuracy for source prediction. Consequently, Model 7 was selected, with 47 top predictor genes identified as being most informative for distinguishing between the two groups of strains (Excel File S2). A heatmap illustrating the clustering of these 47 genes based on their bitscore matrix values is shown in Fig. 3.

Function analysis of top predictor genes

The 47 top predictor genes were assigned to 17 COG categories based on functional annotation. Apart from those with function unknown (S), a large proportion of them was involved in five COG categories, namely coenzyme transport and metabolism (H), lipid transport and metabolism (I), general function prediction only (R), signal transduction mechanisms (T) (Fig. 4).

The KO (KEGG Orthology) pathway analysis showed that the top predictor genes are mainly involved in three pathways: Metabolic pathways, biosynthesis of secondary metabolites, and Two-component system, starch and sucrose metabolism, biosynthesis of cofactors and Microbial metabolism in diverse environments (Table 1). The KO Brite analysis indicated that the predictor genes are primarily associated with two functional classifications: KO and enzymes (Table 2).

Performance evaluation of CNN-based genomic classifier using stratified five-fold cross-validation

The CNN-based genomic classifier demonstrated robust performance through stratified five-fold cross-validation, achieving a mean AUC of 0.964 (±0.046) with consistent discriminative capacity across all folds (Fig. 5). Where the high area under the curve indicates excellent separability between classes across all classification thresholds. The steep initial ascent of the curves reflects strong true positive rate (sensitivity) at low false positive rates, clinically valuable for minimizing false positives. Performance remained consistent across all folds, with high accuracy (95.8 ± 3.4%) and perfect recall (100.0 ± 0.0%). Precision (91.0 ± 7.4%) and F1-scores (95.1 ± 4.0%) showed minimal variability between partitions. The composite score (0.921 ± 0.038), calculated as a weighted average of key metrics, confirmed balanced performance without evidence of overfitting. Based on cross-validation stability and biological coherence of the extracted features, the final model was trained on the complete dataset and saved as final_cnn_bcg_classifier.h5 for deployment.