GapFiller: a de novo assembly approach to fill the gap within paired reads

Definitions

Let Σ be an alphabet and Σ* be the set of the words from Σ. For every S ∈ Σ^∗ we will denote with |S| the number of characters of S and with S[p, . . ., p + l - 1] the sub-sequence of S starting in p ∈ {0, . . ., |S| - 1} and of length l ∈ {0, . . ., |S| - p}. We will refer to S[p, . . ., p + l - 1] as prefix if p = 0, suffix if p + l = |S|, and as the p- th character of S if l = 1, and we will simply write S[p].

In order to quickly identify overlaps between the contig under construction and the reads' tips, we use an approach closely related to the one presented in [20] based on an Hamming-aware hash function. The idea is that, by representing a string of length l as a base-|Σ| number, one can often replace expensive char-by-char comparison by fast integer (or bit-string) comparison. However, for practical values of l, the integers to be compared would not fit in a memory word. For this reason, as in the classical Karp-Rabin exact string matching algorithm [21], we can work with numbers modulo q considering equality modulo q only as an indication (necessary condition) that pairs of strings may be the same (i.e., operating with the strings' fingerprints). Policriti et al. in [22] proposed an extension of the approach by Karp and Rabin, introducing a technique to deal with mismatches, based on the idea of replacing simple fingerprints comparison with a more articulated test. In particular they noticed that, by choosing q to be a Mersenne (prime, when possible) number (i.e., q = 2 ^w - 1, for some w ∈ ℕ), to check whether two strings align against each other at a small Hamming distance can be implemented in average linear time.

where q is a (prime) number of the form q = 2 ^w - 1, for some w ∈ ℕ. The value f _H (S) is called the fingerprint of the sequence in S ∈ Σ* coded with s.

contains at least all the reads $r\in ℛ$ such that the hamming distance between r[x, . . ., x + b - 1] and $S\left[y,...,y+b-1\right]$ is not greater than k. False positives can be present but, as showed in [22], their amount is limited. On this ground the search for reads overlapping $\mathrm{ S}$ can be restricted to those belonging to $\mathcal{R}\left(S,x,y\right)$, for some x, y ∈ ℤ.

As far as GapFiller is concerned, we set k = 0 as default, meaning that we search for exact b-length substrings in the reads (i.e., $r\left[x,...,x+b-1\right]=S\left[y,...,y+b-1\right]$, for some x and y). As a consequence, better quality output will be obtained if we select a position x in r such that the average base quality is expected to be the highest possible. This point will be further discussed in the section specifically addressing data structures' design and implementation.

Dataset preparation

In order to avoid the generation of wrong contigs, it is of utmost importance to use only correct reads over the entire extension phase. Several tools are available to perform error correction on Illumina data using the so-called "read spectrum" (consider QUAKE [23], Hammer [24], and Allpaths [25] just to mention the most recent ones). Other tools discard reads or try to improve their reliability using quality information (rNA [20] and QSRA [16]).

Our approach, when we are given raw data, is to first trim (and possibly filter) the reads on the ground of quality information using a specific rNA option (refer to [20] for details), and to subsequently correct them with an error correction tool like QUAKE [23].

Another important way to assess a dataset's global quality is to plot the reads' k-mers distribution. This can be easily done using Jellyfish [26]. If the genome has been sequenced tens of times, then two peaks are expected: one in correspondence of the expected coverage and one in correspondence of coverage one. k-mers composing this second peak are likely to be sequencing errors. As a rule of thumb, a low number of k-mers occurring only once suggests that the dataset has a good global quality.

Contig extension

In the contig extension phase, each read is selected in a loop and used as seed in order to create a new contig. Once a seed read is selected, the suffix-prefix overlaps with other reads are computed and, if a sufficiently high level of global similarity is reached, they are clustered in a consensus string, which is subsequently used to perform further extensions. The procedure continues while some overlapping reads exist and the consensus string is highly representative of the clustered reads. If either one of the previous two conditions is not met, the extension phase stops, the current sequence is returned in output, and the loop continues.

Before the extension phase some parameters are set: the minimum overlap length L and the maximum shift Δ: an overlap between the current contig's suffix and the read's prefix is considered only if the overlap length l belongs to the interval [L, L + Δ].

GapFiller builds a cluster every time a contig is to be extended with the overlapping reads. In particular, GapFiller uses only those reads aligning against the contig's suffix with at most δ mismatches (where δ = δ(l) is a function of the overlap length l) and requires at least m reads in order to compute a consensus string. Notice that b ≤ L ≤ l holds, hence suffix-prefix overlaps might occur with more than k = 0 mismatches (see section Definitions).

Let $\mathrm{ ℛ}$ be the set of the input reads for GapFiller and $r_0\in ℛ$ be a seed read. At step i = 0 the current sequence is initialized with the seed S₀ := r₀. Denoting by S _i the current contig at the generic i-th step of the algorithm, the procedure to build S_i+1is described below:

Step1 Reads are selected according to their similarity with the current contig S _i (see Figure 2a). At this point, every read overlapping S _i for l ∈ [L, L + Δ] characters with at most δ mismatches is selected.

Step2 The reads are clustered and a consensus string is computed. Every character of the consensus string is assigned a flag indicating how it is representative of the reads from which it is built. More precisely, for every position j, GapFiller selects the most occurring character in the considered reads, and the majority consensus string C is computed (see Figure 2b). Depending on two parameters T₁ and T₂ such that T₁ <T₂, we say that a position j is non-represented, low-represented, or high-represented if the representation rate of the corresponding character in C is lower than T₁, lower than T₂, or higher than T₂, respectively.

Step3 The reads used to build the consensus C are filtered and trimmed, depending on the presence of low-represented and non-represented positions, respectively. The idea is that on low-represented positions we need a minimum percentage of reads matching the consensus string, and that on non-represented positions the extension is considered to be unsafe. Reads differing from C in correspondence of low-represented positions are discarded and the remaining ones are also trimmed if a non-represented position occurs (see Figure 2c).

Step4 A new consensus string C _new is computed, considering only the reads obtained at Step 3, and possibly the current contig is extended (see Figure 2d). The extension is done only if the number of reads is at least m and the consensus C _new exceeds S _i 's right end: in this case, a new contig S_i+1is built and the procedure restarts. Otherwise the algorithm stops and the contig S _i is returned.

The adopted strategy is aimed at either avoiding errors and overcoming the problems arising when GapFiller attempts to cluster reads that are different from each other. In the last part of this section we will discuss in more detail how the algorithm works. The reader who is not interested in the technical formalism might skip this part and move directly to the Subsection Stop criteria.

Stop criteria

The algorithm described in the previous section may potentially extend a contig for an arbitrarily large number of times, without checking any "global" properties of the current sequence. With our method the extension phase halts if at least one of the following conditions is met: (i) the available overlapping reads for the consensus C are less than m; (ii) the available overlapping reads for the new consensus C _new are less than m; (iii) contig's length exceeds the maximum length; (iv) the seed-mate has been found.

Let S _i be the contig obtained at the i-th step, starting from the seed read r₀. Criterion (i) applies when the consensus string C does not exceed the current contig. This means that there are no more than m - 1 overlapping reads, or that they are too short. In such a case, the contig produced is labelled as NO_MORE_EXTENSION.

Criterion (ii) applies when the consensus may have been produced as consequence of the presence of reads belonging to different genomic locations. More precisely, this situation is likely to appear when the consensus extension is "trying" to exit from a repeat. In this case, either too many reads are discarded (due to the presence of low-represented positions) or a significant trimming of them has been performed (as some non-represented positions occur far before the end of the consensus). In such a situation, the extension is halted and the contig is labelled as REPEAT_FOUND.

Criterion (iii) is satisfied as |S_i+1| >L_max, where L_max is fixed at the beginning of the algorithm and is usually set to the maximum insert size, plus a tolerance value. In such a situation, we could have been able to continue the extension but, however, we could not find the seed-mate. This suggests that the contig produced may be wrong or, at least, that it contains a high number of unreliable bases. When the maximum allowed length is exceeded, the computation is halted and the contig, labelled as LENGTH_EXCEED, is returned.

where M is the maximum number of mismatches allowed between ${\tilde{r}}_0$ and S _i . Inequality (7) is satisfied if and only if the mate is found in S _i at position p with no more than M mismatches. This control is performed on-the-fly and hence the positions already checked at the i-th step will not be re-checked. The mate-check criterion is used as a guarantee of correctness of the whole contig. This is in contrast to previous criteria, which are used to detect and prevent errors introduced in the extension phase. From this point of view, criteria (i) and (ii) can be seen as strictly local, since no information collected during previous steps is used. In this last case the contig returned is labelled as MATE_FOUND.

Data structures

In this section we will take a closer look at the data structures designed for our algorithm and at their implementation. GapFiller's core is the module working during the extension phase. At this point, we assume that the set $\mathrm{ ℛ}$ has already been trimmed and possibly filtered.

The basic idea is to pre-compute as much as possible of the useful information on the reads, in order to speed up the computation of the overlaps needed to perform the extension phase. Suppose that GapFiller is working at the (i+1)-th step of an extension, with i ≥ 0, and let S _i be the current contig. When constructing the consensus string C (see Figure 2a) we are always interested in obtaining overlaps between suffixes of S _i and prefixes of reads belonging to $\mathrm{ ℛ}$.

In order to compute overlaps, GapFiller employs a hashing schema based on the one implemented in rNA [20]; in particular, a data structure similar to the one proposed in [22] is built. A simplified schema of GapFiller's data structure is presented in Figure 3. The basic idea behind GapFiller is the possibility to obtain in a fast and efficient way the set of reads whose prefixes overlap a suffix of the partial contig under construction. Therefore we used the rNA hash function to find reads that are likely to overlap a suffix of S _i ; those reads are subsequently checked to see if they actually overlap S _i or not.

Obviously, all the data must be stored in the main memory, thus requiring a careful data structures' engineering. It is clear that, since overlaps between reads and the current contig can take place on both strands, reads must be stored together with their reverse complement.

The overall memory requirement for GapFiller depends on the size of HASHcounter and on the number of reads. As for rNA, a reasonable value for q is 2³⁰ -1. Such a number guarantees a reduction of the number of false positives (i.e., reads reported to align with the contig suffix, even though they do not overlap with it). As far as the number of reads is concerned, we can limit q, without loss of generality, to 2³¹: with state-of-the-art Illumina technology, such a number of reads represents approximately a 70× coverage of the human genome. An Illumina read of length 100 bp requires two memory locations in HASHvalues of 4 bytes each (31 bits to access array Reads and one bit to store the overlap orientation) and one entry in Reads of 9 bytes (7 bytes to store the read's numerical representation, one to store the mate position in Reads, and one more byte to store several useful informations about read status). In total the amount of memory required is $4q+2*4\left|ℛ\right|+9\left|ℛ\right|=4q+17\left|ℛ\right|$ bytes.

The reads' fingerprint is computed on a precise substring of length b (see (2)). As pointed out in section Definitions, the fingerprint of $r\in ℛ$ should be computed on the position x such that the (expected) average base quality is as high as possible and the substring r[x, . . ., x + b - 1] falls into the contigs' suffix, independently on the overlap length l. For these two reasons, having the Illumina error-profile in mind, we choose x = 0 if r is considered on its original strand, x = L - b if r has been reverse-complemented (see Figure 4).

In order to compute the overlaps between the current contig S _i and the reads, one has to compute the fingerprints of the substrings of length b starting from y, for every y ∈ {|S _i | - L - Δ, . . ., |S _i | - L} if original-stranded reads are searched, and for every y ∈ {|S _i | - Δ - b, . . ., |S _i | - b} if reverse-complemented ones are to be extracted. Let us indicate with s _y the fingerprint computed from S _i [y, . . ., y + b - 1] (see Figure 4). GapFiller uses this number to retrieve reads whose l-length prefix (l = |S _i | - y for original-stranded reads, l = |S _i | - y + L - b for reverse-complemented ones) is likely to match a substring of S _i close to the sequence's end. In particular GapFiller accesses all HASHvalues positions between HASHcounter[s _y ] and HASHcounter[s _y + 1] and, subsequently, accesses Reads to identify the set of candidate overlapping sequences $ℛ(S_i,l)$ (in Figure 3 GapFiller scans all positions between k and r - 1 of HASHvalues). Finally, the set $ℛ(S_i,l)$ is used to compute $\widehat{\mathcal{R}}\left(S_i,l\right)$, the set of real overlapping reads. This is done by checking all candidate reads singularly. Due to the fact that only a limited number of mismatches is allowed in this phase and that the employed hash function guarantees a low false positive rate, this step is extremely fast.

Dataset

As far as the experiments on real data are concerned, it is important to notice that the datasets provided by GAGE, together with the assembly results described in [9], represent the first available benchmarks that can be used to evaluate new instruments like GapFiller.

Using a specific rNA option, each simulated dataset was filtered to prune and trim reads on the basis of their quality information. For the real datasets, instead, we chose to use the Allpaths error-corrected reads, hence there was no need to trim them.

Design of experiments

We used simulated data in order to evaluate GapFiller's ability to correctly reconstruct the gap between two paired reads and to assess the reliability of the output classification (NO_MORE_EXTENSION, REPEAT_FOUND, LENGTH_EXCEED, and MATE_FOUND). In particular we used these datasets--easy to build and validate--to explore how coverage affects GapFiller's extension phase. Results on real datasets have been used instead to evaluate GapFiller's potential when its output is used as an input dataset for an assembler for long reads. However, the capability of producing correct contigs is a fundamental feature when GapFiller is used in this context.

GapFiller's performances rely on the choice of three crucial parameters: the minimum overlap length L, the slack Δ, and the length b of the substring on which the fingerprint is computed. We decided to set L = 50 and Δ = 40, as reads' length is approximately 100 bp for every library used for the experiments. The value of b identifies the length of a substring on which we (almost always) require an exact matching between read and contig (see Figure 4), due to the fact that the employed hash function has a low false-positives rate (see (2)). We set b = 20 because we observed that a greater value of b (i.e., close to L) dramatically prevents GapFiller to find even few-mismatch-affected overlaps.

The parameters T₁ and T₂, necessary to discern among high/low/non-represented positions in the consensus string (see Subsection Implementation-Contig extension), are set to T₁ = 0.6 and T₂ = 0.9. Recall that when a position in the consensus string has a representation rate lower than T₁, all the reads are trimmed on that position; instead, if the representation rate is lower than T₂, only the reads not matching the consensus string are dropped. The value of m, the minimum number of reads required in order to compute the consensus string, has always been set to 2. We chose not to let m depend upon coverage, since the number of reads after Step 3 strongly depends on the parameters used (say, T₁ and T₂).

We set the maximum length of a contig to be much greater than the expected mean insert size, i.e., 1800 bp for simulated data, 550 bp and 4500 bp for GAGE fragment and short jump libraries, respectively (see also Table 1 and Table 2).

We allowed for the presence of mismatches when looking for the seed-mate in the contig being constructed with parameter M. In all the performed experiments we set M = 10 (i.e., approximately 10% of the reads' length). This choice is justified by two reasons: the first one lies in the fact that the data simulated with SimSeq have a quite high amount of low-quality bases even far from the rightmost positions within the reads; the second one is that, on real datasets, lower values of M (e.g., 5 or 2) do not increase output quality. The value of δ, representing the maximum number of mismatches allowed when computing overlaps, depends on the overlap length l and was set to Ml / |r|, where |r| is the average read length.

Analysis

The post-processing phase of GapFiller's output is aimed at both quantitative and qualitative analysis. The first is focused on evaluating the amount of trusted contigs our tool is able to produce, the second on results' validation. The main goal is to compare the performances on different input datasets and coverages.

Due to their nature, experiments on simulated data allow to precisely estimate correctness by aligning a contig in the exact place where it is supposed to occur in the reference genome. More precisely, we used the Smith-Waterman alignment algorithm [29], assigning a score of 1 to a match, -1 to a mismatch, and -2 to an indel. For instance, let us consider a contig S generated by extending a seed read r, and suppose that r has been extracted from the genome G at position x, on the forward strand. To test its correctness, S is aligned against G[x, ..., x + |S| + g - 1], where g is the maximum number of allowed indels, depending on a user-defined threshold for the alignment score. We say that S is correctly aligned if and only if the ratio between the best alignment score of S against G[x, ..., x + |S| + g - 1] and |S| is at least 0.95 (for instance, we allow up to 5 mismatches, 1 indel and 1 mismatch, or 3 indels every 200 bp, on average). For this particular choice of the alignment score, g is fixed to be ⌈3|S|/200⌉.

Thanks to the presence of theoretical optimal assemblies for the two real datasets (see [9]) we evaluated the performances of GapFiller with respect to other assemblers. In particular, we extracted a set of contigs corresponding to a fixed coverage (10× for Staphylococcus aureus and 15× for Rhodobacter sphaeroides datasets, respectively) and assembled it with PHRAP [30], a well known Overlap-Layout-Consensus assembler. We produced a set of statistics representing the correctness of our assembly using the same scripts used in [9] and available for download at [28].