[earonesty]

Well, I've never tried this with short adapters like that. Looks like the thing I did to fix your case made other things worse. I need to put the code back the way it was. The original behavior was more correct. The problem is the use of a very short adapter sequence, which is shorter than the scan length the sytem uses (15 bases). Wtih something that short, it will find them everywhere,and it isn't a realistic test. I panicked and gave you a version that was a lot worse.

Solution: Use the original version ( i put it back... v 148), and use a sequence longer than 15 BP for your test adapter.

( If you really need it to work with a very short string like that... I can get it to work... just the code that does the subsampling needs to change... not the rest. )

[fabrice]

earonesty,

I have tried your new release. The output is fellow.

Threshold used: 1 out of 4. What is the mean here?

Why both Adapter date_day and Adapter adaptor1 counted 4 at the 'start' of 'test.fastq'?

I also noticed that the "N" at the begin/end of reads that are not trimed.

_____________________output___________________________

fastq-mcf test.fa test.fastq -q 0 -l 1 -o test.trim
Scale used: 2.2
Threshold used: 1 out of 4
Adapter date_day (CCTTAAGG): counted 4 at the 'start' of 'test.fastq', clip set to 1
Adapter adaptor1 (AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT): counted 4 at the 'start' of 'test.fastq', clip set to 1
Files: 1
Total reads: 4
Too short after clip: 0
Clipped 'start' reads: Count: 3, Mean: 8.00, Sd: 0.00

[earonesty]

They do deserve some elaboration in the help:

# -s N.N Log scale for clip pct to threshold (2.5)

Scaling factor which causes more frequent occurrences to clip heavily. The negative log base "scale" of the ratio of adapters found becomes the minimum match-length for the adapter.

# -t N % occurance threshold before clipping (0.25)

Minimum number of times an adapter needs to be found before clipping is considered necessary


# -m N Minimum clip length, overrides scaled auto (1)

Replaces scaling algorithm (above) with hardcoded limit

# -p N Maximum adapter difference percentage (20) This is for mismatch?

Yes this is the max # mismatch

# -k N sKew percentage causing trimming (2)

If one of the bases is only 2% of the reads, then that cycle is considered "skewed".... if it's at the edge of a read, it's trimmmed off.

[fabrice]

In my opions, most of adaptor is at the start of reads. Is it right?

In my RNA-seq data, I just see a pattern at the start, not the end.

In fact, I cannot understand the means of these options.

-s N.N Log scale for clip pct to threshold (2.5)
-t N % occurance threshold before clipping (0.25)
-m N Minimum clip length, overrides scaled auto (1)
-p N Maximum adapter difference percentage (20) This is for mismatch?
-k N sKew percentage causing trimming (2)

[earonesty]

You're right, I'd introduced a bug in detecting starts. I'm sending a release and adding it to my tests. The problem was it should have been saying at 'start' of sequence... not 'end'. By thinking it was an 'end' adapter, it was trimming from the right... hence removing the whole sequence.

(REVISION: The problem was very short adapter sequences in the input file. fastq-mcf can't be used that way. It uses subsampling to get the adapter thresholds, and requires a minimum of 15 character matches during the subsampling pass)

[fabrice]

But when I run it on the test data. It seems fastq-mcf remove the whole reads. Do I set something wrong?

fastq-mcf adapters.fa test.fastq -q 0 -l 1 -o test.trim
Scale used: 2.2
Threshold used: 1 out of 4
Adapter date_day (CCTTAAGG): counted 3 at the 'end' of 'test.fastq', clip set to 1
Files: 1
Total reads: 4
Too short after clip: 3

----------------adapters.fa----------------
>date_day
CCTTAAGG
>adaptor1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT

------------------test.fastq------------------------
@test1
CCTTAAGGAAAAAAAAAAGGGGGGGGGG
+test1
HHHHHHHHHHHHHHHHHHHHHHHHHHHH
@test2
CCTTAAGGAAAAAAAAAGGGGGGGGGGG
+test2
HHHHHHHHHHHHHHHHHHHHHHHHHHHH
@test3
AGAGAGAGAGAGAGAGAGAGAGAGAGAG
+test3
HHHHHHHHHHHHHHHHHHHHHHHHHHHH
@test4
CCTTAAGGTTGACGTGATCGACACCTGG
+test4
[[[[[[[[[[[[[[[[[[[[[[[[[[[[

[earonesty]

It only removes the whole read if the remaining length is less than the minimum length threshold, which, by default is 15:

-l N Minimum remaining sequence length (15)

15's pretty low, but you can lower it.

[fabrice]

I have tryed your ea-utils. But it seems as the same FASTXtoolkit adapter trimming. ea-utils also remove the whole read which contained adapter.

We rewrote a lot of fastx's toolkit stuff, and posted it here: https://code.google.com/p/ea-utils/. It attempts to do things like adapter removal, trimming, etc... without as much configuration by detecting presence of adapters located in a common file.

[earonesty]

I think the barcode in the middle was making it odd. Also, I think your solution is great.

[angelawu]

Oh, the adapter sequences are not identical. If you look closely at the middle portion of the sequences, there is a barcode in the middle that is different for each sequence. But I also do not think this would make any difference...

In any case, I think I have a solution to my particular application, so I don't know how much time I want to spend debugging this, but thanks for reminding me of the -d option, which will surely come in handy later on as well. The -e option may be the trick, since the barcode only begins in the middle of the adapter sequence?

Thanks once again!

[earonesty]

Hi earonesty,

Thanks for getting back to me!
Here is an example of the output I received:

Scale used: 2.2
Threshold used: 101 out of 40000
Adapter TruSeq-Adapter1 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter2 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter3 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTAGGCATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter4 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTGACCAATCTCGTATGCCGTCTTCTGCTTG): counted 10159 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter5 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACACAGTGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter6 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGCCAATATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter7 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCAGATCATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter8 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACACTTGAATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter9 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGATCAGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter10 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTAGCTTATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter11 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGGCTACATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter12 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCTTGTAATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Files: 1
Total reads: 21964185
Too short after clip: 35672
Clipped 'start' reads: Count: 13283064, Mean: 1.58, Sd: 1.15
Trimmed 394023 reads by an average of 7.15 bases on quality < 10


So far, what I understand is that my samples probably have a lot of adapter-pair ligations in there without any genomic insert. This leads to the entirety of my 36bp read being a portion of the index/adapter. And since the adapter is much longer than 36bp, i think those reads are not being removed. e.g.:

My sample DNA: ADAPTER1adapter2
read: dapte

I say this because when I put the cleaned up reads through FastQC again, I see that all the "Over-represented Sequences" that are TruSeq Indexes are still present in my file.

I've managed to resolve this particular issue by basically copying in the sequence given by FastQC as the overrepresented sequence, and using those in a fasta file as the adapter sequences. It works well for my case, so maybe there is nothing wrong with the toolkit, and it's just my particular sample?

Yes, I understand that Illumina reads use Phred64, but I always convert directly to Sanger Phred33 as soon as I get my files, which is why I put the -P 33 option in there.

Thanks,

Angela

- Your adapter file seems to have the same sequence over and over? I'm not sure how that will affect things. TruSeq-Adapter2 is the same as TruSeq-Adapter1.... etc. Try just using 1 per unique sequence. This probably won't help.

- Out of 40000 reads, 10000 had an exact match for 15 base pairs of adapter sequence. That's a lot. So when it says "clip set to 1" it will clip any matching subsequence.

- It only discarded 35672 reads and only a few bases. That's surprising to me considering the number of sequences it found in the subsample with exact matches. I would expect a higher rate of discards, and a higher number of mean bases clipped.

- This is a situation where I wish I could see about 100K reads from your sample and just run it a few times to see what happened why it did that. It should be walking the adapter along the sequence looking for the best match. It seems to be stopping early on....or perhaps the sequences that match the adapter are somewhere else (at the end...?) and it guessed wrong (you can force -e)

- There's also an undocumented "-d" option that spits out lots of debug info that I find useful.

[angelawu]

1. It does remove partial matches. It searches only from one "end" of the file. The default settings are very conservative, so if it's removing 10 million reads, that's an enormous number - you may want to change the settings to be more aggressive for that data.

2. Can you post the summary output... it should say why sequences were removed/clipped, and why etc.

3. Until very recently, GAII's output base-64 by default, not 33, so you may want to double-check that.

EXAMPLE OUTPUT:

Scale used: 2.2
Threshold used: 101 out of 40000
Adapter ILMN RT_primer_rc (TCGTATGCCGTCTTCTGCTTG): counted 193 at the 'end' of 'example.fastq', clip set to 6
Adapter FLUIDIGM Index-SP (AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCG): counted 1063 at the 'end' of 'example.fastq', clip set to 4
Files: 1
Total reads: 250000
Too short after clip: 53
Clipped 'end' reads: Count: 16612, Mean: 18.12, Sd: 17.44
Trimmed 24474 reads by an average of 10.81 bases on quality < 10

Hi earonesty,

Thanks for getting back to me!
Here is an example of the output I received:

Scale used: 2.2
Threshold used: 101 out of 40000
Adapter TruSeq-Adapter1 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter2 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter3 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTAGGCATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter4 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTGACCAATCTCGTATGCCGTCTTCTGCTTG): counted 10159 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter5 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACACAGTGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter6 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGCCAATATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter7 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCAGATCATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter8 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACACTTGAATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter9 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGATCAGATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter10 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACTAGCTTATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter11 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACGGCTACATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Adapter TruSeq-Adapter12 (GATCGGAAGAGCACACGTCTGAACTCCAGTCACCTTGTAATCTCGTATGCCGTCTTCTGCTTG): counted 10158 at the 'start' of './sanger-fastq/s_2.fastq', clip set to 1
Files: 1
Total reads: 21964185
Too short after clip: 35672
Clipped 'start' reads: Count: 13283064, Mean: 1.58, Sd: 1.15
Trimmed 394023 reads by an average of 7.15 bases on quality < 10


So far, what I understand is that my samples probably have a lot of adapter-pair ligations in there without any genomic insert. This leads to the entirety of my 36bp read being a portion of the index/adapter. And since the adapter is much longer than 36bp, i think those reads are not being removed. e.g.:

My sample DNA: ADAPTER1adapter2
read: dapte

I say this because when I put the cleaned up reads through FastQC again, I see that all the "Over-represented Sequences" that are TruSeq Indexes are still present in my file.

I've managed to resolve this particular issue by basically copying in the sequence given by FastQC as the overrepresented sequence, and using those in a fasta file as the adapter sequences. It works well for my case, so maybe there is nothing wrong with the toolkit, and it's just my particular sample?

Yes, I understand that Illumina reads use Phred64, but I always convert directly to Sanger Phred33 as soon as I get my files, which is why I put the -P 33 option in there.

Thanks,

Angela

[earonesty]

1. It does remove partial matches. It searches only from one "end" of the file. The default settings are very conservative, so if it's removing 10 million reads, that's an enormous number - you may want to change the settings to be more aggressive for that data.

2. Can you post the summary output... it should say why sequences were removed/clipped, and why etc.

3. Until very recently, GAII's output base-64 by default, not 33, so you may want to double-check that.

EXAMPLE OUTPUT:

Scale used: 2.2
Threshold used: 101 out of 40000
Adapter ILMN RT_primer_rc (TCGTATGCCGTCTTCTGCTTG): counted 193 at the 'end' of 'example.fastq', clip set to 6
Adapter FLUIDIGM Index-SP (AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCG): counted 1063 at the 'end' of 'example.fastq', clip set to 4
Files: 1
Total reads: 250000
Too short after clip: 53
Clipped 'end' reads: Count: 16612, Mean: 18.12, Sd: 17.44
Trimmed 24474 reads by an average of 10.81 bases on quality < 10

[angelawu]

Question about fastq-mcf

Hi,

I encountered an issue when using fastq-mcf on my GA2 generated 1x36 reads, and wondering if you could shed some light.

So I made my fasta file with all the TruSeq adapter sequences in there, and ran fastq-mcf using that file, -P Phred scale set to 33 for my files are in Sanger fastq format. All other parameters were left as default.

After trimming was completed, the outfile reports removing about 10 million reads out of 24 million.

I run the trimmed file through FastQC, and under the "over-represented sequences" tab, I see that partial adapter sequences (e.g. starting from bp #2) are still over-represented in my file, which suggests that they were not trimmed.

My question is, does fastq-mcf remove partial matches to adapter sequences provided, as well as full? If so, am I doing something wrong with the way I am using the tool?

I am pretty new to bioinformatics, so sorry if this is a stupid question...

Thank you!

[earonesty]

Note: I made a change recently to properly use the "better quality base" in the overlapping region... there was a bug in it that someone pointed out. If you're using it, you'll want the newer version.