Yes .. BBMap can do that!

[GenoMax]

I think @duane.hassane is asking to remove primers from amplicons with a common core sequence but with variable (length and/or sequence?) primers on the ends of that core sequence for each group.

[duane.hassane]

Thanks. To be more clear, I am after removing the primers themselves (either from the fastq or BAM) and preserving the intervening sequence only. In this use case, there are hundreds of different primer pairs with known mapping locations (because they are used for targeted enrichment of different exons). Thus, for each primer pair, the intervening sequence (target being enriched) is different.

[luc]

Have you tried to remove these primer sequences with bbduk?

[duane.hassane]

Yes, bbduk with restrictright and/or restrictleft set to constrain the search to the ends where primers are located. This is actually extremely fast and efficient even with hundreds of primers. Unfortunately, we have encountered some specificity issues with this approach (over- and under-trimming) leading to small gaps in coverage ... with no one set of parameters being optimal for all amplicons.

[Brian Bushnell]

Hmm, I don't currently have a better solution than BBDuk or the msa/cutprimers pair. You might get a better result if you first bin the amplicons by which primers they use, though, so that then you can trim them only using the correct primers using aggressive settings, rather than with all possible primers. You can try using Seal for that. For example:

seal.sh in=amplicons.fq ref=primer1.fa pattern=out%.fq k=(whatever)

That will produce one file per sequence in the primer file, with all the reads that best match it (meaning, share the most kmers). For that purpose you should only use the left primers or right primers as reference, not both at the same time. Or, ideally, the reference would look like this:

Primer pairs: {(AAAT, GGCC), (TTTT, GGGG)}

Reference file:

>1
AAATNGGCC
>2
TTTTNGGGG

...etc, where each primer pair is concatenated together with a single N between the sequences.

[duane.hassane]

Ok. Will take a look at this approach. Thank you!

[fdts]

Just a comment, not sure if this is the place:
It would be helpful for downstream analysis if the refstat files that can be output using bbsplit produced 0s (zeros) for references that have no reads mapped to them instead of omitting data for these references. I was working on generating a table for visualization when I noticed missing data.

[Brian Bushnell]

Zero-count lines are suppressed by default, but they should be printed if you include the flag "nzo=f" (nonzeroonly=false). That's not documented in BBSplit's shellscript, just BBMap's; I'll add it to the documentation.

[fdts]

Wonderful, thanks!

[Thowell]

I have been using bbmap (and other bbtools) for some of my analyses, and so far they are a great set of tools! Unfortunately I have come across a problem, and some preliminary searching hasn't turned up an answer.

I am trying to build bbmap into a galaxy workflow, using it to remove reads that map to a reference of a contaminant (E.coli in this case). The problem is that I just want the unmapped reads in fastq format, but because galaxy names everything with the .dat extension, I cannot tell bbmap to output in fastq format and it defaults to sam. Is there another way to specify the output format?

I'm sure I could convert the sam back to fastq, but the picard tools sam to fastq tool is throwing an error saying that there are unpaired reads in the sam. I know I can probably do this another way, but it seems silly to use another tool when bbmap can already output in the format I want!

I'm also confused as to why there would be unpaired reads, since the bbmap documentation for "outu" says:

Write only unmapped reads to this file. Does not include unmapped paired reads with a mapped mate.

Sounds to me like all the reads should be paired?


Thank you in advance, and thanks for the great tools!

[Brian Bushnell]

Thanks for noting this. Some tools (like reformat) support the "extin" and "extout" flags which let you override the default, so you could do this:

reformat.sh in=file.dat out=file2.dat extin=.sam extout=.fq

But, BBMap doesn't support that right now. I'll add it. And I don't particularly recommend sam -> fastq conversion because the names change, since in sam format read 1 and read 2 must have identical names, whereas in fastq format they will typically have "/1" and "/2" or similar to differentiate them. Though you can do that conversion if you want.

I have not used Galaxy and don't know what's possible, but until I make this change, I would suggest one of these:

1) Use BBDuk for this filtering; its default output format is fastq and it's probably faster than BBMap anyway in this case. The syntax is very similar. On the command line, it would be something like "bbduk.sh in=reads.fq outu=clean.fq ref=ecoli.fasta".

2) Tell BBMap "outu=stdout.fq" and pipe that to a file, if Galaxy supports pipes.

As for your question about pairing, the normal behavior in paired-mapping mode is:

"out=" will get everything.
"outm=" will get all pairs in which either of the reads mapped to the reference.
"outu=" will get all pairs in which neither read mapped to the reference.

For BBDuk, it's slightly different but essentially the same:
"out=" is the same as "outu=".
"outu", aka "out", will get all pairs in which neither had a kmer match to the reference.
"outm" will get all pairs in which either had a kmer match to the reference.
For BBDuk, this behavior can be changed with the "reib" (removeIfEitherBad) flag. The assumption of that flag's name is that the reference is contaminants being filtered against, so the default "reib=true" means any pair where either matches the contaminant is removed.

So, for both tools, if the input data is paired, the output data will also be paired - pairs are always kept together in all streams.

[Thowell]

Thank you for the quick reply Brian.

I was able to get things working with a pipe. I'm guessing the reads have to be interleaved with this method, but that will work fine for me until you can implement the alternate output flag.

Thanks again!

[mendezg]

How would you like us to cite bbduk in papers?

[GenoMax]

Originally Posted by mendezg

How would you like us to cite bbduk in papers?

In the past Brian has asked people to link to the project site on source forge.

[Brian Bushnell]

Hi - sorry I somehow missed this question! Yes, as Genomax stated, please just cite it something like this (altered according the format of the journal):

"BBDuk - Bushnell B. - sourceforge.net/projects/bbmap/"

[vmikk]

Hello!
Is it possible to use cutprimers.sh to cut the sequence AND to preserve the primer sites around?

[Brian Bushnell]

Not currently... but I'll plan to add a flag for that.

[Brian Bushnell]

I added the "include" flag to cutprimers. Default is "include=f". If you set "include=t" the primers will be retained for the output.

[vmikk]

Hello Brian! Thanks a lot for the implementation of this feature!
Meanwhile I thought to modify sam files from msa.sh, but the out of the box functionality is much more convenient!
Thanks again!

[dkainer]

Brian,

is there a way with the BB Suite to demultiplex paired-end reads based on inline barcodes, like Flexbar does?

I can see it can be done one barcode at a time by outputting matching reads based on the first 6 left bases. But can it be done in one command to demultiplex for multiple barcodes?

cheers
DK