[jweger1988]

Genomax,

Awesome, changing CIGAR fixed the issue. Thanks.

About Clumpify vs Dedupe, that's a good idea. I remember Brian saying somewhere that clumpify is not quite as good as Dedupe? What would be your recommended settings?

Thanks

[GenoMax]

Not sure if this is applicable but BBMap produces SAM v.1.4 CIGAR strings by default. If your variant calling program expects v.1.3 format then you can either add sam=1.3 option when you align (or reformat.sh in=your.bam out=new.bam sam=1.3).

Dedupe can be used with PE files but only if reads are interleaved. You may want to use clumpify instead of dedupe.

[jweger1988]

Hi Brian,

I have a couple of questions.

I've been comparing different alignment programs and then some different variant calling software. For some reason when I use a bam file produced with bbmap and try to call variants with lofreq I get an error message, but not with bam files produced with other programs (bowtie, mosaic, etc...). Even after sorting the bam I get the same message. Any thoughts on why this might be?

Also, I'm wondering if it's possible to use dedupe with the input as either paired fastq files (i.e. in1 and in2) or a bam file?

Thanks,
James

[Brian Bushnell]

Thanks a lot for your helpful replies. It is now working perfectly . One more question to the "xstag". I see that if the gene is on the + strand, R2 maps on+ and R1 maps on -. If I use "xstag=ss" I get for spliced R1 and R2 reads XS=+. So I think this is the correct library type. Am I right?

I think so. I've always found Tophat's description of those terms to be very confusing, so I have trouble remember for more than a few days which is the correct one to use and when.

Will you be incorporating the soft-clipped base trimming in this next update?

Yes, I have added that as of v37.23 which I just released. The flag is "trimclip", for example:

bbduk.sh in=mapped.sam out=trimmed.sam trimclip ordered

You can also output them as fastq if you want.

[jweger1988]

Brian,

Thanks for the response. My issue was exactly that. It seems to be working fine now by only using one at a time. The other issue was a hamming distance thing that I also figured out.

Will you be incorporating the soft-clipped base trimming in this next update?

Thanks again.

[Brian Bushnell]

Can you post an example of a read that is not getting trimmed with the command you posted?

Also, "ktrim=rl" with BBDuk won't work; it only allows "ktrim=r" or "ktrim=l". I've changed it now so in the next release it will print an error message if you use other values.

[jweger1988]

Hi Brian,

I'm trying to count some Kmers in a specific location in a read. I was hoping I could use bbduk or bbduk2 to trim around the sequence I want to count. I can't seem to make bbduk trim all the reads at these sites though. I'm trying many iterations of something like

bbduk2.sh in=my.fq.gz out=my_trim.fq.gz lliteral=CCACGCACTGATAACATCAC rliteral=CGGGGCACACTGCTTGTGGC k=20 hdist=2 tbo tpe

I've also tried this with bbduk specifying ktrim=rl and separating the sequences with a comma.

Is there a way to make bbduk(2) actually trim all of the sequences? I'm ending up with a bunch of the reads still having that sequence within a hdist=2.

You had previously suggested to reformat to a specific location in the reads but the problem with this is that even though they were all amplified with the same primers some additional nucleotides have managed to stay on the end of a subset of the reads. I can filter these out with size parameters but then I lose a lot of potentially valuable reads.

Any help is greatly appreciated.

Thanks!

[Mchicken]

Thanks a lot for your helpful replies. It is now working perfectly . One more question to the "xstag". I see that if the gene is on the + strand, R2 maps on+ and R1 maps on -. If I use "xstag=ss" I get for spliced R1 and R2 reads XS=+. So I think this is the correct library type. Am I right?

[Brian Bushnell]

As Genomax said, you need to set the "intronlen" flag, maybe to about 10. Arabidopsis has very short introns as I recall, so you shouldn't need to increase maxindel.

As for the xstag, you should not need to set that unless you plan on using the Tuxedo pipeline for analysis, but in that case... reads get XS:A:+ if R1 is mapped to the plus strand in firststrand mode. But, only spliced reads get an XS tag, so without setting intronlen, none of them will get it.

[GenoMax]

You should take a look at the maxindel= and intronlen= settings.

Sometimes it is not apparent, but if you did not specify a setting, it means that the program is using the default values for it.

This is what help for those two parameters says

maxindel=16000          Don't look for indels longer than this. Lower is faster.
                        Set to >=100k for RNAseq with long introns like mammals.
intronlen=999999999     Set to a lower number like 10 to change 'D' to 'N' in 
                        cigar strings for deletions of at least that length.

So those are the two values you used for this run.

[Mchicken]

Hi,

I try to align paired-end reads of Pisum sativum (pea) to the genome of Arabidopsis thaliana, as there is no reference genome for pea. I first tried HÌSAT2 but had alignment rates of ~4%. Therefore I wanted to give BBMap a try but have some questions:

1) I know we sequenced strand-specific but don't know the protocol but I see that the R2 read (second in pair) aligns always on the strand on which the gene is located and the R1 read (first in pair) on the opposite strand. How do I set the "xstag" parameter? xstag=fs or xstag=ss?

2) Regarding question 1 I tried both xstag=fs and xstag=ss but don't see a XS flag in the SAM file? Why?

3) When I aligned the first 1,250,000 read pairs (alignment rate ~12%) I saw from the CIGAR format (and IGV) that spliced reads (so normally sth. like 120N) are reported as having deletions (e.g. 120D). I thought BBMap was aware of spliced alignments?

Here is how I call BBMap:

sh /opt/bbmap/bbmap.sh \
mappedonly=t \
xstag=fs \
in=/.../test_R1.fastq \
in2=/.../test_R2.fastq \
out=/.../test_BBMap_fs.sam \
ref=/.../TAIR10_full.fa

I am really thankful for answers.

Best Regards
Mario

[sk8bro]

Thank you so much Brian!!! I will download it

[Brian Bushnell]

OK, I just released the latest version with the "deterministic" flag in 37.22.

[Brian Bushnell]

I just looked at the code, and it's trivial, there was already a variable:

protected static final boolean DYNAMIC_INSERT_LENGTH=true;

I just need to remove the word "final" and parse a new argument (in fact, I just did that). So that will be in the next release, probably Monday. Then you'll need to add these flags:

deterministic averagepairdist=100

...where instead of 100, you use the expected average inner distance between pairs (that's the insert size minus the sum of the read lengths, or the length of the unsequenced portion). You don't need the number to be correct, and you don't need to set it at all (the default is 100), but setting it accurately will yield slightly better pairing.

[sk8bro]

Hey Brian,

Do you have any sort of time estimate on when you might consider adding the deterministic mode to bbmap in a future release? It would help me decide what course of action to take if you had a rough (but realistic) time estimate!