Massively parallel read mapping on GPUs with the q-group index and PEANUT

Designing for efficient GPU usage

We use the terminology of NVIDIA™, while the general concepts are also applicable to competitors like AMD™. A GPU is partitioned into streaming multiprocessors (SMs), each of which has its own on-chip memory, cache and processing cores. By adjusting the thread block size it can be controlled how threads are distributed among the SMs. One thread block is executed on one SM and continues resident until all threads in the block are completed. Once a thread block has finished, another will be scheduled to the SM if any blocks are left. An SM can execute 32 threads in parallel (restricting the thread block size to be a multiple of 32); such a group of threads is called a warp or wavefront. At any time, each of these threads has to execute the same instruction in the code, but may do so on different data, a concept which is called single instruction, multiple threads (SIMT). Hence, conditionals with diverging branches should be avoided, since threads taking an if-branch have to wait for threads taking the corresponding else-branch to finish and vice versa. All SMs may access a slow common global memory (about ≤3 GB on most of today’s GPUs) in addition to their fast on-chip cache and memory. While the size of the fast cache is extremely limited, accesses to the global memory are slow and should be minimized. However, the memory latency can be reduced by coalescing the access, i.e., letting threads in a warp access contiguous memory addresses, such that the same memory transaction can serve many threads. In addition, an SM can execute a different warp while waiting on a transaction to finish, thereby hiding the latency. For the latter, threads should minimize their register usage such that the number of warps that can reside on an SM is maximized. Finally, data transfers from the main system memory to the GPU global memory are comparatively slow. Hence it is advisable to minimize them as well.

A useful programming pattern that is used extensively in the algorithms presented in this work are parallel prefix scans, a special case of which is the computation of a cumulative sum, which at first appears to be a serial process, or filtering an index set. Parallel prefix scans are used to solve these problems in a data parallel way with a minimum amount of branching, thus nicely fitting the above considerations; Blelloch (1990) provides a comprehensive introduction. We use the PyOpenCL implementation of prefix scans (Klöckner et al., 2012).

The q-group index presented in this work is designed with above considerations in mind: The data structure supports simple access patterns that minimize the register usage. This allows to hide latency caused by not coalesced memory accesses. We evaluate this in “GPU resource usage”. The parallel algorithms presented in this work further minimize the data transfers necessary between their steps to constant amounts of data (e.g., single integers). Where appropriate, we discuss architecture specific considerations in the text.

The q-group index

A classical DNA q-gram index of a text T stores, for each string of length q over the alphabet, at which positions in T the q-gram occurs and allows retrieving these positions in constant time per position. It is commonly implemented via two arrays that we call the address table A and the position table P.

Q-grams are encoded as machine words of appropriate size with two subsequent bits encoding one genomic letter (i.e., A = 00, C = 01, G = 10, T = 11). Unknown nucleotides (usually encoded as N) are converted randomly to A, C, G or T, and larger subsequences of Ns are omitted from the index. Hence, a q-gram needs 2q bits in hardware and is represented (encoded) as a number g∈{0, …, 4^q−1}. The address table provides for each (encoded) q-gram g a starting index A[g] that points into the position table such that P[A[g]], P[A[g] + 1], …, P[A[g + 1]−1] are the occurrence positions of g.

Deciding about the q-gram length q entails a tradeoff between specificity of the q-grams and the size of the data structure. Array A needs 4^q integers and thus grows exponentially with q, while array P needs |T| integers, independently of q. Larger values of q lead to fewer hits per q-gram that need to be validated or rejected in later stages. Further, the choice of q determines the sensitivity or error-tolerance of the search via the pigeonhole principle (q-gram lemma): only if there are e < (n + 1)/q−1 errors, we can guarantee that a q-gram match exists.

Idea

The idea of the q-group index is to have the same functionality as the q-gram index (i.e., retrieve all positions where a given q-gram occurs in constant time per position), but with a smaller memory footprint for large q. This is achieved by introducing additional layers in the data structure. In the following, we always consider a q-gram as its numeric representation g∈{0, …, 4^q−1}.

We divide all 4^q q-grams into groups of size w, where w is the GPU word size (typically w = 32). The q-gram with the number g is assigned to group number ⌊g/w⌋. Thus the ith group is the set G_i = {g|⌊g/w⌋ = i} of w consecutive q-grams according to their numeric order. The set of all q-groups is ${\mathcal{G}}_q≔\left\{G_0,G_1,\dots ,G_{\lceil \frac{4^q}{w}\rceil -1}\right\}$. We write g_ij for the jth q-gram in G_i.

For a given q and text T, the q-group index is a tuple of arrays ${\displaystyle {\mathcal{I}}_{T,q}≔\left(I,S,S^{\prime },O\right).}$ Array I consists of $\left|{\mathcal{G}}_q\right|$ words with w bits each (overall 4^q bits), and bit j of I[i] indicates whether g_ij occurs at all as a substring in the text, i.e., ${\displaystyle I{\left[i\right]}_j=\left\{\begin{array}{ll}1\hfill & \text{if }g\text{ is a substring of }T,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\right.}$

The array O corresponds to the position table P of a regular q-gram index: it is the concatenation of all occurrence positions of each q-gram in sorted numeric q-gram order. To find where the positions of a particular q-gram g begin in O, we first determine the group index i and the q-gram number j within the group, such that g = g_ij. With the bit pattern of I[i], we determine whether q_ij occurs in T. If not, there is nothing else to do. If yes, i.e., I[i]_j = 1, we determine the j′ such that bit j is the j′th one-bit in I[i].

The address array S contains, for each group i, an index into another address array S′, such that S′[S[i] + j′] is the starting index in O where the positions of g_ij can be found; see Fig. 1 for an illustration. All occurrence positions are now listed as ${\displaystyle O\left[S^{\prime }\left[S\left[i\right]+j^{\prime }\right]\right],\dots ,O\left[S^{\prime }\left[S\left[i\right]+j^{\prime }+1\right]-1\right].}$

Similar to a plain q-gram index, access is in constant time per position: For a given q-gram g, to determine (i, j) such that g = g_ij, we simply compute i = ⌊g/w⌋ and j = g−wi = jmodw. To compute j′, i.e., find how many one-bits occur up to bit j in I[i], we use the population count instruction with a bit mask: ${\displaystyle j^{\prime }=\text{Popcount}\left(I\left[i\right]\hspace{0.17em}\&\hspace{0.17em}\left(2^j-1\right)\right).}$ The population count Popcount(x) returns the number of 1-bits in x.

Size

To determine the worst case size of the q-group index, we note that both I and S consist of ⌈4^q/w⌉ words, S′ contains an index for each occurring q-gram and hence of min{4^q, |T|} words in the worst case, and O is a permutation of text positions consisting of |T| words. Thus the q-group index needs up to 2/w⋅4^q + min{4^q, |T|} + |T| words, compared to a conventional q-gram index with 4^q + |T| words. To evaluate how the worst case size of the q-group index behaves compared to the size of the q-gram index, we consider the ratio K between the possible q-grams and the text size, i.e., 4^q = K|T|.

If 4^q⪅|T|, the conventional q-gram index has a small advantage because each q-gram occurs (even multiple times). In fact, assuming w = 32, if K ≤ 1, the size ratio between q-group index and q-gram index is $\left(K/16+K+1\right)\left|T\right|/\left(\left(K+1\right)\left|T\right|\right)=1+\frac{K}{16\left(1+K\right)}$. For K = 1 (or |T| = 4^q), this means a small size disadvantage of 3% for the q-group index.

If, on the other hand, q becomes larger for fixed text size (such that q-grams become sparse), the q-group index saves memory, up to a factor of 16: For K > 1 the size ratio is (K/16 + 1 + 1)|T|/((1 + K)|T|) = (2 + K/16)/(1 + K), which tends to 1/16 for large K. The break-even point is reached for K = 16/15.

In practice, we use q = 16 because a bit-encoded DNA 16-mer just fills a 32-bit word and q = 16 offers reasonable error tolerance and high specificity. In this regime, with current GPU memory size, we use |T| = 10⁸ (processing 100 million nucleotides at a time), so the size ratio is K = 42.95, and the q-group index needs only 10% of the memory of the conventional q-gram index. Figure 2 shows the behaviour of the index size ratio depending on K.

Without noticable increase in access time, we can reduce the memory usage of the q-group index further by sampling every second (say, even) position of S and adding another population count and addition instruction to the code for odd positions.

The PEANUT algorithm

PEANUT uses the filtration plus validation approach: It first quickly detects exact matches of a given length q between each read and each reference (so-called q-gram matches) and only computes alignment scores where such matches are found. In total, the PEANUT algorithm for read mapping consists of the three steps (1) filtration, (2) validation and (3) postprocessing. The first two steps, filtration and validation, are handled on the GPU, while postprocessing is computed on the CPU.

The steps are conducted on a stream of reads. Reads are collected until buffers of configurable size are saturated. Then, any computation is done in parallel for all buffered reads (see Fig. 3).

In the filtration step, potential hits between the reference sequence and the reads are detected using the q-group index. Next, the potential hits are validated using a variant of Myers’ bit-parallel alignment algorithm (Myers, 1999; Hyyrö, 2003). The validated hits undergo a postprocessing that annotates them with a mapping quality and calculates the actual alignment. The postprocessed hits are streamed out in SAM format (Li et al., 2009). Because of memory constraints on the GPU, all steps are performed per chromosome instead of using the reference as a whole.

GPU resource usage

In order to maximize utilization of the GPU hardware, idle cores have to be avoided. The two most important reasons for idle cores are branching and memory latency (see “Designing for efficient GPU usage”). The latter can be hidden if the SM is able to execute a different warp while waiting on the transaction. The capability to do so can be measured with the occupancy, that is the fraction of active warps among the maximum number of warps on an SM. The more active warps exist on an SM, the higher is the chance that latency can be hidden by executing another warp. Figure 4 shows the occupancy patterns of the implemented OpenCL kernels, as measured with the NVIDIA™ CUDA command line profiler depending on the used thread block size (see “Designing for efficient GPU usage”). The thread block size influences the occupancy by limiting the number of potentially active warps and determining the amount of used registers and shared memory on the SM. Since the latter are limited, a bigger thread block size does not necessarily lead to a higher occupancy. As can be seen, the occupancy for all steps is high. For building of the q-group index (create_queries_index) and the filtration step (filter_reference), it even reaches 1.0 which illustrates the benefit of the q-group index being tailored towards the GPU architecture.

Sensitivity

Here, we strive to evaluate the sensitivity of PEANUT in terms of its ability (and hence that of the q-group index) to detect all alignments up to a given error rate. For this, we use the Rabema (Holtgrewe et al., 2011) benchmark that allows to compare mapping results based on a formalized framework. With Rabema, the genomic origin of a read itself defines how tight a mapping has to be using equivalence classes of ambiguous alignments. Further, it supports fair benchmarking of both all-mapping and best-mapping. First, 10,000 Illumina³ reads of length 100 were simulated using the read simulator Mason (Holtgrewe, 2010), which is part of the Seqan project (Döring et al., 2008), with the Saccharomyces cerevisiae genome,⁴ default parameters and error rates. Second, the simulated reads were mapped to the genome using RazerS 3 (Weese, Holtgrewe & Reinert, 2012) with full sensitivity and different error tolerances. In this configuration, RazerS 3 guarantees to report all alignments of a read up to a given error rate. For above notion of sensitivity, it would not be sufficient to consider the true origins of a read as they are known from the simulation with Mason. The mapped reads were used to generate gold standards for Rabema to test against. In “Comparison with other read mappers”, we assess the accuracy of PEANUT in comparison to other mappers on larger datasets.

We analyze the sensitivity of our algorithm using q-grams of length 16, because it is computationally optimal on the current GPU hardware. We configure PEANUT to provide all semi-global alignments of a read and leave all other parameters at their default values. Sensitivity is assessed by the Rabema measure “Normalized found intervals” (Holtgrewe et al., 2011) and all alignments of a read are considered (all-mode). With a percent identity threshold of 60 (see “Validation”) our algorithm provides 100% sensitivity for error rates below 5%, 99.86% sensitivity for error rates up to 10% and still 98.86% sensitivity with an error rate up to an unrealistically high 20%. With a stricter threshold of 80, PEANUT still reaches 98.81% sensitivity for the latter.

In general, we advise to set the percent identity threshold to be slightly more permissive than the expected error rate. This is because the replacement of N-characters in the reads and the reference with random bases can introduce additional mismatches. Above rates are far better than the worst case sensitivity that can be expected by applying the pigeonhole principle (i.e., with reads of length 100 and q-grams of length 16, we can expect to find at least one perfectly matching q-gram for all alignments with at most 6 errors), such that using 16-grams appears to be a reasonable default choice in practice.

Comparison with other read mappers

We compare the PEANUT algorithm with other state of the art read mapping algorithms in terms of run time and accuracy. The evaluation is conducted on 4 datasets:

1. 5 million simulated Illumina HiSeq 2000 reads;

2. 5 million real Illumina HiSeq 2000 reads from the human exome;

3. 10 million real paired-end Illumina HiSeq 2000 reads from the human exome;

4. 50 million real paired-end Illumina HiSeq 2000 reads from the whole human genome.

The simulated reads (dataset 1) were created from the ENSEMBL human reference genome⁵ version 37 with Mason (see “Sensitivity”). The read length is set to 100 and all other parameters of Mason are left at their default values, such that reads with a typical error profile and mutation rate are generated. The second and third datasets are generated from real paired-end exome sequencing reads⁶ (Martin et al., 2013) of length 100 obtained from a patient suffering from uveal melanoma (a cancer of the eye) sequenced with an Illumina HiSeq 2000 sequencer. Dataset 2 consists of the first 5 million forward reads. Dataset 3 consists of both the first 5 million forward and backward reads, i.e., 10 million reads in total. Dataset 4 is generated from real paired-end whole genome reads of length 200 obtained from an african male.⁷ The reads are part of the Illumina Platinum Genomes.⁸ The first 25 million forward and backward reads were chosen, i.e., 50 million reads in total.

The benchmark was conducted on an Intel Core i7-3770™ system (4 cores with hyperthreading, 3.4 GHz, 16 GB RAM) with an NVIDIA Geforce™ 780 GPU (12 SMs, 3GB RAM) and a 7200 rpm hard disk. We evaluated two modes of PEANUT. First, PEANUT was configured to find the best stratum of semi-global alignments (best-stratum mode) for each read. Second, PEANUT was configured to find all semi-global alignments (all mode) for each read. For comparison, we benchmarked the newest generation of BWA (BWA-MEM, version 0.7.5; Li (2013)), Bowtie 2 (version 2.0.2; Langmead & Salzberg (2012)), CUSHAW 3 (version 3.0.3; Liu, Popp & Schmidt (2014)), CUSHAW2-GPU (version 2.1.8; Liu & Schmidt (2014)), NextGenMap (version 0.4.11; Sedlazeck, Rescheneder & von Haeseler (2013)), RazerS 3 (version 3.2; Weese, Holtgrewe & Reinert (2012)) and MrFast (version 2.6.0.1; Alkan et al. (2009)). All tools were configured to use 8 threads (the reasonable choice in case of 4 cores with hyperthreading). For MrFast, which does not support multithreading directly, this was achieved by partitioning the input files containing the reads into 100 equally sized chunks and running 8 parallel instances of MrFast with the Unix command parallel (the time for merging the resulting output was not included into the run time). NextGenMap was used in GPU mode, such that it makes maximum use of the available hardware. All read mappers were configured to output alignments in SAM format (Li et al., 2009) directly to the hard disk.

We outline the reasons for excluding several available read mappers from the benchmarks. At the time of writing (07/2014), no working installations of SOAP3 (Liu et al., 2012), SOAP3-dp (Luo et al., 2013) and BarraCUDA (Klus et al., 2012) could be obtained. A binary compiled against the NVIDIA CUDA setup of the test system (NVIDIA CUDA 6) was not available for SOAP3 and SOAP3-dp. The compilation of SOAP3-dp-r177 and SOAP3-r146 fails on our test system. The Barracuda read mapper compiles but refuses to run on CUDA 6. Finally, read mappers specialized on RNA-Seq (Dobin et al., 2013; Trapnell, Pachter & Salzberg, 2009) were excluded, as this exceeds the scope of PEANUT. In principle, STAR (Dobin et al., 2013) could be applicable to DNA reads, and the authors claim significant speedups compared at least to other RNA-Seq focused read mappers. However, this comes at the cost of extensive memory usage by an uncompressed suffix array, which exceeds the capacity of the used test system.

In the following, we distinguish between all-mappers and best-mappers (PEANUT occurs in both categories, using the all-mode and the best-stratum mode). Best-mappers only strive to find the single origin of a read on the reference sequence of a single organism. All-mappers provide all alignments of a read down to a given error rate. Hence, all-mapping is computationally more intensive.

Run time performance is measured three times as the overall wall clock time for processing a dataset on the test system. Table 1 shows the run times for PEANUT and its competitors on all datasets. First, PEANUT in best-stratum mode outperforms all best-mappers (including the other GPU based mappers NextGenMap and CUSHAW2-GPU) on all datasets. On the biggest (and therefore most realistic) dataset, PEANUT is 2 times faster than the best competitor (BWA-MEM). Second, PEANUT in all-mode is 4 to 10 times faster than the all-mapper RazerS 3 and 4 to 6 times faster than the all-mapper MrFast. While Bowtie 2 provides an all-mode, too, it did not terminate in competetive time due to extensive memory requirements exceeding the capabilities of the test system.

The accuracy of the obtained alignments is assessed using Rabema (see “Validation”). For best-mappers, using the simulated dataset 1, we measure precision and recall with Rabema as defined by Siragusa, Weese & Reinert (2013): Recall is the fraction of reads correctly mapped to their original location. This location is known from the read simulation with Mason (see above). Precision is the fraction of correctly mapped reads among all reads that were mapped unambiguously (i.e., for which the mapper only provided exactly one hit). Figures 5 and 6 shows the results with increasing maximum error rate of the reads. Except when restricting to reads with zero errors (there, CUSHAW3 is slightly better), PEANUT slightly outperforms all other best-mappers in this benchmark. For best-mapping, the GCAT⁹ benchmark provides an alternative, less formal approach to measure the accuracy of a mapper by testing whether a mapped read is within 5 base pairs of its known origin in given datasets simulated with DWGSIM.¹⁰ For comparison, we provide GCAT results in Table 3 for the four datasets (1) 100bp-se-small-indel,¹¹ (2) 100bp-se-large-indel,¹² (3) 100bp-pe-small-indel¹³ and (4) 100bp-pe-large-indel,¹⁴ consisting of about 12 million 100 basepair reads each. The comparison is limited to those of our considered read mappers for which public results are available in GCAT (i.e., BWA-MEM, Bowtie 2, CUSHAW2 and CUSHAW3). The full reports are available via the GCAT web service under the URLs provided in the footnotes.

All-mappers are compared by their ability to find all alignments of a given error rate. Here, we again use the sensitivity provided by Rabema (see “Sensitivity”). This involves creating a gold standard with RazerS 3 configured to full sensitivity, which is computationally expensive. We therefore perform this on only 1,000 reads simulated with the same parameters as dataset 1. The gold standard is calculated for an error rate of at most 15%. Figure 7 shows that up to an error rate of 4% PEANUT provides a sensitivity of almost 100%, similar to RazerS 3 and nearly as good as MrFast. Beyond 5%, the sensitivity of PEANUT is supreme to RazerS 3 and MrFast. This is due to RazerS 3 and MrFast being restricted to low error rates per default to achieve acceptable performance. In summary, at default settings, PEANUT provides similar or even better sensitivity than RazerS 3 and MrFast while being 4–10 times faster (see above).

While PEANUT requires at least 2.5 GB of GPU memory for filtration and validation, it is not restricted to running on high-end GPU models like the Geforce 780 used above: Table 2 shows that an advantage can be maintained when benchmarking on a different test system with an Intel Core i7-2600 (3.4 GHz, 16 GB RAM) and a four year old NVIDIA™ Geforce 580 GPU. We see that the older system is about 9% slower and still faster than the best competitor in Table 1 on the newer test system. In Table 2, we further provide run times for the same system equipped with an NVIDIA Tesla K40c GPU with 12 GB RAM. While the Tesla is supreme in theory, here, it cannot outperform the Geforce. This is because a faster CPU (or more CPU cores) would be needed in order to also speed up the postprocessing which runs in parallel to the GPU based filtration and validation (see “The PEANUT algorithm”).

Evaluation of mapping qualities

Finally, we evaluate the mapping qualities provided by PEANUT. For each PHRED-scaled mapping quality q, the corresponding probability $10^{\frac{-q}{10}}$ equals the expected false positive rate. Here, true positives are the correctly identified true sampling positions of the reads, whereas false positives are reported mapping locations that may have the same alignment score but are not the true origins of a read. Figure 8 shows the measured and expected false positive rate (i.e., the fraction of false positives among all hits) at increasing mapping qualities for the hits reported by PEANUT in all-mode on dataset 1. As can be seen, the measured false positive rate rapidly decays when increasing the mapping quality from zero. A mapping quality above 10 (i.e., a probability of 0.1) already guarantees almost no false positives. Comparison with the dashed line, depicting the expected false positive rate for each mapping quality, shows that the mapping qualities provided by PEANUT are indeed conservative in the sense of underestimating the true probabilities. The parameters λ and C of the approximation (see “Postprocessing”) can be configured to different values to achieve a closer fit for a particular reference or dataset. Despite that, the defaults provide a reasonable way to distinguish between true and false positives without the computational overhead needed for calculating concrete alignments of many suboptimal hits (see “Postprocessing”).