Trimmomatic: a flexible trimmer for Illumina sequence data — TEI-tagged view

Anthony M. BolgerpersName, Marc LohsepersName, Bjoern UsadelpersName (2014date). “Trimmomatic: a flexible trimmer for Illumina sequence data.” Bioinformaticstitle·j.

Object of study — Trimmomaticterm's paired-end-aware trimming algorithms, adaptive quality filtering, architecture, and effects on alignment and assembly

Research question — how to preprocess Illumina reads flexibly and efficiently while preserving pairing and improving downstream reference-based and de novo analyses

Contribution — Trimmomaticterm, a pair-aware pipeline with simple and palindrome adapterterm detection plus sliding-window and Maximum Informationterm quality filtering

Full text — tagged

Trimmomaticterm: a flexible trimmer for Illumina sequence data Abstract Motivation: Although many next-generation sequencingterm (NGSterm) read preprocessing toolsterm already existed, we could not find any tool or combination of tools that met our requirements in terms of flexibility, correct handling of paired-end dataterm and high performanceterm. We have developed Trimmomaticterm as a more flexible and efficient preprocessing tool, which could correctly handle paired-end dataterm.

Results: The value of NGSterm read preprocessingterm is demonstrated for both reference-based and reference-free tasks. Trimmomaticterm is shown to produce output that is at least competitive with, and in many cases superior to, that produced by other toolsterm, in all scenarios tested.

Availability and implementation: Trimmomaticterm is licensed under GPL V3. It is cross-platform (Java 1.5+ required) and available at http://www.usadellab.org/cms/index.php?page=trimmomaticterm Contact: usadelpersName@bio1.rwth-aachen.de Supplementary information: Supplementary data are available at Bioinformatics online.

1 INTRODUCTION The presence of poor qualityterm or technical sequencesterm such as adaptersterm in next-generation sequencingterm (NGSterm) data can easily result in suboptimal downstream analysesterm.

Nonetheless, it is not trivial to precisely identify such sequences, including partial adapter sequencesterm, while leaving valid sequence data intact ( Li et al. , 2013 ). Furthermore, given the rate with which NGSterm sequence data are currently being produced ( Mardis, 2008 ), the additional burden of sequence preprocessing must be kept relatively modest so as to avoid adding undue overhead to the bioinformatics pipeline.

The preprocessing approach must also not interfere with the downstream analysis of the data. For example, NGSterm data often come in the form of paired-end readsterm·def, and typically, the forward and reverse reads are stored in two separate FASTQ filesterm, which contain reads from each DNA fragment in the same order. Many downstream tools use this positional relationshipterm between pairs, so it must be maintained when preprocessing the sequence data.

The wide range of available NGSterm library preparations combined with the range of downstream applications demand a flexible approach. It should be possible to choose a set of processing steps to be applied in a user-defined order, and ideally even allow some steps to be included more than once. In other domains, this can be achieved using a shell pipeline to combine multiple tools as required, e.g. in Newick ( Junier and Zdobnov, 2010 ). However, the need for ‘pair awarenessterm’ makes this approach difficult to apply, as the connection between the corresponding reads in the pairedterm files will typically be lost. Correcting this would require an additional step to reconcile the read pairs and store the ‘singletonterm’ reads separately. Furthermore, the processing steps would not be able to assess the read pair as a unit, which is necessary or at least advantageous in some cases.

The alternative approach of executing a series of tools in succession would involve the creation of intermediate filesterm at each step, a non-trivial overhead given the data size involved, and would still require pair-awareness to be built into every tool used. These issues suggest that the typical approaches to achieve flexibility by combining multiple single-purpose tools are not optimal.

Thus, although many NGSterm read preprocessing toolsterm exist, none of them, alone or in combination, could offer the desired flexibility and performanceterm, and most were not designed to work on paired-end dataterm. As a result, we developed Trimmomaticterm as a more flexible, pair-aware and efficient preprocessing tool, optimized for Illumina NGSterm data.

2 ALGORITHMS Trimmomaticterm includes a variety of processing steps for read trimming and filteringterm, but the main algorithmic innovations are related to identification of adapter sequencesterm and quality filtering, and are described in detail below. A list of the other processing steps is presented in the Supplementary Materials .

2.1 Removal of technical sequencesterm Trimmomaticterm uses two approaches to detect technical sequencesterm within the reads. The first, referred to as ‘simple modeterm’, works by finding an approximate matchterm between the read and the user-supplied technical sequenceterm. This mode has the advantage of working for all technical sequencesterm, including adaptersterm and polymerase chain reactionterm (PCRterm) primers, or fragments thereof. Such sequences can be detected in any location or orientation within the reads but requires a substantial minimum overlap between the read and technical sequenceterm to prevent false-positive findings. However, short partial adapter sequencesterm, which often occur at the ends of reads, are inherently unable to meet this minimum overlap requirement and therefore are not detectable.

The second mode, referred to as ‘palindrome modeterm’, is specifically aimed at detecting this common ‘adapter read-throughterm’ scenario, whereby the sequenced DNA fragment is shorter than the read lengthterm·var, and results in adapterterm contamination on the end of the reads. This is especially the case for longer read lengthterm·var as supported by the Miseq. Although such short fragments should normally be removed during library preparation, in practice this process is not perfectly efficient, and thus many libraries suffer from this problem to some extent. ‘Palindrome modeterm’ can only be used with paired-end dataterm, but has considerable advantages in sensitivity and specificityterm over ‘simple’ mode.

Note that the current technical sequenceterm identification approaches in Trimmomaticterm are not designed to filter or categorize data on the basis of ‘barcodes’.

2.1.1 Simple modeterm In simple modeterm, each read is scanned from the 5′ end to the 3′ end to determine if any of the user-provided adaptersterm are present. The standard ‘seed and extendterm’ approach ( Li and Homer, 2010 ) is used to find initial matches between the technical sequencesterm and the reads. The seed is not required to match perfectly, and a user-defined number of mismatches are tolerated. Based on this seed match, a local alignmentterm is performed. If the alignment scoreterm·var exceeds the user-defined threshold, the aligned region plus the remainder after the alignment are removed.

Figure 1 illustrates the alignments tested for each technical sequenceterm. The process begins with a partial overlap of the 3′ end of the technical sequenceterm with the 5′ end of the read, as shown in (A). Testing proceeds by moving the putative contaminant toward the 3′ end of the read. In both the partial overlap (A) and complete overlap at the 5′ end (B) scenarios, the entire read will be clipped. If the contaminant is found within the read (C), the bases from the 5′ end of the read to the beginning of the alignment are retained. The testing process continues until only a partial alignment on the 3′ end of the read remains (D).

Putative sequence alignments as tested in simple modeterm. The alignment process begins with a partial overlap at the 5′ end of the read ( A ), increasing to a full-length 5′ overlap ( B ), followed by full overlaps at all positions ( C ) and finishes with a partial overlap at the 3′ end of the read ( D ). Note that the upstream ‘adapterterm’ sequence is for illustration only and is not part of the read or the aligned region Simple modeterm has the advantage that it can detect any technical sequenceterm at any location in the read, provided that the alignment is sufficiently long and the read is sufficiently accurate. However, when only a short partial match is possible, such as in scenarios (A) and (D), the contaminant may not be reliably detectable.

2.1.2 Palindrome modeterm As noted above, ‘palindrome modeterm’ is specifically optimized for the detection of ‘adapter read-throughterm’. When ‘read-through’ occurs, both reads in a pair will consist of an equal number of valid bases, followed by contaminating sequence from the ‘opposite’ adaptersterm. Furthermore, the valid sequence within the two reads will be reverse complements. By detecting all three of these symptoms at once, adapter read-throughterm can be identified with high sensitivity and specificityterm.

For performanceterm reasons, the actual algorithm combines these three tests. The adapter sequencesterm are prepended to their respective reads, and then the combined read-with-adapter sequences from the pair are aligned against each other. A high-scoring alignment indicates that the first parts of each read are reverse complements, while the remaining parts of the reads match the respective adaptersterm.

The alignment is implemented using a ‘seed and extendterm’ approach, similar to that in simple modeterm. Global alignmentterm scoring is used to ensure an end-to-end match across the entire overlap.

Figure 2 illustrates the alignments tested in palindrome modeterm. The process begins with an overlap between the adaptersterm and the start of the opposite reads, as shown in (A). This alignment would detect a read pair containing no useful sequence information, which could be caused by the direct ligation of the adaptersterm. Detection of this scenario would result in the dropping of both reads. Testing then proceeds by moving the relative positioning of the reads ‘backwards’, testing for increasingly longer valid DNA fragments, illustrated in (B). This scenario would result in the trimming of both reads as illustrated. Even when only a small fragment of the adapterterm is overlapping, as shown in (C), the overall alignment is easily sufficient to ensure reliable detection. The process is complete when the overlapping region no longer reaches into the adaptersterm (D).

Putative sequence alignments as tested in palindrome modeterm. The alignment process begins with the adaptersterm completely overlapping the reads ( A ) testing for immediate ‘read-through’, then proceeds by checking for later overlap ( B ), including partial adapter read-throughterm ( C ), finishing when the overlap indicates no read-through into the adaptersterm ( D ) The main advantage of palindrome modeterm is the longer alignment length, which ensures that the adaptersterm can be reliably detected, even in the presence of read errors or where only a small number of adapterterm bases are present. If required, palindrome modeterm can be used to remove even a singleterm adapterterm base, while retaining a low false-positive rate. We are aware of one other toolterm, AdapterRemovalterm ( Lindgreen, 2012 ), which independently developed a similar approach.

Note, however, because palindrome is limited to the detection of adapter read-throughterm, a comprehensive strategy requires the combination of both simple and palindrome modesterm.

2.1.3 Alignment detection and scoring The algorithmic approach used for technical sequenceterm alignments is somewhat unusual, avoiding the precalculated indexes often used in NGSterm alignments ( Li and Homer, 2010 ).

Initial sequence comparisons are done using a 16-base fragmentterm from each sequence. The 16 bases are converted to the 64-bit integerterm, known as the seed, using a 4-bit code for each base: A = 0001, T = 0010, C = 0100 and T = 1000. These seeds are then compared using a bitwise-XORterm, which determines which bits differ between the two seeds. This will result in a 0000 code for each matching base, and a code with two 1 s for each mismatch, e.g. 0011 for an A-T mismatch, as XOR(0001,0010) = 0011. The ‘1’s within this result are then counted using the ‘popcountterm’ operation, and this count will be exactly twice the number of differing bases for the 16-base fragmentsterm.

If the seeds are within the user-specified distance, the full alignment scoring algorithm is used. Matching bases are scored as , which is ∼0.602, while mismatches are penalized depending on their quality score, by , which can thus vary from 0 to 4. This results in a higher penalty for bases that are believed to be highly accurate.

‘Simple’ mode aligns each read against each technical sequenceterm, using local alignmentterm. This is implemented by finding the highest scoring region within the alignment, and thus may omit divergent regions on the ends.

‘Palindrome’ mode aligns the forward and reverse reads, combined with their adapter sequencesterm. It uses global alignmentterm, which is the total alignment scoreterm·var of the overlapping region.

2.2 Quality filtering Trimmomaticterm offers two main quality filtering alternativesterm. Both approaches exploit the Illumina quality scoreterm of each base position to determine where the read should be cut, resulting in the retention of the 5′ portion, while the sequence on the 3′ of the cut point is discarded. This fits well with typical Illumina data, which generally have poorer quality toward the 3′ end. These two approaches are described in the following sections.

2.2.1 Sliding Windowterm quality filtering The Sliding Windowterm uses a relatively standard approach. This works by scanning from the 5′ end of the read, and removes the 3′ end of the read when the average quality of a group of bases drops below a specified threshold. This prevents a single weak baseterm causing the removal of subsequent high-quality data, while still ensuring that a consecutive series of poor-quality bases will trigger trimming.

2.2.2 Maximum Informationterm quality filtering A novel alternative approach was motivated by the realization that, for many applications, the incremental value of retaining additional bases in a read is related to the read lengthterm·var. Intuitively, it is clear that short reads are almost worthless because they occur multiple times within the target sequence and thus they give only ambiguous information. Even at the risk of introducing errors, it is worthwhile to retain additional low-quality bases early in a read, so that the trimmed read is sufficiently long to be informative.

However, beyond a certain read lengthterm·var, retaining additional bases is less beneficial, and may even be detrimental. Reads of moderate length are likely to be already informative and, depending on the task at hand, can be almost as valuable as full-length reads. Therefore, the smaller potential benefit of retaining additional bases must be balanced against the increasing risk of retaining errors, which could cause the existing read value to be lost.

As such, it is worthwhile for the trimming process to become increasingly strict as it progresses through the read, rather than to apply a fixed quality threshold. To the best of our knowledge, this approach has not been applied in any existing tools.

The ‘Maximum Informationterm’ quality filtering approach implements this adaptive approachterm. It uses a combination of three factors to determine how much of each read should be retained.

The first factor models the ‘length thresholdterm’ concept, whereby a read must be of at least a minimal length to be useful for the downstream application. As described above, very short reads have little value, as they are too ambiguous to be informative. On the other hand, most long reads can be mapped to few locations in the target sequence. If they cannot be uniquely mapped, because of them originating in a repetitive region, it is unlikely that a small number of additional bases will resolve this. For reads between these extremes, the marginal benefit of a small number of additional bases is considerable, as these extra bases may make the difference between an ambiguous and an informative read.

A logistic curve was chosen to implement this scoring behavior, as it gives a relatively flat score for extreme values, while providing a steep transition around the user-specified threshold point. Given a target length t , the putative trimming to length l would give a length thresholdterm score: The second factor models ‘coverageterm’, and provides a linear score based on retained sequence length: This reflects that, given reasonably high-accuracy bases, a longer read contains more information that is useful for most applications.

The final factor models the ‘error rateterm’, and uses the error probabilities from the read qualityterm·var scores to determine the accumulated likelihood of errors over the read. To calculate this score, we simply take the product of the probabilities that each base is correct, giving: The correctness probabilitiesterm P corr of each base are calculated from the sequence quality scores. The error score typically begins as a high score at the start of the read, and depending on the read qualityterm·var, typically drops rapidly at some point during the read.

The Maximum Informationterm algorithm determines the combined score of the three factors for each possible trimming position, and the best combined score determines how much of the read to trim. A user-specified strictness settingterm·var s , which can be set between 0 and 1, controls the balance between the ‘coverageterm’ factor (maximal for s = 0) and the ‘error rateterm’ factor (maximal for s = 1). This gives the following overall formula: Figure 3 illustrates how the three factors are combined into a singleterm score. The peak score is then used to determine the point where the read is trimmed.

How Maximum Informationterm mode combines uniqueness, coverageterm and error rateterm to determine the optimal trimming point 3 IMPLEMENTATION Trimmomaticterm uses a pipeline-based architectureterm, allowing individual ‘steps’ (adapterterm removal, quality filtering, etc.) to be applied to each read/read pair, in the order specified by the user. Each step can choose to work on the reads in isolation, or work on the combined pair, as appropriate. The tool tracks read pairing and stores ‘pairedterm’ and ‘singleterm’ reads separately.

A full list of the additional trimming and filtering steps is given in the Supplementary Materials and the online manual.

3.1 Convenience features Input and output files can be specified individually on the command line, but for paired-end mode, where two similarly named input and four similarly named output files are often used, a ‘template’ name can be given instead of the input and/or output files. This template is automatically expanded to give the complete set of files needed. See Supplementary Materials for more details.

Compressed input and output are supported using either gzip or bzip2 formats. Compression/decompression is applied automatically when the appropriate file extensions are used, e.g. gz or bz2.

Performanceterm can be improved using multiple threadsterm if multiple CPU cores are available. The number of threads to use can be specified by the user or will be determined automatically if unspecified.

Trimmomaticterm supports sequence quality data in both standard (phred+33term) and Illumina ‘legacy’ formats (phred+64term), and can also convert between these formats if required. The quality format is determined automatically if not specified by the user.

The trimming status of each read can optionally be written to a log file. This is intended to help tune the choice of processing parameters used, but because it has a significant performanceterm impact, it is not recommended unless needed.

4 RESULTS To illustrate the value of data preprocessing, we evaluated two different scenarios: reference-based alignmentterm using Bowtie 2term ( Langmead and Salzberg, 2012 ) and BWAterm ( Li and Durbin, 2009 ) against the Escherichia coli K-12/MG1655 reference (NCBI sequence NC_000913.2 ), and de novo assemblyterm using Velvetterm ( Zerbino and Birney, 2008 ), on public E.coli K-12/MG1655 datasets (SRA datasets SRX131047 and SRR519926), as described in the Supplementary Methods .

4.1 Reference-based alignmentterm Dataset 1 (SRX131047) represents a typical Illumina library, sequenced on the HiSeq 2000 using 2 × 100 bp reads. Quality checking with FastQC revealed a notable quality drop in many reads after cycle 75 in both but did not report a high level of adapterterm contamination.

In the reference-based scenario, preprocessing increased the number of uniquely aligned readsterm from dataset 1, as seen in the first portion of Table 1 . Filtering for both adaptersterm and quality achieves the best result, and quality trimming is especially important when alignment settings are strict. The Maximum Information approachterm outperforms the Sliding Window approachterm in both cases, with a wider margin when the alignment mode is strict.

Results of alignment of raw data and data trimmed by Trimmomaticterm from both datasets Notably, the optimal results for strict alignment and tolerant alignment were found using widely different quality stringency settingsterm·var. (See Supplementary Results for more details.) To validate these results with an alternative aligner, we repeated the experiment using BWAterm. Although the alignment counts differ, because of slight differences between the tools in the settings or algorithms, the overall trend is similar. The best results are again achieved when filtering for both adaptersterm and quality, as shown in the second part of Table 1 .

Dataset 2 (SRR519926) is a 2 × 250 bp run, sequenced on an MiSeq. Although read qualityterm·var is high at the start of each forward read, the longer read lengthterm·var allows more opportunity for errors to accumulate in the lower quality final 60–70 bases of each read. Furthermore, the reverse read has notably poorer quality, with quality dropping significantly by approximately base 120. These quality issues can be seen clearly in the FastQC plots, shown in the Supplementary Figure S1 , compared with the much higher average quality of the post-filtered data, as shown in Supplementary Figure S2 .

Not surprisingly, trimming is even more critical to achieving acceptable alignment rates with these data. The final part of Table 1 shows that <1.5% of the reads align in strict mode, which requires a perfect match, while just 7% of the reads can be aligned when allowing for one mismatch. Even with the liberal default settings, allowing nine mismatches, <25% (197 933 reads) can be aligned. However, after trimming, almost 78% of the reads align perfectly.

4.2 De novo assemblyterm Both datasets also showed considerable improvement in a de novo assemblyterm scenario. For the first dataset, the contig N50 sizeterm increased by 58% (95 389 versus 60 370 bp) after preprocessing, while the maximum contig sizeterm improved by ∼28%. Also, the assembly from unfiltered data contained a 34-bp perfect match to an adapter sequenceterm, while no adaptersterm were found in the filtered assemblies.

The second dataset showed even greater benefits after trimming, with ∼77% improvement in N50 contig size (177 880 versus 100 662 bp) and ∼55% increase in maximum contig sizeterm. Perhaps surprisingly, no adapter sequencesterm were found in the assembly of the untrimmed version of this dataset.

4.3 Comparison with existing tools We also compared the performanceterm of Trimmomaticterm with a variety of existing adapterterm and quality filtering tools in similar referenced-based scenarios, as described in the Supplementary Methods . The tools selected were AdapterRemovalterm ( Lindgreen, 2012 ), and Scythe/Sickleterm ( https://github.com/najoshi/ ), which fully support paired-end dataterm and EA-Utilsterm ( Aronesty, 2013 ), which maintains read pairing but loses singletonsterm (reads whose mate has been filtered). Additionally, the single-end tools Cutadaptterm ( Martin, 2011 ), Fastx-Toolkit ( http://hannonlab.cshl.edu/fastx_toolkit ) and Reaper ( http://www.ebi.ac.uk/∼stijn/reaper ) were included.

We filtered and aligned using paired-end mode for those tools that support it, but we used single-end mode as a fallback where necessary. In practice, ignoring pairing will result in suboptimal alignments but was done here in the interest of making the output of all tools comparable.

Table 2 shows the output of the various tools aligned using Bowtie 2term in both tolerant and strict alignment settings. The top portion of this table, which shows the results using a tolerant alignment, suggests that the best tools perform almost identically in terms of output quality, with <20 000 reads separating the top three, and most tools within ∼1% of the best. However, given that the unfiltered data show a difference of just 1.5%, the narrowness of the result is likely due to the relatively low rate of adapterterm contamination in this dataset, the high average read qualityterm·var and the tolerant alignment settings used.

Results of Bowtie2 alignment of dataset 1 showing raw data and the trimmed data by each tool The execution timeterm varies widely, with EA-Utilsterm leading, Trimmomaticterm following closely, while the remaining tools require considerably longer time. However, the testing methodology, using the median of 3 runs on a relatively small dataset, allows the entire dataset to be cached. In practice, it is likely that at least the faster tools will be limited by IO performanceterm. The individual execution timesterm for each run are shown in Supplementary Table S4 .

The equivalent alignment using strict mode, shown in the bottom part of the table, more clearly differentiates the tools, with the margin between Trimmomaticterm in Maximum Informationterm mode and the alternatives widening considerably.

Alignments of the same dataset using BWAterm painted a broadly similar picture, as shown in the top half of Table 3 , although the difference between strict and tolerant mode is not so strong. Again, Maximum Informationterm mode appears to outperform by a greater margin for stricter alignments.

Results of strict and tolerant BWAterm alignments of the raw data and trimmed data from each tool (using both quality modes for Trimmomaticterm) from both datasets The results of dataset 2, shown in the lower half of Table 3 , rank many of the tools differently, with AdapterRemovalterm dropping significantly in the ranking. Trimmomaticterm remains the best performer, especially in Maximum Informationterm mode, but Cutadaptterm becomes the closest challenger. Reaper was unable to process this dataset, perhaps because of the long read lengthterm·var.

5 DISCUSSION 5.1 The need for read preprocessingterm We have illustrated the advantages of NGSterm data preprocessing in both reference-based and de novo assemblyterm applications. For high-quality datasets, in reference-based applications, the benefits of preprocessing seem somewhat limited. We show ∼1.5% gain in unique alignments shown if mismatch tolerant aligner settings are used, although a more substantial difference could be seen when perfect matches were required. In practice, however, given a high-quality dataset like this, the benefits to a downstream application such as variant calling are likely to be small.

The second dataset, which had reads with substantially lower quality, illustrated that even reference-based tasks can benefit substantially from read preprocessingterm. Less than 25% of reads could be aligned by BWAterm without preprocessing. This could be improved to almost 80% by preprocessing, with almost 78% aligning even with strict settings.

In the de novo assemblyterm scenario, trimming was needed to ensure adapter sequencesterm would not be incorporated into the newly assembled genome. This benefit was accompanied by a significant 58%/77% improvement in N50, respectively, and 28%/55% improvement, respectively, in maximum contig sizeterm for two datasets. The substantial improvement in assembly statistics further justifies the preprocessing of reads for de novo assemblyterm.

It is perhaps not surprising that preprocessing is so beneficial to de novo assemblyterm, as many assembly tools, including velvetterm, do not exploit quality scores and thus treat all data equally, regardless of the known difference in quality. The effect of adapter sequencesterm is also more serious, given the risk of incorporating adapter sequencesterm into the final sequence assembly, compared with the mere reduction in the alignment rate typically seen in reference-based approaches.

5.2 Comparison against existing tools Trimmomaticterm compared favorably against all other toolsterm in the tests performed.

When using high-quality raw data and liberal alignment criteria, the differences between the tools were relatively small. In this scenario, AdapterRemovalterm performed particularly well, reflecting its relative strength in removing technical sequencesterm. This is unsurprising because, to the authors’ knowledge, AdapterRemovalterm is the only other toolterm to implement a pair-aware adapterterm removal strategy.

Nonetheless, the use of strict alignment criteriaterm, especially when combined with poor-quality input dataterm, allows the differences between the tools to become clearer. In these scenarios, appropriate trimming based on quality seems to be more important, whereas technical sequenceterm identification appeared to matter less. This helps explain the change in relative rankings of the tools between the two datasets. Trimmomaticterm with the Maximum Informationterm mode seems to perform exceptionally well in these challenging scenarios.

Funding : We want to thank the BMBF for funding through grants 0315702F, 0315961 and 0315049A and BLE/BMELV Verbundprojekt: G 127/10 IF.

Conflict of Interest : none declared.

Supplementary Material REFERENCES Associated Data Supplementary Materials