[nlubock]

BBMap Error

Hi Brian,

I'm getting the following error when using BBMap (Version 36.11) for some of my alignments:

java -Djava.library.path=/usr/local/stow/bbmap_36.11/bin/jni/ -ea -Xmx8g -cp /usr/local/stow/bbmap_36.11/bin/current/ align2.BBMap build=1 overwrite=true fastareadlen=500 in=1_nonDoped.adapters.fastq vslow=t k=8 ref=../ref-seq.fasta mdtag=t -Xmx8g outm=bbmap.sam nodisk
Executing align2.BBMap [tipsearch=150, minhits=1, minratio=0.25, rescuemismatches=50, rescuedist=3000, build=1, overwrite=true, fastareadlen=500, in=1_nonDoped.adapters.fastq, k=8, ref=../ref-seq.fasta, mdtag=t, -Xmx8g, outm=bbmap.sam, nodisk]

BBMap version 36.11
Set MINIMUM_ALIGNMENT_SCORE_RATIO to 0.250
Retaining first best site only for ambiguous mappings.
Executing dna.FastaToChromArrays2 [../ref-seq.fasta, 1, writeinthread=false, genscaffoldinfo=true, retain, waitforwriting=false, gz=true, maxlen=536670912, writechroms=false, minscaf=1, midpad=300, startpad=8000, stoppad=8000, nodisk=true]

Set genScaffoldInfo=true
Set genome to 1

Loaded Reference:	0.001 seconds.
Loading index for chunk 1-1, build 1
Indexing threads started for block 0-1
Indexing threads finished for block 0-1
Generated Index:	0.013 seconds.
Analyzed Index:   	0.031 seconds.
Started output stream:	0.022 seconds.
Cleared Memory:    	0.123 seconds.
Processing reads in single-ended mode.
Started read stream.
Started 40 mapping threads.
Exception in thread "Thread-25" java.lang.ArrayIndexOutOfBoundsException: 136
	at align2.BBIndex.getOffsetArray(BBIndex.java:3065)
	at align2.BBIndex.shrink2(BBIndex.java:830)
	at align2.BBIndex.find(BBIndex.java:443)
	at align2.BBIndex.findAdvanced(BBIndex.java:396)
	at align2.AbstractMapThread.quickMap(AbstractMapThread.java:736)
	at align2.BBMapThread.processRead(BBMapThread.java:405)
	at align2.AbstractMapThread.run(AbstractMapThread.java:495)
Detecting finished threads: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39

**************************************************************************
Warning!  1 mapping thread did not terminate normally.
Check the error log; the output may be corrupt or incomplete.
Please submit the full stderr output as a bug report, not just this message.
**************************************************************************




   ------------------   Results   ------------------

Genome:                	1
Key Length:            	8
Max Indel:             	16000
Minimum Score Ratio:  	0.25
Mapping Mode:         	normal
Reads Used:           	570587	(170899784 bases)

Mapping:          	34.755 seconds.
Reads/sec:       	16417.62
kBases/sec:      	4917.33


Read 1 data:      	pct reads	num reads 	pct bases	   num bases

mapped:          	 99.9855% 	   570504 	 99.9840% 	   170872487
unambiguous:     	 99.9855% 	   570504 	 99.9840% 	   170872487
ambiguous:       	  0.0000% 	        0 	  0.0000% 	           0
low-Q discards:  	  0.0000% 	        0 	  0.0000% 	           0

perfect best site:	 62.1351% 	   354535 	 62.2356% 	   106360500
semiperfect site:	 62.1351% 	   354535 	 62.2356% 	   106360500

Match Rate:      	      NA 	       NA 	 99.7089% 	   170661054
Error Rate:      	 37.4881% 	   213882 	  0.2887% 	      494106
Sub Rate:        	 28.9387% 	   165105 	  0.1164% 	      199217
Del Rate:        	 10.8304% 	    61791 	  0.1676% 	      286801
Ins Rate:        	  1.1807% 	     6736 	  0.0047% 	        8088
N Rate:          	  0.6976% 	     3980 	  0.0024% 	        4128
Exception in thread "main" java.lang.AssertionError:
The number of reads out does not add up to the number of reads in.
This may indicate that a mapping thread crashed.
If you submit a bug report, include the entire console output, not just this error message.
570504+0+0+82+0 = 570586 != 570587
	at align2.AbstractMapper.printOutput(AbstractMapper.java:1867)
	at align2.BBMap.testSpeed(BBMap.java:489)
	at align2.BBMap.main(BBMap.java:33)

Any insights would be greatly appreciated!

[Brian Bushnell]

Hi SC,

No publication yet, so for now, please cite it like this:

"BBMap - Bushnell B. - sourceforge.net/projects/bbmap/"

...or similar, depending on the style guidelines.

[schalivendra]

BBMap

Hi,

How I can cite BBMap? Is there a publication that I can list under references?

thanks very much,
SC

[Brian Bushnell]

Hi Kate,

Yes, please do that - it would be very helpful.

-Brian

[sk8bro]

Thanks so much! I can try to replicate the issue and email you a small dataset/reference if that would help? Just let me know. BBsplit is otherwise really fantastic!

Kate

[Brian Bushnell]

Hi Kate,

That's odd... I will investigate and see if I can figure out why that's happening. minid and idfilter are slightly different, but certainly with idfilter=0.95 you should not get any alignments with under 95% identity.

-Brian

[sk8bro]

hey Brian,

I am still trying to get rid of these alignments, and others that don't make sense coming through at 95% identity of both reads. Changing kmer size, pairlen, rescuedist, ambiguous2 and maxindel back to default as well as removing vslow had no effect. However, changing minid from unspecified/default of 0.76 to 0.95 did change the output such that these alignments are no longer outputted.

I didn't think minid should effect output with idfilter set, only speed. Or in other words, that idfilter trumps at the very end.

Do you think there could be a bug? Thanks! Kate

[sk8bro]

idfilter?

Hey Brian,

I have some unexpected BBsplit mapping results. Do you have any insights?

Here is my command:

bbsplit.sh k=8 in=1P.fastq\
in2=2P.fastq\
basename=blah\ outu=clean_1P.fastq\
outu2=clean_2P.fastq\
ref=path_to_fasta\
ambiguous2=split pairlen=200000 rescuedist=200\ refstats=stats.txt\
idfilter=0.95 maxindel=20 nzo=f vslow pairedonly 2>> BBsplit_stats.txt


I don't understand why, for example, this read pair would align to the reference that BBsplit reports, because neither read actually has 95% identity to the reference. See attached Geneious alignment of both Reads to Reference sequence. Note: these reads are a little weird in that there isn't an insert, but still this is unexpected behavior that I'd like to understand and remedy. Thanks!




------------------ Results ------------------

Genome: 1
Key Length: 8
Max Indel: 20
Minimum Score Ratio: 0.56
Mapping Mode: normal
Reads Used: 2215982 (265861975 bases)

Mapping: 38.699 seconds.
Reads/sec: 57261.41
kBases/sec: 6869.93


Pairing data: pct reads num reads pct bases num bases

mated pairs: 0.0001% 1 0.0001% 210
bad pairs: 0.0000% 0 0.0000% 0
insert size avg: 51.00


Read 1 data: pct reads num reads pct bases num bases

mapped: 0.0001% 1 0.0001% 105
unambiguous: 0.0001% 1 0.0001% 105
ambiguous: 0.0000% 0 0.0000% 0
low-Q discards: 0.0000% 0 0.0000% 0

perfect best site: 0.0000% 0 0.0000% 0
semiperfect site: 0.0000% 0 0.0000% 0
rescued: 0.0000% 0

Match Rate: NA NA 49.4444% 89
Error Rate: 100.0000% 1 50.5556% 91
Sub Rate: 100.0000% 1 8.8889% 16
Del Rate: 100.0000% 1 41.6667% 75
Ins Rate: 0.0000% 0 0.0000% 0
N Rate: 0.0000% 0 0.0000% 0


Read 2 data: pct reads num reads pct bases num bases

mapped: 0.0001% 1 0.0001% 106
unambiguous: 0.0001% 1 0.0001% 106
ambiguous: 0.0000% 0 0.0000% 0
low-Q discards: 0.0000% 0 0.0000% 0

perfect best site: 0.0000% 0 0.0000% 0
semiperfect site: 0.0000% 0 0.0000% 0
rescued: 0.0000% 0

Match Rate: NA NA 49.1713% 89
Error Rate: 0.0096% 1 50.8287% 92
Sub Rate: 0.0096% 1 9.3923% 17
Del Rate: 0.0096% 1 41.4365% 75
Ins Rate: 0.0000% 0 0.0000% 0
N Rate: 0.0000% 0 0.0000% 0

[Brian Bushnell]

Hi telia,

The command depends on the operating system you are using. Assuming you unzipped the package to /home/BWA_sam/countnew - meaning there is a file in that directory called "pileup.sh" (which is normally in a directory called "bbmap") - for linux, it would be:

/home/BWA_sam/countnew/pileup.sh in=tube1_dedup.sam stats=stats.txt -Xmx30g

If you unzip the package normally to "/home/BWA_sam/countnew/", pileup.sh should be at "/home/BWA_sam/countnew/bbmap/pileup.sh"

You can also skip the shellscript and use this command, which works in any operating system (as long as the -cp argument points to the location of the directory named "current"):

java -ea -Xmx30G -cp /home/BWA_sam/countnew/current jgi.CoveragePileup jgi.Pileupin=tube1_dedup.sam stats=stats.txt

These are equivalent, but the shellscript is shorter if you are using Linux/bash.

[GenoMax]

Dear Brain,
thank you for your attention. I am new in Bioinformatic. I have some sam file which I made using BWA. I am going to count unique mapped reads for gene expression. I have tried several times to work with Pileup.sh, but I had an error like this: Could not find or load main class pileup.sh
after that I used this command :
java -ea -Xmx30G -cp /home/BWA_sam/countnew/Pileup.sh in=tube1_dedup.sam stats=stats.txt

but unfortunately I have this error....
Error: Could not find or load main class in=tube1_dedup.sam

any help is much appreciated
Best regards

Did you move pileup.sh script out of the BBMap folder? Try running it from its original location.

[telia]

error in pileup.sh

Dear Brain,
thank you for your attention. I am new in Bioinformatic. I have some sam file which I made using BWA. I am going to count unique mapped reads for gene expression. I have tried several times to work with Pileup.sh, but I had an error like this: Could not find or load main class pileup.sh
after that I used this command :
java -ea -Xmx30G -cp /home/BWA_sam/countnew/Pileup.sh in=tube1_dedup.sam stats=stats.txt

but unfortunately I have this error....
Error: Could not find or load main class in=tube1_dedup.sam

any help is much appreciated
Best regards

[skbrimer]

Hi Brian,

Thanks for the input, I'm glad to hear that I was interpreting the outputs correctly as a biologist and not a computer engineer. I agree with our interpretation of the results as well. I understand that the way I uses this comparison was not overall fair, but mostly I was looking to see if there was a large difference between the final alignment of the type types of aligners and I think in that manner this test worked well.

I was afraid that because most of the viral species we are studying make their different proteins by differential spice events or from subgenomic start sites that TMAP made have been making false positive calls for variation. From this test I feel that it was not but I also feel that using a splice-aware aligner would be more prudent for this work. It looks like BBMap was able to find splice sites (real or not because some of them are in areas that there is no information on them yet) and completed the alignment using the longer read segments, whereas TMAP overcame this issue using the smaller segments or clipping the longer ones until it matched the one area.

Thank you for getting back to me, this makes me feel more confident about my interpretation of the data.

[Brian Bushnell]

It's hard to directly compare a splice-aware versus non-splice-aware aligner directly on data that may contain splice-like features since you don't really know which is right. Also, viral data is very noisy due to the rapid mutation rate, or something; maybe the RNA forms catalytic structures and the reads edit each other in solution. And finally, the coverage is tricky to calculate using a output from a splice-aware aligner, because the splices may or may not be counted - if they are counted, the coverage is artificially inflated.

BBMap has a tool called Pileup for calculating coverage. You can use it like this:

pileup.sh in=mapped.bam delcov=f covstats=stats.txt basecov=perbasecov.txt bincov=binnedcov.txt

delcov=f will ignore deletions (such as splice sites) in coverage calculation; delcov=t will include them.

The main factor in the difference in mapping rates is probably clipping. Note that TMAP clipped 25% of the reads, which allows possibly chimeric reads to map; BBMap does not do clipping by default, but you can turn it on by adding the "local=t" flag.

The mapping quality is not important in this case. You can't compare it between aligners since the quality values are subjective and the implementations are different. In fact, TMAP's average mapq of over 70 is fairly absurd; normally they don't go above 50 (by definition, a mapq of 70 means a 1/10,000,000 chance of being wrong, which would be difficult to achieve in practice). But anyway, that's just the aligner's opinion of how likely it is that an alignment is correct. Mapq is sometimes useful for filtering reads, but I don't think average mapq is ever useful.

The error rate here is being miscalculated due to the "splices" - the long deletions are being calculated as errors. Also, the fact that BBMap is not clipping reads that have poorly-matching ends inflates the error rate.

Perhaps the best way to do an apples-to-apples comparison is to turn off BBMap's long-indel detection and turn on soft-clipping:

bbmap.sh local=t maxindel=20 [other arguments]

That will make their outputs more comparable. It won't tell you which is correct, though.

Anyway... it's possible that those things that look like splices, really are splices - not traditional ones defining introns like in a eukaryote, but the RNA cutting itself or other pieces and randomly recombining, or inserting itself into other nearby pieces of viral RNA. I'm not a biologist and I have not studied viral behavior, but I have seen a lot of results of viral mapping that looked like this (with what appear to be large splices all over), and assembled isolate viral genomes that were supposed to be uniform and exactly like the reference, but turned out to have tons of different versions with hypervariable regions, and so forth. They're often really ugly.

[skbrimer]

Hi Brian,

I have been playing with BBmap to see if there is any difference between mapping virus cDNA (RNAseq) with a splice aware aligner vs the default TMAP (Ion Torrent) and I was hoping you can help me understand what I am seeing.

I can see that BBmap is saying I have reads that could be splices, how do I check this? When compared to the TMAP output the overall output looks the same. Same shape in the coverage histogram on IGV, mostly the same SNP calls.

When I evaluated them with Qualimap, that is where I can see a lot more differences. BBmap mapped less reads but had more coverage and more uniform coverage then TMAP. However TMAP had a higher mapping quality and lower overall error rate. I think these last two indexes inflated in the BBmap report because of the splicing but that is why I am asking. Are these two mappers equivalent in this case?

Also I wanted to include a SNP report for you to look as well but freebayes was unable to locate the index for the ref when trying to call SNPs with the BBMap output. Since it makes its on ref index I think I should use that but I'm not sure how to call it.

I would appreciate any insight from you when you get a chance.

Thank you,

Sean

[horvathdp]

Ahh...OK, thanks for using the small words and simple pictures for my old brain :-). I think I understand the issue. I'll play with this a bit and see what I can teach myself.