BACKGROUND
high-throughput sequencing  technology has recently shown a rapid and impressive development and this has led to the production of gigabases of sequence in a few hours for only a fraction of the former cost  <cit> .
hts has produced an explosion of knowledge in genetics and genomics thanks to the development of specific applications such as genome re-sequencing .
this technological evolution was paralleled by the development of new algorithms to deal with the quantity and the quality of reads produced.
a fundamental analysis steps in re-sequencing approaches is the mapping of the reads onto a reference genome.
this step, which involves the accurate positioning of reads onto a reference genome sequence, is highly important because it determines the global quality of downstream analyses.
the algorithms used for this step are called mappers.
mappers have to be sensitive and accurate and, if possible, fast and not too computationally demanding.
they should be able to find the true position of each read on a reference genome and ideally distinguish between technical sequencing errors and natural genetic variations.
in recent years many mappers have been developed and distributed .
two studies  <cit>  have classified mappers using a wide variety of features that include: the type of data, their application, the sequencing platform, the read length, the allowed error rate, parallel implementation, the ability to deal with multi-mapped reads , the input and output formats, and the available parameters.
mappers have multiplied and so has the range of possible settings.
hence, the growing difficulty in selecting a mapper has been raised in recent studies aimed at evaluating mapper performances through a multiplicity of comparison criteria.
some of these studies have focused on mapper sensitivity   <cit> .
schbath et al.
studied the ability of mappers to identify unique versus multi-mapped reads using a well-controlled benchmark containing reads with exactly three mismatches  <cit> .
hatem et al.
introduced a benchmarking suite to analyze mapping tools  <cit> , which consists of tests that cover input properties and algorithmic features.
in addition to the difficulty in determining evaluation criteria, choosing an appropriate evaluation method, i.e.
how to compare mappers according to the evaluation criteria, and using the appropriate metrics, are also problematical.
using real datasets to evaluate mapper performances allows only a rough assessment and classification of mappers by comparing the percentage of mapped reads, but does not reveal the actual accuracy of mappers.
attempts have been made to avoid this pitfall using simulated datasets in which the original read positions are known.
another difficulty lies in the accurate definition of what a correctly mapped read is.
the basic definition is to consider a read as correctly mapped if the original location is retrieved  <cit> .
ruffalo et al.
broadened this definition by adding a condition on the quality score, which had to be superior to a given threshold  <cit> .
in a more recent paper  <cit> , a new definition was introduced in which a read was considered to be correctly mapped if the mapping criteria were not violated, i.e.
contained less errors than the threshold parameter set by the user.
using simulated data allows numerical values to be obtained and compared between a set of mappers.
however, simulated data do not have the same characteristics as real data, even when an error model based on real data is used.
real hts data present biases  <cit>  that can be very difficult to simulate.
additionally, the current definition of the mapping correctness based only on the original start location presents some weaknesses: a read can have several correct positions on the reference sequence and sequencing errors or true genetic variations can lead to a better alignment in a genome position different from the original one.
holtgrewe et al.
introduced the interval definition, rather than the genome position, to describe a read mapping  <cit>  and used a full-sensitivity algorithm to identify all possible matching intervals within a given error rate range for each read.
this method has been implemented in rabema , a tool that evaluates the result of arbitrary read mappers that support the sam output format with real and simulated datasets.
our analysis of the published literature on mapper evaluation led us to conclude that for a complete and robust comparison of mappers, real and simulated datasets should be used.
using real datasets avoids simulation biases and gives a real picture of mapper behavior, whereas simulated datasets are benchmarks from which all parameters can be controlled.
additionally, a sound, more complete definition of what constitutes a correctly mapped read needs to be considered .
in all the previous studies, mapper performance was evaluated using large eukaryotic genomes  and, for the most part, short illumina or illumina-like reads data were used, except in  <cit>  where  <dig> datasets were evaluated with a reduced number of mappers and metrics.
the type of sequencing errors and their rate is inherent to the sequencing technology and more precisely to the nucleotide elongation detection methods used.
for example, life technologies sequencing by oligonucleotide ligation and detection  technology showed a strong bias in its coverage of repetitive elements  <cit> , whereas the illumina reversible dye-terminator sequencing technology  mainly caused substitutions  <cit> .
pyrosequencing on solid support  and ion semiconductor sequencing technology  produced indel errors associated with homopolymer-regions  <cit> .
in the published evaluations, the criteria that were tested and the default parameters of the mappers were usually chosen to address or deal with substitution-type errors and are, therefore, less informative for mapping the reads from new technologies like the ion torrent platform.
furthermore, the analysis of small microbial genomes compared with the analysis of large eukaryotic genomes poses other challenges because microbial genomes contain a wide range of gc content, which is sometimes extreme.
very high or very low gc content means that there is a high probability of encountering homopolymers in a genome sequence and this is known to be a specific problem for pyrosequencing and ion semiconductor sequencers.
a recent development in the hts technologies has made available benchtop sequencers targeted at the quick and inexpensive sequencing of small to moderate-sized genomes, mainly bacteria, viruses, fungi, and parasites.
small microbial genome sequences could be considered to present a simpler, less demanding mapping process compared with the mapping process for larger eukaryotic genomes.
however, this is only partially true because the characteristics of small microbial genomes are not the same as those of eukaryotic genomes.
the questions of interest are also usually different and, consequently, the expected mapping quality criteria are not exactly the same.
whole genome sequencing or re-sequencing is an important application in the new field of microorganism characterization using hts.
for instance, clinical diagnosis and the epidemiological study of microbial strain circulation will be profoundly remodeled in the near future by the use of hts, which should, very soon, be used as a characterization approach for pathogens and which will probably slowly replace the present pcr and biochemical based characterization methods  <cit> .
in this particular context the re-sequencing applications and derived analyses are in the front-line of research and development.
the focus includes the sequencing of the entire length of a microbial genome and the analysis of obtained reads by mapping them onto one or several reference strains to identify potential relevant changes in the studied genome.
the aim is to accurately identify the gain or loss in genetic elements  as well as small changes  to predict a potential new phenotype or a derived new pathogenicity profile.
this requirement poses several challenges, the most important of which is the necessity to distinguish true genetic variations from sequencing errors.
in this paper, we focus on the evaluation of mappers in the context of whole genome sequencing or re-sequencing for small microbial, mainly bacterial genomes.
we tested  <dig> mappers, mostly using their default settings to be in the general context of non-expert users.
we selected four criteria to match this context:  computational resource and time requirements,  robustness of mapping through the evaluation of precision, recall and f-measure,  ability to report positions for reads in repetitive regions, and  ability to retrieve true genetic variation positions.
to evaluate a mapper’s robustness on simulated datasets, we introduced a new definition of a correctly mapped read.
in addition to the original start position  that was used in most previous studies, the end position as well as the numbers of insertions, deletions, and substitutions in the alignment were also used to classify the mapping of a read as correct.
this definition is more stringent than the previous ones because it implies that it is a full-length read alignment and that the error count is correct.
indeed, sequencing errors can mean that the original location of a read is not necessarily the best alignment location.
using mappers tuned to report all possible hits  and to accept a higher error rate than the error rate introduced in simulated reads, it should be possible to retrieve the original location in addition to potential equivalent or better hits.
with the new definition of a correctly mapped read used in this study, we ensured that the mapper was able to retrieve the expected original alignment despite inevitable sequencing errors in the reads, thereby allowing a true evaluation of the mapper’s robustness.
the analysis was applied to data generated by the ion torrent personal genome machine , a newly arrived technology dedicated mainly to small genome sequencing, for which mapper performances have not yet been evaluated.
reads from real datasets and artificially simulated reads were used.
simulated reads were generated using a new customizable read simulator, curesim, which can generate reads of user-determined lengths with insertions, deletions, and substitutions introduced at a controlled rate and with an adjustable error distribution along the read.
curesim and curesimeval, a script that can be used to evaluate mapping quality, were developed in java to run on all operating systems  and are freely available at http://www.pegase-biosciences.com/tools/curesim/.
we have shown that in microbial genome sequencing, some mappers, such as segemehl, present higher robustness than others, especially when the number of sequencing errors was high.
other mappers are more robust for other applications that demand other quality criteria.
for example, bwasw, shrimp <dig>  smalt, ssaha <dig> and tmap, might perform particularly well for sequencing focused on rare variant discovery because they show a robust discrimination of variations.
smalt can localize most of the positions of reads located in repeated regions.
some mappers, such as novoalign, smalt and srmapper, needed very small memory resources , while snap was very fast and required only about two minutes to process the bigger datasets used in this study.
these results emphasize the observation that mapper choice is application dependent and users should carefully consider the targeted aim before choosing a mapper.
the evaluation approach presented here, together with the developed tools  can be considered as a general method to evaluate existing or in-development mappers and could prove interesting in the evaluation of the performances of mappers for the coming third generation of sequencers that may have yet another type and rate of errors.
RESULTS
computational resource requirement and time measurement
all mapping processes involve the alignment of millions of reads onto a reference sequence.
this is true even for small genome sequencing projects where the small size of the reference sequence is generally compensated by the multiplicity of samples to be analyzed.
in clinical microbiology, the time and the computational resources required for the analysis are critical; therefore,  <dig> these factors also need to be evaluated for the different mappers.
all the mappers tested were run with  <dig> threads  and the memory consumption and runtime were recorded for three different ion torrent datasets rd_ <dig>  rd_ <dig>  and rd_ <dig>  these three datasets contain real single-reads with different mean sizes and are described in table  <dig>  the reference genome used was escherichia coli str.
k- <dig> substr.
dh10b  for the rd_ <dig> and rd_ <dig> datasets and escherichia coli str.
k- <dig> substr.
mg <dig>  for the rd_ <dig> dataset.
figure  <dig> shows the memory consumption for each mapper for the real datasets when the indexing and mapping steps were considered together.
novoalign, smalt, and srmapper needed very low memory resources .
it should be noted that srmapper was developed to run on a computer with  <dig> gb of ram for genomes the size of the human genome, but, in such a case, it can be run only in ‘all-best’ mode and does not allow indels in the mapping.
the novoalign version used in this study was the free academic version that has not been implemented in parallel.
a second group comprising bowtie <dig>  mosaik, and segemehl, needed less than  <dig> gb of ram, while a third group, bwa, bwasw, and tmap, needed less than  <dig> gb of memory.
bwa had peak memory usage of  <dig> gb for rd_ <dig> and of more than  <dig> gb for the rd_ <dig> dataset.
bwa was developed to map short reads of up to  <dig> bases, which may explain the high peak usage for 400-base reads.
shrimp <dig>  snap, and ssaha <dig> required more ram  and ssaha <dig> needed about  <dig> gb for the rd_ <dig> dataset.
finally, the gsnap and pass mappers were highly memory-consuming; for the rd_ <dig> dataset, gsnap needed  <dig> gb of ram, with a peak usage of  <dig> gb while pass needed about  <dig> gb of ram with a peak usage of  <dig> gb.
the ram requirement increased proportionally with the dataset size for half of the mappers tested, while for bowtie <dig>  bwasw, mosaik, novoalign, smalt, segemehl, shrimp <dig>  and tmap memory consumption was about the same for all dataset sizes.
these experiments revealed that the computational resource requirements varied considerably among the mappers, from a few megabytes to  <dig> gb.
the datasets all contain only single-reads with different mean sizes.
the time required for the sequencing process relies mainly on the biotechnological part of the protocol  but the runtime of the mapping step could also constitute a bottleneck for some mappers.
figure  <dig> shows runtime measurements for each mapper running with the three real datasets.
the mappers had very different runtimes that were not all proportional to the dataset size.
snap was very quick and needed only about  <dig> minutes to map the rd_ <dig> dataset.
however, this runtime was for the program run in ‘any-best’ mode, which is always quicker than the other modes.
srmapper, ssaha <dig>  pass, bowtie <dig>  tmap, smalt, bwasw, shrimp <dig>  and mosaik needed less than  <dig> minutes to map the rd_ <dig> dataset and between  <dig> minute  and  <dig> minutes  to map the rd_ <dig> dataset.
bwa had quick runtimes of  <dig> and  <dig> minutes for the rd_ <dig> and rd_ <dig> datasets, respectively, but was slower with the biggest dataset rd_ <dig> , probably because bwa is optimized for short reads.
the slowest mappers were novoalign, gsnap, and segemehl.
the novoalign version used in this study could only be run with one thread which explains the long runtimes observed in this study .
gsnap runtimes were  <dig>   <dig> and  <dig> minutes, and segemehl needed  <dig>   <dig> and  <dig> minutes for the rd_ <dig>  rd_ <dig>  and rd_ <dig> datasets, respectively.
for all the mappers, the runtimes for the rd_ <dig> dataset  were longer.
generally speaking, the more bases in the dataset, the longer was the runtime, although the runtimes ranged from one minute to up to five hours.
mapper robustness
the accuracy of the sequencing technology is usually the criterion of first importance in the choice of a sequencer.
nevertheless, the mappers used to analyze the sequencing data must be able to efficiently take into account the inherent and inevitable raw data errors.
a robust mapper will permit compensation for sequencing defects and will contribute to maximizing coverage while limiting noise.
to evaluate mapper robustness, one method is to compute metrics  through a benchmark formed by simulated reads for which their original location in the genome and the number and type of introduced errors are known.
we used simulated datasets with varying error rates to compare mapper robustness.
to avoid simulation biases, we also studied mapper robustness with rabema  <cit>  using real sub-datasets.
these experiments were repeated for simulated datasets containing reads with a mean length of  <dig> and  <dig> bases .
overall, the f-measure values were marginally higher for the shortest reads but the mapper behaviors were similar to the behaviors observed with the dataset containing the reads with a mean length of  <dig> bases.
however, differences were observed for bwa and gsnap for which the f-measures were significantly better for the dataset with the shorter reads.
bwa was designed to map reads up to  <dig> bp long, which explained the better results with short reads.
the f-measures for the dataset of reads of  <dig> bases were lower for all the mappers and a significant decrease was observed for novoalign, bwa , and gsnap.
for the reads of  <dig> bases, novoalign showed an f-measure close to  <dig> even when the reads contained no errors.
this finding can be explained by the fact that novoalign truncates reads before alignment .
the maximum allowed read length is  <dig>  so all reads longer than this are truncated to  <dig> before mapping.
these experiments showed that most of the mappers were less robust when the indel rate increased, probably because most mappers are tuned mainly to deal with substitutions.
in the alignment step, the scoring parameters used by the mappers are often those currently used in bioinformatics, i.e.
from the evolutionary point of view, a substitution is less penalized than an insertion or a deletion.
however, in sequencing, mutations do not follow evolutionary rules; rather, they are dependent on the error model of the sequencing technology.
the ion torrent pgm, for example, is known to introduce more indels than substitutions into homopolymer stretches.
therefore, mapper robustness could probably be improved by modifying the scoring parameters in the alignment step by decreasing the indel penalty.
to test this idea, we changed the gap penalty for two mappers, shrimp <dig> and pass.
for shrimp <dig>  the gap open and extension penalties were set to match the penalty for substitutions .
for pass, a maximum gap of  <dig> bases was allowed with a gap open and extension penalty of  <dig>  the f-measures that were obtained with the adapted scoring parameters behaved in the same way as previously observed for these two methods, but they were globally better for all error rates than the f-measures obtained with the default parameters .
all the simulated datasets described above contained  <dig>  random reads , which could not be mapped onto the reference genome.
all the mappers, except smalt and tmap, returned all the random reads as unmapped.
for smalt and tmap, the longer the read length the higher the number of mapped random reads.
smalt mapped only a small number of the random reads , whereas tmap mapped around 10%, 12%, and 16% of the random reads in the  <dig>   <dig>  and  <dig> bases datasets, respectively, with around  <dig> matches.
these percentages are not negligible and indicated that the tmap strategy  was to map a maximum number of reads even if the mapping was not always relevant.
the reported alignments for the random reads were short and could be filtered out easily, but for non-expert users these reported hits will add to the complexity of the read mapping task.
in conclusion, most of the tested mappers were robust with low error rates.
segemehl showed the best f-measures even for datasets with high error rates and for all read lengths considered in this study.
mosaik, smalt, ssaha <dig>  bowtie <dig>  and shrimp <dig> correctly mapped a major part of the read datasets.
the results also showed that to handle ion torrent reads, mappers need to allow indels in the alignments, as was clear for all tested mappers except for srmapper and pass with their default settings.
we also demonstrated that decreasing the gap penalties could improve the mapping results for ion torrent data.
to avoid simulation biases, rabema was used to evaluate mapper performances with real datasets.
in rabema, a full-sensitivity algorithm was used to identify all possible matching intervals within a given error rate range for each read and the mapper evaluation was based on a metric called normalized found intervals , in which each interval for a read contributed 1/x points, where x is the number of alignments for the read.
the number of points was divided by the number of reads and multiplied by  <dig> to get the percentage.
figure  <dig> shows the percentage of nfi for mappers run in the ‘all’ mode with varying error rates.
only  <dig> mappers were considered because bwasw, snap, and srmapper cannot be run in ‘all’ mode.
all the mappers identified between 100% and 95% of the nfi for datasets with no errors.
however, for datasets with errors, the nfi fell rapidly to below 10% of nfi for some mappers , while others  maintained a high nfi percentage for datasets with up to a 4% error rate and finished at between  <dig> and 20% nfi for an 8% error rate .
only segemehl, shrimp <dig>  and bowtie <dig> maintained nfi above 80%, even at an 8% error rate.
the experiments were repeated with datasets that contained reads  <dig> and  <dig> bases long .
the ranking and behavior of the  <dig> mappers were similar to those obtained with datasets containing read lengths of  <dig>  except novoalign which was significantly better with the shorter reads.
for datasets with reads  <dig> bases long, the behavior of most of the mappers was similar to the behavior observed with 200-base long reads but the nfi percentages were lower.
the novoalign plot with several increases and decreases was atypical and only around 16% of the 400-base reads were mapped, probably because novoalign trims reads to a maximum length of  <dig> bases.
bwa identified around 100% nfi in the 400-base reads dataset with no errors, while with an error rate of 8% the nfi only fell to 40%.
this behavior for bwa was surprising when compared with its behavior in the previous experiments; however, it can be explained by the definition of nfi used by rabema.
in rabema, reads do not have to be aligned over their entire length to be considered as correctly mapped; so, many of the short alignments returned by bwa were classified as correct by rabema, which was not the case with our new definition.
the analysis of mapper performances on real datasets with rabema indicated that bowtie <dig>  segemehl, and shrimp <dig> were better than the other mappers, even for datasets with high error rates and regardless of the read lengths.
similar observations and similar rankings were obtained with the real and simulated datasets.
this double strategy built our confidence in the conclusions drawn from these experiments and confirmed that our simulator generated reads that were similar to sequencer generated reads .
study of repeats
the study and analysis of repeated sequences is as important for small microbial genomes, especially for bacterial genomes, as it is for eukaryotic genomes.
repeats in bacterial genomes represent a smaller proportion of the total genomic dna that they do in eukaryotic genomes, but the repeated elements are usually longer.
mapper behavior when dealing with repetitive regions in a reference genome is, therefore, an important parameter when the dna repeat regions may also be informative regions.
to study the ability of a mapper to report all possible positions for a read in a repeated sequence, we used an artificial genome containing five repeats.
in theory, a mapper, in ‘all’ mode, must report  <dig> hits for each repeat-located read.
figure  <dig> shows the percentage of repeat-located reads correctly reported by the mappers with reads of  <dig> bases, subdivided in classes depending on the number of hits found.
for each of the repeat-located reads, the number of locations in a repeat were counted.
note that bwasw and snap can report only one hit  and srmapper is limited to all-best hits.
most of the mappers were able to map repeat-located reads in at least one repeat , except for bwa and pass.
only two mappers  retrieved a large proportion  of the  <dig> hits and four few others  retrieved an average proportion of the  <dig> hits .
the other mappers performed quite poorly in this task, retrieving only a small percentage or none of the  <dig> hits.
with 100-base and 400-base reads, the mappers gave better and worse global results, respectively, than they did with the 200-base reads .
in conclusion, smalt was very good at retrieving multi-mapped reads whatever the read length, while gsnap, mosaik, and shrimp <dig> also gave correct results.
tmap was better with longer reads and novoalign was better with shorter reads.
mappers that cannot be run in ‘all-mode’ or that are not able to deal with indels  are not suitable for identifying multi-mapped reads.
mutation discovery
distinguishing between sequencing or mapping errors and true genetic variations is a challenge in variant analysis.
exome sequencing and genome re-sequencing require robust mapping results with as little noise as possible to identify a mutation of interest and to limit false positive mutations.
real reads from e. coli dh10b sequencing were mapped onto a genome sequence in which mutations with known positions and types  had been introduced artificially.
freebayes software  <cit>  was used to call variants, and precision and recall values were computed for mutation discovery in a reference genome with varying mutation rates.
figure  <dig> shows the precision and recall values obtained for mutation discovery with real datasets containing reads of  <dig> bases and a theoretical depth of 40x.
generally, precision and especially recall decreased when the mutation rate was increased in the reference genome.
in all the experiments, the precision values were high, indicating that the mutations predicted by the variant caller from the mapping files were mainly correct for all mappers.
most of the tested mappers presented good precision and recall values for all mutation rates; the exceptions were bwa, novoalign, pass and srmapper.
srmapper and pass presented lower precision and recall values than all the other mappers mostly because these two mappers do not allow for indels in the alignments, which decreased the precision of the mapping  and made the variant calling less accurate.
the mutation discovery performances of bwa and novoalign diminished when the mutation rate reached 5%.
it should be noted that for these two mappers, the percentage of mapped reads, and therefore the mean depth, was low compared with the percentage of mapped reads for the other mappers .
this reduced number of mapped reads did not permit the accurate detection of mutations in the reference genome.
roc curves were constructed , which confirmed the mutation discovery results that we obtained.
the experiments were repeated with simulated datasets .
the conclusions that were drawn were similar to those obtained with the real datasets; however, the precision and recall values were lower for all mappers.
we also performed similar experiments with real and simulated datasets for read lengths of  <dig> and  <dig> bases .
mapper behavior was similar regardless of the read length, except for bwa and novoalign.
these two mappers showed better values with reads of  <dig> bases, and showed near zero recall values with reads of  <dig> bases.
these results were not surprising because bwa was designed for short reads and novoalign truncates reads to a maximum length of  <dig> bases.
the behavior of the mappers in variant discovery was coherent with the results obtained in the robustness study and could be deduced from them.
for example, srmapper and bwa show a significant decrease in f-measure values when the error rate increased and similar behavior has been observed when the mutation rate was increased in the reference genomes.
variant discovery is impacted directly by the quality of the mapper alignments, i.e.
position and type of edit operations .
the definition of a correctly mapped read introduced in this study is more stringent than for previous studies, because it takes into account the correctness of the alignment .
these results demonstrated that the method we used to evaluate mapper robustness was efficient.
for the simulated data, similar behavior was observed for all the mappers and for all datasets but with lower precision and recall values than was observed for the real data.
this decrease could be explained by a lower error rate in the real data than in the simulated data.
we performed complementary analyses to observe the precision and recall values obtained with lower sequencing error rates .
when reads were generated without errors, the precision and recall values were close to  <dig>  precision and recall values were closer to the values obtained for the real dataset values when reads were generated with  <dig> % deletions,  <dig> % insertions, and  <dig> % substitutions, suggesting that the real dataset used here contained less than 2% sequencing errors.
these experiments again showed that the data simulated with curesim have characteristics that are similar to the real data produced by the ion torrent pgm.
finally, because we used simulated data, the impact of sequencing depth in mutation discovery could be tested.
we used shrimp <dig> because this mapper behaved well in the variant discovery experiments.
the same procedure was applied with four different read datasets of  <dig> bases with mean depths of 20x, 80x, 160x, and 320x .
the precision and recall values were lower with a mean sequencing depth of 20x and were equivalent for the other tested sequencing depths.
these results showed that a mean sequencing depth of 40x was enough to call variations correctly.
increasing the depth of sequencing did not seem to improve the quality of variant calling.
these experiments showed that most of the tested mappers gave correct results in mutation discovery even when used with their default settings.
the only exceptions were the bwa, novoalign, pass, and srmapper mappers.
srmapper and pass do not allow indels in alignments.
these kinds of mappers should be avoided for variant calling analysis.
discussion
here, a benchmark procedure to compare mappers for hts that can be applied to any sequencing platforms and any applications is described.
the different steps involved in this procedure are shown in figure  <dig>  in step  <dig>  a list of mappers is defined.
depending on the sequencing technology and the application, the most appropriate mapper can be selected for use.
in step  <dig>  real datasets are collected and simulated datasets are generated before being mapped onto the reference genome.
step  <dig> is a comparison step based on four criteria: mapper computational resource and time requirements; mapper robustness; mapper behavior with repetitive regions; and mapper mutation discovery ability.
the benchmark procedure uses simulated and real datasets to provide the user with a robust method for mapper comparison.
the results obtained can be used to answer questions such as: how much ram is required?
how long will it take to map a set of reads?
how does the robustness vary in relation to the error rate?
how does a mapper deal with multi-mapped reads?
could a mapper be used with a distant reference genome?
what is the quality of the reported alignment?
answers to these questions can help users chose a mapper that best fits a particular application and sequencing platform.
this procedure could also be used to evaluate performances of a newly developed mapper or to optimize parameters of already existing mappers.
we also presented a new read simulator, curesim , which generates synthetic hts reads for the major letter-base sequencing platforms.
users can fix the mutation rates, the read lengths, and can generate random reads.
several error distribution modes are available and particular attention was paid to special cases in which several introduced errors in the same read can lower the number of errors because of compensatory changes.
curesimeval is a complementary tool that evaluates the mapping quality from sam files produced by aligning curesim simulated reads with any mapper.
curesim and curesimeval are freely available at http://www.pegase-biosciences.com/tools/curesim/.
the curesim suite has been developed in java and is distributed as jar files to be operating system independent and easy to use by non-expert users.
we used the curesim suite in a mapper comparison with ion torrent data applied to small genomes.
to obtain a robust evaluation procedure, we introduced a new definition for mapping correctness.
this newly introduced definition is more stringent than the previous ones because the end of the alignment and the number of mutations were considered in addition to the start position.
the mapper robustness results obtained with the curesim suite simulated data matched the results obtained with real datasets and rabema, demonstrating that the curesim suite simulated reads with characteristics similar to real reads.
we performed completely independent experiments to evaluate the mutation discovery ability of the mappers and found that the results obtained for mapper robustness can also be used to predict the mutation discovery ability of the mappers.
variant calling efficiency is directly dependent on the alignment quality obtained by the mapping algorithms.
checking whether a mapped read is in its expected position is not sufficient because the position and number of edit operations in the produced alignment must also be as close as possible to the expected alignment.
the sequencing errors in ion torrent reads are mainly indels.
for mappers that are unable to deal correctly with indels, the resulting alignments, even those at the expected positions, can result in biased mapping that could impact the variant calling results.
all our results demonstrated the reliability of our evaluation method.
our benchmark procedure was applied to ion torrent data from small genome sequencing.
mapping algorithms and previous mapper comparison studies focused mainly on short reads with substitutions  and on large reference genomes .
these evaluation studies therefore are poorly informative for mapping new technology sequencer data with different error models.
additionally, some features of bacterial and small genomes, such as possible extreme gc content, make the extrapolation from the previous studies difficult.
for example, very high or low gc content percentages create a higher probability of encountering homopolymers in the genome sequence, which can significantly increase the number of indels in homopolymers.
our benchmark procedure for ion torrent data with bacterial genomes did not reveal a single best mapper but rather indicated several options depending on the particular application and technology.
when only a desktop computer with  <dig> gb of ram is available, users can select a mapper that is not highly memory-consuming; for example, novoalign, smalt, srmapper, bowtie <dig>  mosaik, or segemehl.
novoalign, segemehl, and gsnap require very long runtimes, while snap is very fast .
other mappers had runtimes shorter than  <dig> minutes for the bigger datasets.
concerning mapper robustness with varying error rates, all mappers manage to correctly map reads when the sequencing error rate was low; however, some mappers were clearly not suitable for use with datasets containing high error rates .
segemehl presented good f-measure values with all the tested datasets even at high error rates.
mosaik, shrimp <dig>  and bowtie <dig> also gave correct results.
smalt was well fitted to retrieve all hits for repeat-located reads and gsnap, mosaik, and shrimp <dig> also give correct results in this task.
one of these mappers is therefore suitable for the identification of unique and non-unique reads, whereas pass, bwa, bwasw, snap, and srmapper are not.
mapper behavior for mutation discovery with datasets with varying mutation rates using a close but not identical reference genome is of special interest because often only the genome of a closely-related species is available as a reference.
mutation discovery ability is also important when the genomes of two closely-related strains are compared to detect variants or mutated strains, for example.
in such cases, mappers need to produce accurate alignments so that true mutations can be detected.
all the mappers tested here showed good precision and recall with all tested mutation rates and for all datasets, except bwa, novoalign, pass, and srmapper.
our results show that some mappers dealt correctly with the ion torrent data although they were not initially designed for this technology.
for example, shrimp <dig> which was designed for illumina, solid, and  <dig> reads, showed robust results with ion torrent data.
the mapper default parameters were used deliberately in this study to mimic the general case of a non-expert user; therefore, different results could have been obtained with other parameter settings.
even with the default settings, several mappers that can be used with ion torrent data were identified.
additionally, we showed that the mapping results could be improved by adapting the parameter settings to the error model, for example, by decreasing the indel penalty with shrimp <dig>  for ion torrent data, our study demonstrated that to be efficient a mapper had to allow indels in the alignments and that the results were more reliable when the mapping algorithm allowed multi-mapped reads.
the mutation discovery experiments showed that a sequencing depth of 40x was enough to correctly call variants.
CONCLUSIONS
all the different applications that arise from hts technologies need not have the same mapping characteristics.
some applications may require robust mapping that deals with high error rates while others may require the ability to deal with repeats, for example, when re-sequencing is performed for bacterial variant identification aimed at efficiently detecting mutations and indels.
mappers such as ssaha <dig>  tmap, shrimp <dig>  or bowtie <dig> will support the detection of mutations even at high rates and without the necessity for deep sequencing.
in other applications, such as amplicon sequencing to study of repeated motifs , the ability to map correctly on repeat regions will be essential and a mapper like smalt, which performs such tasks very well even though its robustness is not among the highest could be used.
however, for some specific applications, such as the discovery of mutations in viral genomes, mappers such as bowtie <dig>  segemehl, and shrimp <dig> with strong robustness could be used because accurate mapping of the maximum number of reads, especially the few that bear the mutation, is essential  <cit> .
for some applications, it could be better to use a combination of mappers; for example, in pathogen identification, the strain might be unknown.
in this case, snap can be used to quickly identify a close reference genome among a set of available genomes, then a more robust mapper can be used to identify mutations or unique reads.
the correct choice of mapper is crucial in hts data analysis.
in this paper, we have presented a benchmark procedure to compare mapping algorithms that are used currently in hts.
therefore, we introduced a stringent definition of mapping correctness together with a new read simulator, curesim, to generate simulated reads with controlled type, rate, and/or distribution of errors along the reads.
the read simulator is freely distributed along with a tool to evaluate the mapping quality, curesimeval; both are available at http://www.pegase-biosciences.com/tools/curesim/.
this procedure was applied to small genomes with ion torrent data.
our results do not lead to the selection of a unique, omnipotent mapper but rather show that the choice of mapper has to be application and sequencing technology driven.
our study also demonstrates that a combination of several complementary mappers could significantly improve the mapping step in pipelines.
possible combinations should be tested and evaluated using the same approach.
the benchmark procedure presented here greatly helps in the choice of a good mapper for a given application and dataset.
this procedure could also be used to evaluate a newly developed mapper or to optimize parameters of an already existing one.
an optimized solution for read mapping, adapted to sequencing technology and biological applications, will help compensate for hts defects.

 high-throughput sequencing technology has recently shown a rapid and impressive development and this has led to the production of gigabases of sequence in a few hours for only a fraction of the former cost. hts has produced an explosion of knowledge in genetics and genomics thanks to the development of specific applications such as genome re-sequencing. this technological evolution was paralleled by the development of new algorithms to deal with the quantity and the quality of reads produced. a fundamental analysis steps in re-sequencing approaches is the mapping of the reads onto a reference genome. this step, which involves the accurate positioning of reads onto a reference genome sequence, is highly important because it determines the global quality of downstream analyses. the algorithms used for this step are called mappers. mappers have to be sensitive and accurate and, if possible, fast and not too computationally demanding. they should be able to find the true position of each read on a reference genome and ideally distinguish between technical sequencing errors and natural genetic variations. in recent years many mappers have been developed and distributed. two studies have classified mappers using a wide variety of features that include: the type of data, their application, the sequencing platform, the read length, the allowed error rate, parallel implementation, the ability to deal with multi-mapped reads, the input and output formats, and the available parameters. mappers have multiplied and so has the range of possible settings. hence, the growing difficulty in selecting a mapper has been raised in recent studies aimed at evaluating mapper performances through a multiplicity of comparison criteria. some of these studies have focused on mapper sensitivity. schbath et al. studied the ability of mappers to identify unique versus multi-mapped reads using a well-controlled benchmark containing reads with exactly three mismatches. hatem et al. introduced a benchmarking suite to analyze mapping tools, which consists of tests that cover input properties and algorithmic features. in addition to the difficulty in determining evaluation criteria, choosing an appropriate evaluation method, i.e. how to compare mappers according to the evaluation criteria, and using the appropriate metrics, are also problematical. using real datasets to evaluate mapper performances allows only a rough assessment and classification of mappers by comparing the percentage of mapped reads, but does not reveal the actual accuracy of mappers. attempts have been made to avoid this pitfall using simulated datasets in which the original read positions are known. another difficulty lies in the accurate definition of what a correctly mapped read is. the basic definition is to consider a read as correctly mapped if the original location is retrieved. ruffalo et al. broadened this definition by adding a condition on the quality score, which had to be superior to a given threshold. in a more recent paper, a new definition was introduced in which a read was considered to be correctly mapped if the mapping criteria were not violated, i.e. contained less errors than the threshold parameter set by the user. using simulated data allows numerical values to be obtained and compared between a set of mappers. however, simulated data do not have the same characteristics as real data, even when an error model based on real data is used. real hts data present biases that can be very difficult to simulate. additionally, the current definition of the mapping correctness based only on the original start location presents some weaknesses: a read can have several correct positions on the reference sequence and sequencing errors or true genetic variations can lead to a better alignment in a genome position different from the original one. holtgrewe et al. introduced the interval definition, rather than the genome position, to describe a read mapping and used a full-sensitivity algorithm to identify all possible matching intervals within a given error rate range for each read. this method has been implemented in rabema, a tool that evaluates the result of arbitrary read mappers that support the sam output format with real and simulated datasets. our analysis of the published literature on mapper evaluation led us to conclude that for a complete and robust comparison of mappers, real and simulated datasets should be used. using real datasets avoids simulation biases and gives a real picture of mapper behavior, whereas simulated datasets are benchmarks from which all parameters can be controlled. additionally, a sound, more complete definition of what constitutes a correctly mapped read needs to be considered. in all the previous studies, mapper performance was evaluated using large eukaryotic genomes and, for the most part, short illumina or illumina-like reads data were used, except in where datasets were evaluated with a reduced number of mappers and metrics. the type of sequencing errors and their rate is inherent to the sequencing technology and more precisely to the nucleotide elongation detection methods used. for example, life technologies sequencing by oligonucleotide ligation and detection technology showed a strong bias in its coverage of repetitive elements, whereas the illumina reversible dye-terminator sequencing technology mainly caused substitutions. pyrosequencing on solid support and ion semiconductor sequencing technology produced indel errors associated with homopolymer-regions. in the published evaluations, the criteria that were tested and the default parameters of the mappers were usually chosen to address or deal with substitution-type errors and are, therefore, less informative for mapping the reads from new technologies like the ion torrent platform. furthermore, the analysis of small microbial genomes compared with the analysis of large eukaryotic genomes poses other challenges because microbial genomes contain a wide range of gc content, which is sometimes extreme. very high or very low gc content means that there is a high probability of encountering homopolymers in a genome sequence and this is known to be a specific problem for pyrosequencing and ion semiconductor sequencers. a recent development in the hts technologies has made available benchtop sequencers targeted at the quick and inexpensive sequencing of small to moderate-sized genomes, mainly bacteria, viruses, fungi, and parasites. small microbial genome sequences could be considered to present a simpler, less demanding mapping process compared with the mapping process for larger eukaryotic genomes. however, this is only partially true because the characteristics of small microbial genomes are not the same as those of eukaryotic genomes. the questions of interest are also usually different and, consequently, the expected mapping quality criteria are not exactly the same. whole genome sequencing or re-sequencing is an important application in the new field of microorganism characterization using hts. for instance, clinical diagnosis and the epidemiological study of microbial strain circulation will be profoundly remodeled in the near future by the use of hts, which should, very soon, be used as a characterization approach for pathogens and which will probably slowly replace the present pcr and biochemical based characterization methods. in this particular context the re-sequencing applications and derived analyses are in the front-line of research and development. the focus includes the sequencing of the entire length of a microbial genome and the analysis of obtained reads by mapping them onto one or several reference strains to identify potential relevant changes in the studied genome. the aim is to accurately identify the gain or loss in genetic elements as well as small changes to predict a potential new phenotype or a derived new pathogenicity profile. this requirement poses several challenges, the most important of which is the necessity to distinguish true genetic variations from sequencing errors. in this paper, we focus on the evaluation of mappers in the context of whole genome sequencing or re-sequencing for small microbial, mainly bacterial genomes. we tested mappers, mostly using their default settings to be in the general context of non-expert users. we selected four criteria to match this context: computational resource and time requirements, robustness of mapping through the evaluation of precision, recall and f-measure, ability to report positions for reads in repetitive regions, and ability to retrieve true genetic variation positions. to evaluate a mapper’s robustness on simulated datasets, we introduced a new definition of a correctly mapped read. in addition to the original start position that was used in most previous studies, the end position as well as the numbers of insertions, deletions, and substitutions in the alignment were also used to classify the mapping of a read as correct. this definition is more stringent than the previous ones because it implies that it is a full-length read alignment and that the error count is correct. indeed, sequencing errors can mean that the original location of a read is not necessarily the best alignment location. using mappers tuned to report all possible hits and to accept a higher error rate than the error rate introduced in simulated reads, it should be possible to retrieve the original location in addition to potential equivalent or better hits. with the new definition of a correctly mapped read used in this study, we ensured that the mapper was able to retrieve the expected original alignment despite inevitable sequencing errors in the reads, thereby allowing a true evaluation of the mapper’s robustness. the analysis was applied to data generated by the ion torrent personal genome machine, a newly arrived technology dedicated mainly to small genome sequencing, for which mapper performances have not yet been evaluated. reads from real datasets and artificially simulated reads were used. simulated reads were generated using a new customizable read simulator, curesim, which can generate reads of user-determined lengths with insertions, deletions, and substitutions introduced at a controlled rate and with an adjustable error distribution along the read. curesim and curesimeval, a script that can be used to evaluate mapping quality, were developed in java to run on all operating systems and are freely available at http://www.pegase-biosciences.com/tools/curesim/. we have shown that in microbial genome sequencing, some mappers, such as segemehl, present higher robustness than others, especially when the number of sequencing errors was high. other mappers are more robust for other applications that demand other quality criteria. for example, bwasw, shrimp smalt, ssaha and tmap, might perform particularly well for sequencing focused on rare variant discovery because they show a robust discrimination of variations. smalt can localize most of the positions of reads located in repeated regions. some mappers, such as novoalign, smalt and srmapper, needed very small memory resources, while snap was very fast and required only about two minutes to process the bigger datasets used in this study. these results emphasize the observation that mapper choice is application dependent and users should carefully consider the targeted aim before choosing a mapper. the evaluation approach presented here, together with the developed tools can be considered as a general method to evaluate existing or in-development mappers and could prove interesting in the evaluation of the performances of mappers for the coming third generation of sequencers that may have yet another type and rate of errors.  computational resource requirement and time measurement all mapping processes involve the alignment of millions of reads onto a reference sequence. this is true even for small genome sequencing projects where the small size of the reference sequence is generally compensated by the multiplicity of samples to be analyzed. in clinical microbiology, the time and the computational resources required for the analysis are critical; therefore, these factors also need to be evaluated for the different mappers. all the mappers tested were run with threads and the memory consumption and runtime were recorded for three different ion torrent datasets rd_ rd_ and rd_ these three datasets contain real single-reads with different mean sizes and are described in table the reference genome used was escherichia coli str. k- substr. dh10b for the rd_ and rd_ datasets and escherichia coli str. k- substr. mg for the rd_ dataset. figure shows the memory consumption for each mapper for the real datasets when the indexing and mapping steps were considered together. novoalign, smalt, and srmapper needed very low memory resources. it should be noted that srmapper was developed to run on a computer with gb of ram for genomes the size of the human genome, but, in such a case, it can be run only in ‘all-best’ mode and does not allow indels in the mapping. the novoalign version used in this study was the free academic version that has not been implemented in parallel. a second group comprising bowtie mosaik, and segemehl, needed less than gb of ram, while a third group, bwa, bwasw, and tmap, needed less than gb of memory. bwa had peak memory usage of gb for rd_ and of more than gb for the rd_ dataset. bwa was developed to map short reads of up to bases, which may explain the high peak usage for 400-base reads. shrimp snap, and ssaha required more ram and ssaha needed about gb for the rd_ dataset. finally, the gsnap and pass mappers were highly memory-consuming; for the rd_ dataset, gsnap needed gb of ram, with a peak usage of gb while pass needed about gb of ram with a peak usage of gb. the ram requirement increased proportionally with the dataset size for half of the mappers tested, while for bowtie bwasw, mosaik, novoalign, smalt, segemehl, shrimp and tmap memory consumption was about the same for all dataset sizes. these experiments revealed that the computational resource requirements varied considerably among the mappers, from a few megabytes to gb. the datasets all contain only single-reads with different mean sizes. the time required for the sequencing process relies mainly on the biotechnological part of the protocol but the runtime of the mapping step could also constitute a bottleneck for some mappers. figure shows runtime measurements for each mapper running with the three real datasets. the mappers had very different runtimes that were not all proportional to the dataset size. snap was very quick and needed only about minutes to map the rd_ dataset. however, this runtime was for the program run in ‘any-best’ mode, which is always quicker than the other modes. srmapper, ssaha pass, bowtie tmap, smalt, bwasw, shrimp and mosaik needed less than minutes to map the rd_ dataset and between minute and minutes to map the rd_ dataset. bwa had quick runtimes of and minutes for the rd_ and rd_ datasets, respectively, but was slower with the biggest dataset rd_, probably because bwa is optimized for short reads. the slowest mappers were novoalign, gsnap, and segemehl. the novoalign version used in this study could only be run with one thread which explains the long runtimes observed in this study. gsnap runtimes were and minutes, and segemehl needed and minutes for the rd_ rd_ and rd_ datasets, respectively. for all the mappers, the runtimes for the rd_ dataset were longer. generally speaking, the more bases in the dataset, the longer was the runtime, although the runtimes ranged from one minute to up to five hours. mapper robustness the accuracy of the sequencing technology is usually the criterion of first importance in the choice of a sequencer. nevertheless, the mappers used to analyze the sequencing data must be able to efficiently take into account the inherent and inevitable raw data errors. a robust mapper will permit compensation for sequencing defects and will contribute to maximizing coverage while limiting noise. to evaluate mapper robustness, one method is to compute metrics through a benchmark formed by simulated reads for which their original location in the genome and the number and type of introduced errors are known. we used simulated datasets with varying error rates to compare mapper robustness. to avoid simulation biases, we also studied mapper robustness with rabema using real sub-datasets. these experiments were repeated for simulated datasets containing reads with a mean length of and bases. overall, the f-measure values were marginally higher for the shortest reads but the mapper behaviors were similar to the behaviors observed with the dataset containing the reads with a mean length of bases. however, differences were observed for bwa and gsnap for which the f-measures were significantly better for the dataset with the shorter reads. bwa was designed to map reads up to bp long, which explained the better results with short reads. the f-measures for the dataset of reads of bases were lower for all the mappers and a significant decrease was observed for novoalign, bwa, and gsnap. for the reads of bases, novoalign showed an f-measure close to even when the reads contained no errors. this finding can be explained by the fact that novoalign truncates reads before alignment. the maximum allowed read length is so all reads longer than this are truncated to before mapping. these experiments showed that most of the mappers were less robust when the indel rate increased, probably because most mappers are tuned mainly to deal with substitutions. in the alignment step, the scoring parameters used by the mappers are often those currently used in bioinformatics, i.e. from the evolutionary point of view, a substitution is less penalized than an insertion or a deletion. however, in sequencing, mutations do not follow evolutionary rules; rather, they are dependent on the error model of the sequencing technology. the ion torrent pgm, for example, is known to introduce more indels than substitutions into homopolymer stretches. therefore, mapper robustness could probably be improved by modifying the scoring parameters in the alignment step by decreasing the indel penalty. to test this idea, we changed the gap penalty for two mappers, shrimp and pass. for shrimp the gap open and extension penalties were set to match the penalty for substitutions. for pass, a maximum gap of bases was allowed with a gap open and extension penalty of the f-measures that were obtained with the adapted scoring parameters behaved in the same way as previously observed for these two methods, but they were globally better for all error rates than the f-measures obtained with the default parameters. all the simulated datasets described above contained random reads, which could not be mapped onto the reference genome. all the mappers, except smalt and tmap, returned all the random reads as unmapped. for smalt and tmap, the longer the read length the higher the number of mapped random reads. smalt mapped only a small number of the random reads, whereas tmap mapped around 10%, 12%, and 16% of the random reads in the and bases datasets, respectively, with around matches. these percentages are not negligible and indicated that the tmap strategy was to map a maximum number of reads even if the mapping was not always relevant. the reported alignments for the random reads were short and could be filtered out easily, but for non-expert users these reported hits will add to the complexity of the read mapping task. in conclusion, most of the tested mappers were robust with low error rates. segemehl showed the best f-measures even for datasets with high error rates and for all read lengths considered in this study. mosaik, smalt, ssaha bowtie and shrimp correctly mapped a major part of the read datasets. the results also showed that to handle ion torrent reads, mappers need to allow indels in the alignments, as was clear for all tested mappers except for srmapper and pass with their default settings. we also demonstrated that decreasing the gap penalties could improve the mapping results for ion torrent data. to avoid simulation biases, rabema was used to evaluate mapper performances with real datasets. in rabema, a full-sensitivity algorithm was used to identify all possible matching intervals within a given error rate range for each read and the mapper evaluation was based on a metric called normalized found intervals, in which each interval for a read contributed 1/x points, where x is the number of alignments for the read. the number of points was divided by the number of reads and multiplied by to get the percentage. figure shows the percentage of nfi for mappers run in the ‘all’ mode with varying error rates. only mappers were considered because bwasw, snap, and srmapper cannot be run in ‘all’ mode. all the mappers identified between 100% and 95% of the nfi for datasets with no errors. however, for datasets with errors, the nfi fell rapidly to below 10% of nfi for some mappers, while others maintained a high nfi percentage for datasets with up to a 4% error rate and finished at between and 20% nfi for an 8% error rate. only segemehl, shrimp and bowtie maintained nfi above 80%, even at an 8% error rate. the experiments were repeated with datasets that contained reads and bases long. the ranking and behavior of the mappers were similar to those obtained with datasets containing read lengths of except novoalign which was significantly better with the shorter reads. for datasets with reads bases long, the behavior of most of the mappers was similar to the behavior observed with 200-base long reads but the nfi percentages were lower. the novoalign plot with several increases and decreases was atypical and only around 16% of the 400-base reads were mapped, probably because novoalign trims reads to a maximum length of bases. bwa identified around 100% nfi in the 400-base reads dataset with no errors, while with an error rate of 8% the nfi only fell to 40%. this behavior for bwa was surprising when compared with its behavior in the previous experiments; however, it can be explained by the definition of nfi used by rabema. in rabema, reads do not have to be aligned over their entire length to be considered as correctly mapped; so, many of the short alignments returned by bwa were classified as correct by rabema, which was not the case with our new definition. the analysis of mapper performances on real datasets with rabema indicated that bowtie segemehl, and shrimp were better than the other mappers, even for datasets with high error rates and regardless of the read lengths. similar observations and similar rankings were obtained with the real and simulated datasets. this double strategy built our confidence in the conclusions drawn from these experiments and confirmed that our simulator generated reads that were similar to sequencer generated reads. study of repeats the study and analysis of repeated sequences is as important for small microbial genomes, especially for bacterial genomes, as it is for eukaryotic genomes. repeats in bacterial genomes represent a smaller proportion of the total genomic dna that they do in eukaryotic genomes, but the repeated elements are usually longer. mapper behavior when dealing with repetitive regions in a reference genome is, therefore, an important parameter when the dna repeat regions may also be informative regions. to study the ability of a mapper to report all possible positions for a read in a repeated sequence, we used an artificial genome containing five repeats. in theory, a mapper, in ‘all’ mode, must report hits for each repeat-located read. figure shows the percentage of repeat-located reads correctly reported by the mappers with reads of bases, subdivided in classes depending on the number of hits found. for each of the repeat-located reads, the number of locations in a repeat were counted. note that bwasw and snap can report only one hit and srmapper is limited to all-best hits. most of the mappers were able to map repeat-located reads in at least one repeat, except for bwa and pass. only two mappers retrieved a large proportion of the hits and four few others retrieved an average proportion of the hits. the other mappers performed quite poorly in this task, retrieving only a small percentage or none of the hits. with 100-base and 400-base reads, the mappers gave better and worse global results, respectively, than they did with the 200-base reads. in conclusion, smalt was very good at retrieving multi-mapped reads whatever the read length, while gsnap, mosaik, and shrimp also gave correct results. tmap was better with longer reads and novoalign was better with shorter reads. mappers that cannot be run in ‘all-mode’ or that are not able to deal with indels are not suitable for identifying multi-mapped reads. mutation discovery distinguishing between sequencing or mapping errors and true genetic variations is a challenge in variant analysis. exome sequencing and genome re-sequencing require robust mapping results with as little noise as possible to identify a mutation of interest and to limit false positive mutations. real reads from e. coli dh10b sequencing were mapped onto a genome sequence in which mutations with known positions and types had been introduced artificially. freebayes software was used to call variants, and precision and recall values were computed for mutation discovery in a reference genome with varying mutation rates. figure shows the precision and recall values obtained for mutation discovery with real datasets containing reads of bases and a theoretical depth of 40x. generally, precision and especially recall decreased when the mutation rate was increased in the reference genome. in all the experiments, the precision values were high, indicating that the mutations predicted by the variant caller from the mapping files were mainly correct for all mappers. most of the tested mappers presented good precision and recall values for all mutation rates; the exceptions were bwa, novoalign, pass and srmapper. srmapper and pass presented lower precision and recall values than all the other mappers mostly because these two mappers do not allow for indels in the alignments, which decreased the precision of the mapping and made the variant calling less accurate. the mutation discovery performances of bwa and novoalign diminished when the mutation rate reached 5%. it should be noted that for these two mappers, the percentage of mapped reads, and therefore the mean depth, was low compared with the percentage of mapped reads for the other mappers. this reduced number of mapped reads did not permit the accurate detection of mutations in the reference genome. roc curves were constructed, which confirmed the mutation discovery results that we obtained. the experiments were repeated with simulated datasets. the conclusions that were drawn were similar to those obtained with the real datasets; however, the precision and recall values were lower for all mappers. we also performed similar experiments with real and simulated datasets for read lengths of and bases. mapper behavior was similar regardless of the read length, except for bwa and novoalign. these two mappers showed better values with reads of bases, and showed near zero recall values with reads of bases. these results were not surprising because bwa was designed for short reads and novoalign truncates reads to a maximum length of bases. the behavior of the mappers in variant discovery was coherent with the results obtained in the robustness study and could be deduced from them. for example, srmapper and bwa show a significant decrease in f-measure values when the error rate increased and similar behavior has been observed when the mutation rate was increased in the reference genomes. variant discovery is impacted directly by the quality of the mapper alignments, i.e. position and type of edit operations. the definition of a correctly mapped read introduced in this study is more stringent than for previous studies, because it takes into account the correctness of the alignment. these results demonstrated that the method we used to evaluate mapper robustness was efficient. for the simulated data, similar behavior was observed for all the mappers and for all datasets but with lower precision and recall values than was observed for the real data. this decrease could be explained by a lower error rate in the real data than in the simulated data. we performed complementary analyses to observe the precision and recall values obtained with lower sequencing error rates. when reads were generated without errors, the precision and recall values were close to precision and recall values were closer to the values obtained for the real dataset values when reads were generated with % deletions, % insertions, and % substitutions, suggesting that the real dataset used here contained less than 2% sequencing errors. these experiments again showed that the data simulated with curesim have characteristics that are similar to the real data produced by the ion torrent pgm. finally, because we used simulated data, the impact of sequencing depth in mutation discovery could be tested. we used shrimp because this mapper behaved well in the variant discovery experiments. the same procedure was applied with four different read datasets of bases with mean depths of 20x, 80x, 160x, and 320x. the precision and recall values were lower with a mean sequencing depth of 20x and were equivalent for the other tested sequencing depths. these results showed that a mean sequencing depth of 40x was enough to call variations correctly. increasing the depth of sequencing did not seem to improve the quality of variant calling. these experiments showed that most of the tested mappers gave correct results in mutation discovery even when used with their default settings. the only exceptions were the bwa, novoalign, pass, and srmapper mappers. srmapper and pass do not allow indels in alignments. these kinds of mappers should be avoided for variant calling analysis. discussion here, a benchmark procedure to compare mappers for hts that can be applied to any sequencing platforms and any applications is described. the different steps involved in this procedure are shown in figure in step a list of mappers is defined. depending on the sequencing technology and the application, the most appropriate mapper can be selected for use. in step real datasets are collected and simulated datasets are generated before being mapped onto the reference genome. step is a comparison step based on four criteria: mapper computational resource and time requirements; mapper robustness; mapper behavior with repetitive regions; and mapper mutation discovery ability. the benchmark procedure uses simulated and real datasets to provide the user with a robust method for mapper comparison. the results obtained can be used to answer questions such as: how much ram is required? how long will it take to map a set of reads? how does the robustness vary in relation to the error rate? how does a mapper deal with multi-mapped reads? could a mapper be used with a distant reference genome? what is the quality of the reported alignment? answers to these questions can help users chose a mapper that best fits a particular application and sequencing platform. this procedure could also be used to evaluate performances of a newly developed mapper or to optimize parameters of already existing mappers. we also presented a new read simulator, curesim, which generates synthetic hts reads for the major letter-base sequencing platforms. users can fix the mutation rates, the read lengths, and can generate random reads. several error distribution modes are available and particular attention was paid to special cases in which several introduced errors in the same read can lower the number of errors because of compensatory changes. curesimeval is a complementary tool that evaluates the mapping quality from sam files produced by aligning curesim simulated reads with any mapper. curesim and curesimeval are freely available at http://www.pegase-biosciences.com/tools/curesim/. the curesim suite has been developed in java and is distributed as jar files to be operating system independent and easy to use by non-expert users. we used the curesim suite in a mapper comparison with ion torrent data applied to small genomes. to obtain a robust evaluation procedure, we introduced a new definition for mapping correctness. this newly introduced definition is more stringent than the previous ones because the end of the alignment and the number of mutations were considered in addition to the start position. the mapper robustness results obtained with the curesim suite simulated data matched the results obtained with real datasets and rabema, demonstrating that the curesim suite simulated reads with characteristics similar to real reads. we performed completely independent experiments to evaluate the mutation discovery ability of the mappers and found that the results obtained for mapper robustness can also be used to predict the mutation discovery ability of the mappers. variant calling efficiency is directly dependent on the alignment quality obtained by the mapping algorithms. checking whether a mapped read is in its expected position is not sufficient because the position and number of edit operations in the produced alignment must also be as close as possible to the expected alignment. the sequencing errors in ion torrent reads are mainly indels. for mappers that are unable to deal correctly with indels, the resulting alignments, even those at the expected positions, can result in biased mapping that could impact the variant calling results. all our results demonstrated the reliability of our evaluation method. our benchmark procedure was applied to ion torrent data from small genome sequencing. mapping algorithms and previous mapper comparison studies focused mainly on short reads with substitutions and on large reference genomes. these evaluation studies therefore are poorly informative for mapping new technology sequencer data with different error models. additionally, some features of bacterial and small genomes, such as possible extreme gc content, make the extrapolation from the previous studies difficult. for example, very high or low gc content percentages create a higher probability of encountering homopolymers in the genome sequence, which can significantly increase the number of indels in homopolymers. our benchmark procedure for ion torrent data with bacterial genomes did not reveal a single best mapper but rather indicated several options depending on the particular application and technology. when only a desktop computer with gb of ram is available, users can select a mapper that is not highly memory-consuming; for example, novoalign, smalt, srmapper, bowtie mosaik, or segemehl. novoalign, segemehl, and gsnap require very long runtimes, while snap is very fast. other mappers had runtimes shorter than minutes for the bigger datasets. concerning mapper robustness with varying error rates, all mappers manage to correctly map reads when the sequencing error rate was low; however, some mappers were clearly not suitable for use with datasets containing high error rates. segemehl presented good f-measure values with all the tested datasets even at high error rates. mosaik, shrimp and bowtie also gave correct results. smalt was well fitted to retrieve all hits for repeat-located reads and gsnap, mosaik, and shrimp also give correct results in this task. one of these mappers is therefore suitable for the identification of unique and non-unique reads, whereas pass, bwa, bwasw, snap, and srmapper are not. mapper behavior for mutation discovery with datasets with varying mutation rates using a close but not identical reference genome is of special interest because often only the genome of a closely-related species is available as a reference. mutation discovery ability is also important when the genomes of two closely-related strains are compared to detect variants or mutated strains, for example. in such cases, mappers need to produce accurate alignments so that true mutations can be detected. all the mappers tested here showed good precision and recall with all tested mutation rates and for all datasets, except bwa, novoalign, pass, and srmapper. our results show that some mappers dealt correctly with the ion torrent data although they were not initially designed for this technology. for example, shrimp which was designed for illumina, solid, and reads, showed robust results with ion torrent data. the mapper default parameters were used deliberately in this study to mimic the general case of a non-expert user; therefore, different results could have been obtained with other parameter settings. even with the default settings, several mappers that can be used with ion torrent data were identified. additionally, we showed that the mapping results could be improved by adapting the parameter settings to the error model, for example, by decreasing the indel penalty with shrimp for ion torrent data, our study demonstrated that to be efficient a mapper had to allow indels in the alignments and that the results were more reliable when the mapping algorithm allowed multi-mapped reads. the mutation discovery experiments showed that a sequencing depth of 40x was enough to correctly call variants.  all the different applications that arise from hts technologies need not have the same mapping characteristics. some applications may require robust mapping that deals with high error rates while others may require the ability to deal with repeats, for example, when re-sequencing is performed for bacterial variant identification aimed at efficiently detecting mutations and indels. mappers such as ssaha tmap, shrimp or bowtie will support the detection of mutations even at high rates and without the necessity for deep sequencing. in other applications, such as amplicon sequencing to study of repeated motifs, the ability to map correctly on repeat regions will be essential and a mapper like smalt, which performs such tasks very well even though its robustness is not among the highest could be used. however, for some specific applications, such as the discovery of mutations in viral genomes, mappers such as bowtie segemehl, and shrimp with strong robustness could be used because accurate mapping of the maximum number of reads, especially the few that bear the mutation, is essential. for some applications, it could be better to use a combination of mappers; for example, in pathogen identification, the strain might be unknown. in this case, snap can be used to quickly identify a close reference genome among a set of available genomes, then a more robust mapper can be used to identify mutations or unique reads. the correct choice of mapper is crucial in hts data analysis. in this paper, we have presented a benchmark procedure to compare mapping algorithms that are used currently in hts. therefore, we introduced a stringent definition of mapping correctness together with a new read simulator, curesim, to generate simulated reads with controlled type, rate, and/or distribution of errors along the reads. the read simulator is freely distributed along with a tool to evaluate the mapping quality, curesimeval; both are available at http://www.pegase-biosciences.com/tools/curesim/. this procedure was applied to small genomes with ion torrent data. our results do not lead to the selection of a unique, omnipotent mapper but rather show that the choice of mapper has to be application and sequencing technology driven. our study also demonstrates that a combination of several complementary mappers could significantly improve the mapping step in pipelines. possible combinations should be tested and evaluated using the same approach. the benchmark procedure presented here greatly helps in the choice of a good mapper for a given application and dataset. this procedure could also be used to evaluate a newly developed mapper or to optimize parameters of an already existing one. an optimized solution for read mapping, adapted to sequencing technology and biological applications, will help compensate for hts defects. 
to measure robustness, we introduced a new definition for a correctly mapped read taking into account not only the expected start position of the read but also the end position and the number of indels and substitutions. in this paper, we present a benchmark procedure to compare mapping algorithms used in hts using both real and simulated datasets and considering four evaluation criteria: computational resource and time requirements, robustness of mapping, ability to report positions for reads in repetitive regions, and ability to retrieve true genetic variation positions. a benchmark procedure to compare hts data mappers is introduced with a new definition for the mapping correctness as well as tools to generate simulated reads and evaluate mapping quality. this benchmark procedure can be used to evaluate existing or in-development mappers as well as to optimize parameters of a chosen mapper for any application and any sequencing platform. we applied our benchmark procedure to evaluate mappers in the context of whole genome sequencing of small genomes with ion torrent data for which such a comparison has not yet been established. the application of this procedure to ion torrent data from the whole genome sequencing of small genomes has allowed us to validate our benchmark procedure and demonstrate that it is helpful for selecting a mapper based on the intended application, questions to be addressed, and the technology used. a fundamental step in hts data analysis is the mapping of reads onto reference sequences.
 the rapid evolution in high-throughput sequencing technologies has opened up new perspectives in several research fields and led to the production of large volumes of sequence data. a fundamental step in hts data analysis is the mapping of reads onto reference sequences. choosing a suitable mapper for a given technology and a given application is a subtle task because of the difficulty of evaluating mapping algorithms.  in this paper, we present a benchmark procedure to compare mapping algorithms used in hts using both real and simulated datasets and considering four evaluation criteria: computational resource and time requirements, robustness of mapping, ability to report positions for reads in repetitive regions, and ability to retrieve true genetic variation positions. to measure robustness, we introduced a new definition for a correctly mapped read taking into account not only the expected start position of the read but also the end position and the number of indels and substitutions. we developed curesim, a new read simulator, that is able to generate customized benchmark data for any kind of hts technology by adjusting parameters to the error types. curesim and curesimeval, a tool to evaluate the mapping quality of the curesim simulated reads, are freely available. we applied our benchmark procedure to evaluate mappers in the context of whole genome sequencing of small genomes with ion torrent data for which such a comparison has not yet been established.  a benchmark procedure to compare hts data mappers is introduced with a new definition for the mapping correctness as well as tools to generate simulated reads and evaluate mapping quality. the application of this procedure to ion torrent data from the whole genome sequencing of small genomes has allowed us to validate our benchmark procedure and demonstrate that it is helpful for selecting a mapper based on the intended application, questions to be addressed, and the technology used. this benchmark procedure can be used to evaluate existing or in-development mappers as well as to optimize parameters of a chosen mapper for any application and any sequencing platform. high-throughput sequencingmapping algorithmsmapper comparisonread simulator 