[Brian Bushnell]

You can output fastq to outm and sam to out or outu (or any similar combination), but you can't output both sam and fastq to outm.

[Brian Bushnell]

If I set

 outm=normal_pairs.fq

, I will have interleaved paired end ?

Yes, they will be interleaved. If you want you can use "outm1=r1.fq outm2=r2.fq" to keep them in 2 files but I find interleaved reads more convenient.

[Brian Bushnell]

Basically, you can do something like this:

bbmap.sh in=reads.fq outm=normal_pairs.sam outu=abnormal.fq killbadpairs=t pairedonly=t

Then the proper pairs will be in normal_pairs.sam, and everything else will be in abnormal.fq. You can subsequently apply further filters to abnormal.fq to look at it in more detail. Either file can put output in fastq or sam (depending on the extension) but everything in outu will have its mapping information removed so if you want to see the orientations, you'd need to remap them. You could in a subsequent pass, do this:

bbmap.sh in=abnormal.fq outm=strange_pairs.sam outu=unmapped.fq requirecorrectstrand=f pairedonly=t

That would give you FF and RR pairs in strange_pairs.sam and all the unmapped or half-mapped pairs in unmapped.fq. You can further split them like this:

bbmap.sh in=unmapped.fq outm=half_mapped.sam outu=both_unmapped.fq

[GenoMax]

@moistplus: Take a look at these options for BBMap and their default values.

samestrandpairs=f       (ssp) Specify whether paired reads should map to the
                        same strand or opposite strands.
requirecorrectstrand=t  (rcs) Forbid pairing of reads without correct strand 
                        orientation.  Set to false for long-mate-pair libraries.
killbadpairs=f          (kbp) If a read pair is mapped with an inappropriate
                        insert size or orientation, the read with the lower  
                        mapping quality is marked unmapped.
pairedonly=f            (po) Treat unpaired reads as unmapped.  Thus they will 
                        be sent to 'outu' but not 'outm'.

[Brian Bushnell]

There is no equivalent to -M, since BBMap does not produce split alignments. Also, by default, BBMap does not consider a pair "proper" unless it aligns in the fragment orientation, with reads on opposite strands. if you are trying to quantify the insert size distribution of a long-mate-pair library with a different orientation, you should add the flag "requirecorrectstrand=f".

To generate an insert size histogram, you can do this:

(for interleaved reads)
bbmap.sh ref=ref.fa in=reads.fq ihist=ihist.txt

(for reads in two files)
bbmap.sh ref=ref.fa in1=r1.fq in2=r2.fq ihist=ihist.txt

[GenoMax]

is it possible to estimate insert size of non overlapped library ?

By mapping with bbmap (assuming you have a reference).

[Brian Bushnell]

BBMerge's "minentropy" is an approximation; it generates a score based on counting the number of unique 3-mers in the head and tail of the read to adjust the minimum allowed overlap. It can be disabled with the flag "entropy=f". Specifically, the score is the sum of 4*(kmers occurring at least once)+1*(kmers occurring at least twice).

For normal (complex) sequence, minoverlap will be used as the minimum allowed overlap to merge reads. But if there is a minentropy setting, the head (3' end) and tail (5' end) of the read are examined and the minimum overlap may be increased (never decreased) if the sequence is low-complexity, until the entropy score reaches a minimum value. For example, if one read ends with "...ATCGTTAGTCCCCCCC....CCCCCCC", the poly-C part will contribute only a single unique 3-mer no matter how much there is, so the score will remain at 5 until it gets to the complex part, at which the score will start to rapidly increase due to new unique 3-mers. When the score passes a threshold (default 39, corresponding to roughly 10 unique 3-mers) the minimum overlap length is determined.

This reduces the false-positive merge rate because low-complexity sequence is much more prone to coincidental matches than high-complexity sequence, so it makes more sense to determine the minimum overlap dynamically based on information content than to use the same fixed number for all reads. You can manually set the entropy cutoff with the "minentropy=" flag on the command line, but it's not documented because it's hard to explain concisely. Rather, the different stringencies (strict, loose, etc) each have their own default value for minentropy, with higher values for higher stringency. For your reads, the minentropy flag may be important because they are so low-complexity.

[sganschow]

Ah I see. That was really helpful, thank you Brian.
Since I know the shortest and longest expected amplicon, it was possible to get perfectly merged reads by setting the correct values for minoverlap0, minoverlap, mininsert0, and mininsert.

One last question: what does minentropy stand for? I haven't seen it in the bbmerge usage.

[Brian Bushnell]

Maybe it just depends on the parameters in bbmerge?
Thanks in advance for your help.
Sebastian

Indeed, it does depend on the parameters. They will merge at 311bp if you use these parameters:

bbmerge.sh in=amplicons.fa loose minoverlap0=14

..or anything less stringent, such as

bbmerge.sh in=amplicons.fa xloose minoverlap0=50

...which would yield a much higher overall merge rate for this kind of read pair. The problem is that in this case the last 13 bp of the reads look like this:

AAAGAAAGAAAGA

TCTTTCTTTCTTT

...which, once you reverse-complement read 2, look like this:

AAAGAAAGAAAGA
|||||||||||||
AAAGAAAGAAAGA

They match perfectly, so there's no way to determine that that isn't the correct overlap frame based on sequence information alone, without additional information about the expected insert size or overlap length.

[sganschow]

Thanks for your reply.
I followed your advise to simulate different repeat lengths and generated paired end reads. Indeed, when adding more repeat units, there's a point where the pair is not joined but classified as ambiguous. I could manage to merge some longer amplicons by setting loose=t.

However, it is not clear to my why a pair like the following can't be merged, although the reads have an overlap of ~190nt and the repeat region has a complex structure, i.e. repeat stretches are interspersed by single bases and non-repetitive elements. Merging by hand (aligning R1 and rev-comp R2) revealed that there should be only 1 solution without any mismatches.

Please have a look at the following read pair (2x250nt):

Read 1

>31.2_14/17
AATCTGGGCGACAAGAGTGAAACTCCGTCAAAAGAAAGAAAGAAAGAGACAAAGAGAGTTAGAAAGAAAGAAAGAGAGAGAGAGAGAAAGGAAGAAAGGAAGAAAAAGAAAGAAAAAGAAAGAAAGAGAAAGAAAGAAAGAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAA

Read 2

>31.2_14/17
ACATCTCCCCTACCGCTATAGTAACTTGCTCTTTCTTTCCTTCCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTCTTTCTTTCTTTCTCTTTCTTTCTTTTTCTTTCTTTTTCTTCCTTTCTTCCTTTCTCTCTCTCTCTCTTTCTTTCTTTC

Amplicon (311nt)

AATCTGGGCGACAAGAGTGAAACTCCGTCAAAAGAAAGAAAGAAAGAGACAAAGAGAGTTAGAAAGAAAGAAAGAGAGAGAGAGAGAAAGGAAGAAAGGAAGAAAAAGAAAGAAAAAGAAAGAAAGAGAAAGAAAGAAAGAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGGAAGGAAAGAAAGAGCAAGTTACTATAGCGGTAGGGGAGATGT

Maybe it just depends on the parameters in bbmerge?
Thanks in advance for your help.
Sebastian

[Brian Bushnell]

If reads overlap only in a repetitive region, they will not be merged. BBMerge looks at all the possible overlaps, and keeps track of the two top-scoring ones (based on length and match/mismatch ratio). If those two are close, the pair will be classified as "ambiguous" and not merged. However, that does not mean it will always be correct; say you have a tandem repeat of two copies, like this:

ARRB

...where A and B are unique, and R is repeat. If the reads look like this:

1: ARR
2: RRB

...then there are 2 good overlap frames, forming ARRB or ARRRB, so the merge will be rejected as ambiguous. But if the reads only span a single repeat each, like this:

1: AR
2: RB

...then there will only be one apparently good overlap frame, and the reads will probably be merged incorrectly to form ARB. BBMerge's false-positive merge rate is extremely low, but it's not perfect. With shotgun data you can add the flag "rsem" to greatly reduce the chances of false-positive merges due to short tandem repeats, but that does not really work with amplicon data. You may want to simulate some data using your expected sequence to see what the actual behavior is.

[sganschow]

Hello,

how does BBMerge behave when the reads contain repetitive regions at the right end?

My amplicons are variable in length and derive from STR's, meaning that they are like:
(non-repetitive flanking region) - (tandem repeats) - (non-repetitive flanking region).
If the the amplicon is long enough, I could imagine that there is a case where paired reads overlap only in the repetitive region and thus multiple ways of merging are theoretically correct.
Until now, merging works fine. Can BBMerge ensure that merged reads are always consistent with the "real" amplicon sequence?

Sebastian

[shimingt]

Dear Genomax,

Thanks for your help!

[GenoMax]

Can you try

for x in *_R1_001.fastq;do bbmerge.sh -Xmx20g in1=$x in2=${x%_R1_001.*}_R2_001.fastq out=bbmerge\/${x%_*-*_R1_001.*}_R1R2_bbmerge.fastq outu1=outu1\/$x outu2=outu2\/${x%_R1_001.*}_R2_001.fastq > bbmerge\/${x%_*-*_L001_R1_001.*}_bbmerge_log.txt 2>&1;done

[shimingt]

Log file from bbmerge

Hello GenoMax,

Thanks for your reply.

I tried a command like this:

for x in *_R1_001.fastq;do echo bbmerge.sh -Xmx20g in1=$x in2=${x%_R1_001.*}_R2_001.fastq out=bbmerge\/${x%_*-*_R1_001.*}_R1R2_bbmerge.fastq outu1=outu1\/$x outu2=outu2\/${x%_R1_001.*}_R2_001.fastq > bbmerge\/${x%_*-*_L001_R1_001.*}_bbmerge_log.txt 2>&1 >> bbmerge.sh;done


The log files are generated according to the names, but the log files are empty.

Am I doing something wrong here?