How to take into account alternative bwa mem mapping when computing coverage

Question

when mapping short reads with bwa mem, if a read has alternative mapping positions they are reported by bwa mem in the X0 and XA tags. Now, let's say I want to compute the coverage of my bam file. How would I do for the alternative hits to be taken into accounts. What I mean is that, if a read QUERY has a selected best HIT position and 3 alternate HIT positions, I would like that QUERY read to be counted as 1 unit of coverage for the best HIT position, but also as 1 unit of coverage for each of the alternate HIT position.

My understanding is that samtools depth will only consider the best hit (which makes sense in 99.99% of the cases).

I guess I could edit the bam file with custom scripts so that the alternative hits are also reported as primary hits but I was wondering if there was an already existing tool (as I am not that confident with my ability to mess around with sam files).

Thanks a lot.

Answer 1

You would want to instead use the -a option with bwa mem to output all alignments. You would be best off filtering these with a custom script to only keep those with a close enough score to the primary hit. samtools depth can't be told to include secondary alignments (it's hard-coded), so you'll need to use something else like bamCoverage or samtools mpileup.