Get the coverage after read alignment with bwa-mem2

Question

[maxime.policarpo]

Hi,

I have aligned genomic reads to a reference genome using bwa-mem2 (using bwa-mem2 mem). From what I understood, if a read map to 2 different locations in the genome with the same score, bwa-mem2 will choose at random one of the location and this will have a MAPQ score of 0 (I did not run mem with the -a option).

What I want to do then is to get the mean read coverage at different part of the genome using bedtools (Thats why I absolutely need reads to be mapped to only 1 genomic location) and compare these different coverages.

However, I have several questions :

• Is my first asumption about bwa-mem correct ?
• Will bedtools coverage take into account reads flagged with a MAPQ of 0 ?

Thanks for any help and tips !

Maxime

[Pappu]

You may check out deeptools.

Answer 1

[Pierre Lindenbaum]

1) yes

2)

$ samtools depth -a -r "chr1:123-456" in.bam -q 0  | awk '{D+=int($3)} END{print D/(NR*1.0);}'

[maxime.policarpo]

Thanks a lot for your answer :D