How To Get The Correct Alignment Statistics After Using Bwa Mem With -M Option

Question

[skm770]

I want to get uniquely mapped reads from bam files aligned with bwa-mem using the -M option (to make it compatible with picard-tools). If I try to use the following command with samtools to get uniquely mapping reads( or reads with mapping quality greater than 10. I get a very high number of around 90% uniquely mapping. I know I can decrease this number using a more strong threshold of mapping quality of 30 or so but I want to know what would be the ideal way to do so and if the bwa-mem needs to be run with different set of parameters to get more reliable alignments. Following are the bwa and samtools commands that I am using.

bwa mem -t 4  -M ref_genome.fasta  fastq[0] fastq[1] > samfile
samtools view -c -q10 -f 1 -F 12 -S samfile

Thanks

[brentp]

You can get non-secondary alignments with samtools view -F 0x100 (see: http://samtools.sourceforge.net/SAMv1.pdf)

Answer 1

[Pavel Senin]

I am not sure if I understand the purpose of the procedure aiming on the reliability of an alignment right, but you can make bwa to search for seeds in more rigorous way - exhaustively - by lowering the r paremeters which is set by default to 1.5. Alternatively or in addition, you could adjust A and B parameters to penalize mismatches more aggressively - if your reads are of good quality, then it makes sense not to align a read if there is mismatch, IMO. Then, you can scrutinize what is going on with gaps and tighten that too as well as clipping and unpaired read penalties.