Available E. coli genomes (n = 31,230) were retrieved from the GenBank database during the database construction. Based on the genome sequence completeness (full E. coli genomes were included in priority), the country of isolation and the E. coli pathotype, 1,425 genomes were selected to maximize the diversity. Sixty-eight additional genomes sequenced and assembled in a previous study by Jaudou and colleagues (Jaudou et al., 2023), were added to reach a total of 1,493 genomes. The genome accession numbers and the metadata associated with the selected genomes are reported in Supplementary Table 1. All the genomes were screened against a custom database (Available at https://github.com/fabgenomics/ML_EHEC) containing all the stx subtypes, eae and O-group genes using abricate v1.0.1. (https://github.com/tseemann/abricate). The phylogroup of each genome was determined using the EzClermont phylotyping tool available at https://github.com/nickp60/EzClermont.

The selected genomes were annotated using the rapid prokaryotic genome annotation software prokka v1.13.3 (Seemann, 2014) using the proteins option with the reference genome of E. coli O157:H7 str. Sakai (NC 002695.2) (Hayashi et al., 2001). Resulting General Feature Format (.gff) files were processed through a Pangenome analysis pipeline using panaroo v1.2.7 with –clean-mode strict, –remove invalid-genes and –merge paralogs options (Tonkin-Hill et al., 2020). Panaroo collapses genes into putative families with a family sequence identity level of 70% by default and creates groups. The gene_presence_absence.Rtab table provided by the panaroo output contains all the groups and genes and the presence/absence information from each genome. These groups and genes were renamed using a custom script (Available at https://github.com/fabgenomics/ML_EHEC) in which we used the information in the panaroo output file matrix (gene_presence_absence.csv available at https://doi.org/10.5281/zenodo.7129021) to retrieve the corresponding locus tag number from the E. coli O157:H7 Sakai strain annotation (ECs number) relative to the group or the gene (CDS). Only groups renamed with the locus tag number by our custom script were retained for further analysis. We created a CDS presence/absence matrix (ECs_presence_absence.csv available at https://doi.org/10.5281/zenodo.7129021), by adding a new column based on the pathotype of each genome. Genomes that were found to be stx+/eae+ were assigned to the EHEC pathotype and noted 1 in the table and all the other genomes (STEC for stx+/eae-, EPEC for stx-/eae+ and the abbreviation COM was used for stx-/eae-) were considered non-EHEC pathotype and noted 0. A phylogenetic tree was reconstructed using IQtree v2.0.3 (Minh et al., 2020) with the Generalised Time Reversible (GTR) model on the core genome alignment produced by panaroo. The tree was annotated with the molecular serogroups using the CLC genomic workbench v21 (QIAGEN, Aarhus, Denmark).

Before evaluating the performance of the different classifiers, the CDS presence/absence matrix was filtered on non-informative features. Any CDS with less than 10% variance on its presence/absence vector was removed, as these loci do not contain useful information to test machine learning algorithms. This step removed 3,603 CDS, resulting in a dataset with 1,178 CDS. Then, to avoid possible data leakage between training and testing datasets, we grouped the samples based on their similarity. For each pair of samples, their CDS presence/absence vector was used to compute a hamming distance (i.e. the number of differences). Any two samples with a hamming distance lower than or equal to D were allocated to the same cluster. We considered possible values for D of 5, 10, 50, 100 and 200. Note that the resulting clusters have extremely homogeneous pathotypes (at D = 100, only one cluster consists of a mixture of EHEC and non-EHEC samples). The cluster table is available at https://github.com/fabgenomics/ML_EHEC. The main dataset (n = 1,493) was then subsampled to keep only one genome from each cluster. Then, each subsampled dataset was randomly split using 80% of the samples for training/validation and the remaining 20% for testing. The train_test_split module from Sklearn was used, with stratify option to control the proportion of EHEC in both datasets. Eight classification algorithms were trained on each dataset using 10-fold Cross-Validation: Decision Tree, Extra Tree, Gradient Boosting, LGBMClassifier, Logistic Regression, Random Forest, XGBClassifier and Support Vector Machine. The evaluation metric used was the function cross_val_score from the sklearn library. For all the cross-validation scores (10 folds i.e. 10 scores), the mean accuracy was calculated (Supplementary Table 2 and Figure 1). Further analysis were performed using a distance D = 100 for clustering (dataset Cluster-D100). This implies that two samples in this dataset differ by at least 8.4% (100/1178 CDS) of their gene content. A module from Sklearn library v0.23.1 (RandomUnderSampler) was used to select randomly the same amount of non-EHEC genomes to be equal to the number of remaining EHEC genomes from the cluster analysis. The Cluster-D100 dataset was randomly split with a ratio of 80/20% for training and testing datasets respectively and the stratify option. Eight classification algorithms were used to select the most important features with the SelectFromModel library from Sklearn. The most important CDS are listed in Table 1. We arbitrarily chose to select the six most important features to create a new reduced dataset. With this resulting matrix, hyper-parameter tuning was done on each of the eight models using RandomizedSearchCV and GridSearchCV (scoring on roc_auc metric) and cross-validation steps (n = 5) when the option was available. Finally, each classifier was retrained with its best hyper-parameter and evaluated on the testing dataset previously set aside (accuracy, precision, recall and F1-score are obtained using the classification_report from Sklearn, see Supplementary Table 3). To understand which gene combination led to the prediction of the EHEC pathotype, we generated all 2⁶ = 64 combinations of the six genes presence/absence and computed, for each ML model, the probability of the EHEC pathotype. We then kept the cases where the probability was ≥0.7 and transformed the set of gene combinations into simplified boolean expressions (using a boolean Algebra Solver). The results are reported on Figure 2 Charts of the training pipeline and the prediction pipeline are presented in Figures 3A, B, respectively.

From a previous study conducted by Jaudou et al. (2023), raw ONT MinION reads from two STEC (ECA279 (O174:H2) = SRR18191627 and 97HMPL652 (O110:H9) = SRR18191587), one stx-negative eae-positive E. coli (2142-O103 (O103:H25) = SRR18191529), one eae-positive STEC (E. coli 12-1 (O157:H7) = SRR18191640) and one commensal E. coli [i.e., negative for both stx and eae (NC809 (O41:H7) = SRR18191621)] were collected. One stx-negative eae-positive E. coli strain (KK072/05 - O156:H8) was newly sequenced during this study following the same protocol described by Jaudou et al. (Jaudou et al., 2023) and the raw ONT MinION reads were deposited to the NCBI database under the number SRX18376762. Raw reads were subsampled using Rasusa v0.6.0 (Hall, 2022) with 5.5Mb for the target number of bases and 10, 20, 30, 40 and 50x for the coverage. From these subsampled reads, twenty-five in silico mixtures were generated by concatenating into a 1:1 ratio the same coverage of subsampled reads. Details of the different mixtures are presented in the Table 2. From these mixtures, an assembly was generated using metaFlye v2.9-b1768 (Kolmogorov et al., 2020) with nano-raw and meta parameters. The resulting assemblies were annotated with the same parameters as described in the Genome annotation and classification paragraph (Section 2.2). The produced GFF file of each annotation was integrated individually in the pangenome graph generated with the 1,493 genomes using the panaroo-integrate command from the panaroo program. From the new gene_presence_absence.csv and the gene presence absence.Rtab generated by the panaroo-integrate command, the row corresponding to the added mixture was extracted using a newly developed python script (Available at https://github.com/fabgenomics/ML_EHEC) to reconstruct the CDS presence/absence table. We extracted only the features required for the tested model and when a feature was absent, we created it and introduced a 0 value. The data extracted were used to perform the predictions (Figure 3B). The predict_proba method from all the algorithms was used to estimate the probability that the sample is an EHEC.

Eight metagenomes from artificially contaminated raw milks described in a previous study (Jaudou et al., 2022) were downloaded from the Genbank public database (Table 3). The estimated level of contamination was 0.5x10³, 0.5x10² and 0.5x10¹ CFU.mL⁻¹ of EHEC O26 plus one EHEC-free milk. Raw reads were processed using the STECmetadetector pipeline developed by Jaudou et al. (2022) available at https://gitlab.com/Bfr_bioinformatics/STECmetadetector and the extracted E. coli reads were assembled using metaFlye v2.9-b1768 with the same parameters as described in Section 2.4. The resulting assemblies were annotated with the same parameters as described in the Genome annotation and classification paragraph (Section 2.2). The resulting GFF file was treated with the same process than the in-silico mixture GFF file and the EHEC predictions were performed, as described in Section 2.3 (Figure 3B).

To take advantages of ML to find patterns and preserve its generalization potential, we constituted an E. coli database for which we selected in priority complete E. coli genomes with a minimum of required metadata (origin, isolation date and location) and verified their pathotypes (stx and eae genes presence). A total of 1,493 genomes were downloaded from the GenBank database. The geographic origin of the strains was 33, 26, 23, 7, and 3% from Europe, Asia, America, Africa and Oceania, respectively and the 8% remaining were missing the country of origin. During the genome selection, we tried to respect an equal proportions of stx/eae-positive strains (i.e., EHEC) and non-EHEC strains based on the metadata provided. Genomes simultaneously positive for at least one stx gene and the eae gene (n = 632) were assigned to the EHEC group. The other genomes (n = 861) were assigned to the non-EHEC group (Available at https://doi.org/10.5281/zenodo.7129021). In addition, the custom database was reporting all the O-group sequences so that the serogroup information, in particular the most frequent EHEC serogroup, was available (Supplementary Table 1). The top seven most represented serogroups were O26 (n = 160), O157 (n = 126), O103 (n = 66), O121 (n = 61), O145 (n = 60), O111 (n = 36), and O45 (n = 9). These serogroups comprised mostly EHEC strains with 126 EHEC O26 strains, 124 EHEC O157 strains, 57 EHEC O103 strains, 55 EHEC O121 strains, 59 EHEC O145 strains, and 31 EHEC O111 strains. The phylogroup analysis showed that the diversity of the species is well represented. The phylogeny of the final dataset is illustrated in Figure 4.

From 37,380 groups generated by panaroo, 13,952 remained with an ECs annotation. Some ECs were duplicated in the table due to the genetic diversity of some genes and the panaroo identity level threshold. We aggregated the results for these genes, which resulted in 4,780 unique CDS and kept the presence/absence information. Before splitting the ECs presence/absence table (Available at https://doi.org/10.5281/zenodo.7129021), we first selected CDS with enough variation between the samples (see Section 2). This filtering step removed 3,602 CDS resulting in a dataset containing 1,178 CDS. To avoid optimistic performance estimates that would result from near duplicate samples both in training and testing sets, we performed a clustering based on gene content similarity (see Section 2). Clustering at distance D ensures that two clusters in the dataset have at least D genes that are different. We evaluated the change in accuracy for increasing values of D, starting from 5 and up to 200 (Figure 1). Up to D = 100 (8.4% of the CDS), the accuracy of all classifiers remains very high (above 97%). This indicates that robust information can be extracted from gene profiles to predict EHEC pathotype with high accuracy. The performances drop by around 5–10% for D = 200 for an accuracy of around 93%. While still high, this decrease could be attributed to the low number of clusters that remain for prediction (n = 137, 15 EHEC and 122 non-EHEC). From these results, a value of D = 100 was chosen to build the dataset for further analysis, as it combines both good performance with a sufficient sample size (n = 756). However, this dataset was imbalanced. To avoid overfitting problems, we subsampled the non-EHEC class to be equal to the EHEC class and ended with a balanced dataset (see Section 2). The matrix was then randomly split into a training dataset containing 80% of the genomes used in the study (n = 139) with 69 non-EHEC and 70 EHEC. The testing dataset contained the remaining 20% (n = 35) with 18 non-EHEC and 17 EHEC.

We conducted a comprehensive analysis of the top features extracted from the training of eight classification models, and the results are summarized in Table 1. The table shows the top six most important features, ranked based on the number of times they were used by the models. The most important feature was found to be ECs_1056, which corresponds to a phage excisionase gene. This feature was used by eight of the models. The second most important feature was ECs_1812, which corresponds to the nleA/espI gene, coding for a type III secretion system (T3SS) secreted effector protein, and was used by seven of the models. The third most important feature was ECs_1824, which corresponds to the nleG gene, which encodes another T3SS secreted effector protein, and was used by five of the models. The remaining features, namely ECs_3858, ECs_1815, and ECs_1561, were used by four models each, and correspond to the nleE, nleF, and espN genes, respectively, all of which encodes T3SS secreted effector proteins. Multiple classifiers achieved near perfect performance when evaluated on the D100 dataset. To better understand which combinations of genes contribute to the prediction of the EHEC pathotype, we generated all 64 (2⁶) genes presence/absence profiles and recorded in which case one of the models predicted an EHEC pathotype with confidence (Figure 2 and Section 2). This confirms that some genes, such as ECs_1056 (phage excisionase) and ECs_1824 (nleG) have an importance (they are needed in the majority of the predictions). Decision Tree and Gradient boosting learned the same decision rule: “ECs_1812 (nleA) and (ECs_1056 or (ECs_1815 (nleF) and ECs_1824))”. The SVM classifier predicts EHEC for any two combination of ECs_1056, ECs_1815, and ECs_1824. All those classifiers can predict an EHEC pathotype with as little as two genes. On the other hand, Extra tree and logistic regression make predictions involving three to four genes, showing that they can have different sensitivity.

Supplementary Table 3 shows the evaluation metrics obtained by training the eight classifiers using the selected six features: ECs_1056, ECs_1812, ECs_1824, ECs_3858, ECs_1815, ECs_1561. All classifiers achieved high accuracy scores, ranging from 0.97 to 1.00, indicating a good performance in predicting the target variable. Logistic Regression, Extra Trees XGBCLassifier, LGBMClassifier, Decision Tree, and SVM achieved perfect accuracy scores of 1.00, while Random Forest and Gradient Boosting achieved a slightly lower accuracy score of 0.97. Extra Tree and Random Forest achieved a precision of 0.94 on the EHEC class and a recall of 0.94 on the non-EHEC class. All the other classifiers achieved perfect precision, recall, and F1-scores, indicating that the selected features were informative and sufficient to discriminate perfectly between the two classes.

We first tested the ability of the different models to predict the presence of an EHEC strain in a mixture of E. coli strains. For this purpose, in silico mixtures of raw MinION reads were assembled. Assemblies produced using meta-Flye ranged from 5.63 to7.94 Mb (mean = 6.57 Mb) and the number of contigs from 80 to 155 (mean = 106) (Supplementary Table 4). In all mixtures of E. coli strains, the assembly size was longer than a normal E. coli assembly (4.8–5.5 Mb). Shorter assemblies were produced by meta-Flye with the STEC-COM mixture (5.63–6.16 Mb). On the contrary, the EHEC-COM mixture produced longer assemblies (7.35–7.94 Mb). The eight models were then used to perform predictions on the in silico mixtures (Table 2). Predictions of the pathotype ranged from 0 to 1 and the average prediction from 0.01 to 0.99. The cut-off for binary classification is usually set to 0.5. Above or equal to this cut-off value the presence of an EHEC is predicted, and below a non-EHEC is predicted. This cut-off of 0.5 was used to report the results presented in this study. The eight classifiers were able to predict the correct class 22 times over 25 predictions (88%). However, all models incorrectly predicted the presence of an EHEC three times for the STEC-EPEC mixture with the strain 97HMPL652 and 2142-O103 for a coverage of 30x, 40x and 50x, respectively. For all non-EHEC containing mixtures, the higher value was 0.47 for the same STEC-EPEC mixture with the Extra Tree classifier (Table 2). For the EHEC class the lower value was 0.75 for the EHEC-COM mixture and the STEC-EHEC mixture with the Decision Tree classifier and the XGBClassifier. Taken together, these data indicate that all the classifiers were able to predict with high confidence the presence of an EHEC in E. coli mixtures that combine different E. coli pathotypes. Only three false positives were predicted for the most difficult mixture combining a STEC and an EPEC.

We then tested the performance of the eight models on complex mixtures using artificially contaminated raw milk. A bioinformatic pipeline called STECmetadetector developed by Jaudou et al. (2022), was used to specifically extract E. coli reads from raw milk samples artificially contaminated with an O26 EHEC strain (from 0 to 500 CFU.mL⁻¹). The E. coli reads were assembled using the same method as described in Section 2 (paragraph 2.5). Assemblies ranged from 5.2 Mb to 6.83 Mb (mean = 5.95 Mb) and the numbers of contigs from 5 to 62 (mean = 15). The eight classifiers were used to predict the pathotype of artificially contaminated raw milks (Table 3). All models were able to accurately predict the EHEC pathotype in the sample with high confidence, at the three contamination levels tested, regardless of the strain used for the artificial contamination. Notably, the raw milk used for spiking with the 6423-O26 strain was naturally containing commensal E. coli of serotype O185:H2 and O8:H19 (Jaudou et al., 2022). Predictions of the pathotype ranged from 0.95 to 1 for the class EHEC for all contamination levels. The negative control (a non-contaminated raw milk), was classified accurately as non-EHEC by all eight classifiers with the higher value of 0.21 for the SVM (Table 3).