ABSTRACT
In principle, whole-genome sequencing (WGS) can predict phenotypic resistance directly from a genotype, replacing laboratory-based tests. However, the contribution of different bioinformatics methods to genotype-phenotype discrepancies has not been systematically explored to date. We compared three WGS-based bioinformatics methods (Genefinder [read based], Mykrobe [de Bruijn graph based], and Typewriter [BLAST based]) for predicting the presence/absence of 83 different resistance determinants and virulence genes and overall antimicrobial susceptibility in 1,379 Staphylococcus aureus isolates previously characterized by standard laboratory methods (disc diffusion, broth and/or agar dilution, and PCR). In total, 99.5% (113,830/114,457) of individual resistance-determinant/virulence gene predictions were identical between all three methods, with only 627 (0.5%) discordant predictions, demonstrating high overall agreement (Fleiss' kappa = 0.98, P < 0.0001). Discrepancies when identified were in only one of the three methods for all genes except the cassette recombinase, ccrC(b). The genotypic antimicrobial susceptibility prediction matched the laboratory phenotype in 98.3% (14,224/14,464) of cases (2,720 [18.8%] resistant, 11,504 [79.5%] susceptible). There was greater disagreement between the laboratory phenotypes and the combined genotypic predictions (97 [0.7%] phenotypically susceptible, but all bioinformatic methods reported resistance; 89 [0.6%] phenotypically resistant, but all bioinformatics methods reported susceptible) than within the three bioinformatics methods (54 [0.4%] cases, 16 phenotypically resistant, 38 phenotypically susceptible). However, in 36/54 (67%) cases, the consensus genotype matched the laboratory phenotype. In this study, the choice between these three specific bioinformatic methods to identify resistance determinants or other genes in S. aureus did not prove critical, with all demonstrating high concordance with each other and phenotypic/molecular methods. However, each has some limitations; therefore, consensus methods provide some assurance.
INTRODUCTION
Staphylococcus aureus causes both superficial infections (such as boils) and life-threatening disease, including septicemia (1). There were 11,405 S. aureus bacteremias in England in 2015 and 2016 (2); 7.2% were methicillin-resistant S. aureus (MRSA), which has increased costs and poorer patient outcomes (3). Fast accurate resistance prediction is key to managing S. aureus infections. Molecular-based methods directed at detecting specific genes, e.g., through rapid multiplex PCR and microarrays, can reduce the time to identify resistance determinants and the time on broad-spectrum antibiotics (4–6). However, they require specific primers that impact sensitivity and specificity.
In principle, whole-genome sequencing (WGS) has the potential to predict phenotypic resistance directly from the genotype, replacing laboratory-based phenotypic tests (7). Several studies report high concordance between genotypic predictions based on known or novel resistant determinants and phenotypic methods (8–13). However, these studies used various sequence-processing pipelines and bioinformatics methods to identify in silico resistance determinants. Without formal comparisons between the various methods, it is unclear whether the underlying differences affect results or whether differences in methodology could cause some of the observed discrepancies between genotypic predictions and phenotype.
Therefore, we compared three WGS-based bioinformatics methods (Genefinder [read based], Mykrobe [de Bruijn graph based], and Typewriter [BLAST based]) in terms of predictions of the presence/absence of different resistance determinants and the overall prediction of antimicrobial susceptibility and the presence/absence of virulence genes from short-read Illumina WGS.
MATERIALS AND METHODS
Three sets of S. aureus isolates with known high-quality phenotypes were analyzed: derivation, n = 501, and validation, n = 491, sets (denoted Oxford derivation/validation) from blood cultures and nasal swab isolates at the Oxford Radcliffe Hospitals NHS Trust and Brighton and Sussex University Hospitals NHS Trust, spanning a period of 13 years, sequenced for an initial assessment of genotypic prediction of susceptibility phenotype in S. aureus (9, 10) and 397 isolates that had been referred to the Public Health England (PHE) reference laboratory for investigation (denoted Colindale 397; NCBI BioProject number PRJNA445516). The Oxford derivation set had previously been used in the development of Typewriter and Mykrobe but not Genefinder; the former methods were then applied to the Oxford validation set.
The phenotypes for Oxford derivation/validation isolates were determined by disc diffusion and/or automated broth diffusion (BD Phoenix), with discrepancies between phenotype and genotype resolved as described previously (11). All PHE isolates (n = 397) were subjected to MIC testing by the PHE Staphylococcal Reference Laboratory using the agar dilution method (14). In addition, the mecA-mecC status and virulence gene profile of the PHE isolates was determined by PCR or microarray testing as described previously (15). The European Committee on Antimicrobial Susceptibility Testing (EUCAST) thresholds were used to determine the sensitivity or resistance for each phenotype (http://www.eucast.org/clinical_breakpoints).
All Oxford derivation/validation isolates were sequenced using the Illumina HiSeq 2000 platform as previously described (16). PHE samples were sequenced in an Illumina HiSeq 2500 platform as described previously (17) (both 150-bp reads). Samples determined as mixed based on WGS were excluded from further analysis. For quality control of sequences at PHE, trimmomatic software was used (Illumina adapter removed, leading and trailing quality threshold set to 30, and minimum length of read set to 50 bases) (18). Isolates from Oxford analyzed by Typewriter were mapped and de novo assembled with exclusion parameters of <70% coverage of the reference genome for mapping and <50% of the genome in contigs >1 kb (10). Mykrobe processes raw sequence data with no prior cleaning of the data. The isolates came from 111 sequence types, including 29 new sequence types (STs)/alleles, covering the range of S. aureus genomic diversity as previously described in Oxfordshire.
Three programs, Genefinder (M. Doumith, PHE, not published), Mykrobe (P. Bradly, version v0.3.13-2-gd5880fa, open-source at https://github.com/iqbal-lab/Mykrobe-predictor), and Typewriter (T. Golubchik, MMM group, Oxford University; version 2.0, https://github.com/tgolubch/typewriter) (Table 1), were compared to determine presence/absence of resistance determinants (genes or variants) and toxin genes (Tables 2, 3, and 4). Mykrobe is part of the automated processing with the Complete Pathogen Software Solution (COMPASS) developed at the University of Oxford. This returns quality and depth of sequence metrics, maps against a reference (MRSA 252, GenBank accession no. BX571856.1) using Stampy (19) and performs de novo assembly using Velvet v1.0.18 (20). These de novo assemblies formed the basis for the Typewriter program, whereas Genefinder used the raw sequencing reads.
Overview of Genefinder, Mykrobe, and Typewriter methods and requirements
Predicted antibiotic susceptibility phenotypes from WGS by Genefinder, Mykrobe, and Typewriter
Predicted phenotype for antimicrobial susceptibility
Predicted genotype for virulence genes, ccr genes, and mecA-mecC
Although all three methods search for matches to a predefined list of alleles, they have different approaches to their identification (further details below). Genefinder and Mykrobe required fastq files, whereas Typewriter used BLAST use de novo assemblies. All used preset thresholds to detect genes. Thresholds are adapted for certain genes (e.g., blaZ, which can be chromosomally integrated or carried on plasmids) to improve the prediction and for quality control. Both Typewriter and Mykrobe identified the presence or absence of each target singly, whereas Genefinder identified which of closely related homologs is most plausibly present. Genefinder and Mykrobe were very fast, between 1 and 3 min, and can be used on a standard desktop computer (specification of a 2.3-GHz processor and 16-GB memory). Typewriter, as it requires de novo assembly, took up to 3 h and used cloud computing or high-capacity servers.
Genefinder was written by M. Doumith. It used a mapping approach (similar to SRST2, https://github.com/katholt/srst2) to detect the presence or absence of predefined genes or variations in predefined genes using Bowtie. The thresholds were defined at 90% overall, but amended where required in order to distinguish between both variants where genes were represented with multiple reference sequences and the level of diversity expected for each gene sought. Genefinder also checked for premature stop codons and compared the average depth of read coverage to identify any potential sequence contamination.
Mykrobe was written by P. Bradley and Z. Iqbal (9). A threshold frequency was generated for each gene (K minimum percentage) based on the empirical level of diversity observed in the training set described by Bradley (K = 0.3 for blaZ, K = 0.6 for fusB and fusC, K = 0.8 otherwise). The maximum likelihood from 3 models (gene absent, gene present in minor proportion, gene present) was chosen. The models took into account the expected proportion of kmers based on the depth of coverage and empirical level of diversity (described in reference 9). Mutations were genotyped by choosing the maximum likelihood model from 3 Poisson models comparing the depths of coverage across 63-bp reference and alternate alleles while demanding 100% coverage across the allele, also described in reference 9.
Typewriter was developed by T. Golubchik (described in reference 10). It considered BLAST results over a query reference (blastn for sequence identity, tblastn for mutations). It used a “relative coverage” to determine presence/absence of a gene, a metric that gives equal weights to coverage and sequence identity. Typewriter reported this value for each query gene of interest, and cutoffs were adjusted to optimize the specificity/sensitivity for different genes. In this study, a relative cutoff of 90% for resistance and toxin genes was used except blaZ, for which a cutoff of 80% was used. For variant reporting, mutations were reported above a given threshold of relative coverage (e.g., 90%); however, this could be changed or set to 0% to report all identified differences from the query sequence. Stop codons were predicted, as were novel mutations.
Eighty-four genes were included in the analysis; 46 acquired resistance genes, 5 sets of chromosomal variants within resistance-associated genes, 5 cassette chromosome recombinase genes (ccr), and 28 virulence genes (Tables 2, 3, and 4). Acquired resistance genes were classified as present (p, P) or absent (a, A), setting 3 missing Genefinder predictions (“ND” or “X”) to absent. Chromosomal resistance variants were those listed in Table S4 in the supplemental material; 23 other mutations were reported in the relevant genes but were not compared, as they are not considered resistance-determinants (Table S4 in the supplemental material). For all methods, genotype predictions of susceptibility phenotypes were based on the presence of any relevant resistance determinant as shown in Tables 2, 3, and 4 (as described in reference 10 with minor modifications and updates from reference 9). Intermediate phenotype results were excluded from analysis (80 cases [0.5%]).
RESULTS
Short-read Illumina WGS data were available from 1,389 samples; 992 from a collection held in Oxford (previously described by Gordon et al. [9, 10]) and 397 from Public Health England (PHE) Staphylococcus Reference Service, Colindale. Ten samples were excluded due to mixed/contaminated WGS results, leaving 1,379 for analysis. Samples were analyzed by Genefinder and Typewriter (Table 1) after sequence mapping and variant calling and by Mykrobe from raw fastq reads.
Eighty-four genes were included: 46 acquired resistance genes, 5 sets of chromosomal variants within genes associated with resistance, 3 cassette chromosome recombinase, ccrA, ccrB, and ccrC, including three variants of ccrC [ccrC(a), ccrC(b), and ccrC(c)], and 28 virulence genes (see Table S1 in the supplemental material). Overall, 99.5% (113,830/114,457) of the individual resistance determinant/virulence gene predictions were identical between all three methods (Fig. 1; Table S1), with only 627 (0.5%) discordant predictions, demonstrating high overall agreement (Fleiss' kappa = 0.98, P < 0.0001). Overall, one method disagreed with both other methods in 0.23% for Typewriter (263/114,457 predictions), 0.16% for Mykrobe (183/114,457), and 0.16% for Genefinder (181/114,457). The three most common discrepancies for Typewriter were the nondetection of virulence genes identified by other methods (seu, 57 samples; chp, 46 samples; sei, 33 samples). Similarly, for Genefinder, the three most common discrepancies were nondetection of resistance genes (qacB, 44 samples; dfrC, 34 samples) or other genes (ccBb, 22 samples) identified by other methods. Genefinder reported the presence of dfrA, qacA, or ccrC(b) genes in these samples. In contrast, Typewriter and Mykrobe reported the presence of two dfr, two qac, and three ccrC genes, where the detected variants for each of these three genes shared more than 90% nucleotide identity. The most common discrepancies for Mykrobe were identifying resistance/other genes as present when the other two methods called them absent (aadE–ant(6)-Ia, 28 samples; blaZ, 19 samples; ccrCB, 22 samples). No gene was ever identified as present by Typewriter alone. Fourteen of the 84 genes had >1% discrepancies (maximum, 4.3% for seu), but the majority of discrepancies were in only one method for all genes except ccrC(b).
Determinant-by-determinant disagreements between methods. Each panel shows percentage differences in proportions of the detected presence of each determinant between the first method and the second.
Discrepancies were similar in acquired resistance genes (0.3% [221/63,434]) and chromosomal resistance genes (0.1% [8/5,516]), but slightly larger for ccr genes (1.8% [123/6,895]) and virulence genes (0.7% [275/38,612]) (see Table S2). The percentage discrepancies varied modestly across the different sample sets, being higher for the PHE set (1.1% [349/32,928]; particularly for ccr genes with 4.2% [83/1,960] discrepancies), intermediate for the Oxford derivation set (0.6% [233/42,084]) and lowest for the Oxford validation set (0.1% [45/40,824]) (Table S2).
Genotypic predictions of antimicrobial susceptibility were also identical in 99.6% of cases (16,477/16,548 predictions) (Table 2). Of the 71 discrepancies in susceptibility prediction between the methods, 42% (30/71) occurred with Typewriter reporting susceptible when Genefinder and Mykrobe reported resistant, and 49% (35/71) occurred with Mykrobe reporting resistant where Genefinder and Typewriter reported susceptible.
Comparing genetic predictions to laboratory phenotypes (restricted to samples either phenotypically resistant or susceptible), in 98.3% (14,224/14,464) of cases, all three bioinformatics methods and the gold standard laboratory results agreed completely (2720 [18.8%] resistant, 11,504 [79.5%] susceptible) (Table 3 and Fig. 2). There was greater disagreement between the laboratory phenotypic results and the combined genotypic predictions than within the three bioinformatics methods. In 97 (0.7%) instances, the laboratory phenotype was susceptible, but all bioinformatic methods reported resistance. Of these, 33% (32/97) were for penicillin, 23% (22/97) for clindamycin, and 11% (11/97) for erythromycin, with smaller numbers for fusidic acid (7), tetracycline (6), mupirocin (6), methicillin (5), ciprofloxacin (4), gentamicin (3), and rifampin (1), and none for trimethoprim. In 89 (0.6%) instances, the laboratory phenotype was resistant, but all three bioinformatics methods reported susceptible, most commonly to gentamicin (21% [15/89]), ciprofloxacin (17% [15/89]), and fusidic acid (15% [13/89]). The remaining 54 (0.4%) cases (16 phenotypically resistant, 38 phenotypically susceptible) had different genotypic predictions made from the different methods. However, in 36/54 (67%), the consensus genotype (predicted by two of the three methods) matched the laboratory phenotype.
Antimicrobial susceptibility genotypic predictions compared to phenotype.
PCR/array results were available for some virulence genes (15) and mecA-mecC for all 397 PHE isolates. Compared with genetic predictions, in 96.8% (3,983/4,115) of cases, all three bioinformatics methods and the PCR/array results agreed completely (3,364 [81.7%] absent, 619 [15.0%] present) (Table 4; see also Fig. S1). As for antimicrobial resistance, there was greater disagreement between the laboratory PCR/array results and the combined genotypic predictions than within the three bioinformatics methods, with 81 (2.0%) cases where all three methods called a gene present that had not been detected by PCR/array and 12 (0.3%) where no method called a gene present that had been detected by PCR/array, in comparison with 39 (0.9%) discrepant predictions between the methods. In 20/39 (51%) cases, the consensus genotype matched the PCR/array result.
The sensitivity and specificity of all three bioinformatics methods compared to laboratory phenotypic methods in predicting antimicrobial susceptibility were very similar. Across the 14,464 genotypic predictions, Typewriter had the lowest overall sensitivity (0.964 [95% CI, 0.956 to 0.970]), but the highest specificity (0.992 [0.990 to 0.993]), while Mykrobe had higher sensitivity (0.967 [0.960 to 0.974]) and the lowest specificity (0.989 [0.987 to 0.990]). Genefinder's performance fell between that of Mykrobe and Typewriter for specificity (0.990 [0.988 to 0.992]), with a sensitivity equal to that of Mykrobe (0.967 [0.960 to 0.973]). Specificity and sensitivity varied across the different antibiotics (Fig. 3), but were broadly similar between the three methods, overall and within each data set (see Table S3). There were no vancomycin-resistant isolates identified by either phenotyping or bioinformatics methods. Similarly, specificity and sensitivity to identify PCR-detected virulence and other genes varied across the different genes, but were broadly similar between the three methods (see Fig. S2).
Sensitivity and specificity of genotypic predictions of antimicrobial susceptibility.
DISCUSSION
While WGS is increasingly used to detect antibiotic resistance and virulence determinants, to our knowledge, this is the first study that compares three methods for predicting genotypes of large numbers of isolates. As discussed in the recent European Committee on Antimicrobial Susceptibility Testing (EUCAST) report (21), discordance can occur between phenotypic and genotypic resistance due to inadequate limits of detection for WGS methods, incomplete understanding of the genotypic basis of phenotypic resistance, flaws with the phenotypic or molecular (e.g., PCR) methods currently used to detect resistance, and/or WGS failures, including the lack of assembly caused by multiple operons or similar sequences, incomplete gene coverage, nonfunctional genes (e.g., due to the presence of stop codons/indels), or cropped contigs.
Here, we found that three different approaches for identifying genetic determinants of resistance and virulence (Genefinder, Mykrobe, and Typewriter) agreed in 99.5% of predictions. Genefinder and Mykrobe were fast, taking under 5 min, whereas Typewriter, while also taking a few minutes per sample, required an initial genome assembly that increased the turnaround time by up to 3 h. Mykrobe and Typewriter are freely available (https://github.com/iqbal-lab/Mykrobe-predictor and https://github.com/tgolubch/typewriter, respectively); Genefinder is not, but the underpinning methods are relatively straightforward, and the freely available SRST2 (https://github.com/katholt/srst2) follows an analogous mapping approach (22), which would likely provide very similar results with the same catalogue. Previous comparisons of bioinformatics methods relevant to the microbiology community are limited. Bradley et al. (9) found good concordance between Mykrobe and SeqSphere (23), an allele-based method that detects the presence/absence of a limited number of resistance and virulence markers. SeqSphere took longer than Mykrobe as, like Typewriter, it uses Velvet assemblies. Other previous studies have shown 100% concordance between the resistome and toxome in 14 MRSA isolates (24), 98.6% concordance across 5,288 susceptibility predictions in 308 S. aureus isolates (both MRSA and MSSA) (25), 100% concordance for selected resistance and toxin gene presence/absence in 18 MRSA strains (23), and 97% and 97% sensitivity and specificity, respectively, for Typewriter and 99.1% and 99.6% sensitivity and specificity, respectively, for Mykrobe for predicting phenotypic resistance in the Oxford validation samples used here (9, 10). A comparison between microarray and WGS in 154 isolates reported 1.7% discordancy in detecting resistance and virulence genes (26), mainly due to the failure of WGS to detect enterotoxins and super antigens (similar to that for Typewriter in this study).
Individually, the three programs demonstrated high concordance, but interestingly, in almost all genes, only one of the three bioinformatics methods did not identify a determinant that the other two methods did identify, or vice versa. The most common discrepancy with Typewriter was the failure to identify virulence genes identified by Mykrobe and Genefinder (namely, seu, chp, and sei). Two of these genes, sei and seu, are located on the enterotoxin gene cluster (egc) (27, 28), referred to as an enterotoxin gene nursery (29), and the other, chp, is located on a prophage (30). Such regions may be particularly susceptible to recombination (31, 32) and paralogs. As Typewriter uses BLAST, it may have a higher chance of detecting one of multiple closely related genes than the other two methods.
Similarly to Typewriter, the most common discrepancy with Genefinder was a failure to identify genes reported by Typewriter or Mykrobe, particularly, ccrB, qacB (quaternary ammonium compound B, conferring resistance to chlorhexidine [33] via an efflux drug pump, but differing from another gene, qacA, by only seven nucleotides [34]), and dfrC (a dihydrofolate conferring resistance to trimethoprim believed to be the origin of the more common transposon-associated drfA gene). The fact that Genefinder identified only one variant of acquired dfr and qac may indicate that the other two methods were misidentifying paralogs (35). Alternatively, as Genefinder detects predetermined alleles, the recombination of partial genes or differences in flanking sites or genomic variation alone may reduce its ability to detect some genes. One advantage of Genefinder is its ability to detect variations in multicopy genes such as the rRNA encoding genes associated with linezolid resistance in staphylococci.
In contrast, Mykrobe most commonly identified a determinant that other methods did not, particularly, aadE(ant6′)-Ia, an adenyltransferase encoding resistance to aminoglycosides. This gene is associated with small plasmids flanked by direct repeats of staphylococcal insertion sequence IS257 (36). Although Mykrobe is kmer based, it requires a high match across the whole gene, not just flanking sequences, and so the reason for this is unclear. Mykrobe also had a higher false-positive rate in blaZ, as reported previously (9). Although this was previously attributed to phenotypic errors, the fact that neither Genefinder nor Typewriter identified blaZ in these isolates suggests the algorithm/threshold may need adjusting for this gene. Mykrobe also had a high false-positive rate for the ccrCB gene, which is part of the cassette chromosome recombinase (ccr) associated with SSCmec (37). As all ccrC genes share >87% similarity and were not included in the original Mykrobe implementation, further investigation and modification of sequence identity thresholds may be required to accurately classify this gene, whose different alleles can have 60% to 82% sequence identity.
Overall, the comparison highlights key challenges inherent in all methods. First are the trade-off between specificity and sensitivity to detect specific genes/variants and the need for adjustment based on specific features, such as the proximity to repetitive elements or similarity with other alleles. Specific genes may also require different approaches, e.g., the ccr genes were the most discordant overall in the study. These genes were more often present in the staphylococcal reference laboratory isolates, increasing the overall error rates for this sample set. Reference libraries of genes/variants also require frequent updating with new alleles, and appropriate thresholds must be set to enable separate copies of closely related genes (e.g., qacA and qacB) to be detected if genuinely present. Taking the consensus prediction across the three different bioinformatics methods is one strategy for balancing these different trade-offs. As error rates were low overall, this only improved the genetic predictions slightly, but in samples where the susceptibility is unknown, it could be valuable, particularly if the two fast implementations (Genefinder and Mykrobe) are used, followed by the slower assembly-based method only if they disagree.
Our main findings were that the largest discordance occurred between phenotype and genotype regardless of the method used to predict genotype and that the “consensus” genotypic prediction agreed with the phenotype in two-thirds of the small number of cases where bioinformatics methods made different predictions. Where bioinformatics methods are concordant but disagree with phenotype, the unresolved question is which is “correct,” in terms of a drug achieving clinical cure in a patient infected with this strain. Penicillin and clindamycin/erythromycin were most likely to be called resistant by all methods but susceptible by phenotyping. Previous studies of erythromycin and clindamycin resistance have reported positive ermC PCR results from nondetectable resistance phenotypes (38) and have suggested that plasmids conferring resistance to these antibiotics may be lost in subculture (9, 39). Sensitivity to penicillin by phenotypic methods where genotype methods predict resistance has been reported previously (40, 41), and the evidence suggests that phenotyping underreports resistance. The EUCAST guidelines illustrate the challenges in distinguishing between penicillin-resistant and -susceptible isolates based on fuzzy versus sharp zones (42). Overall, therefore, it is plausible that the genetic detection of resistance may reflect more closely the impact of the strain on a patient.
An interpretation where phenotyping reports resistance but WGS methods predict susceptibility is more difficult. One possibility is small colony variants (SCV) being present phenotypically but overgrown in WGS culture and thus not represented in the sequence. Resistance associated with gentamicin, fusidic acid, and ciprofloxacin, the main antibiotics where this phenomenon was observed, is observed with SCV phenotypes (43, 44). An alternative explanation is novel resistance mechanisms, for example, for ciprofloxacin (45), leading to false-negative WGS predictions. The need for a continuously updated curated database is a key challenge for WGS methods. As more sequencing occurs, novel mutations will be identified in resistance genes that may or may not confer phenotypic resistance, but these can at least be identified and tested; identifying entirely new resistance-conferring genes is more complex, and prediction software that can recognize new clinically important genes a priori would be a valuable addition to an analysis pipeline. However, we observed similar differences between concordant genotypic predictions and both phenotypic antimicrobial susceptibilities and single gene PCR results, suggesting that the underlying causes may not necessarily be related to resistance. As previously noted, the agreement between WGS and phenotyping is higher (98.6%) than between phenotyping undertaken by two separate laboratories (97.6%) (25); thus, at least some discrepancies are probably due to incorrect phenotyping results. In contrast, the concordance between genotypic predictions made using a single method but based on WGS generated from 5 different laboratories was recently shown to be >99.8% (46).
Limitations.This comparison was based on a prespecified set of resistance- or virulence-associated genes: some genetic traits previously associated with resistance were omitted (e.g., IleS mutations linked to low-level mupirocin resistance). Despite this, we found good agreement between genotypic predictions and phenotype. Typewriter used Velvet de novo assemblies: other newer assemblers (e.g., SPADES [47]) might have improved predictions further. We included data which had been used in the development of two of the methods compared, which could potentially have led to overfitting, although the performances of all three methods were in fact similar on this data set (see Table S2 in the supplemental material). All analyses were undertaken on short-read Illumina data. The increasing use of long-read sequences will require further software testing, although Mykrobe has been successfully used for initial resistance calling in Mycobacterium tuberculosis from Nanopore sequencing in a small number of samples (48). However, it has not been comprehensively tested, nor have Typewriter and Genefinder, with long-read sequences generated using Nanopore or PacBio technology. The greatest differences detected in this study were between phenotype and genotype, which could be partly due to the method of phenotypic testing and recognized issues with reproducibility. We did not have resources to rephenotype all or a subset of the isolates; well-characterized sets of repeatedly phenotyped isolates would be useful for further studies. We found no suggestion that missing calls in one program were associated with scores just below a threshold but did not undertake a more detailed assessment of specific sequence coverage and quality around discrepant genetic predictions.
Conclusion.In summary, in this study, the choice between three specific bioinformatics methods to identify resistance determinants or other genes in S. aureus did not prove critical. All demonstrated a high concordance with each other and with phenotypic methods and can be recommended for genotype prediction. However, each has some limitations; therefore, consensus methods provide at least some assurance. Due to computational speed, Mykrobe (de Bruijn graph based) and Genefinder (or equivalent mapping-based program such as SRST2 [22]) are a sensible combination to use as an initial consensus method, followed by Typewriter (BLAST based) if these two methods disagree. As a set of 34 diverse bacteria have been made available for whole-genome sequencing validation (49), the study strains and genotypic predictions are available as a resource for other studies investigating different bioinformatics analysis methods, which will become increasingly important as this technique is more widely used to inform clinical management, through bacterial identification, antimicrobial susceptibility prediction, and virulence profiling. External quality control of clinical laboratory performance in predicting antibiotic resistance is provided by UK proficiency testing schemes such as the United Kingdom National External Quality Assessment Service for Microbiology (UK NEQAS) (50); a similar set of standards will need to be created to accredit whole-genome sequencing methods.
ACKNOWLEDGMENTS
This research was supported by the National Institute for Health Research Health Protection Research Unit (NIHR HPRU) in Healthcare Associated Infections and Antimicrobial Resistance at Oxford University in partnership with Public Health England ([PHE] grant HPRU-2012-10041) and the NIHR Oxford Biomedical Research Centre; D.C. and T.P. are NIHR senior investigators.
This report presents independent research funded by the National Institute for Health Research and the Department of Health. The views expressed in this publication are those of the authors and not necessarily those of the NHS, the National Institute for Health Research, the Department of Health, or Public Health England.
FOOTNOTES
- Received 20 November 2017.
- Returned for modification 4 December 2017.
- Accepted 9 May 2018.
- Accepted manuscript posted online 20 June 2018.
- Address correspondence to Dona Foster, dona.foster{at}ndm.ox.ac.uk.
↵* Present address: Amy Mason, Department of Mathematics and Department of Statistics, University of Oxford, Oxford, United Kingdom; Tanya Golubchik, Wellcome Trust Centre for Human Genetics, University of Oxford, Oxford, United Kingdom; N. Claire Gordon, KEMRI-Wellcome Trust Collaborative Research Programme, Kilifi, Kenya.
A.M., D.F., P.B., T.G., and M.D. contributed equally to this article, as did A.S.W., A.K., and T.P.
Citation Mason A, Foster D, Bradley P, Golubchik T, Doumith M, Gordon NC, Pichon B, Iqbal Z, Staves P, Crook D, Walker AS, Kearns A, Peto T. 2018. Accuracy of different bioinformatics methods in detecting antibiotic resistance and virulence factors from Staphylococcus aureus whole-genome sequences. J Clin Microbiol 56:e01815-17. https://doi.org/10.1128/JCM.01815-17.
For a commentary on this article, see https://doi.org/10.1128/JCM.00813-18.
Supplemental material for this article may be found at https://doi.org/10.1128/JCM.01815-17.
REFERENCES
- Copyright © 2018 American Society for Microbiology.
