SRST2: Rapid genomic surveillance for public health and hospital microbiology labs

Implementation

Given a read set and database of reference allele sequences, SRST2 is designed to perform two key tasks: (1) detect the presence of a gene or locus; and (2) determine the precise or closest matching allele for that locus, among a set of possible reference allele sequences. The approach is illustrated in Figure 1. A database of reference sequences must be provided in fasta format, in which the fasta headers indicate both the locus (so that alleles of the same locus can be compared) and a unique name for each allele. In the case of MLST data an additional database of ST profiles is provided as tab-delimited text, which assigns STs to unique combinations of alleles. Current MLST data (allele sequences and profile definitions) can be downloaded from pubmlst.org automatically using the getmlst.py script supplied with SRST2. For most MLST schemes, these files are compatible with SRST2 and can be used without modification. Improperly formatted databases, or other private sequence databases, can be easily formatted for use with SRST2 using the scripts supplied with the program. Any number of sequence databases can be analysed in a single run, allowing for simultaneous typing of MLST, resistance genes and virulence genes.

For each input database, reads are aligned using bowtie2[29] v2.1.0 or above with the `--very-sensitive-local and `-a settings, and all alignments are reported to a file in SAM format. Mapping sensitivity can be fine-tuned by specifying to SRST2 any of the parameters available within the bowtie2-align command or a maximum number of mismatches per read (default 10 mismatches allowed). Flags in the resulting SAM file are modified so that each read is included in the pileup for every allele to which it is aligned. Pileups are generated using SAMtools v0.1.18 [30]mpileup and parsed by SRST2 to determine percent coverage, divergence, and mismatches as well as to calculate a score for each possible allele.

Allele scoring

An overview of the scoring approach is given in Figure 1. We begin with an alignment of reads from sample s to a reference sequence r. At each position i in the reference sequence r (r _i ), let s _i be the set of reads in sample s that align to r _i . Let a _i be the total number of reads in s _i , and let b _i be the number of reads in s _i in which the aligned base does not match the reference base at r _i . If sample s contains the precise sequence r, then the probability of a mismatched base at any position in an aligned read is equal to the per-base error rate of the sequencing technology e _i , which for Illumina is taken to be 0.01, although this can vary depending on what preprocessing steps are implemented [31],[32].

To quantify the evidence against the presence of the reference sequence r in s, we perform a Binomial test at each position r _i , to generate a one-sided P value P _i to assess the probability of observing a _i -b _i successes in a _i trials, with a probability of success of 1-e _i . Any change at position r _i - including a base substitution, an insertion of any size or a deleted base - is treated as a mismatch, incrementing b _i by 1. For large deletions that result in an absence of any aligned reads (including truncations of the end of the sequence), a _i = 0 and no binomial test is possible. In this case, the evidence for the deletion is provided by the reads which align adjacent to the deletion but do not align across the deletion. Hence we calculate the average number of reads aligned to the two bases preceding the deletion, d _i , and conduct the binomial test with a _i ?=?b _i ?=?d _i .

We then compare the sorted -log₁₀(P _i ) values versus those of the theoretical distribution of -log₁₀(x _j /n) where n?=?length(r) and x _j = 1,2,n, analogous to a quantile-quantile (QQ) plot (Figure 1). A linear model is fitted to the two probability distributions and the resulting slope is taken as the score for reference sequence r, score _r . Here we leverage a common criticism of linear models to our advantage: the susceptibility to outliers at the tails of the distribution. In this case, outliers are typically SNPs or indels relative to the sequence r which, because they result in low P values in the binomial test and thus very high values of -log₁₀(P), are at the end of the observed distribution (Figure 1). Thus when a linear model is fitted, its slope increases with the number of well-supported SNPs and indels compared to the reference. As a result, among reference alleles of the same locus, the sequence r with the lowest score _r (flattest slope in the QQ plot) is the most likely match for sample s.

Reporting outputs

For each sample s and each locus or gene cluster, SRST2 output tables report the lowest scoring allele sequence r, the average read depth of s across r and indicators of any evidence against a precise match with r (including mismatches supported by >50% of aligned reads, or read depth falling below a cutoff). Only matches passing the user-set coverage and divergence cutoffs (by default, >90% coverage and <10% divergence) are reported. For MLST data, STs are calculated according to the MLST profiles database provided, based on the closest matching alleles at each locus.

Normally, an exact match between r and s would be assigned if (a) r has the lowest score _r among the set of alleles of the same locus or gene cluster, and (b) there are no SNPs or indels between r and s. If (a) holds but (b) does not, this is indicative of a novel allele and SRST2 will flag the result in output tables. In such cases, we recommend that users who are interested in defining novel alleles should inspect the raw sequence data (which may be assisted by the alignments, pileups and consensus fastq files generated by SRST2).

Optionally, SRST2 can report the full details of scoring s against all reference sequences r, to enable users to parse and interpret the results to suit specific needs. These include average depth of s across r, average depth across the first and last two bases of r, the number of positions in r in which the majority of aligned reads in s show a mismatch against r (with SNPs, insertion/deletions and truncations reported separately), the depth of bases neighbouring truncations and, for the position with the greatest proportion of mismatching reads, the total aligned reads, total mismatching, proportion mismatching, and binomial P value.

Major differences between SRST and SRST2

SRST2 is new code and takes an entirely different approach to read mapping, scoring alignments and reporting results than SRST [28], which was designed solely for MLST and is unsuitable for detection of acquired genes. In SRST, bwa was used for global alignment of reads to MLST loci and their flanking sequences; in SRST2, bowtie2 is used for local alignment of reads to any locus, without need for flanking sequences, allowing detection of acquired genes as well as MLST. SRST scores were calculated in an entirely different way and were not designed to take into account deletions/truncations or the relative weight of evidence provided by each position in the alignment (differences in read depth). SRST2 allows finer control of mapping and scoring parameters and provides more detailed reports than SRST. SRST2 is also faster (2 to 5) and slightly more reliable than SRST for MLST analysis (see below).

Bacterial isolates and sequencing

A total of 231 Listeria monocytogenes isolates were analysed in this study, at the Microbiological Diagnostic Unit (MDU) Public Health Laboratory in Victoria, Australia. MDU is the national reference laboratory for L. monocytogenes and the isolates analysed include several from recent outbreaks as well as from the laboratory's reference collection. Ethical approval was not required for the use of reference laboratory isolates in this project. Cultures of L. monocytogenes isolated from food, environmental or clinical specimens were purified by two successive single colony selections after streaking onto horse blood agar (HBA) incubated for 18 to 24 h at 37 °C. Resultant bacterial growth on the surface of HBA medium was aseptically collected and resuspended in a cryotube (Nalgene) containing 1mL of sterile glycerol storage broth (1.6% w/v Tryptone, Oxoid Pty Ltd, LP0042 containing 20% v/v glycerol) prior to storage at -70 °C. Cultures were retrieved from storage as required and freshly grown (HBA, 18 to 24 h at 37 °C) in preparation for DNA extraction. DNA was extracted from each isolate using QIAmp DNA Mini Kit (Qiagen) and eluted in EB buffer (Qiagen) (Tris buffer, no EDTA).

DNA samples were subjected to traditional L. monocytogenes MLST analysis [33],[34], with a minor modification to the annealing temperature for the bglA PCR (52 °C not 45 °C). The PCR products were purified with FastAP Thermosensitive Alkaline Phosphatase (Thermo Scientific) and Exonuclease I (Thermo Scientific). The purified PCR products were sequenced using BigDye Terminator v3 chemistry followed by capillary sequencing using a 3130xL Genetic Analyzer (Applied Biosystems). Trace analysis was conducted using BioNumerics version 6.6 with MLST Online plugin version 2.13 and Batch Sequence Assembly plugin version 1.34.

DNA was subjected to multiplex library preparation using Nextera XT followed by sequencing using an Illumina MiSeq. DNA was quantified by Qubit dsDNA HS Assay Kit (Invitrogen) and normalised to 0.2ng/?L. Total 1ng of DNA was used for Nextera XT DNA Sample Preparation Kit (Illumina). Tagmentation of genomic DNA, PCR amplification with dual index primers, PCR clean-up using Agencourt AMPure XP (Beckman Coulter), DNA libraries normalization, library pooling and MiSeq sample loading were performed according to the manufacturer 's instruction with minor modifications. For longer than 2250bp runs on the MiSeq, 25L of AMPure XP beads was added to each PCR-amplified product during the PCR purification step otherwise 30L of AMPure XP beads was added. For some samples, after PCR purification, DNA fragment size and library concentration was analysed by 2100 Bioanalyzer (Agilent Technologies) and Qubit dsDNA HS Assay Kit (Invitrogen). DNA libraries were normalized manually to 4 nM and libraries with unique indexes were pooled in equal volumes. Each resulting pooled library was denatured and diluted with 0.2N NaOH and pre-chilled HT1 (Illumina) to produce a 20 pM denatured library in 10mM NaOH. Prior to the MiSeq run, the denatured library was further diluted with pre-chilled HT1 to approximately 12 to 13.5 pM. A total of 600L of library including 2% (v/v) 20 pM denatured PhiX library (Illumina) was loaded together with MiSeq reagent kit v3 (Illumina) according to the manufacturer 's instructions.

Publicly available short read data used in this study

Subsampling of read sets

To explore accuracy at low read depths, 10 genomes each of S. aureus and E. faecium were selected for random subsampling of reads to simulate genomes sequenced to low read depth. To do this, we used the mean read depth across MLST loci to calculate the sampling fraction required to achieve approximately 10, 20, 100 mean read depth. We randomly sampled reads from the forward reads file at the required sampling fraction, and extracted the corresponding reverse reads, using Perl scripts. Ten random samples were generated from each read set at each depth level, generating a total of 1,000 read sets for each species.

Sequence databases used in this study

MLST databases for Staphylococcus aureus, Streptococcus pneumoniae, Salmonella enterica, Escherichia coli, Enterococcus faecium, Listeria monocytogenes and Enterobacter cloaceae were downloaded from pubmlst.org using the getmlst.py script included with SRST2 (June 2014).

Antimicrobial resistance gene detection was performed using the ARG-Annot database of acquired resistance genes [18]. Allele sequences (DNA) were downloaded in fasta format [43] (May, 2014). Sequences were clustered into gene groups with ?80% identity using CD-hit [44] and the headers formatted for use with SRST2 using the scripts provided (cdhit_to_csv.py, csv_to_gene_db.py). A copy of the formatted sequence database used in this study is included in the SRST2 github repository [35].

Representative sequences for 18 plasmid replicons were extracted from GenBank using the accessions and primer sequences specified by Carattoli et al.[45]. A copy of the formatted sequence database used in this study is included in the SRST2 github repository [35].

Simulation of expanded S. aureus MLST database

As more genomes are sequenced and as bacteria continue to evolve, novel alleles will continue to be discovered and thus the size of allele databases will increase. To explore the impact of database size on accuracy of allele detection with SRST2, we simulated expansion of the current S. aureus MLST database from 2,161 alleles (mean 309 per locus) to 5,578 alleles (mean 797 per locus). The additional 500 alleles (approximately) per locus were generated using netrecodon v6.0.0 [46]. Sequences derived from the true MLST database were used to seed the simulation at each locus as follows. Existing alleles were translation-aligned between start (alignment start) and stop (alignment end) codons, those containing a frameshift or stop codon were removed, and the modal consensus sequence was exported. The best-fit DNA substitution model of each true alignment was determined using the AIC in MrModeltest v2.3, as implemented in PAUP* v4.0b. In netrecodon, the modal sequences were forward evolved under the coalescent, using the parameters of the best-fitting model for each locus, mutation rate 1E-7 and recombination rate 1E-7/15 (based on reported r/m of 1/15 [47]). A total of 100 independent replicates of forward evolution were performed per locus, retaining 2,000 sequences per replicate (N = 200,000 simulated sequences per locus). The first 500 unique simulated sequences at each locus were added to the MLST database, and duplicate sequences were removed.

Assembly-based analysis

Assemblies were generated using the de novo assembler Velvet v1.2.10 [24], with optimal kmer choice for each read set refined through iterative calls to VelvetOptimiser v2.2.5 [48]. Briefly, each read set was assembled using a call to VelvetOptimiser with kmers from 29 up to 89, in steps of 12. The optimal kmer, k₁, was extracted and a second call to VelvetOptimiser was made using kmers from k₁-12 up to k₁?+?12, in steps of 4. A final call to VelvetOptimiser was run using kmers from k₂-4 up to k₂?+?4, in steps of 2. The final assembly was that output from the third and final call to VelvetOptimiser.

For MLST analysis from assemblies, a nucleotide BLAST?+?(v2.2.25) search was performed for each locus and each contig set. In this BLAST search, the contig set was used to query the database containing all known allele sequences for a given locus, and the top BLAST hit was reported. If this hit had ?90% nucleotide identity across ?90% of the length of the reference allele sequence, an allele call was recorded. If the hit was an exact match to a known allele (that is, 100% nucleotide identity across 100% of the length of the allele sequence), this was considered a precise allele call. The Python code used is available within the SRST2 distribution. Where the hit was not an exact match (n = 42), an alternative nucleotide BLAST analysis was performed using the allele sequences as query and the contig set as database, and the results manually inspected to determine whether it was possible to identify the correct allele from the assembly. For gene detection analysis from assemblies, a nucleotide BLAST search was performed in which the set of reference sequences (sequence database, that is antimicrobial resistance gene database) was used to query the database of all contigs for that assembly.

SRST (v1) analysis

The 543 read sets used for validation of SRST2 allele calling (Table 1) were also analysed using SRST [28], run with default settings.

Web-based analysis with MLST Typer and ResFinder

The 44?S. enterica read sets (Table 1) were analysed using the MLST Typer [21] and ResFinder [19] websites. This data set was chosen to begin with as it is the smallest of those used for validation in the manuscript (n = 44). Each read set took 3 to 4 h to upload and analyse using these websites, and had to be done in serial as attempting to run multiple jobs crashed Java. Therefore it was not feasible to test all read sets for comparison.

Analysis runs and time calculations

All SRST2, SRST, assembly and BLAST analysis was run on a Linux cluster (iDataplex86 system, `Barcoo cluster at VLSCI [49]). SRST2 was run with default parameters. Details of Velvet assembly and BLAST analysis are given below. Run times were calculated from time stamps extracted from log files for SRST2 and Velvet Optimiser assembly runs.

Statistical analysis

All statistical analysis and data plotting was performed in R. Allele calling performance of SRST2 and assembly?+?BLAST was assessed via three metrics: (1) call rate?=?total number of allele calls made, for SRST2 this was a call with ?90% coverage and no uncertainty recorded (that is, with 2 read depth at both ends and also neighbouring any truncations or deleted bases), for BLAST this was a call with ?90% coverage and ?90% nucleotide identity; (2) false positive rate?=?total number of correct allele calls as a proportion of all calls; and (3) proportion of all tests resulting in a call with a correct allele, equal to (call rate) * (1 - (false positive rate)). As these metrics are proportions, the significance of differences in performance metrics was calculated using a two-sided test for equality of proportions (prop.test function in R). Resistance gene detection was assessed using a cutoff of ?90% coverage and ?90% identity to define the presence of a gene.

Validation of allele calling

To assess the accuracy of allele identification with SRST2, we analysed publicly available Illumina data from 543 bacterial genomes of five different species for which independent MLST data were available (Table 1). With seven loci in each MLST scheme, this yielded 3,801 allele calls across 35 loci to assess call rate and false positive rate. The read sets represented a wide range of average read depths, with 90% in the range 12 - 130 and 50% between 20 - 60 (Table 1). For each species, we used SRST2 to download the latest MLST database from pubmlst.org and subsequently ran SRST2 using default parameters. Median run time was 6min per sample (interquartile range, 4 to 10min) and increased linearly with number of reads (Figure 2). Efficiency can be easily improved or standardised, without data preprocessing, by instructing SRST2 to map the first N reads only.

SRST2 call rates and true positive rates increased with average read depth, stabilizing with depths ?15 (Figure 3a, Additional file 2). For comparison, we also assembled each read set using Velvet[24] and VelvetOptimiser[48] and used nucleotide BLAST to identify MLST alleles (assembly?+?BLAST method; Methods). At read depths 15, SRST2 made significantly more allele calls than assembly?+?BLAST (call rates 99.9% vs. 95.9%, respectively; P <110^-15), with equivalent accuracy (0.46% vs. 0.22% of allele calls incorrect; P = 0.16). The heuristic information provided by SRST2 (that is, confident mismatches, insertions, deletions or truncations reported from read mapping) was a strong indicator of accuracy in the result: where an exact match was reported (98% of calls with depth ?15), 0.2% of allele calls were incorrect; where an inexact match was reported, 11.7% of allele calls were incorrect. Hence, the key difference between the two methods was the ability of SRST2 to make correct calls where assembly?+?BLAST could not make any call: for read depths ?15, SRST2 made a call with the correct allele 99.4% of the time, compared to only 95.7% for assembly analysis (P <110^-15 for difference in frequencies of correct allele calls). At sequence type (ST) level, the difference was even greater: SRST2 achieved accurate ST assignment for 98% of isolates with average depth ?15, whereas assembly?+?BLAST correctly identified only 79%. SRST2 also performed better than SRST and the MLST typer website (which implements an alternative assembly?+?blast approach), see Tables 4 and 5.

To assess performance at low read depths (?15), 10?S. aureus read sets were subsampled to low depths (Methods). This confirmed that an average depth of only 10 was required for SRST2 to achieve >90% call rate and <0.5% false positives (Figure 3a, Figure 4). MLST databases can be expected to grow indefinitely due to increasing diversity and broader sampling. However simulations (Methods) indicated that doubling the size of the S. aureus MLST database had no impact on SRST2 accuracy (Figure 3a, Figure 4).

Validation of gene detection using the vanA-B resistance gene

In addition to reliably distinguishing alleles of a given gene, SRST2 can also accurately determine the presence or absence of genes of interest, such as those encoding antimicrobial resistance or virulence. To evaluate this, we used 43 E. faecium genomes (Table 1), previously screened for vancomycin susceptibility and presence of the VanB vancomycin resistance operon vanABHSXY[41],[50]. Seventeen isolates were vancomycin resistant (VRE), and all were PCR positive for the vanA-B gene. These genomes were sequenced to approximately 1,000 depth and SRST2 correctly detected vanA-B in 17/17 VRE. In five vancomycin sensitive (VSE) isolates PCR negative for vanA-B, SRST2 detected VanA-B sequences at very low depths (<0.2% of average depth), probably caused by minor but easily identifiable contamination during VRE-VSE multiplexed sequencing. SRST2 also confirmed the presence of the entire VanB operon, which is strongly predictive of the VRE phenotype. For comparison, assembly?+?BLAST identified full-length vanA-B sequences in just 7/17 VRE genomes, with multiple smaller hits spanning the full-length gene in five VRE and <50% coverage of the gene identified in the remaining five VRE. To investigate the effect of sequencing depth on gene detection, we randomly selected five VRE and five VSE read sets for subsampling at <10 average read depth. VanA-B was only ever detected in confirmed VRE genomes, and sensitivity of detection with SRST2 reached 100% for read sets with ?5 average read depth (Figure 3c).

To further explore the relative sensitivity of gene detection with SRST2, we screened all the read sets used for MLST validation (Table 1) for antimicrobial resistance genes in the ARG-Annot database of acquired resistance genes [18] (Methods). SRST2 's detection of whole genes was more sensitive than detection of whole or partial gene sequences by assembly?+?BLAST (Figure 5): 6.8% of genes detected at ?90% coverage by SRST2 at depths ?15 were not found at ?90% coverage in assemblies. For most of these genes, smaller fragments were detected by BLAST (Figure 5); however, SRST2 has the advantage of sensitive detection and confident allele-calling across the full length of genes, even at low depths (Figures 3c and 5). SRST2 also performed substantially better than the ResFinder website, which implements an alternative assembly?+?blast approach (Table 3).

Validation of SRSTin a public health laboratory

To validate SRST2 in a public health laboratory setting, we analysed 231 clinical isolates of Listeria monocytogenes and compared MLST data obtained from gold-standard PCR and amplicon sequencing with those obtained from SRST2 or assembly?+?BLAST analysis of Illumina MiSeq data (Figure 3b). Sequencing and analysis were performed by the Microbiological Diagnostic Unit Public Health Laboratory in Melbourne, Australia, the national reference laboratory for L. monocytogenes. For average read depths ?15, SRST2 had a substantially higher call rate than assembly-based analysis (99.6% vs. 95.7%; P <110^-12), with similar low false positive rates (0.7% vs. 0.6%; P = 0.9). Hence, for samples with ?15 data, a total of 99% of all alleles were called correctly by SRST2, a significantly higher proportion than the 95% achieved by assembly?+?BLAST (P <110^-12). At <15 read depths, SRST2 also performed better than assembly-based analysis (87% vs. 72% of alleles correctly called, respectively, P <110^-3; Figure 3b).

Further, SRST2 is already being assessed for routine MLST analysis of Streptococcus pneumoniae at Public Health England (Anthony Underwood, personal communication), and the open-source SRST2 code has been adapted by Public Health Ontario, Canada to perform specialist emm typing of Group A Streptococcus[51].

Identification of antimicrobial resistant clones

In a hospital setting, the combination of MLST and gene detection can provide rapid and powerful insights for infection control without specialist bioinformatics knowledge. SRST2 analysis of 69?K. pneumoniae and 74 E. coli genomes from a UK hospital [42] revealed that each was dominated by a single ST with a high rate of the extended-spectrum beta-lactamase (ESBL) gene CTX-M-15 (K. pneumoniae ST490 comprising 25% of total, 71% of ESBL; E. coli ST131 comprising 40% of total, 77% of ESBL; Figure 6). Routine SRST2 surveillance of ESBL infections could be indicative of hospital outbreaks and used to identify which isolates should be investigated via transmission analysis.

Using the E. faecium genome data, collected as part of a 12-year hospital study of vancomycin resistance [41], SRST2 took approximately 30min to generate the results and visualisations shown in Figure 7, indicating: (1) increasing vancomycin resistance over time; (2) a shift in dominant ST during the same period; and importantly (3) that this was not attributable to the introduction nor transmission of a new resistant clone, as the resistance rates were steady (approximately 50%) across all dominant STs. Similar conclusions typically require many days of labour and specialised assays in the diagnostic laboratory [52] and have been confirmed by detailed WGS analysis showing frequent acquisition of VanB transposons by diverse circulating strains [41].

Investigation of outbreaks and carbapenem resistance mechanisms

We next applied SRST2 to analyse data from real-world small-scale infection control investigations [15]. SRST2 took 5min to generate results for suspected outbreaks of VRE and E. cloaceae (Figure 8), in which suspected outbreak isolates were readily distinguishable from epidemiologically unrelated isolates, consistent with WGS phylogenies and manual analysis of antimicrobial resistance markers [15]. SRST2 typing of 18 plasmid replicons [45] also indicated specific plasmid replicons (IncHI2, IncA/C) associated with two of the resistance profiles. The authors also reported use of a complex hybrid of assembly, mapping and manual inspection to determine carbapenem resistance mechanisms in five Gram-negative bacteria isolated in close proximity [15]. SRST2 analysis of these five read sets identified the acquired beta-lactamases OXA-23 in AB223; IMP, SHV-12 and TEM-1 in EC1a; CTX-M-15 and TEM-1 in Eco216; CTX-M-15 and SHV-133 in KP652; and no acquired carbapenemase genes in EC302. These results are consistent with those reported from manual analysis [15].