[GenoMax]

Because your literal adapter is also matching the beginning of the read entire sequence to the right is removed (ktrim=r). If you change the sequence at the beginning by one base you will see that the initial part is retained as you expect.

[Stimpsky]

Ah OK, that explains it, nevermind. I thought ktrim=r would go from right to left and stop at the first hit. Thanks for the quick reply!

[cuencam]

Hi Brian,
I would like to switch from solexaqa to BBduk for my read trim and filtering option, however since we are working with bacterial strains typification we would like to have the option to only keep trimmed reads where no individual base has a quality lower than a defined threshold (instead of average region quality). Could this be done with BBduk?

Thanks!

[Brian Bushnell]

Hi cuencam,

It is currently not possible to do this, other than discarding all reads that have any undefined (quality 0) bases with "maxns=0". I never saw a reason to discard all reads with a single base below a specified cutoff. It would be simple enough to add (or implement via a custom script), but can you explain why you're doing it? The "average region quality" method is the best at maximizing coverage while minimizing the total number of errors.

Edit - anyway, it was quick to add, so it will be in the next release as the "mbq" ("minbasequality") flag.

[cuencam]

Hi Brian,
Thanks for such a quick response, and implementing this so fast!
Our interest is that during SNV calling in low coverage regions or low abundance taxa (in metagenomics) base quality can be more important than coverage. This way we can assess properly different alleles and avoid the creation of artifacts
Cheers!

Edit - Do you have an estimated next release date?

[cuencam]

Hi Brian,
In the same lines of my previous question, what is the rationale of using maq=10? We are interesting in de novo assembly of metagenomic data and we were worried that low quality bases at the ends of the reads might feed artificial k-mers in to the assembler (SPADES). I read that you recommend read normalization, but since our coverage is highly unequal (due to unequal species abundance, not because sequencing artifacts) we are worried that this might introduce more biases than the ones it solves.

We were thinking on using your newly implemented option "mbq" to secure that all bases have 20 as minimum quality. Do you believe that this is a good alternative?

[Brian Bushnell]

"maq=10" is to throw away really junky reads. The only way to really verify whether a setting is beneficial is to actually test it, unfortunately. But personally, I think "mbq=20" would be too aggressive (particularly if your sequencing run had a single low-quality cycle, in which case it would discard all of the data)... if you really want to get rid of the low-quality trailing bases, I'd suggest quality-trimming instead (qtrim=r trimq=14 or something like that). Spades is pretty robust with respect to low-quality data anyway; the biggest problem is that it low quality reads balloon the kmer-space which can make it run out of memory.

The main advantage of normalization with metagenomes, in fact, is that it removes a lot of data which allows Spades to run on datasets that it can't otherwise handle. It's not strictly beneficial and if you can assemble a metagenome without normalization, that may be better - sometimes normalization improves the assembly, sometimes it doesn't.

[cuencam]

Thanks for this response! I'm pretty sure that your excellent user support is only comparable to the high quality of your tools!

I will implement quality-trimming at a higher threshold and then test. I do agree that mbq=20 is hard for assembly (but probably useful for SNV).
Cheers

[EssigSchurke]

Hi Brian,

I tried to filter reads longer 10bp. I used the following command:

bbduk.sh -in=input.fq -out=output.fq -maxlength=10

However, nothing happens, I get the same amount of reads as in the input. But all reads are longer 10bp.
I used the latest version of bbduk 37.53

Test Input:

@test
ACTGGACTTGGAGTCAGAAGGC
+
b\\[\ZZ[][a]_]]cbbbabc

Input:                  	1 reads 		22 bases.
Total Removed:          	0 reads (0.00%) 	0 bases (0.00%)
Result:                 	1 reads (100.00%) 	22 bases (100.00%)

[jazz710]

The BBDuk commands don't have '-' before them. Your command should read:

bbduk.sh in=input.fq out=output.fq maxlength=10

Give that a shot?

[EssigSchurke]

With or without "-" does not matter, I get same results.

[cuencam]

Hi EssigSchurke
The flag is minlength=10

The whole command is
bbduk.sh in=input.fq out=output.fq minlength=10

Edit:

I misread your question. The command provided by jazz710 is the appropriate, and works on my computer. You want to remove the big reads, correct?

[EssigSchurke]

Hi cuencam,

minlength=10 filters only reads shorter 10bp. I want to filter reads longer 10bp, whereas 10bp is only a dummy for my test case.

[EssigSchurke]

Yes, I want to exclude large reads, but I tested the command provided by jazz710. It produces the same result, the test read is still in the output.

[Brian Bushnell]

Actually, all the BBTools strip off the leading "-" so you can put as many of them as you want

This is a bug. Thanks for the report! It looks like BBDuk only removes reads under minlen or over maxlen if they were trimmed; untrimmed sequences will pass regardless of their length. Sorry about that! Reformat actually works correctly in this case, though:

reformat.sh in=x.fq out=y.fq minlen=A maxlen=B

I'll fix BBDuk ASAP. Thanks again!