[GenoMax]

You can try reformatting your alignment file using "reformat.sh in=you.bam out=new.bam sam=1.4" to see if that helps.

[kimng]

Bug in call variants with relation to sam files created from minimap2

Hello, I've been using different parts of the BBMap suite and mixing them with other softwares for WGS of bacterial samples (primarily for the purpose of automated QC). I'm a fan of callvariants.sh as it works quickly and doesn't require me to sort .sam files before hand and it's simple to clearfilters to return all values. Unfortunately when using callvariants.sh with .sam files created from minimap2 (https://github.com/lh3/minimap2) using illumina nextseq reads. I receive the error seen below. I was curious if this was a minor bug that could be corrected?

I'm running callvariants.sh from bbmap v37.90 obtained from conda (https://anaconda.org/bioconda/bbmap). Appreciate any assistance given.

-Kim

Error:
java -ea -Xmx48779m -Xms48779m -cp /home/ssi.ad/kimn/.conda/envs/env_kim/opt/bbmap-37.90/current/ var2.CallVariants in=aln.sam ref=contigs.fasta vcf=minimap.vcf ploidy=1 clearfilters
Executing var2.CallVariants [in=aln.sam, ref=contigs.fasta, vcf=minimap.vcf, ploidy=1, clearfilters]

Loading reference.
Time: 0.304 seconds.
Processing input files.
Exception in thread "Thread-3" Exception in thread "Thread-4" Exception in thread "Thread-5" java.lang.AssertionError: TODO: Encountered a read with 'M' in cigar string but no MD tag.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
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.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
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.
at stream.SamLine.toShortMatch(SamLine.java:1200)
at stream.SamLine.toRead(SamLine.java:1972)
at stream.SamLine.toRead(SamLine.java:1832)
at stream.SamReadStreamer$ProcessThread.makeReads(SamReadStreamer.java:206)
at stream.SamReadStreamer$ProcessThread.run(SamReadStreamer.java:135)

[mcmc]

mcmc- what are you setting the idfilter value to? Because reads still take a hit to their alignment score for soft-clipping. That is to say a read that's aligned with no mismatches but it was soft-clipped should have a lower score than a read that was mapped perfectly without soft-clipping. This strategy doesn't make sense when you're read is being soft-clipped due to adapter sequence at the end of the read...but for biological purposes it does make sense. Maybe a more appropriate question should be how to get rid of those "couple bases of adapter remaining" after trimming them.

yes I know the scores may be different... but I'm simply wondering whether the idfilter and subfilter post-filtering options work on the soft-clipped alignments or the non-clipped alignments. I want to be able to filter with, say, 2 MM or 98% ID *after* the soft clipping.

Thanks,
MCMC

[sdriscoll]

mcmc- what are you setting the idfilter value to? Because reads still take a hit to their alignment score for soft-clipping. That is to say a read that's aligned with no mismatches but it was soft-clipped should have a lower score than a read that was mapped perfectly without soft-clipping. This strategy doesn't make sense when you're read is being soft-clipped due to adapter sequence at the end of the read...but for biological purposes it does make sense. Maybe a more appropriate question should be how to get rid of those "couple bases of adapter remaining" after trimming them.

[mcmc]

bbmap local + idfilter

Brian et al., I am using the 'local' option with bbmap because I sometimes have a couple bases of adapter remaining (after bbduk mink=6) and my understanding is that 'local' will trim the ends to achieve a better score.

My question: does the idfilter option work on the local or global alignment? I have reads with 100% ID after trimming ends but they are not passing the idfilter. I'm curious whether there is another way to filter based on the local alignment percent id?

thanks!
MCMC

[jweger1988]

Hi Brian,

I'm using BBtools to call indels from small RNA virus genomes using both bbmap and callvariants. It's doing a great job of calling these. We are interested in the possibility of also calling inversion events from this data. I know this is somewhat common in transcriptomics data. Do these programs have this functionality?

To be more specific, is it possible to identify reads that are aligning to both the sense and antisense of the given reference?

Thanks in advance!

James

[santiagorevale]

So, I've spoken with 10X support. For now they are not doing or planning to do any filtering step on the rawdata. However, thanks to my email they've opened a software feature request. Running "cellranger count" without the Index file will just skip the metric on that file, nothing else.

In summary, I'm so sorry for making you waste your time! If appropriate, we could even delete all this posts because we have spoken more about 10X than bbduk.

Thanks GenoMax for your feedback! Cheers!

[santiagorevale]

Yes, that's how I'm running mkfastq. The sequencing is currently producing 2x75 bp, but by running mkfastq with the above mention cycle pattern: Read 1 is being trimmed to 26 bp, and for Read 2, we are only getting 75 bp, which is still enough data for doing alignment. Even though this is not the recommended way of sequencing this 10X libraries, 10X support team tell us they had other customers doing it this way with good results.

Regarding the N thing, that's not the case. These high N proportion reads are not being filtered at all by the first step (not by the bcl2fastq nor by the mkfastq wrapper).

I have contacted 10X and they are saying that this reads are not being filtered, but they are not being mapped also thus downstream analysis is not being affected, which is something we already knew.

However, they told me that for v2 chemistry the Index read is not being used any more after the demultiplexing step. I replied them asking how are they calculating "Q30 Bases in Sample Index" metric then without that file. In the meanwhile, I'm running a test so I'll be back later with both the results and the reply. So lets see. If this works, then I won't be needing any change on bbduk. Lets hope for the best!

Cheers!

[GenoMax]

So how are you doing the base calling then?

10x recommends that you run 26x8x98 but you are running this as a 75x8x75 run instead? If that is the case my recommendation above would cut the read 1 down to 26 as expected by 10x. The index read is correct length and the last read, which should be 98 bp is only 75 in your case. Is this how you are running cellranger mkfq?

With above base mask you should get reads that don't have any N's. Or so I would think.

[santiagorevale]

Oh, sorry if I was not clear enough. Those are the sequencing cycles not how I'm doing the base calling.

I've been working on 10X projects for several months now and evaluating them in different platforms. That's why I also included the additional question above.

[GenoMax]

Yikes! You are certainly doing something completely off-label here :-)

In the cellranger mkfastq run, you could use a base-mask such as --use-bases-mask=Y26n*,I8,Y*. That should produce reads in the format 10x wants. If this is a 2D, 4000 run then you may have to do --use-bases-mask=Y26n*,I8,n*,Y*, Your second read is going to be short by about 20 bases. I am not sure if the software will like that.

@Brian would have to comment on your original request. For him, being a programmer, anything would be possible to implement.

[santiagorevale]

Regarding the issue that produces the Ns, it's not supported by them, because it occurs when sequencing on a HiSeq 4000 using 75,8,0,75 bp cycle pattern instead of their recommended 26,8,98 bp pattern.

Regarding the N filtering (or any other quality filtering) I agree with you that they should be doing it in advance and incorporate it into their pipeline. I know that so far it's not being done and I have just sent them an email about it.

I'll see what their answer is and be back to you with it. Let's see what their feedback is.

In the meantime, do you think this would be something complicated to implement? If not as a permanent thing on bbduk, could you walk me on how to make it work for me? Because I was unable to match your scripts performance (I tried coding things in perl, python and bash).

Thanks for this quick replies!

[GenoMax]

While I have not worked with cellranger (I have only used longranger for WGS data) this sounds like an odd behavior. 10x should be accounting for presence of N's in the reads in their software.

If you have not talked with their software tech support then I would suggest that you give that a try. They may have some specific suggestions or will need to implement a fix in their software. Support has been responsive to me in past.

[santiagorevale]

cellranger pipeline is splitted in two steps: "cellranger mkfastq" which produces the three above mentioned files (R1, R2, I1) and "cellranger count" which uses this 3 files as input.

The reason why I need to do the filtering first is that some platforms fail sometimes when calling reads from the TOP surface of the flowcell for Read 2, resulting in reads which will be fully or partly composed by Ns. Then, "cellranger count" will use all of the reads (including those that are 100% Ns) to do their calculations, ending up with aberrant metrics (like "mean reads per cell", which is determined as the raw number of reads divided by the number of called cells).

[GenoMax]

If you are using the cellranger pipeline then you should not need to do anything to the data outside the pipeline. Cellranger pipeline takes in raw Illumina flowcell data folder and does the entire analysis (including demultiplexing/alignments/counting etc).