SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BBMap (aligner for DNA/RNAseq) is now open-source and available for download. Brian Bushnell Bioinformatics 537 Yesterday 05:01 AM
BBMap for BitSeq dietmar13 Bioinformatics 1 04-30-2015 08:40 AM
Please help my BBMap cannot remove Illumina adapter TofuKaj Bioinformatics 4 04-28-2015 08:53 AM
BBMap Error Phage Hunter Bioinformatics 5 01-14-2015 04:34 AM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 09:37 AM

Reply
 
Thread Tools
Old 06-11-2015, 09:24 AM   #21
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,409
Default

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.
GenoMax is offline   Reply With Quote
Old 06-11-2015, 09:44 AM   #22
duane.hassane
Junior Member
 
Location: New York, NY

Join Date: Jun 2012
Posts: 5
Default

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.
duane.hassane is offline   Reply With Quote
Old 06-11-2015, 09:58 AM   #23
luc
Senior Member
 
Location: US

Join Date: Dec 2010
Posts: 282
Default

Have you tried to remove these primer sequences with bbduk?
luc is offline   Reply With Quote
Old 06-11-2015, 10:18 AM   #24
duane.hassane
Junior Member
 
Location: New York, NY

Join Date: Jun 2012
Posts: 5
Default

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.
duane.hassane is offline   Reply With Quote
Old 06-11-2015, 10:41 AM   #25
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

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.

Last edited by Brian Bushnell; 06-11-2015 at 11:17 AM.
Brian Bushnell is offline   Reply With Quote
Old 06-11-2015, 10:47 AM   #26
duane.hassane
Junior Member
 
Location: New York, NY

Join Date: Jun 2012
Posts: 5
Default

Ok. Will take a look at this approach. Thank you!
duane.hassane is offline   Reply With Quote
Old 07-01-2015, 05:50 AM   #27
fdts
Junior Member
 
Location: DE

Join Date: May 2015
Posts: 4
Default

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.
fdts is offline   Reply With Quote
Old 07-01-2015, 09:38 AM   #28
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

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.
Brian Bushnell is offline   Reply With Quote
Old 07-01-2015, 11:39 PM   #29
fdts
Junior Member
 
Location: DE

Join Date: May 2015
Posts: 4
Default

Wonderful, thanks!
fdts is offline   Reply With Quote
Old 07-11-2015, 06:35 PM   #30
Thowell
Junior Member
 
Location: California

Join Date: Nov 2011
Posts: 3
Default

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:
Quote:
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!
Thowell is offline   Reply With Quote
Old 07-11-2015, 07:26 PM   #31
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

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.
Brian Bushnell is offline   Reply With Quote
Old 07-11-2015, 08:32 PM   #32
Thowell
Junior Member
 
Location: California

Join Date: Nov 2011
Posts: 3
Default

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!
Thowell is offline   Reply With Quote
Old 09-16-2015, 08:32 AM   #33
mendezg
Junior Member
 
Location: Maryland

Join Date: Oct 2014
Posts: 2
Default

How would you like us to cite bbduk in papers?
mendezg is offline   Reply With Quote
Old 09-18-2015, 04:23 PM   #34
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,409
Default

Quote:
Originally Posted by mendezg View Post
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.
GenoMax is offline   Reply With Quote
Old 09-20-2015, 07:18 PM   #35
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

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/"
Brian Bushnell is offline   Reply With Quote
Old 09-24-2015, 04:42 AM   #36
vmikk
Junior Member
 
Location: RF

Join Date: Jun 2015
Posts: 3
Default

Hello!
Is it possible to use cutprimers.sh to cut the sequence AND to preserve the primer sites around?
vmikk is offline   Reply With Quote
Old 09-24-2015, 08:57 AM   #37
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

Not currently... but I'll plan to add a flag for that.
Brian Bushnell is offline   Reply With Quote
Old 09-28-2015, 09:56 AM   #38
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,637
Default

I added the "include" flag to cutprimers. Default is "include=f". If you set "include=t" the primers will be retained for the output.
Brian Bushnell is offline   Reply With Quote
Old 09-28-2015, 10:37 PM   #39
vmikk
Junior Member
 
Location: RF

Join Date: Jun 2015
Posts: 3
Default

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!
vmikk is offline   Reply With Quote
Old 09-30-2015, 06:36 PM   #40
dkainer
Junior Member
 
Location: Australia

Join Date: May 2015
Posts: 9
Default

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
dkainer is offline   Reply With Quote
Reply

Tags
bbmap

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:27 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO