[sk8bro]

Hi Brian,

Thanks for letting me know. That would be really awesome... keep me posted! My email is [email protected]. Til then, I will probably need to stick with the single thread.

[Brian Bushnell]

BBSplit is based on BBMap, so it is also nondeterministic in paired mode with multiple threads. BBDuk and Seal (which can be used similarly to BBSplit) are always deterministic. I could add a deterministic mode to BBMap without too much effort.

[sk8bro]

Thank you for the explanation and suggested resolution Brian,

How about BBduk and BBsplit? Can they similarly behave non-deterministically in paired-end mode?

[Brian Bushnell]

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.

It sounds like related strains are being binned together. Ribosomes are always hard to assemble anyway. But yes, you can try minidentity=0.98 idfilter=0.98. That might knock out a bunch of your 2x300bp reads since those are usually low quality, but you can add the flags "qtrim=r trimq=15 untrim" to trim the reads before mapping, then undo the trimming after mapping, so that the low-quality right end won't contribute to the identity calculation. Alternatively, you can try error-correction via Tadpole and/or BBMerge prior to mapping:

bbmerge.sh in=reads.fq out=ecco.fq ecco mix rsem strict k=62
tadpole.sh in=ecco.fq out=ecct.fq ecc k=62

[Brian Bushnell]

Pre-filtering with reformat.sh in the mean time then?

Yep, looks like it

[Brian Bushnell]

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

Hi Kate,

BBMap is always nondeterministic when run in paired-end mode with multiple threads, because the insert-size average is calculated on a per-thread basis, which affects mapping; and which reads are assigned to which thread is nondeterministic. The only way to avoid that would be to restrict it to a single thread (threads=1), or map the reads as single-ended and then fix pairing afterward:

bbmap.sh in=reads.fq outu=unmapped.fq int=f
repair.sh in=unmapped.fq out=paired.fq fint outs=singletons.fq

In this case you'd want to only keep the paired output.

[GenoMax]

I'll add the ability to handle paired reads too.

Pre-filtering with reformat.sh in the mean time then?

[Brian Bushnell]

Does this work with paired end reads?

I get this error: Read of length 36 outside of range 50--1. Paired input is incompatible with 'breaklength'

Yep, looks like that's only for single reads, since I intended it for PacBio/Nanopore (breaking up long reads to a fixed length and discarding short dangling pieces)... I'll add the ability to handle paired reads too.

[darthsequencer]

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).

Does this work with paired end reads?

I get this error: Read of length 36 outside of range 50--1. Paired input is incompatible with 'breaklength'

[sk8bro]

It looks like the input is identically sorted by the way

[sk8bro]

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

[sebastian_1]

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.

[darthsequencer]

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.

[Brian Bushnell]

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.

[darthsequencer]

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?