{'fraction_non_alphanumeric': 0.030049831378668145, 'fraction_numerical': 0.04263351286052247, 'mean_word_length': 4.871158392434988, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pan-genomic and Polymorphic Driven Prediction of Antibiotic Resistance in Elizabethkingia

The Elizabethkingia are a genetically diverse genus of emerging pathogens that exhibit multidrug resistance to a range of common antibiotics. Two representative species, Elizabethkingia bruuniana and Elizabethkingia meningoseptica, were phenotypically tested to determine minimum inhibitory concentrations for five antibiotics. Ultra-long read sequencing with Oxford Nanopore Technologies and subsequent de novo assembly produced complete, gapless circular genomes for each strain. Alignment based annotation with Prokka identified 5,480 features in E. bruuniana and 5,203 features in E. meningoseptica, where none of these identified genes or gene combinations corresponded to observed phenotypic resistance values. Pan-genomic analysis, performed with an additional 19 Elizabethkingia strains, identified a core-genome size of 2,658,537 bp, 32 uniquely identifiable intrinsic chromosomal antibiotic resistance core-genes and 77 antibiotic resistance pan-genes. Using core-SNPs and pan-genes in combination with six machine learning algorithms, binary classification of clindamycin and vancomycin resistance achieved f1 scores of 0.94 and 0.84 respectively. Performance on the more challenging multiclass problem for fusidic acid, rifampin and ciprofloxacin resulted in f1 scores of 0.70, 0.75 and 0.54 respectively.


Introduction 20
Emerging antimicrobial resistance (AMR) is a global crisis. A recent report has predicted that by 21 2050 antimicrobial resistance will lead to 10 million deaths annually and cost the world's respectively (Gordon et al., 2014). Another study investigating 681 Neisseria gonorrhoeae 39 isolates achieved an acceptable major error rate (93% accuracy within one MIC doubling 40 dilution and 98% for two) by regressing MIC phenotypes on genetic mutations of known AMR 41 genes (Eyre et al., 2017). While the results are comparable to routine antimicrobial testing, they 42 rely on the prior knowledge of single-gene products that relate to observable AMR phenotypes. 43 The wide range of expressed AMR phenotypes in the bacterial domain suggests that the genetic 44 basis for the evolution and transmission of AMR is driven by a complex interplay of several 45 factors. These include the rate at which resistance genes and mutations arise, the level of 46 resistance contributed by the acquired genetic variants, and the relative fitness of the resistant 47 mutants while under selective pressure from the drugs (Sommer et al., 2017). In addition, there 48 are a large number of intrinsic chromosomal genes found in many bacteria, such as the marRAB 49 operon (Vinue et al., 2013), which are involved in complex epistatic interactions regarding 50 antimicrobial resistance expression (Gambino et al., 1993;Martin and Rosner, 2002;Lim et al., 51 2018). This can further render single AMR gene-based interpretations ineffective for predicting 52 resistance phenotypes (Vinue et al., 2013). 53 Machine learning (ML) has led the way in cutting-edge prediction accuracy for vision tasks 54 (Voulodimos et al., 2018), time-series problems (Ahmed et al., 2010), and the clustering of high-55 dimensional data (Assent, 2012). These methods have seen success in genomics for gene-finding 56 (Libbrecht and Noble, 2015), predicting the functional consequences of protein missense 57 mutations (Shihab et al., 2013), and genetic structure discovery using Markov clustering 58 (Kopelman et al., 2015). ML is capable of dealing with high-dimensional interactions and 59 nonlinear relationships in data, and has shown promise in using SNPs for predicting phenotypes 60 that have a complex genetic architecture (Ban et al., 2010) and using k-mer counts (Nguyen et  other non-Elizabethkingia isolates were included to construct core-SNPs and pan-gene 80 presence/absence matrices, with which AMR predictability was evaluated for several ML 81 algorithms. In this report, we detail the methods used to produce efficient AMR prediction for 82 both binary and multiclass resistance profiles.  Antibiotic resistance profile evaluation 92 Bacterial isolates were grown and maintained as described previously (Johnson et al., 2018).   Per-read quality filtering consisted of a multi-step procedure to maximize read length and read 143 quality for assembly. The reads, in FASTQ format, were subjected to a quality filter pass with a 144 minimum Phred score of 12 using PRINSEQ (Schmieder and Edwards, 2011) with the -145 min_qual_mean switch. Reads with a length of 1,000 bp or lower were also discarded with 146 PRINSEQ's -min_len switch.

147
De novo assembly of each organism's reads was completed with Canu v1.5 , 148 a Celera Assembler successor designed to generate high-quality assemblies from Nanopore or 149 PacBio long-reads. Canu was chosen because it provides higher assembly sequence identity than  After producing the initial de novo assembly with Canu, Nanopolish v7.1 (Loman et al., 2015) 155 was used to improve the overall assembly quality for each sequence using a hidden Markov 156 model. All original base-called, signal-level reads were re-extracted to tag the reads with 157 identifying information for Nanopolish; these tagged reads were then aligned to their respective 158 assemblies using BWA-MEM 0.7.15 (Li, 2013) using the -x ont2d switch. This command switch 159 reduces the initial seed lengths and uses a relaxed scoring matrix, which allows the effective 160 mapping of ONT's noisy reads to the reference assemblies without producing large-scale 161 fragmentation. After alignment, the produced SAM file was converted into the corresponding 162 binary format (BAM file) using samtools (Li et al., 2009) using the view command and the -sB 163 switch. Nanopolish was run in parallel to produce the consensus sequence for each assembly.       Finally, a third approach using the precise HMMER3 (Eddy, 1998)  AMR. Identification of AMR gene clusters was done by filtering the output for regions with four 199 or more AMR genes that had at most three non-AMR genes between any given gene pair. These  For each antibiotic studied (vancomycin, clindamycin, fusidic acid, ciprofloxacin, and rifampin),  (Table 2). Separately, statistics were generated for the core-genome of Elizabethkingia strains 224 only.

225
To generate a core genome for each group, a multiple sequence alignment of the assemblies for 226 the strains within that group was first completed. The progressiveMauve algorithm (Darling et 227 al., 2010) in Mauve was used to create six different alignments of the assigned bacterial groups.  245 To prepare the input data for downstream predictive algorithms, a predictor X matrix was 246 constructed by directly loading the data from the Mauve SNPs file. Genotypic information from 247 SNPs was represented as the encoded additive value for all polymorphisms on that site. The 248 smallest value (starting at zero) represented the major allele. For the minor alleles, the encoded 249 value increases as the frequency of the allele at that polymorphic site decreases. Effectively, the 250 first minor allele will be encoded as 1, and in the case of multi-allelic polymorphisms, the next 251 most frequency allele will be encoded as 2. The matrix was then transposed to adhere to the 252 traditional structure of an input X matrix, in which all rows represent individuals and columns 253 contain additively encoded SNPs.

254
While SNPs often function as genetic predictors of phenotypic traits, models utilizing only core-255 SNPs do not consider the presence of genes that are out of the core genome. To include these 256 genes in the models, a list for all genes that exist within the strains of each group (pan-genes) 257 was generated. A second, separate X predictor matrix was constructed using a presence/absence 258 value. Each column represents the presence or absence of one putative pan-gene with a unique 259 UniProt ID. A value of 1 is given for strains that contain the gene, while a value of 0 is given 260 when the gene is absent. We consider this gene-centric approach appropriate in cases when  Each set of predictors (pan-genes and core-SNPs) was used in evaluating AMR prediction. A 265 final hybrid method was also evaluated by appending both matrices (SNPs and genes) together to 266 form a combined, third predictor X matrix.

267
In order to assess the predictability of AMR with gene/SNP predictors, strains were labeled 268 either "resistant" or "susceptible" to a particular antibiotic. Members of the Elizabethkingia 269 genus have no standardized breakpoints so other species breakpoints are used for reference. the corresponding row in the X matrix (X i ). In the Y vector, a value of 0 represents a susceptible 281 strain, while a value of 1 represents a strain that is labeled as resistant, with each respective 282 group being assigned its own Y vector for that AMR category.

283
Higher resolution AMR prediction is possible by training models to predict a multiclass 284 "resistance level" instead of a binary resistant/susceptible label. This was accomplished by 285 assigning each strain a resistance level based on their MIC score. Resistance levels are therefore 286 represented categorically as a sliding scale of AMR for the purposes of classification. Raw MIC 287 results for each phenotype were collapsed into these resistance levels (see Supplementary Table   288 3 for resistance level assignments). The selection of ranges for binning MIC values for each 289 phenotype was determined to maximize the uniformity of the categorical phenotypic distribution   Decision tree and random forest algorithms 314 The decision tree is a non-parametric algorithm that can be used for classification or regression. 315 Unlike the Naïve Bayes model, decision trees allow modeling of variable interactions and 316 perform well on samples that cannot be linearly separated. In this tree structure, each node splits 317 up samples based on a determined variable rule. This can be a threshold for continuous variables 318 or a categorical value for discrete and categorical variables. To effectively split samples based on 319 variables, the Gini impurity metric measures how well a particular node split in the tree will 320 separate samples based on output category. Here, this calculation is given by where t is the SNP or gene predictor and tree (Breiman, 2001). This is done through an ensemble-based learning approach, similar to 336 bootstrap aggregation, also known as "bagging" (Dietterich, 2000). Bagging follows the  category for that sample.

417
A key hyper-parameter for SVMs is the constant C, which acts within the soft margin cost 418 function; this C term controls the tightness of the two margins used to produce the hyperplane.

419
Larger C's will produce tighter margins, resulting in less misclassified training samples. A 420 smaller C will produce larger margins, allowing for the misclassification of some training  3, and yield rapidly diminished to 3 mb per-hour by hour 20.

470
Following the results of read filtering, Figure 1

494
The assemblies of the two newly assembled Elizabethkingia genomes have been deposited in a 495 NCBI BioProject (for early access, contact charles.chen@okstate.edu).   The vancomycin group contained a total of 11,066 SNP predictors. With the gene-centric 549 approach, the pan-genome consisted of 4,865 genes.

550
Binary classification of vancomycin resistance revealed that genes are a superior predictor for 551 this particular task. With the exception of the decision tree, Naïve Bayes, and AdaBoost, the 552 gene-centric approach produced consistently superior f1 micro-scores with every algorithm. The 553 highest scoring algorithm, the support vector machine, was able to achieve a mean 0.84 f1 micro-554 score (standard deviation 0.152) with gene predictors (Figure 3). Both the linear SVM and a 555 radial basis kernel SVM performed similarly with large C parameters. In particular, the radial 556 SVM dramatically improved in classification accuracy when increasing C from 0.1 to 10000 557 (Figure 4).

558
Outside of support vector methods, k-NN with k=1 performed with similar accuracy to the linear 559 SVM, in the binary case, being very sensitive to the neighbor parameter: as shown in Figure 4, k-560 NN performed best when k=1, with subsequent degradation in accuracy. Conversely, greater 561 consistency can be achieved with larger k values. For example, k=9 and k=11 reduce the total f1 562 micro-score compared to k=1, but have a smaller standard deviation when computing their mean 563 f1 micro-score over many iterations (when k=1 std.=0.15 and when k=11 std.=0.11 for SNPs).

564
With SNP predictors, accurate binary AMR classification with a decision tree produced 565 comparable results to gene-based SVMs and k-NN and required no hyper-parameter tuning. This 566 method performed the best of all SNP predictor based methods.

567
With genes, ensemble based methods are less accurate than SVMs and neighbor methods in 568 binary classification. Ensemble techniques did not show significant improvements when 569 increasing the number of trees; although, adding more trees improved f1-scores slightly ( Figure   570 4).

571
In the case of the multiclass classification, the core-SNP approach with a radial SVM performed 572 the best, with an f1-score of 0.71 (Figure 3). k-NN (k=1) with genes performed roughly the 573 same. Multiclass classification with tree based methods (random forest and decision tree) 574 showed a 0.06 -0.07 improvement to the f1 micro-score by using SNP variables.

576
The clindamycin group of 28 individuals contained an initial total of 1,996 SNP predictors. With 577 the gene-centric approach, the pan-genome consisted of 6,949 genes.

578
Prediction of binary AMR classification with respect to clindamycin was most effective using 579 SNPs ( Figure 3). An f1 micro-score of 0.94 was achieved using a linear SVM and SNP 580 predictors, compared to a lower micro-score of 0.89 using gene predictors. The radial SVM 581 performed slightly worse than the linear variant and required a carefully chosen C parameter. As 582 seen with other AMR classifications, larger C values were necessary for the radial SVM to 583 perform adequately (Figure 4). k-NN with SNPs produced an f1 micro-score of 0.92 with k=3. With gene predictors, this score 585 decreased dramatically (0.85 score, k=1 was best with gene predictors). When using SNP 586 predictors, any k > 3 marked a consistent degradation in algorithm prediction performance. A 587 similar trend was seen in the variation of the f1-scores during cross-validation; k=3 yielded a 588 lower standard deviation (std. = 0.10) when compared to most other k values.

589
The random forest, with either genes or SNPs as predictors, performed competitively with SVM 590 methods and neighbor methods in terms of f1 micro-score. The other ensemble method,

591
AdaBoost, also generated similar performance between SNPs and genes as predictors but is 592 significantly outperformed by random forests. Changing the number of estimators (trees) in 593 either ensemble algorithm did not affect accuracy; although, having less than 20 trees in a 594 random forest negatively impacted performance.

595
Using SNPs, and a large alpha value, Naïve Bayes achieved a 0.93 f1 score for clindamycin 596 AMR prediction ( Figure 3B). A large discrepancy between the gene predictor and SNP predictor 597 methods was also found in Naïve Bayes AMR classification accuracy where using gene 598 predictors resulted in a ~22% reduction in performance.

599
In the multiclass AMR classification, SVMs (with SNPs) again outperformed all other 600 algorithms with a 0.71 f1 micro-score (for both kernels). Using the other described algorithms With a total of 29 individuals in the rifampin group, a total of 2,044 core SNPs were generated.

616
The pan-genome consisted of 7,805 unique genes.

617
Multiclass prediction of rifampin resistance levels produced the most successful results among 618 the multiclass tests, with a mean f1 micro-score 0.75 (0.173) with a decision tree and SNPs. revealed that support vector machines, k-NN, and random forests can be trained to high accuracy 641 with less than 35 samples, using thousands or tens of thousands of predictors (Table 2). Using  data, with respect to its usage in algorithms, must be applied for appropriate predictive modeling.

689
Prediction accuracy is data-dependent and relies on the quality of the predictor variables used as 690 input. Also, some algorithms, like SVMs and k-NN, require careful hyper-parameter tuning to 691 function properly (Figure 4). In this study, we tested different biological predictor types and 692 found that, depending on the AMR phenotype, either genes or SNPs reliably produce better 693 results. For example, the multiclass fusidic acid classification greatly favors genes over SNPs, 694 while the vancomycin multiclass classification performs better using SNP predictors. Unlike the 695 core-SNPs, our gene-centric model can produce predictors from genomic regions outside of the 696 core-genome, and gene presence/absence predictors do not suffer from small sequencing errors 697 due to their reliance on alignments, which can tolerate some errors using score matrices. This is 698 also advantageous when the strain exhibits an AMR phenotype that is based around a gene 699 product, such as fusB-type fusidic acid resistance. In these genic AMR cases, the effect of 700 mutational changes in the core genome might be too insignificant to be reliably explanatory for 701 AMR phenotypes.

702
A common problem seen when analyzing biological data from individuals that are not part of a 703 controlled population is the relative lack of available samples compared to the wealth of genomic 704 data. The curse of dimensionality characterizes this problem as requiring a tremendous amount 705 of sample data to guarantee that each potential combination of SNP/gene predictor exists within 706 the dataset. As shown in Table 2, the number of variables can far outnumber the sample size and,   707 the ratio of samples to SNPs can be as low as 0.002 (Table 2), depending on the core genome.

728
Together with this research, we suggest that this sample size dilemma can be significantly 729 lessened by capitalizing on the wealth of information stored on cloud services. There are 730 currently more than 150,000 prokaryotic genome assemblies available on the NCBI. We used 731 this community-driven "cloud knowledge" to increase our sample population by 57% and 33% 732 for the vancomycin and clindamycin groups respectively (Table 2), compared to just using in-

733
house Elizabethkingia strains and also produced a more uniform phenotypic distribution which 734 made prediction more feasible.

735
Decreasing sequencing costs and the subsequent increase in availability of genomic information 736 stored on cloud services like NCBI and the European Nucleotide Archive means that "cloud 737 knowledge" will be an effective means of improving sample sets for AMR prediction, provided 738 proper phenotyping is performed. Practical use of a computational AMR prediction pipeline in a 739 clinical or hospital setting will also depend on a computationally efficient predictor generation 740 method. Also, given the large, diverse collection on the NCBI, constructing a core-genome could 741 only be advantageous for species/strains that share common ancestry; the same approach can be prediction. To conclude, we expect ML predictive pipelines, in combination with metagenomics, 754 will shift the paradigm from phenotype-based diagnostics to data-driven prediction for AMR 755 detection and outbreak prevention.

758
The authors declare that there is no conflict of interest.