A large scale evaluation of TBProfiler and Mykrobe for antibiotic resistance prediction in Mycobacterium tuberculosis

Dataset constitution

A dataset of 6,616 genomes with reads available on the sequence read archive (SRA) and resistance phenotypes was gathered from several studies (Coll et al., 2018; Bradley et al., 2015; Schleusener et al., 2017; Vincent et al., 2012) and from the ReSeqTB website (https://platform.reseqtb.org/), as of June 2017. After removal of duplicated genomes, for which the raw reads had the same nucleotidic content, the final dataset included 6,574 genomes. All sequencing experiments were based on the Illumina technology and their great majority involved Illumina HiSeq or Miseq instruments.

The resulting dataset contains a variable number of phenotypes for 20 antibiotics. The current study focused on the 10 antibiotics addressed by TBProfiler and/or Mykrobe, namely amikacin, capreomycin, ethambutol, ethionamide, fluoroquinolones, isoniazid, kanamycin, pyrazinamide, rifampicin, and streptomycin. Although addressed by TBProfiler, para aminosalicylic acid was not considered due to the limited number of strains with phenotypes available (45 resistant for 368 susceptible ones). Moreover, since TBProfiler and Mykrobe both provide a prediction for the entire family of fluoroquinolones, a composite phenotype was created by aggregating the phenotypes of ofloxacin, moxifloxacin and ciprofloxacin: it was set as the corresponding phenotype when they agreed, and set as missing when they did not. For each antibiotic, the number of strains ranges from 635 (ethionamide) to 6,464 (isoniazid), with a small fraction of resistant strains ranging from 14.7% to 37.8% (see Table 1 for further details).

Command-line versions of TBProfiler (version 0.3.4) and Mykrobe (v0.3.3-0-gc211bf2) were run on these 6,574 samples, using their default configurations (see Supplementary Materials for further details). A resistance genotype matrix, encoding the presence/absence of all the mutations detected in the entire panel, was then built for each software. In total, 469 distinct mutations were detected by TBProfiler out of the 1,195 candidate mutations embedded in its catalog. In contrast, Mykrobe, which considers a smaller catalog detected 272 distinct mutations. Table 1 compares, for each software and each antibiotic, the number of detected and candidate markers. The difference can be important in some cases. Given the large number of samples considered in this study, the predictive power of the undetected mutations is most likely to be limited.

We noted that TBProfiler provides the frequency with which a given marker was detected, defined as the fraction of reads presenting the resistance marker among reads mapped at this position. A similar information can be computed from the results of Mykrobe, which reports the number of hits obtained for a given marker and its corresponding reference allele. This frequency information allowed us to investigate whether considering a minimum frequency threshold to call a marker present could be beneficial in terms of prediction or not.

Interestingly, as a by-product of its genotyping process, TBProfiler detects novel mutations within the 33 resistance loci it considers, with respect to their reference sequence in the H37Rv genome. As can be seen from Table 1, the number of these novel mutations can be relatively high (22,883 in total, but more than 50% observed in a single genome). While not used to predict resistance, their predictive power can be empirically measured, as will be described later on.

As an ending remark, we note that TBProfiler and Mykrobe did not provide results for three and two samples, respectively.

Multivariate modeling

Motivated by the observation that many markers embedded in TBProfiler and Mykrobe catalogs seem to individually have a limited predictive power, we aimed to evaluate whether alternative multivariate machine learning approaches operating from the set of markers detected by these softwares could improve over the direct-association strategy. For this purpose, we represented each genome by a vector of binary variables encoding the presence or absence of the mutations detected by each tool on the entire panel of genomes, and relied on the Lasso-penalized logistic regression to learn a resistance prediction model for each antibiotic. The logistic regression model is appealing in this context since it provides a probabilistic prediction, thereby allowing to measure the confidence of the prediction, and to control the trade-off between sensitivity and specificity by adjusting its decision threshold. Coupled with the Lasso penalty, it leads to sparse solutions allowing to identify key (combinations of) mutations, hence to interpretable models.

To estimate the predictive performance of the Lasso-penalized logistic regression model, we relied on a nested 10-fold cross-validation procedure, that is, a 10-fold cross-validation procedure with an inner optimization of the regularization parameter, based on the AUC criterion. The final models were similarly obtained by applying a standard cross-validation procedure, designed to maximize the AUC, on the entire dataset. To test whether a difference observed in terms of AUC was significantly different from zero, we relied on the non-parametric approach proposed by DeLong, DeLong & Clarke-Pearson (1988), which is implemented in the R package pROC (Robin et al., 2011). In particular, the statistical significance was assessed at the 0.01 level.

Data availability

All data used and obtained in this study are available as Supplementary Materials. They contain in particular, (i) a table containing the reference phenotypes and the SRA accessions, (ii) the resistance genotype matrices obtained by TBProfiler and Mykrobe, (iii) the inferred lineages, (iv) the predictive power of individual markers, and (v) the final models by the Lasso-penalized logistic regression models (i.e., mutations selected and beta coefficients in the corresponding multivariate logistic model). Note that accessing and using data originating from the ReseqTB platform is subject to the following terms and conditions: https://platform.reseqtb.org/main/acceptTerms.html.

Inferred lineages

A very high level of consistency between lineages inferred by both softwares was observed, consistently with a previous study (Schleusener et al., 2017). We noted however that TBProfiler provided “mixed” lineage calls for 228 samples out of the 6,571 successfully processed. These ambiguities were especially observed among the West-African, bovis and bovis/African lineages, for which TBProfiler systematically provided mixed results. Still, between 1.8% and 5.8% of ambiguous results were observed among the four major lineages (lineages 1–4). Further analyzes revealed however that in many cases one of the lineages was called from a marker detected at a low frequency, and that including a minimum frequency threshold to detect a lineage allowed to significantly decrease the number of ambiguous results observed among main lineages. The above rates decreased indeed respectively to 0.5–2% and 0.1–0.8% when a minimum threshold of 0.1 and 0.2 was considered. Among the 6,342 genomes unambiguously classified by both softwares, they disagreed on 21 samples only, among which 19 were called “unknown” by one software or the other. Moreover, more than 99% of the samples originate from four main lineages: lineages 1 (East-African/Indian Ocean—9.5%), 2 (Beijing/East-Asia—15%), 3 (Delhi/Central-Asia—16%), and mainly 4 (European/American—59%). This is illustrated in Figs. S1–S3.

Resistance genotyping agreement

To assess the agreement of the softwares in terms of resistance, we first evaluated their agreement at the marker level, that is, their ability to detect the same markers within the same panel of genomes. For this purpose, we extracted a list of 116 mutations embedded in both catalogs and compared the number of calls made by both softwares. We noted that the number of calls differed quite significantly in some cases, as illustrated in Fig. 1. For instance, three markers detected in a single strain by TBProfiler are detected in 10 strains or more by Mykrobe, with one being detected in almost 30 strains. Likewise, a marker detected in around 20 strains by TBProfiler is detected in more than 70 strains by Mykrobe. Conversely, however, a marker detected in around 400 strains by TBProfiler is detected in around 35 strains less by Mykrobe. It is worth noting that these results were obtained without considering any minimum frequency threshold for calling a marker present. Similar results were also obtained when considering a threshold of 0.5 (see Fig. S4), indicating that the observed discrepancies did not result from calling present markers observed at a low frequency. This analysis therefore revealed that the discrepancies in R/S predictions made by both softwares may not only be due to differences in the catalogs of mutations they embed, but also to their underlying genotyping pipelines.

We then compared their overall agreement in terms of R/S prediction. For this purpose, we worked from the 6,570 samples and the nine antibiotics for which both tools provided a prediction (Mykrobe does not predict resistance to ethionamide). Overall, TBProfiler and Mykrobe agreed on 95.4% of the 28,426 reference phenotypes available, with 92.4% of agreement for the 6,451 resistant reference phenotypes, and 94.1% for the 21,975 susceptible ones.

Predictive performance

Table 2 gives a summary of TBProfiler and Mykrobe predictive performance, measured in terms of the following indicators:

• sensitivity: the fraction of resistant strains predicted as such (also called recall).

• specificity: the fraction of susceptible strains predicted as such.

• precision: the fraction of resistant strains among strains predicted as resistant (also called positive predictive value).

We first noted that performance differed between both tools, TBProfiler offering in general a higher sensitivity for a lower specificity. This is mainly due to the fact that TBProfiler contains more resistance-defining mutations than Mykrobe, and that the prediction strategy amounts to calling a strain resistant whenever one of them is detected. Considering the average of sensitivity and specificity, referred to in the following as macro-accuracy, TBProfiler achieves a better performance than Mykrobe for five drugs (capreomycin, ethambutol, kanamycin, pyrazinamide, and fluoroquinolones). Interestingly however, Mykrobe outperformed TBProfiler in terms both of sensitivity and specificity for streptomycin (by 4.9 points of macro-accuracy), the sole antibiotic for which it detected more mutations than TBProfiler, especially within the gid gene (41 mutations instead of 1).

We also noted that the performance varied across antibiotics. Sensitivity was only around 80% for streptomycin and capreomycin, and even lower for pyrazinamide (59.1% and 34.4% with TBProfiler and Mykrobe, respectively). Specificity generally exceeded 90%, but precision was quite low in some cases (e.g., around 57% and 68% using TBProfiler for ethionamide and ethambutol, respectively). For ethambutol and TBProfiler, for instance, this means that while 92.5% of the R strains were indeed identified as such, only 67.8% of the strains predicted to be resistant were actually resistant. This might seem contradictory with the fact that the specificity was quite high (92.5%, meaning that 92.5% of the S strains were recognized as such), but was simply a consequence of the fact that the dataset involved around five times more S than R strains. While the precision (or positive predictive value) is an indicator of great clinical interest, since it is helpful to interpret the results actually provided by the software, it intrinsically depends on the prevalence of R strains in the panel studied hence must be interpreted with care (Zignol et al., 2018).

Strikingly, we observed major differences between our performance estimations and that reported in previous studies. For instance, streptomycin sensitivities differed by more than 20 points from the ones reported in Schleusener et al. (2017), 57% for both TBProfiler and Mykrobe, instead of 77.4% and 81.4%, respectively, in the current study. Likewise, performance differed by more than 10 points for isoniazid (84% specificity for TBProfiler in Schleusener et al. (2017) instead of 97%, and 79% sensitivity for Mykrobe instead of 88.5%). These differences can be explained by the small sample size, 91 strains considered in Schleusener et al. (2017), leading to uncertain performance prediction estimations.

We noted finally that the minimum frequency threshold considered to call a marker present sometimes had an impact on the predictive performance. This was especially the case with TBProfiler for ethambutol, pyrazinamide and rifampicin, where an improvement of up to three points of macro-accuracy was obtained for rifampicin, as shown in Fig. S5. Mykrobe, on the other hand, showed a more stable behavior, with marginal improvements obtained for streptomycin and pyrazinamide. Accordingly, the results presented in Table 2 were obtained by optimizing this minimum frequency threshold, for each software and each antibiotic, to maximize the corresponding macro-accuracy. An identical table was included in Table S1 when no such minimum frequency thresholds were considered. In practice, we recommend to use the thresholds provided in Table 2.

Predictive performance by lineage

The size of the dataset considered allowed then to evaluate how TBProfiler and Mykrobe behaved across lineages. For this purpose, we focused on the four main lineages, that is lineages 1–4, and computed the same performance indicators per lineage. Table 3 summarizes the results obtained for TBProfiler. The global entries correspond to the performance obtained from the 6,342 samples coming from these four major lineages, and excluding all samples with ambiguous lineage calls. These values can slightly differ from the values in Table 2, obtained from the entire dataset of 6,571 samples, including samples of other lineages and ambiguous lineage calls. As before, the minimum frequency threshold to call a marker present was fine-tuned, for each drug, to maximize the macro-accuracy indicator on the entire dataset (i.e., of the global entry).

We noted a large discrepancy across lineages for some antibiotics. This was notably the case of capreomycin, ethambutol, ethionamide, pyrazinamide, and streptomycin, where a difference of more than 10 points of macro-accuracy was observed between the lineages of highest and lowest performance. This difference reached 46.4 points for capreomycin and 34.1 points for pyrazinamide, where no resistant strain from lineage 1 was detected in both cases, and exceeded 17 points for ethionamide and streptomycin. Interestingly, performance observed within lineage 4 was not systematically better than within other lineages, despite the fact that it accounted for more than 50% of the samples for all drugs but ethionamide, where it represented 44% of the samples. We also noted that a lower performance could be due either to a lack of sensitivity, as it was the case for amikacin/lineage 4, capreomycin/lineage 4 and streptomycin/lineage 4, or to a lack of specificity, for ethambutol/lineage 2 and ethionamide/lineage 2. Similar observations could be made from the results obtained by Mykrobe, shown in Table S2.

Individual markers performance

As mentioned above TBProfiler and Mykrobe consider a pre-defined list of resistance-defining mutations and predict a strain resistant whenever one of these mutations is detected. The number and the quality of these markers therefore intrinsically define the overall performance of the softwares. Analyzing individually the predictive power of the markers is interesting in at least two respects. First, to get a better understanding of the discrepancies observed previously in terms of the overall predictive performances of the softwares. Moreover, as described by Miotto et al. (2017) and suggested earlier by Farhat et al. (2016), to appreciate the prediction provided by the algorithm: the prediction may be deemed more confident when it is based on a marker with a higher predictive power.

We first noted that a limited number of markers had a high level of sensitivity. Considering all drugs, 80.6% and 70.6% of the markers detected by TBProfiler and Mykrobe, respectively, had a sensitivity smaller than 1%. Only 19 and 14 of TBProfiler and Mykrobe markers had a sensitivity above 10%, and only seven above 40% in both cases. In contrast, the great majority of markers were highly specific: 98.6% of the TBProfiler and all Mykrobe markers had a specificity greater than 95%. This is illustrated in Fig. S6.

Figure 2 shows the individual markers performance obtained for amikacin and streptomycin to illustrate how they can help better understanding the trade-offs achieved by both softwares. For amikacin, we noted that the predictive power was mostly driven by a specific mutation in both cases. Table S3 indicates that it corresponded to the mutation rss_1401A > G, which was captured by both softwares, albeit in a different number of strains (339 for TBProfiler and 311 for Mykrobe). This mutation taken individually would lead to a sensitivity of 86.7% and a specificity of 98.2% for TBProfiler, and of 80.7% and 98.6% for Mykrobe. All the other mutations considered by Mykrobe are highly specific but not very sensitive (2.7% at best, see Table S3), hence only allowing for a marginal improvement over the sensitivity brought by the rss_1401A > G mutation. TBProfiler, on the other hand, includes a mutation with an intermediate level of sensitivity, not belonging to the Mykrobe catalog. It corresponds to a mutation at position 514 in rss, and while it indeed allows improving the overall sensitivity of TBProfiler, it has a specificity lower than any mutation considered by Mykrobe for amikacin (96.8% instead of 98.6%, see Table S3). Including this mutation improves TBProfiler sensitivity but at the cost of a reduced specificity. Choosing to include or not this mutation amounts to favoring sensitivity over specificity, and is therefore a matter of choice regarding the target performance one wants to achieve. In the case of streptomycin, Mykrobe is both more sensitive and specific than TBProfiler. Interestingly, Fig. 2 suggested a similar behavior to that of amikacin, TBProfiler including a marker of intermediate sensitivity and lesser specificity. While this can explain the lower specificity of TBProfiler, this marker does not allow TBProfiler to outperform Mykrobe in terms of sensitivity in this case. This therefore suggests that the Mykrobe catalog is more thorough for this antibiotic. We note from Table 1 that Mykrobe has detected 51 mutations for this antibiotic, instead of 27 by TBProfiler. It turns out that Mykrobe contains much more mutations in gid than TBProfiler (41 instead of one). We note however, that some novel mutations detected by TBProfiler may have a good predictive power. For streptomycin, this is in particular the case of a highly specific mutation (99.7%) having a sensitivity of 7.7% which, if not already explained by other markers, may bring additional predictive power. As can be seen on Fig. 2, only three markers of the native TBProfiler catalog actually exhibit a greater sensitivity. These novel mutations are not taken into account by TBProfiler in its prediction process, but can readily be considered in multivariate machine learning algorithms, as will be described below.

A complete table providing the individual performance of the known and novel mutations is available as Supplementary Materials to allow a finer exploration.

Multivariate modeling strategies

Models based on the original TBProfiler and Mykrobe catalogs of markers

We first built models from the markers included in the TBProfiler and Mykrobe catalogs, separately, using the nested cross-validation procedure described above to estimate the generalization performance of the models. Note that we considered only the markers related to a given antibiotic to build the corresponding prediction model.

Table 4 summarizes the results obtained in terms of AUC and support size, defined as the number of markers involved in the final model. Performance was comparable for most antibiotics, but better AUCs were obtained by the model learned from TBProfiler markers for kanamycin and pyrazinamide (significantly at the 0.01 significance level, see Methods). Interestingly, the number of markers involved in the models was relatively limited, with often less than 20 markers. This represents a drastic reduction with respect to the size of the original catalogs. ROC curve analysis was performed to objectively compare the performance of these models and that of TBProfiler and Mykrobe. The resulting ROC curves for four antibiotics are illustrated in Fig. 3, the curves obtained for the remaining ones are shown in Fig. S9. For all antibiotics, the performance of TBProfiler and Mykrobe lied on the ROC curves obtained for the corresponding logistic regression models. This means that a similar sensitivity/specificity trade-off can be achieved by adjusting the threshold on the probability returned by the logistic regression model, indicating that both strategies achieve the same level of performance. We noted moreover that TBProfiler and Mykrobe performance often lied where the ROC curves were the closest to the optimal point of the ROC space (which lies at x = 0 and y = 1 and denotes perfect sensitivity and specificity). The most notable exception was observed for amikacin and TBProfiler, where the logistic regression model allowed to trade a bit of sensitivity for a greater gain in specificity.

Table 4:

Performance obtained using Lasso-penalized logistic regression models operating from the markers detected by TBProfiler or Mykrobe.

	TBProfiler markers		Mykrobe markers
	AUC	Support	AUC	Support
Amikacin	92.4	6/9	89.5	2/7
Capreomycin	87.1	2/9	85.8	3/9
Ethambutol	92.4	21/72	91	11/24
Ethionamide	82.4	12/23	–	–
Fluoroquinolones	92.3	11/33	90.6	8/32
Isoniazid	94.2	18/100	92.7	7/44
Kanamycin	93.7	5/11	89.1	3/8
Pyrazinamide	76.1	60/154	65.3	20/56
Rifampicin	95	21/53	95.1	21/55
Streptomycin	87.7	9/27	88.3	45/51

DOI: 10.7717/peerj.6857/table-2

Note:

AUC was estimated by the nested cross-validation procedure described in the Materials and Methods section. The support corresponds to the number of markers included in the model, over the number of candidate markers (see Table 1). AUC values shown in bold are significantly different between the two models. Statistical significance is assessed at the 0.01 level.

Overall, while it did not improve the performance of TBProfiler and Mykrobe, the Lasso-penalized logistic regression strategy allowed to reach the same level of performance while providing a probabilistic prediction, hence allowing to adjust, to some extent, the desired levels of sensitivity and specificity, and relying on much shorter list of mutations.