Rapid and precise alignment of raw reads against redundant databases with KMA

Heuristic k-mer mapping

When mapping against databases, a query sequence can fall within two categories: the query does not look like anything in the database or the query looks like something in the database. Depending on the database, it is therefore necessary to determine at an early stage whether alignment of a given sequence is feasible or if the sequence should be discarded from further analysis. Computationally, it is much faster to map raw reads against template sequences and determine their approximate origin than it is to perform alignment giving the exact differences between sequences. In order to exploit this, KMA maps k-mers between the query sequence and template database. The matching k-mers are then used to produce a mapping score where a matching k-mer is rewarded with a score of k while k-mer mismatches are penalized with − 1. For a query sequence to be accepted for alignment, it must have a positive score based on this mapping, equivalent to a minimum identity of (k + 1)^− 1 in k-mer space. In the case where a query sequence is accepted, the potential template candidates are limited by linking the matching k-mers to their respective template sequences.

Fine alignment

KMA uses a hash map of indexed k-mers from which to start alignment, where this indexing is computed for templates passing the heuristic k-mer mapping. The indexed k-mers are then used as seeds from which to start alignment, where each matching seed is extended to give an optimal score. To enable a high resolution of gaps and mismatches, KMA utilizes the Needleman-Wunsch algorithm [9] to align the pieces of query and template sequence between seed-extends. The tails of the alignment, before the first seed and after the last, are aligned with a specialized version of the Needleman-Wunsch algorithm allowing for early stopping so that leading and trailing gaps can be skipped.

ConClave scoring

As KMA is designed to map sequences against redundant databases, the possibility of a tie for best matching template should be considered more as rule than an exception. For KMA to be able to find the one most likely template for each query sequence, a novel scoring scheme, ConClave, was developed.

The ConClave scoring scheme requires that all best matching templates are initially accepted for all best matching query sequences. The alignment scores obtained are then summed together for each template candidate in order to form the ConClave score, in this way reflecting a maximum score for each template. For alignments where ties do not occur, a second score is conducted as the unique ConClave score, which are used to reflect the unique alignments.

Ties for best matching templates can now be resolved by choosing the most likely template based on the ConClave score. In case of a tie on the ConClave score, the unique ConClave score is used to choose a best matching template. For samples were two templates are equally well-suited candidates, reflected by identical ConClave scores, the algorithm will choose the first sequence added to the database.

Reference assembly

The ConClave scoring and alignment of reads enables assembly guided by a reference, resulting in a consensus sequence for the given reference / template. From the alignment, the differences between each read and the associated template are known, and this gives information about which nucleotides are called at each position in the template. Each nucleotide is then determined by majority voting, where the strength of the called nucleotide is tested with a McNemar test, giving a robust base calling performance across different sequencing platforms [10]. In case a nucleotide is not significantly overrepresented according to the McNemar test (according to a user specified α, default 0.05), it is reported in lower case.

After reference assembly, statistics for each template are collected in a single file, containing features such as template identity, coverage and depth.

Performance evaluation and comparison

The performance of KMA was evaluated on a dataset of 13 Escherichia coli from Grad et al. (2012). The dataset included 11 E. coli from the German / French outbreak in 2011 and two historical E. coli [11]. Phenotypic susceptibility tests were known for all of the 13 E. coli, and the presence of the beta-lactamase blaCTX-M-15 were verified by PCR [11].

The 13 E. coli were mapped to the ResFinder database of known resistance genes [12], where the associated phenotype of each gene were compared with phenotypic susceptibility tests.

In order to measure the performance on larger databases, a core genome MLST (cgMLST) database for E. coli was downloaded from EnteroBase (available at: http://enterobase.warwick.ac.uk [Accessed 18 January 2018]. The cgMLST scheme contains 2447741 closely related genes, which together matches the size of the human genome.

KMA was compared with six different methods: SRST2 (v. 0.1.8) [3], MGmapper (v. 2.7) [7], BWA-MEM (v- 0.7.15) [5], Bowtie2 (v. 2.2.4) [4], MiniMap2 (v. 2.6) [6] and Salmon (v. 0.9.1) [8], each capable of mapping raw reads directly against entire databases of sequences.

SRST2 uses Bowtie2 to map against a homology-reduced database, clustered at 90% identity using CD-hit [13]. From each cluster, one allele is chosen, as the one with the highest alignment score from all alignments towards that cluster. After mapping and alignment, SRST2 performs SNP calling and reports the differences between query and template sequences [3].

For MGmapper to accept a template, it has to contain at least one unique mapping read when mapped with BWA-MEM. After mapping the reads, MGmapper reports the depth and template coverage of the accepted genes [7].

Salmon was developed for quantifying RNA transcript levels using the EM-algorithm. Which partially matches the redundant database problem, as only a subset of the templates in the database are likely to be present. Salmon can be used as an extension to BWA-MEM, Bowtie2 and Minimap2, when the alignment method is set to report several alignments per read [8].

SAMtools and BEDTools were used to estimate depth and template coverage of mappings performed with BWA-MEM, Bowtie2 and Minimap2 so that these could be used for a more direct comparison as well [14, 15].

To prove the difficulty of mapping raw reads directly against redundant databases, a simulated dataset of single-end and paired-end reads was created, where each gene in the ResFinder database was split into raw reads with a length 100 bp, and an insert size of 250 for the paired end set. Error and indel rates were set to zero so that the simulated reads would represent the perfect case scenario of short-read sequencers. Meaning, that any inconsistencies from the mapping methods, would reveal each of their weaknesses in mapping against redundant databases, and give a clear view of their assumptions. Simulated reads were created with chopDB, which is included alongside the commands used to compare KMA with existing methods and the ResFinder database.

Simulated datasets

For the simulated datasets, the comparison was straight forward as the location of every read was known.

For both the single and paired-end datasets, KMA mapped everything correctly, where every read in each sample was mapped to the correct template only (see Table 1). For all samples, SRST2 mapped almost all reads to the correct template, only differing by a read in most cases. However, since SRST2 allows for one read to map towards several templates and even several places within each template, a large number of false positives were achieved, thus giving a lower correlation score (see Table 1). For a few genes, SRST2 accepted up to 21 alignments (see Fig. 2). MGmapper did not assign any false positives, but underestimated the depth of the correct templates. This was especially the case for the paired-end dataset where MGmapper seemingly had trouble with coupling of the reads. This resulted in MGmapper only being able to assign 19.6% of the actual reads. BWA-MEM and MiniMap2 follow a similar pattern compared with MGmapper, where the methods are driven towards false negatives (FN). They do however perform better on the paired end reads. In contrast to the other methods tested, Bowtie2 has a more equal distribution between false positives (FP) and FN. Only BWA-MEM had improved predictions when post processed with Salmon, but only fro the single end reads, whereas Bowtie2 had a performance decrease, on both the single and paired end datasets (see Table 1 and Fig. 2).

Method / Performance	Single end read set			Paired end read set
	MCC	Sensitivity	PPV	MCC	Sensitivity	PPV
KMA	1.000	1.000	1.000	1.000	1.000	1.000
SRST2	0.591	0.999	0.350	0.659	0.999	0.436
MGmapper	0.676	0.457	1.000	0.443	0.196	1.000
BWA-MEM	0.585	0.342	1.000	0.580	0.337	1.000
Bowtie2	0.480	0.480	0.480	0.577	0.577	0.577
Minimap2	0.591	0.353	0.988	0.671	0.455	0.991
BWA-MEM / Salmon	0.720	0.720	0.720	0.500	0.353	0.707
Bowtie2 / Salmon	0.390	0.389	0.398	0.250	0.177	0.368

MCC Matthews correlation coefficient, PPV Positive Prediction Value

German / French E. coli outbreak

Ninety percent template coverage, 90% template identity and a minimum depth of 5 were the thresholds for predicting genes with KMA and SRST2. The identity threshold was skipped for MGmapper and methods post processed by SAMtools and BEDTools because introducing an identity threshold would require further post-processing, while only giving lower correlation scores. Only the depth threshold was used for Salmon, because Salmon only estimates the transcript levels of each template. Meaning that the exact location of each read is lost, making the analysis of template coverage and identity impossible.

MGmapper crashed when mapping the paired end reads, but succeeded in single end mode when only the forward reads were used. Paired end reads were converted to interleaved format for BWA-MEM because it crashed when the reads were split into two files.

Predicting antimicrobial resistance

Full concordance between phenotypic susceptibility tests and the predicted antimicrobial resistance genes were achieved for both KMA and SRST2. Results for nalidixic acid resistance were excluded, since it is caused by chromosomal resistance genes that were not included in this study [16]. The only disagreement between KMA and SRST2 for this dataset were that SRST2 predicted the presence of dfrA17 and dfrA32 in addition to the dfrA7 gene predicted by KMA, in all of the outbreak strains (see Additional file 1: Table S1). A closer examination of this result revealed that the two additional genes predicted by SRST2 were due to multi mapping reads between them and dfrA7.

As for the simulated dataset, MGmapper had no FP predictions, but it missed most of the phenotypes. Even though the outbreak strains are thought to be highly similar, the results from MGmapper showed a highly diverse pattern between the samples (see Additional file 1: Table S1). Neither BWA-MEM, Bowtie2 nor Minimap2 predicted any genes exceeding the thresholds, no matter the settings of SAMtools and BEDTools (see Additional file 1: Table S2).

Mapping towards cgMLST alleles

The performance of each method was measured as their ability to find one allele within each locus of the cgMLST database, which together make the core genome of the bacteria.

As for the antimicrobial resistance genes, BWA-MEM, Bowtie2 and Minimap2 had trouble assigning genes exceeding the thresholds, thus giving a low correlation coefficient (see Table 2). The performance of MGmapper was slightly better than that of BWA-MEM used directly, giving a correlation of 0.062 compared with BWA-MEM with a correlation of 0.021 (see Table 2). As for the resistance gene prediction, the correlation of BWA-MEM was improved when using Salmon. In this case, however, Salmon did take up the majority of the computational resources when compared with the mapping. When estimating the allele abundances, Salmon had a peak memory consumption of 104.2 GB, and used more time on the post-processing than BWA-MEM used on mapping (see Table 2). KMA outperformed all of the methods with a correlation of 0.998, only missing two genes and detecting three duplicates on average.

SRST2 and Bowtie2 combined with Salmon were stopped after 100 h, where neither method had completed a single sample.

Runtime analysis

CPU time and peak memory requirements were measured with GNU time (v1.7), mapping and post-processing were measured separately for all methods, unless included in the method. KMA outperformed all of the existing methods on speed when mapping both resistance genes and cgMLST genes. The memory usage of KMA was above that of Bowtie2, but just below that of MGmapper and BWA-MEM, and less than a fourth of MiniMap2 (see Table 2). The memory consumption and speed of KMA is thus feasible to carry out on a laptop, with results accurate enough to be used in surveillance.