calctruequality in bbmap

Question

[coyk]

I'm trying to recalibrate Q scores of a NextSeq run using MiSeq contigs assembled with Tadpole. The commands I used to map the reads to the reference were as follows:

bbmap.sh in=concatABC.fastq.gz outm=mapped.sam ref=./Lpe09_06TdpAssemblies/contigs09_06.fa ignorequality maxindel=100 minratio=0.4 ambig=toss qahist=qahist_raw.txt qhist=qhist_raw.txt mhist=mhist_raw.txt

The command I used to get the calibration matrices was:

calctruequality.sh in=mapped.sam

I get the following output

java -ea -Xmx57992m -Xms57992m -cp /home/me/bbmap/current/ jgi.CalcTrueQuality in=mapped.sam
Executing jgi.CalcTrueQuality [in=mapped.sam]

Exception in thread "Thread-2" Exception in thread "Thread-1" java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag and no ScafMap loaded.
          at stream.SamLine.toShortMatch(SamLine.java:1212)
          at stream.SamLine.toRead (SamLine.java:2015)
          at stream.SamLine.toRead (SamLine.java:1875)
          at stream.SamReadStreamer$ProcessThread.makeReads (SamReadStreamer.java:206) 
          at stream.SamReadStreamer$ProcessThread.run (SamReadStreamer.java:135) 
java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag and no ScafMap loaded.
          at stream.SamLine.toShortMatch(SamLine.java:1212)
          at stream.SamLine.toRead (SamLine.java:2015)
          at stream.SamLine.toRead (SamLine.java:1875)
          at stream.SamReadStreamer$ProcessThread.makeReads (SamReadStreamer.java:206) 
          at stream.SamReadStreamer$ProcessThread.run (SamReadStreamer.java:135) 
Exception in thread "Thread-3" java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag and no ScafMap loaded.
          at stream.SamLine.toShortMatch(SamLine.java:1212)
          at stream.SamLine.toRead (SamLine.java:2015)
          at stream.SamLine.toRead (SamLine.java:1875)
          at stream.SamReadStreamer$ProcessThread.makeReads (SamReadStreamer.java:206) 
          at stream.SamReadStreamer$ProcessThread.run (SamReadStreamer.java:135)

The program never quit so I don't know if it had stopped or was hanging but I let it sit there for over an hour. I didn't see /ref/qual/ files being generated. I'm not sure how to proceed. Since the mapped sam file was generated with bbmap using contigs from Tadpole, I'm not sure why there are Ms in the cigar strings. The infile I used for bbmap was made by first concatenating all Read1 from 3 libraries then all Read2 from 3 libraries then those two files concatenated.

[GenoMax]

Are you using the latest bbmap available?

Note: Reason I asked is this error had come up before in a different tool in BBMap and @Brian had released an update to fix the problem.

[coyk]

Thanks for responding so quickly and I apologize for not including that information in my original post. I swear I googled this and looked on SEQanswers for this error but didn't find it. Ah well. Anyway, It looks to be version 37.95, which is the version that was supposed to correct that error. However, I do see I'm a few versions behind. I suppose I should update and retry.

coyk

[GenoMax]

Do that and let us know if that fixes the problem.

[coyk]

I installed BBMap_38.0.0 and got a very similar result. There's still an "exception in thread" message about encountering a read with an M in cigar string. I added the ref=file command as the stout suggested and I got the same result. There does seem to be a little more output to stout that looks like reads were beginning to be processed, but I didn't see any files being generated. The link is to a screen shot of the output, I hope. I'm not sure how to use the image button to generate a picture.

coyk

[GenoMax]

Can you try one more thing:

reformat.sh in=your.sam out=new.sam sam=1.4

That should change all the M tags to sam v.1.4.

Note: I will make @Brian aware of this error. Unfortunately getting a hold of him has become rather difficult of late so let us see if we hear from him.

Note 2: Also post this as an issue on SF here. That should put it on @Brian's radar.

[coyk]

I should say, there were three outputs generated immediately: qhist_raw.txt, qahist_raw.txt, m_hist.txt.

[coyk]

I tried reformat.sh in=mapped.sam out=new.sam sam=1.4 and got similar results. The files that are supposed to be generated in /ref/qual/ appeared though they are small. There's still the message about the M in cigar string but it managed to process 3400 reads before the exception ReformatReads terminated in an error state; the output may be corrupt. A truncated new sam file was created. When I added ref= myref file to the command, the whole reformatting did take place and I got a new 90 GB sam file.

But alas.

calctruequality.sh in=mapped2.sam

and

calctruequality.sh in=mapped2.sam ref=myreference.fa

are still giving me the same M in cigar string error. And no /ref/qual files.

I went to the source forge site, but it's not clear to me where to report a bug. I see there's a discussion board for general topics, but the last entry was 7 months ago with no reply. I originally posted about this problem on GitHub and SeqAnswers.

https://postimg.cc/image/8g7wlo7av/

Answer 1

[Andreas]

I had the same issue just now with calctruequality.sh and a related one with reformat.sh. Fixing the BAM file as follows did the trick (adding MD, removing all non-aligned secondary etc. reads):

samtools view -h -F 3844 in.bam | samtools calmd  -  ref.fa | samtools view - -o out.bam 

Hope this helps, Andreas