Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads

General information

We simulated 10 million 100 bp + 100 bp PE Solexa reads from the Arabidopsis thaliana genome using ART, a NGS read simulator [15], with some revision on the source codes for simulating adapter-contaminated reads (http://sourceforge.net/projects/skewer/files/Simulator/). The trained profile was from the real sequencing data of A. thaliana where about 36% of the reads were contaminated with adapters. We compared Skewer with mainstream adapter trimmers that can handle PE reads as well as four representative adapter trimmers that can handle only SE reads.

To assess trimming quality, we defined the following metrics: FP (false positive) as the number of reads that were over-trimmed, either for trimming non-contaminant reads (false trimming), noted as F P_f t, or for over-trimming contaminant reads, noted as F P_o t; FN (false negative) as the number of reads that were under-trimmed, either for not trimming contaminant reads (false retaining), noted as F N_f r, or under-trimming contaminant reads, noted as F N_u t; and TN (true negative) as the number of untrimmed non-contaminant reads.

Primary result

Each method was run with its default parameters, except that the minimum output fragment length and the thread number (if applied) were set to 1.

The results of the runs listed in Table 2 show that Skewer outperformed all the mainstream tools in terms of mCC for both SE and PE trimming (0.9291 and 0.9989 respectively), although Skewer was only marginally better than Cutadapt and Scythe in SE trimming. Furthermore, Skewer was substantially faster (one times faster for SE and more than 12 times faster for PE trimming) than the tools that had comparative performances.

Trimmomatic and AlienTrimmer both used above 2G bytes peak memory. Most of the other trimmers used less than 50M bytes memory except SeqTrim, which used 115.7M bytes. Although Skewer did not have the least memory usage, its memory consumption (less than 35M bytes, see Table S2 of Additional file 1 for details) was far from a bottleneck on a 64-bit computer. In fact, Skewer uses additional memory to facilitate the processing of IUPAC (International Union of Pure and Applied Chemistry) characters and for parallel computing.

Scalability for parallel computing

We ran various adapter trimmers that support multi-threading to compare their scalability under a parallel computing environment. SeqTrim was excluded because it is too slow to gain comparative speed even with tens of threads. In addition, Cutadapt, AdapterRemoval, and other tools were not included because they currently lack a multi-threading function. In the eight threads case, for both the uncompressed and compressed inputs, Skewer achieved the highest speedup among the adapter trimmers tested (7.87 for uncompressed input, and 4.54 for compressed input) (see Table S2 of Additional file 1 for details).

Receiver operating characteristic (ROC) curves

ROC curves for various adapter trimmers under different stringencies were plotted and are shown in Figure 1 and Figure 2 (see Table S3 of Additional file 1 for details).

In high stringency (left) regions, both Trimmomatic and Cutadapt performed well in that they had low FPRs (false positive rates) and high TPRs (true positive rates); however, as the stringency decreases, both their performance degrade gradually (Figure 1). A similar trend was seen for TrimGalore where the ROC curve shifted to the upper-right region. This implies that these trimmers greedily picked up the first candidate that met the stringency rather than select the optimal one. The ROC curve for AlienTrimmer was similar to the above ones, but with a worse performance. FastX may adopt some optimization technique, however, its performance was worse than those of Trimmomatic and Cutadapt within all the stringency range. Other adapter trimmers showed advantages on a specific metric; e.g. AdapterRemoval, Flexbar, and EA-Tools were the most sensitive, while TagCleaner, Btrim, and Scythe were the most conservative. SeqTrim appeared only as a dot in the ROC curves plot, because it does not provide a stringency threshold. Skewer outperformed the other adapter trimmers in that it had the least FPR to gain a specific TPR, when TPR > 95%.

ROC curves for the adapter trimmers that are aware of PE information were plotted and are shown in Figure 2. Other trimmers that can process PE reads have worse ROC curves than corresponding ROC curves for processing SE reads since the second reads usually have lower sequencing qualities. From Figure 2, we can see that Skewer had a nearly perfect ROC curve close to the upper-left corner. For example, it achieved a TPR of 99.951% with a FPR of 0.001%.

sRNA sequencing data for Caenorhabditis elegans

A recently published real sRNA data set (short read archive [SRA:SRR014966]) [16], which includes 14,251,981 reads of small non-coding RNA (ncRNA) from C. elegans, was used to evaluate the adapter trimmers. Because it is hard to recover all the underlying sRNA fragments for sequencing, we aligned the trimmed reads to the reference genome and used delta of the number of uniquely aligned reads relative to the number of uniquely aligned raw reads, noted as TT (true trimming), as a substitute for true positive. We also used delta of number of non-uniquely aligned reads relative to the number of non-uniquely aligned raw reads, noted as FT (false trimming), as a substitute for false positive. The rationale was that correct-trimming tends to change unaligned fragments to uniquely aligned fragments (true positive), while over-trimming tends to change uniquely aligned fragments to non-uniquely aligned fragments (false positive). Note that these metrics tolerate tiny mistakes that can be rescued by the alignment software and are useful for practical evaluation.

To evaluate the performances under various trimming stringencies, all the tools were used to trim adapter sequences from the C. elegans data set using various trimming stringency. Next the processed reads were aligned to the C. elegans genome [17] (version 10) using Bowtie2 [18] (version 2.1.0). We then used the above metrics for final plotting, with higher FT representing lower stringency.

The results are presented in Table S4 and Table S5 of Additional file 2, and illustrated in Figure 3. AdapterRemoval and Flexbar exhibited similar performance curves, while AdapterRemoval was slightly better than Flexbar within all the tested stringency range; and TrimGalore and Cutadapt had similar curves, while TrimGalore was slightly better than Cutadapt at all the stringencies. Under high stringency, EA-Tools, Skewer, and TrimGalore shared the first rank in terms of low FT and high TT; Trimmomatic, AdapterRemoval, and Schythe were ranked second under middle stringency, middle low stringency, and low stringency respectively; and Skewer ranked first at all the stringencies.

Paired-end RNA sequencing data for Drosophila simulans

A real RNA-Seq data set with 27,005,344 pairs of 101 bp reads (short read archive [SRA:SRR330569]) from the gonads and carcasses of D. simulans was used to compare the performances of the adapter trimmers in trimming artificial contaminants from PE reads.

For the evaluation, we first used each of the tools that can deal with PE reads with default setting to trim adapters from the reads, with the exception that the minimum output fragment length was set to 20 and quality trimming was inhibited. We then used TopHat [19, 20] (version 2.0.10) to align the processed reads to the reference genome of D. simulans[21] (dsim revision 1.4). Finally the number of uniquely and concordantly aligned pairs was used as the performance metrics.

The results are presented in Table S6 of Additional file 3 and illustrated in Figure 4. Skewer outperformed the other adapter trimmers in terms of the number of uniquely and concordantly aligned pairs of the trimmed PE reads. Trimmomatic and AdapterRemoval, both of which performed well in processing the sRNA data, performed poorly in processing the long PE data. This finding implies that these tools may be tuned specifically for trimming adapters from sRNA data. Similarly, Btrim also performed less well with the PE data in this experiment. After investigating the processed data, we found that Btrim could recognize only the occurrence of the whole adapter sequence with a limited tolerance for insertions and deletions. It should be noted that all quality trimming was inhibited from these experiments to compare the adapter trimming performance alone. However, in real applications, quality trimming, which is outside of the scope of this paper, has been reported to improve the mapping rate and facilitate downstream data analysis [22].

Nextera long mate pair (LMP) data for Arabidopsis thaliana

A 5-kb insert size Nextera LMP library of A. thaliana Col-0 with 6,602,426 pairs of 251bp reads (European Nucleotide Archive [ENA:ERA264981]) and a 400-bp insert size Illumina HiSeq PE library of the same species with 17,341,797 pairs of 100 bp reads [ENA:SRR519624] were sequenced previously and used to demonstrate the utility of NextClip [14], a dedicated tool for trimming adapters from Nextera LMP libraries.

To compare Skewer with NextClip, we followed a validation procedure similar to the one described for NextClip [14]. Briefly, the LMP library was first trimmed using the adaptor trimmer. Then the trimmed LMP reads and the PE reads were de novo assembled using ABySS [23]. The time needed for the adapter trimming and the N50 lengths of the scaffolds were used as the metrics for the evaluation.

Notation

where δ_a,b=1 if a≠b, otherwise δ_a,b=0. When m i n(|A|,|B|)=0, ∥A,B∥ _{l e v} :=m a x(|A|,|B|).

k-difference problem

Given a sequence S=s₁s₂…s _n , a query pattern P=p₁p₂…p _m , and a threshold k(0≤k<m), search all substrings of S, noted as {P^′}, such that ∥P^′,P∥ _{l e v} ≤k.

Extended k-difference problem

Given a sequence S=s₁s₂…s _n , a query pattern P=p₁p₂…p _m , and a threshold e(0≤e<1,⌊n×e⌋=k), search all substrings of S, noted as {P^′}, such that ∥P^′,P∥ _{l e v} ≤k; and all suffixes of S, noted as {S^′}, such that ∃ a prefix of P, noted as P^′, and ∥S^′,P^′∥ _{l e v} ≤|S^′|×e

Lemma 1

In the dynamic programming matrix C for tackling the k-difference problem, the values of elements along each diagonal are monotonically non-decreasing.

The proof is provided in Additional file 5 Appendix A.

Theorem 1

All the matched elements of the query pattern and sequence are equal to their precursors in the diagonal and do not need to be updated in the dynamic programming process.

Proof.

This theorem is a direct consequence of Lemma 1 and the dynamic programming recurrence, when δ _{i j} =0.

In other words, only mismatched elements need to be updated in the dynamic programming process. If bit-vectors that denote mismatched positions of comparison between the adapter sequence and each of the four nucleotide characters are pre-computed and a bit-vector that marks all positions of the queue elements that exceed the k-difference constraint is maintained, then unnecessary computations in updating the column vector can be inhibited. This is the key point that led to the main improvement of our algorithm over Ukkonen’s algorithm.

As listed in Algorithm 1, the bit-masked k-difference matching algorithmhas the following characteristics:

This algorithm uses a queue of size m and several bit vectors of size ⌈m/w⌉, where w is the word length of the computer (for example w equals 64 for a 64-bit machine), and hence has a space of O(m). For each of the n characters in a target sequence, the character enters the queue once and exits from the queue at most once. For a random sequence, the expected size of the queue is O(k); hence, generally the algorithm has O(k n) expected time. However, because it is restricted by the bit-mask operations, each element in the queue usually updates at most k+1 times. Because bit operations are negligible compared with element update operations, this algorithm achieves O(k n) worst-case time in practice, which is better than the O(k n) expected time for Ukkonen’s algorithm.

Algorithm 1 can be improved further by avoiding all unnecessary updates through constant time bit operations within each iteration cycle of the target nucleotide. The basic principle is that when an element in a diagonal of the original dynamic programming matrix has a value that is derived from an adjacent diagonal (i.e. an indel occurs in the corresponding path), the score and associated index of the element will remain unchanged if the precursor remains unchanged.

Although the above improvement can reduce the theoretical time complexity from expected O(k n) to worst-case O(k n), experiments on large volumes of real data showed that the reduced element update operations did not compensate for the additional bit operations.

For the extended k-difference problem, an additional step bounded by O(m) is performed to check all the elements remained in the queue if no hit was found in previous steps.