BACKGROUND
sequencing by capillary electrophoresis, known as the sanger method  <cit> , has been employed in many historically significant large-scale sequencing projects and is regarded as the gold standard in terms of both read length and sequencing accuracy  <cit> .
several massively parallel dna sequencing  technologies have recently emerged, including the roche/ <dig> gs flx system, the illumina/solexa genome analyser, and the ab solid system, which are able to generate a few orders of magnitude more bases per instrument run, being considerably less expensive than the sanger method  <cit> .
these technologies are enabling researchers and practitioners to efficiently sequence genomes, leading to very significant advances in biology and medicine.
however, the huge volume of data produced by mpds technologies creates important computational challenges  <cit> .
moreover, the different platform-specific data characteristics require different algorithmic approaches.
for instance, some applications may use the  <dig> titanium platform to produce reads  <dig> bases long, some other studies may employ a solid system set to produce short reads of  <dig> bases, and yet other projects may use the illumina system to produce  <dig> ×  <dig> bases paired-end reads.
given their large variety, it would be rather difficult for a single algorithm to handle all kinds of data optimally.
when sequencing a new organism, one is usually faced with the problem of assembling the sequence fragments  together from scratch.
however, when a sufficiently close sequence is already known, one may choose to use it as a reference and proceed by first mapping the reads to this reference and then determining the new sequence by extracting the consensus from the mapping results.
the former strategy is called de novo sequencing, while the latter is known as re-sequencing.
several tools have recently been developed for generating assemblies from short reads, e.g  <cit> .
similarly, several methods have been proposed to address the problem of efficiently mapping mpds reads to a reference sequence, like  <cit> , to cite a few.
as referred before, the sheer volume of data generated by mpds technologies , and the need to align reads to large reference genomes limit the applicability of standard techniques.
indeed, in a typical application, we may have to align hundreds of millions of reads to a reference genome that can be as large as few gigabases, a job that cannot be efficiently achieved through standard dynamic programming procedures.
one way to speed up the read alignment task is to resort to approximate indexing techniques.
a first generation of aligners was based on hash tables of k-mers.
some of them, like ssaha <dig>  <cit> , build tables of k-mers of the target sequence, whilst others, like newbler  <cit> , index the reads, thus presumably requiring re-indexing for each new run.
recent developments in the field of compressed approximate indexes have led to a new family of alignment algorithms such as segemehl  <cit> , which uses an enhanced suffix array , and bwa-sw  <cit> , which uses a fm-index  to accelerate smith-waterman alignments.
yet the number of aligners that support gs flx pyrosequencing data is, as of today, relatively scarce compared to other technologies, most notably illumina.
moreover, some of these tools find their origins in the days before the advent of the new sequencing technologies and only later were adapted to cope with new kinds of data  <cit> , and some others target multiple kinds of data  <cit>  being not optimized for pyrosequencing data.
given this state of affairs, we argue that there is still room for improvement in the realm of publicly available aligners specifically designed for high-throughput pyrosequencing data.
in this paper we present a new method for the alignment of pyrosequencing reads, like those produced by the  <dig> gs flx platform.
by focusing on this specific technology, our procedure manages to explore its data characteristics to achieve improved performance over other mainstream methods.
like many of those methods, ours also builds an index of the target  sequence to accelerate the alignment.
it then employs a multiple seed heuristic to anchor the best candidate alignments.
contrary to other seed-based alignment tools, our strategy adds more flexibility by dispensing with the need of determining the number and length of the seeds beforehand.
our heuristic relies on some assumptions that can be reasonably expected to hold true for re-sequencing projects based on pyrosequencing data, namely, that the optimal alignments are mostly composed of relatively large chunks of exact matches interspersed by small, possibly gapped, divergent regions.
a banded dynamic programming is used to finish up the candidate multiple seed alignments considering user-specified error constraints.
a detailed description of the algorithm and data structures is given in the "implementation" section.
in the "results and discussion" section, we present a comparison between our method and a set of tools of widespread use for the local alignment of pyrosequencing reads.
we base our discussion on results obtained with both real and simulated data.
implementation
compressed indexes
the main data structure for sequence pattern matching is an index.
indexes reduce the time for matching a pattern because they restrict the search to the positions where it may occur instead of scanning the whole text.
one of the most, if not the most, popular index structures is the suffix tree , which is obtained by identifying common prefixes of the different suffixes of the represented text to nodes of a tree  for an example).
in such a structure, a pattern can be searched following edges with matching labels down from the root.
each leaf of the suffix tree represents a suffix of the text and, more generally, each node represents the subset of suffixes corresponding to the leaves of the subtree rooted at that node.
the downside of indexes is that they need to be constructed a priori and have a bad reputation of using too much space.
despite providing for fast searching algorithms, suffix trees are particularly known for this bad characteristic.
a popular alternative to suffix trees are suffix arrays , that require asymptotically the same space, o computer words for a text of size n, but with a smaller proportionality factor.
suffix arrays are obtained by ordering the suffixes of a text lexicographically.
a correspondence can be established between nodes of the suffix tree and contiguous intervals of the suffix array.
an example of a suffix array is shown in figure  <dig>  detailed descriptions of string matching and indexes, including the ones mentioned here, are widely available  <cit> .
recent research on indexes has focused on the fact that pointer representations require o bits whereas the original text  requires only n log σ bits, where σ is the alphabet size, e.g.
<dig> for dna and  <dig> for proteins.
in an effort to reduce this gap, new indexes have been designed which became collectively known as compressed indexes  <cit>  due to the fact that they rely heavily on data compression techniques.
in spite of their reduced space, compressed indexes can be made to allow for an even broader range of operations than classical indexes, like generalized branching, that combines blocks of letters instead of just one letter at a time  <cit> .
our method uses an implementation of the fmindex  <cit>  optimized for the dna alphabet.
the fmindex is a compressed index based on the burrows-wheeler transform   <cit>  requiring only o bits of memory space.
the bwt of a text t is obtained by appending an extra symbol $ to t, and then sorting all cyclic permutations  of t$ according to the lexicographic order, with $ being the lowest symbol.
thus a bwt of a text is essentially equivalent to its suffix array.
in fact, they are related through the formula bwt = t - 1], for i =  <dig>  .
.
.
, | t |.
an example of a bwt is given in figure  <dig>
the seed-based search approach
due to the relatively large size of the gs flx reads, it is not practical to use a plain index-based exact matching algorithm.
some sort of backtracking strategy could be used to allow for errors but, since the number of possible comparisons increases exponentially with their number, this becomes rapidly inefficient.
instead, we propose a seed-based search heuristic that explores the characteristics of the pyrosequencing data and of the bona fide alignments that are likely to arise in the context of re-sequencing applications.
since the error rates are usually low  <cit> , and the prevalent type of pyrosequencing errors are small indels, with mismatches being much less common, we conjecture that the optimal alignments are expected to be formed of large chunks of exact matches interspersed by divergent gapped regions.
moreover, since the read lengths are of a few hundred bases, we can expect the exact match regions to be large enough so we can use segments of them, called seeds, as a backbone to pin down the position of the alignment on the reference sequence, or at least reduce the amount of candidate positions to a manageable number of possibilities that can be tested individually.
in this case, the optimal alignments can be obtained by expanding the candidate multiple seed matches into alignments of the whole read by filling up the remaining regions and selecting those with best overall scores.
our strategy for choosing the seeds consists in approximately partitioning the read into maximal exact match blocks in a greedy fashion.
more precisely, let r = r <dig> ⋯ rm be the read.
the procedure starts at the first position of the read and uses the index to find the largest prefix of the read with exact occurrences in the reference sequence, say r = r <dig> ⋯ rl.
in practice we obtain the equivalent of an interval of the bwt which contains the positions of the reference sequence g at which r occurs.
obviously, by maximality, r does not occur in g. this happens because none of the occurrences of r is followed by rl+ <dig> or, put another way, because there is a mismatch between rl+ <dig> and the letter following each occurrence of r .
if rl+ <dig> ≠ rl, then we set r as the first seed and proceed as above to find the next seed starting from position l +  <dig>  if, however, we have rl+ <dig> = rl, this means that the difference occurred in the middle of a homopolymer , most likely due to an insertion sequencing error.
in this case, we set r as the first seed as before, but advance the cursor to the start of the next homopolymer, i.e.
to the smallest l' >l s.t.
rl' ≠ rl' -  <dig>  we repeat this process until the end of the read is reached.
once we have the set of seeds and their individual positions in the reference sequence, we need to identify subsets of occurrences of distinct seeds that are in accordance with their original order and spacing in the read, which can then serve as a support for the final alignments.
more formally, let g be the reference sequence, and r be a read.
let also s <dig>  .
.
.
, sk be k substrings of r such that r = s1a1s2a <dig> ⋯ ak-1skak, where, for i =  <dig>  .
.
.
, k, si denote the seeds and ai denote the substrings in between them.
for each i =  <dig>  .
.
.
, k, let oi ≥  <dig> be the number of exact occurrences of si in g. we then have  ways to choose a set containing one occurrence for each of the k substrings.
let p =  be one of such o tuples of distinct seed occurrences.
if, for i =  <dig>  .
.
.
, k - <dig>  we have pi ≤ pi+ <dig> and |ai| ∈i pi+ <dig> -  |ai| + ∈i for some given ∈i ≥  <dig>  then we say that these occurrences are coherent.
hence, a coherent set of seed occurrences is composed of positions which respect the relative order of the corresponding seeds in the read and such that, for any two consecutive seeds, the distance between their occurrences lies within a restricted interval around the actual distance between those seeds in the read.
for the sake of flexibility, we do not restrict ourselves to coherent sets containing occurrences of all the seeds.
we also take partial sets containing occurrences of only some of those seeds as good candidates for further expansion.
the set of seeds  and their occurrence positions in g can be obtained with the index in linear time.
however, the number of combinations of occurrences of different seeds can be rather large, especially if some of them are short.
this makes it impractical to test all possibilities for coherence.
nonetheless, most of these combinations will typically be non-coherent, and if we care to previously sort the set of occurrences of each seed, we can efficiently search for coherent combinations using, again, a greedy strategy, simply by scanning the seed matches in g from left to right, partitioning them into maximal non-overlapping sets of consecutive coherent matches.
although this might seem a rough approach at first glance, in fact this strategy has shown to be adequate because of the relatively large size of the seeds and small separation between them, which makes it difficult for occurrences of two consecutive seeds to be interspersed with an occurrence of a third one.
once the potential read occurrences indicated by coherent multiple-seed matches are found, the algorithm runs a banded needleman-wunsch dynamic programming procedure with gotoh's modifications to align the non-seed segments of the read to their counterparts in the genome.
that is, if we have a coherent set of occurrences of the seeds sj <dig>  .
.
.
, sjq , s.t.
the read decomposes into r = b0sj1b1sj <dig> b <dig> bq-1sjq bq, and the reference sequence into g = c0sj1c1sj <dig> c <dig> ⋯ cq 1sjqcq, then we align each pair , for i =  <dig>  .
.
.
, q. of course, for both ends,  and , we perform semi-global alignments.
the largest candidate coherent multiple seed matches are extended this way and accepted as a read occurrence if either the overall alignment identity stays above a given threshold percentage t or, alternatively, if the sum of the errors in-between the seeds and at the extremes of the read do not exceed a pre-established number e. the algorithm can be chosen to report all the accepted occurrences or only the one with the least errors.
the strategy described above is illustrated in figure  <dig>
synthetic data generation
in order to evaluate the algorithms in a controlled setting, we generated artificial data sets with a procedure inspired by empirical studies on gs flx data  <cit> , and designed to yield reads with characteristics similar to real data.
in our procedure, n random contiguous subsequences are extracted from a given 'source' sequence g. the lengths of these initial subsequences are drawn from a normal distribution with mean μl and standard deviation σl, also provided as input.
next, these subsequences are modified to simulate sequencing errors as follows.
in the gs flx high-throughput pyrosequencing procedure, the template molecules are sequenced one maximal homopolymer at a time , as opposed to one base at a time in the traditional sanger method.
hence, the most common type of error in pyrosequencing consist in the misinterpretation of the intensity of the signal that determines the length of the homopolymer being read, leading to an insertion or deletion of identical consecutive bases in the read, relative to the actual template sequence.
miscalled base errors  also occur but they are comparatively much less frequent.
sequence quality is known to be non-uniform along the read, being lower at the extremes, particularly towards the 3' end.
also, errors tend to affect long homopolymers more than short ones.
however, for the sake of simplification, we consider that errors are uniformly distributed along the read and that the prevalence and size of indels are not affected by the length of the homopolymers.
more precisely, the procedure takes in three parameters psub, pins and pdel which correspond to the probabilities of having a substitution, an insertion or a deletion in any given homopolymer, regardless of its length and position in the read.
moreover, these events are considered to be mutually exclusive, that is, we assume that for any homopolymer being sequenced, there can either be a substitution error with probability psub, an insertion with probability pins, a deletion with probability pdel or it can be correctly sequenced with probability  <dig> - .
whenever a mismatch takes place, the miscalled base is randomly chosen according to substitution probabilities indicated in a matrix m, given as input.
each row/column of m corresponds to a nucleotide and the element m indicates the probability for a to be miscalled as  b. as for indels, the lengths of the gaps are drawn from zipfian distributions, which are discrete power-law distributions with mass function .
in our case, ω is a positive integer parameter that corresponds to a maximum allowed gap size, and γ >  <dig> controls the shape of the distribution: the greater its value the higher the prevalence of small gaps.
we use specific exponent parameters, γins and γdel, for insertion and deletion operations, respectively.
RESULTS
we evaluated tapyr against other mainstream mapping tools which are also able to deal with high-throughput pyrosequencing reads, namely bwa-sw  <cit> , ssaha <dig>  <cit> , segemehl  <cit> , gassst  <cit> , and newbler  <cit> .
our analyses were performed with real and simulated data sets, with the objective of assessing the efficiency and accuracy of the aforementioned tools in the context of re-sequencing projects.
results on real data
the biological data sets we used, summarized in table  <dig>  encompass a reasonable variety of organism types, including two bacteria , one protozoan , one nematode , one insect , and one human chromosome.
they also cover re-sequencing applications with reads from individuals of the same species , different and mutated strains of the same species , and different species .
s. pneumoniae
○ srr <dig> s. pneumoniae sp195
○ srr <dig> e. coli k-12
○ srr <dig> e. coli etec ws3080a
○ srr <dig> p. falciparum 3d7
○ srr <dig> p. falciparum 3d7
○ srr <dig> p. falciparum 3d7
c. elegans
wormdb rel.
○ srr <dig> d. pseudoobscura bogotana er
○ srr <dig> d. pseudoobscura bogotana er
○ srr <dig> human individual na15510
○ srr <dig> human individual na15510
○ srr <dig> human individual na15510
○ srr <dig> human individual na15510
biological data sets used for the evaluation of the algorithms.
the read data sets were downloaded from the sequence read archive  public repository.
in this experiment, we wanted to analyze the ability of the algorithms to produce high coverage mappings, which directly relates to the proportion of reads that can be successfully mapped.
high coverage is essential to the successful completion of a re-sequencing project, with about 20-25× coverage being required for optimal results with the gs flx technology  <cit> .
attaining such high levels depends naturally on the amount of available data, but equally on the capacity of the alignment tool to map the reads correctly, especially in the presence of inevitable sequencing errors and natural variation.
the other aspect we wanted to assess was the efficiency of the algorithms in terms of computation time.
efficiency is a critical aspect for any algorithm in modern high-throughput data processing pipelines, given the rapid increase in the volume of data being produced.
the results of our tests are shown in table  <dig>  in that comparison, we included two lines corresponding to tapyr being set to report alignments with at least 50%  and 85%  identity.
these illustrative values match the default options of other tools: 85% for segemehl, and 50% for ssaha <dig>  as can be seen, in almost all direct comparison scenarios, tapyr has shown to be several times faster than the other tools.
as for the number of successfully aligned reads , we notice first that the other algorithms display quite similar figures, with no tool consistently aligning more reads that the others.
with the minimum identity threshold set to 85% , our method aligns a smaller quantity of reads.
however, if we investigate the number of errors  of the reported alignments by computing the ratio between the number of base-pair matches and the number errors , we see that tapyr is using a more conservative heuristic which tends to produce alignments of a higher identity level at the expense of dropping a slightly larger number of reads.
indeed, if we lower minimum identity requirement to 50% , then our tool aligns more reads than all the others in all data sets at comparable average error rates, and with a minimal time overhead.
results of the experiments performed with real biological data.
these tests were run on a linux server with  <dig> gb of ram.
tapyr was tested under two configurations, requiring alignments with at least 50% and 85% identity, designated by tapyr  <dig> and tapyr  <dig>  respectively.
the other tools were run with their default options for  <dig> data, except for the following modifications.
ssaha <dig> and segemehl were set to report only the best alignment for each read.
for gassst, we set the minimum identity to 85%  to match segemehl.
newbler was set not to generate large files  and to load the index into the main memory .
reported times refer to the total number of cpu-seconds that the process used directly, as given by the linux command time -f "%u".
the average base-pairs-per-error rates  are computed based on the best reported alignment of each read only.
results on synthetic data
we also performed tests using simulated data produced according to the procedure described in the methods section.
we generated three data sets of n =  <dig>  synthetic reads from the  <dig> mbp sequence of the human chromosome  <dig> .
these data sets are supposed to mimic the data obtained in a typical run of the gs flx instrument at different sequencing error levels.
hence, the first data set, hereafter referred to as hs <dig>  was generated with the read generator parameters set as μl =  <dig>  σl =  <dig>  pins = pdel = psub =  <dig> , ω =  <dig>  γins = γdel =  <dig>  and equiprobable substitution rates m = 1/ <dig>  for a ≠ b. in this setting, we have a 1% chance of each of the three kinds of error when reading a homopolymer.
for the second data set, hs <dig>  we increased the error levels of each kind to 5% by setting pins = pdel = psub =  <dig> , and for the third data set, we added a considerable amount of noise by setting pins = pdel = psub =  <dig>  .
the purpose of this experiment was mainly to test accuracy of the procedures by computing the fraction of reads mapped back to their original positions, as well as to assess the robustness of the heuristic to different error levels.
the results of the tests are shown in table  <dig>  we notice that all algorithms give quite accurate results in all the tested conditions.
in any case, tapyr behaved among the best in terms of accuracy, mapping virtually all reads correctly, showing thus resilience to noise up the tested levels.
moreover, as in the previous experiments, our method has confirmed to be fastest by a comfortable margin.
results of the experiments performed with synthetic data.
the data sets were generated as described in the text.
the tests were run under the same conditions as those with real data.
shown is the percentage of reads of each data set that were successfully aligned to their original positions in the reference sequence.
memory usage
we also measured the memory requirements of the evaluated tools in the tests with real data discussed above.
the figures presented in table  <dig> show the sizes of the index files on disk, when they exist, and the peak usage of main memory for the different data sets.
as it can be seen, bwa displayed the smallest requirements in absolute terms, followed closely by tapyr.
the other tools, especially those based on k-mer tables, demand substantially more memory.
as expected, tapyr's index files scale linearly with the size of the indexed genomes .
apart from the index, which is loaded into main memory, tapyr uses only a small additional amount of space , so that the total amount of required ram also scales linearly with the indexed genome .
these modest and predictable requirements make tapyr suitable for large genomes with moderately-sized machines.
memory requirements for the real biological data sets of table  <dig>  shown are the index file size  and peak main memory usage as measured by the tool tstime
http://bitbucket.org/gsauthof/tstime, except for newbler, whose overall memory demand was estimated through the linux htop tool.
CONCLUSIONS
the combination of state of the art indexing techniques and a seed-based search approach led to the development of a new read mapping method for high-throughput pyrosequencing data.
by using an effective heuristic which explores the characteristics of this particular kind of data in the context of typical re-sequencing applications, our method manages to achieve convincing performance in terms of speed and in terms of the number and precision of aligned reads, as demonstrated by our tests with real and simulated data.
in fact, our proposed solution has displayed class-leading cpu-time performance and excellent use of input reads in comparison to other mainstream tools.
an added-value of our procedure comes from the fact that it requires almost no external parameterization.
as a matter of fact, the main user options are end-of-the-chain cutoff parameters that concern the quality of the reported alignments in terms of minimal identity or maximal number of errors, having no consequence on the accuracy of the heuristic and only marginal impact on the overall execution time.
memory requirements are also on par with the best in this category of tools, being not only small in absolute terms but, more importantly, linearly proportional to the size of the input reference sequence by a small factor.
based on these results, we propose that tapyr constitutes an advantageous alternative for re-sequencing projects based on pyrosequencing data.
availability and requirements
project name: tapyr--tool for the alignment of pyrosequencing reads
project home page: http://www.tapyr.net
operating system: multiple
programming language: c
other requirements: none
license: gnu gpl
restrictions to use by non-academics: none additional
authors' contributions
af and ao conceived the project.
all authors have participated in the design and refinement of the method and in the analysis of the results.
ff wrote most of the code, helped by lr and pf.
ff and pf collected the data and performed the experiments.
the authors collectively drafted the manuscript.
pf, af, and ff wrote most of the text.
all authors revised and approved the final version.

 sequencing by capillary electrophoresis, known as the sanger method, has been employed in many historically significant large-scale sequencing projects and is regarded as the gold standard in terms of both read length and sequencing accuracy. several massively parallel dna sequencing technologies have recently emerged, including the roche/ gs flx system, the illumina/solexa genome analyser, and the ab solid system, which are able to generate a few orders of magnitude more bases per instrument run, being considerably less expensive than the sanger method. these technologies are enabling researchers and practitioners to efficiently sequence genomes, leading to very significant advances in biology and medicine. however, the huge volume of data produced by mpds technologies creates important computational challenges. moreover, the different platform-specific data characteristics require different algorithmic approaches. for instance, some applications may use the titanium platform to produce reads bases long, some other studies may employ a solid system set to produce short reads of bases, and yet other projects may use the illumina system to produce × bases paired-end reads. given their large variety, it would be rather difficult for a single algorithm to handle all kinds of data optimally. when sequencing a new organism, one is usually faced with the problem of assembling the sequence fragments together from scratch. however, when a sufficiently close sequence is already known, one may choose to use it as a reference and proceed by first mapping the reads to this reference and then determining the new sequence by extracting the consensus from the mapping results. the former strategy is called de novo sequencing, while the latter is known as re-sequencing. several tools have recently been developed for generating assemblies from short reads, e.g. similarly, several methods have been proposed to address the problem of efficiently mapping mpds reads to a reference sequence, like, to cite a few. as referred before, the sheer volume of data generated by mpds technologies, and the need to align reads to large reference genomes limit the applicability of standard techniques. indeed, in a typical application, we may have to align hundreds of millions of reads to a reference genome that can be as large as few gigabases, a job that cannot be efficiently achieved through standard dynamic programming procedures. one way to speed up the read alignment task is to resort to approximate indexing techniques. a first generation of aligners was based on hash tables of k-mers. some of them, like ssaha, build tables of k-mers of the target sequence, whilst others, like newbler, index the reads, thus presumably requiring re-indexing for each new run. recent developments in the field of compressed approximate indexes have led to a new family of alignment algorithms such as segemehl, which uses an enhanced suffix array, and bwa-sw, which uses a fm-index to accelerate smith-waterman alignments. yet the number of aligners that support gs flx pyrosequencing data is, as of today, relatively scarce compared to other technologies, most notably illumina. moreover, some of these tools find their origins in the days before the advent of the new sequencing technologies and only later were adapted to cope with new kinds of data, and some others target multiple kinds of data being not optimized for pyrosequencing data. given this state of affairs, we argue that there is still room for improvement in the realm of publicly available aligners specifically designed for high-throughput pyrosequencing data. in this paper we present a new method for the alignment of pyrosequencing reads, like those produced by the gs flx platform. by focusing on this specific technology, our procedure manages to explore its data characteristics to achieve improved performance over other mainstream methods. like many of those methods, ours also builds an index of the target sequence to accelerate the alignment. it then employs a multiple seed heuristic to anchor the best candidate alignments. contrary to other seed-based alignment tools, our strategy adds more flexibility by dispensing with the need of determining the number and length of the seeds beforehand. our heuristic relies on some assumptions that can be reasonably expected to hold true for re-sequencing projects based on pyrosequencing data, namely, that the optimal alignments are mostly composed of relatively large chunks of exact matches interspersed by small, possibly gapped, divergent regions. a banded dynamic programming is used to finish up the candidate multiple seed alignments considering user-specified error constraints. a detailed description of the algorithm and data structures is given in the "implementation" section. in the "results and discussion" section, we present a comparison between our method and a set of tools of widespread use for the local alignment of pyrosequencing reads. we base our discussion on results obtained with both real and simulated data. implementation compressed indexes the main data structure for sequence pattern matching is an index. indexes reduce the time for matching a pattern because they restrict the search to the positions where it may occur instead of scanning the whole text. one of the most, if not the most, popular index structures is the suffix tree, which is obtained by identifying common prefixes of the different suffixes of the represented text to nodes of a tree for an example). in such a structure, a pattern can be searched following edges with matching labels down from the root. each leaf of the suffix tree represents a suffix of the text and, more generally, each node represents the subset of suffixes corresponding to the leaves of the subtree rooted at that node. the downside of indexes is that they need to be constructed a priori and have a bad reputation of using too much space. despite providing for fast searching algorithms, suffix trees are particularly known for this bad characteristic. a popular alternative to suffix trees are suffix arrays, that require asymptotically the same space, o computer words for a text of size n, but with a smaller proportionality factor. suffix arrays are obtained by ordering the suffixes of a text lexicographically. a correspondence can be established between nodes of the suffix tree and contiguous intervals of the suffix array. an example of a suffix array is shown in figure detailed descriptions of string matching and indexes, including the ones mentioned here, are widely available. recent research on indexes has focused on the fact that pointer representations require o bits whereas the original text requires only n log σ bits, where σ is the alphabet size, e.g. for dna and for proteins. in an effort to reduce this gap, new indexes have been designed which became collectively known as compressed indexes due to the fact that they rely heavily on data compression techniques. in spite of their reduced space, compressed indexes can be made to allow for an even broader range of operations than classical indexes, like generalized branching, that combines blocks of letters instead of just one letter at a time. our method uses an implementation of the fmindex optimized for the dna alphabet. the fmindex is a compressed index based on the burrows-wheeler transform requiring only o bits of memory space. the bwt of a text t is obtained by appending an extra symbol $ to t, and then sorting all cyclic permutations of t$ according to the lexicographic order, with $ being the lowest symbol. thus a bwt of a text is essentially equivalent to its suffix array. in fact, they are related through the formula bwt = t - 1], for i =..., | t |. an example of a bwt is given in figure the seed-based search approach due to the relatively large size of the gs flx reads, it is not practical to use a plain index-based exact matching algorithm. some sort of backtracking strategy could be used to allow for errors but, since the number of possible comparisons increases exponentially with their number, this becomes rapidly inefficient. instead, we propose a seed-based search heuristic that explores the characteristics of the pyrosequencing data and of the bona fide alignments that are likely to arise in the context of re-sequencing applications. since the error rates are usually low, and the prevalent type of pyrosequencing errors are small indels, with mismatches being much less common, we conjecture that the optimal alignments are expected to be formed of large chunks of exact matches interspersed by divergent gapped regions. moreover, since the read lengths are of a few hundred bases, we can expect the exact match regions to be large enough so we can use segments of them, called seeds, as a backbone to pin down the position of the alignment on the reference sequence, or at least reduce the amount of candidate positions to a manageable number of possibilities that can be tested individually. in this case, the optimal alignments can be obtained by expanding the candidate multiple seed matches into alignments of the whole read by filling up the remaining regions and selecting those with best overall scores. our strategy for choosing the seeds consists in approximately partitioning the read into maximal exact match blocks in a greedy fashion. more precisely, let r = r ⋯ rm be the read. the procedure starts at the first position of the read and uses the index to find the largest prefix of the read with exact occurrences in the reference sequence, say r = r ⋯ rl. in practice we obtain the equivalent of an interval of the bwt which contains the positions of the reference sequence g at which r occurs. obviously, by maximality, r does not occur in g. this happens because none of the occurrences of r is followed by rl+ or, put another way, because there is a mismatch between rl+ and the letter following each occurrence of r. if rl+ ≠ rl, then we set r as the first seed and proceed as above to find the next seed starting from position l + if, however, we have rl+ = rl, this means that the difference occurred in the middle of a homopolymer, most likely due to an insertion sequencing error. in this case, we set r as the first seed as before, but advance the cursor to the start of the next homopolymer, i.e. to the smallest l' >l s.t. rl' ≠ rl' - we repeat this process until the end of the read is reached. once we have the set of seeds and their individual positions in the reference sequence, we need to identify subsets of occurrences of distinct seeds that are in accordance with their original order and spacing in the read, which can then serve as a support for the final alignments. more formally, let g be the reference sequence, and r be a read. let also s..., sk be k substrings of r such that r = s1a1s2a ⋯ ak-1skak, where, for i =..., k, si denote the seeds and ai denote the substrings in between them. for each i =..., k, let oi ≥ be the number of exact occurrences of si in g. we then have ways to choose a set containing one occurrence for each of the k substrings. let p = be one of such o tuples of distinct seed occurrences. if, for i =..., k - we have pi ≤ pi+ and |ai| ∈i pi+ - |ai| + ∈i for some given ∈i ≥ then we say that these occurrences are coherent. hence, a coherent set of seed occurrences is composed of positions which respect the relative order of the corresponding seeds in the read and such that, for any two consecutive seeds, the distance between their occurrences lies within a restricted interval around the actual distance between those seeds in the read. for the sake of flexibility, we do not restrict ourselves to coherent sets containing occurrences of all the seeds. we also take partial sets containing occurrences of only some of those seeds as good candidates for further expansion. the set of seeds and their occurrence positions in g can be obtained with the index in linear time. however, the number of combinations of occurrences of different seeds can be rather large, especially if some of them are short. this makes it impractical to test all possibilities for coherence. nonetheless, most of these combinations will typically be non-coherent, and if we care to previously sort the set of occurrences of each seed, we can efficiently search for coherent combinations using, again, a greedy strategy, simply by scanning the seed matches in g from left to right, partitioning them into maximal non-overlapping sets of consecutive coherent matches. although this might seem a rough approach at first glance, in fact this strategy has shown to be adequate because of the relatively large size of the seeds and small separation between them, which makes it difficult for occurrences of two consecutive seeds to be interspersed with an occurrence of a third one. once the potential read occurrences indicated by coherent multiple-seed matches are found, the algorithm runs a banded needleman-wunsch dynamic programming procedure with gotoh's modifications to align the non-seed segments of the read to their counterparts in the genome. that is, if we have a coherent set of occurrences of the seeds sj..., sjq, s.t. the read decomposes into r = b0sj1b1sj b bq-1sjq bq, and the reference sequence into g = c0sj1c1sj c ⋯ cq 1sjqcq, then we align each pair, for i =..., q. of course, for both ends, and, we perform semi-global alignments. the largest candidate coherent multiple seed matches are extended this way and accepted as a read occurrence if either the overall alignment identity stays above a given threshold percentage t or, alternatively, if the sum of the errors in-between the seeds and at the extremes of the read do not exceed a pre-established number e. the algorithm can be chosen to report all the accepted occurrences or only the one with the least errors. the strategy described above is illustrated in figure synthetic data generation in order to evaluate the algorithms in a controlled setting, we generated artificial data sets with a procedure inspired by empirical studies on gs flx data, and designed to yield reads with characteristics similar to real data. in our procedure, n random contiguous subsequences are extracted from a given'source' sequence g. the lengths of these initial subsequences are drawn from a normal distribution with mean μl and standard deviation σl, also provided as input. next, these subsequences are modified to simulate sequencing errors as follows. in the gs flx high-throughput pyrosequencing procedure, the template molecules are sequenced one maximal homopolymer at a time, as opposed to one base at a time in the traditional sanger method. hence, the most common type of error in pyrosequencing consist in the misinterpretation of the intensity of the signal that determines the length of the homopolymer being read, leading to an insertion or deletion of identical consecutive bases in the read, relative to the actual template sequence. miscalled base errors also occur but they are comparatively much less frequent. sequence quality is known to be non-uniform along the read, being lower at the extremes, particularly towards the 3' end. also, errors tend to affect long homopolymers more than short ones. however, for the sake of simplification, we consider that errors are uniformly distributed along the read and that the prevalence and size of indels are not affected by the length of the homopolymers. more precisely, the procedure takes in three parameters psub, pins and pdel which correspond to the probabilities of having a substitution, an insertion or a deletion in any given homopolymer, regardless of its length and position in the read. moreover, these events are considered to be mutually exclusive, that is, we assume that for any homopolymer being sequenced, there can either be a substitution error with probability psub, an insertion with probability pins, a deletion with probability pdel or it can be correctly sequenced with probability -. whenever a mismatch takes place, the miscalled base is randomly chosen according to substitution probabilities indicated in a matrix m, given as input. each row/column of m corresponds to a nucleotide and the element m indicates the probability for a to be miscalled as b. as for indels, the lengths of the gaps are drawn from zipfian distributions, which are discrete power-law distributions with mass function. in our case, ω is a positive integer parameter that corresponds to a maximum allowed gap size, and γ > controls the shape of the distribution: the greater its value the higher the prevalence of small gaps. we use specific exponent parameters, γins and γdel, for insertion and deletion operations, respectively.  we evaluated tapyr against other mainstream mapping tools which are also able to deal with high-throughput pyrosequencing reads, namely bwa-sw, ssaha, segemehl, gassst, and newbler. our analyses were performed with real and simulated data sets, with the objective of assessing the efficiency and accuracy of the aforementioned tools in the context of re-sequencing projects. results on real data the biological data sets we used, summarized in table encompass a reasonable variety of organism types, including two bacteria, one protozoan, one nematode, one insect, and one human chromosome. they also cover re-sequencing applications with reads from individuals of the same species, different and mutated strains of the same species, and different species. s. pneumoniae ○ srr s. pneumoniae sp195 ○ srr e. coli k-12 ○ srr e. coli etec ws3080a ○ srr p. falciparum 3d7 ○ srr p. falciparum 3d7 ○ srr p. falciparum 3d7 c. elegans wormdb rel. ○ srr d. pseudoobscura bogotana er ○ srr d. pseudoobscura bogotana er ○ srr human individual na15510 ○ srr human individual na15510 ○ srr human individual na15510 ○ srr human individual na15510 biological data sets used for the evaluation of the algorithms. the read data sets were downloaded from the sequence read archive public repository. in this experiment, we wanted to analyze the ability of the algorithms to produce high coverage mappings, which directly relates to the proportion of reads that can be successfully mapped. high coverage is essential to the successful completion of a re-sequencing project, with about 20-25× coverage being required for optimal results with the gs flx technology. attaining such high levels depends naturally on the amount of available data, but equally on the capacity of the alignment tool to map the reads correctly, especially in the presence of inevitable sequencing errors and natural variation. the other aspect we wanted to assess was the efficiency of the algorithms in terms of computation time. efficiency is a critical aspect for any algorithm in modern high-throughput data processing pipelines, given the rapid increase in the volume of data being produced. the results of our tests are shown in table in that comparison, we included two lines corresponding to tapyr being set to report alignments with at least 50% and 85% identity. these illustrative values match the default options of other tools: 85% for segemehl, and 50% for ssaha as can be seen, in almost all direct comparison scenarios, tapyr has shown to be several times faster than the other tools. as for the number of successfully aligned reads, we notice first that the other algorithms display quite similar figures, with no tool consistently aligning more reads that the others. with the minimum identity threshold set to 85%, our method aligns a smaller quantity of reads. however, if we investigate the number of errors of the reported alignments by computing the ratio between the number of base-pair matches and the number errors, we see that tapyr is using a more conservative heuristic which tends to produce alignments of a higher identity level at the expense of dropping a slightly larger number of reads. indeed, if we lower minimum identity requirement to 50%, then our tool aligns more reads than all the others in all data sets at comparable average error rates, and with a minimal time overhead. results of the experiments performed with real biological data. these tests were run on a linux server with gb of ram. tapyr was tested under two configurations, requiring alignments with at least 50% and 85% identity, designated by tapyr and tapyr respectively. the other tools were run with their default options for data, except for the following modifications. ssaha and segemehl were set to report only the best alignment for each read. for gassst, we set the minimum identity to 85% to match segemehl. newbler was set not to generate large files and to load the index into the main memory. reported times refer to the total number of cpu-seconds that the process used directly, as given by the linux command time -f "%u". the average base-pairs-per-error rates are computed based on the best reported alignment of each read only. results on synthetic data we also performed tests using simulated data produced according to the procedure described in the methods section. we generated three data sets of n = synthetic reads from the mbp sequence of the human chromosome. these data sets are supposed to mimic the data obtained in a typical run of the gs flx instrument at different sequencing error levels. hence, the first data set, hereafter referred to as hs was generated with the read generator parameters set as μl = σl = pins = pdel = psub =, ω = γins = γdel = and equiprobable substitution rates m = 1/ for a ≠ b. in this setting, we have a 1% chance of each of the three kinds of error when reading a homopolymer. for the second data set, hs we increased the error levels of each kind to 5% by setting pins = pdel = psub =, and for the third data set, we added a considerable amount of noise by setting pins = pdel = psub =. the purpose of this experiment was mainly to test accuracy of the procedures by computing the fraction of reads mapped back to their original positions, as well as to assess the robustness of the heuristic to different error levels. the results of the tests are shown in table we notice that all algorithms give quite accurate results in all the tested conditions. in any case, tapyr behaved among the best in terms of accuracy, mapping virtually all reads correctly, showing thus resilience to noise up the tested levels. moreover, as in the previous experiments, our method has confirmed to be fastest by a comfortable margin. results of the experiments performed with synthetic data. the data sets were generated as described in the text. the tests were run under the same conditions as those with real data. shown is the percentage of reads of each data set that were successfully aligned to their original positions in the reference sequence. memory usage we also measured the memory requirements of the evaluated tools in the tests with real data discussed above. the figures presented in table show the sizes of the index files on disk, when they exist, and the peak usage of main memory for the different data sets. as it can be seen, bwa displayed the smallest requirements in absolute terms, followed closely by tapyr. the other tools, especially those based on k-mer tables, demand substantially more memory. as expected, tapyr's index files scale linearly with the size of the indexed genomes. apart from the index, which is loaded into main memory, tapyr uses only a small additional amount of space, so that the total amount of required ram also scales linearly with the indexed genome. these modest and predictable requirements make tapyr suitable for large genomes with moderately-sized machines. memory requirements for the real biological data sets of table shown are the index file size and peak main memory usage as measured by the tool tstime http://bitbucket.org/gsauthof/tstime, except for newbler, whose overall memory demand was estimated through the linux htop tool.  the combination of state of the art indexing techniques and a seed-based search approach led to the development of a new read mapping method for high-throughput pyrosequencing data. by using an effective heuristic which explores the characteristics of this particular kind of data in the context of typical re-sequencing applications, our method manages to achieve convincing performance in terms of speed and in terms of the number and precision of aligned reads, as demonstrated by our tests with real and simulated data. in fact, our proposed solution has displayed class-leading cpu-time performance and excellent use of input reads in comparison to other mainstream tools. an added-value of our procedure comes from the fact that it requires almost no external parameterization. as a matter of fact, the main user options are end-of-the-chain cutoff parameters that concern the quality of the reported alignments in terms of minimal identity or maximal number of errors, having no consequence on the accuracy of the heuristic and only marginal impact on the overall execution time. memory requirements are also on par with the best in this category of tools, being not only small in absolute terms but, more importantly, linearly proportional to the size of the input reference sequence by a small factor. based on these results, we propose that tapyr constitutes an advantageous alternative for re-sequencing projects based on pyrosequencing data. availability and requirements project name: tapyr--tool for the alignment of pyrosequencing reads project home page: http://www.tapyr.net operating system: multiple programming language: c other requirements: none license: gnu gpl restrictions to use by non-academics: none additional authors' contributions af and ao conceived the project. all authors have participated in the design and refinement of the method and in the analysis of the results. ff wrote most of the code, helped by lr and pf. ff and pf collected the data and performed the experiments. the authors collectively drafted the manuscript. pf, af, and ff wrote most of the text. all authors revised and approved the final version. 
we present an efficient method for the local alignment of pyrosequencing reads produced by the gs flx system against a reference sequence. an evaluation performed using real and simulated data shows that our proposed method outperforms a number of mainstream tools on the quantity and quality of successful alignments, as well as on the execution time. over the past few years, new massively parallel dna sequencing technologies have emerged. among the most critical problems is that of efficiently and accurately mapping reads to a reference genome in the context of re-sequencing projects. however, these techniques also raise important computational difficulties mostly due to the huge volume of data produced, but also because of some of their specific characteristics such as read length and sequencing errors.
 over the past few years, new massively parallel dna sequencing technologies have emerged. these platforms generate massive amounts of data per run, greatly reducing the cost of dna sequencing. however, these techniques also raise important computational difficulties mostly due to the huge volume of data produced, but also because of some of their specific characteristics such as read length and sequencing errors. among the most critical problems is that of efficiently and accurately mapping reads to a reference genome in the context of re-sequencing projects.  we present an efficient method for the local alignment of pyrosequencing reads produced by the gs flx system against a reference sequence. our approach explores the characteristics of the data in these re-sequencing applications and uses state of the art indexing techniques combined with a flexible seed-based approach, leading to a fast and accurate algorithm which needs very little user parameterization. an evaluation performed using real and simulated data shows that our proposed method outperforms a number of mainstream tools on the quantity and quality of successful alignments, as well as on the execution time.  the proposed methodology was implemented in a software tool called tapyr--tool for the alignment of pyrosequencing reads--which is publicly available from http://www.tapyr.net. 