URMAP, an ultra-fast read mapper

Abstract

Mapping of reads to reference sequences is an essential step in a wide range of biological studies. The large size of datasets generated with next-generation sequencing technologies motivates the development of fast mapping software. Here, I describe URMAP, a new read mapping algorithm. URMAP is an order of magnitude faster than BWA with comparable accuracy on several validation tests. On a Genome in a Bottle (GIAB) variant calling test with 30× coverage 2×150 reads, URMAP achieves high accuracy (precision 0.998, sensitivity 0.982 and F-measure 0.990) with the strelka2 caller. However, GIAB reference variants are shown to be biased against repetitive regions which are difficult to map and may therefore pose an unrealistically easy challenge to read mappers and variant callers.

Keywords: Next generation sequencing; Read mapping.

Figure 1. Schematic of the URMAP algorithm.

Words in the plus strand of the reference sequence are indexed using a hash table with 5 bytes per row comprising a tally byte and a 32-bit pointer. Pins, i.e., words with a hash value that is unique across both strands of the reference, are indicated by a reserved tally value (pin flag). URMAP searches for a brace, i.e., a pair of pins close in the reference, one in the forward read (R1) and one in the reverse read (R2). If a brace is found, it is almost certain to be the correct location. Words found more than once in the reference are indexed using a linked list with forward pointers which are stored in tally bytes if they fit into 7 bits, otherwise in the 32-bit pointer field. The first bit of the tally is set if the row is in a list but not the head.

Figure 2. Design of the Urbench benchmark.

For each locus L in a source genome (NA12878 or GRCh38), ten simulated reads pairs are generated (five shown in figure) such that either R1 or R2 contains the locus. This enables systematic errors to be identified where a majority of reads of a given locus are mapped to the same incorrect location. With NA12878 a locus is the position of an experimentally determined variant (SNP or indel) in one of the parental chromosomes, with GRCh38 a locus is a randomly-chosen position. Base call substitution errors are introduced with probabilities given by quality scores in sequencing run SRR9091899.

Figure 3. Speed on Urbench.

Speed is measured relative to BWA with file i/o overhead minimized.

Figure 4. Mapping accuracy on Urbench.

Accuracy metrics are sensitivity and error rate with MAPQ ≥10, expressed as percentages.

Figure 5. Pair-wise method comparisons on Urbench.

Methods are sorted by decreasing total improvement (TI) (see Methods). Cells are colored according to mean improvement. A pairwise comparison of the method in row X vs. the method in column Y is given using the notation described in Methods; e.g., BWA >5(2.0) Bowtie2 means that BWA has five of eight metrics that are better than Bowtie2 with a mean improvement of 2.0. The symbols > > and < < indicate that all metrics are better or worse, respectively, e.g., URMAP > >(4.4) FSVA means that URMAP is better than FSVA by all metrics with a mean improvement of 4.4.

Figure 6. Scatterplot of reported MAPQ vs. measured MAPQ.

For each integer value of MAPQ, the measured MAPQ is determined by the frequency of incorrectly mapped reads in the Urbench benchmark. Hisat2 is not shown because it reports only three distinct MAPQ values (see main text).

Figure 7. Mapping accuracy on wgsim test.

Accuracy metrics are sensitivity and error rate with MAPQ ≥10, expressed as percentages. Tests were performed with three simulated read lengths: (A) 150, (B) 250 and (C) 300, respectively.

Figure 8. Mapper agreement on unmappable regions in the human reference genome.

Venn diagram showing agreement of BWA, Bowtie2 and URMAP on unmappable regions with 2 ×150 reads of GRCh38 with MAPQ ≤3. These mappers agree that 49.6M bases (intersection of the three regions) are not mappable.