`quickbin.sh` yields 0 coverage for all contigs when inputting `bwa mem` SAM files

[ucassee]

Hi QuickBin team,

I ran into an issue with quickbin.sh yields 0 coverage for all contigs when I directly input SAM files generated by bwa-mem2. However, the SAM files generated by bbmap.sh work perfectly.

Could you explain the key formatting differences between these two SAM outputs that might cause this? Also, how can I make the bwa SAM files compatible with quickbin.sh?

I attached the two types of SAM files.

250714-H3-PC3-2_1k_bbmap_sam.txt
250714-H3-PC3-2_1k_bwa_sam.txt

Thanks for your help!

[btpsj]

Are the sam files you added only a small portion of the actual output? Most of the reads look unmapped. I can't reproduce with any quickbin.sh commands since I don't have access to the reference, so I tried it on some simulated data, instructions below.

One potential issue, is that bwa-mem2 outputs CIGAR strings wtih M (alignment match) but quickbin.sh only processes = or X . So in the below example, I do get similar a similar pathology but hard to say its exactly what you are running into since the few reads you sent are unmapped and don't have CIGAR strings.

Assuming that is the problem, and you absolutely must use the bwa-mem2 SAM outputs, a quick fix you can do is write a script that goes through each alignment and just swaps M for =. I tried that and it will "fix" the problem with the caveat that swapping it can potentially cause downstream issues if any other tools you use need take into account the = vs X.

Notes for reproducing

randomgenome.sh len=10000 chroms=5 out=bbcontigs.fa
randomreadsmg.sh bbcontigs.fa reads=100 out=bbreads.1.fq.gz out2=bbreads.2.fq.gz
bwa-mem2 index bbcontigs.fa
bwa-mem2 mem bbcontigs.fa bbreads.1.fq.gz bbreads.2.fq.gz >bbreads.bwamem2.sam
bbmap.sh ref=bbcontigs.fa in=bbreads.1.fq.gz in2=bbreads.2.fq.gz out=bbreads.bbmap.sam

quickbin.sh in=bbcontigs.fa reads=bbreads.bwamem2.sam covout=bwamem2.cov.txt ow=t
quickbin.sh in=bbcontigs.fa reads=bbreads.bbmap.sam covout=bbmap.cov.txt ow=t

Outputs

% cat bbmap.cov.txt 
#Contigs        5
#Depths 1
#ShortName      ID      Size    Cov_0   Edge    Weight
contig1 0       2000    3.30
contig2 1       2000    3.34
contig3 2       2000    3.32
contig4 3       2000    3.31
contig5 4       2000    3.34
% cat bwamem2.cov.txt 
#Contigs        5
#Depths 1
#ShortName      ID      Size    Cov_0   Edge    Weight
contig1 0       2000    0
contig2 1       2000    0
contig3 2       2000    0
contig4 3       2000    0
contig5 4       2000    0

Swapping M for =

% awk 'BEGIN{OFS="\t"} /^@/{print; next} {gsub(/M/, "=", $6); print}' bbreads.bwamem2.sam > swapped.sam
% quickbin.sh in=bbcontigs.fa reads=swapped.sam covout=swapped.cov.txt ow=t
% cat swapped.cov.txt 
#Contigs        5
#Depths 1
#ShortName      ID      Size    Cov_0   Edge    Weight
contig1 0       2000    3.30
contig2 1       2000    3.34
contig3 2       2000    3.32
contig4 3       2000    3.31
contig5 4       2000    3.34

[ucassee]

Hi @btpsj ,

Thank you so much for the detailed analysis and reproduction! The SAM files I attached earlier were randomly subsampled to 100 reads (the original BAM files are huge), which explains why they contained mostly unmapped reads.

You've pinpointed the exact issue — the M vs = in CIGAR strings. After swapping M to = with awk, quickbin.sh works perfectly.

However, for my metagenomic data (hundreds of gigabytes per BAM file), using an awk script to replace M with = on the entire file is very slow and, more importantly, requires huge temporary storage for the converted file. Is it possible to make quickbin.sh natively accept the M operator? Supporting M would make quickbin.sh directly compatible with bwa-mem2 (and many other mainstream aligners) without any preprocessing.

Thank you again for your support!

[btpsj]

A recent fix (46607ed) should now add support for the different CIGAR strings. I can now run quickbin.sh and it outputs a coverage text file that is the same as the output for bbmap.sh in the example I posted above. Try out the latest version and let me know if that works for you!