skewer: A fast and sensitive adapter trimmer for paired-end reads

[relipmoc]

Hello everyone,

We've implemented a novel tool named skewer for adapter trimming. It is aimed for preprocessing Illumina Paired-end/single-end reads at the moment. The main features are as follows:
* Allow full-length adapter sequence trimming for higher specificity;
* Allow indel errors when finding adapter sequence;
* Very fast: Internally it uses a novel local alignment algorithm that has not been published. In single thread mode, it can process a pair of compressed files, whose uncompressed sizes were about 12G bytes each, in about 30 minutes; it is even faster in multi-thread mode, but the speedup is limited due to the parallelism is only made in the sequence alignment part.
* Quality values aware. It evaluates alignments based on sequence qualities.
* Paired information aware. It is more accurate in case of processing paired-end reads.

If you are interested in using it, please download it from
https://sourceforge.net/projects/skewer/

Any feedback or feature requests are welcome!

[BhariD]

Hi,

I have to trim full-length adapter sequences with zero number of mismatches. I do not want to trim reads on any other criteria at this point.

I am using the following command line:
./skewer-0.1.99-linux-x86_64 -x ACACTCTTTCCCTACACGACGCTCTTCCGATCT -y GATCGGAAGAGCGGTTCA
GCAGGAATGCCGAG -r 0 -d 0 -o exact_trim_15 -t 8 read_1.fastq paired_read2.fastq

Log file includes:
Parameters used:
-- 3' end adapter sequence (-x): ACACTCTTTCCCTACACGACGCTCTTCCGATCT
-- paired 3' end adapter sequence (-y): GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG
-- maximum error ratio allowed (-r): 0.000
-- maximum indel error ratio allowed (-d): 0.000
-- minimum read length allowed after trimming (-l): 18
-- file format (-f): Sanger/Illumina 1.8+ FASTQ (auto detected)
-- number of concurrent threads (-t): 8
Tue Jan 14 02:18:14 2014 >> started

Tue Jan 14 02:19:33 2014 >> done (78.699s)
47656840 read pairs processed; of these:
0 ( 0.00%) short read pairs filtered out after trimming by size control
0 ( 0.00%) empty read pairs filtered out after trimming by size control
47656840 (100.00%) read pairs available; of these:
3202 ( 0.01%) trimmed read pairs available after processing
47653638 (99.99%) untrimmed read pairs available after processing

Length distribution of reads after trimming:
length count percentage
97 1 0.00%
98 4 0.00%
99 3197 0.01%
100 47653638 99.99%


My questions are:

1) The 3197 read pairs trimmed, given the input parameter settings, are they really trimmed just based on exact full-length adapter sequence match? any default parameter that I should be aware of?
2) What is the overlap length for adapter detection in paired-end mode? is it like initial 17 bp of the total length? Is there a way I can change this?
3) How can I change the number of mismatches to detect the adapter region in the read? Let's say if I want to allow only 2 mismatches (instead of zero mismatches) in the full-length adapter sequence?
4) How can I specify multiple adapter sequences for read 1 and read 2 data files?

I would appreciate your help! Thank you!

[relipmoc]

Thank you so much for your feedback!

Quick answers to your questions:
1) The searching process is based on exact full-length adapter sequence, but for the 3197 read pairs, only the last nucleotides were identified as the first nucleotides of corresponding adapter sequences. In current implementation, adapter sequence longer than 64 nt will be cut to 64 nt before processing.

2) There's no need to specify the overlap length in paired-end mode. The program knows how to do it correctly.

3) The program only provides a parameter of error ratio (by -r) and detect the most possible adapter location by a statistical scheme which takes into account the quality values. If you just want to specify the number of maximum allowed mismatches in the full-length adapter sequence, you can use fq2fa.sh to transfer the FASTQ files to FASTA files, and specify the maximum allowed error ratio (-r) as 2/33=0.06. For small RNA adapter trimming, it is something like the following command:
$ fq2fa.sh srnaReads.fq | skewer -x TCGTATGCCGTCTTCTGCTTGAAAAAAA -L 30 -r 0.06 -o trimmed -

4) For multiple adapter sequences, you just need to specify two FASTA files which contain adapter sequences, and input something like:
$ skewer -x adapters1.fa -y adapters2.fa flowcell1_lane7_pair1.fastq.gz flowcell1_lane7_pair2.fastq.gz

[BhariD]

Thank you for your prompt response!

I am sorry, I couldn't quite get the "In current implementation, adapter sequence longer than 64 nt will be cut to 64 nt before processing"? I don't think I have adapter more than 62 bp so then why its looking for last few nucleotides (3 I guess here?)?

[BhariD]

Also, what is the base quality value threshold used by the tool to be considered as a mismatch? in "3) The program only provides a parameter of error ratio (by -r) and detect the most possible adapter location by a statistical scheme which takes into account the quality values"

Thanks!

[relipmoc]

As I said, "there's no need to specify the overlap length in paired-end mode", actually there's no parameter or default parameter for the overlap length in paired-end mode.

The 64 nt statement is irrelevant to your question. I just misunderstood your question "any default parameter that I should be aware of". ^_^

"why its looking for last few nucleotides (3 I guess here?)". Unfortunately your guess is not the truth. It's by chance that you got this result.

Originally Posted by BhariD

Thank you for your prompt response!

I am sorry, I couldn't quite get the "In current implementation, adapter sequence longer than 64 nt will be cut to 64 nt before processing"? I don't think I have adapter more than 62 bp so then why its looking for last few nucleotides (3 I guess here?)?

[relipmoc]

There's no base quality value threshold. That's all integrated into the statistical scheme. Since we have not published the paper, I can not tell you the details at the moment. Sorry for that!

Originally Posted by BhariD

Also, what is the base quality value threshold used by the tool to be considered as a mismatch? in "3) The program only provides a parameter of error ratio (by -r) and detect the most possible adapter location by a statistical scheme which takes into account the quality values"

Thanks!

[roryk]

Hi relipmoc,

I have a couple of questions:

1) How does skewer handle partial matches? For example if I have a sequence that goes SEQUENCE-ADAPTER-BARCODE, and I just input ADAPTER, will I end up with SEQUENCE?

2) Why is this sequence not being trimmed? Does skewer only match the entire adapter sequence?

@test_truseq/1
CGATGATCAAGACCCAAGTGTGAGATTACGGAGATCGGAA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@test_truseq/2
CGATGATCAAGACCCAAGTGTGAGATTACTCAGATCGGAA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

~/tmp/skewer-0.1.104-linux-x86_64 -x AGATCGGAAGAG -y AGATCGGAAGAG test_cutadapt_1.fastq test_cutadapt_2.fastq

Thanks! I've been looking around for a faster trimmer and was hoping skewer would be the solution.

[frozenlyse]

Originally Posted by relipmoc

There's no base quality value threshold. That's all integrated into the statistical scheme. Since we have not published the paper, I can not tell you the details at the moment. Sorry for that!

That's not a good way to get people to use your software!

[relipmoc]

Originally Posted by roryk

1) How does skewer handle partial matches? For example if I have a sequence that goes SEQUENCE-ADAPTER-BARCODE, and I just input ADAPTER, will I end up with SEQUENCE?

The answer is Yes. However, if you want an improved specificity, you'd better use ADAPTER-BARCODE as the adapter sequence. Furthermore, if you want to demultiplex the reads, you can specify the --barcode option.

Originally Posted by roryk

2) ... Does skewer only match the entire adapter sequence?

The answer is No. skewer can detect partially matched adapter sequence at the 3' end (or 5' end if '-e 5' is specified).

Originally Posted by roryk

2) Why is this sequence not being trimmed? ...
@test_truseq/1
CGATGATCAAGACCCAAGTGTGAGATTACGGAGATCGGAA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@test_truseq/2
CGATGATCAAGACCCAAGTGTGAGATTACTCAGATCGGAA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

The sequences are not being trimmed because they are not as skewer expected. Is this from real data? Or could you explain why the paired sequences before adapter sequences are not reverse complementary to each other? Are they from mate-pair sequencing instead of paired-end sequencing?

Originally Posted by roryk

~/tmp/skewer-0.1.104-linux-x86_64 -x AGATCGGAAGAG -y AGATCGGAAGAG test_cutadapt_1.fastq test_cutadapt_2.fastq

There's no need to specify -y, if pair1 and pair2 share the same adapter sequence.

Originally Posted by roryk

Thanks! I've been looking around for a faster trimmer and was hoping skewer would be the solution.

My pleasure! Hope it will make your work easier.

[relipmoc]

Originally Posted by frozenlyse

That's not a good way to get people to use your software!

For those people who want to know technique details, I have to say sorry to them. However, you can't wait too long. I'll inform you once our submission is accepted. Thanks!

[hartmaier]

Originally Posted by relipmoc

For those people who want to know technique details, I have to say sorry to them. However, you can't wait too long. I'll inform you once our submission is accepted. Thanks!

Any chance you can release OSX binaries?

[kidaaaa]

I have a question:

How does skewer handle the situation only one read of a pair survives and the other one does not? I didn't see singletons in the output files.

[kidaaaa]

I have a question:

How does skewer handle the situation only one read of a pair survives and the other one does not? I didn't see singletons in the output files.

[relipmoc]

Originally Posted by kidaaaa

I have a question:

How does skewer handle the situation only one read of a pair survives and the other one does not? I didn't see singletons in the output files.

Thank you for your question! Now skewer only output those pairs that are concordantly trimmed or untouched.

[luc]

Hi replimoc,

I just found the answer in another thread:
http://seqanswers.com/forums/showthread.php?t=41976

My use case should be no problem, I guess.

(can skewer trim multiple potential contaminant sequences in the same run? - it seems it is focused on adapter pairs similar to trimmomatic's palindrome trimming mode?)

[roryk]

Thanks replimoc, with 0.1.114 skewer trims off the test sequences i posted correctly. Thanks!

[luc]

Hi replimoc,

it seems to me several sequences in the reverse reads are escaping removal when using the length threshold on paired reads. The filtering for length works fine for the forward reads.

Does the "paired information aware" trimming option work when providing a single (-x) adapter file containing several adapter sequences?

[relipmoc]

Hi Luc,
Thank you for your question! My answer goes as follows:

Originally Posted by luc

it seems to me several sequences in the reverse reads are escaping removal when using the length threshold on paired reads. The filtering for length works fine for the forward reads.

The length threshold -k does not influence the trimming result of paired-end data. Other length thresholds such as -l and -L do influence the results. Could you explain your case with more details?

Originally Posted by luc

Does the "paired information aware" trimming option work when providing a single (-x) adapter file containing several adapter sequences?

The answer is YES. However, you need to pay special attention on the trimming efficiency. The semantics of your case is to try n * n adapter combinations in adapter trimming, where n is the number of adapter sequences provided in the adapter file. If the adapter sequences share most of their content but differ in some region, e.g. 6-bp region for indexing, you may use degenerative characters in this region and specify one representative adapter sequence. For instance, if the content of the adapter file is:

>Index 1, ATCACG
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG
>Index 2, CGATGT
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG
>Index 3, TTAGGC
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACTTAGGCATCTCGTATGCCGTCTTCTGCTTG
>Index 4, TGACCA
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACTGACCAATCTCGTATGCCGTCTTCTGCTTG
>Index 5, ACAGTG
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACACAGTGATCTCGTATGCCGTCTTCTGCTTG
>Index 6, GCCAAT
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACGCCAATATCTCGTATGCCGTCTTCTGCTTG

you may specify -x TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG to improve trimming efficiency.

[relipmoc]

Originally Posted by roryk

Thanks replimoc, with 0.1.114 skewer trims off the test sequences i posted correctly. Thanks!

Hi roryk, thank you for your feedback!