RAMICS: trainable, high-speed and biologically relevant alignment of high-throughput sequencing reads to coding DNA

The challenge presented by high-throughput sequencing necessitates the development of novel tools for accurate alignment of reads to reference sequences. Current approaches focus on using heuristics to map reads quickly to large genomes, rather than generating highly accurate alignments in coding regions. Such approaches are, thus, unsuited for applications such as amplicon-based analysis and the realignment phase of exome sequencing and RNA-seq, where accurate and biologically relevant alignment of coding regions is critical. To facilitate such analyses, we have developed a novel tool, RAMICS, that is tailored to mapping large numbers of sequence reads to short lengths (<10 000 bp) of coding DNA. RAMICS utilizes profile hidden Markov models to discover the open reading frame of each sequence and aligns to the reference sequence in a biologically relevant manner, distinguishing between genuine codon-sized indels and frameshift mutations. This approach facilitates the generation of highly accurate alignments, accounting for the error biases of the sequencing machine used to generate reads, particularly at homopolymer regions. Performance improvements are gained through the use of graphics processing units, which increase the speed of mapping through parallelization. RAMICS substantially outperforms all other mapping approaches tested in terms of alignment quality while maintaining highly competitive speed performance.


INTRODUCTION
The issue of accurate pairwise sequence alignment is common to many fields in bioinformatics, whether as a core tool in fields such as reference-guided genome assembly (1)(2)(3)(4)(5)(6) or as the seed for the generation of a progressive multiple se-quence alignment (7)(8)(9). Ideally, analysis of a pairwise alignment should proceed in the knowledge that no inherent bias exists as a result of the alignment approach used. This is particularly challenging in the era of high-throughput sequencing, where every platform produces systematic errors (10)(11)(12)(13)(14)(15)) that should be considered in constructing an alignment.
The alignment of coding DNA, in particular, presents a unique challenge as it is critical that the final alignment takes into account the correct reading frame. Preservation of the reading frame ensures correct calling of gene structure (16)(17)(18) and SNPs, whether in, for example, the realignment phase of an exome sequencing pipeline (19), single-nucleotide polymorphism (SNP) calling from existing RNA-seq data (20,21) or amplicon-based analyses such as human immunodeficiency virus (HIV) drug resistance genotyping (22). When aligning coding DNA generated using high-throughput sequencing platforms, it is critical that codons present in the open reading frame remain intact and that codon-sized insertions and deletions are recognized and called correctly, as distinct from both genuine frameshifts and single indels created through sequencing error.
The landscape of pairwise alignment and reference mapping tools for high-throughput sequencing data is broad. Tools such as BOWTIE 2 (2) and BWA-MEM (3), while well-suited to mapping the location of query sequence reads within a complete reference genome, lack the subtle nuances required to correctly distinguish spurious indels from genuine codon-sized mutations in coding DNA. Other tools, such as MOSAIK (6) and SSAHA2 (5), perform full mapping and realignment using a Smith-Waterman approach, with a focus on correcting next-generation sequencing errors. However, basic Smith-Waterman alignment, even when taking into account quality scores as in MOSAIK and BWA-MEM, is not appropriate for reference mapping of coding DNA as it fails to maintain the intactness of codons. Finally, the Genome Analysis Toolkit (GATK) (23,24) performs the realignment step for tools such as BWA-MEM in, for example, the 1000 Genomes e106 Nucleic Acids Research, 2014 PAGE 2 OF 12 project pipeline (25), but does not consider the biological context (coding or non-coding) in its alignment.
One approach to ensure codons are maintained for downstream analysis, taken by the RevTrans program (26), is to initially translate query sequences into their corresponding amino acid sequences, align them to a translated reference sequence at the amino acid level and then 'back-align' the nucleotide sequences based on the amino acid alignment (26). While effective in some cases, this approach is useless in the presence of indels in the open reading frame, which results in mistranslation to amino acid space. Tools such as transAlign (27) repetitively translate, align, backtranslate and correct multiple sequence alignments, which to some extent addresses frameshift errors. This approach is, however, time-consuming, as in order to align robustly it requires a full amino acid multiple sequence alignment to create a back-translated DNA profile to which low-scoring sequences (assumed to be frameshifted) are compared. An alternative is to use an approach that undertakes pairwise alignment in 'codon space'--weighting the likelihood of each nucleotide substitution within the context of the codon to which it belongs. This, thus, facilitates the identification of indels, whether genuine frameshift mutations or the result of sequencing error.
Initial models of coding DNA such as those used by Pair-Wise (16) and EST Genome (28) were designed specifically for gene prediction, comparing sequences against all six translated frames simultaneously using a classical dynamic programming approach to discover exons. However, this approach was found to be time-consuming and overly specific (29). The HMMER tool (30) took a different approach, using profile hidden Markov models to align and compare protein sequences. GeneWise and GenomeWise (29) use a hidden Markov model approach to allow for the alignment of a nucleotide sequence, potentially including frameshifts, to an existing protein sequence, and thus represent the closest approximation to a codon-codon alignment tool. However, because the alignment is to a protein sequence, it is impossible to penalize synonymous mutations. Protein mismatches are all scored equally, as are single-nucleotide insertions resulting from sequencing error. This prevents the use of GeneWise and GenomeWise to correctly identify and account for homopolymer errors in next-generation sequencing platforms.
Here, we present a novel reference mapper (RAMICS: rapid amplicon mapping in codon space) that uses a hidden Markov model approach to align large numbers of coding DNA sequences to a coding region of a reference sequence in an accurate and biologically relevant manner, accounting for the inherent error biases of the sequencing platform used. RAMICS is developed to be executed on graphics processing units (GPUs) and, thus, provides significant speed improvements over other, similarly complex, mapping tools. RAMICS is, thus, ideal for any situation where a large number of sequence reads from a known, coding location in a genome must be mapped to a coding reference sequence.

The open reading frame
RAMICS is designed specifically to align coding DNA. The tool first discerns the optimal open reading frame of the reference sequence by selecting a frame that would minimize the number of stop codons. RAMICS then performs a pairwise alignment to each query sequence in 'codon space' using a hidden Markov model (31) (Figure 1 and Supplementary information). Single-and double-nucleotide indels are identified and accounted for, thereby ensuring that biologically relevant complete codons are kept intact. This enables a correct and relevant comparison of a query sequence to the reference, as the concept of a codon-sized insertion that begins within one codon and ends within another (as would be output by most conventional alignment approaches) is structurally meaningless in a biological context ( Figure 2A).
RAMICS uses models of codon substitution (32-34) to give synonymous codon substitutions higher transition probabilities as well as favoring more likely nonsynonymous substitutions. RAMICS can also be set to use a simple non-coding substitution model to align non-coding DNA to a non-coding reference sequence. However, RAM-ICS cannot concurrently align sequence reads to coding and non-coding regions of the same reference sequence (e.g. a reference sequence containing both introns and exons), and, as such, is unsuitable for use as the first-pass mapping tool in full-genome association studies. It can, however, be subsequently used for targeted realignment of coding regions as a second-pass mapping tool.

The sequencing platform
RAMICS uses reported average error rates for given sequencing platforms (10,12,13,15) to set the transition probabilities to single-and double-nucleotide indel states in the . Sequences X and Y differ from each other by a three-nucleotide insertion. Conventional nucleotidebased approaches generate a statistically optimal pairwise alignment matching as many nucleotides as possible. RAMICS, on the other hand, considers the biological importance of codons and generates a biologically relevant alignment preserving the codon structure of the coding DNA being aligned. (B) 'Correctly accounting for homopolymer errors'. Comparing alignments of an HIV pol gene sequence containing a homopolymer error (additional adenosine, green arrow) to a reference sequence shows that RAMICS correctly identifies and accounts for this error. The codon of interest in the reference sequence is marked with a black box. Conventional nucleotide alignment approaches would call this codon as AAA (flagging the G as an insertion) which encodes for lysine (sensitive form), while RAMICS would call the codon as AGA (arginine, the resistant form) correctly identifying and discounting the extra adenosine introduced as a result of the sequencing error at the homopolymer region. profile hidden Markov model. In addition, the model is more likely to place a single-nucleotide indel on either side of a homopolymer region when the sequencing platform is set to Roche/454 (35) or Ion Torrent Personal Genome Machine (36), as homopolymer miscalls are a common sequencing error on these platforms (13).
To illustrate the necessity of considering the sequencing platform used to generate the sequence data, we take an example from HIV-1 antiretroviral drug resistance testing. The drug resistance profile of a sequence is determined by characterizing mutations at various amino acid positions throughout the polymerase gene mapped in relation to the HXB2 HIV reference sequence (51). Drug resistant mutations include the mutation of lysine (K) to asparagine (N) at codon 103 (K103N), which confers drug resistance to certain non-nucleoside reverse transcriptase inhibitors (NNR-TIs) (22).
The K103N codon is located in a homopolymer region and some high-throughput sequencing platforms such as Roche/454 (35) and Ion Torrent PGM (32), as well as to a lesser extent the Illumina platform (37), are prone to increased error rates in regions containing homopolymers, as they cannot successfully identify the correct number of nucleotides in these regions (12,13). Incorrect assignment of a e106 Nucleic Acids Research, 2014 PAGE 4 OF 12 homopolymer indel during the mapping step would result in a viral sequence being falsely called as susceptible when it is actually resistant (54).
A standard Needleman-Wunsch alignment to the HXB2 reference sequence of a sequence containing both the K103N mutation and a spurious extra adenosine directly following homopolymer region would result in the virus in question being falsely classified as susceptible to NNRTIs such as efavirenz ( Figure 2B). Mapping of the same data using RAMICS, however, weighs the possibility of a false tyrosine (T) insertion in the center of a homopolymer region against a false adenosine (A) insertion at the end of the region, assigns the latter a higher probability and, thus, correctly calls the presence of the K103M mutation in spite of the hompolymer error ( Figure 2B). Thus, a major strength of RAMICS is in its application to situations where the correctness of alignment in a coding region is crucial, such as alignment optimization in RNA-seq and exome sequencing (19)(20)(21)38,39).

Training
It can be a common case in reference mapping, particularly for amplicon-based analyses, that a set of query sequence reads will be more closely related to each other than they will to a given reference sequence--this is true in reference-guided assembly (40), drug-resistance genotyping (41) and applications in metagenomic function prediction (42). However, a comparison to the reference sequence is still required for these applications. RAMICS circumvents any potential bias introduced through a distantly related reference sequence by allowing the use of the reference sequence as a guide for the initial mapping and then iteratively training the underlying hidden Markov model to the resulting pairwise alignments. The result is a set of emission and transition probabilities that will correctly align outlier sequences (potentially with many sequencing errors) in the query set, based on the substitutions in sequences more similar to the reference. This approach generates a global reference alignment that approaches a reference-seeded multiple sequence alignment, and which can be used as such.

Output formats
RAMICS produces results in various formats with the choice of output format very much driven by the end-users requirements:

Evaluation of RAMICS
RAMICS is designed to map large numbers of highthroughput sequencing reads from a coding DNA region to a similar reference sequence, both quickly and accurately. RAMICS also has potential as an alignment refiner for the coding regions of whole-genome mapping or as the seed pairwise aligner for a profile HMM multiple sequence aligner in the style of the protein sequence aligner Clustal Omega (44). Both as a standalone amplicon mapper and as a component of these larger toolchains, it can align large numbers of reads quickly as well as aligning reads accurately in regions that have similar biological features but contain many synonymous or functionally similar substitutions. Although RAMICS has the ability to remove frameshifts in reads, under the assumption that they are sequencing errors, we did not remove any indels before evaluating the alignments.
We first compared the quality of alignments with RAM-ICS in non-coding mode, coding mode without training and coding mode with training, respectively. We then compared RAMICS for quality and speed against Bowtie-2 (2), BWA-MEM (3), BWA-MEM post-processed by the Indel-Realigner from the GATK (23), MOSAIK (6), the basic Needleman-Wunsch aligners EMBOSS 'needle' and 'water' (45) ('needle' is used where a global alignment is more appropriate, 'water' for local alignment), NextGenMap (46), SSAHA-2 (5) and SHRiMP2 (1). These tools represent a mix of global, local and 'glocal' alignment techniques (47), as for the mapping of short reads to a relatively short amplicon all methods are potentially valid.

Quality performance
The RABEMA benchmark. Several benchmarks for the mapping of simulated next-generation sequencing data exist (48,49). The SEAL benchmark (49) is not easily extensible to new alignment tools, so we first performed a simple benchmark of RAMICS' relative ability to correctly dealign spurious indels using the Rabema benchmark (48). Rabema uses the MASON read simulator (50) to simulate artificial short reads containing sequencing errors from a user-provided reference sequence and to output a 'gold standard' SAM files containing each of these reads mapped to the reference sequence. The alignments generated by a mapping tool are then compared to the 'gold standard' alignments.

PAGE 5 OF 12
Nucleic Acids Research, 2014 e106 The tool then bins reads by the number of errors found, ranging from perfect reads to highly erroneous outlier reads containing 4% or more errors. This allows us to test alignment performance not only of reads exhibiting the average error rates reported by manufacturers but also of rare outlier reads that exhibit a higher error rate than the reported average. The Rabema output gives the percentage of correctly mapped reads as a function of percentage sequencing error in the read. We ran each benchmark 10 times with 10 different random seeds to ensure statistical significance.
The AFAB benchmark. While Rabema is a good benchmark of a mapping tool's ability to detect sequencing errors such as spurious indels, it only determines whether an alignment is correct or incorrect and, thus, does not allow quantification of how incorrect the alignment is from the 'gold standard'. Further, because all of the reads generated by Rabema originate from the same reference sequence, the levels of diversity of the query sequences being mapped are lower than what would be generated in many genuine sequencing experiments. Thus, we developed a flexible alignment benchmark (AFAB) that could evaluate not only whether a correct alignment was found, but also how similar this alignment was to a gold standard alignment, in the case where very divergent sequences made perfect mapping impossible. This was achieved using a combination of tools.
We began by selecting a 'gold standard' alignment of protein coding DNA which was the highly manually curated HIV-1 group M subtype reference alignment of the HIV envelope gene from the LANL HIV reference sequence database (51). This alignment was chosen for two reasons: (i) it encompasses the broad spectrum of diversity present within the HIV-1 group M subtypes (pairwise genetic diversity measured with the generalized time-reversible substitution model (52) ranged from 0.15 to 0.24 substitutions per site, median 0.218) and (ii) the envelope gene comprises both conserved and variable regions. This ensured that the benchmarking was performed on a complex example of 'real-life' data to which mapping tools are applied.
Using one representative sequence per subtype we extracted all possible pairwise alignments resulting in 55 unique alignments. For each alignment we used ART (48) to generate artificial Roche/454 and Illumina reads (1000 reads per pairwise alignment) from the second sequence of the pair. Further, ART generated SAM files containing the correct alignment of each simulated read to the reference sequence (the first sequence), which were used as the 'gold standard' for evaluation of the accuracy of the various mapping tools. The choice of ART has the added benefit of creating greater test coverage, ensuring that RAMICS is not just incidentally biased to MASON's model of read simulation. For both the Roche/454 and Illumina data the qscore tool (8) was used to evaluate each of the mapping tools, giving the Q Score and Cline shift similarity scores of each tool's generated alignment to the gold standard.
We then ranked the Q Score and Cline shift for each tool that chose to map each read, without penalizing tools that chose not to map the read as Rabema does. When multiple tools achieved the same score (most often 1.0 in the case of trivially mappable reads), we assigned each tool the av-erage score they would receive, to avoid artificially inflating the scores of tools that chose to map only trivially mappable reads. This follows the ranking scheme used in the Friedman test (53) and removes the assumption of normality across easy and hard to map reads. We then report the average score obtained across all reads mapped, with a score closer to 1 indicating a consistently higher ranking of a mapping approach. We repeated each benchmark five times with five different random seeds to ensure statistical significance.

Speed performance
We benchmarked the wall clock time of RAMICS, and other tools, in mapping amplicon-style data. We mapped one million Roche GS Junior FASTQ reads (average length 237 base pairs) from the HIV reverse transcriptase (RT) gene to a short, 312 base pair region of the HXB2 HIV-1 reference sequence (54). All benchmarks were performed on an Intel R Core TM i5-3210M CPU @ 2.50 GHz and an NVIDIA R GeForce R GT 525M graphics card.

Evaluation of the RAMICS algorithm
Initially we evaluated the various components of the RAM-ICS algorithm to explore whether the addition of complexity improves the quality of the resulting alignment. We mapped both the Roche/454 and Illumina versions of the AFAB benchmarking data using the simplest incarnation of RAMICS (without codon alignment or training), RAM-ICS with codon alignment and the full RAMICS algorithm (with both the codon-based alignment and training parameters invoked). The addition of codon-based alignment increased the number of perfectly aligned sequence reads by 7.9 and 10.5% for the Roche/454 and Illumina data, respectively ( Figure 3). Further, including three rounds of training of the HMM to the query data provided marginal improvements in the number of perfectly mapped reads by 1.5% (Roche/454) and 0.9% (Illumina). The number of discarded reads remained mostly unchanged regardless of the parameters used ( Figure 3) showing that while alignment accuracy increases with algorithmic complexity, RAMICS is consistently capable of distinguishing between mappable and unmappable reads. The lesser impact of the training parameter can be explained by its intended use case: the proper alignment of extreme outlier reads with large numbers of sequencing errors. Thus, in most instances we suggest that RAMICS should be used with both the codonbased alignment and training parameters invoked. For all further benchmarking analysis we used the full RAMICS algorithm including both the codon-based alignment and training parameters.

Comparison of RAMICS to other tools
The use case for which RAMICS is designed is rather unlike that for which most other mapping tools are tailored, because of its focus on coding DNA. Despite this, many studies have used various mapping tools to analyze sequence reads in ways similar to that which RAMICS was designed for (55-58) and, thus, we present comparisons of RAMICS to a number of commonly used mapping tools on the basis of alignment quality and speed.

Comparison on the basis of alignment quality
The Rabema benchmark. Initially we compared RAM-ICS' mapping quality against other tools using the Rabema benchmark (39) for both Roche/454 and Illumina data. Reads were simulated from a single reference sequence using the reported error rates of the sequencing platforms as encoded by the Mason read simulator (50). These resulting sequence reads are then binned by their percentage difference from their seed sequence. For example, if the seed sequence used to generate a read was 100 nucleotides long and the resulting read had three differences relative to its seed, then that read would be added into the 3% error bin. This allowed the evaluation of the mapping tools using realistic error rates for both typical and outlier sequences. As would be expected, the observed error rate for vast majority of sequence reads for the Illumina data set fell around 1% with 88% of sequences exhibiting an error rate of ≤1%. For the 454 data, 49% of the sequence reads exhibited an error rate of ≤1% while the error rate of the remaining reads was in excess of 1%. The higher-than-expected number of sequences with an error rate of >1% most likely results from the fact that the single reference sequence used to seed the simulations contained a number of homopolymer regions, which are known to increase the overall observed error rate of the Roche/454 platform (12,13).
For the Roche/454 benchmark, the performance of all of the tools tested was comparable for reads containing zero errors, ranging from 100% to 97.9% of reads correctly mapped for RAMICS and SSAHA, respectively ( Figure  4A). Such high performance is to be expected, as these sequence reads are identical to the reference sequence. As the percentage of erroneous nucleotides in a read increased, the accuracy of each of the mapping tools, with the exception of RAMICS, reduced ( Figure 4A). SHRiMP (1) was excluded from the Roche/454 benchmark as it discarded over 70% of the reads containing 4% erroneous nucleotides with the recommended settings.
A similar result was observed for the Illumina benchmarking data with all methods correctly mapping 100% of reads at zero error rate and the accuracy of the various approaches, with the exception of RAMICS, reducing as the percentage error rate increased ( Figure 4B). For the Illumina benchmark data we observed that compared to the Roche/454 data, the rate of reduction of accuracy is much sharper as the error rate increases. We suggest this effect is due to the shorter lengths of Illumina reads, which leave fewer correct nucleotides on either side of an error with which to anchor an alignment.
RAMICS outperforms all of the other approaches for both the Roche/454 and Illumina benchmarking data sets. This is because RAMICS renders sequencing error in reads when aligning to a reference sequence almost irrelevant, even in outlier reads with a large number of errors, by correctly identifying likely sequencing errors. MOSAIK was the next best performing tool (albeit with a reduction of 8% in accuracy at an error rate of 4% in the Illumina data), while SHRiMP and SSAHA2, both of which conservatively choose not to align reads, quickly deteriorate in their mapping ability with increasing error rate (Figure 4). The GATK IndelRealigner (23) made no statistically significant difference to the test results, so we have removed it from Figure 4 for clarity.
The AFAB benchmark. The Rabema benchmark clearly shows that RAMICS outperforms all other mapping approaches on the outlier sequences with a higher-thanaverage error rate. However, these data are not a true reflection of the type of data for which RAMICS has been optimized, namely mapping a diverse set of sequence reads to a related reference sequence. Rabema data only differ from the reference sequence through the insertion of sequencing errors and, thus, we developed a benchmarking tool to enable the assessment of mapping tools at aligning coding sequence reads containing sequencing errors that are diverse from the reference sequence. In contrast to the Rabema benchmarking tool, this approach enabled the use of two robust methods for evaluating relative alignment accuracy enabling the identification of perfectly aligned reads and misaligned reads. Reads discarded during mapping were not penalized in the final ranking thereby allowing tools to discard ambiguous reads without incurring unnecessary bias in the benchmarking process.
Extraction of all possible pairwise alignments of subtypes A through G from the 'gold standard' HIV alignment resulted in 55 sequence pairs, each of which were evaluated using this approach. Firstly, we considered the percentage of reads that each mapping approach was capable of aligning regardless of the quality of the alignment. Needle (45) does not possess the capability to discard any reads, resulting in 100% of reads being mapped for both the Illumina and Roche/454 data ( Figure 5). The percentage of reads retained by the other approaches ranges from 38% (SHRiMP) to 94% (MOSAIK) for the Roche/454 data and 20% (BWA-MEM) and 62% (RAMICS) for the Illumina data ( Figure  5).
As expected, because these data contained sequencing errors and were diverse from the reference sequence, the percentage of perfectly mapped reads was lower for all approaches than was observed for the Rabema benchmarking data. For the Illumina data, RAMICS significantly outperforms all of the other mapping approaches by aligning 62% of the total number of reads perfectly, with the next best approach (needle) mapping 48% of all reads correctly (Figure 5). For the Roche/454 data RAMICS performed best, correctly mapping 36% of reads perfectly, while the performance of the other methods ranges between 19% (MO-SAIK) and 34% (SSAHA2). The GATK IndelRealigner made no difference to the outcome for BWA-MEM, so we have removed standalone BWA-MEM from Figure 5 for clarity.
It is also interesting to consider the read discarding strategies employed by the various approaches. Ideally, a mapping tool would discard all reads that it could not align perfectly, while still aligning as many reads as possible. Considering that all reads in this benchmark data set are simulated from a gold standard alignment, Bowtie and BWA-MEM's strategy of discarding almost 80% of the reads leads to a high rate of perfect alignment, but an unacceptable loss of coverage depth in situations where many variants are being aligned or deep coverage is necessary.
For all of the reads that were aligned by one or more mapping tools we used two approaches for evaluating alignment quality, the Q Score (59) and Cline shift score (60), to rank each of the mapping tools for each of the 5500 pairs of sequences reads that were tested.
The perfect mapping tool would rank first or joint first compared to all other tools for every single comparison. When two or more tools received equal ranking ties were broken according to the Friedman test (53) by assigning, for example, a ranking of 1.5 to each tool if two tools were equally first, preventing a bias toward tools that aligned more reads. The average Q Score ranks for the methods ranged from 2.512 to 4.846 and 2.978 to 4.562 for the Illumina and Roche/454 data, respectively ( Table 1). The average Cline shift ranks range from 2.522 to 4.565 and 3.039 to 4.817 for the Illumina and Roche/454 data, respectively. For both benchmarking data sets RAMICS was the highest ranked mapping approach tested, with BWA-MEM performing worst over both data sets ( Table 1). The addition of the GATK (23) to the BWA-MEM toolchain, the approach used in the 1000 Genomes Project (25), had a negligible effect on the quality of alignment for this use case.
Thus, for both benchmarking approaches used across both Illumina and Roche/454 data, RAMICS outperforms all other methods tested here in terms of alignment quality. We must reiterate, however, that all approaches here are being tested against data that have been simulated to test the use case for which RAMICS has been specifically designed. While RAMICS outperforms all of the mapping approaches tested here, it is context based (focused on coding DNA that is divergent from the reference sequence) and, thus, it should not be assumed that RAMICS outperforms all methods in all scenarios, e.g. whole-genome mapping.

Speed performance
While alignment accuracy for mapping tools is essential, the scale of the data sets being generated by high-throughput sequencing approaches means that mapping speed is also a critical factor for an optimal method. With this in mind we have designed RAMICS to harness the computing power of graphic processing units (GPUs) for optimal speed. We have also developed a CPU-based version for researchers who do not have access to GPU hardware. Here, we have compared the mapping speeds of both the CPU and GPU versions of RAMICS against the other mapping approaches in mapping one million Roche/454 reads from the HIV RT gene to a reference sequence. We found that, on average, the speeds for all approaches scale linearly as the data set size increases (data not shown) and find that in most instances NextGen-Map is the fastest of the tools tested while MOSAIK is the slowest (Figure 6). Not surprisingly, we find that the GPU version of RAMICS outperforms the CPU version with as much as a 3-fold speedup. Nonetheless we find that RAM-ICS performs in line with most of the read mappers tested.
The majority of the other mapping tools are built to align reads quickly to a large genome using various heuristics Illumina reads that were misaligned (i.e. not discarded by the mapping tool but with an incorrect alignment when compared to the 'gold standard') as well as the percentage of reads that resulted in a perfect alignment when compared with the 'gold standard' alignment. We have excluded BWA-MEM without GATK as it performed identically to BWA-MEM + GATK. For both test data sets, the highest ranked mapping approach is marked in bold.
e106 Nucleic Acids Research, 2014 PAGE 10 OF 12 Figure 6. Comparing mapping tools in terms of both speed (reads mapped per second) and accuracy (average Q Score of the Roche/454 reads generated using the AFAB benchmarking tool). The error bars on each data point reflect the range for both the speed (y-axis) and accuracy (x-axis). RAMICS' speed performance remains competitive with some of the fastest mappers, while its quality performance outstrips them significantly. and, thus, much of the overhead necessary to hash a large genome is redundant in these analyses. RAMICS, on the other hand, is designed to perform a full dynamic programming alignment to a shorter reference sequence quickly and accurately.
Comparing the RAMICS GPU version to the EMBOSS tool needle (45), which performs a full dynamic programming alignment, shows that RAMICS is roughly twice as fast ( Figure 6). Further, the efficiency of RAMICS is evident in that the codon-based approach implemented means it must perform eight comparisons to needle's three, and makes extensive use of floating point arithmetic to achieve the fine-grained distinctions required for its more complicated alignment model. Indeed, RAMICS non-coding, which performs a Needleman-Wunsch style alignment but takes into account homopolymer errors, is itself twice as fast as the full RAMICS algorithm (data not shown).
There should, however, be no trade-off in quality for speed or vice versa and, thus, the optimal mapping approach is one that can map the greatest number of reads with the highest accuracy in the shortest time. By evaluating each mapping approach in terms of speed and accuracy together (Figure 6), we find that the GPU-based version of the full RAMICS codon-based algorithm outperforms all other approaches suggesting that, for amplicon-style analyses at least, many of the other approaches exhibit such a trade-off.
To our knowledge, RAMICS is the first mapping tool that undertakes accurate mapping of coding DNA to a reference sequence at speed, while concurrently identifying and correcting sequencing platform-induced errors. RAM-ICS' accuracy in generating biologically correct alignments of coding DNA means that it can be applied wherever the generation of high-quality alignments for sensitive detection of SNPs and novel variants are required.