Geographic heterogeneity impacts drug resistance predictions in Mycobacterium tuberculosis

The efficacy of antibiotic drug treatments in tuberculosis (TB) is significantly threatened by the development of drug resistance. There is a need for a robust diagnostic system that can accurately predict drug resistance in patients. In recent years, researchers have been taking advantage of whole-genome sequencing (WGS) data to infer antibiotic resistance. In this work we investigate the power of machine learning tools in inferring drug resistance from WGS data on three distinct datasets differing in their geographical diversity. We analyzed data from the Relational Sequencing TB Data Platform, which comprises global isolates from 32 different countries, the PATRIC database, containing isolates contributed by researchers around the world, and isolates collected by the British Columbia Centre for Disease Control in Canada. We predicted drug resistance to the first-line drugs: isoniazid, rifampicin, ethambutol, pyrazinamide, and streptomycin. We focused on the genes which previous evidence suggests are involved in drug resistance in TB. We called single-nucleotide polymorphisms using the Snippy pipeline, then applied different machine learning models. Following best practices, we chose the best parameters for each model via cross-validation on the training set and evaluated the performance via the sensitivity-specificity tradeoffs on the testing set. To the best of our knowledge, our study is the first to predict antibiotic resistance in TB across multiple datasets. We obtained a performance comparable to that seen in previous studies, but observed that performance may be negatively affected when training on one dataset and testing on another, suggesting the importance of geographical heterogeneity in drug resistance predictions. In addition, we investigated the importance of each gene within each model, and recapitulated some previously known biology of drug resistance. This study paves the way for further investigations, with the ultimate goal of creating an accurate, interpretable and globally generalizable model for predicting drug resistance in TB. Author summary Drug resistance in pathogenic bacteria such as Mycobacterium tuberculosis can be predicted by an application of machine learning models to next-generation sequencing data. The received wisdom is that following standard protocols for training commonly used machine learning models should produce accurate drug resistance predictions. In this paper, we propose an important caveat to this idea. Specifically, we show that considering geographical diversity is critical for making accurate predictions, and that different geographic regions may have disparate drug resistance mechanisms that are predominant. By comparing the results within and across a regional dataset and two international datasets, we show that model performance may vary dramatically between settings. In addition, we propose a new method for extracting the most important variants responsible for predicting resistance to each first-line drug, and show that it is to recapitulate a large amount of what is known about the biology of drug resistance in Mycobacterium tuberculosis.

Tuberculosis (TB) is a serious infectious disease caused by Mycobacterium tuberculosis 2 (MTB) and one of the 10 leading causes of death in the world, with roughly 1.4 million 3 people dying from the disease every year [1]. Based on the latest WHO TB report, TB 4 incidence is decreasing at approximately 2% each year, and about 16% of the cases of 5 the disease result in death. More efforts are needed to reach the targets of at most a 5% 6 case fatality ratio and at most a 10% TB treatment failure rate adopted by the End TB 7 Strategy [2]. In order to achieve these milestones, a lot of research has focused on 8 addressing one of the most pressing issues in TB, drug resistance. MTB is susceptible to 9 acquiring resistance, which poses challenges for an effective treatment of the disease. 10 600,000 new cases of resistance to rifampicin, one of the key first-line drugs, were 11 reported in 2016, and 490,000 of these cases were multi-drug resistant, defined as 12 resistance to rifampicin and isoniazid, the two most effective first line-drugs [1]. Drug 13 resistance arises and continues to spread largely due to two reasons: improper TB drug 14 treatment, whereby poor quality drugs, premature treatment discontinuation, and 15 inappropriate drug regimens apply selective pressure on the MTB bacteria within 16 patients that leads to drug resistance, and person-to-person transmission, whereby these 17 drug-resistant MTB strains are spread to new patients [1]. 18 The gold standard in detecting drug resistance is phenotypic drug susceptibility 19 testing via bacterial culture, whereby MTB is cultured in a Lowenstein-Jensen medium 20 and exposed to a particular drug, and its response to the drug is then observed. This 21 process is slow and can require up to two months before results are obtained. A faster, 22 also commonly used tool for drug resistance detection is commercial genotypic line 23 probe assays such as XPert MTB/RIF and GenoType MTBDRplus. However, such 24 assays only predict resistance to isoniazid and/or rifampicin; they also only consider 25 common mutations associated with resistance, and therefore have lower sensitivity. Due 26 to recent advancements in whole-genome sequencing (WGS) technology, researchers are 27 able to sequence the whole genome of bacteria within a few hours at a much lower cost 28 while extracting more information about genotypic resistance markers. Mutation 29 libraries and databases based on single-nucleotide polymorphisms (SNPs) have been 30 constructed to detect drug resistance; these are included in tools such as TBProfiler [3], 31 TBDreaMDB [4], PhyResSE [5] and Mykrobe [6]. However, these libraries are not able 32 to capture the interactions between different resistance markers in predicting resistance; 33 moreover, they are highly dependent on the diversity of the sequenced samples.
Machine learning and deep learning have been successfully applied to different areas 35 of research such as natural language processing, computer vision, and more recently 36 bioinformatics, including GWAS and cancer prognosis prediction [7,8]. Leveraging the 37 abundance of WGS data, researchers have also applied machine learning tools to 38 predicting antibiotic resistance in MTB; for instance, in [9] they applied machine 39 learning models such as logistic regression, support vector machines, and naive Bayes 40 classifiers to predict resistance to 8 drugs from WGS data; in [10], the authors presented 41 a deep learning model, a multi-task wide and deep neural network for drug resistance 42 prediction, where they combined logistic regression and various deep learning techniques. 43 Both of these approaches have demonstrated promising results relative to curated 44 libraries in terms of specificity and sensitivity.

45
Previous work has focused exclusively on the prediction aspect and has not fully 46 explored the insights which could be obtained from the trained machine learning models 47 themselves, namely, a feature importance analysis. Such an analysis is especially 48 important for biological problems, where practitioners such as clinicians or public health 49 professionals may be interested in studying the interpretable factors, such as the genes 50 or SNPs, which are involved in determining resistance. Moreover, previous work has not 51 directly demonstrated the robustness of these models in making predictions on datasets 52 of different diversity.

53
In this work, we investigate the robustness of machine learning models for MTB 54 drug resistance prediction by training them on multiple datasets of different diversity 55 and composition. We demonstrate the pitfalls with regards to the generalizability of 56 models trained on homogeneous datasets with low diversity and regional data.

57
Moreover, we investigate more in depth those genes that are highly associated with 58 determining resistance by using different methods to determine which genes are of 59 particular importance for each model in determining drug resistance.

61
In this work, we focus on three geographically distinct datasets: the Relational 62 Sequencing TB Data Platform (ReSeqTB) [11], the Pathosystems Resource Integration 63 Center (PATRIC) database [12], and a collection of isolates from the British Columbia 64 Centre for Disease Control (BCCDC). ReSeqTB is a consortium dedicated to the 65 development of next-generation TB diagnostics, and the dataset comprises global 66 isolates from 32 different countries. The PATRIC database contains isolates from 67 around the world. The BCCDC dataset is considered to be homogeneous, since it 68 originates from only one location, while the ReSeqTB and PATRIC datsets are 69 heterogeneous. Since all three datasets contain phenotypes for five first-line drugs -  In this work, we predict drug resistance to single drugs and focus on SNPs within 78 two lists of candidate genes, one with 23 genes and one with 179 genes, both of which 79 were previously reported to be associated with drug resistance in MTB [13,14]. The   Figure 3 illustrates the workflow of our machine learning approach.

95
Calling single nucleotide polymorphisms (SNPs) 96 For all three datasets, SNPs were called using Snippy [16], with Mycobacterium 97 tuberculosis H37Rv used as the reference genome. We translated the SNPs and the 98 phenotypic data into a feature matrix X and a list of label vectors y D indexed by the 99 drug D, respectively. X is a binary matrix, where x ij = 0 indicates that sample i does 100 not have the j-th SNP, while x ij = 1 indicates that it does. y D is a binary vector with 0 101 indicating susceptibilty and 1 indicating resistance to the given drug D.  For the cross-regional experiments, we used the fully available ReSeqTB genotype 103 data with calls using the UVP pipeline [17]. This was done to ensure more data was 104 available than just the samples committed to a public repository and for which the 105 corresponding accession number was available. This allowed us to obtain SNPs from a 106 total of 6495 isolates through ReSeqTB.

107
Cross-validation protocol 108 We performed a 10-fold cross-validation to choose hyperparameters. In probability threshold which yielded the highest accuracy on the training set. We 114 measured performance using the area under the Receiver Operating Characteristic curve 115 (AUC). We picked the best parameters for the model based on the average AUC across 116 the 10 folds on the validation set. Once the set of best parameters was found, we 117 retrained the model on the training and validation sets, then predicted on the test set 118 and reported the AUC of the test set for each fold. Machine learning classifiers 120 We considered the main models used in [9], namely, Logistic Regression with L 1 121 and L 2 penalty (LR L1 and LR L2), support vector machine with linear kernel and L 2 122 penalty (SVM L2), support vector machine with radial basis kernel (SVM RBF) and 123 random forests (RF). In our work we additionally investigated the ElasticNet (eNet) 124 model. eNet refers to a logistic regression model which combines L 1 and L 2 penalties. 125 Resistance-associated genes 126 We investigated resistance-associated genes for each drug based on the 127 machine learning models using the following two methods:

128
• Submatrix permutation test: For each gene g in the list of (23 or 179) candidate 129 genes, we randomly permuted the rows of the submatrix of the input feature 130 matrix X whose columns were SNPs within g, while leaving the rest of X 131 unchanged. Denoting this permuted matrix as X g , we then input X g to the 132 machine learning models for retraining. Let AU C(X) and AU C(X g ) denote the 133 test AUC using X and X g as the feature matrix, respectively. We computed 134 ∆ g = AU C(X) − AU C(X g ). If ∆ g is positive and large in magnitude, we may 135 conclude that gene g is relevant to resistance determination for the drug of 136 interest.

137
• Gene-wise partial least square (PLS) regression: PLS is a regression model which 138 decomposes X into a set of components explaining as much as possible of the 139 covariance between X and y [18]. We performed PLS on X gene-by-gene and kept 140 only 1 component for each gene. Hence, the PLS-transformed matrix P LS(X) is 141 an n × m matrix, where n is the number of samples and m is the number of genes 142 considered. We then trained and predicted resistance using the LASSO 143 (LR L1) [19] and the random forest (RF) models with X := P LS(X) as the 144 input matrix. By investigating the resulting models, we are able to estimate the 145 importance of each gene for predicting resistance to a particular drug. deletions, which can be determinant for PZA resistance [20].

170
Any given probability threshold for assigning a prediction (drug-sensitive or 171 drug-resistant) to each sample determines a sensitivity and a specificity for the model. 172 Rather than using a fixed threshold, we take a holistic view and make use of the 173 1 Strains from the same lineage are generally thought to have originated from the same geographical location; they tend to share some unique mutations. does not take into account insertions and deletions, which can be determinant for PZA 186 resistance [20].

187
The performance ranking of the heterogeneous datasets is identical, whereas BCCDC 188 has PZA and RIF as the easiest drugs to predict, and SM as the worst one. Moreover, 189 RIF and PZA are the easiest to predict for BCCDC, with AUC's greater than 0.95,  By training our models on one dataset and testing them on another, we were able to 197 observe how generalizable our models are. We use the short-hand notation a → b to 198 denote a model with training set a and testing set b. Overall, prediction performance is 199 higher when training on a heterogeneous dataset (i.e. ReSeqTB or PATRIC) and   other drugs. There is an increase of 11% in the AUC of PZA while there is an increase 205 of only 2% for RIF and even a decrease of 8% for ReSeqTB→BCCDC compared to its 206 opposite. Nearly identical results are obtained for PATRIC→BCCDC and 207 BCCDC→PATRIC. It is interesting to note that for PZA prediction, training on a 208 heterogeneous dataset is better, while for SM, training on a homogeneous dataset is 209 better. 210 We also found that training with 179 genes instead of 23 can in some cases make the 211 model perform better in cross-dataset evaluation. For example, for ReSeqTB→BCCDC, 212 the AUC for RIF increased from 73% to 83% when training on 179 genes. However, it is 213 also possible for the performance to become worse, which can be seen for SM.
214 Cross-regional evaluation shows the pitfalls of applying a model 215 trained on regional data on a different region 216 We also hypothesized that regional variation can have an effect on the predictive power 217 of models, and thus in this section, we further investigate the regional impact of a can range from 50% to 90%, while RIF prediction performance is stable around 90%. It 225 may be possible that mutations encoding for RIF and INH resistance are homogeneous 226 across strains, meaning there is not much variation in the resistance mechanism across 227 strains, while regional variation may exist for the resistance mechanism of other drugs 228 such as PZA or SM.

229
Trends were also observed when training on different regional subsets. The North 230 American model was much better at predicting resistance for African samples than 231 Asian samples. Similarly, the African model was much better at predicting resistance on 232 North American than Asian samples. The Asian model performed better on the North 233 American than African samples, although the training data was limited. It appears that 234 prediction of Asian samples may be a harder task from both the North American and 235 African samples. It is possible that North American and African strains are 236 genotypically closer to each other than Asian strains. However, this hypothesis would 237 require further investigation and a much more comprehensive dataset, as each regional 238 subset may not be representative of real-world strains present in each region.

239
Furthermore, we did not use European data because the European subset appears to 240 contain many different types of MTB due to the presence of many different groups of 241 immigrants.

242
Which genes are associated with drug resistance?

243
Submatrix permutation tests show that including more features may dilute 244 the importance of key features 245 We considered different machine learning models for estimating gene importance using 246 random permutation of the rows of the feature matrix, to provide an unbiased 247 perspective. We chose the LR L1, LR L2, SVM L2, SVM RBF and RF models for this 248 study and posited that a gene was important for resistance prediction if all the models 249 assigned it a relatively high importance. We define Delta as the largest positive 250 difference in test AUC between all models before and after permuting. The results are 251 shown in Supplementary Table 1. We consider genes that are predicted to be important 252 through machine learning as well as associated with resistance for the corresponding 253 drug in [23] Table 2. Top 5 genes with the highest Delta in the submatrix permutation tests in BCCDC. u indicates upstream regions and * indicates a gene known to be resistance determinant.
Gene-wise partial least squares using the LASSO model reveals the cause  Gene-wise partial least squares using Random Forests shows the difficulty 288 of working with limited sample size heterogeneous datasets 289 We used the same reduced matrices as for the LASSO model, but as an input for RF. A 290 gene's importance was evaluated by investigating the model's feature importance values. 291 The results are presented in Supplementary Table 3.

292
The difference between using 23 and 179 genes is not as evident, but there is still   In summary, we predicted TB drug resistance using different machine learning models 304 on 3 large datasets, ReSeqTB, PATRIC and BCCDC, and performed the study of 305 resistance-associated genes in a machine learning context. We also performed a 306 cross-dataset evaluation and cross-regional evaluation.

307
The first contribution of our work is demonstrating the robustness of models when 308 training and testing on datasets of different diversity. The test AUC for BCCDC is 309 much higher than the ones for ReSeqTB or PATRIC. This demonstrates the ease of 310 training and testing on a homogeneous dataset compared to a heterogeneous one.

311
Furthermore, the performance difference between training on SNPs within 23 genes and 312 SNPs within 179 genes was not significant, so increasing the size of the feature set does 313 not consistently result in better performance when the sample size remains constant.

314
Increasing the number of features may introduce too much noise that results in the 315 model not performing any better. We also see that a model trained on a homogeneous 316 dataset may not generalize well to a heterogeneous dataset. We demonstrate the 317 potential difficulty of achieving good prediction results when a model trained on a 318 homogeneous dataset is used to make predictions on a heterogeneous one. We also 319 demonstrate the potential challenges of training a model using data from one region. 320 Such a model may only be suitable for that region, but may not be generalizable to 321 other regions.

322
Apart from achieving good results on prediction, it is equally important to study the 323 trained models and understand resistance-associated markers, something that seems to 324 be a limitation of most previous work. Our second contribution is a rigorous method for 325 identifying putative resistance-associated genes by permuting the rows of the submatrix 326 of the corresponding gene as well as performing gene-wise PLS with LASSO. The 327 gene-wise PLS method tends to give a better view of which genes are important, while 328 the submatrix permutation testing method tends to be noisier and less consistent.

329
However, our analysis still showed that some of the highly ranked genes were previously 330 reported as resistance-determinant, supporting the use of our methods. Some other 331 identified genes have not been extensively studied in conjunction with drug resistance, 332 and thus may contain true signal or just noise. Therefore, one key challenge is to 333 construct a reliable feature importance extraction method, preferably from a machine 334 learning perspective, but possibly involving expertise-based feature engineering, to 335 eliminate SNPs that may contribute to irrelevant genes being identified as important.

336
By effectively eliminating noise, we are also able to include more relevant SNPs, which 337 will potentially amplify the advantage of training models on a diverse dataset.

338
Our work can be seen as an exploratory study. To the best of our knowledge, it is 339 the first one to investigate the issues of dataset homogeneity and heterogeneity and the 340 identification of drug resistance markers at different resolutions by analyzing machine 341 learning models trained on 3 major datasets and two sets of genes. We hope that this 342 study will inform future work on drug resistance in pathogenic bacteria and the 343 application of machine learning to the drug resistance problem.