Efficient alignment of pyrosequencing reads for re-sequencing applications
© Fernandes et al; licensee BioMed Central Ltd. 2011
Received: 16 December 2010
Accepted: 16 May 2011
Published: 16 May 2011
Skip to main content
© Fernandes et al; licensee BioMed Central Ltd. 2011
Received: 16 December 2010
Accepted: 16 May 2011
Published: 16 May 2011
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 (454) 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.
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 (MPDS) technologies have recently emerged, including the Roche/454 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 [2, 3]. 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 454 Titanium platform to produce reads 400 bases long, some other studies may employ a SOLiD system set to produce short reads of 35 bases, and yet other projects may use the Illumina system to produce 2 × 75 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 (reads) 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 [5, 6]. Similarly, several methods have been proposed to address the problem of efficiently mapping MPDS reads to a reference sequence, like [7–12], to cite a few. As referred before, the sheer volume of data generated by MPDS technologies (to the order of hundreds of gigabases per run), 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 SSAHA2 , 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 (see Implementation), and BWA-SW , which uses a FM-index (see Implementation) 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 454 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 (reference) 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.
Recent research on indexes has focused on the fact that pointer representations require O(n log n) bits whereas the original text (the target genome, in our case) requires only n log σ bits, where σ is the alphabet size, e.g. 4 for DNA and 20 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 (BWT)  requiring only O(n log σ) 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 (rotations) 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[i] = t[SA[i] - 1], for i = 1, . . . , | t |. An example of a BWT is given in Figure 1(c).
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 1 ⋯ r m 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[1 ⋯ l] = r 1 ⋯ r l . In practice we obtain the equivalent of an interval of the BWT which contains the positions of the reference sequence g at which r[1 ⋯ l] occurs. Obviously, by maximality, r[1 ⋯ l + 1] does not occur in g. This happens because none of the occurrences of r[1 ⋯ l] is followed by r l+1 or, put another way, because there is a mismatch between r l+1 and the letter following each occurrence of r[1 ⋯ l] (or because it occurs at the very end of g). If r l+1 ≠ r l , then we set r[1 ⋯ l] as the first seed and proceed as above to find the next seed starting from position l + 2. If, however, we have r l+1 = r l , this means that the difference occurred in the middle of a homopolymer (contiguous subsequence of identical bases), most likely due to an insertion sequencing error. In this case, we set r[1 ⋯ l] 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. r l' ≠ r l' - 1. 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 1, . . . , s k be k substrings of r such that r = s 1 a 1 s 2 a 2 ⋯ a k-1 s k a k , where, for i = 1, . . . , k, s i denote the seeds and a i denote the substrings in between them. For each i = 1, . . . , k, let o i ≥ 0 be the number of exact occurrences of s i in g. We then have ways to choose a set containing one occurrence for each of the k substrings. Let p = (p 1, . . . , p k ) be one of such o tuples of distinct seed occurrences. If, for i = 1, . . . , k -1, we have p i ≤ p i+1 and |a i | ∈ i p i+1 - (p i + |s i |) |a i | + ∈ i for some given ∈ i ≥ 0, 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 (s 1, . . . , s k ) 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.
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 [20, 21], 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 (formally, a homopolymer can consist of a single base), 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 (substitutions) 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 p sub, p ins and p del 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 p sub, an insertion with probability p ins, a deletion with probability p del or it can be correctly sequenced with probability 1 - (p sub + p ins + p del). 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[a, b] indicates the probability for a to be miscalled as (replaced by) 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 γ > 0 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 , SSAHA2 , 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.
Real data sets
Reference genome Source
Reference genome size
SRA accession and species
Average read length
ATCC 700669 [GenBank: FM211187]
○ SRR001327 S. pneumoniae CDC1873-00
○ SRR001328 S. pneumoniae SP195
○ SRR001329 S. pneumoniae CDC0288-04
E. coli 0127:H6 E2348/69 [GenBank: FM180568.1]
○ SRR000868 E. coli K-12
○ SRR000870 E. coli K-12
○ SRR031369 E. coli ETEC WS3080A
○ SRR031370 E. coli ETEC TW03576
P. falciparum 3D7 PlasmoDB rel 7.0
○ SRR006911 P. falciparum 3D7
○ SRR006912 P. falciparum 3D7
○ SRR006913 P. falciparum 3D7
○ SRR006914 P. falciparum 3D7
○ SRR006915 P. falciparum 3D7
○ SRR022943 C. elegans Lynch MA41 mutation-accumulation line derived from N2.
D. pseudoobscura FlyBase rel. 2.14
○ SRR003807 D. pseudoobscura Flagstaff 1993
○ SRR014458 D. pseudoobscura bogotana ER (white)
○ SRR014459 D. pseudoobscura bogotana ER (white)
○ SRR014460 D. miranda strain Mather 1993
H. sapiens Chr. 15 ENSEMBL ver. GRCh37
○ SRR014420 Human individual NA15510
○ SRR014421 Human individual NA15510
○ SRR014422 Human individual NA15510
○ SRR014423 Human individual NA15510
○ SRR014424 Human individual NA15510
○ SRR014425 Human individual NA15510
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.
Experimental results with real 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 = 300,000 synthetic reads from the 250 Mbp sequence of the human chromosome 1 [GenBank: NC_000001.10]. 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 HS1, was generated with the read generator parameters set as μ l = 250, σ l = 50, p ins = p del = p sub = 0.01, ω = 10, γ ins = γ del = 3, and equiprobable substitution rates M[a, b] = 1/3, 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, HS2, we increased the error levels of each kind to 5% by setting p ins = p del = p sub = 0.05, and for the third data set, we added a considerable amount of noise by setting p ins = p del = p sub = 0.10 (only 70% of chance for each homopolymer to be sequenced correctly). 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.
Experimental results with synthetic data
Memory requirements for real data
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.
Project name: TAPyR--Tool for the Alignment of Pyrosequencing Reads
Project home page: http://www.tapyr.net
Operating system(s): multiple (requires a C compiler only)
Programming language: C
Other requirements: none (see Results section for an idea of memory usage)
License: GNU GPL
Restrictions to use by non-academics: none additional
We thank Ana Tereza Vasconcelos and her staff from LNCC, Brazil, for sharing their expertise with the GS FLX data and for granting us access to their platform.
Funding This work was supported by FCT (INESC-ID multiannual funding) through the PIDDAC Program funds, by the FCT PhD grant SFRH/BD/45586/2008 (FF), and by the PTDC program (project PTDC/EIA-EIA/112283/2009).
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.