SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bowtie, an ultrafast, memory-efficient, open source short read aligner Ben Langmead Bioinformatics 513 05-14-2015 03:29 PM
Introducing BBMap, a new short-read aligner for DNA and RNA Brian Bushnell Bioinformatics 24 07-07-2014 10:37 AM
Miso's open source joyce kang Bioinformatics 1 01-25-2012 07:25 AM
Targeted resequencing - open source stanford_genome_tech Genomic Resequencing 3 09-27-2011 04:27 PM
EKOPath 4 going open source dnusol Bioinformatics 0 06-15-2011 02:10 AM

Reply
 
Thread Tools
Old 04-22-2017, 12:55 PM   #461
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 37
Default

Hi Brian,

I'm wondering if there is any way you can get bbmap (or another of your tools) to give a consensus sequence of an alignment? I went to map to a ref sequence and then use that to create a consensus sequence and then use that to map again and then call variants.

Thanks!
jweger1988 is offline   Reply With Quote
Old 04-24-2017, 11:08 AM   #462
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Hi Jweger,

Sorry, I don't have anything that will do that. Clumpify allows something sort of similar; you can create consensus sequence from raw reads, and map those. But that loses per-base depth information which is important for variant-calling, so I don't think it's what you want.
Brian Bushnell is offline   Reply With Quote
Old 04-26-2017, 10:43 PM   #463
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

I was testing something with BBMap today and I realized I had forgotten something about how to use it and couldn't figure it out. I was mapping RNA-Seq and I thought that it would report spliced alignment cigars (using sam 1.3) as /[0-9]+M[0-9]+N[0-9]+M/ but it was reporting the introns as deletions with the [D] cigar value. Is that right or is there a way to get it to not do that?

Thanks-
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 04-27-2017, 01:32 AM   #464
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

You can use the flag "intronlen" to control this. For example, "intronlen=100" will change the reporting so that all deletions of at least 100bp will reported using 'N' instead of 'D' in the cigar string.
Brian Bushnell is offline   Reply With Quote
Old 04-27-2017, 10:08 AM   #465
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

Oh yes! Thanks Brian.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 05-02-2017, 04:52 AM   #466
jweger1988
Member
 
Location: Paris, France

Join Date: Apr 2017
Posts: 37
Default

Hi Brian,

Is there anyway to use bbmap (or any other of your tools) to map a read to a reference file and then trim anything to the left of the reference sequence?

For example

My reference is
XXXXXXXXXXXXXXXXXX

And my reads are

NNNNXXXXXXXXXXXXXXXXXX
NNNXXXXXXXXXXXXXXXXXX
NNXXXXXXXXXXXXXXXXXX

I want them to be
XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX

I basically want to just trim anything of the left of the reads that doesn't match my reference? Thanks in advance.
jweger1988 is offline   Reply With Quote
Old 05-02-2017, 11:02 AM   #467
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

To do that, you'd need to trim all soft-clipped bases. I don't have any programs that will do so, but it looks like I could add it to Reformat without too much difficulty.
Brian Bushnell is offline   Reply With Quote
Old 05-02-2017, 12:31 PM   #468
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Hi Brian,
Is there a way to set a minimum input fragment length when mapping with BBMap?
darthsequencer is offline   Reply With Quote
Old 05-02-2017, 12:34 PM   #469
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,015
Default

You could pre-filter input data with
Code:
reformat.sh in=file.fq out=filt.fq minlength=N
GenoMax is offline   Reply With Quote
Old 05-02-2017, 12:51 PM   #470
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Actually, while it's not documented, the flag "minlength" also works with BBMap. Reads shorter than that will be discarded completely (they won't be output as unmapped).
Brian Bushnell is offline   Reply With Quote
Old 05-02-2017, 01:00 PM   #471
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Quote:
Originally Posted by Brian Bushnell View Post
Actually, while it's not documented, the flag "minlength" also works with BBMap. Reads shorter than that will be discarded completely (they won't be output as unmapped).
Thanks for the tip - save me some space so I don't have to make additional files. I deal with a lot of libraries with different read lengths.

Is there a way to make BBMap require the smallest read to be some fraction of the longest read length? I know that's a niche use but BBMap always suprises me with it's built in functions.
darthsequencer is offline   Reply With Quote
Old 05-02-2017, 01:23 PM   #472
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Ahh... sorry, but nope That would require reading the file twice, so it would not be easy to implement.
Brian Bushnell is offline   Reply With Quote
Old 05-02-2017, 01:25 PM   #473
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Quote:
Originally Posted by Brian Bushnell View Post
Ahh... sorry, but nope That would require reading the file twice, so it would not be easy to implement.
Oh yeah - duh
darthsequencer is offline   Reply With Quote
Old 05-02-2017, 01:30 PM   #474
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Quote:
Originally Posted by Shini Sunagawa View Post
Dear Brian,

I have been looking for a tool that would quickly dereplicate (100% containments) nucleotide sequences and track for each unique sequence the identifiers of the removed duplicates.

Something like:

dedupe.sh in=in.fa out=out.fa outd=outd.fa mid=100 mop=100

where:

in.fa:
seq1
seq2 (contained in seq1)
seq3 (contained in seq1)
seq4

out.fa:
seq1
seq4

outd.fa:
seq2
seq3

I am interested in:
seq1<tab>seq2,seq3
seq4

dedupe.sh does a fantastic job in returning out and outd, but I cannot find any option that would return the information I am interested in. Is this something that I am missing? Otherwise, I believe this could be a great feature, since compared to other tools that return this information, dedupe is so much faster.

Best,
Shini
Did anything like this get added to BBMap? Would be really helpful for me too!
darthsequencer is offline   Reply With Quote
Old 05-02-2017, 01:53 PM   #475
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Hmmm... no, not yet, though I did add into Clumpify's dereplication step the ability to count the number of duplicate reads and add "count=3", for example, to the name of a read representing 3 total reads (itself and 2 duplicates). It would not be difficult to modify that to report read identifiers. I'll add it to my todo list.
Brian Bushnell is offline   Reply With Quote
Old 05-02-2017, 06:13 PM   #476
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Quote:
Originally Posted by Brian Bushnell View Post
Hmmm... no, not yet, though I did add into Clumpify's dereplication step the ability to count the number of duplicate reads and add "count=3", for example, to the name of a read representing 3 total reads (itself and 2 duplicates). It would not be difficult to modify that to report read identifiers. I'll add it to my todo list.
Well my real use it to group sequences that >99% similar. Whereas clumpify find exact matches? Although I suppose that clumpify could be run in error correct mode with "midid" and it should be the same as dedupe with minidentity?
darthsequencer is offline   Reply With Quote
Old 05-02-2017, 06:56 PM   #477
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Clumpify can consider sequences as duplicates if they have at most X substitutions, but it's not as flexible as Dedupe. For example, Clumpify requires duplicates to overlap 100% with neither overhanging, while Dedupe allows containments (this only matters when using variable-length sequences) and also allows indels. What I actually added to my todo list was to update both of them with that capability, since it seems useful.
Brian Bushnell is offline   Reply With Quote
Old 05-02-2017, 08:14 PM   #478
darthsequencer
Member
 
Location: california

Join Date: Feb 2012
Posts: 35
Default

Quote:
Originally Posted by Brian Bushnell View Post
Clumpify can consider sequences as duplicates if they have at most X substitutions, but it's not as flexible as Dedupe. For example, Clumpify requires duplicates to overlap 100% with neither overhanging, while Dedupe allows containments (this only matters when using variable-length sequences) and also allows indels. What I actually added to my todo list was to update both of them with that capability, since it seems useful.
Great! Yeah I'd like to know the duplicate membership for containments too.
darthsequencer is offline   Reply With Quote
Old 05-03-2017, 05:25 AM   #479
sebastian_1
Junior Member
 
Location: Czech republic

Join Date: Mar 2017
Posts: 3
Default

Hello Brian,

I would like to ask you a suggestion for bbmap.
I am trying to reassemble a bin from a metagenomic data, hoping that I will get better assembly if I use just the mapped reads.
I tried bbmap.sh on normal parameters and used the outm to take all the aligned reads, then I normalized the reads with bbnorm.sh, and reassembled with SPAdes. I want to note that the initial metagenomic assembly was done on normalized reads, and SPAdes does error correction, but I did not use these libraries, I used the adapter and quality trimmed libraries but not normalized and error corrected (this is why I do normalization after mapping).

I got better assembly, (some longer scaffolds, and slightly larger N50) but checking briefly the SSU I noticed that some "contaminants" were present. Also, the amount of SSU sequences was much higher than expected. (I expect 4, 3 complete ones and one near complete).
In the metagemic data (assembled using all the data) I have 10 SSUs, but here I got a lot of them (15+) and most of them are really partial.

What I am thinking is that bbmap, includes in the output some(not all) reads from other bacterial SSU which map to a certain degree to my reference (since it can have very conserved regions) then SPAdes is somehow confused by these, and fragments my SSU sequences in multiple places due to these reads. Sorry if this sounds strange, but I am just speculating, I am not sure if this is the case, and unfortunately I am not an expert bioinformatician.

I was thinking to add to the command line the parameters minidentity=0.98 idfilter=0.98 hoping that I will somehow avoid the mapping of "non-specific" reads.

I avoid using the perfectmode=t due to the fact that SPAdes does error correction, and this would somehow cause some small mismatches, and thus I would lose some reads.

Would you have any suggestion for better parameters of bbmap?

Thank you!

PS: I am using PE 2X250bp and PE2x300bp libraries for the mapping.
sebastian_1 is offline   Reply With Quote
Old 05-03-2017, 09:23 AM   #480
sk8bro
Member
 
Location: boston

Join Date: Feb 2012
Posts: 30
Default Non-Deterministic BBmap results... how ensure deterministic?

Hi Brian,

It is important that my pipeline's results can be perfectly reproduced. I notice that the non-deterministic behavior is coming from the human read removal from my datasets... here are the parameters I have specified in this call:

bbmap.sh\
-Xmx23g\
minid=0.9\
idfilter=0.9\
maxindel=3\
bwr=0.16\
bw=12\
minhits=2\
printunmappedcount=t\

The numbers are close (ie. 672510 vs 672492)

I run the PE reads through BBduk to remove low-quality pairs first. It looks like that output is sorted the same and deterministic upon re-runs.

Thanks for your thoughts, Kate

Last edited by sk8bro; 05-03-2017 at 09:32 AM.
sk8bro is offline   Reply With Quote
Reply

Tags
bbmap, metagenomics, rna-seq aligners, short read alignment

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 03:33 PM.


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